EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 639 (July 2020) A&A, 639 (2020) A38 Full HTML
Free Access
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/0004-6361/202038225
Published online 03 July 2020

© 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 made-to-measure 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. Price-Whelan 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 N-body simulations. We note that N-body 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 N-body simulations is that the potential models are self-consistent; 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 density-potential pairs using basis-function expansion techniques (e.g. Rojas-Niño et al. 2016). The Poisson’s equation can be solved in spherical coordinates (e.g. Clutton-Brock 1973; Hernquist & Ostriker 1992; Zhao 1996; Rahmati & Jalali 2009; Lowing et al. 2011; Lilley et al. 2018) or in cylindrical coordinates (e.g. Clutton-Brock 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; Price-Whelan et al. 2016; Pearson et al. 2017) and for N-body 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 spherical-harmonic 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 basis-function in spherical coordinates

ρ ( r , θ , ϕ ) = ∑ l = 0 l max ∑ m = 0 l ∑ n = 0 n max A nlm ρ nl Y lm ( θ , ϕ ) , $$ \begin{aligned}&\rho (r,\theta ,\phi ) = \sum _{l=0}^{l_{\rm max}}\sum _{m=0}^{l}\sum _{n=0}^{n_{\rm max}}A_{nlm}\rho _{nl}Y_{lm}(\theta ,\phi ), \end{aligned} $$(1)

Φ ( r , θ , ϕ ) = ∑ l = 0 l max ∑ m = 0 l ∑ n = 0 n max B nlm Φ nl Y lm ( θ , ϕ ) , $$ \begin{aligned}&\Phi (r,\theta ,\phi ) = \sum _{l=0}^{l_{\rm max}}\sum _{m=0}^{l}\sum _{n=0}^{n_{\rm max}}B_{nlm}\Phi _{nl}Y_{lm}(\theta ,\phi ), \end{aligned} $$(2)

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. Clutton-Brock 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 three-dimensional 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

ρ ( R , z , ϕ ) = ∑ m = 0 m max ρ m ( R , z ) exp ( i m ϕ ) , $$ \begin{aligned}&\rho (R,z,\phi ) = \sum _{m=0}^{m_{\rm max}}\rho _m(R,z)\exp (im\phi ), \end{aligned} $$(3)

Φ ( R , z , ϕ ) = ∑ m = 0 m max Φ m ( R , z ) exp ( i m ϕ ) , $$ \begin{aligned}&\Phi (R,z,\phi ) = \sum _{m=0}^{m_{\rm max}}\Phi _m(R,z)\exp (im\phi ), \end{aligned} $$(4)

and each term in the expansion is given by

ρ m ( R , z ) = 1 2 π ∫ 0 2 π d ϕ ρ ( R , z , ϕ ) exp ( − i m ϕ ) , $$ \begin{aligned}&\rho _m(R,z) = \frac{1}{2\pi }\int _0^{2\pi }\mathrm{d}\phi \rho (R,z,\phi )\exp (-im\phi ), \end{aligned} $$(5)

Φ m ( R , z ) = − G ∫ − ∞ + ∞ ∫ 0 ∞ d z ′ d R ′ 2 π R ′ ρ m ( R ′ , z ′ ) Ξ . $$ \begin{aligned}&\Phi _m(R,z) = -G\int _{-\infty }^{+\infty }\int _{0}^{\infty }\mathrm{d}z^{\prime }\mathrm{d}R^{\prime }2\pi R^{\prime }\rho _m(R^{\prime },z^{\prime })\,\Xi . \end{aligned} $$(6)

and

