LYAPUNOV AND DIFFUSION TIMESCALES IN THE SOLAR NEIGHBORHOOD

Chaotic dynamics due to interaction of nonlinear resonances in Hamiltonian systems is studied in a broad range of application areas, from plasma physics to celestial mechanics (e.g., Chirikov 1979; Lichtenberg & Lieberman 1992; Abdullaev 2006). The characteristic time of predictability of any motion is nothing but the Lyapunov time (the inverse of the maximum Lyapunov exponent) of the motion. Generally, the estimation of the Lyapunov exponents is one of the most important tools in the study of chaotic motion (Lichtenberg & Lieberman 1992), in particular in celestial mechanics. The Lyapunov exponents characterize the mean rate of exponential divergence of trajectories close to each other in phase space. Non-zero Lyapunov exponents indicate chaotic character of motion, while the maximum Lyapunov exponent equal to zero signifies regular (periodic or quasi-periodic) motion.

The development of methods of numerical computation of the Lyapunov exponents has more than a 30 year history (e.g., Froeschlé 1984; Lichtenberg & Lieberman 1992). In contrast, methods of analytical estimation of the Lyapunov exponents started to be developed only recently (Holman & Murray 1996; Murray & Holman 1997; Shevchenko 2000a, 2002, 2007, 2008b).

In studies of the dynamics of the Milky Way, analysis of the Lyapunov exponents has not yet been widely used, but nevertheless there are important achievements. Fux (2001) used Lyapunov exponents as a tool to find bar-induced manifestations of chaos in the local disk stellar kinematics. Taking into account the perturbations from the bar solely (for some particular bar strengths), he constructed "Lyapunov diagrams," presenting the Lyapunov timescales of the orbits in the u–v velocity plane at fixed space positions, and identified regular and chaotic domains in velocity space as a function of space position with respect to the bar. The Lyapunov exponents were calculated on a Cartesian grid of planar velocities. The fraction of chaotic orbits was demonstrated to obviously increase with the bar strength. However, the diagrams in Fux (2001) can hardly be used to estimate real values of Lyapunov times, because the saturation of the computed values of the Lyapunov exponents was not controlled, while the adopted computation time, corresponding to three Hubble times (equivalently, ∼100 Galactic years), might not be enough for the saturation. In other words, the obtained values characterize not the whole chaotic regions of the phase space, but only rather small vicinities of the initial data. Therefore, the computed values are the "local" values of the Lyapunov exponents. This is testified by the fact that the variation of the Lyapunov exponents in the diagrams of Fux (2001) is smooth, while there must be a sharp distinction between the chaotic regions (with non-zero Lyapunov exponents) and regular regions (with zero Lyapunov exponents) in the divided phase space.

In connection with the problem of estimation of Lyapunov timescales in the solar neighborhood, Quillen (2003) noted that, according to Holman & Murray (1996), for a fully overlapped system, the chaotic zone should have a Lyapunov time ∼2π/ν (where ν is the frequency of perturbation), corresponding to the separatrix pulsation period. In what follows, we shall consider the model set of Quillen (2003) for the stellar dynamics in the solar neighborhood and show that the heuristic estimate ∼2π/ν severely underestimates the real Lyapunov time. Besides, we shall see that it is rather seldom that the considered dynamical systems, modeling the dynamics in the solar neighborhood, can be called "fully overlapped."

Besides obtaining the Lyapunov times, we shall estimate the diffusion times in the chaotic domain of the phase space, in the same model set. This will allow one to judge on the possibility for migration of the Sun from the inner regions of the Milky Way to its current location. Such a possibility, arising due to the overlapping of resonances in the phase space, was advocated and studied in detail by Minchev & Famaey (2010) and Minchev et al. (2011).

Many problems on nonlinear resonances in mechanics and physics are described by the perturbed pendulum Hamiltonian:

$$\begin{equation} H = {{{\cal G} p^2} \over 2} - {\cal F} \cos \varphi + a \cos (\varphi - \tau) + b \cos (\varphi + \tau) \end{equation} \tag{ 1 }$$

(e.g., Shevchenko 2000b). The first two terms in Equation (1) represent the Hamiltonian H₀ of the unperturbed pendulum, where φ is the pendulum angle (the resonance phase angle) and p is the momentum. The periodic perturbations are given by the last two terms, where τ is the phase angle of perturbation: τ = Ωt + τ₀, where Ω is the perturbation frequency and τ₀ is the initial phase of the perturbation. The quantities ${\cal F}$, ${\cal G}$, a, b are parameters.

