© The Authors 2025

1 Introduction

A popular method to identify debris from past mergers in our Galaxy is via clustering in the space of energy E, and two components of the angular momentum, L_z and L_⊥ (Helmi & de Zeeuw 2000). The argument behind the method first put forward by these authors is that in an integrable, static potential, these quantities are conserved, and therefore a satellite system that is initially compact in position and velocity space, will also be compact in (E, L_z, L_⊥), and remain so even after the system has fully disrupted and phase-mixed. Thus far clustering analyses have been applied relatively successfully to samples of stars near the Sun (see the reviews by Helmi 2020; Deason & Belokurov 2024), where our Galaxy has been assumed to be axisymmetric and (very close to) equilibrium. In this case, most stars would have orbits with at least two integrals, E and L_z, and while the total angular momentum, and hence L_⊥, is not an integral, it does not vary wildly. Time dependence due to cosmological growth does not seem to have a very significant effect on these simple assumptions (Gómez et al. 2013), although the debris’ complexity is larger, especially for more massive progenitors (e.g. Khoperskov et al. 2023).

There is, however, an additional time dependence and deviation from axisymmetry that is often neglected in such analyses, and that is the presence of the Galactic bar. Whereas its effect on the kinematics of disk stars has been well-established since HIPPARCOS (e.g. Dehnen 1999), its impact on the halo has only recently become clear (Dillamore et al. 2023). Through time-dependence, energy is not conserved, and furthermore the bar introduces resonances, which can induce kinematic substructure (Dillamore et al. 2024). The bar’s evident triaxial shape implies that none of the angular moment components are conserved, which may introduce chaoticity (Manos & Athanassoula 2011). This will obviously be important for (i) stars in the disk and (ii) stars on strongly radial orbits that penetrate the innermost regions of the Galaxy.

While it has been widely known that the stellar halo velocity ellipsoid is radially anisotropic, the extent of it became really apparent only with Gaia data (Gaia Collaboration 2018b). The halo velocity distribution has a significant population of stars that are on extremely radial orbits, the Gaia-Sausage (Belokurov et al. 2018), which is associated with the Gaia-Enceladus merger (Helmi et al. 2018). These stars should be strongly affected by the bar. Stars on the red sequence, (the hot thick disk, or splash population, see Koppelman et al. 2018; Belokurov et al. 2019, respectively) should also be affected by the bar, but differently due to their differing orbital properties.

Here we show that the influence of the bar on the stellar halo near the Sun cannot be neglected. We demonstrate that a large fraction of the halo stars are on chaotic orbits and that therefore the use of E, L_z, and L_⊥ needs some reassessment, particularly for stars on radial and/or very bound orbits. By integrating the orbits of a local sample of halo stars, we studied both the degree and effects of chaoticity introduced by a rotating bar on their distribution in (E, L_z, L_⊥) space.

This paper is structured as follows. Sect. 2 introduces the data and Galactic potential used. In Sect. 3, we describe the orbital characterisation, which includes orbital frequency analysis and the application of chaos indicators. In Sect. 4, we present the results of our analyses. We end with a general discussion in Sect. 5 and present our conclusions in Sect. 6.

2 Data and method

2.1 Data

We used a sample of stars extracted from the Gaia DR3 RVS set (Gaia Collaboration 2022). The distances to the stars were obtained by inverting their parallaxes after correcting them for a zeropoint offset following Lindegren et al. (2021), and retaining only those for which 𝟉_corr/ε_𝟉 > 5. We corrected for extinction using the Lallement et al. (2022) dust extinction maps. We also removed likely globular cluster (GC) members (with probabilities >90%) from the sample using the membership criteria derived by Vasiliev & Baumgardt (2021). See Dodd et al. (2023) for more details on the quality criteria applied.

To convert the observables to Galactocentric Cartesian and cylindrical coordinates¹ we also followed Dodd et al. (2023). Hence, we assumed a right-handed coordinated system in which the Sun is located in the mid-plane (z_☉ = 0 kpc) at x_☉ = −8.2 kpc from the Galactic centre (McMillan 2017). To correct for the motion of the Sun, we used (U, V, W)_☉ = (11.1, 12.24, 7.25) km s⁻¹ (Schönrich et al. 2010) and |V_LSR| = 232.8 km s⁻¹ (McMillan 2017). The z-component of the angular momentum vector, L_z, is defined such that it is positive for prograde orbits, meaning that its sign is flipped.

To obtain a halo sample we selected stars using the Toomre cut |V−V_{LS R}| > 210 km s⁻¹ and considered only those stars within a distance of 1 kpc from the Sun, and with total velocities smaller than 525 km s⁻¹. This velocity cut (well below the estimated local escape velocity, e.g. Koppelman & Helmi 2021) ensures that all stars are bound in our Galactic potential model (as E = E_pot+E_kin ∼−5000 km² s⁻² for this velocity). This leaves us with a sample containing 27 885 stars.

2.2 Milky Way model

We used AGAMA (Vasiliev 2019) for orbit integrations and to compute dynamical quantities, and employed the AGAMA implementation of a barred Galactic potential model as described in Hunter et al. (2024)². This is a basis function expansion (BFE) of the analytic model, which allows a higher computational efficiency. This is our fiducial potential. It uses the Sormani et al. (2022) bar model, which is the analytic approximation of the bar model from Portail et al. (2016). It consists of an X-shaped boxy bulge/bar component, a short bar, and a long bar, and has a total mass of M_bar = 1.83 · 10¹⁰ M_☉. Further, it employs the gas disks from McMillan (2017) and an exponential thin and thick stellar disk with a density dip in the centre with parameters that were fit to the Eilers et al. (2019) and Mróz et al. (2019) circular velocity curves. It includes a central black hole, nuclear star cluster, nuclear stellar disk, and a spherical Einasto (1969) dark matter halo with a mass of M = 1.1 · 10¹² M_☉. We did not include a potential describing the spiral arms in this work. We set the bar pattern speed equal to Ω_b = −37.5 km s⁻¹ kpc⁻¹ and assume a bar angle ϕ_b = −25 degrees, in agreement with recent estimates (see Hunt & Vasiliev 2025).

