Globular clusters and bar: captured or not captured?

Studies of the dynamics of globular clusters assume different values of bar parameters (mass, velocity, size) and analyse the results of orbit classifications over the range of the chosen values. It is also a usual thing that a spherical bulge component is converted into the bar to obtain a non-axisymmetric potential from an axisymmetric one. The choice of bar parameters and the way the bar is converted from the bulge introduce systematics into the orbit classifications that we explore in the present study. We integrate orbits of 30 bulge globular clusters residing in the inner area of the Galaxy ( 𝑅 ≲ 5 kpc) backwards in time for three different potentials, two of which are obtained by fitting the rotation curve, and one is taken from the surrogate 𝑁 -body model representing our Galaxy. We analyse each orbit in terms of dominant frequencies obtained from its coordinate spectra. We find that the bar pattern speed is a key factor in orbital classification. With an increase of it, frequencies deviate more and more from the “bar” frequency ratio 2:1. The bar-to-bulge mass ratio (assuming the total mass of the bar plus the bulge is fixed) and size of the bar play a smaller role. We also find that, in the 𝑁 -body potential, the fraction of orbits that follow the bar is higher than in those obtained from fitting the rotation curve.


INTRODUCTION
Several physical components co-exist within the area of about 5 kiloparsecs from the centre of our Galaxy.These components are a bar, its vertically thick part, which is usually referred to as the boxy/peanut-shaped (B/P) bulge (McWilliam & Zoccali 2010;Nataf et al. 2010;Wegg & Gerhard 2013;Mosenkov et al. 2021), and possibly an another bulge, commonly referred to as the classical one.The existence of the latter has came under the question in the past few years due to various indicators pointing out that bulge stars exhibit cylindrical rotation (Kunder et al. 2012;Ness et al. 2013Ness et al. , 2016)), i. e. support the B/P bulge rather than the classical one, although there are some exceptions (Kunder et al. 2016, also see the review by Bland-Hawthorn & Gerhard 2016).We do not mention here the most inner part subsystems, such as the nuclear disc and the nuclear star cluster (Becklin & Neugebauer 1968), since they are not relevant to the present work and are important for considering on much smaller spatial scales than those considered here.
Globular clusters (GCs) are tracers of the secular evolution of bar and bulge components, since GCs include a large bulk of stars that reflect how these components form/evolve in their metallicity and stellar populations.However, the question of whether a particular GC belongs to a certain component (e.g. a bulge, a bar, a disc, or a halo) is not easy to answer.On the contrary, determining the origin of a globular cluster is a rather difficult task, which requires reliable knowledge of the clusters' proper motions, their radial velocities, positions, and metallicity (Côté 1999;Bica et al. 2016;Massari et al. 2019;Pasquato & Chung 2019;Ortolani et al. 2019a,b;Pérez-Villegas et al. 2020;Bajkova et al. 2020a;Bajkova & Bobylev 2021;Sun et al. 2023).For example, Ortolani et al. (2019a) recently found that the CGs Terzan 10 and Djorgovski 1 have typical halo orbits, while their orbits are contained within the bulge volume.Another illustrative example is that Pérez-Villegas et al. (2020) and Ortolani et al. (2019b) showed that several GCs, while do not belong to either the disc or the halo and appear to belong to the bulge, nevertheless do not follow the bar.Meaning that these GCs move either faster or slower than the bar, but not synchronously with it.The ambiguity in the classification of GCs stems from the fact that several physical components of the Galaxy mentioned above overlap in physical space and, at the same time, the observations of the inner part of the Galaxy are affected by heavy extinction and crowding (Bland-Hawthorn & Gerhard 2016).
An additional problem, which especially concerns the dynamics of the GCs of the inner Galaxy and the classifications based on it, is that the parameters of the bar itself are also not set in stone.Bar pattern speed estimates range from about 30 km/s/kpc to 40 km/s/kpc (Portail et al. 2017;Bovy et al. 2019;Sanders et al. 2019;Binney 2020;Asano et al. 2020;Kawata et al. 2021;Chiba & Schönrich 2021;Li et al. 2022;Clarke & Gerhard 2022), while some authors provide an even higher value of about 50 km/s/kpc (Minchev et al. 2007;Antoja et al. 2014).Naturally, the centrifugal force that influence the motion of GCs depends on the bar pattern speed.It is also important that the changes in the pattern speed force the resonances to move, and, thus, orbits will differ depending on how close the GC to a particular resonance.Therefore, the classifications of the orbits should differ depending on the bar pattern speed, and one should consider a set of bar pattern speed values, as it was done, for example, in Ortolani et al. (2019b) and Pérez-Villegas et al. (2020).In Pérez-Villegas et al. (2020), the authors calculated the probability that an orbit belongs to one or another component separately for each of the pattern speeds considered there.
The uncertainty in the existence of classical bulge mentioned above can also implicitly affect the results of GCs' classification.One of the approaches to modelling of GCs' orbits is to transform the spherical central component into a bar.This means the the central spherical bulge in the originally axisymmetric model of the Milky Way is replaced by an elongated bar with exactly the same mass as the bulge.This approach has been used in recent studies by Ortolani et al. (2019a,b); Pérez-Villegas et al. (2020) and many previous ones.At the same time, various -body studies showed that the inclusion of even a small classic bulge component can drastically change the overall evolution of the model, leading to the formation of the socalled barlenses (Salo & Laurikainen 2017;Smirnov & Sotnikova 2018) or preventing the bar buckling (Smirnov & Sotnikova 2018) altogether.
In the present work, we want to address the mentioned issues in the context of the capturing of CGs by the bar.We want to explore how the choice of the bar parameters (pattern speed, mass, size) affects the state of the CGs relative to the bar, i. e. is there any systematics in the frequency rations of GCs orbits  R /   (see definition is Section 3) depending on the bar parameters.
To this aim, we study the motion of GCs in three different instances of the Milky Way potential.Two of them are based on observational data from Bajkova & Bobylev (2016, 2017) and McMillan (2017) and one is based on the -body model from Tepper-Garcia et al. (2021), which was specifically prepared to represent the mass distribution of the Milky Way and has a spatial resolution of about 30 pc.This -body model also contains a classical bulge and a naturally formed bar, thus providing an opportunity to study the GCs kinematics in case of a self-consistent model, obtained without transforming one component into another.
The article is structured as follows.In Section 2, we describe our sample of GCs.In Section 3, we provide details on the potentials considered in the present work and how the classification and integration of the orbits backwards in time was carried out.In Section 4, we analyse the systematics in the classification of orbits introduced by changing bar parameters using one GC, NGC 6266, as an example.Section 5 presents the results of the classification for all GCs in the sample.We compare our results with those of previous works in Section 6.In Section 7, we give our conclusions.

