Issue 
A&A
Volume 639, July 2020



Article Number  A38  
Number of page(s)  14  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202038225  
Published online  03 July 2020 
Basis function expansions for galactic dynamics: Spherical versus cylindrical coordinates
^{1}
Key Laboratory of Computational Astrophysics, National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100012, PR China
email: wangyg@bao.ac.cn
^{2}
Aix Marseille Univ., CNRS, CNES, LAM, Marseille, France
^{3}
Department of Astronomy, Tsinghua University, Beijing 100084, PR China
Received:
22
April
2020
Accepted:
9
May
2020
Aims. The orbital structure of galaxies is strongly influenced by the accuracy of the force calculation during orbit integration. We explore the accuracy of force calculations for two expansion methods and determine which one is preferable for orbit integration.
Methods. We specifically compare two methods, one was introduced by Hernquist & Ostriker (HO), which uses a spherical coordinate system and was built specifically for the Hernquist model, and the other by Vasiliev & Athanassoula (CylSP) has a cylindrical coordinate system. Our comparisons include the Dehnen profile, its triaxial extension (of which the Hernquist profile is a special case) and a multicomponent system including a bar and disk density distributions for both analytical models and Nbody realizations.
Results. For the generalized Dehnen density, the CylSP method is more accurate than the HO method for nearly all inner powerlaw indices and shapes at all radii. For Nbody realizations of the Dehnen model, or snapshots of an Nbody simulation, the CylSP method is more accurate than the HO method in the central region for the oblate, prolate, and triaxial Hernquist profiles if the particle number is more than 5 × 10^{5}. For snapshots of the Hernquist models with spherical shape, the HO method is preferred. For the Ferrers bar model, the force from the CylSP method is more accurate than the HO method. The CPU time required for the initialization of the HO method is significantly shorter than that for the CylSP method, while the HO method costs subsequently much more CPU time than the CylSP method if the input corresponds to particle positions. From surface of section analyses, we find that the HO method creates more chaotic orbits than the CylSP method in the bar model. This could be understood to be due to a spurious peak in the central region when the force is calculated with the HO expansion.
Conclusions. For an analytical model, the CylSP method with an inner cutoff radius of interpolation R_{min} as calculated by the AGAMA software, is preferred due to its accuracy. For snapshots or Nbody realizations not including a disk or a bar component, a detailed comparison between these two methods is needed if a density model other than the Dehnen model is used. For multicomponent systems, including a disk and a bar, the CylSP method is preferable.
Key words: methods: analytical / Galaxy: general / galaxies: kinematics and dynamics / galaxies: formation
© ESO 2020
1. Introduction
Stellar orbits and their families are the fundamental building blocks for the dynamics of galaxies. They are also critical in constructing dynamical models of galaxies using the Schwarzschild method (e.g. Schwarzschild 1979; Rix et al. 1997; Thomas et al. 2004; van de Ven et al. 2006; van den Bosch et al. 2008; Wang et al. 2008; Wu et al. 2009, 2017; Breddels et al. 2013; Zhu et al. 2018; Magorrian 2019; Jin et al. 2019) or the madetomeasure method (e.g. Syer & Tremaine 1996; de Lorenzi et al. 2007, 2009; Long & Mao 2012; Hunt et al. 2013; Long et al. 2013; Zhu et al. 2014; Long 2016; Portail et al. 2017; Bovy et al. 2018). Moreover, orbits can be used to explain tidal streams (e.g. PriceWhelan et al. 2016; Pearson et al. 2017), the kinematics in the Solar neighborhood (e.g. Dehnen & Binney 1998a; Fux 2001), the origin of hypervelocity stars (e.g. Zheng et al. 2014; Kenyon et al. 2014; Erkal et al. 2019), the flow of gas in bars (Athanassoula 1992a,b) and the formation and structure of bars, including their boxy and peanut bulges (e.g. Patsis et al. 1997, 2002, 2003).
The studies mentioned above often require orbit integrations. For simple analytic potentials, the accuracy of orbit integration is high. However, most of these models are not sufficiently realistic for galactic studies. In more realistic models, for example for the bar in the Milky Way, the potential and force cannot be given directly in an analytical form, even for an analytic density profile. Therefore, methods for calculating the potential and force are needed even for analytical models.
A direct way to obtain the orbit properties is by using Nbody simulations. We note that Nbody simulations have successfully reproduced many properties of galaxies, but it is difficult to control the accuracy of the force calculation in a simulation. The advantage of Nbody simulations is that the potential models are selfconsistent; therefore, the orbit properties can be studied in a frozen simulation potential as long as it varies slowly. One accurate method to calculate the force in such a potential is to sum the force contribution for each particle; however, this method is not realistic if the particle number is large. An alternative is to calculate densitypotential pairs using basisfunction expansion techniques (e.g. RojasNiño et al. 2016). The Poisson’s equation can be solved in spherical coordinates (e.g. CluttonBrock 1973; Hernquist & Ostriker 1992; Zhao 1996; Rahmati & Jalali 2009; Lowing et al. 2011; Lilley et al. 2018) or in cylindrical coordinates (e.g. CluttonBrock 1972; Kuijken & Dubinski 1995; Earn 1996; Dehnen & Binney 1998b; Cohl & Tohline 1999).
For each coordinate system, there are many different basis functions. One such expansion technique was proposed by Hernquist & Ostriker (1992, hereafter HO92). Hereafter we refer to this method as the HO method. This method has been widely used for galactic bars (e.g. Zhao 1996; Wang et al. 2012; PriceWhelan et al. 2016; Pearson et al. 2017) and for Nbody potentials (e.g. Weinberg & Katz 2002; Bryan et al. 2012; Wang et al. 2013, 2019; Meiron et al. 2014; Zhu et al. 2017). For the expansion in cylinder coordinates, Vasiliev & Athanassoula (2015) have improved the speed and accuracy of the force calculation by using a spline interpolation technique. We refer to their method as the CylSP method.
Each of these expansion methods has its advantages and disadvantages, and some recent studies already gave comparisons between the sphericalharmonic expansion method and the CylSP method (Vasiliev & Athanassoula 2015; Vasiliev 2019). There are many parameters in each method, and the accuracy of each method depends on the values chosen for these parameters, as well as on the density distribution and the particle number. In this paper, we give here a more detailed comparison between the HO and CylSP methods. The accuracy of the force calculation can affect the orbit shape, and the resulting orbit families in the systems (e.g. Carpintero & Wachlin 2006). In this paper we also study how the expansion methods affect the orbit families. The outline of the paper is as follows. In Sect. 2 we describe the two expansion methods for the potential and accelerations. In Sect. 3 we compare their accuracy for a triaxial extension of the Dehnen models (Dehnen 1993). In Sect. 4 we present our comparisons for a multicomponent model which includes a Ferrers bar. In Sect. 5 we present the summary and discussions.
2. Potential and accelerations
Given a model density distribution, the potential can be calculated using Poisson’s equation. Then the accelerations of the model can be obtained from the derivatives of the potential. For complicated density distributions, the potential and density can be expanded on a set of simple orthogonal basisfunction in spherical coordinates
where n is the radial expansion parameter, and l and m are the usual angular parameters. There are many basis functions and different expansion formats in the spherical coordinates (e.g. CluttonBrock 1973; Hernquist & Ostriker 1992; Zhao 1996; Weinberg 1999; Wang et al. 2012; Vasiliev 2013). Here we focus on the HO method. The potential and accelerations can be written as Eqs. (3.13) and (3.21)–(3.23) in HO92, respectively. In the HO method, the expansion coefficients are determined by threedimensional integrations (see Eq. (3.14) in HO92). For a given system, the expansion coefficients are calculated only once, and the potential and force in any arbitrary position can be obtained using these coefficients.
For a nearly disk system, the density and potential pairs are naturally expanded in cylindrical coordinates (R, z, ϕ). The density and potential can be represented by the Fourier series
and each term in the expansion is given by
and
where m_{max} is the truncated number of azimuthal Fourier harmonic terms. For an axisymmetric model, only the m = 0 term is needed, while m_{max} ≲ 12 is enough for typical asymmetric models (Vasiliev 2019). We note that J_{m} is an mthorder Bessel function of the first kind, and can be expressed analytically by the halfinteger degree Legendre function of the second kind Q_{m − 1/2} (Gradshteyn & Ryzhik 1965; Cohl & Tohline 1999; HolleyBockelmann et al. 2005)
Equation (6) shows that the potential can be obtained by a twodimensional integration, and this integration depends on the (R, z) coordinates of the system. Therefore, the calculation is time consuming if the input model is analytical. In order to increase the speed for the calculations of potential and forces, Vasiliev & Athanassoula (2015) improved this expansion using a spline interpolation. Therefore, the potential and forces can be obtained only in some fixed grids (N_{R} × N_{z}) within a fixed region R ∈ [R_{min}, R_{max}], z ∈ [z_{min}, z_{max}], and then the potential and forces in any position can be derived by the spline interpolation. There are two advantages in using the spline interpolation. First, the speed for the orbit integration can be improved by several orders of magnitude. Second, the unphysical noise for the force can be smoothed if the input model is a snapshot of an Nbody simulation.
3. Generalized Dehnen models
3.1. General density
The density models we consider here are triaxial and their spherical analogs were first discussed by Dehnen (1993). The density model is (Merritt & Fridman 1996)
with
Above, M is the total mass of the model and γ is the inner density profile index. We refer to this kind of density distribution as the generalized Dehnen models or simply as Dehnen models. Specifically, the γ = 1 case is the same as the Hernquist (1990) model. The isodensity contours are ellipsoids with axis ratios of a:b:c. We adopt four combinations of a:b:c = 1:1:1, 1:1:0.2, 1:0.2:0.2, and 1:0.79:0.5, corresponding to spherical, oblate, prolate, and triaxial models, respectively. Merritt & Fridman (1996) analytically reduced the triple integrals of the potential and force calculations to a single integral each. This is of course a considerable improvement as the corresponding calculations are much easier to perform numerically and, moreover, are much more accurate. We thus use them here, labeling them as “true” values. For a spherical model, the potential and forces have analytical expressions (see Eqs. (4) and (6) in Merritt & Fridman 1996). We have checked the accuracy for the true values obtained using onedimensional integrals relative to analytic values in the spherical model, and we found that they match each other within a relative accuracy of 10^{−15} for the model with γ = 1 and better than 10^{−8} for the other values of γ. It is important to notice that the xaxis, is the major axis while the zaxis is the minor axis. Henceforth, we use units such that the total mass M, the major axis length a, and the gravitational constant G are equal to unity.
3.2. Comparisons for analytic potentials
In Fig. 1, we show the potential and accelerations of the Hernquist models for the true values and the results from the HO method. It is seen that the HO method can approximate the potential well for a spherical system. For oblate and prolate systems, the HO method does not provide accurate forces in the central region. Naturally, larger values of n_{max} and l_{max} give more accurate expansions than smaller ones for most models and in most regions. For a spherical system, we find that large values of n_{max} and l_{max} generate larger errors in the force along the minor axis (a_{z}) in the central region than those only using the lower order terms of coefficients, this is reasonable and consistent with that in HO92. In the HO method, the zeroth order terms in the expansion of the density and potential are the Hernquist distribution and its corresponding potential. Therefore, if only the zeroth terms are adopted in the spherical model, then the errors in both the forces and potential are zero. In order to give a consistent comparison for different shapes, we also kept the highorder terms in the expansion. Also, for the triaxial system, if n_{max} = 16 and l_{max} = 12, then the HO method can give accurate expansions.
Fig. 1. Potential and forces for the Hernquist models. From left to right: the input density model is spherical, oblate, prolate, and triaxial, respectively. In each panel, the black line is the result for the model given by Merritt & Fridman (1996). The red, blue, and green lines are the results for the HO method with (n_{max}, l_{max}) = (16,12), (12,12), and (6,4), respectively. From top to bottom: the radius r was measured along the x, x, y, and zaxes, respectively. The xaxis is the major axis and the zaxis is the minor axis of the density distribution. 