Generally, Equation (1) describes a triplet (triad) of resonances: there are three trigonometric terms, each corresponding to a particular resonance. In the following sections, the case of a resonance duad is considered (i.e., a or b equals zero).

It is convenient to describe the motion in the vicinity of the separatrices of the Hamiltonian (Equation (1)) by means of the so-called separatrix map (Chirikov 1979). It is a two-dimensional, area-preserving map given by

$$\begin{eqnarray} w_{i+1} &=& w_i - W \sin \tau _i, \nonumber\\ \tau _{i+1}&=& \tau _i + \lambda \ln {32 \over \vert w_{i+1} \vert } \ \ \ (\mbox{mod } 2 \pi). \end{eqnarray} \tag{ 2 }$$

These equations give the classical separatrix map (Chirikov 1979), valid in the symmetric (a = b) case of perturbation. The variable w of the map denotes the relative (with respect to the unperturbed separatrix value) pendulum energy $w \equiv {H_0 \over {\cal F}} - 1$, and τ retains its meaning of the phase angle of perturbation. The parameter λ is the ratio of |Ω|, the absolute value of the perturbation frequency, to $\omega _0 = | {\cal F G} |^{1/2}$, the frequency of the small-amplitude pendulum oscillations. The parameter W in the case of a = b has the form (Shevchenko 1998)

$$\begin{equation} W = \varepsilon \lambda \left(A_2(\lambda) + A_2(-\lambda) \right) = {4 \pi \varepsilon \lambda ^2 \over \sinh {\pi \lambda \over 2}}, \end{equation} \tag{ 3 }$$

where $\varepsilon = a / {\cal F}$, and A₂(λ) is the value of the Melnikov–Arnold integral as defined in Chirikov (1979):

$$\begin{equation} A_2(\lambda) = 4 \pi \lambda {\exp ({{\pi \lambda } / 2}) \over \sinh (\pi \lambda)}. \end{equation} \tag{ 4 }$$

Formula (3) differs from that given in Chirikov (1979) and Lichtenberg & Lieberman (1992) by the term A₂(− λ), which is small for λ ≫ 1. However, its contribution is significant for small values of λ (Shevchenko 1998), i.e., in the case of adiabatic chaos. The kind of chaos (adiabatic or non-adiabatic) in model (1) is identified by the value of λ: if λ < 1/2, it is slow (adiabatic), otherwise it is "fast" (non-adiabatic; Shevchenko 2008a).

One iteration of the separatrix map corresponds to one period of the pendulum rotation or a half-period of its libration. The applicability of the theory of separatrix maps for description of the motion near the separatrices of the perturbed nonlinear resonance in the full range of the relative frequency of perturbation, including its low values, was considered and shown to be legitimate in Shevchenko (2000b).

The separatrix map in the case of asymmetric (a ≠ b) perturbation is different from that in the symmetric (a = b) case, because the energy increments are different for the prograde and the retrograde motions of the model pendulum (Shevchenko 1999). (The motion is called prograde or retrograde if the variation of φ with time is, respectively, positive or negative.) The algorithm, taking this difference into account, constitutes the separatrix algorithmic map (Shevchenko 1999):

$$\begin{eqnarray} & & \mbox{if } w_i < 0 \mbox{ and } W = W^- \mbox{ then } {W}\,{:=}\, W^+, \nonumber \\ & & \mbox{if } w_i < 0 \mbox{ and } W = W^+ \mbox{ then } {W}\,{:=}\, W^-, \nonumber \\ & & w_{i+1} = w_i - W \sin \tau _i, \nonumber \\ & & \tau _{i+1} = \tau _i + \lambda \ln {32 \over \vert w_{i+1} \vert } \ \ \ (\mbox{mod } 2 \pi), \end{eqnarray} \tag{ 5 }$$

where W⁺ and W⁻ denote the values of the W parameter for the prograde and retrograde motions, respectively. In the case of asymmetric perturbation these values are different.

Equations (5) as well can be written in a shorter way (Shevchenko 2000b):

$$\begin{eqnarray} & & \mbox{if } w_i < 0 \mbox{ and } W = W^\pm \mbox{ then } W := W^\mp, \nonumber \\ & & w_{i+1} = w_i - W \sin \tau _i, \nonumber \\ & & \tau _{i+1} = \tau _i + \lambda \ln {32 \over \vert w_{i+1} \vert } \ \ \ (\mbox{mod } 2 \pi). \end{eqnarray} \tag{ 6 }$$

The sign in the upper index of W alternates at each iteration if w_i < 0 (i.e., at librations); the quantity W^± denotes W⁺ or W⁻, while W^∓ denotes a corresponding value of W⁻ or W⁺.