DATA
To study the kinematics of GCs in different barred potentials, we first selected 30 CGs, which were previously identified in Bajkova et al. (2020a) as those that belong to the bar/bulge.These GCs were selected from a catalogue of 152 GCs from Bajkova & Bobylev (2021) based on the following criteria.First, a geometric criterion was applied to retain only those GCs whose apocentric distance  apo is less than 3.5 kpc (Massari et al. 2019;Bajkova et al. 2020a).This reduces the sample to 39 members.Then, nine GCs were found to belong to the disc based on the angular momentum and eccentricity of the corresponding orbits (see details in Bajkova et al. 2020a) and, thus, were removed from the sample.Table 1 and Table 2 list the chosen GCs, as well as their observational parameters and Cartesian coordinates and velocities, used below to integrate orbits backwards in time for 5 Gyr.Coordinates and velocities are obtained from equatorial coordinates ( 2000 ,  2000 ), line-of-sight velocities from the catalogue of Vasiliev (2019b), distances from Baumgardt & Vasiliev (2021), and proper motions from Vasiliev & Baumgardt (2021).The catalogue of Vasiliev (2019b) is compiled based on the Gaia DR2 data, while the catalogues of Vasiliev & Baumgardt (2021); Baumgardt & Vasiliev (2021) contain new proper motions and refined distances based on Gaia EDR3 data, Hubble Space Telescope (HST) data, and some literature estimates.The transformation from angular coordinates and velocities is performed using the values obtained by Bajkova & Bobylev (2016, 2017) from rotation curve fitting, i.e. under the assumption that the distance from the Galaxy centre to the Sun  ⊙ = 8.3 kpc, the height of the Sun above the disc plane  ⊙ = 17 pc (Bobylev & Bajkova 2016), the velocity of local standard of rest (LSR)  ⊙ = 244 km/s.The peculiar velocity of the Sun relative to LSR ( ⊙ ,  ⊙ ,  ⊙ ) = (−11.1,12.2, 7.3) km/s is taken from Schönrich et al. (2010).For the bar viewing angle, the value 23 deg was taken from Mosenkov et al. (2021), where it was estimated from fitting the boxy/peanut bulge intensity profile for different viewing angles.