Since the bar is rotating, the Galactic potential is time dependent, and hence a star’s total energy is not conserved along its orbit. The Jacobi integral or Jacobi energy E_J may be used instead as it takes this rotation into account. E_J is defined as $$E_J=\frac{1}{2}(v_{x,\mathrm{CR}}^2+v_{y,\mathrm{CR}}^2+v_{z,\mathrm{CR}}^2)+\Phi (x,y,z,t)-\frac{1}{2}{\Omega }_b^2{R}^2,$$(1) where the subscript CR denotes that we are considering quantities in the corotating frame. Positions and velocities in the corotating frame are obtained using $$\begin{array}{cc}& {\mathit{x}}_{\mathrm{CR}}=M_R·\mathit{x}\hfill \\ & {\mathit{v}}_{\mathrm{CR}}=M_R·\mathit{v}-{\mathit{\Omega }}_{\mathit{b}}\times {\mathit{x}}_{\mathrm{CR}}\hfill \end{array}$$(2) with $$M_R=\left(\begin{array}{ccc}\cos (\varphi )& -\sin (\varphi )& 0\\ \sin (\varphi )& \cos (\varphi )& 0\\ 0& 0& 1\end{array}\right),$$(3) Ω_b = (0, 0, Ω_b), and ϕ = −(ϕ_b+Ω_b t). Hence, E_J can also be expressed as E_J = E−Ω_b L_z, which shows that its sign does not indicate whether a star is bound or not.

To ensure that the Jacobi integral is conserved up to at least the order of 10⁻⁵ over the integration time we increased the accuracy parameter in AGAMA if necessary. The limiting factor for the conservation of the Jacobi energy in our integrations is the finite smoothness of the BFE interpolator.

We defined the orbital period for each stars, T_orb,r, as the median time between consecutive apocentres over an integration time of 30 Gyr. The median orbital period for our sample is 120 Myr, with 5% of the stars having orbital periods smaller than 90 Myr and 5% having orbital periods larger than 320 Myr, demonstrating that the stars in our sample have a large range of orbital periods. To make sure that in our analysis each orbit is sampled at a high enough time frequency Δ t for a long enough period of integration time T_int, we split our sample in three subsets:

• regime 1: T_orb,r<0.2 Gyr, Δ t = 0.001 Gyr, T_int = 15 Gyr, 17 741 stars;

• regime 2: 0.2<T_orb,r<0.5 Gyr, Δ t = 0.003 Gyr, T_int = 50 Gyr, 9192 stars;

• regime 3: 0.5 Gyr<T_orb,r, Δ t = 0.005 Gyr, T_int = 100 Gyr, 951 stars.

3 Orbital characterisation: frequency analysis, bar resonances, and chaos indicators

A regular orbit has three integrals of motion (IoM). In the case of (static) axisymmetric systems, the total energy E, L_z, and I₃ are used to denote these integrals, although I₃ is often not known explicitly. Therefore, $L_{\perp }=\sqrt{L_x^2+L_y^2}$ is used as a proxy instead, even though this is not a truly conserved quantity. In a time-dependent triaxial potential such as the one we use here to represent the Galaxy, none of these quantities are conserved. This is illustrated in Fig. 1, which shows (E, L_z) and (L_z, L_⊥) computed at t = 0 and evaluated after 10 Gyr of integration in the fiducial potential, in the top and bottom panels, respectively. Although the distribution of stars globally does not change dramatically³, we do see especially large differences for highly bound, low angular momentum orbits. In particular, a fraction of the orbits which initially had E<−13000 km² s⁻², have become more bound and mildly retrograde.

Arguably a regular orbit can be best characterised using its three fundamental orbital frequencies, which are constant in time. Generally these cannot be computed analytically but may instead be found by performing a Fourier transform of the complex time series of the orbit’s phase-space coordinates x(t)+i v(t). In this work, we considered the orbital frequencies in the corotating frame in Cartesian coordinates (i.e. Ω_{x, CR}, Ω_{y, CR}, Ω_{z, CR}) and in a variation of cylindrical coordinates (i.e..Ω_{R, CR}, Ω_{ϕ, CR}, Ω_{z, CR})⁴. We used Poincaré symplectic coordinates (Papaphilippou & Laskar 1996; Valluri et al. 2012), in which case the complex time series are of the form $R+i{v}_R,\sqrt{2L_z}[\cos (\varphi )+i\sin (\varphi )]$ and z+i v_z. To determine the orbital frequencies, we used SuperFreq (Price-Whelan 2015) with parameters min_freq_diff = 10^**−6, nintvec =15, and break_condition = None.

If two or three of the orbit’s fundamental frequencies are a small integer ratio of one another, the corresponding orbit is said to be resonant. In a rotating barred potential, an orbit can also be in resonance with the bar, and this occurs when $$m({\Omega }_{\varphi }-{\Omega }_b)+l{\Omega }_R=0$$(4) in the inertial frame, or $$m\left({\Omega }_{\varphi ,\mathrm{CR}}\right)+l{\Omega }_{R,\mathrm{CR}}=0$$(5) in the corotating frame, with l and m being an integer and where Ω_b is the bar’s pattern speed. l/m = 0 corresponds to the corotation resonance (an orbit has the same azimuthal orbital frequency as the bar), l/m = −1/2 corresponds to the inner Lindblad resonance (ILR), and l/m = 1/2 corresponds to the outer Lindblad resonance (OLR).

Fig. 2 shows examples of orbits in resonance with the bar. These were selected using $\frac{{\Omega }_{\varphi ,\mathrm{CR}}}{{\Omega }_{R,\mathrm{CR}}}=l/m\pm 0.005$ for a range of integers l and m. Their present-day L_⊥ and L_z values range from 0 to 900 kpc km s⁻¹, and from −1000 to 1500 kpc km s⁻¹, respectively. This figure shows that even orbits with large apocentres (e.g. 20 or 30 kpc) can be in resonance with the bar.

Chaotic orbits have fewer than three IoM and their orbital frequencies thus change over time. An example of a chaotic orbit is shown in Fig. 3. After 10 Gyr of integration time, this orbit seems to be librating around the corotation resonance, and it changed its apocentre but also E, L_z, and L_⊥ by a substantial amount⁵. Clearly, in the case of this orbit, the bar is able to cause chaotic behaviour that manifests itself after only a few Gyr.

Fig. 1 Distribution of nearby halo stars in E−Lz (top row) and L⊥−Lz (bottom row) at t = 0 and after 10 Gyr of evolution in our fiducial Milky Way potential including a rotating bar. Especially highly bound stars with Lz ∼ 0 (and small L⊥) change their distribution significantly.

Fig. 2 Examples of orbits of stars on bar resonance, identified as having $\frac{{\Omega }_{\varphi ,\mathrm{CR}}}{{\Omega }_{R,\mathrm{CR}}}=l/m\pm 0.005$ and λ = 0 and fdr<−1.9. The orbits are shown in the corotating frame, and the resonance they are associated to is indicated in the top right corner. Isodensity contours of the fiducial potential in the mid-plane (z = 0) are shown in the background, the black cross indicates the Galactic Centre.