Open with DEXTER 
In Fig. 2, we show the potentials and accelerations of the Hernquist models, comparing the true values and the results from the CylSP method. Here we consider different values of R_{min} in this method. In the AGAMA software implementation (Vasiliev 2019) of the CylSP method that we are using, the code automatically chooses the R_{min} value and in the case of our models gives a few times 10^{−4}. It is clear that this method can reproduce the potential well for spherical, oblate, prolate, and triaxial systems. For the forces from the CylSP method, the value of R_{min} that is given by the software can give accurate values in all cases except in the centralmost region (r < 2 × 10^{−3}) for the spherical, prolate, and triaxial models. For the forces from R_{min} = 0.01 and R_{min} = 0.1, obvious wiggles appear in the central region (r < 0.1).
Fig. 2. Same as in Fig. 1, but for the results of the CylSP method. In each panel, the red, and blue lines are the results of the CylSP method with R_{min} = 0.1, and 0.01, respectively, while the green lines correspond to the value calculated by the AGAMA software (see Sect. 3.2). The black line is the analytical result given by Merritt & Fridman (1996). Here N_{z} = N_{R} = 40, z_{min} = R_{min}, and R_{max} = z_{max} = 100, m_{max} = 12. 

Open with DEXTER 
We also consider different values of m_{max} and N_{R}, N_{z} for the CylSP method. It should be noted that N_{R} = 40 can give good accuracy as shown in Vasiliev & Athanassoula (2015). We then searched for the optimum m_{max} value and found that there is hardly any improvement beyond m_{max} = 12, while the necessary CPU time is considerably increased. For m_{max} ≥ 12, the accuracy for the force and potential is almost unchanged even when m_{max} is further increased together with CPU time. Therefore, we adopt N_{R} = N_{z} = 40 and m_{max} = 12 in the CylSP method in this paper unless stated otherwise.
In Fig. 3, we show the relative error in the potential and forces for both the HO and CylSP methods for the analytical Hernquist models. We note that the CylSP method is more accurate than the HO method for nearly all cases, all shapes, and all radial ranges.
Fig. 3. Relative error in the potential and forces for both the HO (red lines) and CylSP (blue lines) methods for the analytical Hernquist models. For the HO method, (n_{max}, l_{max}) = (16,12). For the CylSP method, R_{min} is given by the AGAMA software, N_{R} = 40, N_{z} = 40, z_{min} = 0, R_{max} = z_{max} = 100. 

