Issue 
A&A
Volume 642, October 2020



Article Number  A207  
Number of page(s)  10  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/202038959  
Published online  23 October 2020 
Diffusion of radial action in a galactic disc
LUPM, Univ Montpellier, CNRS, Montpellier, France
email: herve.wozniak@umontpellier.fr
Received:
17
July
2020
Accepted:
1
September
2020
Context. The stellar migration of the galactic disc stars has been invoked to explain the dispersion of stellar metallicity observed in the solar neighbourhood.
Aims. We seek to identify the dynamical mechanisms underlying stellar migration in an isolated galaxy disc under the influence of a bar. Our approach is to analyse the diffusion of dynamical quantities.
Methods. We extend our previous work by exploring Chirikov’s diffusion rate (and derived timescale) of the radial action J_{R} in an idealised Nbody simulation of an isolated disc galaxy. We limit our study to the evolution of the disc region well after the formation of the bar, in a regime of adiabatic evolution.
Results. The J_{R} diffusion timescale T_{D}(J_{R}) is less than 3 Gyr for roughly half the galaxy mass. It is always much shorter than the angular momentum diffusion timescale T_{D}(L_{z}) outside the stellar bar. In the disc, ⟨T_{D}(J_{R})⟩ ∼ 1 Gyr. All nonaxisymmetric morphological structures that are characteristic of resonances and waves in the disc are associated to particles with T_{D}(J_{R}) < 3 Gyr and T_{D}(L_{z}) > 10 Gyr. Short T_{D}(J_{R}) can be explained by the gradual decircularisation of initially circular orbits (J_{R} = 0) under the effect of intermittent. Inner Linblad resonance scattering by wave trains propagating in the disc, well beyond the outer Lindblad resonance of the bar (OLR). This leads to a moderate secular heating of the disc beyond the bar’s OLR for 7 Gyr, which is comparable to solar neighbourhood observations. The complex multiwave structure, mixing permanent and intermittent modes, allows for multiple resonance overlaps.
Key words: Galaxy: disk / Galaxy: evolution / Galaxy: kinematics and dynamics / Galaxy: structure
© H. Wozniak 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
In Wozniak (2020), we used, for the first time, the formulation of the diffusion rates introduced by Chirikov (1979), applied to both specific energy E and angular momentum L_{z} in selfconsistent Nbody experiments of isolated galactic discs. Using the same definition, we extend our previous work by now focusing on the radial action J_{R}.
The radial action has been introduced in the galactic problem by Freeman (1966) and Kalnajs (1971). In the context of the epicycle approximation, J_{R} = E_{R}/κ, where κ is the epicycle radial frequency and E_{R} is the specific radial kinetic energy. This approximation is only valid if J_{R} ≪ L_{z} or for nearcircular orbits (Kalnajs 1971).
Additionally, J_{R} has been considered as a thermometer for measuring the disc heating in the radial direction. Whether or not the stellar disc heats during the secular exchange of angular momentum is an open issue. In particular, it is a question of being able to determine the relative importance of two phenomena, churning and blurring, for stellar migration (Schönrich & Binney 2009; Halle et al. 2015). In the context of epicycle approximation, the effect of blurring is to increase the amplitude of epicycle motion around a fixed guidingcentre radius (thus generating no net radial migration in principle), whereas that of churning is to move the guidingcentre radius inwards or outwards.
According to Sellwood & Binney (2002), any variation of J_{R} in a stellar disc excited by a constant nonaxisymmetric perturbation of pattern speed Ω_{p}, can be related to ΔL_{z} as:
where Ω is the orbit and particle angular frequency. At corotation, Ω = Ω_{p} so that there is no J_{R} variation; whereas at Lindblad resonances with an marmed perturbation,
This applies to any kind of disturbance, such as a bar or spiral(s). It can be generalised to higher order resonances, so that ΔJ_{R}/ΔL_{z} = ∓l/m for any lthorder resonance, where l = 0 is the corotation.
To complete the picture, it is necessary to add a last relation: J_{R} is proportional to the square of the epicycle amplitude (Binney & Tremaine 2008). Therefore, variations in L_{z} that do not produce significant variations in J_{R} can be interpreted as the mark of the churning mechanism. Corotation scattering is thus expected to play a major role in this mechanism (Sellwood & Binney 2002; Roškar et al. 2012). On the contrary, variations in J_{R} that are not correlated with variations in L_{z} via Eq. (1) are typical of blurring. When both L_{z} and J_{R} vary, the situation is much more complex: Churning and blurring are mixed in timedependent proportions (Halle et al. 2015), and scattering can be nonresonant (e.g. with giant molecular clouds) or made by Lindblad resonances (e.g. the inner Lindblad Resonance (ILR) Sellwood 2012). Moreover, the overlap of resonances between the bar and the spirals (Minchev & Famaey 2010; Minchev et al. 2011) is expected to increase the radial energy so that stars that were originally on nearly circular orbits (J_{R} ≈ 0) can move to more eccentric orbits. Additional complexity comes from the presence of a strong bar. Indeed, the congestion of n/1 resonances near the corotation generates a high degree of stochasticity. The diffusion capacity of a bar’s corotation could therefore be more important than that of a spiral.
The introduction of the Chirikov diffusion rate makes it possible to quantify the impact of the accumulation of small fluctuations in energy, angular momentum, and radial action over time. This quantity takes all variations weighted by all timescales into account. In addition to other methods, it helps to identify the predominant dynamical mechanism(s), through the associated timescales. On the other hand, it does not contain any spatial information. But we have shown (Wozniak 2020) how to get around this obstacle.
The paper is organised as follows. After reintroducing some basic notations and concepts in galactic dynamics and detailing how J_{R} is estimated in our simulations (Sect. 2), we present the results on the Chirikov diffusion timescale of J_{R} in Sect. 3. Section 4 contains a discussion on the stellar distribution function as a function of L_{z} and J_{R}, and its time evolution. Section 5 describes how circular orbits of the outer disc absorb angular momentum and gain eccentricity. A wave analysis is reported in Sect. 6. Our results are discussed in Sect. 7, whereas Sect. 8 summarises our conclusions.
2. Computing J_{R} diffusion
2.1. Dynamical concepts
The Chirikov diffusion rate of J_{R} for individual particles can be defined as follows:
Although the original definition only deals with E (Chirikov 1979, Eq. (4.6)), we have extended the definition to compute D_{n}(L_{z}) and, in fact, any measurable quantity evolving during a timedependent simulation (Wozniak 2020). In Eq. (3) for n = 2, is the value of radial action averaged over a period of Δt_{2} = 10^{2} (in the unit of time of the system). We can also estimate a diffusion timescale T_{D} by renormalising D_{2} by , where the timeaverage is now computed over the longest timescale (e.g. the experiment length):
Several methods exist to estimate actions numerically (see Sanders & Binney 2016, for a comprehensive review). Vasiliev (2019) provides a stateoftheart software package for building dynamical models based on an actionangle description, including positionvelocity and actionangle conversion routines. However, the accurate and fast estimation of actionangle variables remains an open problem, especially for nonaxisymmetric dynamic systems in fast rotation, such as barred galaxies. The approach by Sanders & Binney (2014) requires sampling all orbits in such a way that all periods of each orbit can be represented. In the case of large Nbody simulations, this objective is still extremely difficult, if not impossible, to achieve for strictly computational reasons. Therefore, we have decided to resort to a more qualitative than quantitative approach, estimating the radial action as if the orbits of the dynamical system could be described by the epicycle approximation at each stage of the evolution of a galaxy. This is certainly the best solution, especially as many averages are used in the calculation of the Chirikov diffusion rate.
For a nearly circular orbit in the symmetry plane of an axisymmetric potential in the form Φ(R, z) = Φ_{R}(R)+Φ_{z}(z), J_{R} reduces to
where . In the limit of very small radial excursions around the guiding centre (X = (R_{max} − R_{min})→0), J_{R} reduces even further to , where κ, E_{R}, and X are the epicycle radial frequency, the radial kinetic energy, and the epicycle amplitude, respectively (Binney & Tremaine 2008). The integral must be taken around a full radial period, a quantity largely unknown, and not accessible, in an Nbody simulation. Since we need in Eq. (3) rather than the exact instantaneous value, which is difficult to compute, the following averaging procedure allows us to obtain such an approximate expression:
It is convenient to timeaverage here over Δt_{2} = 10^{2} time units (i.e. 105.4 Myr). Other time windows are used throughout this article and will be clarified in due course. As we use Eq. (6) to study the diffusion of J_{R}, we never have access to an instantaneous value for the radial action.
No assumption that leads to the approximation J_{R} = E_{R}/κ is exact when a galaxy is barred. Indeed, in the bar, the gravitational potential is not vertically separable at any point; additionally, κ, which was derived from the axisymmetric part of the potential, is not correct since isopotentials are no longer circular. Clearly, the epicycle approximation is only correct away from the influence of the bar. Therefore, in order to test the robustness of our conclusions, we also used the generalised κ formulation due to Pfenniger (1990) and we tried to estimate as the ratio of averages . The results quantitatively differ but this does not drastically change our conclusions. In the following, we mention the results most strongly impacted by our choice for estimating .
2.2. Numerical implementation
In the following, we use the RunC simulation of Wozniak (2020) as a reference run. The results are qualitatively identical for the other simulations. Initial stellar populations are set up to reproduce an idealised disc galaxy. Scale lengths and scale heights have been chosen so as to shape an initial axisymmetric disc galaxy with a small but significant bulge. The total mass is M_{tot} = 2 × 10^{11} M_{⊙}. All 4 × 10^{7} particles have the same individual mass so that the plots expressed in a relative particle number or mass fraction are equivalent. The initial disc size is 40 kpc for a scalelength of ≈4 kpc. Initial velocity dispersions are anisotropic solutions of Jeans equations, with σ_{R} = σ_{z} and , where σ_{R}, σ_{θ}, and σ_{z} are three components of the velocity dispersion along the radial, azimuthal, and vertical directions, respectively.
We only study the phase after the formation of the bar (t > 3.16 Gyr, until the end of the simulation at t = 10.54 Gyr), in a regime that can be considered as adiabatic. Therefore, the bar formation mechanism has no direct influence on our results. This obviously does not mean that the bar has no influence on the rest of the disc at all.
Unless otherwise stated, we have excluded, from the analyses, particles escaping the 3D logpolar grid (R > 100 kpc or z > 7.8 kpc) as soon as they were out by even a single timestep between t = 3.16 and 10.54 Gyr. During the Nbody computation, they were tracked by ballistic approximation until they reentered the grid in order to ensure the best possible conservation of momenta. These particles amount to 0.11 M_{tot} at the end of RunC. The other particles are named “neverescaped” in the following.
3. Diffusion timescale (T_{D}(J_{R})) results
The normalised relative frequency distribution (akin to a probability density function) of T_{D}(J_{R}) is plotted in Fig. 1. For the sake of clarity, we have restricted this figure to the range 10^{−2} − 90 Gyr. Distributions for T_{D}(J_{R}) and T_{D}(E) are overplotted for comparison purposes (Wozniak 2020). The distribution shapes are different, especially for T_{D} ≲ 0.9 Gyr. On the contrary to E and L_{z}, there are few particles with T_{D}(J_{R}) < 100 Myr. The J_{R} diffusion timescale is never very short. The maximum of T_{D}(J_{R}) is located at ≈0.9 Gyr, followed by a plateau up to ≈3 Gyr. Additionally, T_{D}(J_{R}) = 3 Gyr is also the median in mass or number of particles. On both sides, the frequency of T_{D}(J_{R}) decreases sharply. In other words, the characteristic diffusion timescale of J_{R} is remarkably of the same order of magnitude as dynamical timescales in the galactic disc. This result contrasts with what was obtained for L_{z} and E, for which T_{D} distributions are dominated by short timescales, outside the local maxima at T_{D}(E) ∼ 10 and T_{D}(L_{z}) ∼ 1 Gyr. On the other hand, the frequencies are quite similar beyond 10 Gyr.
Fig. 1. Normalised particle relative frequency per Gyr as a function of T_{D}(J_{R}) (black histogram). A binsize of 0.01 Gyr has been used. One particle thus represents a probability density of 2.5 × 10^{−6} Gyr^{−1}. Distributions for T_{D}(E) (red) and T_{D}(L_{z}) (blue) are plotted for comparison purposes. 
Figure 2 shows T_{D}(J_{R}) averaged over sets of particles (designated by ⟨T_{D}(J_{R})⟩) sampled by ranges, where is now timeaveraged over ≈7.4 Gyr (i.e. from t = 3.16 to 10.54 Gyr). It must not be confused with an instantaneous L_{z} that is not conserved. The range in occupied by the bar Lindblad resonances during the evolution is delimited by the innermost position of the bar’s ultra harmonic resonance (UHR_{B} at t = 3.16 Gyr) and the outermost bar’s outer Lindblad resonance (OLR_{B} at t = 10.54 Gyr), whereas the range covered by the corotation (CR_{B}) is approximately represented by the shaded area. Again, for the sake of comparison, ⟨T_{D}(E)⟩ and ⟨T_{D}(L_{z})⟩ obtained in Wozniak (2020) are overplotted.
Fig. 2. ⟨T_{D}(J_{R})⟩ (black line), ⟨T_{D}(E)⟩ (red line), and ⟨T_{D}(L_{z})⟩ (blue line) as a function of . The shaded area delimits the region occupied by bar corotation (CR_{B}) over ≈7.4 Gyr. Vertical lines show L_{z} of circular orbits at the innermost bar’s UHR (UHR_{B} at t = 3.16 Gyr, dotdashed) and the positions reached by the outermost bar’s OLR (OLR_{B} at t = 10.54 Gyr, longdashed). 
As shown in Wozniak (2020), ⟨T_{D}(E)⟩ decreases overall from the centre () to the outermost regions, with a strong depression in the area occupied by the bar’s Lindblad resonances, while ⟨T_{D}(L_{z})⟩ increases monotonically to values that can be considered as slowdiffusion. This is not the case for ⟨T_{D}(J_{R})⟩, which decreases regularly from values similar to ⟨T_{D}(L_{z})⟩ in the centre (∼10 Gyr) to values around 1 Gyr in the disc. Beyond UHR_{B}, ⟨T_{D}(J_{R})⟩ shows a plateau between 0.6 and ≈1 Gyr, which explains the bump around 0.9 Gyr in Fig. 1. Apart from the very centre, ⟨T_{D}(J_{R})⟩ is always smaller than ⟨T_{D}(L_{z})⟩.
Particles that are retrograde, on average (0.1 M_{tot}, typical of a barred galaxy, and an identical proportion among neverescaped particles), have shorter diffusion times for L_{z} and J_{R} than for E. The values remain compatible with the timescales in the disc for L_{z} and J_{R}. The scattering of retrograde particles is a subject in and of itself because the relative velocity of these particles with respect to any prograde waves is very large, so that their interaction can only take place over a very short time. We do not address this point in this article.
The particle or mass density distribution in the T_{D}(L_{z})−T_{D}(J_{R}) plane can be studied at different times during the simulation. Nevertheless, for this first exploratory study, we found it more interesting to focus on long times. Therefore, this distribution is analysed for the final snapshot (t = 10.54 Gyr) and plotted in Fig. 3. It shows many structures. They correspond to dynamical sets of particles with a similar behaviour. Projected on the T_{D}(J_{R}) axis, we can recover the peak and plateau identified in Fig. 1. The most remarkable structure extends along the bisector T_{D}(J_{R}) = T_{D}(L_{z}). As is seen in the following, it clearly belongs to the bar. Substructures can be linked to families of orbits, but that work is beyond the scope of this article. By integrating along the T_{D}(L_{z}) axis, this large structure explains both the plateau and the quasilinear decay observed in Fig. 1 when T_{D}(J_{R}) > 3 Gyr. The other noticeable structure with the same density and long T_{D}(L_{z}) (> 10 Gyr) is related to the disc. It contributes to a large extent to the peak at T_{D}(L_{z})) = 0.9 Gyr, although the two structures are blended by integrating along the T_{D}(L_{z}) axis.
Fig. 3. Normalised particle relative frequency (Gyr^{−2}) in the T_{D}(L_{z})−T_{D}(J_{R}) plane in log scale for neverescaped particles only. Contours are spaced by 0.15 dex. The white line divides the domain into two sets (named Subset A and Subset B), which are used for Figs. 4–6 (see text for details). 
It can be expected that these dominant structures correspond to marked morphological counterparts in physical space. It is quite difficult to isolate each of the substructures to determine to which morphological element, orbit family, or wave type they correspond. For the sake of simplicity, we have decided to isolate only one particular subset (named B in Fig. 3 and Subset B in the text). After a few trials and by focusing on the disc properties, Subset B has been defined as T_{D}(L_{z}) > 10 Gyr (as in Wozniak 2020) and T_{D}(J_{R}) < 3 Gyr (roughly the end of the plateau in Fig. 1). Other thresholds have been tested, but these ones roughly isolate particles with “long” T_{D}(L_{z}) and “short” T_{D}(J_{R}). Mass fractions are calculated with respect to the total mass of RunC: 0.65 M_{tot} for Subset A, 0.24 M_{tot} for Subset B, and the rest being excluded particles.
Figure 4 shows the mass surface density projected in the x–y plane for the two particle sets defined above. Subset A, which is the most massive (0.65 M_{tot}), mainly contains particles with a highly symmetrical mass distribution, especially the stellar bar. Beyond the corotation, the distribution in the disc only shows very small deviations from axisymmetry. Inside UHR_{B}, the morphology is ellipticallike, which is the signature of strong bars (Skokos et al. 2002; MichelDansac & Wozniak 2006). The case of Subset B is significantly different. When T_{D}(L_{z}) > 10 Gyr and T_{D}(J_{R}) < 3 Gyr, the mass distribution (0.24 M_{tot}) shows many morphological substructures associated to the presence of waves and resonances, such as spiral arms and rings. Several structures also exist inside the bar (such as ansea, Buta 2019), but they extend well beyond OLR_{B}. As with Subset A, there appear to be several subpopulations, which have not been separated because of our approximate criterion for defining Subset B. However, since Subset B corresponds to a physical region involved in stellar migration, we analyse it in greater detail.
Fig. 4. Projected mass surface density (in M_{⊙} pc^{−2}) at t = 10.54 Gyr inside ±50 kpc. Particle sets are defined in Fig. 3 and in the main text. The log colourscale is common to both figures. Black isodensities are spaced by 0.4 dex for Subset A (left panel) and 0.1 for Subset B (right panel). The spatial scale is in kpc. 
4. Distribution functions (DFs)
4.1. DF in L_{z} and J_{R}
In Fig. 5, the DFs for “neverescaped” particles and the two particular selections (Subset A and Subset B) are plotted as a function of and timeaveraged between t = 3.16 and 10.54 Gyr. The “averaged” DFs include the signatures of all temporal events. This can be compared to Sellwood (2012), for instance. However, we should not expect to find exactly the same results since our initial stellar disc is not such a stable disc. The DF of RunC is obviously much more structured than in Sellwood’s experiments. In particular, both L_{z} and J_{R} bear the stigma of the bar and its formation (occurring during the first Gyr). Both integrals have been largely redistributed, especially in Subset A.
Fig. 5. Mass density per (kpc km s^{−1})^{2} in the plane in log scale for neverescaped particles only. This is a timeaveraged 2D representation of DF over 7.4 Gyr. Contours are spaced by 0.2 dex. White isodensity contours represent Subset A; black contours correspond to Subset B. 
Subset B is identifiable through several substructures. A density peak is present for (circular orbits) and (beyond OLR_{B} at t = 10.54 Gyr). This region is bordered by two almost vertical tails (centred at 4000 and 4500 kpc km s^{−1}) which, similarly to Sellwood (2012), may have been formed by resonant scattering at a Lindblad resonance. That assumption has to be challenged, but we can already claim that it cannot be a resonance with the bar in this region of the disc.
Between 2300 and 3600 kpc km s^{−1}, the large vertical tail with high values seems to include a significant fraction of the socalled hot population (Sparke & Sellwood 1987). These orbits spend most of their time outside the bar and sometimes enter inside the bar from the L_{1, 2} Lagrangian points. The tail width is likely to be related to corotation shifts over time. We note that can reach very high values there. The mass distribution is very sensitive to how is calculated. Using the ratio as an estimator of , the mass would have been concentrated around a maximum density located at , which is far from being representative of typical trajectories in this region. Indeed, particles of the “hot” population explore large portions of the disc, resulting in large variations in κ for each of them, as well as radial kinetic energy E_{R}. Therefore, this invalidates the epicycle approximation for hot population orbits.
For kpc km s^{−1}, Subset B exhibits a ridge that bridges Subset B and Subset A. This part is linked to ansae identified in Fig. 4 and likely associated to UHR_{B}. A component separation based on T_{D} alone is not identical to a separation based on morphological criteria. So it is not surprising that a fraction of Subset B belongs to what we identify as the stellar bar. A slightly smarter component separation algorithm could probably separate this contribution, which mainly seems to be due to the bar from the rest of the disc.
4.2. DF time evolution
Onedimensional DF(J_{R}) or DF(L_{z}) can be recovered by the projection of the density map on the axes. We note that DF(L_{z}) is then similar to typical profiles obtained by a wealth of 3D Nbody simulations, for example, those of Zang & Hohl (1978), Sparke & Sellwood (1987), or Pfenniger & Friedli (1991). The typical DF profile of barred galaxies has been explained by a superposition of various families of orbits (Sparke & Sellwood 1987; Wozniak & Pfenniger 1997), including the abovementioned hot population.
In order to identify the time evolution of some identifiable structures in DF(L_{z}), we have calculated over only 1 Gyr at three different moments of the simulation. The intervals were centred at , 6.8 and 10.0 Gyr. The criteria used to define the selection of Subset A, Subset B, and “neverescaped” particles remain identical as at t = 10.54 Gyr. For the sake of comparison, we have normalised all DF(L_{z}) to the maximum DF(L_{z} = 0) at Gyr.
In Fig. 6 (bottom panels), the contribution of the two subsets is clearly separated. The whole stellar bar forms the peak of DF and contributes mainly to Subset A. The disc, both in a large axisymmetric fraction and the whole resonant structures, forms Subset B. This region contains the hot population bump ( kpc km s^{−1}). Unsurprisingly, particles that spend some time outside the simulation grid, preferentially belong to the disc. They make a significant contribution to the hot population bump. The secondary smaller bump for kpc km s^{−1} has been previously identified as a ridge overlapping with Subset A. It is linked to bar structures (cf. Fig. 4).
Fig. 6. Distribution functions (DFs) for all particles (green triple dotteddashed lines), “neverescaped” particles (dotted lines), and the two subpopulations defined in Fig. 4 (red and black lines) as a function of (top) and (bottom). We note that and have been averaged over ≈1 Gyr, starting at t = 3.2 (left), 6.3 (middle), and 9.5 Gyr (right). DFs have been normalised to DF and DF for the interval t = 3.2−4.2 Gyr, respectively. 
Other evidence is that the number of particles with L_{z} = 0 decreases over time, which means that the number of particles close to the centre or on radial orbits decreases. As the particle number of Subset A is constant by definition, this also means that the redistribution of L_{z} within the bar continues. For Subset B, resonant structures are visible as small oscillations along the hot population bump.
The projected distribution (Fig. 6 top panels) blends all the components discussed so far. Most particles that temporarily escaped have kpc km s^{−1}. The mode is mainly due to Subset A. Its position shifts from kpc km s^{−1} at Gyr to ≈90 kpc km s^{−1} at Gyr. In this time frame, the distribution width increases by ≈23% for Subset A, while a substructure appears in the Subset B distribution, leading to a small bump around kpc km s^{−1}. This behaviour is symptomatic of moderate but regular radial heating of the disc.
5. Evolution of circular orbits (J_{R} = 0)
Let us take a closer look at what may be one of the causes of radial heating in RunC. The spread of DF increases over time as the bin depopulates. Since DF has been normalised by DF at Gyr, Fig. 6 clearly shows that the number of particles with decreases with time, even long after the bar has been formed. This strongly points to an increase in epicycle amplitude as . The evolution of these particles on nearcircular orbits thus deserves a particular analysis since the variation of J_{R} is discriminating with respect to the blurring and churning issue.
In order to take the numerical uncertainties inherent to this type of simulation into account, let us define hereafter circular orbits as kpc km s^{−1}. If we select only the particles with , the evolution of DF and DF can be extracted at time and Gyr (Fig. 7). This circular orbit population contains about 0.039 M_{tot} (i.e. 7.7 × 10^{9} M_{⊙}). Almost all particles are located in the tail in the outermost part of the disc (), and they are mostly well beyond OLR_{B} (located at 20.8 kpc, i.e. L_{z} ≈ 3790 kpc km s^{−1}, at t = 10.54 Gyr). Very few of these particles belong to the hot population bump, which is centred on .
Fig. 7. Particle number as a function of (top) and (bottom), for particles selected as kpc km s^{−1} (black lines). The same particle selection is drawn at t = 6.3−7.3 (in green) and t = 9.5−10.5 Gyr (in red). The contribution of particles with strictly increasing J_{R} to DF() is shown by a dashed line. 
Figure 7 shows that increases significantly for ≈57.4% of particles initially on nearcircular orbits. It leads to an increase in σ_{R} of the order of ≈10 km s^{−1} beyond OLR_{B}. This is comparable with observations of the solar neighbourhood (Soubiran et al. 2008; Mackereth et al. 2019, for instance); although, RunC does not reproduce the properties of Milky Way. The corresponding impact on DF is perfectly identifiable by spikes in the bump of this population. This is clearly the signature of a coherent mode of L_{z} exchange, coupled with the increase in J_{R}, presumably in the form of waves. Well beyond OLR_{B}, it is clearly the impact of propagating density waves with Ω_{p} < Ω_{B} that is illustrated here. According to Eq. (1), if the particle scattering is resonant with a wave, it cannot be at corotation since ΔJ_{R} ≠ 0.
The mass surface density at t = 3 Gyr of particles selected at Gyr as having is plotted in Fig. 8. As expected, this set of particles has a very axisymmetric distribution, mostly beyond OLR_{B}, with the notable exception of those inside the bar, which are concentrated around the Lagrange L_{4, 5} points (the bar is oriented at about 45°). At t = 10.54 Gyr, the distribution of the same particle set shows many wavelettype structures. Since J_{R} of more than half of the particles have increased significantly, then it is normal that the initial axisymmetry has been broken. In anticipation of Sect. 6, we have plotted some resonances identified from three patterns detected in the disc: the bar, an intermediate spiral structure, and a set of external waves. The structures that appear progressively between t = 3 and t = 10.54 Gyr are mainly beyond OLR_{B} and do not exceed the corotation of the outermost waves.
Fig. 8. Projected mass surface density for particles selected as kpc km s^{−1} at t = 3 Gyr (left panel) inside ±50 kpc. Their distribution at t = 10.54 Gyr (right panel). The log colourscale is common to both figures (in M_{⊙} pc^{−2}). Black isodensities are spaced by 0.2 dex. Circles show the position for the ILR (dotdashed line), the UHR and outer m = + 4 resonance (shortdashed), the CR (full line), and the OLR (longdashed). In red, the following resonances for the stellar bar are shown: UHR_{B}, CR_{B}, and OLR_{B}, as defined in Sect. 6. Using the same linestyle, the white circles represent the intermediate spiral: ILR_{iS} (close to UHR_{B}), CR_{iS} (close to OLR_{B}), outer m = +4 resonance, and OLR_{iS}. The black circles represent the outer wave ILR_{oW} (close to intermediate spiral outer m = +4 resonance), UHR_{oW}, and CR_{oW}. 
6. Wave analysis
6.1. Fourier spectrograms
In order to give a coherent overview of the dynamical mechanisms at work in this simulation, we have been looking for waves in the disc that could be associated with variations in angular momentum and radial action. Figure 9 shows the classical m = 2 and m = 4 spectrograms cumulated in time windows t = 3.16−4.24, 6.32−8.48, and 8.38−10.54 Gyr, respectively. The first window (1075 Myr wide) is shorter than the others, which are twice as large, since the slowdown rate of Ω_{B} is higher when the bar is young.
Fig. 9. Top row: m = 2 power spectra in log scale as a function of radius for RunC. The time windows are 3.16−4.24 Gyr (left panel), 6.32−8.48 Gyr (middle), and 8.38−10.54 Gyr (right). The vertical scales give values of Ω in Myr^{−1} (left) and in km s^{−1} kpc^{−1} (right). The radial scale (in kpc) is converted in L_{z} using the circular rotation curve. The mean curves Ω ± κ/2 are drawn as black shortdashed lines, Ω ± κ/4 as dotdashed lines, Ω as a solid line (for the CR), and Ω + κ as a tripledotdashed line. The uppermost full horizontal line represents the mean bar pattern speed Ω_{B} = 21.9−14.2 (full line), which was determined directly from the time variation of the bar positionangle. The lowest one is an estimated intermittent waves packet Ω_{oW} ≈ 4.4−2.6 km s^{−1} kpc^{−1}. The intermediate wave at Ω_{iS} ≈ 13.2−8.4 km s^{−1} kpc^{−1} is computed as a beat mode. Bottom row: same, but for m = 4. 
We have identified at least three kinds of waves that might impact the evolution of the disc in a significant way. The first mode identified as the bar covers a large domain in Ω_{p} = Ω_{B}, from 23.2 down to 13.5 km s^{−1} kpc^{−1}. Since the integration time window is large and the bar is secularly slowing down, it is normal that this mode is spread between the extreme values of Ω_{B}. Beyond the bar corotation (CR_{B}), the bar permanently excites a m = 2 mode (thus of an identical speed pattern), which is visible beyond the bar corotation, but also well beyond the 1:1 resonance (Ω + κ = Ω_{B} resonance).
Secondly, an intermediate spiral wave (named “iS” hereafter), whose maximum power is located at Ω_{iS} ∈ [15−8] km s^{−1} kpc^{−1}. This mode has a spatial extension that goes roughly from UHR_{B} to its own 1/1 resonance, that is, between ≈20 and ≈40 kpc (cf. also Fig. 8). This mode gains power over time. It is clearly more visible at the end of the simulation, but it seems to be present as soon as t = 3 Gyr.
Finally, at lower values of Ω (below ≈10 km s^{−1} kpc^{−1}), other waves appear, which are not permanent. They reappear regularly at slightly decreasing Ω_{p} values, giving an average contribution that is spread out in Ω. Nevertheless, cumulated over 7.4 Gyr, their signature (in terms of power) is at least equivalent to the bar. The behaviour is similar to wave packets that carry angular momentum outwards in a finite time and not like standing waves. The integration over about 2 Gyr shows a cumulative power that exaggerates the comparison with the stellar bar, which is permanent. These wave packets are nevertheless indispensable to evacuating the angular momentum towards the outer edge of the disc. For the most powerful of these intermittent waves, we have roughly determined Ω_{oW} and plotted it in Fig. 9 (Ω_{oW} ∈ [4.4−2.6] km s^{−1} kpc^{−1}, where “oW” stands for “outer wave”). These values, estimated by hand, are very approximate because they are slightly different according to m. We can reasonably approximate that this recurrent wave structure, whose successive values of Ω_{oW} decrease, is equivalent to a permanent wave that would slow down over time, as the bar does. This structure is complex and difficult to analyze because it is the only one whose power spectrum does not cancel for m > 4 or for odd m. Its trace is perfectly visible up to m = 8, which is the limit we imposed on ourselves in our study; whereas for m > 4, the other inner structures have almost no contribution.
The nature of the intermediate spiral wave raises a question since its normal mode is very close to what was obtained by beating the bar mode with the averaged outer waves discussed above. In a linear approach, all the waves present in the disc evolve independently and do not interact. However, if higher order terms of kinetic equations are considered, this is no longer true and waves can exchange energy and angular momentum (Sygnet et al. 1988). Selection rules on wave numbers and frequencies then apply. Therefore, applying these selection rules, ω_{beat} = 2Ω_{B} + 2Ω_{oW} decreases from ≈53 to ≈34, giving Ω_{iS} = ω_{beat}/4 ≈ 13.3−8.5 km s^{−1} kpc^{−1}. This is approximately the location of the intermediate wave in m = 2 and m = 4 spectrograms. Therefore, this could correspond to the mode coupling as illustrated by Masset & Tagger (1997) in an Nbody simulation.
6.2. Resonance overlaps
With at least three patterns in the stellar disc, a number of resonance overlaps are unavoidable. Over 1 Gyr, resonance radii can increase by up to ≈1 kpc due to the changes in Ω_{p} and in Ω(R) and κ(R) curves. Moreover, resonances unavoidably have a width that must be expressed in frequency units. For analytical dynamical systems, this width is often computed in the pendulum approximation and is typically proportional to the squareroot of the perturbation amplitude. For a galaxy, we can make the reasonable assumption that the width depends both on the Ω_{p} bandwidth and on the local slope of Ω(R)+κ(R)/m. The latter dependency results in narrower resonances when only a single bar is involved; the ILR being probably the narrowest of all and the OLR being the widest. On the other hand, resonances with spiral structures are much wider. As for the width of Ω_{p}, it depends on its time derivative and therefore on possible fluctuations.
Therefore, the notion of overlap, which is expressed in the spatial domain, should not be taken literally. The above margin of the order of 1 kpc can be applied. In which case, between t = 6.32 and 8.48 Gyr (Fig. 9 middle panel), OLR_{B} and CR_{iS} are close to each other, as well as UHR_{B}−ILR_{iS}. Between t = 8.38 and 10.54 Gyr (Fig. 9 right panel), UHR_{B} is still close to ILR_{iS}, while OLR_{B} and CR_{iS} are now separated by ≈1.5 kpc. Within the large Ω_{oW} uncertainties, ILR_{oW} might also be close to both the bar 1/1 resonance and the outer m = +4 intermediate wave resonance. Between these two extreme time windows, any other type of overlap may occur.
The bar and the intermediate wave seem to be locked as UHR_{B}−ILR_{iS} overlap is constant within less than 0.5 kpc, while OLR_{B}−CR_{iS} do the same within a slightly wider range. If the intermediate spiral is a beat mode, it would imply that the pattern of outer waves is also locked to the bar one. In view of the uncertainties in determining the pattern speed of the outer wave, we cannot firmly confirm this.
The frequency analysis thus suggests that the dynamical particle and wave interactions have many sources in the disc. This is, in particular, one of the reasons why we claimed to be unable to confirm that OLR_{B} might be a barrier to radial migration (Wozniak 2020). Indeed, OLR_{B} occurs at a rather small radius with respect to the whole disc extension. Any single wave (as a bar) cannot efficiently carry angular momentum over a large radial span. Therefore, spiral waves take over the bar in angular momentum exchanges, at least up to the UHR_{oW} of the lowest frequency outer waves (Fig. 8). The UHR_{oW} resonance seems to mark the end of the set of external waves (as for Patsis et al. 1994) whose properties, both morphological and temporal, are different from the intermediate spiral.
7. Discussion
First of all, let us recall that we focus on the evolution of the disc once the bar formation phase is over, so 3.16 < t < 10.54 Gyr for RunC. The disc is then in a state of adiabatic evolution. It remains subject to the gravitational influence of the bar and to its own selfgravitating instabilities.
The diffusion timescale, in Chirikov’s sense, is shorter for J_{R} than for L_{z} in the disc. On average, T_{D}(J_{R}) is even slightly shorter in the disc than in the bar. The set of Subset B particles, selected such that T_{D}(J_{R}) < 3 and T_{D}(L_{z}) > 10.54, is associated with all morphological structures related to the secular evolution of the disc, excluding the bar. In the context of the epicycle approximation, J_{R} diffusion could be interpreted as X diffusion and thus secular radial heating. In other words, blurring (Schönrich & Binney 2009) would be favoured over churning for Subset B in the time window 3.16 < t < 10.54 Gyr. Alternatively, one could also imagine that J_{R} diffusion could be more strongly linked to that of κ. The fact that the κ frequency can diffuse is mainly due to its dependence on R, which would translate into L_{z} diffusion. However, as Halle et al. (2015) show, the disc evolution is complex since blurring and churning coexist and their relative importance evolves over time. The magnitude of churning decreases with time, unlike blurring. In our case, starting our analysis well after the bar formation, the intensity of churning may have strongly decreased.
The stellar disc is not only made of Subset B particles. On the one hand, there are disc particles in the complementary subpopulation Subset A for which T_{D}(J_{R}) < 3 and T_{D}(L_{z}) < 10.54 (lowerleft corner of Fig. 3). Thus, a fraction of the disc shows L_{z} variations on short timescales, that is < 10.54 Gyr. One should keep in mind that the numerical values of the boundaries have not been firmly established yet. On the other hand, the DF of these Subset A particles also evolves: Both DF(J_{R}) and DF(L_{z}) widen by ≈20% (according to their full width at half maximum), while their maxima decrease. The widening of DF(L_{z}) for Subset A is mainly due to an increase in L_{z} for particles initially with L_{z} ≈ 0. This coevolution of J_{R} and L_{z} suggest that the scattering mechanism is also at work in this simulation, but on a longer timescale that is compatible with a decreasing importance. This dynamical mechanism cannot only be the scattering by corotation(s) because ΔJ_{R} ≠ 0. Inside the stellar bar, ILR scattering is likely to be the most efficient mechanism.
In the framework of a particlemesh code, the waveparticle exchanges shape the evolution of dynamical properties. Several waves have been identified by their power spectrum. Some are highly timedependent. Intermittence has not been studied exhaustively yet, but it deserves special attention as it certainly has a role, as shown by Sellwood & Binney (2002). By simply looking at the evolution of these waves over time windows of 1–2 Gyr, we can nevertheless qualitatively deduce the impact of the resonances that these waves introduce into the disc. Several resonance species are at work in the simulation.
Changes in the distribution function DF(, ) are perfectly identifiable. The two vertical tails, located at L_{z} ≈ 4000 and 4500, and J_{R} ≈ 100, which are similar in shape to Sellwood (2012), could be the signature of an ILR scattering of each of the two waves present, that is to say the intermediate and the external. However, as these tails appear progressively (Fig. 5), we cannot exclude that the origin of this scattering is a single wave, which would reappear at a lower frequency and/or greater radius, thus explaining the duplication. In addition, a lowlevel inspection of the DF(, ) map shows other signatures of the same type but at much lower levels, supporting the latter hypothesis. The depopulation of circular orbits (J_{R} = 0), accompanied by a redistribution in L_{z} (cf. Sect. 5 and Fig. 7), is associated with this ILR scattering by intermittent waves. If the ILR of the intermediate spiral and/or the outer waves are involved, the energy contained in the waves is transferred to the particles. However, since the entire disc is not heated by this mechanism, it remains cold enough to allow for the regular resurgence of the mode(s).
With three sets of waves, we suspect that resonance overlaps in physical space, such as OLR_{B}−CR_{iS} and UHR_{B}−ILR_{iS} or +1/4 iS−ILR_{oW}, may play an important role. For a full efficiency of these couplings, Sygnet et al. (1988) have shown that two patterns must overlap over a radial range. But the theory does not predict anything about the duration of the overlap(s). As far as we know, there is nothing to prevent intermittent waves from interfering nonlinearly. Sygnet et al. (1988) suggested that the corotation of the inner wave must coincide with the ILR of the outer wave, which is the case for CR_{iS}−ILR_{oW} at t = 3.16−4.24 Gyr. They showed that this kind of resonance overlap would make the nonlinear interaction between the two patterns much more efficient. This is the process advocated by Minchev & Famaey (2010) for amplifying radial migration.
It is questionable whether an (m > 0)−(m < 0) overlap (e.g. intermediate spiral outer m = + 4 −ILR_{oW}, or any OLR−ILR overlap as for a doublebarred system in Wozniak 2015), would facilitate the transfer of angular momentum to the outer regions. Indeed, it has long been shown (LyndenBell & Kalnajs 1972) that particles at any m > 0 resonance (as well as at the corotation) absorb L_{z} (and E) from the wave, while those at any m < 0 give L_{z} (and E) to the wave. When two resonances with opposite signs overlap, the angular momentum and energy acquired by the particles at the m > 0 resonance of the inner highΩ_{p} wave can be transferred to the outer lowΩ_{p} wave through any m < 0 resonance.
The wide vertical tail visible in DF(, ) around , which is approximately the value for the bar corotation at t = 10.54 Gyr, has yet to be understood. The large values of J_{R} that are reached are the signature of the hot population. This population should not see a significant variation of J_{R} since it is supposed to be scattered by the bar corotation. However, the approximate calculation of J_{R} (Eq. (6)) prevents a more detailed analysis of this population. Indeed, is very different from in this region. This is easily explained because these particles explore a large fraction of the bar and the disc, where κ values are very different. The epicycle approximation is therefore not valid here.
It would have been interesting to compare the results obtained with real data, notably those concerning the Milky Way obtained by GaiaRVS (Trick et al. 2019), or recent models also expressed in the actionangle domain (e.g. Frankel et al. 2020, and reference therein). However, as mentioned several times, RunC does not have the characteristics of the Milky Way. Therefore, we cannot explain the similarities or differences yet with these recent semianalytical modelling works. We postpone the analysis of ongoing simulations that are actually dedicated to the Milky Way to a later article.
8. Conclusions
Using the epicycle formulation of the radial action J_{R} to calculate the Chirikov diffusion rate and the associated characteristic timescale, in addition to the results already obtained so far for L_{z} and E by Wozniak (2020), we have shown the following.

