Radial migration in a stellar galactic disc with thick components

We study how migration affects stars of a galaxy with a thin stellar disc and thicker stellar components. The simulated galaxy has a strong bar and lasting spiral arms. We find that the amplitude of the churning (change in angular momentum) is similar for thin and thick components, and of limited amplitude, and that stars of all components can be trapped at the corotation of the bar. At the exception of those stars trapped at the corotation, we find that stars that are far from their initial guiding radius are more likely so due to blurring rather than churning effects. We compare the simulation to orbits integration with a fixed gravitational potential rotating at a constant speed. In the latter case, stars trapped at corotation are churned periodically outside and inside the corotation radius, with a zero net average. However, as the bar speed of the simulated galaxy decreases and its corotation radius increases, stars trapped at corotation for several Gyrs can be churned outwards on average. We study the location of extreme migrators (stars experimenting the largest churning) and find that extreme migrators come from regions on the leading side of the effective potential local maxima.


Introduction
Observations and theoretical and/or numerical studies have shown that radial migration is of interest to understand some metallicity and chemistry observations such as the agemetallicty scatter in the solar neighbourhood (e.g. Haywood 2008;Schönrich & Binney 2009), stellar metallicity distributions at different radii in the Milky Way (Loebman et al. 2016) or the upturn of mean stellar age in the outskirts of local galaxies (e.g. Roškar et al. 2008;Bakos et al. 2008;Ruiz-Lara et al. 2017).
The guiding radius of a star can change because of dynamical interactions during mergers, with clumpy structures in the discs or with non-axisymmetric patterns like bars or spiral arms (e.g. Lynden-Bell & Kalnajs 1972;Sellwood & Binney 2002). A number of studies have focused on radial migration generated by resonance with such non-axisymmetric patterns (e.g. Minchev & Famaey 2010;Brunetti et al. 2011;Minchev et al. 2011Minchev et al. , 2012bRoškar et al. 2012;Grand et al. 2012Grand et al. , 2014Di Matteo et al. 2013;Kubryk et al. 2013;Vera-Ciro et al. 2014;Halle et al. 2015;Daniel & Wyse 2015). This change of guiding radius is often named churning as opposed to blurring that is due to epicyclic oscillations around the guiding radius (Schönrich & Binney 2009). The churning can be oscillatory in the case of an interaction with a lasting non-axisymmetric pattern such as long-lived spiral arms or long-lived bars (Sellwood & Binney 2002;Ceverino & Klypin 2007;Binney & Tremaine 2008). This has led to consider that bars, that seem to be long-lived (even if they can be destroyed by gas infall (Bournaud & Combes 2002) and grow again), could not drive substantial radial migration (e.g. Aumer et al. 2016) because stars are periodically churned back and forth in a region around corotation. In contrast, stars corotating with transient spiral arms are able to change their angular momentum permanently when transient spiral vanish (Sellwood & Binney 2002). However, the churning driven by a bar mixes significantly the stars around the corotation of the bar, and, if the bar has a non-constant speed, the portion of the disc affected by churning can shift towards the inside or outside of the disc, in case of an increase (e.g. Ceverino & Klypin 2007;Halle et al. 2015) or decrease of the bar-speed, respectively.
Radial migration has been studied in some simulations including thick discs (e.g. Solway et al. 2012;Aumer et al. 2017) or also in the context of its potential thickening effect on discs (e.g. Schönrich & Binney 2009;Loebman et al. 2011;Minchev et al. 2012a;Roškar et al. 2013;Vera-Ciro et al. 2014Kubryk et al. 2015;Grand et al. 2016;Schönrich & McMillan 2017). Some studies however suggest that, at least in the Milky Way, the characteristics of the thick disc do not seem to require any significant migration and can be well explained by formation from a disc rich in gas and turbulent (e.g. Noguchi 1998;Brook et al. 2004;Haywood et al. 2013Haywood et al. , 2015Lehnert et al. 2014). It has been argued that stars migrating from the inner disc parts to outer regions should heat the outer regions because they come from higher velocity dispersion regions (Schönrich & Binney 2009;Roškar et al. 2013), but this may be balanced by the provenance bias of the migrating stars (Vera-Ciro et al. 2014; see also Minchev et al. 2012a;: the stars that are more likely to migrate are the stars remaining close to the disc mid-plane, hence stars with a globally low velocity dispersion. Solway et al. (2012) found the churning in a thick disc is only mildly more important than the churning of thin disc stars.
In this work, we study how stars that remain trapped around the bar corotation can be globally churned outwards in the case of a slowing-down bar whose corotation radius increases with time, in a disc galaxy with thick components. We use an N-body Thin disc 2.6 × 10 10 4.7 0.3 1 × 10 7 Intermediate disc 1.5 × 10 10 2.3 0.6 6 × 10 6 Thick disc 1.0 × 10 10 2.3 0.9 4 × 10 6 Dark matter halo 1.6 × 10 11 10 5 × 10 6 simulation of an isolated disc galaxy with three disc components of different scale heights, embedded in a live dark-matter halo, and compare the radial migration to the case of a fixed gravitational potential rotating at a constant speed. Section 2 presents the numerical simulation used in this work, its initial conditions and the dynamical evolution of the stellar disc and its non-axisymmetric patterns. Section 3 focuses on radial migration of the thin and thicker disc components with a comparison to the radial migration in the fixed potential rotating at a constant speed case. Section 4 focuses on the extreme migrators churned outwards, comparing their orbits and spatial location to the fixed potential rotating at a constant speed case.

Initial conditions
The simulated galaxy contains a stellar disc with three components of different vertical scale heights, and a dark matter halo. No stellar bulge is included. Masses of the different components, number of particles and parameters are shown in Table 1.
The discs are called thin, intermediate, and thick in reference to their relative average thickness. They have Miyamoto-Nagai density profiles: with height parameters h thin < h inter < h thick and radial parameters a thin and lower a inter = a thick , as detailed in Table 1. The three disc components of different scale heights aim at representing a more realistic total disc than a case with only a thin and a thick component (this is suggested by some studies of stellar populations in the Milky Way indicating a continuous variation of disc properties with scale height; e.g. Bovy et al. 2012). In the following analyses, however, we often show results for the three individual components because the intermediate case is an interesting transition between the thinner and thicker components for the study of radial migration. The dark matter halo has a Plummer profile: The initial conditions are set with an iterative method allowing the components to be dynamically relaxed. Radial migration due to resonance with non-axisymmetric patterns in isolated discs may sometimes be overestimated because of radial expansion of the stellar discs from unrelaxed initial conditions, which is avoided in this work. The low increase in radial extent of the discs is seen on Fig. 1

Dynamical evolution
The simulated time span is of 5 Gyrs. We use the Tree code of Semelin & Combes (2002) with a softening length set to 50 pc. Figure 2 shows maps of the surface density of the whole disc and its three components after 1, 2, 3, 4 and 5 Gyr of evolution. Overplotted in white are the positive isocontours of Σ(R, θ) − Σ(R) Σ(R) , the relative deviation of surface density at radius R and azimuth θ, Σ(R, θ), from the azimuthally averaged value Σ(R). The disc develops a bar and spiral features. The bar is the strongest non-axisymmetric perturbation at almost all times. Spiral arms have the same angular speed as the bar. Figure 3 shows the time evolution of the bar strength as estimated from a Fourier decomposition of the surface density (from stars at |z| < 500 pc) of either each individual disc component or the whole components. We note that this simulation does not aim to reproduce the bar properties of the Milky Way. The bar in the thick disc is significantly weaker (as seen in e.g. Combes et al. 1990;Athanassoula 2003;Fragkoudi et al. 2017). At initial times, the bar of the intermediate disc is slightly stronger than the thin disc bar because of initial heating of the thin disc that makes it more stable than the intermediate disc, until 0.5 Gyr. The bar angular speed decreases with time by angular momentum transfer from the disc to the dark matter halo as seen on Fig. 3 (e.g. Debattista & Sellwood 1998;Athanassoula 2002;Di Matteo et al. 2014). The largest angular momentum transfer occurs from the thin disc to the halo (as in e.g. Fragkoudi et al. 2017). This slowing-down allows the bar to radially extend. It develops a buckling instability associated with decreases in the bar strength visible on Fig. 3 at t = 1.9 and 3 Gyr. This buckling can be seen in the edge-on views of Fig. 2. The bar speed is estimated by a Fourier method as in Halle et al. (2015) and resonances radii (radii at which stars on nearly circular orbits resonate with the bar) are obtained by finding the radii where Ω(R) − Ω p and κ(R) are commensurable, with Ω(R) and κ(R) the angular speed and epicyclic frequency (respectively) of particles     Fig. 4. Left first row: determination of resonances radii at different times from a Fourier analysis on the whole disc. Horizontal lines show the pattern speed Ω p obtained from the maximum of the spectrogram, the solid curve is Ω(R), dashed curve is Ω(R) + κ 2 , and dash-dotted curve is Ω(R) − κ 2 . Vertical lines are the radii of the ILR (blue), corotation (orange), and OLR (green). Left second row: strength of π-periodic non axisymmetries of the whole disc as a function of radius at the same times as in the first row. Right plot: resulting bar speed and resonance radii estimation as a function of time. on almost circular orbits at R. Figure 4 shows the determination of the radii of the inner Lindblad resonance (ILR) where Ω(R) − Ω p = κ 2 , corotation where Ω(R) − Ω p = 0 and of the outer Lindblad resonance (OLR) Ω − Ω p = − κ 2 at different times, and their time evolution.

Blurring and churning
Radial migration is often separated into blurring, due to epicyclic oscillations around a guiding radius with time, and churning that is the change of the guiding radius (Sellwood & Binney 2002;Schönrich & Binney 2009). We compute average radii (hereafter referred to as guiding radii) as in Halle et al. (2015), by a local average of the radial oscillations of a stellar particle using the closest local radial minima and maxima. Figure 5 shows the distributions of the changes in galactocentric radius and guiding radius as a function of the initial values, for different time intervals. The black and red curves show the average and the dispersion (respectively) as a function of initial radius or guiding radius, and the total rms value is indicated in each panel. During these time intervals, the resonance radii grow (see Fig. 4), and the corresponding encompassed radial ranges are shown in the shaded areas in Fig. 5. A mean corotation radius (average of the values of the corotation radius at the beginning and end of the time interval) is represented by a thicker vertical line, together with a line of slope −1 intersecting the x-axis at the mean corotation radius, allowing to select stars crossing the mean corotation outwards or inwards, and a line of slope −2 with the same x-axis intersection at the mean corotation radius, along which stars exchange their position with respect to the mean corotation radius (as in Sellwood & Binney 2002). Some diagonal features are visible in the distributions, they are located around some resonances of the bar such as corotation, with outwards migration of stars below the resonance radius and inwards migration of stars beyond the resonance radius. For each pair of panels corresponding to the same time interval, the global dispersion and maximum change is lower for guiding radius because the epicyclic oscillations around the guiding radii are removed.
As in many other studies (e.g. Sellwood & Binney 2002), we see the main churning occurs around corotation, with other patterns located at the ILR or, for example, the ultra-harmonic resonance Ω − Ω p = κ 4 at R 10 kpc from t = 2 to t = 3 Gyr (only the ILR, corotation and OLR radii are shown on the panels). the bar, implying a larger and larger corotation radius, allows for a large part of the disc to be affected by migration around an outwards-shifting corotation.
Comparison of panels of change in guiding radius during the same duration also show that the amplitude of churning at corotation depends on the strength of the bar as in Halle et al. (2015): it is, for example, stronger from t = 1 to t = 2 Gyr than from t = 3 to t = 4 Gyr (global rms value of 1 kpc and 0.6 kpc, respectively, and churning signal at corotation of a lower amplitude), as the bar is weaker at late times (see Fig. 3).
The top right panel of Fig. 5 shows the change in guiding radius of thin disc stars on a time interval of 3 Gyr (from t = 1 to t = 4 Gyr). It can be seen that stars with initial guiding radii close to the initial corotation radius are the most extreme outwards migrators (their guiding radius can increase by as much as almost 10 kpc). These extreme migrators are studied in greater detail in Sect. 4.

Migration in thick components
Stars in the thick components have larger radial velocity dispersions than thin disc stars by construction (see Fig. A.2) and thus larger radial excursions. The changes in galactocentric radius in a given time interval as from t = 1 to t = 2 Gyr represented on the two top-right panels of Fig. 6 can thus be significantly larger than for the thin disc (top left panel of Fig. 5). The largest migration in terms of galactocentric radius occurs at the OLR (for relatively few stars as this resonance is located near the outermost parts of the intermediate and thick discs). This migration signal at the OLR corresponds to stars migrating of as much as 10 kpc in the intermediate disc, and as much as 20 kpc in the thick disc.
However, comparison of the bottom panels of Fig. 6 to the corresponding plots on this time interval in Fig. 5 shows that the changes in guiding radius are similar, both in their amplitude and root mean square value (global ones or as a function of initial guiding radius (red curves)). This is consistent with results of (Solway et al. 2012) who find similar changes in vertical angular momentum for thin and thick discs of simulated galactic discs with a bar. The relatively weak churning at OLR implies that the large migration in terms of galactocentric radius is due to the shape of orbits of particles at the OLR that allows them to have very large radial excursions.
The higher radial velocity dispersion of the thick component is also associated with a larger asymmetric drift effect, implying that stars resonating with the bar are on average at a lower radius in the thick components than in the thin disc. The migration signal around corotation is thus slightly shifted to lower radii or guiding radii for thick components as can be seen, for example, on the bottom-right panel of Fig. 6 compared to the corresponding panel of Fig. 5.
As for the thin disc, during the 3 Gyrs time interval from t = 1 to t = 4 Gyr, some stars with an initial guiding radius close to the initial corotation radius can be churned outwards by almost 10 kpc, which is also discussed in more detail in Sect. 4. At the exception of those stars, especially in the thick components, stars that are far from their initial guiding radius are more likely to be so due to blurring rather than churning effects.

Comparison to a fixed potential
To compare the simulation to a case in which the bar has a constant strength and length, we took the gravitational potential of the simulation at a time t 0 and integrated orbits by simply rotating this fixed gravitational potential of an angle Ω P,t 0 (t − t 0 ) at each time t, with Ω P,t 0 the pattern speed at t 0 , determined from the previous Fourier analysis of Sect. 2. We used a Kick-Drift-Kick time integration with a timestep of 1 Myr. Orbits were integrated from t 0 = 1 to t = 4 Gyr and we computed their radius and guiding radius as a function of time as for the simulation particles. Figure 7 shows a comparison of radial migration in the simulation and for the orbits integration for the different disc components and the whole disc from t = 1 to t = 2 Gyr. During this time interval, the corotation radius grows by a few Fig. 6. First row: distributions of the change in radius R f − R i as a function of R i from t = 1 to t = 2 Gyr (left panels) or from t = 1 to t = 4 Gyr (right panels) in the intermediate and thick discs. Second row: distributions of the change in guiding radius R g f − R g i as a function of R g i from t = 1 to t = 2 Gyr (left panels) or from t = 1 to t = 4 Gyr (right panels) in the intermediate and thick discs. kpcs (from 10 kpc to 13 kpc), which broadens the churning signal at corotation in the simulation in comparison to the orbits integration. More stars are churned in the shifting corotation region, which makes the global rms values of churning larger in all components. The amplitude of churning is only mildly larger in the simulation than for orbits integration.
Differences between the simulation and the integrated orbits in the fixed potential are more visible on a longer time interval. Figure 8 shows a comparison between the simulation and orbits integration for a time interval from t = 1 to t = 3 Gyr. Unlike the previous time interval of 1 Gyr, for 3 Gyr a significant difference can be seen for both changes in galactocentric radii and guiding radii around corotation. Stars can be churned outwards by almost twice as much a distance in the simulation compared to the orbits integration, for all disc components.
In the case of a fixed pattern rotating at a fixed speed, stars can remain trapped at corotation in a libration movement around the local maxima of the effective potential φ eff = φ g − 1/2Ω 2 p R 2 , where φ g is the gravitational potential and Ω p the pattern speed (see Sect. 3.3 of Binney &Tremaine 2008 andBinney 2002). The time evolution of the galactocentric radius of stars trapped at corotation consists of oscillations around a more slowly oscillating guiding radius. Here, the changes in guiding radius in the orbits integration are fairly similar in both considered time intervals, both in amplitude and rms values because of these churning oscillations of trapped stars (the first time interval is already similar to or larger than the average half-period of oscillations of guiding radius around corotation). For an evolving pattern-speed, the situation is more complicated than periodic churning as stars can remain trapped, but also be liberated along with the shifting of the resonance location, and new stars can become trapped. In the next section we focus on the extreme migrators, stars churned outwards the most.

Stars at corotation
We first looked for stellar particles corotating with the bar at t = 1 Gyr, so as to study their migration. Determining which stars resonate with a non-axisymmetric pattern is possible by extracting their individual angular Ω and radial κ frequencies from a Fourier analysis of their orbits (see Binney & Spergel 1982;Athanassoula 2002;Ceverino & Klypin 2007). We analysed orbits integrated with the simulation gravitational potential of t = 1 Gyr rotating at the t = 1 Gyr bar speed as described in Sect. 3.3, to obtain stars resonating with the bar at t = 1 Gyr. We computed the radial frequency κ from a Fourier transform of the time evolution of the radius, removing the (usually) low frequency of the angular momentum oscillations, and estimate the angular frequency Ω by a Fourier transform of the Cartesian x or y coordinates divided by the radius. Figure 9 shows the distribution of the ratio Ω−Ω p κ , with Ω p the bar speed at t = 1 Gyr. Several peaks can be seen, the most prominent ones being at the ILR ( Ω−Ω p κ = 0.5) and corotation ( Ω−Ω p κ = 0) in all disc components. There is also a complex structure of peaks between corotation and ILR concerning a large fraction of the stars, especially in the thickest component for which a high and broad peak is visible. This is likely due to the beginning of the buckling instability and is similar to the patterns obtained in Martinez-Valpuesta et al. (2006). The OLR peak ( Ω−Ω p κ = −0.5) is also visible in all three components. The surface density declines approximately exponentially with radius and the bar is strong, hence the larger fraction of stars at the ILR peak for all three components. Figure 10 shows in colour the density of particles found at corotation by the Fourier analysis, that is, the particles of the Fig. 7. First and second rows: distributions of the change in radius R f − R i as a function of R i in the simulation (first row) and for the integration of orbits (second row) from t = 1 to t = 2 Gyr. Third and fourth rows: distributions of the change in guiding radius R g f − R g i as a function of R g i in the simulation (third row) and for the orbits integration (fourth row) from t = 1 to t = 2 Gyr. Ω − Ω p κ = 0 peak (making up 7.6% of the thin disc mass, 6% of the intermediate disc mass and 7.8% of the thick disc mass). The black contours represent some isocontours of the effective potential in the frame rotating at the bar speed. The purple contours are positive isocontours of the overdensity of the whole disc as in the left column of Fig. 2. The bar is aligned with the x-axis. The local maxima of the effective potential (encompassed by closed curves in the top right and bottom left) are not on a line orthogonal to the bar as Lagrange points L4 and L5 of a simple barred potential, because the gravitational potential is gradually tilted by the spiral arms contribution as radius increases. Thin disc stars at corotation are localised around those local maxima. The distribution of stars of thicker components is wider in radius because of the higher eccentricities of their orbits. Figures B.1 and B.2 show the density maps of stars at the ILR and at the OLR (respectively).

4.2.
Orbits of stars at corotation for a fixed potential rotating at a constant speed Figure 11 shows an example orbit of thin-disc stellar particles corotating with the bar: the orbit is represented by its x-and y-coordinates, and by x rot and y rot , its coordinates in a frame rotating around the z-axis at the pattern speed Ω p . In both latter plots, the bar is parallel to the x-axis. In the rotating frame, the libration movement is clearly visible (the amplitude of the azimuthal excursions in the rotating frame vary from star to star). The maximum value of the churning amplitude is expected to increase with the strength of the non-axisymmetric patterns (Sellwood & Binney 2002;Binney & Tremaine 2008). The stellar particle rotates either faster than the bar at a lower radius or slower at a larger radius (and on average at the bar speed). Its guiding radius is plotted with a dashed line, exhibiting the churning oscillations happening at corotation. The evolution of the radius as a function of time (bottom left panel) also shows epicyclic oscillations modulated by a lower frequency guiding radius oscillation around the corotation radius. Stars of thicker components can also have librating orbits as see in Fig. 12. Their radial "epicycle" excursions (the orbits are far from close to circular orbits) have an amplitude similar to the global radial motion amplitude. Because of the asymmetric drift effect discussed in Sect. 3.2, corotating stars are on average at lower radii than thin disc corotating stars, and have thus an average lower guiding radius and vertical angular momentum. Their average position in the rotating frame can also be closer to orthogonal to the bar in the gravitational potential of Third and fourth rows: distributions of the change in guiding radius R g f − R g i as a function of R g i in the simulation (third row) and for the orbits integration (fourth row) from t = 1 to t = 4 Gyr. this orbits integration, because of the decreasing spiral tilting of the potential as radius decreases.

Migration of stars at corotation at an initial time
We now compare the fate of corotating stars at t = 1 Gyr in the orbits integration case (with a constant bar speed) to the simulation case with a decreasing bar speed. Figure 13 shows the distribution of these stars at t = 2 and t = 4 Gyr in the orbits integration and in the simulation. Density maps are each rotated so that the bar is parallel to the x-axis.
In the orbits integration, the distribution remains almost the same: stars just evolve on their orbits filling the same spatial area.   Fig. 11 for a thick disc stellar particle.
In the simulation, the distribution is however torn out into a spiralling shape, with stars reaching higher radii, similar to the corotation radius at the end of the time interval for the ones migrating the most (or even a few kpcs higher in the thick disc), and some stars remaining close to the 1 Gyr corotation radius. A fraction of the stars thus remains trapped at corotation during the time interval and migrates outwards, while the rest of the stars remain close to their initial radius. We determined which stars corotate with the bar at t = 2 Gyr and t = 4 Gyr by the same method as in Sects. 3.3 and 4.1, and find that the fraction of stars at corotation at t = 1 Gyr remaining trapped at t = 2 or t = 4 Gyr decreases with the thickness of the disc component. 19% of the thin disc stars corotating with the bar at t = 1 Gyr still corotate with the bar at t = 2 Gyr, while this fraction is of 16% for the intermediate disc and 13% for the thick disc. Between t = 2 and t = 4 Gyr, some stars do not remain in corotation, only 4.8% of the thin disc stars corotating with the bar at t = 1 Gyr still corotate with the bar at t = 4 Gyr, while this fraction is of 3.9% for the intermediate disc and 1.2% for the thick disc. Figure 14 shows the location of the disc stars whose guiding radii increase the most in the orbits integration (the top 1% of the distribution of change in guiding radius for each disc component) between t = 1 Gyr and t = 2 Gyr (top three rows) and between t = 1 and t = 4 Gyr (bottom three rows). Density maps are shown at t = 1 Gyr, and at the final time (t = 2 or t = 4 Gyr) in the orbits integrations case and in the simulation case. At t = 1 Gyr, the thin disc stars (with low eccentricities and therefore galactocentric radii close to their guiding radii) consist, as expected, of corotating stars in the part of their orbits with a radius lower than corotation. At the final time (t = 2 or t = 4 Gyr), they are beyond the (fixed) corotation radius in the orbits integration. Stars of the thicker components are distributed in a wider radial range but their average radius is below the corotation radius at t = 1 Gyr, and beyond it at the final time. In the simulation, those stars do not necessarily migrate much: as discussed in the last section, some of them remain at their initial radius while a fraction is churned outwards. Only a very small fraction of stars migrating the most in the orbits integration from t = 1 to t = 4 Gyr reaches radii close to the corotation radius at t = 4 Gyr in the simulation, indicating the extreme migrators in the simulation in this time interval must consist of different stars. Figure 15 shows the location of the disc stars whose guiding radii increase the most in the simulation (the top 1% of the distribution of change in guiding radius for each disc A86, page 9 of 14  . Stars corotating with the bar at t = 1 Gyr shown at t = 1 Gyr (initial time for the orbits integration) and at t = 2 and t = 4 Gyr in the orbits integration case and simulation. Maps are each rotated so that the bar is parallel to the x-axis. The red circle is the corotation radius at t = 1 Gyr. The orange circle is the corotation radius at t = 2 or t = 4 Gyr. Black contours are isocontours encompassing the local maxima of the effective potential with the pattern speed at t = 1 Gyr for the orbits integration panels, or at t = 2 or t = 4 Gyr for the simulation panels. component) between t = 1 Gyr and t = 2 Gyr (top three rows) and between t = 1 Gyr and t = 4 Gyr (bottom three rows). As in Fig. 14, density maps are shown at t = 1 Gyr, and at the final time (t = 2 or t = 4 Gyr) in the orbits integration case and in the simulation case. The spatial distributions at t = 1 Gyr show that particles that remain trapped at corotation (and are thus significantly churned outwards) are located beyond the local effective potential maxima (in the counter-clockwise rotation sense). Some of these stars are corotating with the bar at t = 1 Gyr, but the distribution is deprived of stars behind the effective potential local maxima (in the counter-clockwise rotation sense). The latter stars must be liberated from trapping at corotation, they do not follow the decrease in average angular frequency of the stars trapped at the corotation of a slowing-down non-axisymmetric pattern. The spatial distribution of extreme migrators is even more tilted towards higher azimuths (rotating counter-clockwise) when looking at extreme migrators in the time interval from t = 1 to t = 4 Gyr. Some extreme migrators in the latter time interval do not corotate with the bar at t = 1 Gyr (as can be deduced from the spatial distribution of the orbits integration case at t = 4 Gyr, departing from the spatial distribution of particles at corotation of Fig. 13). These stars are trapped by the corotation at a later time (after t = 1 Gyr).

Extreme migrators
Finally, we compare the orbits of extreme migrators in the simulation to the orbits of the same stars in the orbits integration. Figure 16 shows such a comparison for a thin disc star and Fig. 17 for a thick disc star. Both stars corotate with the bar at t = 1 Gyr as can be seen in the panel showing the libration of the integrated orbits (thin lines) in the frame rotating at the bar speed at t = 1 Gyr. As they are trapped at the corotation of the slowing-down pattern, their angular speed decreases and the orbits of the simulation circulate around the centre of the galaxy in this frame. The thick disc star has a high eccentricity but it is however trapped at the bar corotation just as thin disc stars, as found by Binney (2018), and has a larger and larger guiding radius.

Conclusion
We have studied radial migration in a galactic disc with thick components with a bar and corotating spiral arms as the most prominent non-axisymmetry. The dark matter halo is massive and concentrated, so as to enhance the shifting of resonances radii of interest in this work.
Stars of the thick components can have large radial excursions because of the high eccentricities of their orbits but the churning is limited and very similar to the churning in the thin disc. This is consistent with results of Solway et al. (2012) that found thick disc stars undergo almost the same churning as thin disc stars. We note that during the dynamical evolution, the thin disc -initially less stable than the thicker components -heats up, which reduces the difference between the disc components in velocity dispersion.
We have shown that stars belonging to thick components can be trapped at the bar corotation, as found by Binney (2018). If the bar keeps the same strength and speed, trapped stars are churned periodically outwards and inwards at corotation, but if the corotation radius increases (in the case of a slowing-down bar), they can be churned to larger guiding radii on average. This outwards churning is possible even in disc components of decreasing total angular momentum because the loss of angular  Stars churned outwards the most in the orbits integration (top 1%) from t = 1 to t = 2 Gyr (top three rows) or from t = 1 to t = 4 Gyr (bottom three rows). Left column: spatial distribution of these extreme migrators at t = 1 Gyr. Middle and right columns: spatial distribution at t = 2 or t = 4 Gyr in the orbits integration case and in the simulation. Each map is rotated such that the bar is parallel to the x-axis. The red circle is the corotation radius at t = 1 Gyr. The orange circle is the corotation radius at t = 2 or t = 4 Gyr. Black contours are the isocontours encompassing the local maxima of the effective potential with the pattern speed at t = 1 Gyr for the orbits integration panels, or at t = 2 or t = 4 Gyr for the simulation panels.  Extreme migrators constitute only a fraction of stars at corotation at an initial time, this fraction decreasing as time goes by and being a little smaller for thicker components. These stars can be distinguished from their spatial location with respect to the maxima of the local effective potential extrema at initial time.  Fig. 11 for a thin disc stellar particle evolution in the orbits integration (first row and thin lines of third row) and in the simulation (second row and thick lines of third row) in which the stellar particle is an extreme migrator. The black dash-dotted line shows the time evolution of the corotation radius (bottom left panel) and of Ω p (bottom right panel) in the simulations.   Fig. 11 for a thick disc stellar particle evolution in the orbits integration (first row and thin lines of third row) and in the simulation (second row and thick lines of third row) in which the stellar particle is an extreme migrator. The lack dash-dotted line shows time evolution of the corotation radius (bottom left panel) and of Ω p (bottom right panel) in the simulations.
Potential observational signatures (in kinematics or metallicity and/or chemistry) of these extreme migrators of thin and thick components will be studied in a following paper.  In this appendix, we show initial condition details in addition to Sect. 2.1, as well as the evolution of some quantities: rotation curves and radial and vertical velocity dispersions of the disc components. Figure A.1 shows the contributions of the different disc components, the whole disc, and the dark matter halo in the initial conditions and at later times. The disc dominates the gravitational potential at low radii at all times, while the dark matter halo dominates at high radii. Figure A.2 shows the radial and vertical velocity dispersions of the different disc components in the initial conditions and at later times.

Appendix B: Stars at ILR and OLR
In this appendix, we show stars found at the ILR and the OLR for the different disc components by the analysis of Sect. 4.1. Figure B.1 shows the density maps of stars at the ILR in the different disc components. They mostly belong to the bar. Figure B.2 shows the density maps of stars at the OLR. While the distribution is an annulus for the thin disc stars, it has a wider radial extent for the thick components.