The essence of the separatrix algorithmic map is in taking into account the alternations of the W parameter. It alternates when the direction of motion alternates. This takes place either when rotation changes to libration, or when the motion is librational. Algorithms (5) and (6) do not contain the condition w_i > 0, because the direction of motion does not change when it holds.

In order to find expressions for W⁺ and W⁻, one should integrate the increment of energy per one iteration of the map, following the usual procedure (Chirikov 1979), but making it separately for prograde and retrograde directions of motion. This gives (Shevchenko 1999)

$$\begin{equation} W^+ (\lambda, \eta) = \varepsilon \lambda \left(A_{2}(\lambda) + \eta A_{2}(-\lambda) \right), \end{equation} \tag{ 7 }$$

$$\begin{equation} W^- (\lambda, \eta) = \varepsilon \lambda \left(\eta A_{2}(\lambda) + A_{2}(-\lambda) \right), \end{equation} \tag{ 8 }$$

where $\varepsilon = a / {\cal F}$ and η = b/a. The Melnikov–Arnold integral A₂(λ) is given by Equation (4).

The separatrix map theory can be used for analytical estimation of the maximum Lyapunov exponents (Shevchenko 2000a, 2002, 2007, 2008b). Comparisons of the predictions of this theory with numerical–experimental results can be found in Shevchenko & Kouprianov (2002) and Shevchenko (2000a, 2002, 2007, 2008b, 2009), where it was applied to various problems of celestial mechanics: rotational dynamics of planetary satellites, orbital dynamics of satellite systems, and orbital dynamics of asteroids.

The maximum Lyapunov exponent is defined by the limit

$$\begin{equation} L = \limsup _ {{t \rightarrow \infty } \atop {d (t_0) \rightarrow 0}} {1 \over {t-t_0}} \ln {d (t) \over d (t_0)}, \end{equation} \tag{ 9 }$$

where d(t₀) is the distance (in phase space) between two nearby initial conditions for two trajectories at the initial moment of time t₀, d(t) is the distance between the evolved initial conditions at time t (e.g., Lichtenberg & Lieberman 1992).

According to the general approach used by Shevchenko (2000a, 2002, 2007), the maximum Lyapunov exponent, L, of the motion in the main chaotic layer of system given by Equation (1) is represented as the ratio of the maximum Lyapunov exponent, L_sx, of its separatrix map and the average period T of rotation (or, equivalently, the average half-period of libration) of the resonance phase angle φ inside the layer:L = L_sx/T. In this way, formulas for the Lyapunov time were derived in Shevchenko (2007) for four generic cases of interacting resonances: the fastly chaotic resonance triplet, fastly chaotic resonance doublet, slowly chaotic resonance triplet, and slowly chaotic resonance doublet (called, respectively, "ft," "fd," "st," and "sd" resonance multiplet types).

In the following sections, we shall need formulas for the Lyapunov times for the "fd" and "sd" resonance types only. The formula for the Lyapunov time for the "fd" resonance type is

$$\begin{equation} T_{\mathrm L} = \frac{T_{\mathrm{pert}}}{2 \pi } \cdot \frac{\mu _{\mathrm{libr}} + 1}{ \mu _{\mathrm{libr}} \frac{L_{\mathrm{sx}}(2 \lambda)}{T_{\mathrm{sx}}(2 \lambda, W)} + \frac{L_{\mathrm{sx}}(\lambda)}{T_{\mathrm{sx}}(\lambda, W)}}, \end{equation} \tag{ 10 }$$

where μ_libr ≈ 4 (Shevchenko 2007). The quantity T_pert = 2π/|Ω| is the period of perturbation. The quantities W, L_sx, and T_sx are given by the formulas

$$\begin{equation} W(\varepsilon, \lambda) = \varepsilon \lambda \left(A_2(\lambda) + A_2(-\lambda) \right) = {4 \pi \varepsilon \lambda ^2 \over \sinh {\pi \lambda \over 2}}, \end{equation} \tag{ 11 }$$

$$\begin{equation} L_{\mathrm{sx}}(\lambda) \approx C_h {2 \lambda \over 1 + 2 \lambda }, \end{equation} \tag{ 12 }$$

$$\begin{equation} T_{\mathrm{sx}}(\lambda, W) \approx \lambda \ln {32 e \over \lambda | W |}, \end{equation} \tag{ 13 }$$

where C_h ≈ 0.80 is Chirikov's constant (Shevchenko 2004) and e is the base of natural logarithms.

The formula for the Lyapunov time for the "sd" resonance type is