The distribution of J_{R} diffusion timescales spikes around T_{D}(J_{R}) = 0.9 Gyr, which is mainly due to disc particles. It is followed by a plateau ending at T_{D}(J_{R}) ≈ 3 Gyr, whose main contributor is the stellar bar. Roughly 0.5 M_{tot} has T_{D}(J_{R}) < 3 Gyr and 0.77 M_{tot} has T_{D}(J_{R}) < 10.54 Gyr (i.e. the simulation length). Beyond the bar UHR, the space averaged timescale is ⟨T_{D}(J_{R})⟩ ∼ 1 Gyr.

By selecting particles as T_{D}(J_{R}) < 3 Gyr and T_{D}(L_{z}) > 10.54 Gyr, that is, 0.25 M_{tot}, we identify all the particles that participate in the morphological structures that are characteristic of resonances and waves with the exception of the stellar bar.

Secular radial heating has been identified in the disc, by means of the depopulation of circular orbits. We note that 57% of particles on circular orbits (J_{R} = 0) at t = 3.16 Gyr have J_{R} increased by a few tens of kpc km s^{−1} after 7 Gyr, leading σ_{R} to increase by ∼10 km s^{−1}. This decircularisation is accompanied by L_{z} transfer through coherent waveparticle interactions.

The signature of ILR scattering by disc wave(s) has been detected by ridges in the average distribution function DF(, ), similarly to Sellwood (2012).