Open with DEXTER 
Figure 4 shows the relative errors in the potential and forces for the models with a different inner density profile index γ, using the HO method for the small radius r = 0.01. It is important to notice that the accuracy of the HO method depends on γ. In the spherical, oblate, and triaxial systems, the accuracy of γ = 1 is, in general, higher than for the other γ values, with the exception of the z component of the force in the triaxial model. This is not surprising since the HO method is designed for γ = 1. For models with another value for γ, the HO method is not accurate for either the potential or the force; this is especially the case for a low number of expansion coefficients. The deviation from the exact values is 10%−100% in the force for modes with γ ≠ 1. This deviation is due to the HO method itself, and it was also presented in HO92. For the Hernquist model, the accuracy for the spherical distribution is higher than that for the prolate one. We note that there is a significant difference of the accuracy between γ = 1.9 and 2 for the triaxial model with (n_{max}, l_{max}) = (16,12) in the potential. These relative errors depend on the accuracy of the true value and the HO results, and their accuracy depends on the density distribution of the model; even two models with similar values for γ may exhibit a significant difference. Moreover, we also tested the relative errors in the force for the triaxial models with γ = 1.9, 1.92, 1.94, 1.96, and 2 with (n_{max}, l_{max}) = (16,12), and we found that these errors between two values for γ that are close can differ quite markedly.
Fig. 4. Relative errors of the potential and forces for models with a different γ value for r = 0.01. From the left to right: the input density model is spherical, oblate, prolate, and triaxial, respectively. The red, blue, and green lines are the results for the HO method with (n_{max}, l_{max}) = (16,12), (12,12), and (6,4), respectively. 

Open with DEXTER 
Figure 5 shows the relative errors in the potential and forces for the models with a different value for γ, using the CylSP method for r = 0.01. It is seen that the CylSP method can approximate the potential and forces well with γ = [0, 2] if R_{min} = 0.01, or for the small value given by the AGAMA software. This shows that 0.01 and the value provided by AGAMA are the appropriate values for R_{min} for the Dehnen models, while R_{min} = 0.1 is not. Compared with the accuracy from the HO method, the CylSP method performs considerably better at r = 0.01, that is, close to the center. In the spherical and oblate models, we found that there is a significant difference in the accuracy between γ = 1.9 and 2 for R_{min} provided by AGAMA and 0.1, while this difference is small for R_{min} = 0.01. It is known that forces from the CylSP method are obtained by spline interpolation and that their values are strongly related to the positions of the grids. Moreover, the grids are divided by R_{min}, R_{max} and N_{R}, and each grid cell has roughly the same mass (Vasiliev 2013). In other words, the grids depend on the mass distribution of the model. Therefore, even when R_{min}, R_{max}, and N_{R} are the same for two different models, the grid positions can be quite different. This is the reason why we find significant differences in the forces for models γ = 1.9 and 2 between R_{min} = 0.1 and the value provide by AGAMA. Furthermore, we checked the relative errors in the forces for the spherical and oblate models with γ = 1.9, 1.92, 1.94, 1.96, and 2 with R_{min} = 0.1 and the value given by AGAMA, and we found that the relative errors between two similar γ values can also be different. Although, both the HO and CylSP methods have sudden jumps from models with γ = 1.9 to 2; the amplitude of the jump can be much bigger for the HO method than for the CylSP method. For the HO method, it can reach on the order of 100% for n_{max} = 6, l_{max} = 4; whereas, in the CylSP method, it can go between 10^{−3} and 10^{−4}, or, in the worst case scenario, between 10^{−3} and 10^{−2}. Thus, here as well the CylSP method performs much better than the HO method.
Fig. 5. Same as in Fig. 4, but for the CylSP method. The red, blue, and green lines represent the results with R_{min} = 0.1, 0.01, and the value provided by AGAMA, respectively. 

Open with DEXTER 
The left panels of Fig. 6 show the relative error of the potential and forces for prolate models with a different c/a by using the HO method for the small radius r = 0.01. We find that the HO method can give a better estimation for the potential than for the forces. For the latter, the HO method gives more accurate values for larger values of c/a. The reason is that the model becomes more spherical with an increasing c/a and that a spherical system has a higher accuracy than a prolate one, as is shown in Fig. 4.
Fig. 6. Relative error in potential and forces for both the HO (left) and the CylSP methods (right) for the prolate models with a different c/a at r = 0.01. Left panels: the red, blue, and green lines are the results for the HO method with (n_{max}, l_{max}) = (16,12), (12,12), and (6,4), respectively. Right panels: the red, blue, and green lines represent the results for the CylSP method with R_{min} = 0.1, 0.01, and the value given by AGAMA, respectively. From top to bottom: the radius r is along the x, x, y, and zaxes, respectively. 

Open with DEXTER 
The right panels of Fig. 6 show the relative error of the potential and forces for prolate models with a different c/a by using the CylSP method for the small radius r = 0.01. It is found that the CylSP method can give accurate values for both the potential and forces in the prolate model at r = 0.01 and the accuracy is nearly the same with c/a ≥ 0.4 in the potential, a_{x} and a_{y}. We can thus conclude that the CylSP method is more accurate than the HO method at r = 0.01.
3.3. Comparisons for particle realizations of generalized Dehnen potentials
In order to study the accuracy of the two expansion methods for particle realizations of given potentials, we used the inverse function sampling method to generate the Dehnen models with spherical, oblate, prolate, and triaxial shapes as well as different particle numbers. Figure 7 shows the relative errors of the potential and forces for the Hernquist models with different particle numbers using the HO method at radius r = 0.01. We randomly generated 40 samples with different particle numbers. The filled circles in Fig. 7 are the mean value and the error bars are the standard error from these 40 samples. We also note that the error bar size in each point is different. One possible reason is that there is a high Possion noise level if the particle number is smaller than 10^{6}, particularly in the central region. In the central region r ≤ 0.01, the density is high but the total enclosed mass is low, resulting in the Poisson noise due to the small number of particles sampled. The second reason is that the error bar size also depends on the mean value of the deviation since we use a logarithmic scale for the xaxis.
Fig. 7. Relative errors of the potential and forces for the Hernquist models with different particle numbers at r = 0.01 when using the HO method. In each panel, the red, blue, and green lines are the results for HO method with (n_{max}, l_{max}) = (16,12), (12,12), and (6,4), respectively. The filled circles are the mean value from 40 random samples, and the error bar is the standard error from these samples. 

