Positive Lynden-Bell derivative as a ticket to the bar trap?

We have translated the results of $N$-body simulations of one barred model into the language of action variables and frequencies. Using this language, we analysed the behaviour of all orbits in the model on a large time scale at the stage of a mature bar. We show that the orbits join the bar while preserving their adiabatic invariant, which takes into account the 3D structure of the orbits. This allows us to apply the concept of the Lynden-Bell derivative for each of these orbits and trace how the sign of the derivative changes, i.e. how asynchronous changes in angular momentum $L_z$ and orbital precession rate $\Omega_\mathrm{pr}$ (normal orbital mode) change to synchronous (abnormal mode). The transition to the abnormal mode occurs when $\Omega_\mathrm{pr}$ reaches the angular velocity of the pattern $\Omega_\mathrm{p}$, after which the orbit becomes stuck in the bar trap. All this happens against the background of secular changes in actions ($L_z$ decreases, $J_\mathrm{R}$ and $J_z$ increase). At the same time, corotation particles near two stable Lagrange points are also subject to secular changes in their actions. They increase $L_z$ and drift to the periphery, shifting corotation outwards. We also show that a change in the orbital mode from normal to abnormal and the trapping of orbits in a bar is possible only when the bar speed decreases with time, regardless of what is causing the bar to slow down. Our findings clarify and expand the picture of bar formation and evolution in numerical models.


INTRODUCTION
The question of the formation of bars in galaxies is usually reduced to the analysis of a bar-forming mode instability.But how exactly the bars form in real galaxies is an unresolved problem in galactic dynamics (Sellwood 2013(Sellwood , 2014)).This is surprising in the era of large (all sky) and deep (high redshifts) surveys and a large suit of cosmological simulations that attempt to replicate the evolution of the observable Universe.However, the "quantity" further confuses the issue, since barred galaxies are not much different from their almost equally numerous unbarred counterparts (see, for example, Courteau et al. 2003;Sellwood 2013).In particular, the location of a galaxy in the Tully-Fisher plane is independent of barredness.Moreover, barred and unbarred galaxies have, on average, comparable fractions of luminous and dark matter at a given radius.But even at the level of individual galaxies, the question turns out to be no less difficult.
One can guess what factors can stop the bar formation (see below).However, even for one of the typical spiral galaxies in the local Universe, namely M 33, there is still no explanation for why this particular galaxy is stable to the bar formation (Sellwood et al. 2019).All reasonable -body models of this galaxy, compatible with the mass distribution and those taking into account the morphological features of this galaxy, lead to the formation of a bar, which the real galaxy does not have (Sellwood et al. 2019).The inability to explain why some galaxies do not have bars actually means that we still do not understand all the details of the bar formation and do not see the whole picture of the galaxy secular evolution.At the same time, numerical simulations of the last years gradually expanded the list of factors that can stop or slow down the formation of bars.For example, Collier et al. (2018) found that a high halo spin  can prevent the bar formation.Nevertheless, a general principle of why a bar either forms or does not form in a particular simulation has yet to be understood.There is also an interesting discrepancy between the number of barred galaxies in the local Universe and cosmological up-to-date simulations (Reddish et al. 2022).
If one goes down to the level of specific mechanisms, two mechanisms for the bar formation in isolated galaxies are usually discussed, namely: the swing amplification (Toomre 1981) and the Lynden-Bell mechanism (Lynden-Bell 1979).
The swing amplification mechanism connects the formation and growth of the bar with a special type of interaction between the wave packets and the corotation (CR) region.At the corotation, outgoing leading spiral waves can be super-reflected into amplified ingoing trailing waves.The ingoing trailing wave, in turn, can be reflected from the disc centre, transforming into a leading wave that propagates outward until it reaches the corotation.The cycle of wave travelling and repetitive transformations creates a positive feedback loop.This, in turn, leads to an exponential increase in the amplitude of the perturbation in time.
The most valuable practical yields of the swing amplification theory are feedback-breaking factors.Specifically, Toomre (1981) believed that a dense bulge-like mass component, that could induce inner Lindblad resonance (ILR), might be the main stabilising factor for the cold massive disc.However, up-to-date -body simulations showed that the disc is stable against the formation of a bar only with a very dense and very massive central component (Fujii et al. 2018;Saha & Elmegreen 2018;Kataria & Das 2018).As argued by Polyachenko et al. (2016), compact, but less massive bulges do not prevent bar formation due to the finite thickness of the disc.In such discs, the precession profile when averaged over the height ("effective" precession profile) turns into a curve with a rather low maximum.Moreover, compact bulges, as well as the presence of gaseous central mass concentration, do not prevent the bars from appearing but lead to the formation of another subsystem, the socalled barlens (Athanassoula et al. 2015;Salo & Laurikainen 2017;Smirnov et al. 2021).The destruction of the bar itself does not occur in this case.Observations show that galaxies with barlenses are a large class of objects (Laurikainen et al. 2011).In their structure, they often exhibit compact classic bulges (Laurikainen et al. 2018).But, despite the presence of a compact bulge, such galaxies also have a long elongated bar.This shows that the wave picture drawn by Toomre (1981) is actually much more complex than originally thought.
The mechanism that was suggested by Lynden-Bell ( 1979) is similar to the radial orbits instability (ROI) arising in threedimensional systems (Polyachenko & Shukhman 1972, 1981;Antonov 1973).Merritt (1987Merritt ( , 1999)), Saha (1991), Weinberg (1991), Palmer (1994) interpret the ROI precisely in terms of the Lynden-Bell's mechanism.The Lynden-Bell's mechanism qualitatively describes how the planar eccentric orbits line up in the direction of a rather slowly rotating bar-like perturbation and gradually become trapped that leads to an increase of the potential perturbation.The key idea is that the eccentric orbits in the central regions of the disc precess slowly and usually satisfy the following condition: Ω pr /Ω ≪ 1, where Ω is an angular velocity and Ω pr is a precession rate.In this case, the orbital motion can be averaged over time, and one can consider not the individual stars moving along the corresponding orbits, but slowly precessing orbital ellipses.One can study how these ellipses interact with the gravitational field (and the bar "well" or "trap") then.Studying the normal modes of the stellar disc, Polyachenko & Polyachenko (2004); Polyachenko (2013) argued that, while the orbits themselves precess slowly, the capturing of such orbits in the bar trap can occur relatively fast, i.e. on a dynamical timescale (Sellwood 2014).
Assuming that |Ω pr − Ω p |/Ω ≪ 1 (where Ω p is the pattern speed of the bar), the evolution of the orbit under the influence of the bar-like potential proceeds adiabatically, that is, with the conservation of a specific integral of motion   =  R +   /2, where   is the vertical projection of the angular momentum,  R = ∫  R  is the radial action (the integration is carried out over the trajectory of a star).Lynden-Bell (1979) explicitly showed that, depending on the sign of the partial derivative (Ω pr /  )   (hereinafter, LB-derivative) calculated at the fixed value of the invariant   , the eccentric orbit will be trapped by the bar potential (the positive sign of the LB-derivative) or moves out of it (negative sign).The region of the galaxy, where the mentioned derivative has the negative sign, is called "normal" (usually the disc region), and the region, which has the negative sign of the derivative, is called "abnormal".
The formalism of Lynden-Bell (1979) was extended by Poly-achenko & Shukhman (2020a,b).These authors tried to connect the ideas of Lynden-Bell and the modern view of the bar as a structure consisting of various orbital families (Contopoulos 1980), like orbits elongated along the bar major axis (the so-called 1 family) or elongated along the bar minor axis (the 2 family).For several models of the potential (including the Milky Way-like potential), Polyachenko & Shukhman (2020a,b) studied stationary points of the Hamiltonian where the parenting orbits of the mentioned orbital families reside.Based on their findings, the authors concluded that, in addition to LB-derivative, two other parameters, the precession rate itself and "bar responsiveness" ( in their notation), determine whether a particular orbital family can support the bar structure or not.Given that the orbital dynamics of the bar is found to be rich in different families in recent years (Patsis & Athanassoula 2019;Smirnov et al. 2021), the studies by Polyachenko & Shukhman (2020a,b) highlight the importance of considering the individual orbits (bar "building blocks") to understand the process of bar formation/evolution and not the regions or areas of the disc as was done by Lynden-Bell (1979).
The orbital approach was partially adopted in Petersen et al. (2019), where the authors specifically divided the orbital types on 1 orbits and "other" bar orbits.In particular, Petersen et al. (2019) studied the evolution of angular momenta separately for the 1 orbits and other orbital types.Their findings indicate that the evolution indeed proceeds differently for different orbital types and, at the same time, bars in different models have different fractions of them (see also Smirnov et al. 2021).The change and evolution of orbital morphology and characteristics are very important during the deeply nonlinear stages of a bar's life.In this case, it is important not only to know about the presence of different types of orbits, but also about the quantitative relationship between them and how they transform when interacting with the bar and falling into its trap.At the same time, the Lynden-Bell language (action variables, the LBderivative), although it was initially applied to the early linear stages of bar formation, can also be useful for describing the later stages, if the focus is on the analysis of secular orbital transformation.
In the present work, we want to extend the results of previous studies on the orbital evolution of the bar using up-to-date galactic models.Specifically, using the language of variables originally suggested by Lynden-Bell, namely   ,  R , the precession rate Ω pr , the invariant   , and partial derivatives (Ω pr /  )   expanded to account for the 3D structure of the orbits (see Sec. 3.1), we strive to describe the dynamical structure of a "live" -body bar, which forms, grows, and changes its pattern speed, while gradually trapping orbits that experience a secular evolution of actions.This is the main goal of the present work.To do this, we developed a method for calculating action variables and frequencies directly from the orbits in the non-axisymmetric -body potential of a numerical model at any time moment, with an emphasis on the late stages of the model's evolution.
Studying the structure of a "live" bar, we want to verify whether the abnormal region, where the LB-derivative is positive, exists for a typical barred -body model, and, if it exists, to what degree the region is connected with the bar.We also want to study the orbital trapping at the level of individual orbits, which is rarely done due to how complex this problem is.Specifically, we want to check how the actions of individual orbits change when the orbits experience trapping and how the overall behaviour of the orbits changes in terms of actions and precession rates in this case.In this way, we hope to obtain a more detailed picture of the orbital capturing by a bar, which presumably occurs or has occurred in our and other galaxies.
The paper is organised as follows.In Section 2, we present our -body model and its global properties, including the evolution of the amplitude and pattern speed of the bar.We describe how we calculate actions and orbital frequencies of all particles from the -body model in Section 3. In Section 4, we show how the system under consideration evolves in terms of   ,   , and Ω pr .In Section 5, we identify the particles that are trapped by the bar during the particular time interval and check how their actions and frequencies evolve during this process.In Section 6, we discuss the picture of the bar formation/growth that stems from our analysis.We summarise our conclusions in Section 7.