Fig. 3 Example of a star on a chaotic orbit (selected from regime1). The red and blue markers show the position of the star at the present day and after 10 Gyr of integration time, respectively. The top panels show projections of the orbit in the corotating frame, (xCR, yCR), and in (R, z), respectively. The bottom panels illustrate the trajectory of this star in (E, Lz, L⊥) space over time. This star ends up librating around the bar’s corotation resonance after 10 Gyr, and moves around significantly in (E, Lz, L⊥) space as it changes its orbit. The orbit has λ = 2.1 Gyr−1 and tλ = 1.3 Gyr, and fdr=−0.2 (see Sect. 3.1 and 3.2 for details on how these quantities are computed).

3.1 Chaos indicators: Frequency diffusion rate

To quantify the degree and effect of chaoticity in our halo sample, we made use of the Lyapunov exponent λ, and we calculated the frequency diffusion rate, which we denote as fdr. These are different but complementary measures of chaoticity. To obtain a chaoticity timescale we used the Lyapunov exponent. This timescale gives an indication of how quickly an orbit changes due to chaos.

The fdr is based on the fact that the orbital frequencies change over time for a chaotic orbit. To compute it, we took the difference in orbital frequency between two consecutive periods of equal integration time, 0<t₁<T_int and T_int<t₂<2 T_int, and determined $$\Delta f_i=\left|\frac{{\Omega }_i\left(t_1\right)-{\Omega }_i\left(t_2\right)}{{\Omega }_i\left(t_1\right)}\right|$$(6) for each of the three frequencies Ω_i. We took the geometric average of the three Δ f_i and then its logarithm to obtain $$\mathtt{f}\mathtt{d}\mathtt{r}={\log }_{10}\mathrm{\Delta }f,\text{\hspace{0.17em}where\hspace{0.17em}}\mathrm{\Delta }f={\left[\prod _i\mathrm{\Delta }{f}_i\right]}^{\frac{1}{3}}.$$(7) The larger the value of the fdr, the more chaotic the orbit is. However, in general there is no clear threshold value of the fdr that separates chaotic from regular orbits, and instead, its spectrum is continuous (Valluri et al. 2010). Nonetheless, we find that the fdr correlates well with the Lyapunov exponent λ (see also Vasiliev 2013, and Sect. 4.1). An important caveat is that for resonant orbits the frequencies are not independent, and consequently the determination of the fdr becomes less reliable in such cases.

3.2 Chaos indicators: Lyapunov exponent