Mass models
In the present work, we consider several types of mass models of the Milky Way.The first one was obtained by Bajkova & Bobylev (2016, 2017) (hereinafter, BB2016) via fitting the rotation curve to the kinematic data of a set of different objects with distances up 200 taken from Bhattacharjee et al. (2014).The mass model consists of three distinct components, namely the bulge (Plummer 1911), the disc (Miyamoto & Nagai 1975), and the halo (Navarro et al. 1996): (3) Description of the parameters and their respective values are given in Table 3.The second model is taken from McMillan (2017) (hereinafter, MC2017) and consists of six different components, namely thin and thick stellar discs, dark matter halo, and H I and molecular discs.In this model, the dark halo is also described by a Navarro-Frank-White profile, given in eq. ( 3).The stellar discs are exponential both in the plane and the vertical direction: while gaseous discs are exponential in the plane and isothermal in the vertical direction and have a hole in the centre with the scale of  m : The central component (bulge) is implemented via the following parametric model: where  ′ = √︃  2 + (/ bulge ) 2 .To avoid repeatance, we refer the reader to McMillan (2017) for a description of the parameters and their values.In both models, we introduce a bar component by decreasing the mass of the central component (bulge) by a certain value and then assigning the bar mass to this value.Essentially, this means that, for all models considered below (except the -body one), the total mass of a spherical bulge and the bar is fixed: where  b,0 is the initial bulge mass of the axisymmetric model and  b is the residue mass of the bulge.Hereinafter, we refer to  b,0 simply as  b , since we do not consider the residue bulge mass as an independent parameter at any part of this work.Below, we consider a set of bar mass values, or, more precisely, a number of bar-tobulge mass ratios  bar / b (see Table 3).Ortolani et al. (2019a,b); Pérez-Villegas et al. ( 2020) assigned all the bulge mass to the bar component in their models.Here, we introduce the ratio of the bulge and bar masses as a free parameter to investigate how uncertainty in the classical bulge parameters possibly existing in our Galaxy can affect the results of orbital classification.For the bar density profile, we take a Ferrers profile: where  bar is the bar mass,  bar is the bar major axis, r = √︁  2 + (/) 2 + (/) 2 is the elliptical radius and  and  characterise the flattening of the bar in disc plane and along the vertical axis, respectively.The bar parameters and their description are given in Table 3.
The third type of potential is taken from a recent work by Tepper-Garcia et al. (2021) (hereinafter, TG2021), where a surrogate Milky Way -body model was presented.Time snapshots of the model were made publicly available by the authors.At start of the simulations, the model consisted of two spherical components, NWF-like halo (Navarro et al. 1996) and a stellar bulge (Hernquist 1990), and an exponential disc isothermal in the vertical direction (similar to eq. ( 5), but without the hole).The evolution of the model was followed up to about 4.3 Gyr.There is no need to insert the bar component separately or transform the bulge as the bar in this model is formed naturally (see Fig. 3).For simplicity, we consider here only the last snapshot of Tepper-Garcia et al. ( 2021)'s simulations, neglecting the time evolution of the bar properties.We leave this for future studies.
For the selected time moment, the -body bar has the size  bar about 4.5 kpc and the pattern speed Ω p ≈ 39 km/s/kpc.The mass of the bar was not estimated directly in Tepper-Garcia et al. ( 2021), but the authors provided an overall estimate  bar +  disc +  bulge = 3.5 × 10 10  ⊙ of stellar mass inside the area of  < 5 kpc (bar region), where  disc is the mass associated with the inner area of the disc and  bulge is the mass of the classic bulge originally included in the model.
The number of particles in the -body model is about 7 • 10 7 .To avoid very time-consuming calculations of gravitational force at each time-step when integrating the orbits, we prepared a multipole expansion of the potential using the convenient Multipole subroutine from AGAMA software package (Vasiliev 2019a):  where    are spherical functions of degree  and order .We truncate the series at  max = 6 and  max = 6 and impose a triaxial type of symmetry (only even harmonics are calculated).Isolines of the potential approximations are shown in the right panel of Fig. 3.Note that the potential isolines are rounder than the density isolines, as they should be (Binney & Tremaine 2008, Chapter 2), but still showing the flattening in the bar area.In the very central part, the classic bulge overweighs other components and, thus, the isolines are circular here.
For potentials of BB2016 and MC2017, we consider a range of bar pattern speeds (Ω p ) and sizes ( bar , the half length of the bar major axis), from 30 km/s/kpc to 60 km/s/kpc and from 5.0 kpc to 2.5 kpc, respectively.Fig. 2 shows how the mentioned limits correspond to the main resonances in the potentials of BB2016 and MC2017.The dynamics of the bars is usually characterised by the rotation rate parameter R =  CR / bar (Binney & Tremaine 2008).As can be seen from the figure, we consider both slow (R ≫ 1) and fast bars (R ≲ 1) here.For other galaxies, R spans the range from almost 0 to about 4 (Buta & Zhang 2009;Cuomo et al. 2019;Guo et al. 2019;Garma-Oehmichen et al. 2020, 2022).For our bars, R is from about 0.5 for Ω p = 60 km/s/kpc and  bar = 5 kpc to about 10 for Ω p = 10 km/s/kpc and  bar = 2.5 kpc.For generality, we consider bars with major axes and pattern speeds having the values from the suggested ranges indisciminately, although longer bars tend to have lower pattern speeds (see figure 15 in Garma-Oehmichen et al. 2022).