SIMULATIONS
Our numerical model is taken from Smirnov & Sotnikova (2018).At the start of the simulations, the model consists of the axisymmetric exponential disc (in the radial direction), which is isothermal along the vertical axis: where  d -mass of the disc,  d -its radial scale,  d -the vertical scale height.We assume that the disc scale length  d = 3.5 kpc, the scale height  d = 0.05 d , and the disc mass  d = 5 × 10 10  ⊙ .
The dimensionless system of units is: Then the time unit  u will correspond to 13.8 Myr, the velocity unit will be about 252 km s −1 , and the angular velocity unit is 71 km s −1 kpc −1 .The disc is embedded into a spherical dark halo of the Navarro-Frank-White-like profile (Navarro et al. 1996) with a slightly steeper inner slope of the density profile (see Smirnov & Sotnikova 2018 for details).We do not consider here the model with the gaseous component or the initial classical bulge to simplify the present analysis and avoid possible complex interactions between the subsystems.The physical parameters of the halo are related to the initial parameters of the disc in the following way:  h ( < 4 d )/ d = 1.5, where  h ( < 4 d ) is the halo mass within the sphere of the radius of 4 d ;  ℎ = 6 d , where  ℎ is the halo scale length.An idea of the mass model used in this study is given by its rotation curve (Fig. 1, left plot).
The kinematics of the components is typical for -body calculations of this kind.The disc has an exponential profile of the radial velocity dispersion: (2) In practice, we set the Toomre parameter value  0 at  =  d to 1.2 and from this obtain the normalization value of  0 .The halo has an isotropic velocity dispersion profile of Cuddeford type (Cuddeford 1991).
Numerical procedures like populating the model with particles and integration of the equation of motions were carried out via mkgalaxy (McMillan & Dehnen 2007) and gyrfalcON (Dehnen 2002) program utilities, respectively, incorporated in the toolbox for -body simulations NEMO (Teuben 1995).We used 4  particles for the disc, 4.5  for the halo and followed the evolution of the system up to ≈ 8 Gyr ( ≈ 600) with the adaptive time step.The smallest value of the time step is about ≈ 0.1 Myr.
After the start of the simulation, the disc quickly evolves (≈ 1 Gyr) to a non-axisymmetric state forming the bar.The bar starts small and fast-rotating but rapidly grows in size and gradually slows down (Fig. 1, right plot).Fig. 2 shows the profile of the angular speed of the model for two specific time moments,  = 100 and  = 500.Both curves are calculated for the spherically symmetric part of the potential and differ a little.The horizontal lines define the bar pattern speed.The figure also indicates where the corotation resonance is located for the bar in our model,  CR ≈ 2 and  CR ≈ 4.5 for  = 100 and  = 500, respectively.We estimate the bar size as about  b ≈ 1.5 and  b ≈ 3.0 for these time moments using the location of the drop in the radial dependence of the bar amplitude  2 , which gives R =  CR / b ≈ 1.3 and R ≈ 1.5, respectively.This means that the model bar marginally fits the definition of a fast bar (Debattista & Sellwood 2000;Rautiainen et al. 2008;Aguerri et al. 2015;Géron et al. 2023), although it becomes slower at later time moments.

MEASURING THE ACTIONS
Action variables are a good tool for studying the dynamics of complex systems such as galaxies with bar and spirals.Recently, this approach has been widely used to study the resonance interactions of stars and non-axisymmetric features in the Milky Way, e.g., Trick et al. (2021); Kawata et al. (2021);Trick (2022); Drimmel et al. (2023).An extensive growth of this kind of studies is partially explained by the development of numerical methods for calculating the actions (Binney 2012;Sanders & Binney 2015).The usage of these methods also made it possible to study the action space of the -body models, where up to several millions particles are present (Debattista et al. 2020;Mikkola et al. 2020).Our work is in league with the mentioned works and further tests what fruitful results can be obtained if one analyses the whole action space of the system.

Different time scales in action measurement
There are several pitfalls in the measurement of actions in a nonaxisymmetric potential that we should mention.
The first problem is that in the case of non-axisymmetric potential, the momentum and the radial action  R cease to be "true" integrals of motion for orbits locked in a bar, and the actions oscillate around some specific values (Binney & Tremaine 2008;Binney 2018;Trick et al. 2021).In the non-axisymmetric rotating potential, the orbits librate or circulate around the orbits that satisfy the exact resonance condition: (Ω − Ω p ) =  (Sellwood 2010;Binney 2018Binney , 2020)).If the orbit is being trapped by the bar, these librations occur simultaneously with the secular changes in the actions.For example, the momentum of the particle can periodically change and, at the same time, show a clear tendency towards a decreasing average value (see Fig. 4).This trend is usually observed when the almost circular orbit becomes an elongated one.The existence of middle-term librations significantly complicates the analysis of orbital trapping, since they are always present and one should distinguish between them and secular evolution of actions.Overall, we distinguish three different types of oscillations in the course of the present article: (i) Short-term oscillations, i.e. fluctuations.They appear as bumps and troughs in Fig. 4 and have a time scale of about one or less radial periods.We smooth out oscillations of this type and do not consider them further in the present analysis.
(ii) Middle-term oscillations, i.e. librations, that occur on a longer time scale than fluctuations.As was already mentioned, they appear when the orbit under study is not in the exact resonance but nevertheless experiences the bar gravitational influence.(iii) Long-term changes (i.e. the decrease or increase of the actions) that happen as a consequence of the orbit secular evolution.
To study the secular evolution of the orbits in the present study, we average actions over the corresponding libration periods.Trick et al. (2021);Trick (2022) also considered librations of the orbits around the average resonance line (ARL in their notation), see figure 3 in Trick et al. (2021).
The second problem concerns the chaotic orbits and their impact on the analysis of the ensemble of orbits.Strictly speaking, the actions can be defined only for regular orbits, chaotic orbits cannot be characterized this way (Binney & Tremaine 2008).Here we do not carry out a throughout analysis of the chaoticity of the orbits in our model, since it would require a separate study.However, various studies of bar orbits in -body models (e.g., Wang et al. 2016;Machado & Manos 2016) showed that most of the orbits supporting the bar should be regular, and their fraction should increase with increasing the bar strength.Voglis et al. (2007) 1 also showed that the fraction of chaotic orbits is rather small in the case of a strong bar and at the late stages of its evolution.Moreover, in Parul et al.
(2020), we touched upon the issue of chaotic orbits in the -body model, which is also used in the present work.We proceeded from the fact that the regular orbits are quasi-periodic, thus the Fourier spectra should consist of discrete lines.Many orbits in the bar of our model are like that.A lot of the orbits have only two frequencies of comparable amplitude in the spectrum.Another important feature of the orbits that assembled into a bar in the model under consideration is their alignment along bright stripes and inside bright spots on the frequency maps (see figure 8, upper plot in Smirnov et al. 2021).This suggests that the frequencies of most orbits are expressed as integer linear combinations of  fundamental frequencies, which is true for regular orbits (Wang et al. 2016). In Parul et al. (2020), we also estimated the frequency shifts, which were small for most bar orbits.According to Wang et al. (2016); Contopoulos (2002), this also indicates the dominance of regular orbits in the model.We further support our statements regarding the dominance of the regular orbits in Appendix A, where we study Lyapunov exponents for bar orbits.
The action space of -body and theoretical models with a fixed  2021).However, the action space in these works was studied only in terms of the usual actions:  R (radial action),   (-projection of angular momentum), and   (vertical action) and, moreover, the potentials correspond to the initial axisymmetric state and do not take into account the mass redistribution induced by the bar.In this work, we investigate the evolution of actions and frequencies in the changing -body potential.In addition, we work in terms of the complete adiabatic invariant   =  R +   +   /2, which should be preserved in the three-dimensional bar of our model (see the Appendix B), the precession rates Ω pr = Ω − /2 and LB-derivative (Ω pr /  )   , which we calculate as the total derivative of the precession frequency with respect to the angular momentum Ω pr /  (a change of this derivative already occurs while preserving the invariant   ).