Open with DEXTER 
It is seen that the forces obtained from the HO method deviate from the true value less than 2% if N_{p} ≥ 5 × 10^{5}, and this deviation is nearly independent of the number of expansion coefficients n_{max} and l_{max}. For the oblate model, the force deviations for the HO method with (n_{max}, l_{max}) = (16,12) are smaller than ∼9% and these deviations weakly depend on the particle number. It is also noted that the forces along the zaxis with (n_{max}, l_{max}) = (12,12) are more accurate than those with (n_{max}, l_{max}) = (16,12). As shown in HO92, increasing the number of expansion terms can improve the accuracy of the system as a whole, the reconstructed acceleration with a higher number of expansion coefficients can deviate more from the corresponding exact values at the central region (HO92). This is one possible reason as to why we find the results from a lower number of expansion coefficients at r = 0.01 that are more accurate than those from a higher number of expansion coefficients. Another possible reason is that there is a high Poisson noise level in the central region. The prolate model has the largest force deviation of these four shapes. These deviations are ∼20% for (n_{max}, l_{max}) = (16,12) and these values are nearly constant along the particle numbers. The triaxial model is more accurate than the oblate and prolate models, but less accurate than the spherical one. The force deviations are ∼3% for (n_{max}, l_{max}) = (12,12) at r = 0.01 and they depend only weakly on the particle number.
Figure 8 shows the relative errors of the potential and forces for the Hernquist models with different particle numbers using the CylSP method at radius r = 0.01. It is clear that the force accuracy for the R_{min} provided by the AGAMA software is higher than that for R_{min} = 0.01 and 0.1; additionally, this accuracy increases with the particle number. This is because the CylSP method uses the spline interpolation method, and the accuracy from the interpolation is higher than that from the extrapolation. A larger number of particles can reduce the Poisson noise of the system and reproduce the density distribution of the system better.
Fig. 8. Same as in Fig. 7, but for the CylSP method. The red, blue, and green lines represent the results in the CylSP method for R_{min} = 0.1, 0.01, and the value given by AGAMA, respectively. 

Open with DEXTER 
We see that the force accuracy for the CylSP method with R_{min} as given by the AGAMA software is equivalent to that of the HO method for the spherical and triaxial models, and higher than the HO method in the oblate and prolate models with N_{p} ≥ 5 × 10^{5}. Compared to the analytical model, we can see that the advantages of the CylSP method over the HO method are weakened for the particle samples, as was also found by Vasiliev & Athanassoula (2015). The reason is that, in the analytical models, the HO method expansion coefficients are determined by a threedimensional integration while the expansion coefficients in the CylSP method are obtained from a twodimensional integration. For particle realizations, however, the integration is unnecessary for both methods in calculating the expansion coefficients.
We also compared the accuracy for these two methods at a smaller inner radius, namely r = 0.001. We found that the accuracy for forces for the HO method is higher than that for the CylSP method in the spherical and triaxial models and this trend is independent of the particle number. For the oblate and prolate models, the CylSP method is more accurate than the HO method if N_{p} ≥ 10^{6}.
Except for the force accuracy, the required CPU time is also important for the orbit integration. We consider both the initialization time, which was performed only once, and the subsequent time. In the HO method, the initialization part calculates the expansion coefficients, while both the expansion coefficients and the force on the grid are calculated in the CylSP method. For an analytical model, the expansion coefficients are determined by threedimensional integrations in the HO method, while those are calculated by twodimensional integrations in the CylSP method. If the CPU time for the force on the fixed grid is also considered, the initialization time for the CylSP method can be comparable to that in the HO method for the analytical model. In this paper, we study both the initialization time for an Nbody realization and the time after initialization to calculate the force and acceleration.
On the left side of Fig. 9, we show the dependence of the initialization CPU time on the particle number from the HO method. It is clear that the initialization time is proportional to the particle number and independent of the model shape. This is normal because the expansion coefficients in the HO method are the sum of a function of the particle position (see Eq. (3.17) in Hernquist & Ostriker 1992). In the HO method, the subsequent time is used to sum the expansion terms, which are related to the number of the expansion coefficients . For both the analytical model and Nbody particles, once the expansion coefficients have been obtained, the subsequent CPU time is the same if the same N_{EC} is adopted. The right panel of Fig. 9 shows the dependence of the subsequent CPU time on the number of the expansion coefficients for the HO method. We find that the subsequent time in the HO method increases significantly with the number of the expansion coefficients, which can be explained by Eq. (3.17) in HO92.
Fig. 9. CPU time for the HO method. Left: dependence of the initialization CPU time on the particle number of a realization of the Hernquist model. The red, blue, green, and black lines are the results for the spherical, oblate, prolate, and triaxial models, respectively. The blue, green, and black lines are slightly shifted along the xaxis. Right: dependence of the subsequent CPU time on the number of the expansion coefficients. The subsequent CPU time is estimated by the forces at 10 000 positions. 

Open with DEXTER 
Figure 10 shows the dependence of the CPU time on the parameters of particle number N_{p} (top left), and N_{R} (top right) and m_{max} (two bottom panels) for the CylSP method. The bottom right panel is for the subsequent calculation time while the other three panels are for the initialization CPU time. We find that the initialization CPU time increases significantly with the particle number N_{P} and the grid number N_{R}. It is noted that the initialization CPU time in the oblate and spherical models is significantly smaller than that in the prolate and triaxial models. The reason is that only the m = 0 terms are needed for a disk system, with the spherical system being a special “disk” system. This is the reason why we find the initialization CPU time to be a constant in the spherical and disk models, which can be seen in the bottom left panel of Fig. 10.
Fig. 10. CPU time for the CylSP method. Top left: dependence of the initialization CPU time on the particle number of the Hernquist model realization. Top right: dependence of the initialization CPU time on the parameter N_{R} for an Nbody realization of the Hernquist model with N_{p} = 10^{6}. Bottom left: dependence of the initialization CPU time on the parameter m_{max} for a snapshot of an Nbody realization of the Hernquist model with N_{p} = 10^{6}. Bottom right: dependence of the subsequent CPU time on the parameter m_{max}. The subsequent CPU time is estimated by the forces at 10 000 positions. Excluding the bottom right panel, the blue and black lines are slightly shifted to the right along the xaxis and the red, blue, green, and black lines are the results of the spherical, oblate, prolate, and triaxial models, respectively. 

Open with DEXTER 
Compared with the HO method, the initialization CPU time in the CylSP is high. However, the subsequent time in the CylSP method is much lower than that in the HO method. The reason is that the subsequent time in the HO method is the cost of summing N_{CE} terms of the basis function; whereas for the CylSP method, it is only the time to use a twodimensional interpolation routine (Vasiliev 2019). Therefore, if a large number of orbits need to be calculated, such as millions of orbits, this is something that should be kept in mind.
4. A multicomponent model with a Ferrers bar
Since nearly twothirds of spiral galaxies in the Universe are barred (e.g. Lee et al. 2012; Buta et al. 2010, 2015), it is important to discuss the application of the different expansion methods to barred galaxies. Here we consider a fiducial barred galaxy model, which consists of a Miyamoto (MN) disk (Miyamoto & Nagai 1975), a Plummer halo (Plummer 1911) and a Ferrers bar (Ferrers 1877). The parameters we used are those listed in Table 1 in Skokos et al. (2002) and, for reasons of continuity, we have the same bar orientation, that is, along the y axis. The potential and force for this bar model can be obtained directly by following the appendix in Pfenniger (1984).
Figure 11 shows the potential and forces for the bar models using the HO method. The input model is an analytical model. It is seen that the HO method can reconstruct the potential well, except for the potential in the central region (r < 0.1) of the Ferrers bar and the MN disk. For the Plummer halo, the HO method can give a good approximation if a large number of expansion coefficients are used.
Fig. 11. Potential and forces for the barred galaxy models and their individual components. From left to right: the input density model is the Ferrers bar, the MN disk, the Plummer halo, and the total, i.e., bar + disk + halo, respectively. In each panel, the black line corresponds to the results for the model given by the analytical method directly. The red, blue, and green lines are the results for the HO method with (n_{max}, l_{max}) = (16,12), (12,12), and (6,4), respectively. 