Orbit integration and classification
For each of the described potentials we add the rotation in accordance with the chosen value of the pattern speed Ω p and integrate orbits of GCs backwards in time for a time period of 5 Gyr.In the present work, we are interested in the orbital families which an orbit can belong to potentially.This "property" should not depend on the type of the integration (forward or backaward) for regular orbits1 , since orbital frequencies are "integral" properties of the orbit.Here, we  consider backward integration, because it is in line with our previous studies, e. g.Bajkova & Bobylev (2019); Bajkova et al. (2020b).Integration is carried out using the AGAMA software package.AGAMA performs integration via 8th order Runge-Kutta scheme with an adaptive time-step.We choose an output time step Δ = 1 Myr.The latter is unrelated to the actual time-step of integration, which is determined internally in ODE solver based on the imposed value of the relative accuracy, 10 −8 in our case.For each orbit, we traced the evolution of the Jacobi energy as an indicator of the accuracy of our calculations.A typical example is shown in Fig. 4. In short, the energy is well conserved during integration (up to six decimal places).
To classify orbits, we apply the methods of spectral dynamics pioneered by Binney & Spergel (1982).In this approach, one calculates the coordinate spectra of the orbit, i.e.Fourier transforms of the time series , , , and  taken in a bar rotating frame, and then finds dominant frequencies   ,   ,   , and  R corresponding to the highest spectral lines.The spectra are calculated as follows where   = /Δ, Δ = 5 Gyr,   = Δ, Δ = 1 Myr, 0 ≤  ≤ (  − 1)/2, and   is the length of the time series.To improve the resolution of the peaks, we use a subroutine similar to zero-padding (see details in Parul et al. 2020, where a similar analysis was applied to the study the orbital families of B/PS bulges).For regular orbits, the spectra consist of discrete lines, these lines can be distinguished, and the corresponding frequencies can be studied to understand which orbital group or family the orbit belongs  to (Binney & Spergel 1982).This approach made it possible to obtain many fruitful results on the orbital composition of the bar and the importance of various resonances for the structure of the bar in a number of studies (Gajda et al. 2016;Wang et al. 2016;Portail et al. 2017;Łokas 2019;Parul et al. 2020;Tikhonenko et al. 2021;Smirnov et al. 2021).Pérez-Villegas et al. ( 2020) also calculated the orbital frequencies to determine whether a particular GC follows the bar or not.Here, we use the same approach and assume that if  R /   = 2.0 ± 0.1, then the GC with such a ratio of frequencies is the bar supporting one, i.e. follows the bar.
A typical example of the orbit of NGC 6266, along with its time series of coordinates and their spectra, is presented in Fig. 5. Hereinafter, all orbits presented in the figures are shown in the bar rotating frame if not specified otherwise.Bar parameters are Ω p = 45 km/s/kpc,  bar / b = 0.95,  bar = 5 kpc,  = 2.0, and  = 3.0 in this case.Integration is carried out in the potential of BB2016.We note that, although the orbit has a nice-looking regular profile, it does not actually follow the bar, since  R /   ≈ 3.5.