$$\begin{equation} T_{\mathrm L} \approx {T_{\mathrm{pert}} \over 2 \pi } \ln \left| \frac{32}{\varepsilon \lambda } \sin \left(\frac{\lambda }{2} \ln \frac{8}{ |\varepsilon | \lambda } \right) \right| \end{equation} \tag{ 14 }$$

(Shevchenko 2007).

The parameter λ = |Ω|/ω₀ measures the separation of the perturbing and guiding resonances in the units of one quarter of the width of the guiding resonance, because the resonances separation in frequency space is equal to Ω, while the guiding resonance width is equal to 4ω₀ (Chirikov 1979). Therefore, λ can be regarded as a kind of the resonance overlap parameter. In the asymptotic limit of the adiabatic (slow) case the resonances in the multiplet strongly overlap, while in the asymptotic limit of the "fast" case the resonances are segregated. However, note that the border λ ≈ 1/2 (Shevchenko 2008a) between slow and fast chaos does not coincide with the borderline between the cases of overlapping and non-overlapping of resonances: the latter borderline lies much higher in λ. For example, in the phase space of the standard map, the integer resonances start to overlap (on decreasing λ) at λ ≈ 2π/0.97 ≈ 6.5 (Chirikov 1979).

Quillen (2003, 2009), based on the dynamical model of Contopoulos (1975, 1988), constructed a Hamiltonian of the motion in the solar neighborhood by adding the perturbations due to the bar and spiral arms to the unperturbed Hamiltonian of the motion. Namely, Quillen's model describes interaction of the bar's 2:1 outer Lindblad resonance with the spiral's 2:1 or 4:1 inner Lindblad resonance. The resulting Hamiltonian has the form

$$\begin{equation} H = A[j^2 + \delta j + \beta j^{1/2} \cos \phi + \epsilon j^{1/2} \cos (\phi + \nu t - \gamma)] \end{equation} \tag{ 15 }$$

(Quillen 2003, Equation (23)). Here, j and ϕ are the conjugate momentum and phase variables; A ≈ −5.7; δ, β, , ν, and γ are free unitless parameters. The ranges for numerical values of the parameters were estimated in Quillen (2003) from observational physical and kinematical considerations.

The frequency ν is counted in units of Ω, and, accordingly, time is in units of Ω⁻¹; here Ω is the rotation rate of the epicyclic center. Therefore, one time unit is equal to one Galactic year at a given distance from the center of the Milky Way, divided by 2π.

The first resonant term (that with the coefficient β) in Equation (15) corresponds to the perturbation due to the bar, while the second one (that with the coefficient ) to the perturbation due to the spiral arms. The strength of the second term is much greater than that of the first one (Quillen 2003). Therefore, it is natural to perform a time-dependent shift ϕ = ψ − νt + γ of the origin of the coordinate system. This shift makes the second resonance explicitly the guiding one. The resulting Hamiltonian is

$$\begin{equation} H = A j^2 + (A \delta + \nu) j + A \epsilon j^{1/2} \cos \psi + A \beta j^{1/2} \cos (\psi - \nu t + \gamma). \end{equation} \tag{ 16 }$$

Then we introduce the parameter Δ = Aδ + ν and make a constant shift j = p − Δ/(2A), ψ = φ − γ, reducing Equation (16) to the form

$$\begin{eqnarray} H &=& A p^2 + A \epsilon \left(p - \frac{\Delta }{2A}\right)^{1/2} \cos \varphi \nonumber\\ && + A \beta \left(p - \frac{\Delta }{2A}\right)^{1/2} \cos (\varphi - \nu t), \end{eqnarray} \tag{ 17 }$$

and expand the coefficients in the neighborhood of p = 0, leaving only the lowest order (constant) terms. This gives

$$\begin{equation} H = A p^2 + A \epsilon \left(- \frac{\Delta }{2A}\right)^{1/2} \cos \varphi + A \beta \left(- \frac{\Delta }{2A}\right)^{1/2} \cos (\varphi - \nu t). \end{equation} \tag{ 18 }$$