Open with DEXTER 
Figure 12 also shows the potential and forces for the models using the HO method, but now the input is an Nbody realization of the model. It shows again that the HO method can provide large errors for the force in the central region, both for the Ferrers bar and the MN disk. It is also noted that the potential of the MN disk, as calculated by the HO method, significantly deviates from the true value. One possible reason was that our input Nbody snapshot is sampled using the revised function method. We have compared the density from the sampled particles with the analytical model; the density at r < 1 has a 5% deviation for the sampled particles, as was also found by Vasiliev & Athanassoula (2015).
Fig. 12. Same as in Fig. 11, but the only difference is that the input is an Nbody realization of our barred galaxy model. The total particle number for the bar is 2 × 10^{6}. The particle numbers for the bar, disk, and halo are 200 000, 1640 000, and 160 000, respectively. 

Open with DEXTER 
Figure 13 shows the potential and forces for the analytical barred galaxy models using the CylSP method. Clearly, CylSP can reproduce the potential and forces well for the Ferrers bar, MN disk, and the Plummer halo. We also find that the results from different R_{min} are nearly same; the reason is that the potential and accelerations vary slowly in the central region in these models.
Fig. 13. Same as in Fig. 11, but for the CylSP method. The red, blue, and green lines represent the result in the CylSP method with R_{min} = 0.1, 0.01, and the value provided by AGAMA, respectively. Since the results from three R_{min} are the same as the true values, only the black lines can been seen in most panels. 

Open with DEXTER 
Figure 14 shows the potential and forces for the barred galaxy models by using the CylSP method. The input is the particle realization with particle number N_{p} = 10^{6}. Generally, we find that the CylSP method is more accurate than the HO method, especially in the central region of the system. However, there are some wiggles in the force distribution in the Ferrers bar, MN disk, and Plummer halo, and the results from R_{min} = 0.01 and the R_{min} provided by AGAMA give larger wiggles than those from R_{min} = 0.1. These wiggles affect the orbit shape in the system, especially the box orbits. To determine the reason why these wiggles occur, we increased the number of particles in the Nbody model. As is shown in Fig. 15, the wiggles become weaker and weaker with the increasing particles and the wiggles vanish if the particle number is 10^{7}.
Fig. 14. Same as in Fig. 13, but the only difference is that the input is an Nbody realization of our barred galaxy model. The total number of particles is 2 × 10^{6}. 

Open with DEXTER 
Fig. 15. Potential and forces for the Ferrers bar with different redistribution numbers of particles. In each panel, the black line is the results for the model given by the analytical method. The red, blue and green lines represent the result in the CylSP method with R_{min} = 0.1, R_{min} = 0.01 and the value provided by AGAMA, respectively. From left to right: the results for Nbody realizations with particle number 2 × 10^{5}, 5 × 10^{5}, 10^{6}, and 10^{7}, respectively. 

Open with DEXTER 
For the Ferrers bar model, we have found that the CylSP method is more accurate than the HO method along the radius even when the expansion coefficient terms with n_{max} = 16, l_{max} = 12 are adopted. As is shown in Fig. 16, the accuracy of the HO method can be comparable to the CylSP method at large radii (r > 0.1) if n_{max} = 40, l_{max} = 40 are used. However, the accuracy in the central region is not improved for the HO method.
Fig. 16. Relative errors of potential and forces for the HO (red lines) and CylSP (blue lines) methods for snapshots of an Nbody realization of a Ferrers bar with a particle number N_{p} = 10^{6}. From left to right: the results for the HO with (n_{max}, l_{max}) = (16,12), (20,20), (30,30) and (40,40), respectively. For the CylSP method, R_{min} = 0.1, N_{R} = 40, N_{z} = 40, z_{min} = 0.1, R_{max} = z_{max} = 100. 

Open with DEXTER 
Since the accuracy of the forces from the HO and CylSP methods is different, this will affect the accuracy of the orbit integration, and thus the orbit structure. One useful way to characterize the orbit structure is the surface of section (SOS). In the bar model, we follow Skokos et al. (2002) and consider the yaxis to be the longest major axis. Therefore, we consider the intersection of the orbit with the plane y = 0 and in particular those with v_{y} > 0. For each energy, we generated 25 initial conditions of orbits with x ≠ 0 and v_{y} ≠ 0.
Figure 17 shows the SOS for four Jacobi’s energies from the analytical potential, the HO method and the CylSP method, respectively. In order to focus on studying how the accuracy of the force from the bar affects the SOS, we use analytical force and acceleration calculations for the MN disk and the Plummer halo. We also consider two different particle numbers of the bar to probe the discreteness noise. It is seen that the overall framework of the SOS at each E_{J} from three methods are similar, but the differences in the details are also obvious.
Fig. 17. Surfaces of section (SOS) for the bar model. Left column: results for the analytic potential. Second and third columns: results for the HO method with N_{p} = 2 × 10^{5} and N_{p} = 10^{7}, respectively. Fourth and fifth columns: results for the CylSP method with N_{p} = 2 × 10^{5} and N_{p} = 10^{7}, respectively. From top to bottom: the rows correspond to E_{J} = −0.41, −0.39, −0.37, and −0.35, respectively. The red filled circle in each panel is the position of the x1 orbit. 