A wave analysis identified at least three types of waves: (1) a permanent stellar bar, (2) a nonpermanent intermediate spiral, which could be a beating phenomenon, with (3) a set of intermittent multiarmed wave packets that carry L_{z} towards the edge of the disc. Several resonance overlaps ensure a continuous cascade of waves, which can also contribute to promote L_{z} exchanges. But none of these overlaps are of the ILR−CR type.
The Chirikov diffusion rate allows for the separation of particle sets with similar dynamic behaviour. It is a useful complementary tool for dynamical analysis. Numerous avenues for further investigation will be the subject of subsequent articles: the study of the bar, the diffusion of R, κ, and X, the impact of live dark matter particles, etc. In particular, here, we have calculated D_{2} based on an average of the properties over 100 Myr, but other scales are possible (Wozniak 2020). It is through these other timescales that possible effects due to chaos could appear (Lichtenberg & Lieberman 1992).
Acknowledgments
I thank the anonymous referee for helpful comments. I would like to acknowledge the Meso@LR computing centre of the University of Montpellier for providing access to computing resources. This project has been funded by a grant from the Scientific Council of the University of Montpellier.
References
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press) [Google Scholar]
 Buta, R. J. 2019, MNRAS, 488, 590 [CrossRef] [Google Scholar]
 Chirikov, B. V. 1979, Phys. Rep., 52, 263 [NASA ADS] [CrossRef] [Google Scholar]
 Frankel, N., Sanders, J., Ting, Y.S., & Rix, H.W. 2020, ApJ, 896, 15 [CrossRef] [Google Scholar]
 Freeman, K. C. 1966, MNRAS, 133, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Halle, A., Di Matteo, P., Haywood, M., & Combes, F. 2015, A&A, 578, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kalnajs, A. J. 1971, ApJ, 166, 275 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Lichtenberg, A., & Lieberman, M. 1992, Regular and Chaotic Dynamics (New York: Springer) [CrossRef] [Google Scholar]
 LyndenBell, D., & Kalnajs, A. J. 1972, MNRAS, 157, 1 [Google Scholar]
 Mackereth, J. T., Bovy, J., Leung, H. W., et al. 2019, MNRAS, 489, 176 [CrossRef] [Google Scholar]
 Masset, F., & Tagger, M. 1997, A&A, 322, 442 [NASA ADS] [Google Scholar]
 MichelDansac, L., & Wozniak, H. 2006, A&A, 452, 97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Minchev, I., & Famaey, B. 2010, ApJ, 722, 112 [NASA ADS] [CrossRef] [Google Scholar]
 Minchev, I., Famaey, B., Combes, F., et al. 2011, A&A, 527, A147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Patsis, P. A., Hiotelis, N., Contopoulos, G., & Grosbol, P. 1994, A&A, 286, 46 [Google Scholar]
 Pfenniger, D. 1990, A&A, 230, 55 [NASA ADS] [Google Scholar]
 Pfenniger, D., & Friedli, D. 1991, A&A, 252, 75 [NASA ADS] [Google Scholar]
 Roškar, R., Debattista, V. P., Quinn, T. R., & Wadsley, J. 2012, MNRAS, 426, 2089 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, J. L., & Binney, J. 2014, MNRAS, 441, 3284 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, J. L., & Binney, J. 2016, MNRAS, 457, 2107 [NASA ADS] [CrossRef] [Google Scholar]
 Schönrich, R., & Binney, J. 2009, MNRAS, 396, 203 [Google Scholar]
 Sellwood, J. A. 2012, ApJ, 751, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Sellwood, J. A., & Binney, J. J. 2002, MNRAS, 336, 785 [NASA ADS] [CrossRef] [Google Scholar]
 Skokos, C., Patsis, P. A., & Athanassoula, E. 2002, MNRAS, 333, 847 [NASA ADS] [CrossRef] [Google Scholar]
 Soubiran, C., Bienaymé, O., Mishenina, T. V., & Kovtyukh, V. V. 2008, A&A, 480, 91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sparke, L. S., & Sellwood, J. A. 1987, MNRAS, 225, 653 [NASA ADS] [Google Scholar]
 Sygnet, J. F., Tagger, M., Athanassoula, E., & Pellat, R. 1988, MNRAS, 232, 733 [NASA ADS] [CrossRef] [Google Scholar]
 Trick, W. H., Coronado, J., & Rix, H.W. 2019, MNRAS, 484, 3291 [NASA ADS] [CrossRef] [Google Scholar]
 Vasiliev, E. 2019, MNRAS, 482, 1525 [NASA ADS] [CrossRef] [Google Scholar]
 Wozniak, H. 2020, ApJ, 889, 81 [CrossRef] [Google Scholar]
 Wozniak, H. 2015, A&A, 575, A7 [CrossRef] [EDP Sciences] [Google Scholar]
 Wozniak, H., & Pfenniger, D. 1997, A&A, 317, 14 [Google Scholar]
 Zang, T. A., & Hohl, F. 1978, ApJ, 226, 521 [CrossRef] [Google Scholar]