Thus, we have reduced the perturbed "second fundamental model for resonance" (called so by Henrard & Lemaître 1983), given by the Hamiltonian (Equation (15)), to the perturbed "first fundamental model for resonance," given by the Hamiltonian (Equation (18)), i.e., to the perturbed pendulum model. This has turned out to be possible because the guiding resonance (that emerging due to the spiral arms) is bifurcated in all cases, the librational "crescent," born from the bifurcation, being situated always quite far from the origin of the coordinate system; see the phase space sections in Figures 2–4 in Quillen (2003), Figure 3 in Quillen (2009), and our Figures 1 and 2. The phase space sections in Figures 1 and 2 have been constructed by numerical integration of the equations of motion with the Hamiltonian (Equation (15)) in the same way as described in Quillen (2003); the coordinates in the sections are x = (2j)^1/2cos ϕ and y = (2j)^1/2sin ϕ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 1, the phase space section for model 6 of group C in Quillen (2003) is shown. It corresponds to Figure 4(f) in Quillen (2003), except that we take here a much denser grid of initial data for the trajectories. In the case of Figure 1, the value of γ is equal to π/2, as in Quillen (2003), while in the case of Figure 2 the value of γ is chosen to be equal to zero. In both figures, the central regular island (small crescent) corresponds to the bar's resonance, while the outer regular island (big crescent) corresponds to the spiral's resonance. Setting γ = 0 results mostly in a shift in the angular location of the spiral's resonance in ϕ, that is why the relative angular orientation of the regular islands changes. The angular orientation of the center of the outer island (the fixed point of the spiral's resonance) changes by π/2. Apart from the change of the relative angular orientation of the regular islands, the phase space structure is qualitatively one and the same in the figures. The value of γ does not influence the Lyapunov and diffusion timescales, because this parameter can be removed by a simple canonical transformation, namely, a constant shift in time.

The chaotic domain is pronounced in both figures. It emerges due to interaction of the bar's and spiral's resonances. The estimates of the Lyapunov and diffusion times are made in what follows just for such domains, emerging due to interaction of the resonances, in various model sets.

Equation (18) describes the resonance duad, and, for estimating the Lyapunov time of the motion in the chaotic layer, we can apply formulas (10) and (14) in the cases of fast and slow chaos, respectively.

Comparing the Hamiltonians (Equations (1) and (18)), one has ${\cal F} = - A \epsilon ({-} \frac{\delta }{2} - \frac{\nu }{2 A})^{1/2}$, ${\cal G} = 2 A$, Ω = ν, $a = | A | \beta ({-} \frac{\delta }{2} - \frac{\nu }{2 A})^{1/2}$, b = 0, $\omega _0 = | A | | \varepsilon |^{1/2} ({-}2 \delta - \frac{2\nu }{A})^{1/4}$, $\lambda = \frac{| \nu |}{\omega _0}$, and $\varepsilon = - \frac{\beta }{\epsilon }$. The results of the calculation of the Lyapunov times by formula (10) (when λ > 1/2, i.e., in the majority of cases) and by formula (14) (when λ < 1/2, only in two cases present) are given in Tables 1, 2, and 3 for model groups A, B, and C, respectively. The model groups A, B, and C correspond to the sets of parameter values considered for the construction of the phase space sections in Figures 2–4 in Quillen (2003).

The Hamiltonian (Equation (15)) has five parameters: δ, β, , ν, and γ. All of them are expressed by means of formulas given in Quillen (2003) through original dynamical characteristics of the Galaxy and the solar neighborhood, such as the pattern speeds of the bar and spiral structure, the radius of the solar circle from the Galactic center, the radius at which the bar ends, etc. These five parameters elude straightforward interpretation in terms of the original dynamical characteristics, but they have straightforward meaning in the framework of the dynamical model (15): the parameters β and characterize the strength of the bar's and spiral's resonances, respectively; δ is the "bifurcation control" parameter, as discussed in Quillen (2003); ν measures the separation of the bar's and spiral's resonances in the frequency space; γ is a phase shift, rather unimportant from the dynamical viewpoint, as discussed at the end of Section 3.

The models 1–10 of each model group correspond to the panel of 10 Poincaré sections in a corresponding figure in Quillen (2003): group A corresponds to Figure 2, group B to Figure 3, and group C to Figure 4. The parameters β and are constants inside each model group. Their values are given in the notes to the tables. The parameters δ and ν vary inside each model group.

The obtained values of the Lyapunov time are all (except in one case, corresponding to adiabatic chaos) in the range from ≈40 to ≈80 time units. Inspection of the data given in Tables 1–3 also makes evident that the heuristic estimate ∼2π/ν (Holman & Murray 1996; Quillen 2003) severely underestimates the real Lyapunov time (by an order of magnitude in many cases). Besides, it is rather seldom that the considered dynamical systems, modeling the dynamics in the solar neighborhood, can be called "fully overlapped": the values of λ (playing the role of the resonance overlap parameter; see Section 2) are typically greater than 1/2; so, chaos is non-adiabatic.