Instantaneous actions
In the present work, we study the bar structure in terms of   and Ω pr besides the usual actions, that is, in terms close to suggested by Lynden-Bell (1979).We also verify how orbits evolve and are trapped in the bar in the overall changing potential of the underlying system.To this aim, we calculate all three action variables for each orbit of the system: the radial action  R , the azimuthal action or the angular momentum projection onto the -axis   , the vertical action   , and the triplet of orbital frequencies  = /J: frequency of radial oscillations , frequency of azimuthal oscillations Ω, and frequency of vertical oscillations   .We use several methods for calculating these quantities.One of them is implemented in the AGAMA software package (Vasiliev 2019) and allows one to find approximate values of actions and frequencies directly from coordinates and velocities at any particular moment in time.The utilities provided by the package are based on the so-called Stäckel fudge method (Binney 2012;Sanders & Binney 2016), which evaluates action values analytically if the considered potential is close enough to the Stäckel potential (locally), which is defined as follows: where  and  are hyperbolic coordinates related to Cartesian coordinates as follows: The Stäckel fudge method from AGAMA utilities can only be applied to an axisymmetric type of potential.The potential of the -body model with the bar is not axisymmetric, of course.We obtain its axisymmetric approximation by expanding the potential into a series of multipoles and selecting only specific harmonics, namely those of 0-th order in azimuth and up to 12-th order in cos : where Φ ,0 are taken from the series expansion: where    are spherical functions of degree  and order .The Stäckel fudge method can be applied to this axisymmetric potential Φ 0 to obtain the approximate values of actions and frequencies.To make sure that the axisymmetric approximation of the potential of our model differs slightly from the non-axisymmetric one, we compared both potentials at several time moments in Fig. 3.As can be seen, there is a good agreement between the approximation of the potential and the actual potential.The maximum value of the relative error (∼ 5%) is achieved in the inner parts of the bar (the bar size is about three length units in this model).
One can assume that the actions obtained by the Stäckel fudge method for the axisymmetric approximate potential will differ slightly from their real values given that the axisymmetric and non-axisymmetric (original) potentials differ only slightly.A comparison of different methods for the same bar orbit is presented in Section 3.3.1.
As an additional bonus, the Stäckel fudge method operating with an axisymmetric potential is much faster than others methods used in this article and allows one to obtain estimates of actions and frequencies, if coordinates and velocities are known only at one point in time.We also want to point out that the axisymmetric approximation refers to the evolving potential, and not to the initial one, as, for example, in Debattista et al. (2020).
We should note the Stäckel fudge method itself is not ideal per se.In addition to errors arising from the approximation of the nonaxisymmetric potential by its axisymmetric part discussed above, there would be errors in the estimates of actions even if the real underlying potential were hypothetically axisymmetric.This is especially true for orbits with significant eccentricity, which populate the bar.However, as we will show in Section 4, this method still allows one to quickly separate bar particles from the disc and obtain a rough view of the bar structure in terms of action variables.

Averaging actions
In this subsection, we present a routine that we use to smooth out short-term fluctuations first (Sec.3.3.1)and then medium-term oscillations or libration (Sec.3.3.2).In the last part of this Section (Sec.3.3.3),we show how one can estimate the value of the LB-derivative on the libration half-period scale.