Open with DEXTER 
For example, around the x1 orbit, more chaotic orbits are generated from the HO method than from the true and CylSP method. This can be due to higher numerical errors in the central region, which we already discussed above. This can introduce a discontinuity and therefore chaos. It is also clear that the discreteness noise can be reduced with increasing particle number in both methods. However, even when the particle number is 10^{7}, the HO method still has large errors in the central region, and thus more chaotic orbits appear in the central region than those in the true and CylSP method. For the CylSP method, the Nbody bar model with N_{p} = 10^{7} can reconstruct almost the same SOS as the true potential, which again shows that the CylSP method works better than the HO method for a bar. In Table 1, we show that the stability of the x1 orbit at four different E_{J}. We found that the stability of the x1 orbit at these E_{J} are the same for the three methods. Since we fixed the disk and halo potential, and as the particle mass of the Ferrers bar only contributes 10% of the mass for the whole system, the effect of stability due to the force accuracy for the three methods is not obvious. We also tested the stability of the x1 orbit at E_{J} = −0.37 and E_{J} = −0.35, which could show considerable differences between the three methods if the disk and halo were also represented by particles.
Four Jacobi energies E_{J} in the bar potential.
5. Summary and discussions
The accuracy of the potential and force calculations directly affects the orbit properties in a model. In this paper, we present a detailed comparison of the Hernquist & Ostriker (1992, HO) method in the spherical coordinate system with the CylSP method (Vasiliev & Athanassoula 2015), which uses the twodimensional spline interpolation for Fourier coefficients in the cylindrical coordinate system. We have compared these two methods for two kinds of density distribution; one has a Dehnen profile (Dehnen 1993), which is generalized to nonspherical shapes (see Eq. (10)). If the inner density powerlaw index γ is equal to 1, this is the Hernquist model. The second kind of model that we consider is a multicomponent galaxy model, which consists of a Miyamoto disk, a Plummer halo, and a Ferrers bar.
For the Hernquist model (with γ = 1), we find that the CyLSP method is more accurate than the HO method. This is especially the case in the central region if the input model is an analytical density distribution and a reasonable parameter set is chosen in the CyLSP method. However, if the input density is an Nbody realization, we find that the force accuracy in the CylSP method is more sensitive to the particle number than the HO method. The force accuracy from the CylSP method at the central region is equivalent to the HO method for the spherical and triaxial models, and it is only higher than the HO method in the oblate and prolate models if the particle number is N_{p} ≥ 5 × 10^{5}. If the particle number is smaller than 5 × 10^{5}, the HO method is better than the CylSP method for the spherical and triaxial models.
For the generalized Dehnen model with γ ≠ 1, the force accuracy from the HO method is lower than that for the Hernquist model. The CylSP method is more accurate than the HO method if the input density is the analytical Dehnen model. For the snapshot of an Nbody realization, we also checked the force accuracy for these two methods at the central region. The advantage of the CylSP method to the HO method is significant for the oblate shapes. For the more prolate Dehnen model, the CylSP method is accurate than the HO method. For the spherical and triaxial models, the HO method is significantly better than the CylSP method if the particle number is smaller than 5 × 10^{5}. However, for the spherical and triaxial Dehnen models, the CylSP method is significantly more accurate than the HO method if the particle number is 5 × 10^{6}.
For the multicomponent galaxy model with a Ferrers bar, MN disk and the Plummer halo, the HO method has a significant cusp in the central region of the bar, while the CylSP method gives good results for both the analytical model and the Nbody realization. If we consider the bar component of this Nbody realization separately, the CylSP method can also generate some wiggles in the central region for the force if a small R_{min} is used, but these wiggles vanish with increasing particle numbers. The force error in the central region from the HO method can affect the orbit structure significantly, especially for box orbits, which approach very closely to the center, thus affecting the SOS. We find more chaotic orbits if the HO method is used and also that the stability of the x1 orbits from these two methods may be different. The effect on detailed orbit families and a frequency analysis of orbits from these two methods will be considered in a future study.
If the input model is an analytical model, the CylSP method with R_{min} as provided by the AGAMA software is preferred. For models with snapshots of an Nbody simulation, the initialization CPU time for the HO method is significantly smaller than that for the CylSP method, while the subsequent CPU time for the CylSP method is much lower than that for the HO method. For models with snapshots of an Nbody simulation, if the model is a spherical or nearly spherical Hernquist model, the HO method can be preferentially used.
For a multicomponent model for a barred galaxy, as expected, the CylSP method works better. A detailed comparison between these two methods for the particle realization is still needed if a density model other than the Dehnen model studied here is used. If the CylSP method is selected for snapshots of an Nbody simulation, a relatively large R_{min} (R_{min} = 0.1) is suggested.
Acknowledgments
We thank E. Vasiliev and JeanCharles Lambert for many useful discussions and for help with software. E. Vasiliev also gave us insightful comments, which greatly improved the paper. The calculations for the CylSP method are made using the publicly available AGAMA software by Vasiliev (2019), and the results of for snapshots of Nbody simulations for the HO method is performed by the SCF.F software provided by Lars Hernquist. We acknowledge the support by the National Key R&D Program (No. 2018YFA0404501, 2017YFA0402603), the National Science Foundation of China (Grant No. 11821303, 11773034, 11761141012, 11761131004 and 11633004), the Chinese Academy of Sciences (CAS) Strategic Priority Research Program XDA15020200 and the CAS Interdisciplinary Innovation Team (JCTD 201905). The “Centre de Calcul Intensif d’AixMarseille” is acknowledged for granting access to its high performance computing resources, where the initial part of this work was performed. Most of the computing was performed on the high performance computing cluster at the Information and Computing Center at National Astronomical Observatories, Chinese Academy of Sciences. This work is also supported by Astronomical Big Data Joint Research Center, cofounded by National Astronomical Observatories, Chinese Academy of Sciences and Alibaba Cloud.
References
 Athanassoula, E. 1992a, MNRAS, 259, 328 [NASA ADS] [CrossRef] [Google Scholar]
 Athanassoula, E. 1992b, MNRAS, 259, 345 [NASA ADS] [CrossRef] [Google Scholar]
 Bovy, J., Kawata, D., & Hunt, J. A. S. 2018, MNRAS, 473, 2288 [CrossRef] [Google Scholar]
 Breddels, M. A., Helmi, A., van den Bosch, R. C. E., van de Ven, G., & Battaglia, G. 2013, MNRAS, 433, 3173 [NASA ADS] [CrossRef] [Google Scholar]
 Bryan, S. E., Mao, S., Kay, S. T., et al. 2012, MNRAS, 422, 1863 [NASA ADS] [CrossRef] [Google Scholar]
 Buta, R. J., Sheth, K., Regan, M., et al. 2010, ApJS, 190, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Buta, R. J., Sheth, K., Athanassoula, E., et al. 2015, ApJS, 217, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Carpintero, D. D., & Wachlin, F. C. 2006, Celest. Mech. Dyn. Astron., 96, 129 [CrossRef] [Google Scholar]
 CluttonBrock, M. 1972, Ap&SS, 16, 101 [CrossRef] [Google Scholar]
 CluttonBrock, M. 1973, Ap&SS, 23, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Cohl, H. S., & Tohline, J. E. 1999, ApJ, 527, 86 [NASA ADS] [CrossRef] [Google Scholar]
 de Lorenzi, F., Debattista, V. P., Gerhard, O., & Sambhus, N. 2007, MNRAS, 376, 71 [NASA ADS] [CrossRef] [Google Scholar]
 de Lorenzi, F., Gerhard, O., Coccato, L., et al. 2009, MNRAS, 395, 76 [NASA ADS] [CrossRef] [Google Scholar]
 Dehnen, W. 1993, MNRAS, 265, 250 [NASA ADS] [CrossRef] [Google Scholar]
 Dehnen, W., & Binney, J. J. 1998a, MNRAS, 298, 387 [NASA ADS] [CrossRef] [Google Scholar]
 Dehnen, W., & Binney, J. 1998b, MNRAS, 294, 429 [NASA ADS] [CrossRef] [Google Scholar]
 Earn, D. J. D. 1996, ApJ, 465, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Erkal, D., Boubert, D., Gualandris, A., Evans, N. W., & Antonini, F. 2019, MNRAS, 483, 2007 [NASA ADS] [CrossRef] [Google Scholar]
 Ferrers, N. M. 1877, Quart. J. Pure Appl. Math., 14, 1 [Google Scholar]
 Fux, R. 2001, A&A, 373, 511 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gradshteyn, I. S., & Ryzhik, I. M. 1965, Table of Integrals, Series and Products (New York: Academic Press) [Google Scholar]
 Hernquist, L. 1990, ApJ, 356, 359 [NASA ADS] [CrossRef] [Google Scholar]
 Hernquist, L., & Ostriker, J. P. 1992, ApJ, 386, 375 [NASA ADS] [CrossRef] [Google Scholar]
 HolleyBockelmann, K., Weinberg, M., & Katz, N. 2005, MNRAS, 363, 991 [NASA ADS] [CrossRef] [Google Scholar]
 Hunt, J. A. S., Kawata, D., & Martel, H. 2013, MNRAS, 432, 3062 [NASA ADS] [CrossRef] [Google Scholar]
 Jin, Y., Zhu, L., Long, R. J., et al. 2019, MNRAS, 486, 4753 [CrossRef] [Google Scholar]
 Kenyon, S. J., Bromley, B. C., Brown, W. R., & Geller, M. J. 2014, ApJ, 793, 122 [NASA ADS] [CrossRef] [Google Scholar]
 Kuijken, K., & Dubinski, J. 1995, MNRAS, 277, 1341 [CrossRef] [Google Scholar]
 Lee, G.H., Park, C., Lee, M. G., & Choi, Y.Y. 2012, ApJ, 745, 125 [NASA ADS] [CrossRef] [Google Scholar]
 Lilley, E. J., Sanders, J. L., Evans, N. W., & Erkal, D. 2018, MNRAS, 476, 2092 [CrossRef] [Google Scholar]
 Long, R. J. 2016, Res. Astron. Astrophys., 16, 189 [CrossRef] [Google Scholar]
 Long, R. J., & Mao, S. 2012, MNRAS, 421, 2580 [CrossRef] [Google Scholar]
 Long, R. J., Mao, S., Shen, J., & Wang, Y. 2013, MNRAS, 428, 3478 [NASA ADS] [CrossRef] [Google Scholar]
 Lowing, B., Jenkins, A., Eke, V., & Frenk, C. 2011, MNRAS, 416, 2697 [CrossRef] [Google Scholar]
 Magorrian, J. 2019, MNRAS, 484, 1166 [CrossRef] [Google Scholar]
 Meiron, Y., Li, B., HolleyBockelmann, K., & Spurzem, R. 2014, ApJ, 792, 98 [CrossRef] [Google Scholar]
 Merritt, D., & Fridman, T. 1996, ApJ, 460, 136 [NASA ADS] [CrossRef] [Google Scholar]
 Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533 [NASA ADS] [Google Scholar]
 Patsis, P. A., Athanassoula, E., & Quillen, A. C. 1997, ApJ, 483, 731 [NASA ADS] [CrossRef] [Google Scholar]
 Patsis, P. A., Skokos, C., & Athanassoula, E. 2002, MNRAS, 337, 578 [NASA ADS] [CrossRef] [Google Scholar]
 Patsis, P. A., Skokos, C., & Athanassoula, E. 2003, MNRAS, 342, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Pearson, S., PriceWhelan, A. M., & Johnston, K. V. 2017, Nat. Astron., 1, 633 [NASA ADS] [CrossRef] [Google Scholar]
 Pfenniger, D. 1984, A&A, 134, 373 [NASA ADS] [Google Scholar]
 Plummer, H. C. 1911, MNRAS, 71, 460 [NASA ADS] [CrossRef] [Google Scholar]
 Portail, M., Wegg, C., Gerhard, O., & Ness, M. 2017, MNRAS, 470, 1233 [NASA ADS] [CrossRef] [Google Scholar]
 PriceWhelan, A. M., Sesar, B., Johnston, K. V., & Rix, H.W. 2016, ApJ, 824, 104 [CrossRef] [Google Scholar]
 Rahmati, A., & Jalali, M. A. 2009, MNRAS, 393, 1459 [CrossRef] [Google Scholar]
 Rix, H., de Zeeuw, P. T., Cretton, N., van der Marel, R. P., & Carollo, C. M. 1997, ApJ, 488, 702 [NASA ADS] [CrossRef] [Google Scholar]
 RojasNiño, A., Read, J. I., Aguilar, L., & Delorme, M. 2016, MNRAS, 459, 3349 [CrossRef] [Google Scholar]
 Schwarzschild, M. 1979, ApJ, 232, 236 [NASA ADS] [CrossRef] [Google Scholar]
 Skokos, C., Patsis, P. A., & Athanassoula, E. 2002, MNRAS, 333, 847 [NASA ADS] [CrossRef] [Google Scholar]
 Syer, D., & Tremaine, S. 1996, MNRAS, 282, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Thomas, J., Saglia, R. P., Bender, R., et al. 2004, MNRAS, 353, 391 [NASA ADS] [CrossRef] [Google Scholar]
 van den Bosch, R. C. E., van de Ven, G., Verolme, E. K., Cappellari, M., & de Zeeuw, P. T. 2008, MNRAS, 385, 647 [NASA ADS] [CrossRef] [Google Scholar]
 van de Ven, G., van den Bosch, R. C. E., Verolme, E. K., & de Zeeuw, P. T. 2006, A&A, 445, 513 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Vasiliev, E. 2013, MNRAS, 434, 3174 [CrossRef] [Google Scholar]
 Vasiliev, E. 2019, MNRAS, 482, 1525 [NASA ADS] [CrossRef] [Google Scholar]
 Vasiliev, E., & Athanassoula, E. 2015, MNRAS, 450, 2842 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, Y., Wu, X., & Zhao, H. 2008, ApJ, 677, 1033 [CrossRef] [Google Scholar]
 Wang, Y., Zhao, H., Mao, S., & Rich, R. M. 2012, MNRAS, 427, 1429 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, Y., Mao, S., Long, R. J., & Shen, J. 2013, MNRAS, 435, 3437 [CrossRef] [Google Scholar]
 Wang, Y., Mao, S., Li, H., et al. 2019, MNRAS, 483, 3048 [CrossRef] [Google Scholar]
 Weinberg, M. D. 1999, AJ, 117, 629 [CrossRef] [Google Scholar]
 Weinberg, M. D., & Katz, N. 2002, ApJ, 580, 627 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, X., Zhao, H., Wang, Y., Llinares, C., & Knebe, A. 2009, MNRAS, 396, 109 [CrossRef] [Google Scholar]
 Wu, X., Wang, Y., Feix, M., & Zhao, H. 2017, ApJ, 844, 130 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, H. S. 1996, MNRAS, 283, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Zheng, Z., Carlin, J. L., Beers, T. C., et al. 2014, ApJ, 785, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, L., Long, R. J., Mao, S., et al. 2014, ApJ, 792, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, Q., Hernquist, L., Marinacci, F., Springel, V., & Li, Y. 2017, MNRAS, 466, 3876 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, L., van den Bosch, R., van de Ven, G., et al. 2018, MNRAS, 473, 3000 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1. Potential and forces for the Hernquist models. From left to right: the input density model is spherical, oblate, prolate, and triaxial, respectively. In each panel, the black line is the result for the model given by Merritt & Fridman (1996). The red, blue, and green lines are the results for the HO method with (n_{max}, l_{max}) = (16,12), (12,12), and (6,4), respectively. From top to bottom: the radius r was measured along the x, x, y, and zaxes, respectively. The xaxis is the major axis and the zaxis is the minor axis of the density distribution. 