Only model 9 (and, marginally, model 8) of model group B can be called adiabatic. The values of the problem parameters , β, ν, and δ in this case are such that the value of the separatrix map parameter λ is less than 1/2. Generally, the adiabaticity of chaos in any model can be checked by calculating λ = |ν|/ω₀, given the values of , β, ν, and δ.

As soon as our time unit is equal to one Galactic year divided by 2π, we see that the obtained typical Lyapunov times are approximately equal to 10 Galactic years. How much is it, in comparison with, e.g., Lyapunov times of the solar system bodies, in comparable time units? The usual Lyapunov times for asteroids in the main belt are 3000–10,000 years or more (Shevchenko 2007), i.e., they are ∼1000 asteroidal orbital periods or more, while the usual Lyapunov times for highly eccentric comets are of the order of one cometary orbital period (Shevchenko 2007). So, the Lyapunov times of the chaotic motion in the considered Galactic model are much less (by at least two orders of magnitude) than the Lyapunov times for the main belt asteroids, but an order of magnitude greater than for the comets, if expressed in adequate time units. In general terms of the loss of predictability of the motion, the Galactic dynamical chaos is rather strong.

From inspection of Tables 1–3, one can deduce that the Lyapunov time depends on the model parameters rather weakly being almost of the same order in all the models. Thus, one can expect that choosing different values for the model parameters, such as relative strengths of the bar and spiral structures (e.g., choosing the spiral structure to be weaker), would not qualitatively change the typical Lyapunov time value. However, note that the change of the model may radically reduce the extent of the chaotic domain. Then chaos is not "global" and so the probability that the Sun belongs to the chaotic domain of the phase space might be small.

It is interesting that, as follows from the above estimates of T_L, the age of the Milky Way measured in its Lyapunov times is about 5–10. This means that now it is already practically impossible to restore exact initial conditions for the stellar dynamics in the solar neighborhood from any observational data.

Let us estimate the characteristic times of chaotic transport (called the diffusion times, when the diffusive approximation is used) in the chaotic domain of the phase space in the same model set of Quillen (2003). Knowledge of these times will allow one to judge on the possibility for migration of the Sun from the inner regions of the Milky Way to its current location. Such a possibility, arising due to overlapping of resonances in the phase space, was advocated and studied in detail by Minchev & Famaey (2010) and Minchev et al. (2011). To estimate the typical diffusion times, we shall base on the approach developed initially by Chirikov & Vecheslavov (1986, 1989) for the purposes of studies in cometary dynamics. Chirikov (1999) employed a similar approach in a study of the separatrix map dynamics (see page 11 and, in particular, formulas (20) and (21) in Chirikov 1999).

First of all, a reservation should be made that it is only very approximate that we can consider the chaotic transport in the problem under study as diffusive. The matter is that the value of λ, as follows from Tables 1–3, in the majority of models is rather low: λ ∼ 1. As explained in Chirikov (1999, p. 12–13), when λ ∼ 1, "... the layer width is reduced down to the size of a single kick.... Hence, the diffusion approximation becomes inapplicable. Instead, the so-called ballistic relaxation comes into play which is much quicker. In other words, a slow diffusive motion ... is replaced now by rapid jumps of a trajectory over the whole layer ..." (A "single kick" is the energy increment per iteration of the separatrix map.) Therefore, all the diffusion rate estimates that we make in this section should be regarded as extrapolation of diffusive description. In reality, they represent upper bounds for the real characteristic times of chaotic transport.

According to Chirikov & Vecheslavov (1986, 1989), the diffusion rate (in the energy variable) in the main chaotic layer in the phase space of the Kepler map¹ approximately equals the mean (over time) squared energy increment per iteration of this map. Analogously, in the case of the ordinary separatrix map (2), the diffusion rate (in the energy variable) in the main chaotic layer in the phase space of the map approximately equals the mean squared energy increment, i.e., 〈W²sin ²τ_i〉. Averaging over the interval 0 ⩽ τ_i < 2π, we find the diffusion rate to be D_map ≈ W²/2.

In the case of the separatrix algorithmic map, the chaotic layer components corresponding to prograde rotations, retrograde rotations, and librations of the phase variable should be considered separately. In the two (prograde and retrograde) rotation cases, the diffusion rate D_map in the energy in the map phase space obviously equals ≈(W⁺)²/2 and ≈(W⁻)²/2 for the prograde and retrograde rotations, respectively. Employing the formulas for Ω, ω₀, a, b (b = 0), λ, and ε, given at the beginning of Section 4, one can calculate the parameters W⁺ and W⁻ of the separatrix algorithmic map (5). If λ > 1/2, the equality b = 0 implies |W⁻| ≪ |W⁺|, while a = 0 implies |W⁻| ≫ |W⁺|. The component of the chaotic layer corresponding to reverse (or direct) rotations does not exist, if W⁻ (or, respectively, W⁺) is equal to zero; its measure is zero. The other circulation component with non-zero measure is described by the ordinary separatrix map (2) with the parameters λ and W = W^± (the non-zero value among W⁺ and W⁻); its extent in w is ≈λ|W| (the half-width of the chaotic layer in the case of fast chaos; see Shevchenko 2008a).