Actions and frequencies averaged over short-time scale
To begin with, we focus on calculating actions.To account for shortterm fluctuations, but still calculate all actions and frequencies over relatively short time intervals, we use two independent approaches.One of them involves calculating actions using the AGAMA package first and then averaging the obtained values between apocentres ( max ) or maximums of the vertical excursions of the orbit ( max ).Another method boils down to finding the actions directly from the orbit of a particle in an -body model using the definition of radial and vertical actions in a cylindrical potential: The actions can be found by integrating the radial or vertical velocities of the particle from apocentre to apocentre or from one maximum of  to another, respectively.Regarding the angular momentum, its instantaneous value is obtained directly from the particle velocities and coordinates.The projection of angular momentum onto the -axis   undergoes short-period oscillations (fluctuations), which we smooth out applying the same procedure as we apply to instantaneous  R and   from AGAMA.As it is done for  R , the averaging is carried out between two adjacent apocentres.
The resulting  R and   averaged between apocenters are piecewise functions of time.However, real actions should continuously change over time.To account for this problem, we apply the numerical procedure described in Ruiz-Arias (2022).Ruiz-Arias (2022) used mean-preserving interpolation for averaged data series in solar radiation modelling.This method of data smoothing allows one to transform a piecewise-defined function into a continuous one and, at the same time, preserve either the value of the integral (7) or the mean value of  R from AGAMA.We use a publicly available software implementation of the described procedure (https://github.com/jararias/mpsplines/)for our calculations.
We demonstrate an example of this procedure application to the radial action of a typical bar orbit shown in the upper panel of Fig. 5 in the second panel of Fig. 5.The blue dashed line in the second panel corresponds to  R calculated as the average of instantaneous  R from AGAMA, while the black dashed line corresponds to  R from Eq. ( 7).In turn, the blue solid and black dash-dotted lines represent the interpolation of these actions by the mean-preserving spline.We note that these two curves appear to be quite closely.Additionally, the curve of instantaneous actions from AGAMA (thin light blue curve) also differs slightly from them, which indicates that the interpolation is carried out correctly.
The interpolation results for all actions from Eq. 7 ( R and   ) and   are shown in the third panel of Fig. 5 by lines of average thickness ("Av.Apo." prefix in the legend).As can be seen, the interpolated actions follow instantaneous AGAMA actions ("Agama" prefix) quite close.
As for orbital frequencies, there are several methods for evaluating them suggested in the literature: by tracing time and angle of apocentre and -maxima passages (Sellwood & Gerhard 2020), by Fourier spectrum (Smirnov & Sotnikova 2018), and by least squares approximation (Ceverino & Klypin 2007).Each of these methods has its own advantages and disadvantages.Since we want to calculate frequencies both on small scales (of the order of one revolution) and on large time scales (in order to trace the secular evolution), we use the method based on apocentre tracing.Within this approach, the frequencies of radial and vertical oscillations are found using the following expressions: where  R and   are the time intervals between two successive passages of apocentres or  maxima, respectively.The mean value of angular velocity ⟨Ω⟩ can be obtained through the increase in the polar angle calculated between two adjacent apocentres: where Δ is the change in the angle, where the orbital revolutions are taken into account.We apply this procedure for the bar orbit.
Since we calculate the mean values, the frequencies are piecewise-defined functions again.To obtain their continuous representations, we apply the same interpolation procedure to them as we apply to actions.Note that, in this case, the mean-preserving spline conserves the total amount of the angle changes that occur as the particle moves from one apocenter to another.The interpolated frequencies  and Ω are shown in the bottom panel of Fig. 5 by medium thickness line.In addition, the thin line demonstrates the instantaneous value of an angular velocity   / 2 , the dots mark the moments of orbit apocentres passage.

Actions averaged over librations
As a second step of our analysis, we smooth out the obtained actions and frequencies over the libration period.We use the same procedure as we applied in the previous subsection to approximate the instantaneous   and  R and   from AGAMA.We take a smoothed time series of actions and frequencies, find the extrema there, and then find the mean value between two adjacent extrema.Thus, as in the previous subsection, a piecewise-defined function is obtained, which can be smoothed with the mean-preserving spline.The only difference between the procedure carried out here and in the previous subsection is that we average the values not between the maxima of the actions and frequencies themselves, but between the extrema.This is due to the fact that the libration period can be quite long.The secular evolution of actions and frequencies (averaged over a half-period of libration) is shown by thick lines in two bottom panels of Fig. 5 ("Sec."prefix in the legend).
As a finishing remark of this subsection, we want to note the following.Here, we describe three different types of oscillation in actions and frequencies that occur when the particle moves along its orbit in the realistic -body potential.In general, we smooth out short-term oscillations.The algorithm that we presented allows us to investigate the middle-term and secular evolution of actions and frequencies for each particle-"star" in the -body model.As we showed here, one can estimate actions either using AGAMA routines applying them to the axisymmetric approximation of the potential or directly integrating the velocities along the orbit in the -body potential.In the course of the present work, we mostly analyse the actions obtained directly from the -body evolving potential.Actions from AGAMA are only used for the analysis of the model at  = 0, when the model is axisymmetric, and for one later time moment for the preliminary analysis of the action space.

Estimation of Lynden-Bell derivative
In the course of the present article, we also obtain an estimation of LB-derivative (Ω pr /  )   for individual orbits.We show below that this quantity turns out to be tightly connected with how the orbit behaves under the bar influence.In the 2D case (i.e. planar orbits), the numerical estimate of (Ω pr /  )   can be obtained by replacing the partial derivative for particles that conserve their adiabatic invariant   with the total derivative: where Δ  is the difference between successive extrema of   in the corresponding time series, where the short-term fluctuations are smoothed out, and ΔΩ pr is calculated as the difference between Ω pr = Ω − /2 corresponding to the extrema of   .We choose the time between   extrema as Δ because the time-series of middleterm   is regular enough, and the extrema can be identified relatively easily (see Fig. 5).In contrast, Ω pr extrema are harder to find numerically, since Ω pr is the difference of two splines Ω and /2 and is not as regular as the time series of   .As a consequence of irregularity, the middle-term time series of Ω pr may contain false extrema.For the same reason, we have not considered the derivative found directly from the time series of medium-term values (Δ is much less than the libration period).As an additional note, we should emphasize that orbits in the -body galaxy are 3D in nature, and one should estimate the total derivative of the precession frequency with respect to the angular momentum ΔΩ pr /Δ  .
To summarise, the algorithm was described in this subsection allows us to calculate the value of LB-derivative between two successive extrema of   , and, in particular, to estimate the sign of LB-derivative.A more detailed analysis of the change in the sign of the LB-derivative for individual orbits will be carried out in Section 5.1 and for an ensemble of particles in Section 5.2.

PRECESSION RATE MAPS
The secular evolution of orbits will be studied in the next Section.In this section, we show how the general structure of the -body system with a bar looks like in the action space taken at a single moment of time.Following the Lynden-Bell (1979) approach, we present the 'portrait' of a bar on the (  ,   ) plane.Here, we show the distributions of non-averaged (instantaneous) actions and precession rates, i.e. those calculated by AGAMA without additional post-processing on our side.This is also to show that instantaneous values obtained from AGAMA can be directly used to identify the bar particles.
One of the important notions concerning the maps of precession rate in numerical models is that the maps can be constructed using two different ways.One way is to calculate the actions and precession rates of the particles constituting the particular model and produce the maps from the obtained values by averaging the values of Ω− /2 over small bins in the plane   −   .The problem, which is implicitly introduced in this approach, is that the particles constituting the self-consistent model occupy only a small part of the phase space.In practice, most of the orbits are circular (especially at the start of the simulation) and, therefore, the maps obtained in this way are not representative, i.e. they do not allow one to picture how the evolution of orbits from circular to more radially elongated can proceed in the given potential.The second approach, which we use here along with the first one, is to artificially populate a wide region of the phase space   −   by orbits and calculate the actions and precession rates using a frozen potential of -body model.The maps obtained in this way are more representative compared to the first case since they show where in the phase space the circular orbits can evolve under the influence of perturbation if they conserve their adiabatic invariant.
Maps of the precession rates in (  ,   ) coordinates for  = 0 and  = 400 are shown in Fig. 6.Each map is obtained by averaging the precession rate value over the bins of width 0.005 by   and the same width of 0.005 by   .The line   =   /2 corresponds to the line of circular orbits.As discussed above, we analyse two types of maps.For the first type, we consider instantaneous actions  and frequencies from AGAMA (left plots).Maps of the second type are presented in the right panels of Fig. 6.Here we also consider the instantaneous values of the precession rate but the maps are constructed analytically using an axisymmetric approximation of the potential at the appropriate time moments (see 3.2).In this case, we consider only in-plane motions and the adiabatic invariant is calculated simply as   =  R +   /2 (without the addition of the vertical component).In general, the maps in the right panel of Fig. 6 were built in the same way as it was done by Lynden-Bell (1979) for the isochrone model, but all integrals and derivatives are calculated numerically here.
In all subpanels of Fig. 6, we also show the location of the line (Ω pr /  )   = 0 (green colour).The line was found via the spline fitting of Ω − /2 dependence on   considered for a fixed value of   .In practice, we take a horizontal slice of the map (  ,   ) and identify the location of a zero value (Ω pr /  )   in each of the slices.The typical error resulting from this procedure is about 0.03 in corresponding simulation units.
As the left top panel of Fig. 6 shows, the most of orbits are circular at  = 0, i.e. they adjoin to the line   =   /2 and have a relatively small value of the radial action.This is expected by the construction of the model.But the unexpected thing is that in contrast to the isochrone model built by Lynden-Bell (1979), in which the abnormal region was constrained within the central part of the model, almost the entire exponential -body disc lies in the abnormal region, where (Ω pr /  )   > 0, as the top right panel of Fig. 6 shows.There, it can be clearly seen that the line (Ω pr /  )   = 0 is located at about   ≈ 3.6.Such a value corresponds to the disc boundary of the initial -body system.The location of zero of (Ω pr /  )   for circular orbits is shown by the small green square in the initial rotation curve of the model (Fig. 1).Of course, we consider a potential that differs from the one considered by Lynden-Bell (1979) and this should explain the difference.Nevertheless, the fact that for a typical -body model the normal region turns out to be almost non-existent is intriguing.In his short review of the dynamics of the Milky Way, Weinberg (2001) gives a general rule stating when two orbits should precess toward or away from each other.They precess toward (away from) each other if the rotation curve is rising (flat).Such an interpretation of the Lynden-Bell's 1979 mechanism comes from the widespread belief that the sticking of elongated orbits occurs in the central regions of galaxies, and the mechanism itself generates the so-called slow bars.We see that in realistic -body models of galaxies, the abnormal region, where orbits can be trapped in a bar-like perturbation, covers almost the entire initial disc and the zero of (Ω pr /  )   lies far from the maximum of the rotation curve (see Fig. 1, left plot).This is also consistent with the analytical result by Polyachenko & Polyachenko (2003) who found that the normal region for an infinitely thin exponential disc lies outside four radial scales2 .Moreover, recently Polyachenko & Shukhman (2020a,b) have shown that, for the initial stages of evolution, the role of the normal orbits in the bar weakening/destruction is much smaller than previously believed, since the actual values of (Ω pr /  )   , even if they are negative over the disc, are actually close to zero for several types of galaxy models typically considered.
Lower panels of Fig. 6 show what happens when the system evolves.The left panel shows the precession rates of the particles of -body model at  = 400, while the right panel demonstrates the theoretical map for the same time moment.Both panels show that the system is now divided into two distinct parts: one part, where (Ω pr /  )   < 0, and the second one, where (Ω pr /  )   > 0. It can be seen that the line (Ω pr /  )   = 0 has a slightly different shape depending on the map, but in general, there is a match between the two maps in terms of how the line runs.
We present snapshots taken for the particles occupying the regions to the left and right from the line (Ω pr /  )   = 0 and for the entire disc in Fig. 7 at  = 400.It can be immediately seen that the region, where (Ω pr /  )   > 0, consists of the particles supporting the bar.The second region is populated mostly by disc particles but also contains a portion of particles that form the socalled bar "spurs", or ansae.
As a result of this Section, we conclude that the bar is mostly contained within the region where (Ω pr /  )   > 0. In addition, we found that, for the initial (non-evolved) model, the region, where (Ω pr /  )   < 0, is almost non-existent.Although these facts can be considered as supporting the idea of orbital trapping suggested by Lynden-Bell (1979), in the next Section we will show that the situation is actually much more complex if one considers individual orbits and how they are captured by the bar, rather than the maps obtained by averaging over some bins.

ORBITAL TRAPPING
When analysing -body simulations, we must distinguish between processes occurring on different time scales.The slow secular changes in actions and frequencies occur simultaneously with the rapid oscillations of the corresponding actions and frequencies in the rotating barred potential, and also simultaneously with longerterm librations.Therefore, to trace the longer-term evolution of the actions and precession rates, we should use their averaged-in-time values obtained via the procedure described in Sec.3.3.2.But first, let us analyse the Lynden-Bell maps, free from short-term oscillations, i.e. calculated according to the scheme presented in Sec.3.3.1.In this case, we investigate actions calculated directly in the -body potential but averaged over small-amplitude oscillations.There will be inevitable librations in the dependence of such actions on time, i.e. periodic changes of a moderate amplitude on a time scale of the order of one revolution of the particle.

The Lynden-Bell maps averaged over short-term oscillations
The Lynden-Bell's (1979) idea of the bar formation due to the mutual attraction and alignment of precessing stellar orbits is based on two interconnected assumptions.The first one is that an orbit should precess relatively slowly with respect to the bar.And the second one, which stems from the first, is that the orbit, falling into the trap, must retain its adiabatic invariant   .The last condition means that the capturing process is slow compared to other characteristic times.We understand the "slowness" of all processes precisely in this sense, namely the precession of the orbits that are being trapped is slow if the precession rate is close to the speed of the pattern, and the difference of these two quantities, divided by the angular velocity Ω, is close to zero.In this case, even if the bar pattern speed slowly decreases, the adiabatic invariant will be preserved for the orbits that are being trapped.Changes in the action variables on the libration time scale have no relation to the characteristic time of orbital trapping in the bar.For capturing, what  is important is the time of gradually accumulating secular changes rather than the time of periodic changes in actions, that occur with the preservation of   throughout the entire period of capturing.
Here we investigate whether the above two assumptions hold true for the orbits in the 'abnormal' region and, in particular, for those orbits that are being trapped in the bar in our -body model.
The left panels of Fig. 8 show how the precession rate maps averaged over short-term oscillations (see Section 3.3.1)look like.Building maps in Fig. 8, we do not get rid of librations of actions and frequencies that are inevitably present for orbits that are far from circular or even regular.It is important that, unlike the maps in Fig. 6, maps in Fig. 8 use actions, frequencies, and precession rates determined directly in the non-axisymmetric -body potential.Judging by the colour bar, it can be seen that the particles in the abnormal region and immediately beyond zero of (Ω pr /  )   in the normal region precess with a speed close to that speed of the bar (dark blue line).
The right plot of Fig. 8 shows how the difference between the bar pattern speed and the orbit precession rate normalised on the value of the angular frequency is distributed over   and   .As can be seen, the precession rate of the orbits is very close to the bar pattern speed (up to ∼ 10%, in the left plot of Fig. 8, the blue lines show the contours of 10% and 15% differences from zero, see also the tips on the colour bar.) in the vicinity of the line (Ω pr /  )   = 0, i.e., the condition |Ω pr − Ω p |/Ω ≪ 1 is well satisfied, which is actually equivalent to conservation of the integral of motion   .This is true both for most of the orbits to the left of the zero LB-derivative line and for a number of orbits adjacent to the line to the right.Thus, the orbits precess relatively slowly compared to the bar before trapping.
Filling the maps with the values of   and   in Fig. 6 and Fig. 8 is very similar, which is not surprising, since the smoothed values of actions determined via  max and  max in a non-axisymmetric potential well fit the dependence of instantaneous actions on time (Fig. 5, second panel from the top) determined in an axisymmetric potential for given time moments.At the same time, the precession rates determined via the apocentres and the instantaneous precession rates, that are obtained by a routine from AGAMA, differ.This is precisely what is causing the difference in the position of the lines of zero partial derivative (Ω pr /  )   .On the maps presented in Fig. 8, a long 'corridor', that is formally a part of the normal region, appears.This 'corridor' is adjacent to the line of circular velocities and stretches from   ≈ 0.9 towards very small   values.However, there are no circular orbits in this region.The orbits there are all bar orbits that experience librations of actions.These orbits draw tracks on the plane (  ;  R +   ) with periodic changes It can be seen that particles 1, 2, and 3 move into the bar region, almost exactly following the lines of constant   .Particle 4 is in the disc and, during the time from  = 250 to  = 550, does not experience secular changes in its actions, but only periodic ones.In the next section, we will analyse the tracks of particles averaged over librations, i.e. we will trace the secular evolution of the actions of the particles that are trapped in the bar.

Trapping of particles
As can be seen from the right plot of Fig. 9, if one averages the actions for all particles over librations, then orbits trapped in the bar by  = 400 will lie completely in the abnormal region.Particle 1 is located in the abnormal region by  = 250.Particle 2 will be there at  ≈ 340, and particle 3 at  ≈ 420, after they cross the zero LBderivative line.At the same time, already trapped orbits or orbits that are being trapped demonstrate secular changes in their actions, despite the averaging over librations, and the trapping process of such orbits is rather slow (Fig. 9, right plot, particle 3).Particle 4, which is not trapped in the bar, experiences periodic changes in actions (without secular evolution), and it does not yet have a chance of getting into the bar.
Fig. 10 shows the orbits of particles (left column), the tracks of which are presented on the maps in Fig. 9. Fig. 10 also demonstrates changes with time of actions for these particles, both non-averaged over librations and averaged.For particles 1, 2, and 3, the averaged values of   (yellow lines) continue to systematically decrease even after  = 400.This occurs not only due to an increase in  R but also due to an increase in   , when the orbit, already trapped in the bar, begins to protrude from the mid-plane in the vertical direction.At the same time, for all particles the action   (red lines) is surprisingly preserved, i.e. the process of orbital evolution occurs adiabatically.
Given the above, the criterion for orbital capture by a bar can be formulated in two ways.In the Lynden-Bell language, a particle can be considered to be trapped in a bar if it crosses the Lynden-Bell zero derivative line.plots in Fig. 10) marks the time moment when the changes of   and Ω pr start to go in the same direction (to increase or to decrease simultaneously).In line with Lynden-Bell (1979), we can introduce the concepts of 'normal' and 'abnormal' orbital behaviour.To the left of the red vertical line, we observe normal behaviour, and to the right, there is a transition to abnormal behaviour 3 .
The outline of the orbit in the abnormal mode is also marked in red (left column in Fig. 10).Orbit 1 by the time  = 250 is already trapped in the bar (abnormal mode), so it is all drawn in red.Orbit 2, upon transition to the abnormal mode (practically coinciding with the moment of intersection of the line of zero LB-derivative), becomes elongated and lies completely inside a rectangle elongated along the major axis of the bar.At the same time, orbit 3, upon transition to the abnormal mode, acquires the shape of a regular bar orbit.All these changes correlate well with Ω pr behaviour in the right plots of Fig. 10.For orbit 2, when it goes into an abnormal mode, the oscillations of Ω pr are getting clung to Ω p .For orbit 3, which acquires a regular shape, Ω pr literally sits on the line of 3 The moment when the transition from abnormal behaviour to normal occurs is indicated on the tracks in Fig. 9 (right plot) by a red bold dot.
changing Ω p values after the mode change.As for orbit 4, no changes occur in its behaviour during the time interval under consideration, and Ω pr always remains below the value of Ω p .