Open with DEXTER  
In the text 
Fig. 2. Same as in Fig. 1, but for the results of the CylSP method. In each panel, the red, and blue lines are the results of the CylSP method with R_{min} = 0.1, and 0.01, respectively, while the green lines correspond to the value calculated by the AGAMA software (see Sect. 3.2). The black line is the analytical result given by Merritt & Fridman (1996). Here N_{z} = N_{R} = 40, z_{min} = R_{min}, and R_{max} = z_{max} = 100, m_{max} = 12. 

Open with DEXTER  
In the text 
Fig. 3. Relative error in the potential and forces for both the HO (red lines) and CylSP (blue lines) methods for the analytical Hernquist models. For the HO method, (n_{max}, l_{max}) = (16,12). For the CylSP method, R_{min} is given by the AGAMA software, N_{R} = 40, N_{z} = 40, z_{min} = 0, R_{max} = z_{max} = 100. 

Open with DEXTER  
In the text 
Fig. 4. Relative errors of the potential and forces for models with a different γ value for r = 0.01. From the left to right: the input density model is spherical, oblate, prolate, and triaxial, respectively. The red, blue, and green lines are the results for the HO method with (n_{max}, l_{max}) = (16,12), (12,12), and (6,4), respectively. 

Open with DEXTER  
In the text 
Fig. 5. Same as in Fig. 4, but for the CylSP method. The red, blue, and green lines represent the results with R_{min} = 0.1, 0.01, and the value provided by AGAMA, respectively. 

