Diffusion of radial action in a galactic disc

Context. Stellar migration of the galactic disc stars has been invoked to explain the dispersion of stellar metallicity observed in the solar neighborhood. 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 analyze 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_\mathrm{R}$ in an idealised N$-$body 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_\mathrm{R}$ diffusion timescales $T_\mathrm{D}(J_\mathrm{R})$ is less than 3 Gyr for roughly half the galaxy mass. It is always much shorter than the angular momentum diffusion timescale $T_\mathrm{D}(L_\mathrm{z})$ outside the stellar bar. In the disc, $\langle T_\mathrm{D}(J_\mathrm{R}) \rangle \sim 1$ Gyr. All non-axisymmetric morphological structures characteristic of resonances and waves in the disc are associated to particles with $T_\mathrm{D}(J_\mathrm{R})<3$ Gyr and $T_\mathrm{D}(L_\mathrm{z})>10$ Gyr. Short $T_\mathrm{D}(J_\mathrm{R})$ can be explained by the gradual decircularisation of initially circular orbits ($J_\mathrm{R}=0$) under the effect of intermittent ILR scattering by wave trains propagating in the disc, well beyond the bar OLR. This leads to a moderate secular heating of the disc beyond the bar OLR for 7 Gyr, comparable to solar neighbourhood observations. The complex multi-wave structure, mixing permanent and intermittent modes, allows for multiple resonance overlaps.


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 self-consistent N−body experiments of isolated galactic discs. Using the same definition, we extend our previous work by now focusing on the radial action J R . J R 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 the specific radial kinetic energy. This approximation is only valid if J R L z or for near-circular orbits (Kalnajs 1971). J R has been considered as a thermometer for measuring the disc heating in the radial direction. Whether the stellar disc heats or not 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 guiding-centre radius (thus generating no net radial migration in principle), whereas that of churning is to move the guiding/mean radius inwards or outwards.
According to Sellwood & Binney (2002), any variation of J R in a stellar disc excited by a constant non-axisymmetric perturbation of pattern speed Ω p , can be related to ∆L z as: where Ω is the orbit/particle angular frequency. At corotation Ω = Ω p so that there is no J R variation, whereas at Lindblad resonances with an m−armed perturbation, This applies to any kind of disturbance, such as bar or spiral(s). It can be generalized to higher order resonances, so that ∆J R /∆L z = ∓l/m for any l th -order 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 for 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), scattering can be nonresonant (e.g. with giant molecular clouds) or made by Lindblad resonances (e.g. the Inner Lindblad Resonance Sellwood 2012). Moreover, the bar-spirals resonance overlap (Minchev & Famaey 2010;Minchev et al. 2011) is expected to increase 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 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 over time of small fluctuations in energy, angular momentum and radial action. 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 organized 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. Sect. 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.

Dynamical concepts
The Chirikov diffusion rate of J R for individual particles, can be defined as : Although the original definition deals only 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 time-dependent simulation (Wozniak 2020). In Eq. (3) for n = 2, J R is the value of radial action averaged over a period of ∆t 2 = 10 2 (in time unit of the system). We can also estimate a diffusion timescale T D by renormalizing D 2 by J R 2 , where the time-average is now computed over the longest timescale (e.g. the experiment length): Several methods exist to estimate numerically actions (see Sanders & Binney 2016, for a comprehensive review). Vasiliev (2019) provides a state-of-the-art software package for building dynamical models based on an action-angle description, including position-velocity and action-angle conversion routines. However, the accurate and fast estimation of action-angle variables remains an open problem, especially for non-axisymmetric 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 N−body simulations, this objective is still extremely difficult, if not impossible, to achieve for strictly computational reasons. Therefore, we have decided to remain in 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 least bad 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 of the form Φ(R, z) = Φ R (R) + Φ z (z), J R reduces to where E z = 1 2ż 2 + Φ z . In the limit of very small radial excursions around the guiding centre (X = (R max − R min ) → 0), J R reduces even further to J R → E R /κ = 1 2 κX 2 where κ, E R and X are respectively the epicycle radial frequency, the radial kinetic energy, and the epicycle amplitude (Binney & Tremaine 2008). The integral must be taken around a full radial period, a quantity largely unknown (and not accessible) in an N−body simulation. Since we need J R 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 time-average here over ∆t 2 = 10 2 time units (i.e. 105.4 Myr). Other time windows will be 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 of 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 and κ 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 generalized κ formulation due to Pfenniger (1990) and tried to estimate J R as the ratio of averages E R /κ. The results quantitatively differ but this does not drastically change our conclusions. We will mention in the following results most strongly impacted by our choice for estimating J R .