Ensemble of trapped particles
Now, using knowledge about the sign of (Ω pr /  )   for each orbit at any time moment, it is possible to identify an ensemble of particles trapped in the bar over a certain period of time and to monitor the evolution of these particles.Based on the results of the previous section, we identify such particles as those that change the sign of (Ω pr /  )   from negative to positive.Below we show that such particles are indeed those that are being trapped.
We limit ourselves by selecting only two specific time moments  = 400 and  = 500, analysing the changes in particle actions that happen between these time moments.Fig. 11 (top row) shows the 2D density distribution of particles over (  ,   ) plane.The maps are obtained by averaging the number of particles in each bin of size of 0.005 both in   and   at time moment  = 400 (left plot) and  = 500 (right plot).The action values were also averaged over librations.That is why the 'corridor' of the normal region, adjacent to the line of circular orbits with  small values of   , disappeared on the maps.The bar particles now lie almost entirely to the left of the Lynden-Bell zero derivative line, without going beyond it.There is a narrow 'bottleneck' near the value of   ≈ 1.0 with a small scatter in   , through which the disc particles fall into the trap of the bar.
An additional, perhaps not so clear matter, is that the line of zero partial derivative changes its location with the evolution of the system potential.Therefore, we also mark the location of the line (Ω pr /  )   = 0 at left and right plots for corresponding time moments4 .It can be seen that the line does not change significantly either its position or its shape in time.
On the maps in Fig. 11 (top row), the particles that experience change from 'normal' to 'abnormal' mode from  = 400 to  = 500 are marked in green (superimposed areas with gradations of the green colour).These particles have a negative LB-derivative at the moment  = 400 (left plot) and a positive derivative by the moment  = 500 (right plot).
Fig. 11 (upper left plot) shows that the particles marked in green (with (Ω pr /  )   < 0) are located almost entirely to the right of the Lynden-Bell zero derivative line at  = 400.A very small number of orbits with (Ω pr /  )   > 0 that have not yet fallen into the trap of the bar is located on the left, in the abnormal region.We also want to note that we excluded a certain number of orbits from consideration.These are chaotic orbits in the central regions, which change the sign of the derivative several times and accidentally have a negative derivative at the time  = 400, and a positive one at the time  = 500.We formally cut off such orbits by the condition   > 0.6 (see also Appendices A and C regarding chaotic orbits and their locations).We also would like to note that for an orbit to have a negative LB-derivative value and to be to the right of the zero LB-derivative line at a given time is not completely the same thing5 .However, the differences between these two properties are small and this mainly concerns the central chaotic orbits.
We selected orbits marked in green as those that have (Ω pr /  )   < 0 at  = 400 and change the sign of the derivative before  = 500 ((Ω pr /  )   > 0).By the time  = 500, most of these orbits have moved to the abnormal region (to the left of the line of the zero LB-derivative, Fig. 11, upper right plot).A certain number of orbits continue to remain on the right, in the normal region.Among them, there are about 10% of particles that have high values of   > 1.6 and are adjacent to the line of circular orbits.These orbits have increased the rate of precession Ω pr to values close to Ω p during the studied time interval and are elliptical in shape, with the major axis elongated along the major axis of the bar.We will discuss these orbits in Section 6.2.
In Fig. 11 (bottom row) we depict the dependence of the precession rate over the -projection of the angular momentum.The bar forms a thick straight strip in this figure , while the not yet trapped particles (marked by green colour) lie on the descending branch of the precession rate curve.The latter is somewhat thick due to the differences in eccentricities of orbits.One can see that the bar perturbation rotates faster than the typical precession rate of the orbits it traps.By the time  = 500 (bottom right plot), all trapped particles increase the precession rate to approximately the same value, the angular speed of the bar Ω p = 0.17 at  = 500.Even particles located in the region of almost circular orbits and having   > 1.6 land on the 'shelf' of the bar if their orbits manage to change their behaviour from normal to abnormal.
Fig. 12 presents the evolution of the subsystem of the particles whose orbits change their behaviour from normal to abnormal over the considered time interval and end up in the 'bar' region.Initially, these particles form a ring-like structure (although, with a clearly inhomogeneous density profile, the upper left plot of Fig. 12).In the course of the evolution, the structure gradually changes its shape to ellipsoidal, becoming part of the bar.Interestingly, the characteristic farfalle-like shape is formed simultaneously in a side-on view.This indicates that the particles increase their vertical excursions and start to support the so-called B/PS bulge (a vertically thick part of the bar).

Secular evolution of actions
The distributions of the actions and the adiabatic invariant   for the trapped particles for two time moments,  = 400 and  = 500, are shown in Fig. 13.The histograms show that the particles lose their angular momentum while increasing their radial and vertical actions.This is consistent with an overall change in morphology of snapshots observed in Fig. 12: orbits become more elongated and acquire greater vertical excursions.At the same time,   is very well preserved over this period of time, and the overall shape of the corresponding distribution does not change.Fig. 14 shows the evolution of median action values over time and their deviations by one quartile.It is clearly seen that the median for the entire distribution over   practically does not change.The trapped orbits keep their   with an accuracy of about 4%-5%.The preservation of   is in good agreement with the slow, secular decrease in   and the increase in  R and   .
Summing up the results of the Section 5, we found that orbits from the normal region, where (Ω pr /  )   < 0, eventually lose their angular momentum and become trapped by a bar.To be more precise, those orbits that change their behaviour from normal to abnormal, fall into the trap of the bar.Before being trapped, such orbits systematically decrease   and increase the precession rate Ω pr .After falling into the trap, the orbits begin to change   in accordance with Ω pr .In the abnormal mode,   increases or decreases simultaneously with Ω pr , which oscillates around the value of the angular velocity of the pattern Ω p .In addition, we find that the trapped orbits initially precess slowly with respect to the bar, i.e. the Lynden-Bell condition for particle trapping is certainly satisfied.In this case, the trapping occurs with the preservation of the adiabatic invariant   .