The Lyapunov exponent measures the evolution in time of the infinitesimal 6D deviation vector w(t) from a given orbit. The divergence of nearby orbits in configuration space grows linearly in time for a regular orbit (Helmi et al. 1999; Helmi & Gomez 2007; Vogelsberger et al. 2008), and exponentially in the chaotic case (Maffione et al. 2015; Price-Whelan et al. 2016a), and therefore we expect a similar behaviour for the deviation vector. In this work, we used the set of six deviation vectors computed by AGAMA to derive the Lyapunov exponent and the onset time⁶ of the exponential growth regime, t_onset. To this end, we took the norm of each deviation vector per timestep, and selected the vector with the largest norm over time, max (∥w₁(t)∥, …,∥w₆(t)∥). We then performed a non-linear least squares fit using scipy.curve_fit of the following function, $$\ln \left(\frac{\parallel w_{max}(t)\parallel }{t}\right)=\{\begin{array}{ll}C& \text{\hspace{0.17em}for\hspace{0.17em}}t<t_{\text{onset\hspace{0.17em}}}\\ C+\lambda (t-t_{\text{onset\hspace{0.17em}}})& \text{\hspace{0.17em}for\hspace{0.17em}}t>t_{\text{onset\hspace{0.17em}}},\end{array}$$(8) where C is a constant. If there is no exponential growth, λ = 0 and t_onset is set to be equal to the integration time, while if there is exponential growth, λ is returned in units of Gyr⁻¹ and t_onset will be smaller than the integration time and zero if chaotic behaviour starts immediately. To make sure chaoticity is robustly detected, we set λ = 0 if T_int−t_onset<5 T_orb, implying that the time of exponential growth is detected over a period shorter than five times the orbital timescale, or if (T_int−t_onset) λ<0.1 ln(T_int), that is if the increase in the norm of the deviation vector due to chaoticity is smaller than 10% of the growth expected from the divergence for a regular orbit over the integration time. We checked that the method implemented via Eq. (8) gives similar values for λ as when computed using an independent method (Voglis et al. 2002), and that the derived magnitude of t_onset roughly agrees with when the orbit starts to depict exponential divergence. To ensure that t_onset is reliably determined within a Hubble time, we integrated stars in the subset defined as regime1 for 25 Gyr when computing their deviation vectors.

We defined a timescale of chaoticity t_λ, using both the Lyapunov exponent and the chaos onset time as $$t_{\lambda }=\frac{1}{\lambda }+t_{\text{onset}}.$$(9) We generally find that t_λ gives a lower limit to the timescale over which chaoticity is manifested macroscopically (i.e. a significant change in the orbit that is apparent through visual inspection) and that typically 5/λ+t_onset better reflects the time over which chaoticity becomes visually discernible.

Fig. 4 Distribution of the values of the two chaoticity indicators used in this work for our local halo sample: the Lyapunov exponent λ, normalised by Torb,r, and the fdr. The different panels correspond to three different versions of the potential described in Section 2.1: with a rotating bar (left), with a static bar (middle), and an axisymmetrised bar (right). The side panels show the cumulative distributions of normalised λ for each of the potentials, while the top panels correspond to those of fdr. The fdr distribution for the orbits that have λ = 0 is shown separately. The vertical magenta lines correspond to fdr=−1.9, the 99 th percentile of the orbits with λ = 0 in the fiducial potential.

4 Results

We present here the analysis of the orbits of our local sample of halo stars. We first focus on the distribution of their Lyapunov exponents and fdr, and then identify where in the present-day’s (E, L_z, L_⊥) space and for what types of orbits the effect of chaoticity is largest. We also locate the bar resonances and determine the fraction of stars trapped on or affected by these. Finally we discuss how substructures identified in previous work in the present-day’s space of (E, L_z, L_⊥) may be affected by the bar.

4.1 Degree of chaoticity for halo stars near the Sun

Figure 4 shows the distribution of the values of the two chaoticity indicators introduced in the previous section for our local halo sample: the Lyapunov exponents and the fdr. We show the Lyapunov exponents normalised by T_orb,r to allow for a more direct comparison between orbits with different orbital timescales. The indicators were computed for three variations of the potential described in Section 2.1: the fiducial potential itself, which is triaxial and time-dependent (left panels), a model with a static bar (middle panels), and an axisymmetrised model of the fiducial potential, meaning that in the basis function expansion only the m = 0 term is retained, making it also time-independent (right panels).

As expected, the largest number of stars on regular orbits are found for the axisymmetric potential, where 56% of the stars have λ = 0. A much smaller fraction of the stars have λ = 0 for the rotating and static bar potential, 4% and 8%, respectively. The distributions for these two triaxial cases are relatively similar, with the fdr increasing sharply for fdr >−2 (generally corresponding to normalised λ > 0.1). The similarity of these distributions indicates that the chaoticity is driven mostly by the triaxiality present in the potential and to a smaller degree by the resonances with the rotating bar.

To obtain a robust estimate of the amount of stars on chaotic orbits, we preferred to use both the fdr and Lyapunov exponent. To identify such orbits we required λ > 0 and fdr > −1.9, since this value corresponds to the 99th percentile of the distribution of fdr for orbits that have λ = 0 in the fiducial potential⁷, as can be inferred from the top sub-panels of Fig. 4. With these criteria we find that 60% of the stars in our sample are on chaotic orbits in the rotating bar potential, and a comparable percentage, namely 56%, in the static bar potential. Blind application of this criterion to the axisymmetrised potential would yield that only 7% of the stars are on chaotic orbits. The right panel of Fig. 4 shows this would be an incorrect assessment, because for orbits with λ = 0 the 99 th percentile takes a lower fdr value, namely −3. In fact, by inspecting the various indicators and the chaos onset times we find that the region with λ<0.15 and fdr<−1.9 hosts sticky chaotic orbits with large chaos onset times, which comprise 13% of the stars. Stars on fully chaotic orbits would be those with (λ > 0.15 and −3<fdr<−1.9) or (λ > 0 and fdr>−1.9), and comprise 26% of the sample.

4.2 Chaoticity in the present-day space of E, L_z, and L_⊥

We investigated which region of present-day (E, L_z, L_⊥) space and hence what types of orbits for the stars in our sample are most affected by the presence of a (rotating) bar. Fig. 5 shows the percentage of stars that have Lyapunov exponents higher than 0 and fdr > −1.9 per bin in (E, L_z, L_⊥) space. Clearly, stars on radial orbits, having L_z ∼ 0 kpc km s⁻¹, and/or strongly bound orbits are for the large majority chaotic. About 70% of the stars in our halo sample with |L_z|<500 kpc km s⁻¹ have chaotic orbits, and this percentage increases by a few percent for stars on more bound orbits (E<−140 000 km² s⁻²). Instead, stars on the retrograde side of the halo, having L_z ≤−1000 kpc km s⁻¹, or stars with high inclination orbits, having L_⊥ ≥ 1500 kpc km s⁻¹, are barely affected by the bar.

We explored how the Lyapunov timescale t_λ (Eq. (9)) varies as a function of the present day location of stars in (E, L_z, L_⊥) space. The top row of Fig. 6 shows the median t_λ per bin in (E, L_z, L_⊥) space, while the bottom row shows the median fdr. The two chaos indicators show good agreement in terms of structure in (E, L_z, L_⊥) space, indicating that chaoticity is detected robustly, as the two measures are independent. The top panels show that there are regions where chaos acts on a timescale (much) shorter than the Hubble time, especially for orbits that are strongly bound and/or radial, which can significantly change even within a couple of Gyr or less.

These short timescales are accompanied by a large change in E, L_z, and L_⊥, as also illustrated in Figs. 1 and 3. Fig. 3 also illustrates the usual path that an initially relatively strongly bound, very radial orbit takes over time in (E, L_z, L_⊥) space, namely, along a line of positive slope in (E, L_z) and along a horizontal line with a large amount of variation in (L_z, L_⊥), generally reaching L_⊥ = 0 kpc km s⁻¹ (see also Dillamore et al. 2024). Due to this, within the timespan of a few Gyrs, a star’s orbit can change from being prograde to retrograde (change in L_z), from being confined to the disk to reaching values of z_max of a few kpc (as a result of a change in L_⊥) and from passing through the local volume to not being able to reach that volume at all because the apocentre has shrunk (change in E). As Figs. 5 and 6 show, this mixing in (E, L_z, L_⊥) space is particularly strong for very bound, relatively radial orbits. This has implications for the reliability of the identification of substructure associated to merger debris, and it also implies that substructures located today in such regions of this space could suffer significant contamination. We return to this point in Sect. 4.4.

Fig. 5 Local halo stars in (E, Lz, L⊥) space computed at t = 0 in the fiducial potential, binned into 100 bins along each axis, where the colour coding of each bin shows the percentage of stars with a Lyapunov exponent higher than 0 and a fdr >−1.9, determined after 25 Gyr of integration.

Fig. 6 Similar to Fig. 5 but showing in the top panel the weighted median Lyapunov timescale (Eq. (9)), and in the bottom panel the median fdr (Eq. (7)) determined over an integration time of 25 Gyr for all stars in the sample. We note the good correspondence between both panels, and with Fig. 5.

Fig. 7 Position of local halo stars close to bar resonances in (E, Lz, L⊥) space at the present day, selected to have $\frac{{\Omega }_{\varphi ,\mathrm{CR}}}{{\Omega }_{R,\mathrm{CR}}}=l/m\pm 0.005$ (see Fig. 2 for examples of each resonant orbit). The stars encircled in black are on resonant bar orbits, meaning with λ = 0 and fdr <−1.9.

4.3 The presence of resonances

The stripy diagonal structures seen in Figs. 5 and 6 can be linked to the bar resonances that are present in this region. This is illustrated in Fig. 7, which shows stars that were selected to have $\frac{{\Omega }_{\varphi ,\mathrm{CR}}}{{\Omega }_{R,\mathrm{CR}}}=l/m\pm 0.005$ for a range of l and m. Approximately 13% of the stars in our sample are sufficiently close to a resonance, meaning that they satisfy this condition. We note that this criterion selects stars that are very close but not necessarily exactly on the resonance. Since resonant orbits should be regular, we imposed a second criterion to select these, namely, we also required that λ = 0 and fdr<−1.9. This additional criterion separates these from orbits trapped by resonances temporarily while being intrinsically chaotic.

Fig. 7 shows the position of the stars in today’s (E, L_z, L_⊥) space selected according to their frequencies as described in the previous paragraph, with each bar resonance shown in a different colour. Stars on bar resonances are shown as symbols with black edge-colours, and can be seen to form a more compact set in the space. Clearly, the stripes apparent in Figs. 5 and 6 (which have larger median t_λ and smaller fdr values) correspond to regions hosting resonant bar orbits (which are stable). These results qualitatively agree with the findings by Dillamore et al. (2024), who used a different gravitational potential hosting a bar. We also find that if the bar is non-rotating, such stripes are not present.

The region marking the transition from a specific bar resonance to another orbit is known as separatrix, and orbits located in this region are prone to chaotic behaviour, as is the case for the chaotic orbit shown in Fig. 3. In fact, we find that the vast majority of stars selected according to the orbital frequency criterion to be close to a bar resonance have positive Lyapunov exponents, with roughly 50% having λ > 1 Gyr⁻¹. Such stars essentially move along the inclined lines in (E, L_z) space delineated by the relation E_J = E−Ω_b L_z, which aligns with some of the bar resonances seen in Fig. 7, most notably the corotation resonance.

The above analysis reveals that the highly chaotic regions essentially have two sources. The most important cause for chaotic behaviour is the presence of a triaxial bar in the centre of the Galaxy, which importantly affects stars that pass close to it (those on radial and/or bound orbits), even if this bar is static (see Section 4.1). The second source of chaoticity are the separatrices of the bar resonances as described above, which affect to a larger degree the stars part of the hot thick-disk, as they are on prograde orbits of low inclination.

4.4 Impact on substructures

In Fig. 8, we show the stars associated to the substructures identified by Dodd et al. (2023) colour-coded by their t_λ. If a star’s t_λ is larger than a Hubble time, it is shown in black. The top panel of this figure shows the present-day distribution in (E, L_z, L_⊥) space, while the middle panel shows the distribution after 10 Gyr of integration time in the fiducial rotating bar potential. There, we only show the stars whose final three apocentres were on average larger than 7 kpc, meaning they should still be able to reach or cross a volume of ∼ 1 kpc radius centred on the Solar circle. This requirement removes a bit more than 4% of the stars, which are on radial and bound orbits close to the corotation resonance.

At first glance, the main difference between the top and middle panels is that the distributions of substructures with t_λ smaller than a Hubble time are smeared out. More specifically, the hot thick disk, ED-1 and L-RL3 are now connected in (E, L_z) space. However, we note that the nature of the latter two substructures has never been clear; for instance, their chemistry reveals a mix of populations (Ruiz-Lara et al. 2022; Dodd et al. 2023), so perhaps this is not so unexpected in hindsight. The most bound stars of Gaia-Enceladus form a tail towards lower energy and more negative L_z after 10 Gyr, meaning that they moved along the E_J = E−Ω_b L_z lines described in the previous section, and get closer to the region occupied by Thamnos. The stars associated to the latter substructure cover a range of regimes in t_λ, from short (∼ 3 Gyr) to greater than a Hubble time, which is the case for orbits with L_z<−1000 kpc km s⁻¹.

Substructures that are on less bound orbits (high-energy halo), or have L_z<−1000 kpc km s⁻¹ (retrograde halo), or with significant L_⊥(≳ 1500 kpc km s⁻¹) and hence reaching high above the mid-plane, are generally on regular orbits. These stable regions include structures such as the Helmi Streams, Sequoia, RLR-64/Antheus, ED-2/3/4/5/6, and Typhon.

Given the results presented here, it is clear that a rotating bar affects the distribution, coherence, and purity of substructures in significant parts of (E, L_z, L_⊥) space. Therefore for any clustering analysis, it would be preferable to consider quantities that remain conserved with time. Hence, instead of using the total energy, it might be interesting to use the Jacobi energy E_J (Eq. (1))⁸. The bottom panel of Fig. 8 shows a projection of the substructures in (L_z, E_J),(L_⊥, E_J) after 10 Gyr of integration time in the fiducial potential with the rotating bar. We see that substructures that are not strongly affected by the bar separate from those with short Lyapunov timescales, especially in (L_⊥, E_J) space. There, almost all stars part of substructures with high L_⊥ and/or high E_J have Lyapunov timescales longer than a Hubble time. We note that because of the somewhat squashed appearance of the (E_J, L_z) subspace, it might be advantageous to do some scaling before applying a clustering algorithm in (E_J, L_z, L_⊥) space.

Fig. 8 Top panel: stars part of the substructures identified by Dodd et al. (2023) in (E, Lz, L⊥) space at the present day, colour-coded by tλ. Stars with tλ > 13.7 Gyr are shown in black. The light grey stars in the background correspond to the whole local halo sample. A number of substructures and regions of (E, Lz, L⊥) space are labelled in the leftmost panel. Middle panel: same as the top panel, but after 10 Gyr of integration time in the fiducial potential, and only showing stars whose final three apocentres were on average larger than 7 kpc. Bottom panel: also after 10 Gyr of integration time, but showing (EJ, Lz) and (EJ, L⊥) instead of E. We have indicated the line EJ = E−Ωb Lz for E = −25 000 km2 s−2 (corresponding to the upper limit in energy of the top and middle panel) and E = 0 km2 s−2 (corresponding to unbound stars in our potential).

5 Discussion

5.1 Robustness under different (bar) potential models

The results derived here hold for the fiducial potential, with it its specific bar model, but many bar parameters such as its length, mass, pattern speed and angle are still somewhat uncertain (see Hunt & Vasiliev (2025) for an overview). To quantify how much changing these parameters affects our results, we have explored two different bar models (see Appendix A), and also considered the effect of different bar pattern speeds on the degree of chaoticity as quantified by the Lyapunov exponents (see Appendix B). Since the stars in our sample phase-mix quickly, the bar orientation should not have an effect (we verified that this is indeed the case).

We find that a less massive boxy bar of a similar extent as the fiducial bar potential does not yield significantly different results. A shorter bar, extending only to a few kpc, causes little chaoticity for stars on prograde orbits, and only affects radial orbits with |L_z| ≲ 500 kpc km s⁻¹, but causes a larger and higher degree of chaoticity (see Fig. A.1).

Regarding the dependence on the pattern speed, we find that our results are robust under reasonable variations, with at most 2% differences in the number of stars having λ = 0 (with the number of regular orbits being slightly larger for lower pattern speeds), which is consistent with literature (Manos & Athanassoula 2011). A larger bar pattern speed causes the distribution’s peak to shift towards higher Lyapunov exponents and thus shorter Lyapunov timescales (see Fig. B.1). This is to be expected, as the bar’s time-dependent perturbation literally increases its frequency.

Finally, throughout this work we assumed a spherical dark matter halo, but there is an increasing body of evidence suggesting that it is triaxial (Law & Majewski 2010; Vasiliev et al. 2021; Woudenberg & Helmi 2024). We explored the effect of a triaxial halo by taking the fiducial bar potential and replacing the spherical DM halo by the best-fit halo found by Woudenberg & Helmi (2024). We find no significant changes to our results: the degree of chaoticity is only 2% larger than in the fiducial bar potential, and the resonance structure is roughly similar.

Fig. 9 Percentage of local halo stars close to different bar resonances for the red halo (dashed outlines) and red halo (solid outlines). The percentages on top of each bar corresponds to all selected stars, the white hatched bar indicates the percentage of stars with a Lyapunov exponent <1 Gyr−1 (meaning the non-hatched part of the bar corresponds to stars with Lyapunov exponent >1 Gyr−1), and the doubly hatched bar indicates the percentage of stars with a Lyapunov exponent equal to zero.

5.2 The blue and red halo

It has recently been established (Gaia Collaboration 2018a) that the colour-(absolute) magnitude diagram of halo stars near the Sun shows two sequences: a red sequence, corresponding to stars on “hot” thick disk-like orbits, and a blue sequence, corresponding to accreted stars (Helmi et al. 2018; Gallart et al. 2019). The stars in the two sequences follow different phase-space distributions, and it is interesting to see whether they are affected differently by the bar because of this. Stars in either sequence can be selected by employing a single isochrone that splits the two sequences. In this work, we followed Dodd et al. (2025) and used a BaSTI-IAC alpha-enhanced isochrone with an age of 11.6 Gyr and [M/H]=−0.509, which corresponds to [Fe/H] ∼−0.8 dex. In this way, we find 14 206 stars in the blue sequence and 13 679 in the red sequence.

Using the definitions introduced earlier to identify chaotic orbits, namely λ > 0 and fdr>−1.9, we find that in the fiducial rotating bar potential, 57% of the blue halo stars are on chaotic orbits, versus 63% for those on the red halo sequence. In the static bar potential, these percentages are 53% versus 59% for the blue and red halo sequences, respectively.

Fig. 9 quantifies the percentage of stars in the red and blue halo sequences close to different bar resonances. Although the total percentage for both is similar, adding up to roughly 13%, individual resonances are populated differently due to the difference in the red and blue halo stars’ phase-space distributions. While the red halo, for example, populates the l/m = 0 and l/m = 1/4 resonance more by factors of 2.2 and 2.8, respectively, the blue halo populates the l/m = 1/2 and l/m = 3/2 resonance more by a factor 1.4 and 3.2, respectively. The main reason for the difference in degree of chaoticity, however, seems to be due to the red halo having a greater percentage of stars on highly bound (E<−150 000 km² s⁻²) orbits, where largest number and most chaotic orbits are found, and also where the corotation resonance is located for our potential.

5.3 Limitations and implications

Throughout this work, we have made some simplifying assumptions. Most importantly, we assumed that the bar has fixed parameters. Recent work suggest the bar has slowed down in the past Gyrs (Chiba et al. 2021), which would result in changes in the resonance structure with time. However, while Machado & Manos (2016) find that the degree of chaoticity for stars on bar or disk-like orbits decreases with time in a decelerating bar potential, for dark matter halo particles with average radii >5 kpc this decrease is minor. Indeed, the rate of change in bar parameters inferred is typically larger than the orbital timescales for the stars in our sample and could thus be considered to be an adiabatic perturbation. More importantly, we show in Sect. 4.1 that most of the chaoticity is driven by the triaxial shape of the bar, rather than by its rotation, so most of the results presented here should hold.

We also neglected the perturbations from the infall of the LMC (Garavito-Camargo et al. 2019; Vasiliev 2024), but this will likely mostly affect stars with large apocentres, on which the bar has no significant effect unless they are extremely radial. Lastly, we neglected the perturbations from the spiral arms, but as these might be transient and limited to very low z, they are dynamically not of great importance for our sample of halo stars.

In this work, we considered a relatively small local volume (d<1 kpc) to ensure completeness. This also implies that we obtained a heliocentric perspective, as the effect of the bar will differ depending on location in the Galaxy. We observe that a small fraction of stars, mainly those on highly bound and radial orbits, can shrink their apocentre to values below 7 kpc within a few Gyr due to the bar. Hence, also in this sense, the observed distribution in (E, L_z, L_⊥) space should be taken as a snapshot taken at the present-time.

The high degrees of chaoticity found in this work imply that some structures, depending on their present-day location in (E,.L_z, L_⊥), may have phase-mixed faster than previously thought (Helmi & White 1999). Spatially coherent stellar streams may also have been affected by the bar. In fact, chaoticity has been put forward as a hypothesis to explain the short length of the Ophiuchus stream Price-Whelan et al. (2016b); Yang et al. (2025), and it would be interesting to explore the presence of such signatures for other streams whose orbit might be chaotic. Next, we find that stars on chaotic orbits move along stripes with positive slope in (E, L_z) space defined by the relation E_J = E−Ω_b L_z, particularly when close to a bar resonance. Hence, a clump of stars in (E, L_z, L_⊥) space, if sufficiently large, may spread out along such a stripe over time. This is consistent with the recent results by Dillamore & Sanders (2025).

We also observe the presence of a large number of bar resonances in the stellar halo. Obviously, the exact locations and number of stars close to bar resonances in (E, L_z, L_⊥) space depends strongly on the exact configuration of the potential. Although the percentage of stars on bar resonances is relatively small, of the order of a few percent, these can cause overdensities, whose interpretation should be conducted with care. Only if kinematic features that appear to be introduced by the bar are detected, can these be used to constrain a bar model, as recently attempted by Dillamore et al. (2025).

6 Conclusion

In this work, we explored the effect of a rotating stellar bar on the orbits of halo stars located within 1 kpc from the Sun, and more specifically the effect on their distribution in (E, L_z, L_⊥) space. We characterised their orbits using orbital frequencies and two chaoticity indicators, the Lyapunov exponent λ and the frequency diffusion rate fdr. We selected stars close to bar resonances and mapped their location in (E, L_z, L_⊥) space.

Using as criteria for the identification of truly chaotic orbits that λ>0 and fdr>−1.9, we find that, while in an axisymmetric potential only ∼ 10% of the orbits is strongly chaotic, in a static bar potential this increases to 56%, and in a rotating bar potential this increases to 60%. This fraction is not strongly dependent on the bar pattern speed, bar orientation, or other of its properties. Stars on the red halo sequence are slightly more often on chaotic orbits than those in the blue halo sequence for a range of rotating bar potentials. This is driven by their different phasespace distribution. We find that a large percentage of stars on bound, radial and prograde orbits are chaotic, interspersed with regions of more stable orbits corresponding to the locations of bar resonances. We find that the main driver of chaoticity is the bar’s triaxial shape, and that the importance of bar resonances (and their overlap) as a driver of chaotic evolution appears to be secondary.

We quantified the impact of the chaoticity over a Hubble time by computing the Lyapunov timescale t_λ (including a chaos onset time). We find that most stars with high binding energy and |L_z| ≲ 1000 kpc km s⁻¹ are on chaotic orbits that have Lyapunov timescales that are much shorter than a Hubble time. Stars may move over time in L_⊥ and along a line in (E, L_z) space defined by the relation E_J = E−Ω_b L_z when close to bar resonances. Orbits can change from prograde to retrograde, from disk-like to reaching high above the mid-plane, and their apocentres can grow or shrink. On the other hand, we find that stars with lower binding energy (≳−80 000 km² s⁻²), or which reach high above the Galactic plane (L_⊥ ≳ 1500 kpc km s⁻¹), or which are highly retrograde (L_z ≲−1000 kpc km s⁻¹), have Lyapunov timescales longer than a Hubble time and do not change their (E, L_z, L_⊥) significantly.

As a consequence of the above, the substructures identified as potentially being associated to merger debris and located in these specific regions of (E, L_z, L_⊥) become more spread out. We explored this for the substructures identified by Dodd et al. (2023), and we find important changes that may lead to overlap for substructures that are on low inclination prograde orbits, as well as for those that are strongly bound and on radial orbits. This leads us to propose the use of (E_J, L_z, L_⊥) space in future clustering analysis, as E_J is a conserved quantity in a rotating bar potential (and will change only adiabatically in case the bar slows down) and substructures separate well in this space. We find that substructures on less bound orbits remain fairly coherent even in the presence of the bar.

These results likely also have implications for Galactic globular clusters and their possible associations to different merger events (see for example Massari et al. 2019). Interestingly, in recent work, Bajkova et al. (2025) find comparable fractions of globular clusters on chaotic and on regular orbits as found for our local sample of halo stars, even though they probe rather different regions of the Galaxy.

Acknowledgements

This research has been supported by a Spinoza Award from NWO (SPI 78-411). We thank the referee for constructive comments on the paper. HCW thanks Eugene Vasiliev for sharing his expertise and his help with a range of questions on AGAMA, and Emma Dodd for helpful discussions. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Throughout this work, we have made use of the following packages: astropy (Astropy Collaboration 2013), vaex (Breddels & Veljanoski 2018), SciPy (Virtanen et al. 2020), matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), AGAMA (Vasiliev 2019), SuperFreq (Price-Whelan 2015), tqdm (da Costa-Luis et al. 2024), and Jupyter Notebooks (Kluyver et al. 2016).

