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/0004-6361/202038959 | |
Published online | 23 October 2020 |
Diffusion of radial action in a galactic disc
LUPM, Univ Montpellier, CNRS, Montpellier, France
e-mail: 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 JR 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 JR diffusion timescale TD(JR) is less than 3 Gyr for roughly half the galaxy mass. It is always much shorter than the angular momentum diffusion timescale TD(Lz) outside the stellar bar. In the disc, ⟨TD(JR)⟩ ∼ 1 Gyr. All non-axisymmetric morphological structures that are characteristic of resonances and waves in the disc are associated to particles with TD(JR) < 3 Gyr and TD(Lz) > 10 Gyr. Short TD(JR) can be explained by the gradual de-circularisation of initially circular orbits (JR = 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 multi-wave 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 Lz 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 JR.
The radial action has been introduced in the galactic problem by Freeman (1966) and Kalnajs (1971). In the context of the epicycle approximation, JR = ER/κ, where κ is the epicycle radial frequency and ER is the specific radial kinetic energy. This approximation is only valid if JR ≪ Lz or for near-circular orbits (Kalnajs 1971).
Additionally, JR 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 JR in a stellar disc excited by a constant non-axisymmetric perturbation of pattern speed Ωp, can be related to ΔLz as:
where Ω is the orbit and particle angular frequency. At co-rotation, Ω = Ωp so that there is no JR variation; whereas at Lindblad resonances with an m-armed 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 ΔJR/ΔLz = ∓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: JR is proportional to the square of the epicycle amplitude (Binney & Tremaine 2008). Therefore, variations in Lz that do not produce significant variations in JR 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 JR that are not correlated with variations in Lz via Eq. (1) are typical of blurring. When both Lz and JR 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 (JR ≈ 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 JR is estimated in our simulations (Sect. 2), we present the results on the Chirikov diffusion timescale of JR in Sect. 3. Section 4 contains a discussion on the stellar distribution function as a function of Lz and JR, 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 JR diffusion
2.1. Dynamical concepts
The Chirikov diffusion rate of JR 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 Dn(Lz) 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 Δt2 = 102 (in the unit of time of the system). We can also estimate a diffusion timescale TD by re-normalising D2 by
, where the time-average 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 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), JR reduces to
where . In the limit of very small radial excursions around the guiding centre (X = (Rmax − Rmin)→0), JR reduces even further to
, where κ, ER, 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:
It is convenient to time-average here over Δt2 = 102 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 JR, we never have access to an instantaneous value for the radial action.
No assumption that leads to the approximation JR = ER/κ 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 Mtot = 2 × 1011 M⊙. All 4 × 107 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 Mtot at the end of RunC. The other particles are named “never-escaped” in the following.
3. Diffusion timescale (TD(JR)) results
The normalised relative frequency distribution (akin to a probability density function) of TD(JR) is plotted in Fig. 1. For the sake of clarity, we have restricted this figure to the range 10−2 − 90 Gyr. Distributions for TD(JR) and TD(E) are overplotted for comparison purposes (Wozniak 2020). The distribution shapes are different, especially for TD ≲ 0.9 Gyr. On the contrary to E and Lz, there are few particles with TD(JR) < 100 Myr. The JR diffusion timescale is never very short. The maximum of TD(JR) is located at ≈0.9 Gyr, followed by a plateau up to ≈3 Gyr. Additionally, TD(JR) = 3 Gyr is also the median in mass or number of particles. On both sides, the frequency of TD(JR) decreases sharply. In other words, the characteristic diffusion timescale of JR is remarkably of the same order of magnitude as dynamical timescales in the galactic disc. This result contrasts with what was obtained for Lz and E, for which TD distributions are dominated by short timescales, outside the local maxima at TD(E) ∼ 10 and TD(Lz) ∼ 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 TD(JR) averaged over sets of particles (designated by ⟨TD(JR)⟩) 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 Lz 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 (UHRB at t = 3.16 Gyr) and the outermost bar’s outer Lindblad resonance (OLRB at t = 10.54 Gyr), whereas the range covered by the co-rotation (CRB) is approximately represented by the shaded area. Again, for the sake of comparison, ⟨TD(E)⟩ and ⟨TD(Lz)⟩ 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 |
As shown in Wozniak (2020), ⟨TD(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 ⟨TD(Lz)⟩ increases monotonically to values that can be considered as slow-diffusion. This is not the case for ⟨TD(JR)⟩, which decreases regularly from values similar to ⟨TD(Lz)⟩ in the centre (∼10 Gyr) to values around 1 Gyr in the disc. Beyond UHRB, ⟨TD(JR)⟩ 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, ⟨TD(JR)⟩ is always smaller than ⟨TD(Lz)⟩.
Particles that are retrograde, on average (0.1 Mtot, typical of a barred galaxy, and an identical proportion among never-escaped particles), have shorter diffusion times for Lz and JR than for E. The values remain compatible with the timescales in the disc for Lz and JR. 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 TD(Lz)−TD(JR) 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 TD(JR) axis, we can recover the peak and plateau identified in Fig. 1. The most remarkable structure extends along the bisector TD(JR) = TD(Lz). 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 TD(Lz) axis, this large structure explains both the plateau and the quasi-linear decay observed in Fig. 1 when TD(JR) > 3 Gyr. The other noticeable structure with the same density and long TD(Lz) (> 10 Gyr) is related to the disc. It contributes to a large extent to the peak at TD(Lz)) = 0.9 Gyr, although the two structures are blended by integrating along the TD(Lz) 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 TD(Lz) > 10 Gyr (as in Wozniak 2020) and TD(JR) < 3 Gyr (roughly the end of the plateau in Fig. 1). Other thresholds have been tested, but these ones roughly isolate particles with “long” TD(Lz) and “short” TD(JR). Mass fractions are calculated with respect to the total mass of RunC: 0.65 Mtot for Subset A, 0.24 Mtot 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 Mtot), 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 UHRB, 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 TD(Lz) > 10 Gyr and TD(JR) < 3 Gyr, the mass distribution (0.24 Mtot) 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 OLRB. 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 Lz and JR
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 Lz and JR 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 |
Subset B is identifiable through several substructures. A density peak is present for (circular orbits) and
(beyond OLRB 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 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 L1, 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 ER. 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 UHRB. A component separation based on TD 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(JR) or DF(Lz) can be recovered by the projection of the density map on the axes. We note that DF(Lz) 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(Lz), 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(Lz) to the maximum DF(Lz = 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 dotted-dashed lines), “never-escaped” particles (dotted lines), and the two subpopulations defined in Fig. 4 (red and black lines) as a function of |
Other evidence is that the number of particles with Lz = 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 Lz 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 (JR = 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 JR 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 Mtot (i.e. 7.7 × 109 M⊙). Almost all particles are located in the
tail in the outermost part of the disc (
), and they are mostly well beyond OLRB (located at 20.8 kpc, i.e. Lz ≈ 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 |
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−1 beyond OLRB. 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 Lz exchange, coupled with the increase in JR, presumably in the form of waves. Well beyond OLRB, 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 ΔJR ≠ 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 OLRB, with the notable exception of those inside the bar, which are concentrated around the Lagrange L4, 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 JR 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 OLRB and do not exceed the co-rotation of the outermost waves.
![]() |
Fig. 8. Projected mass surface density for particles selected as |
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−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 co-rotation (CRB), 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−1 kpc−1. This mode has a spatial extension that goes roughly from UHRB 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 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), OLRB and CRiS are close to each other, as well as UHRB−ILRiS. Between t = 8.38 and 10.54 Gyr (Fig. 9 right panel), UHRB is still close to ILRiS, while OLRB and CRiS are now separated by ≈1.5 kpc. Within the large ΩoW uncertainties, ILRoW 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 UHRB−ILRiS overlap is constant within less than 0.5 kpc, while OLRB−CRiS 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 OLRB might be a barrier to radial migration (Wozniak 2020). Indeed, OLRB 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 UHRoW of the lowest frequency outer waves (Fig. 8). The UHRoW 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 JR than for Lz in the disc. On average, TD(JR) is even slightly shorter in the disc than in the bar. The set of Subset B particles, selected such that TD(JR) < 3 and TD(Lz) > 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, JR 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 JR 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 Lz 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 TD(JR) < 3 and TD(Lz) < 10.54 (lower-left corner of Fig. 3). Thus, a fraction of the disc shows Lz 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(JR) and DF(Lz) widen by ≈20% (according to their full width at half maximum), while their maxima decrease. The widening of DF(Lz) for Subset A is mainly due to an increase in Lz for particles initially with Lz ≈ 0. This coevolution of JR and Lz 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 ΔJR ≠ 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 Lz ≈ 4000 and 4500, and JR ≈ 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 (JR = 0), accompanied by a redistribution in Lz (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 OLRB−CRiS and UHRB−ILRiS or +1/4 iS−ILRoW, 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 CRiS−ILRoW 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 −ILRoW, 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 Lz (and E) from the wave, while those at any m < 0 give Lz (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 JR that are reached are the signature of the hot population. This population should not see a significant variation of JR since it is supposed to be scattered by the bar co-rotation. However, the approximate calculation of JR (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 JR to calculate the Chirikov diffusion rate and the associated characteristic timescale, in addition to the results already obtained so far for Lz and E by Wozniak (2020), we have shown the following.
-
The distribution of JR diffusion timescales spikes around TD(JR) = 0.9 Gyr, which is mainly due to disc particles. It is followed by a plateau ending at TD(JR) ≈ 3 Gyr, whose main contributor is the stellar bar. Roughly 0.5 Mtot has TD(JR) < 3 Gyr and 0.77 Mtot has TD(JR) < 10.54 Gyr (i.e. the simulation length). Beyond the bar UHR, the space averaged timescale is ⟨TD(JR)⟩ ∼ 1 Gyr.
-
By selecting particles as TD(JR) < 3 Gyr and TD(Lz) > 10.54 Gyr, that is, 0.25 Mtot, 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 de-population of circular orbits. We note that 57% of particles on circular orbits (JR = 0) at t = 3.16 Gyr have JR increased by a few tens of kpc km s−1 after 7 Gyr, leading σR to increase by ∼10 km s−1. This de-circularisation is accompanied by Lz transfer through coherent wave-particle 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 non-permanent intermediate spiral, which could be a beating phenomenon, with (3) a set of intermittent multi-armed wave packets that carry Lz towards the edge of the disc. Several resonance overlaps ensure a continuous cascade of waves, which can also contribute to promote Lz 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 D2 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]
- Lynden-Bell, 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]
- Michel-Dansac, 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 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 |
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 |
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 |
In the text |
![]() |
Fig. 7. Particle number as a function of |
In the text |
![]() |
Fig. 8. Projected mass surface density for particles selected as |
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 |
Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.