Normal/abnormal regions versus normal/abnormal orbital modes
A crucial fact that we show for the first time is that the mature -body bar is indeed separated from the disc by the Lynden-Bell zero derivative line (Fig. 6 and Fig. 7).The latter divides the model into two regions, the abnormal one (the bar region, where the partial derivative of angular momentum over the precession rate is positive) and the normal region (disc area with (Ω pr /  )   < 0).However, Fig. 10 and Fig. 11 demonstrate that, for understanding the process of adding new particles to the bar, the important concept is not the concept of a normal/abnormal region, but the behaviour of individual orbits, i.e. the mode in which the orbit is, normal/abnormal.In normal mode, the orbit, before joining the bar, decreases its angular momentum while simultaneously increasing its precession rate until the latter approaches the angular speed of the pattern.After this, the orbit switches to abnormal mode and begins to change its angular momentum in accordance with the precession rate (Fig. 10).In turn, the precession rate begins to oscillate around Ω p (orbit 2) or even literally gets glued to Ω p (orbit 3).This prevents the orbit from jumping out of the bar trap.Now the orbit is locked in it.
It is important to note that, at the same time, an orbit in the abnormal mode is not necessarily located to the left of the zero LBderivative line (Fig. 11, upper right plot, green area).It turns out that normal/abnormal regions and normal/abnormal orbital modes are close, but not the same concept.In the normal/abnormal region, pixel-by-pixel averaging of orbital characteristics occurs.The corresponding pixel shows typical orbital behaviour rather than individual behaviour.That is why, in Fig. 11, we find orbits that reside in the normal region, have large values of   , and, at the same time, change the sign of (Ω pr /  )   from negative to positive over a time interval of  = 400-500.

Bar orbits with large values of 𝐿 𝑧
Examples of three orbits with the abnormal mode of oscillations and large values of   are shown in Fig. 15 Their shape is close to elliptical (left column of Fig. 15).The orbit shown in the top panel changes the orbital mode immediately after  = 400 (red vertical line in the middle and right subplots), while, for the orbits from the middle and lower panels, the change occurs near  = 350.It is after the change of sign of the (Ω pr /  )   that the orbits acquire an elliptical shape, and their major axis becomes elongated along the major axis of the bar.For the upper and middle orbits, this shape practically does not change, just as their   , the value of which keeps these orbits to the right of the line of the zero LB-derivative, does not change.Until these orbits reduce   and increase  R and   , they support the thin peripheral part of the bar outside the B/PS bulge (see Fig. 7, bottom middle plot).
The example of the lower orbit shows that in the future the shape of the orbit may change.In this particular example, the orbit becomes more elongated and looks like candy in a candy wrapper, like orbit 3 in Fig. 10 after falling in the bar trap.In this case, the  R increases, and the value of   drops below 1.2, after which the orbit, without changing Ω pr ≈ Ω p and maintaining   , moves to the region to the left of the line of zero LB-derivative.In the future, it can increase   , leaving the mid-plane and populating the B/PS bulge.

Expanded picture of orbital trapping by a bar
The fact that the -body bar growth occurs because the orbits are trapped mainly from the normal region contradicts Lynden-Bell's (1979) idea that the bar should form and grow from the slowly precessing orbits lying within the abnormal region.We should note that this is not the first time that the role of the normal orbits in the building of the bar has been discussed.In particular, Polyachenko & Shukhman (2020a,b) found that normal orbits can facilitate the bar growth under certain conditions.
The key difference between the case considered by Lynden-Bell (1979) and the physical situation considered here is that the bar perturbation rotates faster than the typical precession rate of the orbits it traps.As can be seen, the trapped particles from the disc area should increase their precession rate while decreasing their angular moment.But such behaviour is precisely what one should expect from particles inhabiting a normal region, i.e. such orbits behave like donkeys increasing their precession rate when slowed down and vice versa.A decrease of angular momentum leads to an increase in orbit precession rate up until the moment its precession rate becomes close enough to the bar pattern speed and the particle falls into an abnormal region, or, more precisely, switches to the abnormal mode.
Naturally, orbits can facilitate the bar growth in the described way if and only if they initially precess slower than the bar.In the opposite case, when the bar is slower than the precession rate of orbits, an orbit will run away from the bar possibly weakening it.To sum up, if an orbit possesses some channel which can be used to transfer its angular momentum then such an orbit will end up in the bar since it will increase its precession rate.Thus, fast bars can be fed by orbits from the normal region.
Polyachenko & Shukhman (2020a,b) analysed the bar formation mechanism in the framework of the standard Hamiltonian approach and individual orbit families.These authors found that not only the partial derivative of angular momentum is the deciding factor in the response of orbits to the bar but also their precession rate in the rotating frame.Polyachenko & Shukhman (2020a,b) also introduced a separate parameter (parameter  in notation of Polyachenko & Shukhman 2020a,b), which describes the orbital responsiveness to the bar-like perturbation and depends on the properties of the bar itself (mass and size) and underlying potential.In their figure 6, the authors provided a scheme that sums up their results on which orbits should be trapped depending on the combinations of all three parameters, partial derivative, precession rate, and orbital responsiveness.In particular, Polyachenko & Shukhman (2020a,b) found that the bar should trap the so-called long-axis orbits 6 from the normal region if the latter are slow (i.e. precess slower than the bar rotates, see also Fig. 11, bottom left plot, green area) and the orbital responsiveness is positive.We do not analyse the orbital responsiveness here.Nevertheless, we obtained here that orbits from the normal region can be trapped in the bar changing their orbital behaviour (mode).Thus, our results support conclusions of Polyachenko & Shukhman (2020a,b) from the perspective of self-consistent -body simulations.