Appendix A Different bar models

The length, mass, and mass distribution of the bar are uncertain. The Sormani et al. (2022) bar model is relatively heavy and long, reaching a radius of about 5 kpc. However, recent work suggests that the length of the bar might be overestimated, as stars part of an attached spiral arm are included, which would bring the length of the bar down to about 3.5 kpc (Lucey et al. 2023; Vislosky et al. 2024). As the bar temporarily overlaps with the inner spiral arms, it can appear up to twice its true size (Hilmi et al. 2020). To investigate how having a less massive and shorter bar with a different mass distribution would change our results, we took the fiducial potential from Woudenberg & Helmi (2024) with a spherical NFW halo (for consistency with the fiducial potential) and replaced its bulge by a bar of the same mass (8.9 · 10⁹ M_☉). This bar is a boxy bar of the form $$r={({\left(\frac{x}{a}\right)}^k+{\left(\frac{y}{b}\right)}^k+{\left(\frac{z}{c}\right)}^k)}^{1/k}$$(A.1) $$\rho \left(r\right)={\rho }_0{e}^{-r^{1/n}}$$(A.2) with k = 4. We kept the pattern speed and bar angle the same as in the fiducial potential, and investigated two different bar parametrizations. Firstly, we investigated a boxy bar that has a similar length as the Sormani et al. (2022) bar, with parameters $n=1.1,a=\frac{7}{10},b=\frac{7}{40},c=\frac{7}{80}$ and ρ₀ = 1.4 · 10¹⁰ M_☉ kpc⁻³. Secondly, we investigated a short boxy bar with a much smaller extent than the Sormani et al. (2022) bar, about 2.5 kpc, with parameters $n=0.7,a=\frac{7}{10},b=c=\frac{7}{40}$, and ρ₀ = 2.9 · 10¹⁰ M_☉ kpc⁻³. We checked that the rotation curves of these two models agree with the rotation curve data by Eilers et al. (2019) and Ou et al. (2023).