ORBIT TYPE DEPENDING ON THE BAR PARAMETERS
First of all, we would like to explore how the choice of bar parameters affects the type of orbit.We begin this Section by considering only  one GC, namely NGC 6626.There is no particular reason for this choice, except that this example is illustrative.By a detailed analysis of one orbit, we outline the systematics in the classification of orbits that arise due to changes in the bar parameters.
Fig. 6 shows how frequencies   and  R and their ratio  R /   change with the bar pattern speed, mass, and size for the potential of BB2016.To study dependencies, we first consider the onedimensional case, where one parameter changes, while the rest are fixed.Unless otherwise specified, all bar parameters are fixed at the following values: Ω p = 45 km/s/kpc  bar / b = 0.95,  bar = 5.0 kpc,  = 2.0, and  = 3.0.We present orbital profiles in Fig. 7 to illustrate how they change when the corresponding parameter is changed.As can be seen from the individual subpanels, there are clear systematic shifts in frequencies and, accordingly, frequency ratios: (i) With an increase of the pattern speed, frequency of radial oscillations decreases.This continues up to the point at about 24 km/s/kpc, then the frequency increases abruptly, after which it remains constant.
(ii) Frequency   decreases monotonically with an increase in the pattern speed.
(iii) The frequency ratio  R /   shows an interesting behaviour as a result of changes in individual frequencies.Initially,  R /   = 2 (a typical ratio for orbits following the bar), but at Ω p ≈ 24 km/s/kpc and after that, it deviates more and more from this value.
The described changes of frequencies are reflected in the orbit profile.In the case of  R /   ≈ 2, one observes a very regular orbit captured by the bar.For  R /   ≳ 2, the orbit becomes more "windy" and now oscillates around the bar.
For the bar-to-bulge mass ratio and size (second and third rows of Fig. 6), one can see that changing these parameters affects the orbit profile and the corresponding frequency ratios, but their influence is not so strong compared to the pattern speed.An increase in the bar mass and size leads to a slight decrease in   , which leads to small changes in the frequency ratios, from  R /   ≈ 2.2 − 2.4 at the left boundary of the interval to about  R /   ≈ 3.0 on the right.
Comparing the results for BB2016 (Fig. 6) and MC2017 (Fig. 8), one can see that the trends in changes of frequencies between them are similar, i.e. there is a sudden change in the frequency ratio at a particular value of the bar pattern speed.For the MC2017 potential, this change occurs at a somewhat smaller value of Ω ≈ 20 km/s/kpc.In the case of MC2017, changing the bar-to-bulge mass ratio has almost no effect on the frequency ratio.This can be explained by the fact that the bulge in the model of MC2017 already has a certain degree of flatness (along the vertical direction) and its transformation into an elongated component does not significantly affect the potential.
In Fig. 6 and Fig. 8, we fixed all bar parameters, except for one, which then varied.However, doing so, we did not take into account the possibility that, with a different combination of bar parameters, the observed dependencies may well change or simply disappear.To explore such a behaviour in more detail, we conduct a following suit of simulations.We run Monte-Carlo simulations, choosing a set of bar parameters from the intervals specified in Table 3 uniformly, then we calculate the orbit and the corresponding ratio of its frequencies.We performed 10 5 of such iterations.Fig. 9 show the results in a form of a matrix plot for all parameters, with the average value of frequency ratio for a given pixel highlighted in different colours.Each subplot presents a 2D histogram obtained by averaging the values within 100 bins from minimum to maximum values for each axis.The subplots show qualitatively similar results compared to those presented in the 1D plots (Fig. 6 and Fig. 8).Again, the pattern speed is the most important parameter, i. e. in each subpanel in the first column there is a gradual progression of colours.For other parameters, there is no such correlations, except for a weak correlation of frequency rations with  bar .Thus, changing all other parameters does not strongly affect the frequency ratio.This means that the pattern speed may be very well the most important factor when one is trying to asses orbit families and check whether a particular orbit follows a bar or not.
To understand why frequencies abruptly change with the pattern speed, we calculated the Poincare surface of sections (SoSs) for the range of pattern speeds.For 3D orbits, SoSs are four dimensional objects, i.e. (, ,   ,   ) taken at  = 0 and   > 0 (or any other similar combinations).Here, we plot SoS projections on (,   ) plane taken at  = 0 and   > 0. A similar approach was used in Kalapotharakos et al. (2004); Voglis et al. (2007), where 3D -body orbits were studied.Fig. 10 demonstrate how the SoSs change either with Ω p or with the corresponding frequency ratio  R /   .Note that the SoSs presented are not exactly typical.Usually, the Jacobi energy is a fixed variable and one investigates various orbits for a chosen energy value.In Fig. 10, the pattern speed (and, thus, the corresponding Jacobi energy) changes from orbit to orbit, not the initial velocity or position.Nevertheless, as can be seen from the figures, the family to which the orbit belongs gradually changes with the pattern speed.The orbit starts on the island close to 1 family (they reside in rightmost corner of the plot), then gradually expands to the left side of the diagram.At some point (after Ω p ≳ 30), new islands appear.The orbit clearly ceases to be a member of the 1 family, as indicated by its increasing frequency ratio  R /   .As for the question of which family an orbit ends up, this question is not easy to answer, since a bar can be populated by orbits with multiplicity greater than 2:1, see a recent work by Wang et al. (2022).From frequency ratios, it follows that the orbit considered here gradually changes its multiplicity with an increase of the pattern speed, becoming a 3:1 orbit, then a 4:1 orbit, and so on.

