Astronomy & Astrophysics manuscript no.

revamp ©ESO 2024

October 17, 2024

Orbital dynamics in galactic potentials under mass transfer

Eduárd Illés1, 2 , Dániel Jánosi2, 3 , and Tamás Kovács4, 5

1
Lund University, Division of Astrophysics, Department of Physics, Box 43, Lund SE-221 00, Sweden
2
Department of Theoretical Physics, Eötvös Loránd University, Pázmány Péter sétány 1A, H-1117 Budapest, Hungary
e-mail: [email protected]
3
HUN-REN Institute of Earth Physics and Space Science, Csatkai Endre utca 6-8, H-9400 Sopron, Hungary
e-mail: [email protected]
4
Department of Atomic Physics, Eötvös Loránd University, Pázmány Péter sétány 1A, H-1117 Budapest, Hungary
5
HUN-REN - ELTE Extragalactic Astrophysics Research Group, Pázmány Péter sétány 1A, H-1117, Budapest, Hun-
gary
arXiv:2405.16367v2 [astro-ph.GA] 15 Oct 2024

e-mail: [email protected]

Received -; accepted -

ABSTRACT

Context. Time-dependent potentials are common in galactic systems that undergo significant evolution, interactions, or
encounters with other galaxies, or when there are dynamic processes like star formation and merging events. Recent
studies show that an ensemble approach along with the so-called snapshot framework in dynamical systems theory
provide a powerful tool to analyze the time-dependent dynamics.
Aims. In this work, we aim to explore and quantify the phase space structure and dynamical complexity in time-
dependent galactic potentials consisting of multiple components.
Methods. We apply the classical method of Poincaré surface of section to analyze the phase space structure in a chaotic
Hamiltonian system subjected to parameter drift. This, however, makes sense only when the evolution of a large
ensemble of initial conditions is followed. Numerical simulations explore the phase space structure of such ensembles
while the system undergoes a continuous parameter change. The pair-wise average distance of ensemble members allows
us to define a generalized Lyapunov-exponent, that might also be time dependent, to describe the system stability.
Results. We provide a comprehensive dynamical analysis of the system under circumstances where linear mass transfer
undergoes between the disk and bulge components of the model.
Key words. galaxies: kinematics and dynamics – methods : numerical – chaos

1. Introduction about the mass distribution within the galaxy (Verheijen &
Sancisi 2001; González et al. 2010).
The theory of stellar orbits in a galactic potential is In a galactic potential, the orbits of stars and other ce-
a crucial aspect of galactic dynamics and astrophysics. It lestial objects can be broadly classified into two categories:
involves understanding how stars and other celestial ob- regular orbits and chaotic orbits (Lichtenberg & Lieber-
jects move under the influence of the combined gravita- man 1992; Contopoulos 2002). These classifications depend
tional forces within a galaxy (Binney & Tremaine 2008; on the initial conditions of the orbits and the underlying
Contopoulos 2002). The motion of stars is influenced not dynamics governed by the gravitational forces within the
only by the central mass concentration (such as a super- galaxy. The study of regular and chaotic orbits is a funda-
massive black hole in the case of many galaxies) but also mental aspect of galactic dynamics and provides valuable
by the combined gravitational influence of all the stars and insights into the overall structure, formation, and evolution
dark matter. This leads to complex and chaotic orbital mo- of galaxies. Understanding the stability of orbits is essen-
tions. Numerically solving the equations of motion for stars tial in explaining the observed structures and dynamics of
in this potential allows to model and predict the overall galaxies. (Henon & Heiles 1964; Benettin et al. 1980; Skokos
behavior of the stellar population and interstellar medium 2010; Froeschlé et al. 1997; Cincotta & Simó 2000; Laskar
from a dynamical standpoint. Galaxies can be classified 1990; Gottwald & Melbourne 2004).
based on the symmetry of their mass distribution, usually The understanding of the global structure of galaxies
as axisymmetric (Toomre 1963; Miyamoto & Nagai 1975) has rapidly gained ground with the development of many
and non-axisymmetric potentials (Binney & Spergel 1982; studies of rotation curves and mass density distributions
Contopoulos & Seimenis 1990; Cincotta et al. 1996). The calculated from observational data, which are relatively
shape of the potential significantly affects the stellar orbits. easy to interpret and decompose into different galaxy com-
The study of stellar orbits is closely related to the measure- ponents. The very first models were built from observa-
ment of a galaxy’s rotation curve, which plots the orbital tional data from our own galaxy, the Milky Way (Schmidt
velocity of stars as a function of distance from the galactic 1956; Toomre 1963; Kuzmin 1956), followed by various
center. The rotation curve can provide valuable information milestones in radio astronomy and more by the Hippar-
Article number, page 1 of 13
A&A proofs: manuscript no. revamp