All Figures
Fig. 1. Normalised particle relative frequency per Gyr as a function of T_{D}(J_{R}) (black histogram). A binsize of 0.01 Gyr has been used. One particle thus represents a probability density of 2.5 × 10^{−6} Gyr^{−1}. Distributions for T_{D}(E) (red) and T_{D}(L_{z}) (blue) are plotted for comparison purposes. 

In the text 
Fig. 2. ⟨T_{D}(J_{R})⟩ (black line), ⟨T_{D}(E)⟩ (red line), and ⟨T_{D}(L_{z})⟩ (blue line) as a function of . The shaded area delimits the region occupied by bar corotation (CR_{B}) over ≈7.4 Gyr. Vertical lines show L_{z} of circular orbits at the innermost bar’s UHR (UHR_{B} at t = 3.16 Gyr, dotdashed) and the positions reached by the outermost bar’s OLR (OLR_{B} at t = 10.54 Gyr, longdashed). 

In the text 
Fig. 3. Normalised particle relative frequency (Gyr^{−2}) in the T_{D}(L_{z})−T_{D}(J_{R}) plane in log scale for neverescaped particles only. Contours are spaced by 0.15 dex. The white line divides the domain into two sets (named Subset A and Subset B), which are used for Figs. 4–6 (see text for details). 