Angular momentum sink
One more interesting question is what triggers the secular changes in orbits, i.e. why the orbit starts to lose its angular momentum.
Here we should note that alongside the orbit capturing, there are some global processes that also affect the bar as a whole.In particular, many numerical and theoretical studies discussed an angular momentum redistribution which occurs alongside the bar formation and growth (Hernquist & Weinberg 1992;Weinberg & Katz 2002;Athanassoula 2003;Ceverino & Klypin 2007;Weinberg & Katz 2007a,b;Sellwood 2008;Petersen et al. 2019).These studies proved that the dark halo particles play an important role in the bar life, taking the momentum from the bar particles due to the resonance interaction.At the same time, the bar slows down due to the loss of angular momentum.This is exactly what happens with the bar in the model considered here, as was already shown in the right panel of Fig. 1.Fig. 16, where secular angular momentum and radial action changes for the disc particles from  = 400 to  = 500 are presented, shows this in a more profound way.The maps are constructed by averaging the overall inflow and outflow of the corresponding variables inside the bins of size 0.04 × 0.04 during the mentioned time interval.It can be clearly seen that the bar region as a whole loses the angular momentum, and orbits become more elongated here ( R increases, see Fig. 10, orbits 2 and 3).It is also interesting that, outside the bar region, the disc is inhomogeneous in terms of angular momentum gain.In particular, there is no clear ring-like structure that can be associated with the corotation 7 (green circle), where the recipients of angular momentum should be located (Athanassoula 2003).The momentum increases in two blobs above and below the bar.But the most prominent increase in   occurs at CR near two stable Langrange points.In Fig. 17 (left column) we present examples of two CR orbits that absorb   .The angular velocity Ω of these orbits practically coincides with Ω p over the entire time interval from  = 250 to  = 550 (Fig. 17, right column).At the same time, the angular momentum shows a "stubborn" secular increase, and the orbits themselves drift outwards (compare black and blue loops of the orbit in the left panels of Fig. 17 Ceverino & Klypin (2007) found that, over one period of the bar (∼ 150 Myr), the largest changes in angular momentum occur at CR, where regions of both gain and loss angular momentum are identified.Considering individual orbits, Ceverino & Klypin (2007) showed that particles at the corotation near two stable Lagrange points move between these regions in such a way that losses and gains of angular momentum are compensated.Here, we study longer periods of time ( = 400 − 500, i. e. ∼ 1.4 Gyr) and, thus, see a secular effect.This effect manifests itself both as the average gain in angular momentum by regions near two stable Lagrange points and at the level of individual orbits.The orbits slowly but steadily increase   , drifting outward and pushing CR to the periphery.At 7 The location of the CR corresponds to the CR in Fig. 2 the same time, they maintain their angular velocity equal to the angular velocity of the bar, which gradually decreases.
Ceverino & Klypin (2007) confirmed their result regarding no net torque over the period time of the bar rotation analysing orbits in an analytical model with a bar rotating with constant speed.And this, the constancy of the bar speed, is the key difference between our results and results by Ceverino & Klypin (2007).Only a bar with decreasing angular velocity8 turns regions near stable Lagrange points into sites of angular momentum sink.

Bar slowdown as a condition for maintaining the trap
Is bar slow downing important for orbit capturing?To verify this, we carry out simple yet qualitatively fruitful experiments concerning orbit trapping.We take a probe particle, which is trapped in the bar, and analyse whether the same particle would be trapped if the bar does not slow down, that is, the bar rotates with a fixed pattern speed.Fig. 18 (left plots in each of the top four panels, consisting of three plots) demonstrates orbits of the particles 2 and 3 (see Fig. 10) starting from initial conditions taken from body simulations corresponding to different time moments.Each of the orbits is obtained by integration of equations of motion in a frozen -body potential rotating with a pattern speed that the bar has for each of the corresponding time moments (see legend in the panel where an orbit is presented).For each time moment, we use a potential of the -body model (including the dark halo contribution) taken for this specific time moment and approximated via the procedure described in Section 3. As can be seen, if the bar does not change its pattern speed the orbit does not experience any change in morphology.
In Fig. 18 (middle plots in each of the top four panels, consisting of three plots), we also show the evolution of  R ,   , and   , respectively, for each of the case described above (see legend).It can be seen that if an orbit is not trapped at any specific time moment it does not experience a change in actions and remains untrapped, keeping the precession rate below the pattern speed (plots with the coevolution of   and Ω − /2 in the top row of Fig. 18).In the middle row of Fig. 18 the orbit is shown after it has been captured into the bar.Although the speed of the pattern remains constant, the captured orbit no longer escapes the trap and retains its shape.Only when the bar slows down, the orbit undergoes a change in morphology and actions and experiences gradual trapping into the bar (Fig. 18, bottom row).
In Villa-Vargas et al. (2009), particle trapping is associated with the outward motion of CR, which, in turn, is only possible if the bar slows down (see also Li et al. 2022).Our findings, including those described in Section 6.4, are in line with these results.
Thus, our experiments clearly show that the slow downing of the bar causes changes in the orbit and may well cause the orbits to become trapped in the bar when the mode of orbital motion changes from normal to abnormal.We also want to emphasise here that in our experiments the orbits are changing in response to the bar potential changes, and not due to the interaction with the dark halo.

CONCLUSIONS
In the present work, we studied the structure of the action space for the typical bar in a simulated galaxy.We used the Lynden-Bell language (action variables, the LB-derivative), but applied it to the later stages of bar evolution.We were interested in how the action variables and frequencies change for each specific orbit of our numerical model on a large time scale.We examined correlations in the behavior of these variables.The crucial point was to determine whether Lynden-Bell's idea of normal/abnormal regions can be utilised to explain the bar growth, which occurs after the bar has formed.Those regions are characterised by a different sign of the partial derivative of the precession rates over -projection of the angular momentum for a fixed value of the adiabatic invariant   , (Ω pr /  )   .Lynden-Bell's hypothesis was that, depending on the sign of the derivative and the bar pattern speed, orbits will be either trapped along the bar and facilitate its growth or run away from the bar, thereby, contributing to its destruction.Analyzing our -body model, we moved from the concept of normal/abnormal regions to the concept of normal/abnormal orbital modes.In the normal mode,   and Ω pr of the orbit change asynchronously, and when transitioning to the abnormal mode, these quantities either increase together or decrease simultaneously.
To calculate actions and frequencies directly from the orbit in -body simulations, we have developed several procedures that allow us to obtain values of all quantities at different times on short, medium, and secular time scales.The first procedure (Section 3.2) uses the publicly available software implementation of Stäckel fudge method (Binney 2012;Sanders & Binney 2016) from AGAMA package (Vasiliev 2019) and allows one to compute instantaneous actions and frequencies from the coordinates and velocities at one point in time.The second procedure aims to extract the medium-term evolution of actions and requires knowledge of the evolution of the orbit over a sufficiently large time interval (Section 3.3.1).The mediumterm evolution of actions is found either by directly integrating the velocities along the orbit in the evolving -body potential or by smoothing instantaneous actions from AGAMA (both approaches give very similar results).In both cases, we smooth the values between successive apocentres using a mean-preserving spline (Ruiz-Arias 2022).The secular evolution of actions is found by further smoothing medium-term time series.Medium-term frequencies are found from the time interval between the successive passages of  and  maxima for  and   , respectively, and from the angle between apocentres for Ω.Medium-term and secular frequency evolution is calculated in a similar way to the actions.Most results presented in this paper are based on the analysis of medium-term and secular actions and frequencies calculated directly in the non-axisymmetric -body potential via the integration of velocities along the corresponding orbits, while actions and frequencies from AGAMA are used for the preliminary analysis of the model.
Using the described procedures, we were able to study our -body model in variables suggested by Lynden-Bell (1979).We tracked (Ω pr /  )   both for individual particles (Section 5.2) and for the ensemble of particles (Section 5.3).Our main results are as follows.
1.For the initial time moment, we found that the normal region of orbits, where the LB-derivative is negative, is non-existent while the whole disc is contained within an abnormal region, where the derivative sign is positive (Fig. 6).This contradicts an old idea that the region should be contained within a small central part of a galaxy.A similar fact was found by Polyachenko & Shukhman (2020a,b).These authors showed that the role of the normal orbits in preventing the bar formation or its destruction was overestimated by Lynden-Bell (1979) for realistic potentials.
2. At late stages, the Lynden-Bell zero derivative line can be used to roughly separate the mature -body bar ((Ω pr /  )   > 0) and disc ((Ω pr /  )   < 0) at any given time (Fig. 7).Such a division of the system allows one to distinguish most of the bar particles relatively easily using procedures from AGAMA software package.This is a good practical result, which can be used in future dynamical studies of bars in -body simulations.We also find that orbits near the Lynden-Bell zero derivative line precess relatively slowly, i. e. the difference between their precession rate and the bar pattern speed, divided by the angular velocity, is close to zero (Fig. 8).
3. However, as we show in Fig. 11, this division is not the full story, since one essentially ignores a number of orbits that are being trapped by the bar and reside to the right of the zero LBderivative line.We conclude that the concept of the Lynden-Bell zero derivative line does not provide a complete picture of the secular evolution of the bar and does not allow us to understand exactly how orbits become trapped in the bar.
4. In addition to Lynden-Bell's ideas, we build on the concept of normal/abnormal orbital modes, considering for each orbit how its behaviour changes over a large time scale in terms of the joint change in   and Ω pr .In the normal mode,   and Ω pr change asynchronously, and when transitioning to the abnormal mode, they either increase together or decrease simultaneously (Fig. 10).
5. The concept of normal/abnormal modes has shed light on how orbits become trapped in a bar in -body simulations.Trapping occurs when a particle, while in a normal orbital mode, loses   on a long-term time scale, while increasing Ω pr until Ω pr becomes almost equal to Ω p .After Ω pr reaches values close to Ω p , the sign of (Ω pr /  )   changes (becomes > 0), and   begins to librate (and even experience secular changes) synchronously with Ω pr .In abnormal mode, Ω pr no longer goes beyond the current value of Ω p , which slowly decreases over time.
6. Before switching to an abnormal mode and being trapped in a bar, the orbit retains   .This means that with a secular decrease in   , it increases  R (stretches out) and   (moves out of the mid-plane).The conservation of   and secular changes in other actions are characteristic both for individual orbits and for the entire ensemble of orbits trapped by the bar (Fig. 13 and Fig. 14).
7. The difference between normal/abnormal regions and normal/abnormal orbital modes is clearly seen in orbits with a positive LB-derivative and high   values (Fig. 15).Such orbits have an elliptical shape close to circular.Formally, they are in the normal region, but have already switched to an abnormal mode of behaviour and have a precession speed almost equal to the speed of the pattern.By decreasing   in the future, they will become more elongated and may even go out of the mid-plane, but for now, having fallen into the trap of the bar, they outline its outer boundaries and make up the peripheral thin part of the bar.8.For the transition to the abnormal mode to be possible, the particle must lose angular momentum on a long-term time scale.This means that there must be particles that gain this momentum.These could be, for example, CR particles in the vicinity of two stable Lagrange points that absorb angular momentum.These particles gradually drift outward, moving the position of the CR to the periphery.This shows up well on a large time scale (Fig. 17).9. We have also established that trapping into the bar can not occur without a secular slowdown of the bar.Experimenting with probe particles, we have found that if one freezes the bar pattern speed no orbit trapping occurs, i.e. a stationary bar does not trap orbits (Fig. 18).
Our approach and findings clarify and expand the picture of bar formation and evolution in numerical models.In the future, we plan to use our approach for analysing the evolution of   and to provide answers to a number of unresolved questions regarding how the bar acquires its thick part, the B/PS bulge.

APPENDIX A: LYAPUNOV EXPONENTS
To obtain an additional argument supporting the dominating role of the regular orbits, we supplemented our analysis of orbital chaoticity by calculating the Lyapunov exponents.We consider the subsequent results as auxiliary, since they are based on the tools of the AGAMA package and are obtained by using only the spherically symmetric part of the evolving potential.
The AGAMA software package has a special subroutine for the calculation of the Lyapunov exponent for the orbit at a given potential.For this study, we used the -body rotating potential of our model at  = 400.We then apply the standard expansion of the model potential into spherical functions (with  max = 12).This was done via the corresponding subroutine from the AGAMA package.Then, we traced the evolution of all model particles in the spherically symmetric part of the potential over 100 time units (we study the particle trapping in this time interval in Section 5.3).An estimate of the Lyapunov exponent over a finite time is one of the quantities that can be obtained by calculating the evolution of an orbit.In the case of AGAMA, its subroutine returns zero for a non-chaotic orbit and an estimate of the Lyapunov exponent for a chaotic orbit.
We plot all the chaotic particles on the (), (), and (  ,   ) planes in Fig. A1.The total number of chaotic orbits is found to be approximately 400 thousand out of more than two million bar particles (i.e. about 20% as was also found in Valluri et al. 2016), and most of them are concentrated in the central regions.  of the central orbits will also be in this mode (Fig. C1, top right plot).Orbits in the abnormal mode occupy a wide range of   from   ≈ 0 (in central regions) to   ≈ 1.3 with a small admixture of orbits with even larger   values, which are adjacent to the line of almost circular orbits (Fig. C1, top left plot).The precession speed of the orbits in the abnormal mode is practically glued to the angular speed of the pattern (Fig. C1, bottom left plot).The precession speed of the disc orbits in the normal mode is lower than the angular speed of the pattern, as is the precession speed of the central orbits in the normal mode (Fig. C1, bottom right plot).The central orbits in the normal mode coincide in their position on the   -  plane with the chaotic orbits in Fig. A1.About half of them are chaotic.In further analysis we do not consider them, separating them from bar particles by the condition   < 0.6.Fig. C2 shows the face-on and edge-on views of three subsystems.On the left, there is a bar assembled from orbits in an abnormal mode at  = 400.In the middle, there is a central subsystem, assem-bled from orbits in the normal mode with   < 0.6, and these are predominantly chaotic orbits.On the right, there is a disc identified by the condition   > 0.6 and negative LB-derivative.

Figure 1 .
Figure 1.Properties of the  -body model.Left: the rotation curve with disc and dark halo contributions for  = 0.The green square marks the location of the zero value of (Ω pr /  )   for circular orbits (see Section 4).Right: the decimal logarithm of the normalized amplitude of the bar and the pattern speed of the bar for different time moments.

Figure 2 .
Figure 2. Profiles of angular speed for time moments  = 100 and  = 500.Horizontal lines define pattern speeds for these time moments.Vertical lines indicate the location of the CR.

Figure 3 .
Figure3.The relative difference between the actual potential of the  -body model used and its approximation, Eq.(5), along the major (upper curves) and minor (lower curves) axes of the bar.The radial distance is given in dimensionless units of simulations.

Figure 4 .
Figure 4. Left: An orbit on  -and -planes over the time interval from  = 250 to  = 550 (-coordinate goes along the major axis of the bar).Right: The evolution of   on different time scales (short-term, middle-term and secular).The short marks at the bottom show the moments of passage of the apocentres of the orbit.

Figure 5 .
Figure 5.The bar particle orbit shown for the time intervals from  = 400 to  = 430 (blue) and from  = 430 to  = 500 (black) (two upper panels; ,  coordinates are given in dimensionless units of the simulations,  coordinate goes along the major axis of the bar).The evolution of its actions  R ,   , and   (middle plots) and frequencies  and Ω (bottom plot), calculated using different procedures (see the main text).

Figure 6 .
Figure 6.Maps of the precession rates for  = 0 (top) and  = 400 (bottom) on the plane   -  .Left column: maps, obtained by averaging the precession rate of the particles from the  -body model in bins of size 0.005 × 0.005.Right column: maps, obtained for the precession rate of the particles artificially populated in the phase space, with orbits of different radial actions and angular momentum.A dark green line corresponds to the condition (Ω pr /  )   = 0.A light green line shows a spline approximation of the latter.Black lines are isolines of the precession rate.

Figure 7 .
Figure 7.The face-on (upper row) and the side-on (bottom row) density plots for all particles of the disc (left column), those from the region where the LBderivative has negative values (middle column), and those from the region of the positive values (right column).All plots are shown for the same time moment  = 400.The face-on view is displayed in the square (  ) = (−4, 4) × (−4, 4) while the side-on view is displayed in the rectangle ( ) = (−4, 4) × (−2, 2).

Figure 8 .
Figure 8.The Lynden-Bell maps derived in the  -body potential for  = 400.Left: 2D distribution of the precession rate Ω pr = Ω − /2 on the plane   -  , where   =  R +   /2 +   .The dark blue line is the line of Ω pr = Ω p .Right: 2D distribution of |Ω pr − Ω p |/Ω on the same plane.In both panels, the dark green line corresponds to the condition (Ω pr /  )   =0.The light green line is the same line but smoothed via the spline approximation.

Figure 9 .
Figure 9. Tracks of several orbits on the plane (  ,  R +   ) (left: short-term averaged actions; right: actions averaged over librations).Arrows indicate the position of the beginning of the track.The background is the density distribution at  = 400.The thick green line is the line of zero value of the LB-derivative smoothed via the spline approximation at  = 400.Each track line is plotted in different colours (see the colour bar) to show how actions change with time.The thick red dot on the right plot indicates the moment of the orbital behaviour changes (see Section 5.2).

Figure 10 .
Figure 10.Left: Orbits on  -plane over the time interval from  = 250 to  = 550.Black and red colours correspond to negative and positive values of LB-derivative, respectively.Centre: Evolution of middle-term (thin lines) and secular (thick lines) orbital actions (  R +   ,   ,   ).Right: Coevolution of middle-term   (orange line) and Ω pr = Ω − /2 (blue line).Coloured arrows indicate the direction of changes in the values during the time between two adjacent extrema of   .The red vertical line in the centre and right plots corresponds to the time when the sign of LB-derivative is changed.
Here, one needs to use actions averaged over librations.It is also possible to trace at what moment   and Ω pr of a particular orbit begin to change synchronously in the same direction (increase or decrease simultaneously), without resorting to the procedure of averaging over librations.Since   is well preserved for orbits near the line of zero LB-derivative, such behaviour means that the orbit has individually acquired a positive value of the LB-derivative.The plots in the right column of Fig. 10 demonstrate the changes with time of the   (yellow lines) and Ω pr (grey lines) values non-averaged over librations.The evolution of Ω p (thin black line) is also plotted.The arrows (yellow and grey) point out the direction of   and Ω pr changes.The red vertical line (middle and right

Figure 11 .
Figure 11.Upper row: The 2D density distribution of particles in  -body model over   -  plane for  = 400 (left) and  = 500 (right).The thick black line is the location of the line of (Ω pr /  )   = 0 (smoothed via spline approximation) at the corresponding point in time.Actions are averaged over librations.The green areas show how particles changing the sign of their (Ω pr /  )   from negative to positive are distributed on the maps for times  = 400 and  = 500, respectively.The dashed horizontal line is the line of the median value of   for 'green' particles.Bottom row: Distributions of particle precession rate depending on the angular momentum -projection   .A straight blue line corresponds to the bar pattern speed with the value indicated in the annotating textbox.The green areas mean the same as in the top plots.

Figure 13 .
Figure 13.Distributions of   ,   ,  R , and   for orbits changing the sign of (Ω pr /  )   in the time interval 400 − 500 at the beginning ( = 400) and at the end ( = 500) of particle path.

Figure 14 .
Figure 14.Evolution of median action values (  ,   ,  R , and   ) for orbits changing the sign of (Ω pr /  )   in the time interval 400 − 500.Dashed areas indicate the first quartiles of the corresponding distributions.

Figure 15 .
Figure 15.The same as in Fig. 10 but for orbits with initially large   .

Figure 16 .
Figure 16.2D histograms depicting the changes in angular momentum Δ  and radial actions Δ R of the orbits from  = 400 to  = 500 in bins of size 0.04×0.04.The picture is smoothed with a Gaussian filter with  = 0.08.Black lines mark isodensity contours, and the green circle shows the location of the corotation.

Figure 17 .
Figure 17.Left: Orbits in the vicinity of stable Lagrange points on  -plane over the time interval from  = 250 to  = 550.Red colour marks the orbit at the late time interval from  = 400 to  = 550 (to show how the orbits drift outward).Right: Coevolution of   (orange line) and angular velocity of the orbit Ω (blue line).Thin grey line is Ω p .Vertical red line marks the beginning of the late stage of the orbit (after  = 400).

Figure 18 .
Figure 18.Top row: Orbits on the  -plane and evolution of actions and precession rate in the frozen  -body potential with constant Ω p .Particles are the same as in Fig. 10, i.e. orbit 2 (first three plots) and orbit 3 (second three plots).The initial moments of evolution are labelled in the plot with the orbit.The potential and the pattern angular speed are also taken for this time moment.Short marks on the abscissa axis of the plots with the evolution of the actions show the moments of the passage of the apocentres of the orbit.Coloured arrows in plots with the coevolution of   (orange line) and Ω − /2 (blue line) indicate the direction of change in the corresponding values during the time between two subsequent extrema of   .Thin grey line is the pattern speed.Middle row: the same as in the top row but for the time interval after orbit trapping in the bar.Bottom row: the same as in the top row but for a slow downing bar.Black and red colour in the plots with orbits correspond to negative and positive values of LB-derivative, respectively.

Figure A1 .
Figure A1.Left two panels: Chaotic particles (with a positive Lyapynov exponent) on the   (top left) and  (bottom left) planes.Right panels: Distribution of chaotic particles on the (  ,   ) plane (bottom).A vertical straight blue line separates central normal particles on the   −   plane for  = 400 (see Section 5.3 and Appendix C).

Figure C1 .
Figure C1.Upper row: The 2D density distribution of particles in the body model over   -  plane for particles in abnormal (left) and normal modes (right) at  = 400.Actions are averaged over librations.Bottom row: Distributions of particle precession rates depending on the angular momentum -projection   .A straight blue line corresponds to the bar pattern speed with the value indicated in the annotating textbox.A vertical straight blue line (  = 0.6) separates central orbits in a normal mode.

Figure C2 .
Figure C2.The 2D density distribution of particles in the  -body model on the (  ) (upper plots) and ( ) (bottom plots) planes for particles in abnormal (left) and normal modes.Orbits in a normal mode are divided into two parts by the conditions   < 0.6 (central orbits) or   > 0.6 (disc orbits).