Consider the libration component of the chaotic layer. Then W⁻ and W⁺ alternate (replace each other) at each iteration of the separatrix algorithmic map (5). It is straightforward to show (Shevchenko 2007) that, if W⁻ or W⁺ is equal to zero, the separatrix algorithmic map (5) on the doubled iteration step reduces to the ordinary separatrix map (2) with the doubled value of λ and the value of W equal to a non-zero value of W^±. Taking into account that one iteration of the new map corresponds to two iterations of the original one, the diffusion rate referred to the original map time units is

$$\begin{equation} D_\mathrm{map} \simeq \frac{1}{4}(W^\pm)^2, \end{equation} \tag{ 19 }$$

where W^± is the non-zero value among W⁺ and W⁻.

For the circulation component of the layer, one has D_map ≃ (W⁺)²/2 (if b = 0) or D_map ≃ (W⁻)²/2 (if a = 0).

As soon as the libration motion is reducible to the ordinary separatrix map (2) with the doubled value of λ, the layer's extent in w on the side of librations doubles, it becomes ≈2λ|W|. Note that the parameters λ and W are considered here as independent from each other. Therefore, the chaotic domain corresponding to libration dominates in extent, and for a rough estimate of the diffusion rate over the entire layer it is sufficient to make an estimate for the libration component alone.

In the case of the Hamiltonian (Equation (18)) b = 0 and η = 0, so, W^± = W⁺, where

$$\begin{equation} W^+ = \varepsilon \lambda A_{2}(\lambda) = 4 \pi \varepsilon \lambda ^2 {\exp ({{\pi \lambda } / 2}) \over \sinh (\pi \lambda)} \end{equation} \tag{ 20 }$$

(see Equation (7)), therefore

$$\begin{equation} D_\mathrm{map} \simeq \left(2 \pi \varepsilon \lambda ^2 {\exp ({{\pi \lambda } / 2}) \over \sinh (\pi \lambda)} \right)^2 \end{equation} \tag{ 21 }$$

in the libration case.

To obtain the diffusion rate referred to real time units, it is necessary to transform the map time units into the real ones. This is performed by dividing the diffusion rate referred to map time units by the mean period of phase rotations (half-periods of librations) inside the chaotic layer, because this mean period is nothing but the average time interval corresponding to one map iteration. Consequently, the diffusion rate referred to real time units is D = |Ω|D_map/T_sx, where T_sx is given by formula (13).

We define the characteristic diffusion time across the chaotic layer to be equal to the inverse of the diffusion rate. Therefore, it is just the time needed for the diffusing particle to cover the relative energy interval equal to one. Note that the maximum possible deviation in the relative energy w from zero in the libration case is equal to −2 (Chirikov 1979); therefore, the defined diffusion time gives an appropriate time estimate for the global mixing inside the chaotic layer, of course, if the chaotic layer is broad enough.

In our Hamiltonian (Equation (18)) b = 0, so, W^± = W⁺, and one gets for the diffusion time

$$\begin{equation} T_{\mathrm{d}} = \frac{1}{D} = \frac{T_{\mathrm{sx}}(\lambda, W^+)}{| \Omega | D_\mathrm{map}} \simeq \frac{4 T_{\mathrm{sx}}(\lambda, W^+)}{| \Omega | (W^+)^2}, \end{equation} \tag{ 22 }$$

where

$$\begin{equation} T_{\mathrm{sx}}(\lambda, W^+) \approx \lambda \ln {32 e \over \lambda | W^+ |}, \end{equation} \tag{ 23 }$$

(see Equation (13)), e is the base of natural logarithms, and W⁺ is given by formula (20). Finally one has

$$\begin{equation} T_{\mathrm{d}} \simeq \frac{4 \lambda }{| \Omega | (W^+)^2} \ln {32 e \over \lambda | W^+ |}. \end{equation} \tag{ 24 }$$

In the case of the prograde rotation component of the chaotic layer, the diffusion rate D_map ≃ (W⁺)²/2; therefore, the diffusion time is two times less than that given by formula (24).

The results of the calculation of the diffusion times T_d by formula (24) are given in Tables 4, 5, and 6 for model groups A, B, and C, respectively.