cos project (Bernacca 1985), several updated mass density nian depends on time and this, in turn, affects the system’s
models (Bahcall & Soneira 1980; Caldwell & Ostriker 1981; dynamics (Lichtenberg & Lieberman 1992). Different pa-
Miyamoto & Nagai 1975), based on data from stars in and rameter shifts can lead to various interesting phenomena in
near our galaxy and other galaxies. These models are usu- Hamiltonian systems such as adiabatic and non-adiabatic
ally completed in the form of a potential, describing the parameter change, resonances, and chaos.
gravitational interactions in the interior of galaxies in a
compressed form, thus easily modelling the motions of dif- A substantial work was done in the past on galactic
ferent stars, clusters and other stellar objects. systems with potentials containing components undergoing
In the recent years the re-evaluation and refinement of perturbations. One of such works was the analysis done by
these models was possible due to more sophisticated and Pichardo et al. (2003), who looked at changes in the systems
accurate observational tools, of which one of the most im- where an axisymmetric galactic potential was perturbed by
portant and ambitious projects was the Gaia Survey (Gaia three-dimensional galactic spiral arms, which were added
Collaboration et al. 2016). One of the Gaia Survey goals as a superposition of inhomogeneous oblate spheroids, in-
was to map out the stellar population and obtain valuable troducing a periodic change into the system. This resulted
parameters, including kinematic data, in great accuracy in in the emergence of nonlinear effects in the orbits of test
order to receive an even more realistic picture of our own particles which showed significant changes in the overall
galaxy (Katz et al. 2018). behavior of the system, which can mainly be seen in the
Chaotic orbits arise when the gravitational potential is Poincaré Surface of Sections shown in the work mentioned
more complex, such as in regions with irregular mass dis- above.
tributions or significant gravitational perturbations from
nearby massive objects. Furthermore, chaotic orbits can An active research area in chaos theory is the study
arise due to orbital resonances induced by different grav- of systems subject to monotonic non-adiabatic parameter
itational influences and other kinematic characteristics of changes. In this context, a comprehensive picture of this
the orbits in a given system, which result in orbital periods type of dissipative (i.e. frictional) systems has emerged over
related by simple integer ratios between the objects. In the the last three decades, motivated by climate change. These
grand schemes, this may lead to the formation of structures systems can be considered as simple models of the Earth’s
in the phase space which give valuable information about climate when a parameter (e.g. carbon dioxide concentra-
the motion of an ensemble in certain regions of the system tion) varies monotonically. Here, it was shown that a long-
(Malhotra 1998). time study of a single trajectory is not representative, in-
Motion in a time-dependent galactic potential refers to stead the relevant phase space structure, the so-called snap-
the study of how stars and other celestial objects move shot chaotic attractor, is revealed by tracing an extended
in a galaxy when the gravitational potential is not static trajectory ensemble (Yu et al. 1990; Sommerer & Ott 1993;
but changes with time. Time-dependent potentials describe Chekroun et al. 2011; Kovács 2020). In contrast, the theory
galactic systems that undergo significant evolution, inter- of chaotic Hamiltonian (or conservative) systems subject
actions, or encounters with other galaxies, or when dy- to monotonic parameter changes has been a completely un-
namic processes like star formation and merging events take known field, but in the last years some discoveries have been
place inside certain regions. Analyzing motion in a time- made in this area, developing several methods to describe
dependent galactic potential can be more challenging than such systems (Jánosi & Tél 2019, 2021, 2022), for a review
in a static potential because the gravitational forces expe- see Jánosi & Tél (2024). One of the most important of these
rienced by objects change over time (Caranicolas (1996); discoveries is the study of so-called snapshot tori, which
Caranicolas & Papadopoulos (2003); Manos et al. (2013); are the relevant ensembles in the phase space. These are
Manos & Machado (2014); Zotos (2011, 2012, 2014); Zotos generalizations of the so-called Kolmogorov-Arnold-Moser
et al. (2020)). As a result, the orbits of stars can become (KAM) tori of stationary Hamiltonian systems, and can
more complex and less predictable on larger timescales. be defined as closed curves in the phase space of systems
There are many studies in the recent past that deal with with varying parameters, but their shape is not constant
mass exchange between different parts of a galaxy. For ex- and over time their dynamics can become chaotic through
ample Caranicolas & Papadopoulos (2003); Zotos (2012) a break-up process. Another important aspect is the gener-
investigated stellar dynamics in cases where low angular alization of Lyapunov exponents, through the calculation of
momentum stars are scattered through the nucleus into the so-called ensemble-averaged pairwise distance (EAPD).
the halo. This process transfers mass from the disk into
the nucleus, while the total mass of the galaxy remains un- In this work, we aim to explore and quantify the
changed. Their model involves exponential mass transfer as phase space structure and dynamical instability in time-
time dependent parameter change in potential. Manos et al. dependent galactic potentials simulating mass exchange be-
(2013) modeled the results of self-consistent N-body simu- tween the disk and the central bulge. We also aim to demon-
lations by introducing an analytical time-dependent poten- strate the efficiency of the recently proposed ensemble de-
tial with linear mass transfer between the disk and the bar. scription and snapshot method in time-varying Hamilto-
They showed that an increase in the mass of the bar leads nian dynamics.
to a more chaotic dynamics.
When a Hamiltonian system is subjected to a parame- The paper is organized as follows. Section 2 presents
ter shift, it means that some parameters in the Hamiltonian the model we use and the parameters obtained from radial
function are changed with time. These parameter shifts can velocity fits. Section 3 is devoted to the phase space analysis
occur due to various reasons, such as external influences of the steady-state and time-dependent stellar dynamics. In
or changing physical conditions in the system. The specific Section 4 we summarize our results and draw conclusions
form of the parameter shift will determine how the Hamilto- on the implications.
Article number, page 2 of 13
Eduárd Illés et al.: Orbital dynamics in galactic potentials under mass transfer

2. The model and the NFW profile in terms of Ms can be written in the
following form using the cylindrical coordinate system:
2.1. Potential profiles for the components √
D ln(1 + R2 + z 2 /d)
In this work the analysis was done on Milky Way-like ΦDM (R, z) = − √ , (4)
galactic potential models. These models can be described as ln 2 − 1/2 R2 + z 2
a combination of three separate potentials, each represent- where D = GMs is the enclosed mass multiplied by the
ing a prominent component of a galactic system: a galactic gravitational constant G and d = rs .
disk, a central bulge and a dark matter halo component.
The overall model of a Milky Way-like disk galaxy is
simply the sum of the disk, bulge, and dark matter halo
2.1.1. Galactic disk potentials written in Eqs. (1, 2, 4):
For several decades the most well-known and commonly Φ(R, z) = Φd + Φb + ΦDM . (5)
used galactic potential model was the one described by
Miyamoto & Nagai (1975), which is a mathematically sim-
ple, yet good approximation when it comes to the numerical 2.2. Equations of motion
analysis of such systems. For the disk components for our
The motion of a test particle (star or interstellar object)
constructed models an extended Miyamoto-Nagai poten-
can be described by the following Lagrangian, with Lz = pϕ
tial was chosen, which is derived from the Poisson equation
being a constant representing the angular momentum com-
(Vogt & Letelier 2005), and was used in the analysis of 3-
ponent around the z-axis since the models under consid-
dimensional axisymmetric potential models for the Milky
eration exhibit axial symmetry (more precisely, cylindrical
Way by Barros et al. (2016).
symmetry):
The equation for the extended Miyamoto-Nagai disk po-
tential is as follows:

1
L= 2 2
Ṙ + ż − Φ(R, z), (6)
A

a(a + z̄) 2
Φd (R, z) = − p 1+ 2 −
R2 + (a + z̄)2 R + (a + z̄)2 where (R, z) represent the cylindrical coordinates, and we

1 a2 R2 − 2(a + z̄)2
incorporate the summed potential given by Eq. (5) into
− , (1) the expression for Φ(R, z). The corresponding generalized
3 R2 + (a + z̄)2 2

momenta are

where R and z are the cylindrical coordinates, A is the pR = Ṙ, pz = ż, (7)
mass of the disk component
√ multiplied by the gravitational and thus, the Hamiltonian can be written as
constant G, while z̄ = z 2 + b2 . The parameter a is the
radial scaling length, while b is the scaling height, as it is 1 2
H= (p + p2z ) + Φeff (R, z), (8)
commonly denoted in the Miyamoto-Nagai model. 2 R
where the effective potential Φeff is given by
2.1.2. Bulge
L2z
For the central bulge, which is commonly approximated Φeff (R, z) = + Φ(R, z), (9)
2R2
as a spheroidal component, we adopted the Hernquist po-
tential as a simple alternative (Hernquist 1990). where Lz is the z-component of the angular momentum.
The potential takes the following form: Thus, out of the original 6 generalized coordinates, only
4 remain, allowing us to define 4 first order differential equa-
C tions describing the motion:
Φb (R, z) = − √ , (2)
R2 + z2 + c
where C again is a mass multiplied by the gravitational Ṙ = pR , ż = pz ,
constant G, while c is the scaling radius of the spheroidal L2z ∂Φ(R, z) ∂Φ(R, z) (10)
bulge. p˙R = − , p˙z = − .
R3 ∂R ∂z