FREQUENCY RATIOS FOR THE SAMPLE OF GLOBULAR CLUSTERS
Here we consider in detail how the frequencies change with the pattern speed for all GCs in our sample.Fig. 11 shows the frequencies   and  R for all three potentials, and Table 4 lists the exact values.
For clarity, we investigate only three values of the pattern speed, Ω p = (30, 45, 60) km/s/kpc, while the rest of the bar parameters are fixed at  bar / b = 0.95,  bar = 5.0 kpc,  = 2.0,  = 3.0.For each orbit, we use Monte-Carlo simulations (10 3 iterations) to estimate frequency errors due to uncertainty in GCs' positions and velocities.
It can be seen from the figure that the orbital frequencies of almost all GCs behave in the same way as it was shown earlier for NGC 6266.
As the pattern speed increases, the frequencies ratio  R /   begins to deviate more and more from the resonance line 2:1.In general, Fig. 11 demonstrates that, in analytical potentials (both in BB2016 and MC2017), most GCs are do not follow the bar for all the considered values of the pattern speed.We should also note that one cannot overstep the limits of the pattern speed considered here, since they are motivated by observations.The rightmost panel of Fig. 11 shows orbital frequencies obtained for the same GCs in the -body potential.For the -body model, we do not consider different pattern speeds, since in this case its value 39 km/s/kpc follows from direct and precise measurements of the bar properties in the model (Tepper-Garcia et al. 2021).As can be seen, there is much more orbits with the resonance frequency ratio of 2:1 in such a potential.We have compiled a list of them in Table 5.Based on the orbital profiles, we distinguish the orbits into two types, the wellknown 1 family consisting of orbits elongated along the bar and supporting its structure and 2 orbits which are elongated in the direction perpendicular to the bar major axis and observed in the most central regions (Contopoulos & Papayannopoulos 1980, see also reviews by Contopoulos & Grosbol 1989 and a more recent one by Sellwood 2014).
In Fig. 12 and Table 6,  R /   values are compared for all considered potentials.The orbits themselves are presented in Fig. 13, Fig. 14, and Fig. 15.For a better comparison, we fixed the bar pattern speed in analytical potentials at the value of the pattern speed in the -body simulation, i. e. Ω p = 39 km/s/kpc for all cases.The rest bar parameters are the same as in the previous section:  bar / b = 0.95,  bar = 5.0 kpc,  = 2.0,  = 3.0.We should note that one can try to change the bar-to-bulge mass ratio somewhat to make the BB2016 and MC2017 pontentials resemble the -body potential more, but, in practice, it is hard to estimate the ratio in the -body model itself.For example, if one considers the ratio of masses of the classical bulge and the bar plus the said bulge in the -body model, it is about 0.6.However, at the same time, the total mass of the bar plus the bulge is about the half of the original disc mass (see table 1 in Smirnov et al. 2021), while, for the BB2016 potential, the bulge plus bar is about 20% of the disc.The root of the problem is that, for the -body model, the bar is formed from the disc material, and the disc itself does not go all the way towards the centre (see Smirnov & Savchenko 2020, figure 1 there).This is clearly not the case for the potentials of BB2016 and MC2017 obtained from the velocity curve fitting, where the disc goes all the way towards the centre and, thus, has high contribution in terms of mass there.One can possibly alleviate this issue by reducing the disc mass in the center or by initially considering the disc with the hole in the centre.We leave the solution of this problem for future studies.Here, we stick to the appoach by Ortolani et al. (2019a,b); Pérez-Villegas et al. (2020), where the whole or almost the whole bulge is thought to be a bar.
As can be seen,  R /   in the analytical potentials are shifted towards larger values compared to those in the -body potential, for both BB2016 and MC2017.Although, we should note that the difference between BB2016 and -body is a bit smaller on average than between MC2017 and -body.
In addition to the orbits following the bar, we want to mention some of the interesting ones with frequency ratios above or below 2:1.Liller 1 in all three potentials has a frequency ratio of about 3:2, the orbit itself looks regular, but has circle-like profile, and clearly does not follow the bar.NGC 6380 in BB2016 has  R /   close to 3, which is reflected in its overall trefoil-like shape.We should also mention, that, while most of orbits do not follow the bar in BB2017 and MC2017, some of their profiles look regular and resemble those previously shown for NGC 6266 in Fig. 7.These are NGC 6642, NGC 6558, Terzan 1, Terzan 5 for BB2016 and NGC 6380, NCC 6440, NGC 6522, NGC 6642, Terzan 1, Terzan 4, and Terzan 5 for MC 2017.It is interesting to note that most of these orbits have a rather small error in their frequency ratios (Δ  R /   = 0.1 − 0.2)