In the text 
Fig. 4. Projected mass surface density (in M_{⊙} pc^{−2}) at t = 10.54 Gyr inside ±50 kpc. Particle sets are defined in Fig. 3 and in the main text. The log colourscale is common to both figures. Black isodensities are spaced by 0.4 dex for Subset A (left panel) and 0.1 for Subset B (right panel). The spatial scale is in kpc. 

In the text 
Fig. 5. Mass density per (kpc km s^{−1})^{2} in the plane in log scale for neverescaped particles only. This is a timeaveraged 2D representation of DF over 7.4 Gyr. Contours are spaced by 0.2 dex. White isodensity contours represent Subset A; black contours correspond to Subset B. 

In the text 
Fig. 6. Distribution functions (DFs) for all particles (green triple dotteddashed lines), “neverescaped” particles (dotted lines), and the two subpopulations defined in Fig. 4 (red and black lines) as a function of (top) and (bottom). We note that and have been averaged over ≈1 Gyr, starting at t = 3.2 (left), 6.3 (middle), and 9.5 Gyr (right). DFs have been normalised to DF and DF for the interval t = 3.2−4.2 Gyr, respectively. 

In the text 
Fig. 7. Particle number as a function of (top) and (bottom), for particles selected as kpc km s^{−1} (black lines). The same particle selection is drawn at t = 6.3−7.3 (in green) and t = 9.5−10.5 Gyr (in red). The contribution of particles with strictly increasing J_{R} to DF() is shown by a dashed line. 