2.1.3. Dark matter halo 2.3. Fit of radial velocity profile

For the dark matter halo component a Navarro-Frenk- The rotation curve equation can be obtained by the fol-
White (NFW) profile (Navarro et al. 1996) was chosen, lowing general relation:
written in terms of the scaling radius rs and the enclosed 12
mass Ms (Sanderson et al. 2017), instead of the commonly

∂Φ(R, z)
vc (R) = R . (11)
used form determined by the scale density ρs or virial ra- ∂R z=0
dius rvir and mass Mvir , in order to keep the potential
components mathematically consistent. The equation for This equation is needed for determining the parameters
the scaling mass is as follows: a, b, c, d, A, C and D. Two main sources of data were used
as a combined data set: one is the grand rotation curve
Ms = 4πρs rs3 (ln 2 − 1/2), (3) of the Milky Way compiled by Sofue (2012, 2013), while
Article number, page 3 of 13
A&A proofs: manuscript no. revamp

the other one is the collection of rotational velocity data

obtained using Classical Cepheids by Mróz et al. (2019),
ensuring that the inner parts of the rotation curve around
the range of R = 0 − 20 kpc is well-captured. After that
the parameter space can be explored to find the best-fitting
parameters for constructing a model representing a Milky
Way-like disk galaxy.
The procedure to find the best parameters was done in
two steps:
– First, the rotation curve Eq. (11) was fitted using Trun-
cated Newtonian Method, or TNC, (Dembo & Stei-
haug 1983) for the parameters a, c, d, A, C, D to receive
approximate values in preparation for the next step.
The initial parameter guesses of scaling distances and
masses, along with the lower and upper limits for the
parameter space, for the TNC fit were obtained from the
DR3+ data analysis results of Labini et al. (2023) and
the works of McMillan (2011) and McMillan (2016). The
fits were done using the lmfit Python library (Newville
et al. 2015).
– After that, the fitted parameters from the TNC were
used as initial conditions for an Affine Invariant Markov
Chain Monte Carlo, or MCMC model (Goodman &
Weare 2010). The change in disk scale height has been
shown to give minimal contribution to the main features
of such models, as was shown in McMillan (2011), so it
was chosen as a constant b = 0.25 kpc, being the low- Fig. 1: The MCMC results on the gathered rotational ve-
est used scale height in the aforementioned work. The locity data illustrated in a corner plot. As expected, the
software of choice used for performing the MCMC sam- highest correlations can be seen between the scaling dis-
plings was emcee (Foreman-Mackey et al. 2013). The tances and mass components corresponding to the same
number of walkers were chosen to be 1000, with a sam- potential components. The red dots in the two-dimensional
pling time of tsmp = 106 and a burning time of 5 times projection, along with the red lines in the one-dimensional
the time of highest auto-correlation time tac ≈ 200. The projection, show the Maximum Likelihood estimates for the
same lower and upper constraints were implemented for parameters, along with the black lines showing the median
the parameter space in this analysis as in the previous and 1σ confidence interval results.
step. The results of this analysis can be seen in Fig. 1,
while the obtained parameters are summarised in Ta-
ble 1.
The rotation curve resulting from the fitting is shown
in Fig. 2. The fitting was done on the joint datasets of So-
fue (2012, 2013) (blue) and Mróz et al. (2019) (orange),
and shows quite good agreement with the observations. Al-
though some of the resulting parameters deviate from the
most recent estimates for the Milky Way (Ou et al. 2024;
Labini et al. 2023; Posti & Helmi 2019), they remain within
acceptable bounds for constructing a Milky Way-like galac-
tic potential model. These parameters are adequate for the
upcoming analysis, where the primary goal is to have a
model which may reproduce the rotation curve in the inner
regions of the galactic model as seen on Fig. 2.

3. Results
3.1. Time-independent case Fig. 2: Rotation curve fitted to the observational data in-
The energy of a trajectory initiated from a particular dicated in the legend, using parameters obtained from the
set of initial conditions can be determined by evaluating MCMC analysis of Fig. 1 and are shown in Table 1.
the Hamiltonian-function (8) using the given initial condi-
tions. Assuming a drift-free case, the energy of the trajec-
tory is conserved at all times. Therefore, we can express it Merely tracing trajectories is far from sufficient for an
as follows: extensive examination of a system, where a deeper under-
standing of a specific galaxy model’s morphology and the
structures within it can be obtained. Investigating trajec-
H(R, z, pR , pz ) = H(R0 , z0 , pR 0 , pz 0 ) = H. (12) tories of a system with multiple initial conditions using
Article number, page 4 of 13
Eduárd Illés et al.: Orbital dynamics in galactic potentials under mass transfer

Table 1: Table of the parameters obtained from the TNC and MCMC methods. The final parameters used in the fitting
of the rotation curve are indicated in green. The scaling height b was kept at a constant 0.25 kpc.

MCMC
Initial TNC
Median MLE
a [kpc] 4.5 9.150 9.190+0.042
−0.041 9.179
c [kpc] 0.25 0.253 0.265+0.003
−0.003 0.265
d [kpc] 183 124.9 +13.7
129.0−12.4 126.3
A [kpc km s−1 ] 1.29 · 105 2.29 · 105 (2.23+0.021
−0.021 ) · 10
5
2.22 · 105
C [kpc km s−1 ] 8.6 · 104 6.48 · 104 (6.60+0.053
−0.053 ) · 10
4
6.60 · 104
D [kpc km s−1 ] 6.9 · 106 1.23 · 107 +0.27
(1.34−0.23 ) · 10 7
1.29 · 107

