© H. Wozniak 2020

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 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.

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 near-circular 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 guiding-centre radius (thus generating no net radial migration in principle), whereas that of churning is to move the guiding-centre 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:

(1)

where Ω is the orbit and particle angular frequency. At co-rotation, Ω = Ω_p so that there is no J_R variation; whereas at Lindblad resonances with an m-armed perturbation,

(2)

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 lth-order resonance, where l = 0 is the co-rotation.

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. Co-rotation 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 time-dependent proportions (Halle et al. 2015), and scattering can be non-resonant (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 co-rotation generates a high degree of stochasticity. The diffusion capacity of a bar’s co-rotation 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:

(3)

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 time-dependent simulation (Wozniak 2020). In Eq. (3) for n = 2, is the value of radial action averaged over a period of Δt₂ = 10² (in the unit of time of the system). We can also estimate a diffusion timescale T_D by re-normalising D₂ by , where the time-average is now computed over the longest timescale (e.g. the experiment length):

(4)

Several methods exist to estimate actions numerically (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 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

(5)

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 N-body 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:

(6)

It is convenient to time-average here over Δt₂ = 10² 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¹¹ M_⊙. All 4 × 10⁷ 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 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-body computation, they were tracked by ballistic approximation until they re-entered 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 “never-escaped” 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⁻² − 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 TD(JR) (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 TD(E) (red) and TD(Lz) (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 time-averaged 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 co-rotation (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. ⟨TD(JR)⟩ (black line), ⟨TD(E)⟩ (red line), and ⟨TD(Lz)⟩ (blue line) as a function of . The shaded area delimits the region occupied by bar co-rotation (CRB) over ≈7.4 Gyr. Vertical lines show Lz of circular orbits at the innermost bar’s UHR (UHRB at t = 3.16 Gyr, dot-dashed) and the positions reached by the outermost bar’s OLR (OLRB at t = 10.54 Gyr, long-dashed).

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 slow-diffusion. 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 never-escaped 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 quasi-linear 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 TD(Lz)−TD(JR) plane in log scale for never-escaped 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 co-rotation, the distribution in the disc only shows 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 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 “never-escaped” particles and the two particular selections (Subset A and Subset B) are plotted as a function of and time-averaged 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 never-escaped particles only. This is a time-averaged 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⁻¹) 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⁻¹, the large vertical tail with high 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. The tail width is likely to be related to co-rotation 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⁻¹, 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

One-dimensional 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 N-body 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 above-mentioned 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 “never-escaped” 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⁻¹). 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⁻¹ 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 dotted-dashed lines), “never-escaped” 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⁻¹. The mode is mainly due to Subset A. Its position shifts from kpc km s⁻¹ at Gyr to ≈90 kpc km s⁻¹ 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⁻¹. 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 near-circular 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⁻¹. 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⁹ 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⁻¹, 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 JR to DF() is shown by a dashed line.

Figure 7 shows that increases significantly for ≈57.4% of particles initially on near-circular orbits. It leads to an increase in σ_R of the order of ≈10 km s⁻¹ 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 co-rotation 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 wavelet-type 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 co-rotation 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 (dot-dashed line), the UHR and outer m​ = ​ + 4 resonance (short-dashed), the CR (full line), and the OLR (long-dashed). In red, the following resonances for the stellar bar are shown: UHRB, CRB, and OLRB, as defined in Sect. 6. Using the same linestyle, the white circles represent the intermediate spiral: ILRiS (close to UHRB), CRiS (close to OLRB), outer m = +4 resonance, and OLRiS. The black circles represent the outer wave ILRoW (close to intermediate spiral outer m = +4 resonance), UHRoW, and CRoW.

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 Lz using the circular rotation curve. The mean curves Ω ± κ/2 are drawn as black short-dashed lines, Ω ± κ/4 as dot-dashed lines, Ω as a solid line (for the CR), and Ω + κ as a triple-dot-dashed 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 position-angle. 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⁻¹ kpc⁻¹. 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 co-rotation (CR_B), the bar permanently excites a m = 2 mode (thus of an identical speed pattern), which is visible beyond the bar co-rotation, 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⁻¹ kpc⁻¹. 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⁻¹ kpc⁻¹), 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⁻¹ kpc⁻¹, 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⁻¹ kpc⁻¹. 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.

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 square-root 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 self-gravitating 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 (lower-left 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 co-rotation(s) because ΔJ_R ≠ 0. Inside the stellar bar, ILR scattering is likely to be 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. 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 low-level 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 non-linearly. Sygnet et al. (1988) suggested that the co-rotation 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 non-linear 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 double-barred system in Wozniak 2015), would facilitate the transfer of angular momentum to the outer regions. Indeed, it has long been shown (Lynden-Bell & Kalnajs 1972) that particles at any m >  0 resonance (as well as at the co-rotation) 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 co-rotation 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 co-rotation. 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 Gaia-RVS (Trick et al. 2019), or recent models also expressed in the action-angle 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 semi-analytical 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.

1. 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.

2. 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.

3. Secular radial heating has been identified in the disc, by means of the de-population 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⁻¹ after 7 Gyr, leading σ_R to increase by ∼10 km s⁻¹. This de-circularisation 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(, ), similarly to Sellwood (2012).

5. A wave analysis identified at least three 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 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₂ 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

All Figures

	Fig. 1. Normalised particle relative frequency per Gyr as a function of TD(JR) (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 TD(E) (red) and TD(Lz) (blue) are plotted for comparison purposes.
In the text

	Fig. 2. ⟨TD(JR)⟩ (black line), ⟨TD(E)⟩ (red line), and ⟨TD(Lz)⟩ (blue line) as a function of . The shaded area delimits the region occupied by bar co-rotation (CRB) over ≈7.4 Gyr. Vertical lines show Lz of circular orbits at the innermost bar’s UHR (UHRB at t = 3.16 Gyr, dot-dashed) and the positions reached by the outermost bar’s OLR (OLRB at t = 10.54 Gyr, long-dashed).
In the text

	Fig. 3. Normalised particle relative frequency (Gyr−2) in the TD(Lz)−TD(JR) plane in log scale for never-escaped 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 never-escaped particles only. This is a time-averaged 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 dotted-dashed lines), “never-escaped” 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 JR 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 (dot-dashed line), the UHR and outer m​ = ​ + 4 resonance (short-dashed), the CR (full line), and the OLR (long-dashed). In red, the following resonances for the stellar bar are shown: UHRB, CRB, and OLRB, as defined in Sect. 6. Using the same linestyle, the white circles represent the intermediate spiral: ILRiS (close to UHRB), CRiS (close to OLRB), outer m = +4 resonance, and OLRiS. The black circles represent the outer wave ILRoW (close to intermediate spiral outer m = +4 resonance), UHRoW, and CRoW.

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 Lz using the circular rotation curve. The mean curves Ω ± κ/2 are drawn as black short-dashed lines, Ω ± κ/4 as dot-dashed lines, Ω as a solid line (for the CR), and Ω + κ as a triple-dot-dashed 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 position-angle. 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