In the text 
Fig. 8. Projected mass surface density for particles selected as kpc km s^{−1} at t = 3 Gyr (left panel) inside ±50 kpc. Their distribution at t = 10.54 Gyr (right panel). The log colourscale is common to both figures (in M_{⊙} pc^{−2}). Black isodensities are spaced by 0.2 dex. Circles show the position for the ILR (dotdashed line), the UHR and outer m = + 4 resonance (shortdashed), the CR (full line), and the OLR (longdashed). In red, the following resonances for the stellar bar are shown: UHR_{B}, CR_{B}, and OLR_{B}, as defined in Sect. 6. Using the same linestyle, the white circles represent the intermediate spiral: ILR_{iS} (close to UHR_{B}), CR_{iS} (close to OLR_{B}), outer m = +4 resonance, and OLR_{iS}. The black circles represent the outer wave ILR_{oW} (close to intermediate spiral outer m = +4 resonance), UHR_{oW}, and CR_{oW}. 

In the text 
Fig. 9. Top row: m = 2 power spectra in log scale as a function of radius for RunC. The time windows are 3.16−4.24 Gyr (left panel), 6.32−8.48 Gyr (middle), and 8.38−10.54 Gyr (right). The vertical scales give values of Ω in Myr^{−1} (left) and in km s^{−1} kpc^{−1} (right). The radial scale (in kpc) is converted in L_{z} using the circular rotation curve. The mean curves Ω ± κ/2 are drawn as black shortdashed lines, Ω ± κ/4 as dotdashed lines, Ω as a solid line (for the CR), and Ω + κ as a tripledotdashed line. The uppermost full horizontal line represents the mean bar pattern speed Ω_{B} = 21.9−14.2 (full line), which was determined directly from the time variation of the bar positionangle. The lowest one is an estimated intermittent waves packet Ω_{oW} ≈ 4.4−2.6 km s^{−1} kpc^{−1}. The intermediate wave at Ω_{iS} ≈ 13.2−8.4 km s^{−1} kpc^{−1} is computed as a beat mode. Bottom row: same, but for m = 4. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.