physical and phase space plots, depending on the number with a resolution of dt = 10−4 Gyr, which is sufficient to
of initial conditions, may not be efficient or even feasible. determine the intersections within a certain time interval,
Therefore, it is required to continue the investigations using to obtain detailed Poincaré SoSs .
the Poincaré surface of section method (Ott 1993; Teschl Figure 3 show the phase space plots in the (R, pR ) plane,
2012). with a colormap representing the value of pz at each point,
Each intersection is registered under the conditions z = obtained from Eq. (12) with z = 0. With a chosen refer-
0 and pz > 0. To clarify, from this point on all the units ence of Lz = 100 we illustrate the phase portrait for four
for momenta, angular momenta, Hamiltonian energy and different values of constant H. The red curves are high-
parameters A, B, C will be indicated per unit mass, making lighted KAM tori of initial conditions R0 = 1.34, pR 0 = 0
them specific quantities. (Fig. 3b) and R0 = 2.96, pR 0 = 267.37 (Fig. 3c). The anal-
In order to obtain the initial conditions for the simu- ysis of the time-dependent dynamics to follow in the next
lations we firstly need to find the circular orbit at Rcirc = section will be done on these curves.
8.249 kpc and zcirc = 0.0208 kpc, which is the estimated lo-
cation of our Sun in galactocentric coordinate system from
GRAVITY Collaboration et al. (2021), and calculate the ra-
dial velocity VR ≈ 226.77 km s−1 from the derived rotation
curve, which then can be used to acquire the z-component
of the angular momentum (per unit mass) Lz = VR R. After
that the reference Hamiltonian energy (per unit mass) for
a circular orbit H0 = −5.26 · 105 km2 s−2 can be calculated
using Eq. (8), with pR 0 , pz 0 = 0.
By lowering the Lz parameter it is possible to create an
ensemble of orbits which have their Hamiltonian energies
distributed to the pR , pz velocity components, resulting in
orbits which are misaligned with the galactic plane. Stellar
streams can have such orbits, and in the last couple of years
they have been proven to be useful tools for measuring the
Milky Way’s gravitational potential (Sanderson et al. 2017;
Walder et al. 2024). Figure A.1 of the Appendix illustrates
the effect of different Lz s while keeping H0 fixed.
It is possible to find initial parameter sets by setting
pR0 = 0 and calculating the maximum permitted value for
the velocity in the z-direction pzmax using Eq. (8), with
known H and a chosen Lz , and scanning the phase space Fig. 3: Phase portraits of Eq. (10) for different values of H
with an arbitrary resolution in the range of pz0 = [0; pzmax ]. with constant Lz = 100. The trajectories are plotted on the
By simultaneously going through the range of radial dis- (R, pR ) plane, following the Poincaré cut at z = 0, and the
tances R0 = [0; Rcirc ] it is possible to find all the possible value of pz obtained from (12), indicated at each point by
orbits in the phase space with a given H0 and Lz parameter the colormap. The red KAM tori belong to initial conditions
values. R0 = 1.34, pR 0 = 0 (panel b) and R0 = 2.96, pR 0 = 267.37
(panel c).
A total of 400 initial conditions were launched in
each chosen Lz , H values. The numerical integrations were
performed using the NumbaLSODA integrator (Wogan One can observe that the system possesses a typical di-
2021), which uses the LSODA algorithm with automati- vided phase space, with the coexistence of four types of dy-
cally switching between stiff and nonstiff methods during namics. Quasi-periodic motion occurs on the closed curves,
the integration depending on which of the two can more i.e. the KAM tori, where trajectories densely fill up these
likely solve the equation in most efficient manner, which shapes. These tori form concentric islands, where the ra-
can especially be beneficial in cases where the nature of dius of the shapes depends on the potential energy of the
the problem is unknown or may change drastically during trajectories, which is conserved within a curve. In the cen-
the chosen time intervals (Petzold 1983). The simulations ter of such islands one finds so-called elliptic fixed points,
in the time-independent case were run for t = 5000 Gyr representing periodic trajectories.
Article number, page 5 of 13
A&A proofs: manuscript no. revamp

The third type of dynamics is chaos, which manifests by a strong stretching in their shape. This is called the
in this illustration as dense bands of irregularly mixed tra- break-up of snapshot tori, and it ends with them becoming
jectories, called the chaotic sea, observable along the edge entrained into what is called the snapshot chaotic sea, the
of the phase space, between the large central and several time-dependent image of the stationary chaotic sea, whose
smaller islands. Motion within this region is ergodic, mean- shape is also changing in time.
ing that trajectories uniformly cover the chaotic sea, and It should also be noted that the snapshot picture taken
have equal probability to reach any parts of it. Here the on consecutive Poincaré sections implies intersections at
dynamics is unpredictable, and trajectories, even if they various time instants, which are in different stages of the
were very close initially, tend to diverge at an exponential parameter drift. Consequently, the energy does not remain
rate. constant. Therefore, what we depict is a projection of the
Although the shown phase portraits do not capture (R, pR , pz ) phase space onto the (R, pR ) plane, which has
this explicitly, it is known from chaos theory (see e.g. Ott evolved from a stationary Poincaré section.
(1993)) that a peculiar phenomenon occurs on the border The most straightforward approach to elucidating our
between the island of KAM tori and the chaotic sea. In simulation results is through visualization. In Figs. 4-7, the
a narrow band outside the outermost KAM torus the tra- time-dependent simulations showcase the time evolution of
jectories do not form any comprehensive pattern, they are snapshot tori in the galactic model previously derived and
scattered just like in the chaotic sea, however their dynam- discussed. From such images we can easily acquire geomet-
ics is not yet chaotic, they do not diverge exponentially. rical information about the general behavior the evolution
Orbits are able to detach from this region into the chaotic of the selected ensemble throughout the integration time.
sea (and vice versa), however the transport of trajectories In Fig. 4 we follow the snapshot torus whose initial im-
in this direction was found to follow a power-law (Channon age is the red KAM torus shown in Fig. 3c. Its evolution is
& Lebowitz 1980; Karney 1983; Meiss et al. 1983; Crista- similar to those shown in Jánosi & Tél (2019, 2021, 2022);
doro & Ketzmerick 2008; Altmann et al. 2013), which is Jánosi & Tél (2024): for some time it does not appear to
much slower than the ergodic wandering in the chaotic sea. show any peculiar change in its dynamics, which can per-
Because of this, chaotic trajectories can be "trapped" in- sist for longer times, like we see here at n = 960, although
side these regions and "stick" to the outermost torus for one can observe slight deformations in its shape, particu-
long times. This is the fourth type of motion present in the larly on the left. However, this suddenly changes when the
phase space, and is often referred to as sticky dynamics. torus strongly stretches in multiple directions (n = 1002),
followed by disintegration and in the end fully chaotic dy-
namics (n = 2210).
3.2. Parameter drift and EAPD-curves For a quantitative analysis, we turn to the ensemble-
The simulations conducted in Section 3.1 were calcu- averaged pairwise distance (EAPD) method (Jánosi & Tél
lated by keeping the parameters received from the model 2019, 2021, 2022; Jánosi & Tél 2024), which allows measur-
fittings fixed, assuming that the physical properties associ- ing the dynamical instability of time-dependent chaos. This
ated with these parameters, as well as those derived from method can be seen as a generalized form of the concept
them, do not change over time. The changes occurring in of Lyapunov-exponents, where the fundamental behavior of
these systems are generally considered relatively slow, mak- the system (or significant parts of it) can be deduced from
ing these assumptions reasonable approximations for ex- the ensemble average of distances between pairs of points
trapolating behaviors of certain galaxies, which can be con- in the phase space. In most analysis methods, these dis-
sidered isolated and stable for a longer period of time when tances are evaluated on projections of the phase space. The
it comes to the stage of stellar evolution. definition of EAPD is as follows:
In this section, we assign a time dependence to the A ρn = ⟨ln rn ⟩, (14)
and C parameters by introducing linear drifts as
where the bracket denotes the average over the ensemble,
A(t) = A(0) − µt, C(t) = C(0) + µt, (13)
ln rn represents the natural logarithm at the n-th Poincaré
that is, the decrease in A(t) and the increase in C(t) has section of the phse space distance rn between one point
the same µ rate, i.e. we simulate a mass transfer from the in the snapshot
√ torus, and its assigned pair of initial dis-
disk to the bulge, so that the overall mass of the galaxy is tance r0 = 2 · 10−12 , which is not on the snapshot torus
conserved. The initial values A(0), C(0) are read from Table itself but very close to it, which turns out to be a good
1, and the rate is chosen to be µ = 10 3
A(0) · T1f , with the approximation for our purposes. If the temporal change of
final intersection time Tf = 200 Gyr. the distances follows an exponential trend, we can fit a lin-
Following the methods used in Jánosi & Tél (2019, ear function defined on a natural semi-logarithmic scale,
2021); Jánosi et al. (2021); Jánosi & Tél (2022); Jánosi where the slope of the function is the so-called instanta-
& Tél (2024), we have to accept that simulating individual neous Lyapunov exponent denoted by λ. In the above case
trajectories is no longer reliable. Instead, one has to fol- λ is constant, while in general, the instantaneous Lyapunov
low an ensemble of trajectories in order to obtain sufficient exponent is given as the derivative of the EAPD function,
statistics. Furthermore, in Hamiltonian systems the shape
λn = ρ̇n , (15)
of this ensemble is not arbitrary, they have to be initiated
on e.g. KAM tori of the stationary system. In our case, we where the dot denotes the time-derivative evaluated in dis-
will focus on the two red tori highlighted in Fig. 3. The crete time. This approach allows making assumptions about
analysis includes studying the temporal evolution of the the system in a manner similar to the traditional counter-
tori, whose shapes become time-dependent, and are called part. If the rn distances diverge linearly during the time
snapshot tori. It was found that at some point, the dy- evolution, it indicates regular or periodic behavior, whereas
namics of these snapshot tori becomes chaotic, indicated in chaotic cases, it results in exponential distance growth.
Article number, page 6 of 13
Eduárd Illés et al.: Orbital dynamics in galactic potentials under mass transfer