Open with DEXTER  
In the text 
Fig. 6. Relative error in potential and forces for both the HO (left) and the CylSP methods (right) for the prolate models with a different c/a at r = 0.01. Left panels: the red, blue, and green lines are the results for the HO method with (n_{max}, l_{max}) = (16,12), (12,12), and (6,4), respectively. Right panels: the red, blue, and green lines represent the results for the CylSP method with R_{min} = 0.1, 0.01, and the value given by AGAMA, respectively. From top to bottom: the radius r is along the x, x, y, and zaxes, respectively. 

Open with DEXTER  
In the text 
Fig. 7. Relative errors of the potential and forces for the Hernquist models with different particle numbers at r = 0.01 when using the HO method. In each panel, the red, blue, and green lines are the results for HO method with (n_{max}, l_{max}) = (16,12), (12,12), and (6,4), respectively. The filled circles are the mean value from 40 random samples, and the error bar is the standard error from these samples. 

Open with DEXTER  
In the text 
Fig. 8. Same as in Fig. 7, but for the CylSP method. The red, blue, and green lines represent the results in the CylSP method for R_{min} = 0.1, 0.01, and the value given by AGAMA, respectively. 

Open with DEXTER  
In the text 
Fig. 9. CPU time for the HO method. Left: dependence of the initialization CPU time on the particle number of a realization of the Hernquist model. The red, blue, green, and black lines are the results for the spherical, oblate, prolate, and triaxial models, respectively. The blue, green, and black lines are slightly shifted along the xaxis. Right: dependence of the subsequent CPU time on the number of the expansion coefficients. The subsequent CPU time is estimated by the forces at 10 000 positions. 

Open with DEXTER  
In the text 
Fig. 10. CPU time for the CylSP method. Top left: dependence of the initialization CPU time on the particle number of the Hernquist model realization. Top right: dependence of the initialization CPU time on the parameter N_{R} for an Nbody realization of the Hernquist model with N_{p} = 10^{6}. Bottom left: dependence of the initialization CPU time on the parameter m_{max} for a snapshot of an Nbody realization of the Hernquist model with N_{p} = 10^{6}. Bottom right: dependence of the subsequent CPU time on the parameter m_{max}. The subsequent CPU time is estimated by the forces at 10 000 positions. Excluding the bottom right panel, the blue and black lines are slightly shifted to the right along the xaxis and the red, blue, green, and black lines are the results of the spherical, oblate, prolate, and triaxial models, respectively. 

Open with DEXTER  
In the text 
Fig. 11. Potential and forces for the barred galaxy models and their individual components. From left to right: the input density model is the Ferrers bar, the MN disk, the Plummer halo, and the total, i.e., bar + disk + halo, respectively. In each panel, the black line corresponds to the results for the model given by the analytical method directly. The red, blue, and green lines are the results for the HO method with (n_{max}, l_{max}) = (16,12), (12,12), and (6,4), respectively. 

Open with DEXTER  
In the text 
Fig. 12. Same as in Fig. 11, but the only difference is that the input is an Nbody realization of our barred galaxy model. The total particle number for the bar is 2 × 10^{6}. The particle numbers for the bar, disk, and halo are 200 000, 1640 000, and 160 000, respectively. 

Open with DEXTER  
In the text 
Fig. 13. Same as in Fig. 11, but for the CylSP method. The red, blue, and green lines represent the result in the CylSP method with R_{min} = 0.1, 0.01, and the value provided by AGAMA, respectively. Since the results from three R_{min} are the same as the true values, only the black lines can been seen in most panels. 

Open with DEXTER  
In the text 
Fig. 14. Same as in Fig. 13, but the only difference is that the input is an Nbody realization of our barred galaxy model. The total number of particles is 2 × 10^{6}. 

Open with DEXTER  
In the text 
Fig. 15. Potential and forces for the Ferrers bar with different redistribution numbers of particles. In each panel, the black line is the results for the model given by the analytical method. The red, blue and green lines represent the result in the CylSP method with R_{min} = 0.1, R_{min} = 0.01 and the value provided by AGAMA, respectively. From left to right: the results for Nbody realizations with particle number 2 × 10^{5}, 5 × 10^{5}, 10^{6}, and 10^{7}, respectively. 

Open with DEXTER  
In the text 
Fig. 16. Relative errors of potential and forces for the HO (red lines) and CylSP (blue lines) methods for snapshots of an Nbody realization of a Ferrers bar with a particle number N_{p} = 10^{6}. From left to right: the results for the HO with (n_{max}, l_{max}) = (16,12), (20,20), (30,30) and (40,40), respectively. For the CylSP method, R_{min} = 0.1, N_{R} = 40, N_{z} = 40, z_{min} = 0.1, R_{max} = z_{max} = 100. 

Open with DEXTER  
In the text 
Fig. 17. Surfaces of section (SOS) for the bar model. Left column: results for the analytic potential. Second and third columns: results for the HO method with N_{p} = 2 × 10^{5} and N_{p} = 10^{7}, respectively. Fourth and fifth columns: results for the CylSP method with N_{p} = 2 × 10^{5} and N_{p} = 10^{7}, respectively. From top to bottom: the rows correspond to E_{J} = −0.41, −0.39, −0.37, and −0.35, respectively. The red filled circle in each panel is the position of the x1 orbit. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.