DISCUSSION
The change in the ratio of frequencies with the pattern speed has been indirectly observed in some of the previous works.In particular, Pérez-Villegas et al. (2020) found that the percentage of orbits following a bar decreases with bar rotation, except for NGC 6304, NGC 6342, and NGC 6637, which are not considered in the present work.If we assume the percentage of orbits following the bar should increase as the frequency ratio gets closer to 2:1, which is reasonable, then results of Pérez-Villegas et al. (2020) support the idea that decreasing the pattern speed causes the frequency ratios  R /   get closer to the bar frequency ratio 2:1.
A decrease in the frequency   with the patter speed, which is one of the reasons why the frequency ratio deviates from 2:1, was also observed by Sellwood & Gerhard (2020) for the self-consistent body model.Strictly speaking, what was observed is an increase of   with a decrease in Ω p (i.e.bar slow downing), which is essentially coincides as our result.We should also note that attributing the effect to a change in the pattern speed only in the case of Sellwood & Gerhard (2020) may be somewhat biased, since other properties of the bar (mass and size), were also changing there in accordance with the self-consistent evolution of the model.As for the particular GCs, Pérez-Villegas et al. ( 2020) found that for Liller 1, NGC 6304, NGC 6522, NGC 6528, NGC 6540, NGC 6553, Terzan 5, and Terzan 9, more then 20 percent of orbits follow the bar.Comparing our results to Pérez-Villegas et al. (2020), we find that, for all potentials, Liller 1 and Terzan 9 do not follow the bar, while Terzan 5 follow the bar in the -body model.NGC 6522 follows the bar in the potentials of BB2016 and MC2017, but perpendicular to it in the -body potential (2 family).For NGC 6528, the frequency ratio is close to 2:1, but the orbits themselves have an irregular profile, therefore, it is hard to say that this orbit can support a bar.

CONCLUSIONS
(i) We calculated the evolution of 30 globular clusters located in the inner area of the Galaxy ( ≲ 5) backwards in time for 5    Gyr in a non-axisymmetric galaxy potential using Gaia DR2 data for line-of sight velocities (Vasiliev 2019b) and the newest Gaia DR3 data for proper motions and distances (Vasiliev & Baumgardt 2021;Baumgardt & Vasiliev 2021).Throughout this work, we have compared the results for three potentials, two of which are analytical, obtained by fitting the rotation curve from Bajkova & Bobylev (2016, 2017) and McMillan (2017), and one is taken directly from body simulations recently preparedby Tepper-Garcia et al. ( 2021) ("surrogate Milky Way").
(ii) For all orbits, we calculated their coordinate spectra and determined the corresponding main frequencies,   and  R , for a range of bar parameters (pattern speed, mass, size, shape) in analytical potentials and for a fixed pattern speed for the -body model.
(iii) We distinguish orbits by their frequency ratio  R /   to test whether a particular orbit follows the bar.Most of orbits in both considered analytical potentials do not support the bar in the "usual" sense (either  R /   ≳ 2.1 or  R /   ≲ 1.9) for a physically reasonable values of the pattern speed, while, for the case of the -body potential, 10 GCs follow the bar (  R /  x ≈ 2.0).
(iv) On the example of one orbit (NGC 6266), we verified how the the frequency ratio changes depending on the pattern speed, the mass and size of the bar tracking the changes in a wide range of parameters  Table 5. List of GCs following the bar for Ω p = 39 km/s/kpc.Orbit types are obtained by visual classification.Note that there are no strictly periodic orbits, and we denote here the family of orbits around which an orbit with a regular profile oscillates.using a small relative step.We found the the frequency ratio does not depend much on the mass ratio of the bar and the spherical bulge ("classic" one), bar size, or its shape parameters.Most of the changes occur due to the changes in the pattern speed.For Ω p ≲ 20 km/s/kpc, the orbit perfectly follows the bar (  R /  x ≈ 2.0) for all values of the pattern speed and has a typical "bar"-like profile.Then, at a certain value of the pattern speed depending on the potential, the frequency ratio changes abruptly, becoming either greater or smaller than  R /  x ≈ 2.0.The orbit then begins to oscillate around the bar and no longer supports it.