the standard deviations and means of phase space coordi-

nates are comparable to each other. This also means, that
the distance growths in the following EAPD curves are not
direct indications of the distances on the Poincaré sections,
rather it is a normalized approach which, however, does
capture the overall form of the evolution of distances in the
actual phase space.
The EAPD function evaluated for the snapshot torus of
Fig. 4 is shown in Fig. 5. This means that the average de-
noted by the bracket in (14) was evaluated on the ensemble
representing this torus. One can clearly see the sudden in-
crease in distance starting from n = 1002. The slope of the
linear fit is λ = 0.373±0.006. After this, the curve goes into
saturation, since the available phase space is of finite vol-
ume, thus once all trajectories of the snapshot torus become
chaotic, their average distance cannot grow any further.
This can be thought of as the typical behavior of snapshot
tori in Hamiltonian systems subjected to parameter drift,
at least in cases where the dynamics tends towards stronger
chaos as the consequence of the drift, which seems to be the
case here. For example, notice that in Fig. 3c the red torus
is in the middle of the island, surrounded by other KAM
tori. However, by the end of the parameter drift scenario at
n = 2100 it has dissolved into a chaotic sea that spans the
whole outer region of the phase space, meaning that all of
the tori that were outside the red one had to have broken
up at some point.
Note that during the period when the dynamics is still
regular (i.e. before n = 1002) the distance between point
pairs does grow, although slowly. An explanation for this
could be that since the assigned pairs for each point are
not exactly on the snapshot torus, during the slowly de-
forming phase they drift slightly further apart from the
torus. Also, there is a short jump at the very beginning of
the curve, which likely happens because the initial assigned
pairs technically correspond to different tori (since they
continuously cover the island), possessing different winding
numbers from the one investigated.

Fig. 4: Time evolution of the red torus on Fig. 3c. For

around 1000 iterations it only deforms slightly, but at
n = 1002 it starts to break up, resulting in the dynam-
ics of more and more trajectories of the snapshot becoming
chaotic.

During the calculations of distances in the phase space

at the moments of intersections it is important to take the
differences in units into account, since it is evident that the
velocity components pR , pz would overbear the distances in Fig. 5: EAPD curve of Eq. (14), evaluated for the snapshot
R. In order to remedy this problem all the phase space di- torus of Fig. 4. The Lyapunov exponent obtained from the
mensions are standardized using z-score transform, so that fitting is λnr = 0.373 ± 0.006.

Article number, page 7 of 13

A&A proofs: manuscript no. revamp

The chosen form of the parameter drift in (13) represent-

ing a mass transfer makes for a non-trivial dynamics even
compared to previous studies Jánosi & Tél (2019, 2021,
2022); Jánosi & Tél (2024), where a simple linear drift in
only one parameter was applied, typically obtaining results
similar to Figs. 4 and 5. In our case however, we found that
the inclusion of just one more time-dependent parameter,
even with the same drift rate, leads to more diverse dynam-
ics on a broad scale of parameters H and Lz . We provide
an example of such peculiar dynamics in Fig. 6 by follow-
ing the snapshot torus initiated from the red KAM torus in
Fig. 3b. Another example is shown in Figs. A.2 and A.3 of
the Appendix.
The first significant part of the evolution of the snapshot
torus in Fig.6 begins at n = 50 with its corners extending,
until about pR = ±350, while some trajectories remain on
the previous border of the torus, essentially closing in parts
of the phase space. A similar thing happens around pR = 0
as well. The created three "closed curves" then seem to be
evolving on their own (panels c-f), inside the outer points
of the snapshot torus, until about n = 100, when they ex-
perience stretching, as if they were separate snapshot tori
themselves, and start to dissolve. By around n = 482 these
points have fully scattered into a snapshot chaotic sea. An-
other peculiar feature to mention is the appearance of small
crinkles at n = 255, developing on the left side of the torus,
the shape of which later dissolve, creating a narrow band
(panel h).
The behavior of the snapshot torus between n = 50
and n = 140 can be explained if one takes a look at the
initial phase space of Fig. 3b. Outside of the red torus,
about pR = ±500 there are two smaller islands. Within
these islands there are KAM tori just like in the largest
one, and if we were to choose these tori as initial ensembles,
they would become snapshot tori as well. Along those two
islands we may also see 4 even smaller elongated islands
scattered in the chaotic sea at around pR = ±200, ±350.
In the time-dependent scenario all those islands move and
deform in the phase space, approaching the outside of the
snapshot torus in Fig. 6. However, what ends up happening
is that points of the large snapshot torus wrap around all
the smaller islands. This state holds until n = 100 when
the snapshot tori within the small islands start to break
up, taking with them the points of the large snapshot torus
around them. In the same way, the behavior around n = 255
is the result of its points wrapping around an island chain
consisting of even more smaller snapshot tori.
The EAPD curve for this snapshot torus is shown in
Fig 7, and in this case there is no clear section of exponential
divergence of the trajectories, (like the one we saw in Fig. 5).
The jump in the beginning can again be explained by the
initial difference in winding numbers, but for the rest of Fig. 6: Time evolution of the red torus on Fig. 3b. At n = 50
the curve the distance only grows slowly, never reaching its points start to wrap around the snapshot tori initially
above ρ = −10. This is despite the fact that by the end belonging to the islands around the initial torus, until these
of the simulation (n = 482) the trajectories are spread out break up starting at n = 100. Points of the large snapshot
in an extended region of the phase space. This means that torus are dispersed by instant n = 482.
the dynamics of this snapshot torus is never (fully) chaotic,
there has to be some other mechanisms at play.
Let us remember that in time-independent Hamiltonian even when one chooses a stationary chaotic sea as initial
phase spaces on the border of islands there exist the so- ensemble, when it is evolved in time as a snapshot chaotic
called sticky regions, which exhibit non-chaotic dynamics, sea, some regions in it could exhibit non-chaotic dynamics.
outside of which one finds the chaotic sea. For systems sub- These regions are termed time-dependent non-chaotic re-
jected to parameter drift, it was found (Jánosi & Tél 2024) gions, and trajectories in them can either become chaotic
that the picture is even more intricate. It turns out, that quickly, or stay non-chaotic for longer times. The subsets of
Article number, page 8 of 13
Eduárd Illés et al.: Orbital dynamics in galactic potentials under mass transfer