The distribution of Lyapunov timescales in IoM space is shown in Fig. A.1. The distribution and amount of stars close to and on bar resonances for the boxy bar resemble the results we find for the fiducial potential, though of course their exact location in (E, L_z, L_⊥) space is shifted (for example, the corotation resonance shifted towards L_z ∼ 0 kpc km s⁻¹ and more bound orbits). The percentage of stars with Lyapunov exponent >0 and fdr > −1.9 is 61% (with 63% in the red halo and 60% in the blue halo). 66% of stars with |L_z|<500 kpc km s⁻¹ are chaotic. Instead, for the short boxy bar, the percentage of stars close to bar resonances is halved, and only stars in the range −500 ≲ L_z ≲ ∼ 500 kpc km s⁻¹ are for a large percentage on chaotic orbits (which are generally more chaotic, meaning having higher Lyapunov exponents and fdr, than for the boxy bar or fiducial potential). The most bound stars in the short boxy bar model are for a larger percentage chaotic. We found that the fdr determined using the frequencies in Poincaré symplectic coordinates in the rotating frame was more reliable than in the Cartesian rotating frame. Using that, the percentage of stars with Lyapunov exponent >0 and fdr > −1.9 is 60% (with 63% in the red halo and 57% in the blue halo). 81% of stars with |L_z|<500 kpc km s⁻¹ are chaotic. Although the location of chaotic orbits is slightly different, for both models, the percentage of stars on chaotic orbits agrees with what has been found for the fiducial potential (see Sect. 4.1).