Inspection of the data in Tables 4–6 makes it evident that the diffusion times vary considerably inside the model groups: by as much as four orders of magnitude in group C. Taking into account that our time unit is equal to one Galactic year divided by 2π, i.e., our time unit ≈32 Myr (1 Myr =10⁶ year; 1 Galactic year ≈200 Myr), it is clear that the obtained diffusion times in most of the models with ordinal numbers up to three are greater than 10 Gyr (1 gigayear ≡ 1 Gyr =10⁹ year), making large-scale radial chaotic migration improbable. Such "junior" models all have positive values of δ. However, models 7–9 in group B also have large values of T_d.

Minchev & Famaey (2010) found in detailed numerical experiments that, due to the resonance overlap of the bar and spiral structure, the Galactic disk mixes in about 3 Gyr. From our Tables 4–6 it is obvious (taking into account that our time unit ≈32 Myr) that most of the models with negative values of δ provide chaotic mixing effective enough to be consistent with the numerical-experimental results of Minchev & Famaey (2010). Almost all models with large ordinal numbers, except models 7–9 in group B, permit such migration. Generally, as follows from data in Tables 4–6, the diffusion rate strongly depends on the radial position in the Galaxy.

Note that the Hamiltonian model (Equation (15)) assumes that the guiding radius in the unperturbed dynamics is fixed. If migration is present, the guiding radius varies; so, a more refined model should be used to account for this variation self-consistently. One should also mention that chaotic transport can be ubiquitous in the Galaxy, though its dynamical origin is presumably different from that considered here. Quillen et al. (2010) showed that the diffusion might occur in many regions due to interaction of multiple waves, and so overlapping of resonances should be common; besides, if the bar slows down (e.g., Weinberg & Katz 2007 and references therein) resonances can sweep through the Galaxy (Quillen et al. 2010). Of course, these processes cannot be described in the specific dynamical model considered here, but this model already provides an insight into the possible effectiveness of chaotic transport in the Galaxy.

The most important conclusion of this section is that models that permit large-scale radial chaotic migration of the Sun (from the inner regions of the Milky Way to its current location) do exist. This confirms the dynamical possibility of the migration concept advocated by Minchev & Famaey (2010) and Minchev et al. (2011). What is more, due to the possibility of ballistic flights mentioned in the beginning of this section, the chaotic mixing might be far more effective and quicker than in the case of normal diffusion. Obviously, the effect of ballistic relaxation should be explored in detail in the future.

We have considered how the Lyapunov and diffusion timescales of the stellar dynamics in the solar neighborhood can be estimated. We have used Quillen's (2003) model to describe interaction of the "spiral" and "bar" nonlinear resonances in the phase space of the motion. A method of analytical estimation of the maximum Lyapunov exponents of the orbital motion has been applied to the solar neighborhood dynamics. The analytical treatment has been performed within a framework of the separatrix map theory (Shevchenko 2000a, 2002, 2007), describing the motion near the separatrices of a perturbed nonlinear resonance.

The Lyapunov times turn out to be basically in the range from 6 to 13 Galactic years. In comparison with the Lyapunov times of the solar system bodies (made in adequate time units), the Galactic dynamical chaos is rather strong in general terms of the loss of predictability of the motion. An interesting inference is that, as soon as the age of the Milky Way measured in its Lyapunov times is about 5–10, it is already practically impossible to restore exact initial conditions for the stellar dynamics in the solar neighborhood from any observational data.

We have also estimated the diffusion times, based on the approach developed initially by Chirikov & Vecheslavov (1986, 1989) for the purposes of studies in cometary dynamics. We have found that, in a number of models, the diffusion times turn out to be small enough to permit radial chaotic migration of the Sun from the inner regions of the Milky Way to its current location. In other words, dynamically adequate models that permit large-scale radial chaotic migration do exist. This confirms the dynamical possibility of the migration concept advocated by Minchev & Famaey (2010) and Minchev et al. (2011). Due to the possibility of ballistic flights inside the chaotic layer, arising because λ ∼ 1, the chaotic mixing might be even far more effective and quicker than in the case of normal diffusion.

We have shown that only in a narrow range of possible values of the problem parameters , β, ν, and δ the Galactic chaos is adiabatic because the values of the separatrix map parameter λ, playing the role of the resonance overlap parameter, are typically greater than 1/2; in other words, adiabatic chaos (λ < 1/2) seems to be not characteristic for the dynamics in the solar neighborhood.

The author is thankful to Alice Quillen for advice and comments. It is also a pleasure to thank the referee for helpful remarks.