Ξ ( R , z , R ′ , z ′ ) = ∫ 0 ∞ d k J m ( k R ) J m ( k R ′ exp ( − k | z − z ′ | ) , $$ \begin{aligned} \Xi (R,z,R^{\prime },z^{\prime }) = \int _0^{\infty }\mathrm{d}k J_m(kR)J_m(kR^{\prime }\exp \left(-k \left|z-z^{\prime }\right|\right), \end{aligned} $$(7)

where mmax is the truncated number of azimuthal Fourier harmonic terms. For an axisymmetric model, only the m = 0 term is needed, while mmax ≲ 12 is enough for typical asymmetric models (Vasiliev 2019). We note that Jm is an mth-order Bessel function of the first kind, and can be expressed analytically by the half-integer degree Legendre function of the second kind Qm − 1/2 (Gradshteyn & Ryzhik 1965; Cohl & Tohline 1999; Holley-Bockelmann et al. 2005)

Ξ = 1 π R R ′ Q m − 1 / 2 ( R 2 + R ′ 2 + ( z − z ′ ) 2 2 R R ′ ) , if R , R ′ > 0 , $$ \begin{aligned}&\Xi =\frac{1}{\pi \sqrt{RR^{\prime }}}Q_{m-1/2}\bigg (\frac{R^2+{R^{\prime }}^2+(z-z^{\prime })^2}{2RR^{\prime }} \bigg ),\;\mathrm{if}\ R, R^{\prime }>0, \end{aligned} $$(8)

Ξ = [ R 2 + R ′ 2 + ( z − z ′ ) 2 ] − 1 / 2 δ m 0 , if R = 0 or R ′ = 0 . $$ \begin{aligned}&\Xi =[R^2+{R^{\prime }}^2+(z-z^{\prime })^2]^{-1/2}\delta _{m0}, \; \mathrm{if}\,R=0\, \mathrm{or}\, R^{\prime }=0. \end{aligned} $$(9)

Equation (6) shows that the potential can be obtained by a two-dimensional 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 (NR × Nz) within a fixed region R ∈ [Rmin, Rmax], z ∈ [zmin, zmax], 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 N-body 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)

ρ ( m ) = ( 3 − γ ) M 4 a b c m − γ ( 1 + m ) − ( 4 − γ ) , 0 ≤ γ < 3 , $$ \begin{aligned} \rho (m) = \frac{(3-\gamma )M}{4abc}m^{-\gamma }(1+m)^{-(4-\gamma )}, \quad 0\le \gamma < 3, \end{aligned} $$(10)

with

m 2 = x 2 a 2 + y 2 b 2 + z 2 c 2 , a ≥ b ≥ c . $$ \begin{aligned} m^2=\frac{x^2}{a^2}+\frac{{ y}^2}{b^2}+\frac{z^2}{c^2}, \quad a\ge b\ge c. \end{aligned} $$(11)

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 iso-density 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 one-dimensional 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 x-axis, is the major axis while the z-axis 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 nmax and lmax give more accurate expansions than smaller ones for most models and in most regions. For a spherical system, we find that large values of nmax and lmax generate larger errors in the force along the minor axis (az) 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 high-order terms in the expansion. Also, for the triaxial system, if nmax = 16 and lmax = 12, then the HO method can give accurate expansions.

thumbnail 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 (nmax, lmax) = (16,12), (12,12), and (6,4), respectively. From top to bottom: the radius r was measured along the x-, x-, y-, and z-axes, respectively. The x-axis is the major axis and the z-axis is the minor axis of the density distribution.

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 Rmin in this method. In the AGAMA software implementation (Vasiliev 2019) of the CylSP method that we are using, the code automatically chooses the Rmin 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 Rmin that is given by the software can give accurate values in all cases except in the central-most region (r <  2 × 10−3) for the spherical, prolate, and triaxial models. For the forces from Rmin = 0.01 and Rmin = 0.1, obvious wiggles appear in the central region (r <  0.1).

thumbnail 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 Rmin = 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 Nz = NR = 40, zmin = Rmin, and Rmax = zmax = 100, mmax = 12.

We also consider different values of mmax and NR, Nz for the CylSP method. It should be noted that NR = 40 can give good accuracy as shown in Vasiliev & Athanassoula (2015). We then searched for the optimum mmax value and found that there is hardly any improvement beyond mmax = 12, while the necessary CPU time is considerably increased. For mmax ≥ 12, the accuracy for the force and potential is almost unchanged even when mmax is further increased together with CPU time. Therefore, we adopt NR = Nz = 40 and mmax = 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.

thumbnail 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, (nmax, lmax) = (16,12). For the CylSP method, Rmin is given by the AGAMA software, NR = 40, Nz = 40, zmin = 0, Rmax = zmax = 100.

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 (nmax, lmax) = (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 (nmax, lmax) = (16,12), and we found that these errors between two values for γ that are close can differ quite markedly.

thumbnail 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 (nmax, lmax) = (16,12), (12,12), and (6,4), respectively.

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 Rmin = 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 Rmin for the Dehnen models, while Rmin = 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 Rmin provided by AGAMA and 0.1, while this difference is small for Rmin = 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 Rmin, Rmax and NR, 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 Rmin, Rmax, and NR 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 Rmin = 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 Rmin = 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 nmax = 6, lmax = 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.

thumbnail Fig. 5.

Same as in Fig. 4, but for the CylSP method. The red, blue, and green lines represent the results with Rmin = 0.1, 0.01, and the value provided by AGAMA, respectively.

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.

thumbnail 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 (nmax, lmax) = (16,12), (12,12), and (6,4), respectively. Right panels: the red, blue, and green lines represent the results for the CylSP method with Rmin = 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 z-axes, respectively.

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, ax and ay. 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 106, 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 x-axis.

thumbnail 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 (nmax, lmax) = (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.

It is seen that the forces obtained from the HO method deviate from the true value less than 2% if Np ≥ 5 × 105, and this deviation is nearly independent of the number of expansion coefficients nmax and lmax. For the oblate model, the force deviations for the HO method with (nmax, lmax) = (16,12) are smaller than ∼9% and these deviations weakly depend on the particle number. It is also noted that the forces along the z-axis with (nmax, lmax) = (12,12) are more accurate than those with (nmax, lmax) = (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 (nmax, lmax) = (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 (nmax, lmax) = (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 Rmin provided by the AGAMA software is higher than that for Rmin = 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.

thumbnail 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 Rmin = 0.1, 0.01, and the value given by AGAMA, respectively.

We see that the force accuracy for the CylSP method with Rmin 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 Np ≥ 5 × 105. 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 three-dimensional integration while the expansion coefficients in the CylSP method are obtained from a two-dimensional 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 Np ≥ 106.

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 three-dimensional integrations in the HO method, while those are calculated by two-dimensional 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 N-body 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 N EC = 1 2 ( n max + 1 ) ( l max + 2 ) ( l max + 1 ) $ N_{\mathrm{EC}}=\frac{1}{2}(n_{\mathrm{max}}+1)(l_{\mathrm{max}}+2)(l_{\mathrm{max}}+1) $. For both the analytical model and N-body particles, once the expansion coefficients have been obtained, the subsequent CPU time is the same if the same NEC 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.

thumbnail 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 x-axis. 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.

Figure 10 shows the dependence of the CPU time on the parameters of particle number Np (top left), and NR (top right) and mmax (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 NP and the grid number NR. 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.

thumbnail 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 NR for an N-body realization of the Hernquist model with Np = 106. Bottom left: dependence of the initialization CPU time on the parameter mmax for a snapshot of an N-body realization of the Hernquist model with Np = 106. Bottom right: dependence of the subsequent CPU time on the parameter mmax. 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 x-axis and the red, blue, green, and black lines are the results of the spherical, oblate, prolate, and triaxial models, respectively.

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 NCE terms of the basis function; whereas for the CylSP method, it is only the time to use a two-dimensional 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 two-thirds 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.

thumbnail 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 (nmax, lmax) = (16,12), (12,12), and (6,4), respectively.

Figure 12 also shows the potential and forces for the models using the HO method, but now the input is an N-body 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 N-body 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).

thumbnail Fig. 12.

Same as in Fig. 11, but the only difference is that the input is an N-body realization of our barred galaxy model. The total particle number for the bar is 2 × 106. The particle numbers for the bar, disk, and halo are 200 000, 1640 000, and 160 000, respectively.

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 Rmin are nearly same; the reason is that the potential and accelerations vary slowly in the central region in these models.

thumbnail 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 Rmin = 0.1, 0.01, and the value provided by AGAMA, respectively. Since the results from three Rmin are the same as the true values, only the black lines can been seen in most panels.

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 Np = 106. 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 Rmin = 0.01 and the Rmin provided by AGAMA give larger wiggles than those from Rmin = 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 N-body 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 107.

thumbnail Fig. 14.

Same as in Fig. 13, but the only difference is that the input is an N-body realization of our barred galaxy model. The total number of particles is 2 × 106.

thumbnail 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 Rmin = 0.1, Rmin = 0.01 and the value provided by AGAMA, respectively. From left to right: the results for N-body realizations with particle number 2 × 105, 5 × 105, 106, and 107, respectively.

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 nmax = 16, lmax = 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 nmax = 40, lmax = 40 are used. However, the accuracy in the central region is not improved for the HO method.

thumbnail Fig. 16.

Relative errors of potential and forces for the HO (red lines) and CylSP (blue lines) methods for snapshots of an N-body realization of a Ferrers bar with a particle number Np = 106. From left to right: the results for the HO with (nmax, lmax) = (16,12), (20,20), (30,30) and (40,40), respectively. For the CylSP method, Rmin = 0.1, NR = 40, Nz = 40, zmin = 0.1, Rmax = zmax = 100.

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 y-axis to be the longest major axis. Therefore, we consider the intersection of the orbit with the plane y = 0 and in particular those with vy >  0. For each energy, we generated 25 initial conditions of orbits with x ≠ 0 and vy ≠ 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 EJ from three methods are similar, but the differences in the details are also obvious.

thumbnail 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 Np = 2 × 105 and Np = 107, respectively. Fourth and fifth columns: results for the CylSP method with Np = 2 × 105 and Np = 107, respectively. From top to bottom: the rows correspond to EJ = −0.41, −0.39, −0.37, and −0.35, respectively. The red filled circle in each panel is the position of the x1 orbit.

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 107, 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 N-body bar model with Np = 107 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 EJ. We found that the stability of the x1 orbit at these EJ 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 EJ = −0.37 and EJ = −0.35, which could show considerable differences between the three methods if the disk and halo were also represented by particles.

Table 1.

Four Jacobi energies EJ 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 two-dimensional 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 power-law 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 N-body 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 Np ≥ 5 × 105. If the particle number is smaller than 5 × 105, 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 N-body 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 × 105. 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 × 106.

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 N-body realization. If we consider the bar component of this N-body realization separately, the CylSP method can also generate some wiggles in the central region for the force if a small Rmin 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 Rmin as provided by the AGAMA software is preferred. For models with snapshots of an N-body 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 N-body 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 N-body simulation, a relatively large Rmin (Rmin = 0.1) is suggested.

Acknowledgments

We thank E. Vasiliev and Jean-Charles 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 N-body 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- 2019-05). The “Centre de Calcul Intensif d’Aix-Marseille” 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, co-founded by National Astronomical Observatories, Chinese Academy of Sciences and Alibaba Cloud.

References

  1. Athanassoula, E. 1992a, MNRAS, 259, 328 [NASA ADS] [CrossRef] [Google Scholar]
  2. Athanassoula, E. 1992b, MNRAS, 259, 345 [NASA ADS] [CrossRef] [Google Scholar]
  3. Bovy, J., Kawata, D., & Hunt, J. A. S. 2018, MNRAS, 473, 2288 [CrossRef] [Google Scholar]
  4. 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]
  5. Bryan, S. E., Mao, S., Kay, S. T., et al. 2012, MNRAS, 422, 1863 [NASA ADS] [CrossRef] [Google Scholar]
  6. Buta, R. J., Sheth, K., Regan, M., et al. 2010, ApJS, 190, 147 [NASA ADS] [CrossRef] [Google Scholar]
  7. Buta, R. J., Sheth, K., Athanassoula, E., et al. 2015, ApJS, 217, 32 [NASA ADS] [CrossRef] [Google Scholar]
  8. Carpintero, D. D., & Wachlin, F. C. 2006, Celest. Mech. Dyn. Astron., 96, 129 [CrossRef] [Google Scholar]
  9. Clutton-Brock, M. 1972, Ap&SS, 16, 101 [CrossRef] [Google Scholar]
  10. Clutton-Brock, M. 1973, Ap&SS, 23, 55 [NASA ADS] [CrossRef] [Google Scholar]
  11. Cohl, H. S., & Tohline, J. E. 1999, ApJ, 527, 86 [NASA ADS] [CrossRef] [Google Scholar]
  12. de Lorenzi, F., Debattista, V. P., Gerhard, O., & Sambhus, N. 2007, MNRAS, 376, 71 [NASA ADS] [CrossRef] [Google Scholar]
  13. de Lorenzi, F., Gerhard, O., Coccato, L., et al. 2009, MNRAS, 395, 76 [NASA ADS] [CrossRef] [Google Scholar]
  14. Dehnen, W. 1993, MNRAS, 265, 250 [NASA ADS] [CrossRef] [Google Scholar]
  15. Dehnen, W., & Binney, J. J. 1998a, MNRAS, 298, 387 [NASA ADS] [CrossRef] [Google Scholar]
  16. Dehnen, W., & Binney, J. 1998b, MNRAS, 294, 429 [NASA ADS] [CrossRef] [Google Scholar]
  17. Earn, D. J. D. 1996, ApJ, 465, 91 [NASA ADS] [CrossRef] [Google Scholar]
  18. Erkal, D., Boubert, D., Gualandris, A., Evans, N. W., & Antonini, F. 2019, MNRAS, 483, 2007 [CrossRef] [Google Scholar]
  19. Ferrers, N. M. 1877, Quart. J. Pure Appl. Math., 14, 1 [Google Scholar]
  20. Fux, R. 2001, A&A, 373, 511 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Gradshteyn, I. S., & Ryzhik, I. M. 1965, Table of Integrals, Series and Products (New York: Academic Press) [Google Scholar]
  22. Hernquist, L. 1990, ApJ, 356, 359 [NASA ADS] [CrossRef] [Google Scholar]
  23. Hernquist, L., & Ostriker, J. P. 1992, ApJ, 386, 375 [NASA ADS] [CrossRef] [Google Scholar]
  24. Holley-Bockelmann, K., Weinberg, M., & Katz, N. 2005, MNRAS, 363, 991 [NASA ADS] [CrossRef] [Google Scholar]
  25. Hunt, J. A. S., Kawata, D., & Martel, H. 2013, MNRAS, 432, 3062 [NASA ADS] [CrossRef] [Google Scholar]
  26. Jin, Y., Zhu, L., Long, R. J., et al. 2019, MNRAS, 486, 4753 [CrossRef] [Google Scholar]
  27. Kenyon, S. J., Bromley, B. C., Brown, W. R., & Geller, M. J. 2014, ApJ, 793, 122 [NASA ADS] [CrossRef] [Google Scholar]
  28. Kuijken, K., & Dubinski, J. 1995, MNRAS, 277, 1341 [CrossRef] [Google Scholar]
  29. Lee, G.-H., Park, C., Lee, M. G., & Choi, Y.-Y. 2012, ApJ, 745, 125 [NASA ADS] [CrossRef] [Google Scholar]
  30. Lilley, E. J., Sanders, J. L., Evans, N. W., & Erkal, D. 2018, MNRAS, 476, 2092 [CrossRef] [Google Scholar]
  31. Long, R. J. 2016, Res. Astron. Astrophys., 16, 189 [CrossRef] [Google Scholar]
  32. Long, R. J., & Mao, S. 2012, MNRAS, 421, 2580 [CrossRef] [Google Scholar]
  33. Long, R. J., Mao, S., Shen, J., & Wang, Y. 2013, MNRAS, 428, 3478 [NASA ADS] [CrossRef] [Google Scholar]
  34. Lowing, B., Jenkins, A., Eke, V., & Frenk, C. 2011, MNRAS, 416, 2697 [CrossRef] [Google Scholar]
  35. Magorrian, J. 2019, MNRAS, 484, 1166 [CrossRef] [Google Scholar]
  36. Meiron, Y., Li, B., Holley-Bockelmann, K., & Spurzem, R. 2014, ApJ, 792, 98 [CrossRef] [Google Scholar]
  37. Merritt, D., & Fridman, T. 1996, ApJ, 460, 136 [NASA ADS] [CrossRef] [Google Scholar]
  38. Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533 [NASA ADS] [Google Scholar]
  39. Patsis, P. A., Athanassoula, E., & Quillen, A. C. 1997, ApJ, 483, 731 [NASA ADS] [CrossRef] [Google Scholar]
  40. Patsis, P. A., Skokos, C., & Athanassoula, E. 2002, MNRAS, 337, 578 [NASA ADS] [CrossRef] [Google Scholar]
  41. Patsis, P. A., Skokos, C., & Athanassoula, E. 2003, MNRAS, 342, 69 [NASA ADS] [CrossRef] [Google Scholar]
  42. Pearson, S., Price-Whelan, A. M., & Johnston, K. V. 2017, Nat. Astron., 1, 633 [Google Scholar]
  43. Pfenniger, D. 1984, A&A, 134, 373 [NASA ADS] [Google Scholar]
  44. Plummer, H. C. 1911, MNRAS, 71, 460 [CrossRef] [Google Scholar]
  45. Portail, M., Wegg, C., Gerhard, O., & Ness, M. 2017, MNRAS, 470, 1233 [NASA ADS] [CrossRef] [Google Scholar]
  46. Price-Whelan, A. M., Sesar, B., Johnston, K. V., & Rix, H.-W. 2016, ApJ, 824, 104 [CrossRef] [Google Scholar]
  47. Rahmati, A., & Jalali, M. A. 2009, MNRAS, 393, 1459 [CrossRef] [Google Scholar]
  48. 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]
  49. Rojas-Niño, A., Read, J. I., Aguilar, L., & Delorme, M. 2016, MNRAS, 459, 3349 [CrossRef] [Google Scholar]
  50. Schwarzschild, M. 1979, ApJ, 232, 236 [NASA ADS] [CrossRef] [Google Scholar]
  51. Skokos, C., Patsis, P. A., & Athanassoula, E. 2002, MNRAS, 333, 847 [NASA ADS] [CrossRef] [Google Scholar]
  52. Syer, D., & Tremaine, S. 1996, MNRAS, 282, 223 [NASA ADS] [CrossRef] [Google Scholar]
  53. Thomas, J., Saglia, R. P., Bender, R., et al. 2004, MNRAS, 353, 391 [NASA ADS] [CrossRef] [Google Scholar]
  54. 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]
  55. 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]
  56. Vasiliev, E. 2013, MNRAS, 434, 3174 [CrossRef] [Google Scholar]
  57. Vasiliev, E. 2019, MNRAS, 482, 1525 [NASA ADS] [CrossRef] [Google Scholar]
  58. Vasiliev, E., & Athanassoula, E. 2015, MNRAS, 450, 2842 [NASA ADS] [CrossRef] [Google Scholar]
  59. Wang, Y., Wu, X., & Zhao, H. 2008, ApJ, 677, 1033 [CrossRef] [Google Scholar]
  60. Wang, Y., Zhao, H., Mao, S., & Rich, R. M. 2012, MNRAS, 427, 1429 [NASA ADS] [CrossRef] [Google Scholar]
  61. Wang, Y., Mao, S., Long, R. J., & Shen, J. 2013, MNRAS, 435, 3437 [CrossRef] [Google Scholar]
  62. Wang, Y., Mao, S., Li, H., et al. 2019, MNRAS, 483, 3048 [CrossRef] [Google Scholar]
  63. Weinberg, M. D. 1999, AJ, 117, 629 [CrossRef] [Google Scholar]
  64. Weinberg, M. D., & Katz, N. 2002, ApJ, 580, 627 [NASA ADS] [CrossRef] [Google Scholar]
  65. Wu, X., Zhao, H., Wang, Y., Llinares, C., & Knebe, A. 2009, MNRAS, 396, 109 [CrossRef] [Google Scholar]
  66. Wu, X., Wang, Y., Feix, M., & Zhao, H. 2017, ApJ, 844, 130 [NASA ADS] [CrossRef] [Google Scholar]
  67. Zhao, H. S. 1996, MNRAS, 283, 149 [NASA ADS] [CrossRef] [Google Scholar]
  68. Zheng, Z., Carlin, J. L., Beers, T. C., et al. 2014, ApJ, 785, L23 [NASA ADS] [CrossRef] [Google Scholar]
  69. Zhu, L., Long, R. J., Mao, S., et al. 2014, ApJ, 792, 59 [NASA ADS] [CrossRef] [Google Scholar]
  70. Zhu, Q., Hernquist, L., Marinacci, F., Springel, V., & Li, Y. 2017, MNRAS, 466, 3876 [NASA ADS] [CrossRef] [Google Scholar]
  71. Zhu, L., van den Bosch, R., van de Ven, G., et al. 2018, MNRAS, 473, 3000 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1.

Four Jacobi energies EJ in the bar potential.

In the text

All Figures

thumbnail 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 (nmax, lmax) = (16,12), (12,12), and (6,4), respectively. From top to bottom: the radius r was measured along the x-, x-, y-, and z-axes, respectively. The x-axis is the major axis and the z-axis is the minor axis of the density distribution.
In the text
thumbnail 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 Rmin = 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 Nz = NR = 40, zmin = Rmin, and Rmax = zmax = 100, mmax = 12.
In the text
thumbnail 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, (nmax, lmax) = (16,12). For the CylSP method, Rmin is given by the AGAMA software, NR = 40, Nz = 40, zmin = 0, Rmax = zmax = 100.
In the text
thumbnail 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 (nmax, lmax) = (16,12), (12,12), and (6,4), respectively.
In the text
thumbnail Fig. 5.

Same as in Fig. 4, but for the CylSP method. The red, blue, and green lines represent the results with Rmin = 0.1, 0.01, and the value provided by AGAMA, respectively.
In the text
thumbnail 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 (nmax, lmax) = (16,12), (12,12), and (6,4), respectively. Right panels: the red, blue, and green lines represent the results for the CylSP method with Rmin = 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 z-axes, respectively.
In the text
thumbnail 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 (nmax, lmax) = (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.
In the text
thumbnail 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 Rmin = 0.1, 0.01, and the value given by AGAMA, respectively.
In the text
thumbnail 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 x-axis. 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.
In the text
thumbnail 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 NR for an N-body realization of the Hernquist model with Np = 106. Bottom left: dependence of the initialization CPU time on the parameter mmax for a snapshot of an N-body realization of the Hernquist model with Np = 106. Bottom right: dependence of the subsequent CPU time on the parameter mmax. 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 x-axis and the red, blue, green, and black lines are the results of the spherical, oblate, prolate, and triaxial models, respectively.
In the text
thumbnail 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 (nmax, lmax) = (16,12), (12,12), and (6,4), respectively.
In the text
thumbnail Fig. 12.

Same as in Fig. 11, but the only difference is that the input is an N-body realization of our barred galaxy model. The total particle number for the bar is 2 × 106. The particle numbers for the bar, disk, and halo are 200 000, 1640 000, and 160 000, respectively.
In the text
thumbnail 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 Rmin = 0.1, 0.01, and the value provided by AGAMA, respectively. Since the results from three Rmin are the same as the true values, only the black lines can been seen in most panels.
In the text
thumbnail Fig. 14.

Same as in Fig. 13, but the only difference is that the input is an N-body realization of our barred galaxy model. The total number of particles is 2 × 106.
In the text
thumbnail 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 Rmin = 0.1, Rmin = 0.01 and the value provided by AGAMA, respectively. From left to right: the results for N-body realizations with particle number 2 × 105, 5 × 105, 106, and 107, respectively.
In the text
thumbnail Fig. 16.

Relative errors of potential and forces for the HO (red lines) and CylSP (blue lines) methods for snapshots of an N-body realization of a Ferrers bar with a particle number Np = 106. From left to right: the results for the HO with (nmax, lmax) = (16,12), (20,20), (30,30) and (40,40), respectively. For the CylSP method, Rmin = 0.1, NR = 40, Nz = 40, zmin = 0.1, Rmax = zmax = 100.
In the text
thumbnail 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 Np = 2 × 105 and Np = 107, respectively. Fourth and fifth columns: results for the CylSP method with Np = 2 × 105 and Np = 107, respectively. From top to bottom: the rows correspond to EJ = −0.41, −0.39, −0.37, and −0.35, respectively. The red filled circle in each panel is the position of the x1 orbit.
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.