tegration. In Fig. 8, along the energy curve of the snapshot

torus shown in Fig. 4, a slight break is visible at the time of
torus break-up. This characteristic becomes even more pro-
nounced when considering the energy dispersion (σ) evalu-
ated over the ensemble. The sudden change in σ around the
disappearance of the torus can be regarded as an alterna-
tive indicator of the dramatic structural evolution of KAM
islands in phase space during parameter changes.

Fig. 7: EAPD curve of Eq. (14) evaluated for the snapshot

torus of Fig. 6.

time-dependent non-chaotic regions where the slow dynam-

ics is persistent are called time-dependent non-hyperbolic
regions, and such latter sets are thought to be the time-
dependent analogs of sticky regions. Since on the EAPD
curve of Fig. 7 no chaotic section was found, and because
in Fig. 6 points of the snapshot torus was observed to wrap
around islands of other snapshot tori, we can speculate
that the time-dependent non-hyperbolic regions described
above, existing around the smaller islands, might be respon-
sible for slowing down the dynamics of the large snapshot
torus.

3.3. Time-dependent energy

Changing the parameter in a Hamiltonian system re-
sults in a violation of energy conservation. The dynamics
of the system can be described by the well-known adia-
batic invariant for sufficiently slow parameter shifts, which
is represented by the volume enclosed by a periodic tra-
jectory in phase space. In this framework, instead of total Fig. 8: Median energy of the ensemble of the snapshot torus
energy conservation, an adiabatic invariant can be consid- shown in Fig. 4 compared to the EAPD curve of Fig. 5. The
ered a constant of motion, provided that the time scale of 1σ dispersion is shown as a light red band around the energy
the periodic orbit is much smaller than that of the parame- curve, while the vertical dashed line marks n = 1002, the
ter change. The conservation of the action variable depends start of the torus break-up.
on the slowness of the parameter variation (Lichtenberg &
Lieberman 1992).
In this study, the speed of parameter variation is cho- For the case of the snapshot torus belonging to Fig. 6,
sen such that the system evolves beyond the adiabatic limit. the energy also monotonically decreases, with the disper-
That is, we do not expect even the action variable to re- sion growing much more widely than in Fig.8, starting at
main constant during the motion. This fact is also evident around n = 50.
from the performed numerical simulations, where not only The Hamiltonian function (12) determines conserva-
the geometric structure but also the size of KAM islands tive dynamics while the parameters remain unchanged over
changes over time in phase portraits. During the mass ex- time. However, if we allow the mass between the disk and
change between the two parts of the galaxy, Eq. (13), the the bulge to evolve, the internal dynamics experience time-
total energy of the system is monitored. It turns out that dependent perturbation, which affects the initial energy of
- see Figures 8 and 9 - the system becomes increasingly the system, as shown here. In other words, manually setting
bound as the median energy decreases. The energy curves the parameters appears as an external disturbance in the
are plotted alongside the respective EAPD curves of Figs. 5 dynamics. Consequently, we model stellar dynamics in a
and 7, in order to compare and highlight relevant parts of time-dependent potential without employing any concrete
the processes. The relative change in total energy averaged physical mechanisms. However, similar results are expected
over the ensemble is less than 1.2% by the end of each in- if we had taken into account an additional term explicitly
Article number, page 9 of 13
A&A proofs: manuscript no. revamp

quasi-periodic, chaotic and sticky dynamics. When the pa-

rameter drift is introduced, as an appropriate ensemble of
initial conditions we chose quasi-periodic KAM tori of the
stationary case. In the time-dependent dynamics they be-
come snapshot tori, which change their shape in the phase
space and in one of the investigated cases their dynamics be-
comes chaotic through a break-up process. The torus break-
up results in an exponential separation of nearby points
on the torus. This is characterized by the so-called EAPD
curve, that describes the evolution of the average distance
of point pairs. The slope of this curve can be thought
of as a new type of finite-time Lyapunov exponent that
characterizes the fate of the torus. In the case mentioned
above, this Lyapunov exponent is approximately constant
and strictly positive through the break-up process. How-
ever, we found that in other cases the EAPD curve does
not contain any exponential parts, the dynamics of the
torus never turns chaotic. This slow dynamics results in
non-trivial Lyapunov-exponents throughout the whole pro-
cess. We note that the two presented cases are far from fully
covering the rich dynamics of the problem, further studies
are needed for its comprehensive discovery.
As it was noted in Sec. 3.3, the snapshot pictures com-
prised of different intersections with various time instants of
the parameter drift can result in changing trajectory ener-
gies. Future studies will be necessary to further explore the
emerging changes induced by different parameter changes in
order to better understand the dynamics of time-dependent
systems. The change in trajectory energies may play a sig-
nificant role in inducing local, or even global, changes in the
morphology of these galaxy models, like the emergence of
more dense star forming regions with distinct structures by
Fig. 9: Median energy of the ensemble of the snapshot torus purely dynamical causes (Daffern-Powell & Parker (2020);
shown in Fig. 6 compared to the EAPD curve of Fig. 7. The Rieder et al. (2021); Arnold et al. (2022)). Classical bulges
1σ dispersion is shown as a light red band around the energy are believed to form through tidal disruption events and
curve, while the vertical dashed lines mark n = 50, 66, 255, galactic mergers (Rodriguez-Gomez et al. 2017), which are
during the process of the snapshot torus wrapping around violent processes taking place on relatively short timescales.
smaller islands. In contrast, pseudobulges are thought to form and evolve
through secular processes (Kormendy & Kennicutt 2004;
Tonini et al. 2016), contributing to the slow buildup and re-
in the effective potential Eq. (5) controlling mass transfer shaping of pseudobulges, in which continuous mass transfer
within the galaxy. between galactic components may also play a role.
The general aim of this paper was to give, through the
4. Conclusions presented example, a universal tool for the description of
problems arising in galactic and celestial mechanics where a
In this paper we have numerically studied ensembles of parameter drift is present. The method originates in chaos
stellar obits for the Milky Way with different potential pro- theory, and is a currently actively developed area of re-
files for the galactic components: the extended Myamoto- search. Some other concepts that could be useful from this
Nagai potential for the disk, the Hernquist potential for method include time-dependent stable and unstable mani-
the bulge and the Navarro-Frenk-White potential for the folds and foliations (Jánosi & Tél 2024), doubly transient
dark matter halo. We simulated a mass transfer from the chaos (Motter et al. 2013), as well as snapshot attractors
disk to the bulge by heuristically introducing a decreasing in dissipative cases (Jánosi & Tél 2019; Jánosi et al. 2021;
and increasing parameter drift of the same rate in the mass Jánosi & Tél 2022; Jánosi & Tél 2024). Some further meth-
parameters of the disk and the bulge, respectively. For sim- ods that could be used for augmenting the presented meth-
plicity, for the form of the drift we chose a linear function, ods include e.g. developing new ways of phase space pro-
but in general cases any monotonous form could be used, jections dependent on trajectory energies, or by expanding
obtained from the observations. The chosen rate results in the EAPD-analysis with mapping methods. Furthermore,
30% of the disk mass transferred to the bulge over time it can be also beneficial to expand the analysis into more
scales which can be deemed reasonably long compared to composite and realistic N-body simulation models, where
the time scales of our Universe (Spergel et al. 1997). the axisymmetric bulge and dark matter halo components
The stellar orbits were illustrated on Poincaré surface and the non-axisymmetric spiral arms are also included in
of sections in the case of constant mass, and were shown the potential (Barros et al. (2016); Pouliasis et al. (2016);
to exhibit a mixed phase space with coexisting periodic, Trick et al. (2017); Michtchenko et al. (2018)).
Article number, page 10 of 13
Eduárd Illés et al.: Orbital dynamics in galactic potentials under mass transfer