Fig. A.1 Local halo stars in (E, Lz, L⊥) space at the present day, binned into 100 bins in both x and y direction, showing the median Lyapunov time in the boxy bar (top row) and short boxy bar model (bottom row).

Appendix B Effect of different pattern speeds and bar orientation on the Lyapunov exponent

To see what influence a different pattern speed has on the presence of chaoticity, we took the fiducial potential and integrated orbits of the 1 kpc volume halo sample as done before. We computed the Lyapunov exponent for various bar pattern speeds, Ω_b = [−32.5, −35,−37.5,−40,−42.5] kpc km s⁻¹ at a fixed bar angle, −25 degrees. The bar angle does not influence the Lyapunov distributions (we ensured this is indeed the case), as the sample we’re considering is phase-mixed. We also compared this to the distribution of Lyapunov exponents when the bar is static (Ω_B = 0). The results are shown in Fig. B.1, which shows for a higher bar pattern speed higher values of the Lyapunov exponent are reached. This is consistent with the findings of Price-Whelan et al. (2016b), see their Fig. 5. When the bar is static, there is a larger percentage of stars on regular orbits (meaning with λ = 0).

Fig. B.1 Lyapunov exponent distributions for the local halo sample in the fiducial potential with various bar pattern speeds, Ωb = [−32.5,−35, −37.5,−40,−42.5] km s−1 kpc−1, at a fixed bar angle, −25 degrees. The grey lines show all cases, while the blue lines show one pattern speed as indicated in the legend. For comparison, the red dotted line shows the distribution for a static bar (Ωb = 0).