Numerical implementation
In the following, we will use 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 such 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 plots expressed in 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 σ 2 θ = σ 2 R κ 2 /(4Ω 2 ), where σ R , σ θ and σ z are three components of the velocity dispersion along respectively the radial, azimuthal and vertical directions.
We only study the phase after the formation of the bar (t > 3.16 Gyr, until simulation end 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 at all on the rest of the disc.
Unless otherwise stated, we have excluded from the analyses particles escaping the 3D log−polar 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 N−nody computation, they have been tracked by ballistic approximation until they re-enter the grid, 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 'never-escaped' in the following.

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 (Wozniak 2020). The distribution shapes are different, especially for T D < ∼ 0.9 Gyr. On the contrary to E and L z , particles with T D (J R ) < 100 Myr are few. 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. 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 that 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. 2 shows T D (J R ) averaged over sets of particles (designated by T D (J R ) ) sampled by L z ranges. L z is now timeaveraged over ≈ 7.4 Gyr (i.e. from t = 3.16 to 10.54 Gyr), and thus must not be confused with an instantaneous L z that is not conserved. The range in L z occupied by the bar Lindblad resonances during the evolution is delimited by the innermost position of the bar UHR (UHR B at t = 3.16 Gyr) and the outermost bar OLR (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.
As shown in Wozniak (2020), T D (E) decreases overall from the centre (L z ≈ 0) to the outermost regions, with a strong depression in the area occupied by the bar Lindblad resonances, while T D (L z ) increases monotonically to values that can be considered as slow-diffusion. This is none of those cases 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 never-escaped particles), have shorter diffusion times for L z and J R than for E. For L z and J R the values remain compatible with the timescales in the disk. The scattering of retrograde particles is a subject in 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/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 analyzed 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 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 ). It will be seen in the following that it clearly belongs to the bar. Substructures can be linked to families of orbits, but this 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 quasi-linear decay observed in Fig. 1 when T D (J R ) > 3 Gyr. The other noticeable structure, of the same density and long T D (L z ) (> 10 Gyr), is related to the disc. It contributes in a large fraction 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. 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 which morphological element, orbit family or wave type they correspond to. 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 some trials, and 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 would have been possible (and 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, the rest being excluded particles. Fig. 4 shows the mass surface density projected in the x−y plane for the two particle sets defined above. Subset A, the most massive (0.65 M tot ), contains mainly particles with a highly symmetrical mass distribution, especially the stellar bar. Beyond the corotation, the distribution in the disc shows only very small deviations from axisymmetry. Inside UHR B the morphology is elliptical-like, which is the signature of strong bars (Skokos et al. 2002;Michel-Dansac & 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 exist also inside the bar (such as ansea, Buta 2019) but extend well beyond OLR B . As with Subset A, there appear to be several subpopulations, which do not have 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 will analyze it in greater detail.

Distribution functions (DF)
4.1. DF in L z and J R In Fig. 5, the DF for 'never-escaped' particles and the two particular selections (Subset A and Subset B) are plotted as a function of L z and J R time-averaged between t = 3.16 and 10.54 Gyr. The 'averaged' DF 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 as our initial stellar disc is not such a stable disc. 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.
Subset B is identifiable through several substructures. A density peak is present for J R ≈ 0 (circular orbits) and L z > 3600 (beyond OLR B at t = 10.54 Gyr). This region is bordered by two almost vertical tails (centered at 4000 and 4500 kpc km s −1 ) which, by similarity 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 J R values seems to include a significant fraction of the so-called '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. Tail width is likely to be related to corotation shifts over time. J R can reach very high values there. The mass distribution is very sensitive to how J R is calculated. Using the E R /κ ratio as an estimator of J R , the mass would have been concentrated around a maximum density located at (L z ≈ 3100, J R ≈ 55), 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 1400 < ∼ L z < ∼ 2300 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 seems mainly due to the bar, from the rest of the disc.

DF time evolution
1D DF(J R ) or DF(L z ) can be recovered by the projection of the L z − J R density map on the axes. DF(L z ) is then similar to typical profiles obtained by a wealth of 3D N−body simulations, e.g. 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 above-mentioned 'hot' population.
In order to identify the time evolution of some identifiable structures in DF(L z ), we have calculated L z over only 1 Gyr at three different moments of the simulation. The intervals were centred at t = 3.7, 6.8 and 10.0 Gyr. The criteria used to define the selection of Subset A, Subset B, and 'never-escaped'   In Fig. 6 (bottom panels), the contribution of the two subsets is clearly separated. The whole stellar bar forms the peak of DF(L z ) 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 (L z > ∼ 2300 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 L z < ∼ 2500 kpc km s −1 has been previously identified as a ridge overlapping with Subset A. It is linked to bar structures (cf. Fig. 4).
Other evidence, 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 J R distribution (Fig. 6 top panels) blends all the components discussed so far. Most particles temporarily escaped have J R > ∼ 150 kpc km s −1 . The mode is mainly due to Subset A. Its position shifts from J R ≈ 70 kpc km s −1 at t = 3.7 Gyr to ≈ 90 kpc km s −1 at t = 10.0 Gyr. In this timeframe, the distribution width increases by ≈ 23% for Subset A while a substructure appears in the Subset B distribution, leading to a small bump around J R ≈ 100 kpc km s −1 . This behaviour is symptomatic of moderate but regular radial heating of the disc.
A&A proofs: manuscript no. diffusion_jr  Fig. 4 (red and black lines) as a function of J R (top) and L z (bottom). J R and L z 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 resp. DF(J R = 0) and DF(L z = 0) for the interval t = 3.2 − 4.2 Gyr.

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(J R ) increases over time as the bin J R ≈ 0 depopulates. Since DF(J R ) has been normalised by DF(J R = 0) at t = 3.7 Gyr, Fig. 6 clearly shows that the number of particles with J R ≈ 0 decreases with time, even long after the bar has been formed. This strongly points to an increase in epicycle amplitude as J R ∼ 1 2 κX 2 . The evolution of these particles on near-circular orbits thus deserves a particular analysis since the variation of J R is discriminating with respect to the blurring/churning issue.
In order to take the numerical uncertainties inherent to this type of simulation into account, let us define hereafter circular orbits as J R < 10 kpc km s −1 . If we select only the particles with J R (t = 3.7) < 10, the evolution of DF(J R ) and DF(L z ) can be extracted at time t = 6.8 and t = 10.0 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 L z tail in the outermost part of the disc (L z 3300), and 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 L z ≈ 3000. Fig. 7 shows that J R increases significantly for ≈ 57.4% of particles initially on near-circular orbits. It leads to an increase of σ R of the order of ≈ 10 km s −1 beyond OLR B . This is comparable with observations of solar neighbourhood (Soubiran et al. 2008;Mackereth et al. 2019, for instance), although RunC does not pretend to reproduce the properties of Milky-Way. The corresponding impact on DF(L z ) 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 t = 3.7 Gyr as having J R ≈ 0 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 wavelet-type structures. J R of more than half of the particles having increased significantly, it is normal that the initial axisymmetry has been broken. In anticipation of Sect. 6, we have plotted some resonances identified Fig. 7. Particle number as a function of J R (top) and L z (bottom), for particles selected as J R (3.2 − 4.2) < 10 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(L z ) is mentioned by a dashed line. from three patterns detected in the disk: 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.

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. Fig. 9 show 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 (twice as large) as the slowdown rate of Ω B is higher when the bar is young.
We have identified at least three kind of waves that might imprint the evolution of the disc in a significant way: -The 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 identical speed pattern), visible beyond the bar corotation, but also well beyond the 1:1 resonance (Ω + κ = Ω B resonance). -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, i.e. 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 seems to be present as soon as t = 3 Gyr. -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 evacuate 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 neither for m > 4 nor for odd m. Its trace is perfectly visible up to m = 8, the limit we imposed on ourselves in our study, while for m > 4 the other inner structures have almost no contribution.
The nature of the intermediate spiral wave raises question as its normal mode is very close to that obtained by beating the bar mode with the averaged outer waves discussed above. In a linear approach, all the waves present in the disk evolve independently and do not interact. However, if higher order terms of kinetic equations are considered, this is no more 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 N−body simulation.

Resonance overlaps
With at least three patterns in the stellar disc, a number of resonance overlap is 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 have unavoidably a width which must be express 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 Ω 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, 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, 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 outer waves pattern 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/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.

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 self-gravitating instabilities.
The diffusion timescale, in Chirikov's sense, is shorter for J R than for L z in the disc. In 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, bar excluded. 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 κ. 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) showed, 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 (lower-left corner of Fig. 3). Thus, a fraction of the disc shows L z variations on short timescales, i.e. < 10.54 Gyr. It should be remembered that the numerical values of the boundaries are not yet firmly established. 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 FWHM) 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 compatible with a decreasing importance. This dynamical mechanism cannot be only the scattering by corotation(s) as ∆J R 0. Inside the stellar bar, ILR scattering is likely the most efficient mechanism.
In the framework of a particle-mesh code, the wave−particle exchanges shape the evolution of dynamical properties. Several waves have been identified by their power spectrum. Some are highly time-dependent. Intermittance 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 to 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(J R ,L z ) are perfectly identifiable. The two vertical tails, located at L z ≈ 4000 and 4500, and J R ≈ 100, similar in shape to Sellwood (2012), could be the signature of an ILR scattering of each of the two waves present, 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 low-level inspection of DF(J R ,L z ) 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 to 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. But, since it is not the entire disc that is heated by this mechanism, it remains cold enough to allow 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 this process, Sygnet et al. (1988) have shown that two patterns must overlapped over a radial range. But the theory does not predict anything about the duration of the overlap(s). As far 1. The distribution of J R diffusion timescales is spiked around T D (J R ) = 0.9 Gyr, mainly due to disc particles. It is followed by a plateau ending at T D (J R ) ≈ 3 Gyr those main contributor is the stellar bar. Roughly 0.5 M tot has T D (J R ) < 3 Gyr, 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. 2. By selecting particles as T D (J R ) < 3 Gyr and T D (L z ) > 10.54 Gyr, i.e. 0.25 M tot , we identify all the particles that participate in the morphological structures characteristic of resonances and waves with the exception of the stellar bar. 3. Secular radial heating has been identified in the disc though the depopulation of circular orbits. 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 to σ R increase by ∼ 10 km s −1 ). This decircularisation is accompanied by L z transfer through coherent wave−particle interactions. 4. The signature of ILR scattering by disc wave(s) has been detected by ridges in the average distribution function DF(J R ,L z ), similarly to Sellwood (2012). 5. A wave analysis identified at least 3 types of waves: 1) a permanent stellar bar, 2) a non-permanent intermediate spiral which could be a beating phenomenon with 3) a set of intermittent multi-armed 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 is of the ILR−CR type.
The Chirikov diffusion rate allows 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 next articles: study of the bar, diffusion of R, κ and X, impact of live dark matter particles, etc. In particular, we have calculated here 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).