Acknowledgements. We are grateful to Tamás Tél and Santi Roca- Malhotra, R. 1998, in Astronomical Society of the Pacific Conference
Fàbrega for their valuable comments, which significantly improved the Series, Vol. 149, Solar System Formation and Evolution, ed. D. Laz-
interpretation of the results. This work was supported by the Hun- zaro, R. Vieira Martins, S. Ferraz-Mello, & J. Fernandez, 37
garian National Research, Development and Innovation Office, under Manos, T., Bountis, T., & Skokos, C. 2013, Journal of Physics A
Grants No. KDP-2023 C2262591, TKP–2021 BME-NVA-02 (D.J.) and Mathematical General, 46, 254017
TKP2021-NKTA-64 (T.K.), financed by the Ministry of Culture and Manos, T. & Machado, R. E. G. 2014, MNRAS, 438, 2201
Innovation of Hungary. McMillan, P. J. 2011, Monthly Notices of the Royal Astronomical
Society, 414, 2446
McMillan, P. J. 2016, Monthly Notices of the Royal Astronomical
Society, 465, 76–94
References Meiss, J. D., Cary, J. R., Grebogi, C., et al. 1983, Physica D: Nonlinear
Phenomena, 6, 375
Altmann, E. G., Portela, J. S. E., & Tél, T. 2013, Reviews of Modern Michtchenko, T., Lepine, J., Barros, D., & S. S. Vieira, R. 2018, As-
Physics, 85, 869 tronomy I& Astrophysics, 615
Arnold, B., Wright, N., & Parker, R. 2022, Monthly Notices of the Miyamoto, M. & Nagai, R. 1975, Publications of the Astronomical
Royal Astronomical Society, 515 Society of Japan, 27, 533
Bahcall, J. N. & Soneira, R. M. 1980, The Astrophysical Journal Let- Motter, A. E., Gruiz, M., Károlyi, G., & Tél, T. 2013, Phys. Rev. Lett.,
ters, 238, L17 111, 194101
Barros, D. A., Lépine, J. R. D., & Dias, W. S. 2016, A&A, 593, A108 Mróz, P., Udalski, A., Skowron, D. M., et al. 2019, ApJ, 870, L10
Benettin, G., Galgani, L., Giorgilli, A., & Strelcyn, J. M. 1980, Mec- Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563
canica, 15, 9 Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2015
Bernacca, P. L. 1985, Ap&SS, 110, 21 Ott, E. 1993, Chaos in dynamical systems, 1st edn. (Cambridge, UK:
Binney, J. & Spergel, D. 1982, ApJ, 252, 308 Cambridge University Press)
Binney, J. & Tremaine, S. 2008, Galactic Dynamics: Second Edition Ou, X., Eilers, A.-C., Necib, L., & Frebel, A. 2024, Monthly Notices
(Princeton, USA: Princeton Univ. Press) of the Royal Astronomical Society, 528, 693
Caldwell, J. A. R. & Ostriker, J. P. 1981, The Astrophysical Journal Petzold, L. 1983, SIAM Journal on Scientific and Statistical Comput-
Letters, 251, 61 ing, 4, 136
Caranicolas, N. D. 1996, Ap&SS, 246, 15 Pichardo, B., Martos, M., Moreno, E., & Espresate, J. 2003, ApJ, 582,
Caranicolas, N. D. & Papadopoulos, N. J. 2003, A&A, 399, 957 230
Channon, S. R. & Lebowitz, J. L. 1980, Annals of the New York Posti, L. & Helmi, A. 2019, A&A, 621, A56
Academy of Sciences, 357, 108 Pouliasis, E., Matteo, P., & Haywood, M. 2016, Astronomy I& Astro-
Chekroun, M. D., Simonnet, E., & Ghil, M. 2011, Physica D Nonlinear physics, 598
Phenomena, 240, 1685 Rieder, S., Dobbs, C., Bending, T., Liow, K. Y., & Wurster, J. 2021,
Cincotta, P. M., Nunez, J. A., & Muzzio, J. C. 1996, ApJ, 456, 274 Monthly Notices of the Royal Astronomical Society, 509
Rodriguez-Gomez, V., Sales, L. V., Genel, S., et al. 2017, Monthly
Cincotta, P. M. & Simó, C. 2000, A&AS, 147, 205
Notices of the Royal Astronomical Society, 467, 3083
Contopoulos, G. 2002, Order and chaos in dynamical astronomy
Sanderson, R. E., Hartke, J., & Helmi, A. 2017, The Astrophysical
(Berlin, Germany: Springer)
Journal, 836, 234
Contopoulos, G. & Seimenis, J. 1990, A&A, 227, 49
Schmidt, M. 1956, Bulletin of the Astronomical Institutes of the
Cristadoro, G. & Ketzmerick, R. 2008, Phys. Rev. Lett., 100, 184101
Netherlands, 13, 15
Daffern-Powell, E. & Parker, R. 2020, Monthly Notices of the Royal Skokos, C. 2010, in Lecture Notes in Physics, Berlin Springer Verlag,
Astronomical Society, 493, 4925 ed. J. Souchay & R. Dvorak, Vol. 790, 63–135
Dembo, R. S. & Steihaug, T. 1983, Mathematical Programming, 26, Sofue, Y. 2012, Publications of the Astronomical Society of Japan, 64,
190 75
Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Sofue, Y. 2013, 985
Publications of the Astronomical Society of the Pacific, 125, Sommerer, J. C. & Ott, E. 1993, Science, 259, 335
306–312 Spergel, D. N., Bolte, M., & Freedman, W. 1997, Proceedings of the
Froeschlé, C., Lega, E., & Gonczi, R. 1997, Celestial Mechanics and National Academy of Sciences, 94, 6579
Dynamical Astronomy, 67, 41 Teschl, G. 2012, Ordinary Differential Equations and Dynamical Sys-
Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, tems, Graduate studies in mathematics (American Mathematical
595, A1 Soc.)
González, G. A., Plata-Plata, S. M., & Ramos-Caro, J. 2010, MNRAS, Tonini, C., Mutch, S. J., Croton, D. J., & Wyithe, J. S. B. 2016,
404, 468 Monthly Notices of the Royal Astronomical Society, 459, 4109
Goodman, J. & Weare, J. 2010, Communications in Applied Mathe- Toomre, A. 1963, The Astrophysical Journal, 138, 385
matics and Computational Science, 5, 65 Trick, W., Bovy, J., D’Onghia, E., & Rix, H.-W. 2017, The Astro-
Gottwald, G. A. & Melbourne, I. 2004, Proceedings of the Royal So- physical Journal, 839
ciety of London Series A, 460, 603 Verheijen, M. A. W. & Sancisi, R. 2001, Astronomy & Astrophysics,
GRAVITY Collaboration, Abuter, R., Amorim, A., et al. 2021, A&A, 370, 765
647, A59 Vogt, D. & Letelier, P. S. 2005, Publications of the Astronomical So-
Henon, M. & Heiles, C. 1964, AJ, 69, 73 ciety of Japan, 57, 871
Hernquist, L. 1990, ApJ, 356, 359 Walder, M., Erkal, D., Collins, M., & Martinez-Delgado, D. 2024,
Jánosi, D., Károlyi, G., & Tél, T. 2021, Nonlinear Dynamics, 106, Probing the dark matter haloes of external galaxies with stellar
2781–2805 streams
Jánosi, D. & Tél, T. 2019, Chaos, 29, 121105 Wogan, N. 2021, NumbaLSODA, https://github.com/Nicholaswogan/numbalsod
Jánosi, D. & Tél, T. 2021, Chaos, 31, 033142 Yu, L., Ott, E., & Chen, Q. 1990, Phys. Rev. Lett., 65, 2935
Jánosi, D. & Tél, T. 2022, Physical Review E, 105, L062202 Zotos, E. E. 2011, Chaos Solitons and Fractals, 44, 501
Jánosi, D. & Tél, T. 2024, Physics Reports, 1092, 1 Zotos, E. E. 2012, New A, 17, 576
Karney, C. F. 1983, Physica D: Nonlinear Phenomena, 8, 360 Zotos, E. E. 2014, Baltic Astronomy, 23, 151
Katz, D., Antoja, T., Romero-Gómez, M., et al. 2018, Astronomy I& Zotos, E. E., Dubeibe, F. L., Steklain, A. F., & Saeed, T. 2020, A&A,
Astrophysics, 616, A11 643, A33
Kormendy, J. & Kennicutt, Robert C., J. 2004, ARA&A, 42, 603
Kovács, T. 2020, J. R. Soc. Interface., 17, 20200648
Kuzmin, G. 1956, Astron. Zh., 33, 27
Labini, F. S., Žofia Chrobáková, Capuzzo-Dolcetta, R., & López-
Corredoira, M. 2023, The Astrophysical Journal, 945, 3
Laskar, J. 1990, Icarus, 88, 266
Lichtenberg, A. & Lieberman, M. 1992, Regular and Chaotic Dynam-
ics, 1st edn. (New York, NY, USA: Springer)