References

All Figures

	Fig. 1 Distribution of nearby halo stars in E−Lz (top row) and L⊥−Lz (bottom row) at t = 0 and after 10 Gyr of evolution in our fiducial Milky Way potential including a rotating bar. Especially highly bound stars with Lz ∼ 0 (and small L⊥) change their distribution significantly.
In the text

Fig. 2 Examples of orbits of stars on bar resonance, identified as having $\frac{{\Omega }_{\varphi ,\mathrm{CR}}}{{\Omega }_{R,\mathrm{CR}}}=l/m\pm 0.005$ and λ = 0 and fdr<−1.9. The orbits are shown in the corotating frame, and the resonance they are associated to is indicated in the top right corner. Isodensity contours of the fiducial potential in the mid-plane (z = 0) are shown in the background, the black cross indicates the Galactic Centre.

In the text

Fig. 3 Example of a star on a chaotic orbit (selected from regime1). The red and blue markers show the position of the star at the present day and after 10 Gyr of integration time, respectively. The top panels show projections of the orbit in the corotating frame, (xCR, yCR), and in (R, z), respectively. The bottom panels illustrate the trajectory of this star in (E, Lz, L⊥) space over time. This star ends up librating around the bar’s corotation resonance after 10 Gyr, and moves around significantly in (E, Lz, L⊥) space as it changes its orbit. The orbit has λ = 2.1 Gyr−1 and tλ = 1.3 Gyr, and fdr=−0.2 (see Sect. 3.1 and 3.2 for details on how these quantities are computed).

In the text

Fig. 4 Distribution of the values of the two chaoticity indicators used in this work for our local halo sample: the Lyapunov exponent λ, normalised by Torb,r, and the fdr. The different panels correspond to three different versions of the potential described in Section 2.1: with a rotating bar (left), with a static bar (middle), and an axisymmetrised bar (right). The side panels show the cumulative distributions of normalised λ for each of the potentials, while the top panels correspond to those of fdr. The fdr distribution for the orbits that have λ = 0 is shown separately. The vertical magenta lines correspond to fdr=−1.9, the 99 th percentile of the orbits with λ = 0 in the fiducial potential.

In the text

	Fig. 5 Local halo stars in (E, Lz, L⊥) space computed at t = 0 in the fiducial potential, binned into 100 bins along each axis, where the colour coding of each bin shows the percentage of stars with a Lyapunov exponent higher than 0 and a fdr >−1.9, determined after 25 Gyr of integration.
In the text

	Fig. 6 Similar to Fig. 5 but showing in the top panel the weighted median Lyapunov timescale (Eq. (9)), and in the bottom panel the median fdr (Eq. (7)) determined over an integration time of 25 Gyr for all stars in the sample. We note the good correspondence between both panels, and with Fig. 5.
In the text

	Fig. 7 Position of local halo stars close to bar resonances in (E, Lz, L⊥) space at the present day, selected to have $\frac{{\Omega }_{\varphi ,\mathrm{CR}}}{{\Omega }_{R,\mathrm{CR}}}=l/m\pm 0.005$ (see Fig. 2 for examples of each resonant orbit). The stars encircled in black are on resonant bar orbits, meaning with λ = 0 and fdr <−1.9.
In the text

Fig. 8 Top panel: stars part of the substructures identified by Dodd et al. (2023) in (E, Lz, L⊥) space at the present day, colour-coded by tλ. Stars with tλ > 13.7 Gyr are shown in black. The light grey stars in the background correspond to the whole local halo sample. A number of substructures and regions of (E, Lz, L⊥) space are labelled in the leftmost panel. Middle panel: same as the top panel, but after 10 Gyr of integration time in the fiducial potential, and only showing stars whose final three apocentres were on average larger than 7 kpc. Bottom panel: also after 10 Gyr of integration time, but showing (EJ, Lz) and (EJ, L⊥) instead of E. We have indicated the line EJ = E−Ωb Lz for E = −25 000 km2 s−2 (corresponding to the upper limit in energy of the top and middle panel) and E = 0 km2 s−2 (corresponding to unbound stars in our potential).

In the text

Fig. 9 Percentage of local halo stars close to different bar resonances for the red halo (dashed outlines) and red halo (solid outlines). The percentages on top of each bar corresponds to all selected stars, the white hatched bar indicates the percentage of stars with a Lyapunov exponent <1 Gyr−1 (meaning the non-hatched part of the bar corresponds to stars with Lyapunov exponent >1 Gyr−1), and the doubly hatched bar indicates the percentage of stars with a Lyapunov exponent equal to zero.

In the text

	Fig. A.1 Local halo stars in (E, Lz, L⊥) space at the present day, binned into 100 bins in both x and y direction, showing the median Lyapunov time in the boxy bar (top row) and short boxy bar model (bottom row).
In the text

	Fig. B.1 Lyapunov exponent distributions for the local halo sample in the fiducial potential with various bar pattern speeds, Ωb = [−32.5,−35, −37.5,−40,−42.5] km s−1 kpc−1, at a fixed bar angle, −25 degrees. The grey lines show all cases, while the blue lines show one pattern speed as indicated in the legend. For comparison, the red dotted line shows the distribution for a static bar (Ωb = 0).
In the text