ID
Overall, our results show that comparing orbital classifications between different potentials is indeed valuable, as the results turn out to be vastly different between them.An interesting question that we could not find an answer to in the present work, is why the -body model demonstrates a lot more bar following orbits compared to the analytic potentials.This can possibly indicate that the self-consistency of the potential should play an important factor in orbital studies of GCs.

DATA AVAILABILITY
The data underlying this article will be shared on reasonable request to the corresponding author.

Figure 1 .
Figure 1.Rotation curves for the three potentials considered in the present work.Note that the velocity is calculated for azimuthally averaged potential in case of TG2021.

Figure 2 .
Figure 2. The co-rotation (thick lines), inner Linblad resonance (thin lines below the co-rotation), and outer Lindblad resonance (thin lines above the co-rotation).The shaded area indicates the range of pattern speeds and bar sizes considered in the present work.

Figure 3 .
Figure 3. Face-on (top) and edge-on (bottom) views of the density distribution in the  -body model (left) and isolines of the potential approximation obtained from eq. 9 (right).

Figure 4 .
Figure 4.The decimal logarithm of the Jacobi energy normalised to its value at the beginning of integration for a typical orbit studied in the present work.

Figure 5 .
Figure 5.An example of the orbit of NGC 6266 in case of Ω p = 45 km/s/kpc,  bar / b = 0.95,  bar = 5 kpc,  = 2.0, and  = 3.0 (left, shown in the bar rotating reference frame).The red dot marks the initial position of the cluster.Second and third columns: coordinate time series (top row) and corresponding spectra (bottom row).

Figure 6 .
Figure 6.Dependence of the orbital frequencies   ,  R , and the ratio  R /   for NGC 6266 on the bar parameters for the potential of BB2016.

Figure 7 .
Figure 7. Evolution of the orbital profiles of NGC 6266 depending on the bar parameters for the potential of BB2016.From the first to the third row: dependencies on the pattern speed, bar-to-bulge mass ratio, and size of the bar, respectively.

Figure 8 .
Figure 8. Dependence of the orbital frequencies   ,  R , and the ratio  R /   on the bar parameters for NGC 6266 in potential of MC2017.

Figure 9 .
Figure 9. Colour coded frequency ratios  R /   for NGC 6266 depending on the bar parameters (pattern speed Ω p , bar-to-bulge mass ratio  bar / b , size  bar , and the ratios of the axes in the disc plane  and vertical direction ) for the potential of BB2016.

Figure 10 .
Figure 10.Surface of sections for the orbit of NGC 6266 with a colorbar indicating the values of the pattern speed Ω p (left) and the frequency ratio  R /   (right).

Figure 11 .
Figure 11.Frequency ratios depending on bar pattern speed for BB2016 and MC2017 potentials (left and middle subpanels, respectively), and frequency ratios in the  -body model.

Figure 12 .
Figure 12.Comparison of frequency ratios obtained for different potentials, BB2016 and  -body (left) and MC2017 and  -body (right).

Figure 13 .
Figure 13.Orbits of globular clusters for Ω p = 39,  bar / b = 0.95,  bar = 5 kpc,  = 2.0,  = 3.0 for the potential of BB2016.For all orbits, (  ) projection in the bar rotating frame is shown in the square area of 5 kpc × 5 kpc.The black line shows the orbit for the middle values of the parameters from Table2.The colour map depicts the probably of finding an orbit according to Monte-Carlo simulations.

Figure 14 .
Figure 14.Same as Fig. 13, but for the potential of MC2017.

Figure 15 .
Figure 15.Orbits of globular clusters for the  -body potential of TG2021.

Table 1 .
Observational parameters of globular clusters considered in the present work.Note that *  is the corrected value of proper motion,  *  =   cos .

Table 2 .
Cartesian coordinates and velocities of globular clusters considered in the present work.

Table 3 .
Description of the model parameters.For bar parameters, a range of values is indicated (symbol "÷" ), which is considered in the course of this work.

Table 4 .
Frequency ratios for different pattern speeds of the bar for the potential of BB2016 and MC2017.The latter are marked with an asterisk.