Article number, page 11 of 13

A&A proofs: manuscript no. revamp

Appendix A: Illustrating another example of

peculiar dynamics
Figure A.1 shows the phase space plots in the (R, pR )
plane, where we illustrate the phase portrait for four differ-
ent values of Lz with the calculated H0 value.

Fig. A.1: Phase portraits of Eq. (10) for different values of

Lz with constant H = H0 . The trajectories are plotted on
the (R, pR ) plane, following the Poincaré cut at z = 0, and
the value of pz obtained from (12), indicated at each point
by the colormap. The red KAM torus belongs to the initial
condition R0 = 0.67, pR 0 = 129.79 (panel a).

A case of even more peculiar dynamics is shown in

Fig. A.2 by following the snapshot torus shown in red in
Fig. A.1a, similarly to the previous two cases in the main
text. In this case an ensemble of orbits is chosen which have
significantly lower angular momentum Lz , hence also hav-
ing significantly low radial velocities compared to the cases
presented before. The first stage of its evolution is similar
to the torus presented in Fig. 4, where the torus develops
stretches, which results in its deformation. However, the
torus does not break up the same way as the one in Fig. 4
does, this time the deformation begins from several small
tongues on the left side of the torus at n = 61. These then
transport points into the snapshot chaotic sea, but, at the
same time, four tongues can be observed to penetrate into
the torus and wrapping around areas that were originally
inside it, see instants n = 80, 98. These structures are sim-
ilar to the islands seen in Fig.6, disappearing by n = 250. Fig. A.2: Time evolution of the red torus on Fig. A.1a.
However, even at this point a significant part of the snap- At n = 61 we can see the initial torus starts to deform
shot torus remains intact, with a large number of points significantly, while a small fraction of trajectories scatter
from the original structure still bordering this new shape, out into the chaotic sea after n = 80. After n = 250 only an
similarly to the case of Fig.6. This property remains true inner structure remains, still separating snapshot tori from
until the end of the simulation at n = 1069. The behav- the snapshot chaotic sea.
ior of the peculiar structures penetrating inside the large
island, present in the snapshots of n = 80, 98, are not yet
ders of magnitude smaller than unity, consistent with the
understood and further, more extended research is needed.
observations of Fig.A.2 at n = 620, 1069.
The EAPD curve for this torus is shown in Fig. A.3, In Fig. A.4, a monotonous growth in energy dispersion
where we see an initial exponential increase from n = can be seen, similarly to the case in Fig. 6, and after n = 51
61, but after the points finished scattering into the outer the mean curve in energy shows a slight break.
chaotic sea the distance goes into saturation at a value or-
Article number, page 12 of 13
Eduárd Illés et al.: Orbital dynamics in galactic potentials under mass transfer

Fig. A.3: EAPD curve of Eq. (14) evaluated for the snapshot
torus of Fig. A.2. The inset shows an oscillating behavior of
small amplitudes after n = 150. A breakup phase is present
after n = 61, for which the Lyapunov exponent obtained
from the fitting is λnr = 0.211 ± 0.003.

Fig. A.4: Median energy of the ensemble of the snapshot

torus shown in Fig. A.2 compared to the EAPD curve of
Fig. A.3. The 1σ dispersion is shown as a light red band
around the energy curve, the vertical dashed lines mark
n = 51, 620.

Article number, page 13 of 13