N-body Self-consistent Stellar-halo Modeling of the Fornax Dwarf Galaxy

Dwarf spheroidal (dSph) galaxies associated with the Milky Way (MW) and M31 galaxies are the best probes for studying the properties of dark matter (DM). This is because these galaxies are largely DM-dominated objects with dynamical mass-to-light ratios of 10–1000 (Mateo 1998; Gilmore et al. 2007; McConnachie 2012). Moreover, in the context of bottom-up formation scenarios based on Λ-cold DM (ΛCDM) theory, these galaxies are the building blocks of more massive galaxies, and thus studying their properties and evolution is important to understanding galaxy formation (e.g., Tolstoy et al. 2009).

The hypothesis on the existence of DM was put forward by Zwicky (1933) to explain the virial paradox in the Coma cluster. Later, Babcock (1939) found the growth of the rotation curve of M31 in its outer parts and concluded that a large amount of invisible mass was present in it. The hypothesis about DM was revived by Einasto et al. (1974) and Ostriker et al. (1974) in their studies on the rotation curves of galaxies, which appeared to be mostly flat. The flat rotation curves could be artificially explained without DM by the surface density in the disk being inversely proportional to the distance from the center, as established by Mestel (1963). Nevertheless, the strongest arguments in favor of the existence of DM come from the need to explain the large-scale structure formation given the small amplitudes of perturbations in the cosmic microwave background radiation (CMB).

ΛCDM theory is a concordance model in modern cosmology that gives an excellent description of the CMB (e.g., Planck Collaboration et al. 2014), large-scale structure (e.g., Tegmark et al. 2004), and the accelerating expansion of the universe (e.g., Riess et al. 1998). On the other hand, the observational studies on the galactic and subgalactic scales have turned up several controversial issues that continue to challenge the ΛCDM paradigm. The core–cusp problem is one of the open questions in the ΛCDM picture: the cosmological ΛCDM-based pure DM simulations predict that the dark halos on all mass scales have cusped DM density profiles. On the other hand, studies of the H i gas rotation curves of low surface brightness galaxies and the stellar kinematics of dwarf galaxies testify in favor of shallower or cored density profiles of galaxies' dark halos (e.g., de Blok et al. 2001; Gilmore et al. 2007).

Another way to explain the lack of visible matter was proposed by Milgrom (1983), who hypothesized that there was a characteristic acceleration below which Newton's law of gravitation was invalid. This is one of the varieties of modified gravity (MG) that was subsequently developed, for example, by Bekenstein (2004, 2010) as the tensor-vector scalar theory, or TeVeS. A variety of MG models has been developed to date (see Casas et al. 2017), but the ΛCDM theory is still the leading paradigm of modern cosmology because the MG is not able to explain the alternative formation paths of the large-scale structure of the universe as well as the entire set of observations for galaxies and their clusters. For a review, see, for example, Read et al. (2019), where different star formation histories cause different central DM densities.

Early DM study of the Milky Way dSphs was done by Lin & Faber (1983). More recent research can be classified as those based on Jeans equations, distribution functions (DFs), action-based DFs, and Schwarzschild's orbit-based methods (for a review, see Battaglia et al. 2013). In this work, we focus on the Fornax dSph galaxy, which is well studied in the electromagnetic range, but whose DM halo structure requires detailed modeling. Moreover, this galaxy is suggested to have a core-like DM density profile, which behaves like ρ_DM ∼ r^α , where α = −0.22 estimated by Hayashi et al. (2016). Thus, the galaxy can assess the core−cusp problem in the ΛCDM models.

The non-self-consistent (stars in the field of DM) spherical Jeans equation approaches were employed taking into account velocity anisotropy for Fornax by Gilmore et al. (2007), Peñarrubia et al. (2008), Strigari et al. (2008), Walker et al. (2009b), Salucci et al. (2012), and Read et al. (2019). Axisymmetric models were constructed by Hayashi & Chiba (2012) and Hayashi et al. (2016). The DF-based modeling for the stellar component in the field of a parameterized spherical DM potential for the Fornax dSph galaxy was done by Wu (2007) and Amorisco & Evans (2011). But all of these models are not self-consistent in contrast to the method of Kuijken & Dubinski (1995), who iteratively solved the Poisson equation, taking into account the stellar and DM DFs in the common gravitational potential. Their method is realized in the framework of the GalactICS Software.

Several methods based on the action-based DFs have been developed (for a review see Sanders & Binney 2016; Jeffreson et al. 2017, models for globular clusters). The implementation for the Fornax galaxy was done by Pascale et al. (2018), but their models were spherically symmetric.

Schwarzschild's modeling (Schwarzschild 1979; Richstone & Tremaine 1984) was applied to Fornax, Sculptor, Carina, and Sextans dSph galaxies by Breddels & Helmi (2013). The nonspherical light distribution in the spherical DM field was taken into account by Jardel & Gebhardt (2012). The spherically symmetric Schwarzschild model of the Fornax dSph galaxy was recently done by Kowalczyk et al. (2019). In this work, we perform two kinds of self-consistent stellar-halo modelings of the Fornax dSph galaxy, which are the DF-based one constructed by Kuijken & Dubinski (1995) and the orbit-based one coded within the AGAMA framework (Vasiliev 2018) without assumptions of spherical symmetry. To this end, we utilize the results of the DM structures in the Fornax galaxy analyzed by Hayashi & Chiba (2012) and Hayashi et al. (2016). In order to trace the N-body dynamical evolution of our models and to check their self-consistency and stability, we used the code by Dehnen (2002) named falcON.

In this work, we shall find the set of parameters such as the visible and DM masses, characteristic radii, and the radial behavior of the DM profile near the center and at the periphery. Our goal is to make the modeled dynamical characteristics satisfy the observed ones, such as the velocity dispersion profile. We will perform our investigation in the framework of the standard CDM model. The construction of a feasible evolutionary stable galaxy model is important to find constraints on possible DM candidates as is shown by González-Morales et al. (2017) and Safarzadeh & Spergel (2020).

N-body simulations began long ago; see Klypin & Shandarin (1983). We construct equilibrium systems (coordinates, masses, and velocities of particles) using two methods, and then we use the falcON code in order to follow the N-body evolution of the systems.

While constructing our systems, we adjust the stellar surface density profiles of the models to the observations, which may be approximated either by the King (1962) or Plummer (1911) profiles. The parameters taken for these observational profiles are listed in the Table 1. In order to compare our kinematic characteristics to the data, we use the observational velocity dispersion projected onto the line-of-sight profiles by Hayashi & Chiba (2012) from the data of Walker et al. (2009a).

Table 1. Observed Parameters of the Fornax Galaxy

Σ0	D	rc	rt	rhalf	M
$\left({L}_{\odot }\,{\mathrm{pc}}^{-2}\right)$	(kpc)	(')	(')	(kpc)	(106 M⊙)
15.7 ± 5.1	147	17.6 ± 0.2	69.1 ± 0.4	0.668	20

Note. The central surface brightness (Σ₀) from Irwin & Hatzidimitriou (1995), the distance to the Fornax galaxy (D) from Pietrzyński et al. (2009); the structural parameters (r_c and r_t ) for the King formulas, Equation (1), from Battaglia et al. (2006); the half-light radius (r_half); and the luminous mass (M) of the Fornax galaxy from McConnachie (2012).

Download table as:  ASCIITypeset image

Our nearly self-consistent stellar-halo models of the Fornax galaxy are constructed using the DF-based method of Kuijken & Dubinski (1995) and the orbit-based Schwarzschild's method by Vasiliev (2018). This is a step forward after previous studies of this object, which were based on the more approximate Jeans equations or spherically symmetric DF or action-based approach. The fitting of the velocity dispersion profile to the data was done only for the mock galaxy by Vasiliev & Valluri (2020) and has not been done for a real astrophysical object. We did not perform the fitting procedure but instead constructed many axisymmetric models with different assumptions about the stellar and DM density distributions that satisfy the velocity data quite well. From the set of these models, we pick up one model developed by the GalaxtICS Software NEMO code and the other two models by the Schwarzchild-based code in the AGAMA framework. We report on these models in order to:

• 1.  
  show the best NEMO and AGAMA models (number 1 and 2 in Table 2);
• 2.  
  compare two different approaches and codes (NEMO model 1 is analogous to AGAMA model 3);
• 3.  
  study the process of relaxation to equilibrium for different values (NEMO models are changed more during this process than the AGAMA ones);
• 4.  
  trace the magnitude of change of the previous Jeans velocity dispersion fitting in the field of the axisymmetric prolate DM halo found by Hayashi et al. (2016; AGAMA model 3).

For the density profile of the stellar component of the galaxy, we use the spherically symmetric King (AGAMA model 2 in Table 2) and axially symmetric Plummer (NEMO model 1 and AGAMA model 3) profiles.

Table 2. Modeling Parameters of the Fornax Galaxy

Model	Code	α	Stellar Profile	Mp	bp	q
				$\left({10}^{6}{M}_{\odot }\right)$	(kpc)
1	N	−0.22	Plummer	20.0	0.668	⋯
3	A	−0.22		20.0	0.668	0.66
				Σ0	rc	rt
				$\left({L}_{\odot }\,{{\rm{pc}}}^{-2}\right)$	$\left(^{\prime} \right)$	$\left(^{\prime} \right)$
2	A	0.0	King	11.0	12.4	38.8

Note. Different visible parameters, DM profile exponent α, and the type of the modeling code: "N"—NEMO DF-based by Kuijken & Dubinski (1995); "A"—orbit-based AGAMA by Vasiliev (2018). For model 1, we have obtained the Plummer-like density profile for the equatorial plane; it is the function of the combined gravitational potential, which is axially symmetric, but the equipotentials are not ellipsoids.

Download table as:  ASCIITypeset image

For the King profile, we employ the following formulas:

$$\begin{eqnarray}&&\begin{array}{l}{\rho }_{\mathrm{King}}(r)\\ \,\,=\,\left\{\begin{array}{ll}{\rho }_{b}{\left(\displaystyle \frac{1}{\sqrt{1+{\left(r/{r}_{c}\right)}^{2}}}-\displaystyle \frac{1}{\sqrt{1+{\left({r}_{t}/{r}_{c}\right)}^{2}}}\right)}^{2}, & r\leqslant {r}_{t};\\ 0, & r\gt {r}_{t},\end{array}\right.\end{array}\end{eqnarray} \tag{ 1 }$$

where r is the spherical radial coordinate, and ρ_b , r_c , and r_t are the central stellar density, the core, and tidal radius, respectively. In order to calculate the density parameter ρ_b , we use formulas (B3) in Appendix B.1. We take into account the central surface brightness Σ₀ = 11.0 L_⊙ pc⁻² and the core and tidal radii as r_c = 0.53 kpc, r_t = 1.66 kpc. The value of Σ₀ is slightly shifted inside the error bar with respect to the data of Irwin & Hatzidimitriou (1995), and the values of the radii are not inside their error bars (see Table 1), but the best fit overall has been obtained with these values. As we shall see further from Figure 1, this stellar density profile does not differ a lot from the profile of Irwin & Hatzidimitriou (1995).

Zoom In Zoom Out Reset image size

For the oblate Plummer profile, we use the following function ρ_Plummer(R, z) of the cylindrical coordinates as in Hayashi et al. (2016):

$$\begin{eqnarray}\left\{\begin{array}{ll} & {\rho }_{\mathrm{Plummer}}(R,z)={\rho }_{p}({m}_{\star })=\displaystyle \frac{3{M}_{p}}{4\pi {{b}_{p}}^{3}}{\left[1+\displaystyle \frac{{{m}_{\star }}^{2}}{{{b}_{p}}^{2}}\right]}^{-5/2},\\ & {{m}_{\star }}^{2}={R}^{2}+\displaystyle \frac{{z}^{2}}{{q}^{2}},\end{array}\right.\end{eqnarray} \tag{ 2 }$$

with the half-light radius b_p calculated by Walker et al. (2010) as half the surface brightness of the King profile. The same value for the Plummer profile was used by Hayashi et al. (2016). As for the mass parameter for the Plummer profile M_p , we use the mass of the Fornax galaxy from McConnachie (2012). These values are listed in Table 1. The oblateness parameter q is calculated from the apparent axial ratio q', and the galaxy inclination i is taken from Hayashi et al. (2016; see Table 3):

$$\begin{eqnarray}&&q=\sqrt{{{q}^{{\prime} }}^{2}-{\cos }^{2}i}/\sin i,\end{eqnarray} \tag{ 3 }$$

with the parameters listed in Table 3. The density profile of the DM halo was also taken from the Hayashi et al. (2016). It is a function ρ(R, z) of the cylindrical coordinates, with parameters also listed in Table 3:

$$\begin{eqnarray}\begin{array}{rcl}\rho (R,z) & = & \rho (m)\\ & = & {\rho }_{0}{\left(\displaystyle \frac{m}{{b}_{\mathrm{halo}}}\right)}^{\alpha }{\left[1+{\left(\displaystyle \frac{m}{{b}_{\mathrm{halo}}}\right)}^{2}\right]}^{-(\alpha +3)/2},\\ {m}^{2} & = & {R}^{2}+\displaystyle \frac{{z}^{2}}{{Q}^{2}},\end{array}\end{eqnarray} \tag{ 4 }$$

where b_halo, α, and Q are the scale length, the inner slope, and the axial ratio of the DM density profile, respectively. We also tried the cored DM profile as was done in Hayashi & Chiba (2015) but for the prolate form of the halo (AGAMA model 2).

Table 3. Hydrodynamical Parameters

	α	${\mathrm{log}}_{10}\left({\rho }_{0}\right)$	${\mathrm{log}}_{10}\left({b}_{\mathrm{halo}}\right)$
		$\left({M}_{\odot }\,{{\rm{pc}}}^{-3}\right)$	(pc)
(1)	$-{0.22}_{-0.22}^{+0.14}$	−1.06 ± 0.15	${2.79}_{-0.15}^{+0.16}$
	Q	i	q'
		(deg)
	${1.11}_{-0.53}^{+0.66}$	${71.85}_{-15.34}^{+11.56}$	0.7

Note. Values from Hayashi et al. (2016)—(1).

Download table as:  ASCIITypeset image

Taking into account all these data and varying the parameters, we constructed three models. The diversity of parameters for the stellar component and two variants of the exponent of the DM profile α are listed in Table 2. The model numbered 1 is implemented using the DF-based mkkd95 code (GalactICS Software, NEMO) designed by Kuijken & Dubinski (1995) and the last two ones by the orbit-based Schwarzschild's method implemented in the AGAMA framework designed for constructing equilibrium models by Vasiliev (2018). For all these methods, we used G = 1 as the gravitational constant, 10⁶ M_⊙ as the mass unit, and 1 kpc as the unit of distance. Then, for the time unit we have:

$$\begin{eqnarray}&&[T]=\displaystyle \frac{{{\rm{kpc}}}^{\tfrac{3}{2}}}{{\left(G\times {10}^{6}{M}_{\odot }\right)}^{\tfrac{1}{2}}}=0.47\,{\rm{Gyr}}.\end{eqnarray} \tag{ 5 }$$

Model 1 relies on the shallow DM profile with α = −0.22 and uses a Plummer-like profile for the visible component of the galaxy. The density distribution depends on the combined gravitational potential, so it is axially symmetric and coincides with the Plummer analytical density profile at the equatorial galaxy plane (z = 0), but the isopycnic surfaces are not ellipsoids. The model does not take into account kinematic constraints for the visible part of the galaxy obtained by Hayashi et al. (2016). The parameter β_z stands for such a constraint. This is a velocity anisotropy parameter:

$$\begin{eqnarray}&&{\beta }_{z}=1-\overline{{{v}_{z}}^{2}}/\overline{{{v}_{R}}^{2}}.\end{eqnarray} \tag{ 6 }$$

For the Fornax galaxy, $-{\mathrm{log}}_{10}(1-{\beta }_{z})={0.28}_{-0.11}^{+0.10}$.

The AGAMA orbit-based code has three different parameters to constrain the velocity anisotropy in the solution: β, β_z , and κ. The β parameter is the spherical anisotropy index:

$$\begin{eqnarray}&&\beta =1-{{\sigma }_{\tau }}^{2}/(2\ast {{\sigma }_{r}}^{2}),\end{eqnarray} \tag{ 7 }$$

where σ_τ is the tangential velocity dispersion, ${{\sigma }_{\tau }}^{2}={{\sigma }_{\phi }}^{2}+{{\sigma }_{\theta }}^{2}$, and σ_r is the dispersion of velocity along the spherical radius. Setting

$$\begin{eqnarray}&&\kappa =1.0\,{\rm{i.e.,}}\,{\sigma }_{\phi }={\sigma }_{R},\end{eqnarray} \tag{ 8 }$$

where σ_R is the dispersion of velocity along the cylindrical radius.

Model 2 assumes the visible King profile and the cored DM profile. Its King profile differs slightly from the observed one, but the model has the best reproduction of the velocity dispersion profile (see the Results section). The velocity anisotropy is expressed by the β parameter. Model 3 describes the cusped DM profile and the visible Plummer profile with the β_z and κ parameters for the kinematic constraints.

For all our N-body models, we used 10⁶ stellar particles and 1.5 × 10⁶ DM halo points. For the AGAMA Schwarzschild's orbit-based modeling, we used 25,000 orbits for each component. For the falcON runs, we used all of the default parameters except ${k}_{\max }=6$ (${\tau }_{\max }={(1/2)}^{{k}_{\max }}$) and softening length = 0.1 (for the code notation, see Appendix E).

3.1. DF-based NEMO Modeling

For the first sample of our nearly self-consistent stellar-halo modeling (Table 2), we use the bulge and DM components of the NEMO code developed by Kuijken & Dubinski (1995). For the visible component in this code, we use a bulge with the King's density profile that has its DF described in Kuijken & Dubinski (1995) by the equation:

$$\begin{eqnarray}{f}_{\mathrm{bulge}}(E)=\left\{\begin{array}{ll}{\rho }_{b}{\left(2\pi {{\sigma }_{b}}^{2}\right)}^{-\tfrac{3}{2}}{e}^{\tfrac{{{\rm{\Psi }}}_{0}-{{\rm{\Psi }}}_{c}}{{{\sigma }_{b}}^{2}}}\left({e}^{-\tfrac{E-{{\rm{\Psi }}}_{c}}{{{\sigma }_{b}}^{2}}}-1\right) & E\lt {{\rm{\Psi }}}_{c},\\ 0 & E\geqslant {{\rm{\Psi }}}_{c},\end{array}\right.\end{eqnarray} \tag{ 9 }$$

and the density distribution in a potential Ψ is described by the equation:

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{\mathrm{bulge}}({\rm{\Psi }}) & = & {\rho }_{b}\left[{e}^{\displaystyle \frac{{{\rm{\Psi }}}_{0}-{\rm{\Psi }}}{{{\sigma }_{b}}^{2}}}{\rm{erf}}\left(\displaystyle \frac{\sqrt{{{\rm{\Psi }}}_{c}-{\rm{\Psi }}}}{{\sigma }_{b}}\right)\right.\\ & & -\left.\displaystyle \frac{1}{\sqrt{\pi }}{e}^{\displaystyle \frac{{{\rm{\Psi }}}_{0}-{\rm{\Psi }}}{{{\sigma }_{b}}^{2}}}\left(2\displaystyle \frac{\sqrt{{{\rm{\Psi }}}_{c}-{\rm{\Psi }}}}{{\sigma }_{b}}+\displaystyle \frac{4}{3}\displaystyle \frac{{\left({{\rm{\Psi }}}_{c}-{\rm{\Psi }}\right)}^{\tfrac{3}{2}}}{{{\sigma }_{b}}^{3}}\right)\right].\end{array}\end{eqnarray} \tag{ 10 }$$

This bulge density distribution follows the equipotential surfaces, so it is neither spherical nor ellipsoidal with the given oblateness. It has three parameters: Ψ_c , ρ_b and σ_b .

The DM component construction is based on the lowered Evans distribution from Evans & Collett (1993) also described in Kuijken & Dubinski (1995). The DF for DM is presented as

$$\begin{eqnarray}&&\begin{array}{l}{f}_{\mathrm{halo}}(E,{{L}_{z}}^{2})\\ \,\,=\,\left\{\begin{array}{ll}\left[\left(A{{L}_{z}}^{2}+B\right){e}^{-\tfrac{E}{{{\sigma }_{0}}^{2}}}+C\right]\left[{e}^{-\tfrac{E}{{{\sigma }_{0}}^{2}}}-1\right] & E\lt 0,\\ 0 & E\geqslant 0,\end{array}\right.\end{array}\end{eqnarray} \tag{ 11 }$$

and the density profile is given by

$$\begin{eqnarray}&&\begin{array}{l}{\rho }_{\mathrm{halo}}(R,{\rm{\Psi }})=\displaystyle \frac{1}{2}{\pi }^{\displaystyle \frac{3}{2}}{{\sigma }_{0}}^{3}\left({{AR}}^{2}{{\sigma }_{0}}^{2}+2B\right){\rm{erf}}\displaystyle \frac{\sqrt{-2{\rm{\Psi }}}}{{\sigma }_{0}}{e}^{\displaystyle \frac{-2{\rm{\Psi }}}{{\sigma }_{0}}}\\ \qquad +\,{\left(2\pi \right)}^{\tfrac{3}{2}}{{\sigma }_{0}}^{3}\left(C-B-{{AR}}^{2}{{\sigma }_{0}}^{2}\right){\rm{erf}}\displaystyle \frac{\sqrt{-{\rm{\Psi }}}}{{\sigma }_{0}}{e}^{\displaystyle \frac{-{\rm{\Psi }}}{{{\sigma }_{0}}^{2}}}\\ \qquad +\,\pi \sqrt{-2{\rm{\Psi }}}\left[{{\sigma }_{0}}^{2}\left(3A{{\sigma }_{0}}^{2}{R}^{2}+2B-4C\right)\right.\\ \left.\qquad +\,\displaystyle \frac{4}{3}{\rm{\Psi }}\left(2C-A{{\sigma }_{0}}^{2}{R}^{2}\right)\right].\end{array}\end{eqnarray} \tag{ 12 }$$

Parameters A, B, and C are expressed by the velocity and density scales σ₀ and ρ₁, the halo core radius R_c , and the flattening parameter q as follows (see Kuijken & Dubinski 1994):

$$\begin{eqnarray}&&A=\displaystyle \frac{8(1-{q}^{2})G{{\rho }_{1}}^{2}}{\sqrt{\pi }{q}^{2}{{\sigma }_{0}}^{7}},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&B=\displaystyle \frac{4{{R}_{c}}^{2}G{{\rho }_{1}}^{2}}{\sqrt{\pi }{q}^{2}{{\sigma }_{0}}^{5}},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&C=\displaystyle \frac{(2{q}^{2}-1){\rho }_{1}}{{\left(2\pi \right)}^{\tfrac{3}{2}}{q}^{2}{{\sigma }_{0}}^{3}}.\end{eqnarray} \tag{ 15 }$$

The density scale ρ₁ is replaced in Kuijken & Dubinski (1995) by R_a :

$$\begin{eqnarray}&&{R}_{a}={\left(\displaystyle \frac{3}{2\pi G{\rho }_{1}}\right)}^{\tfrac{1}{2}}{\sigma }_{0}{e}^{\displaystyle \frac{{{\rm{\Psi }}}_{0}}{2{{\sigma }_{0}}^{2}}}.\end{eqnarray} \tag{ 16 }$$

For the DM NEMO model, we need five parameters:

• 1.  
  q—axial ratio, an optional flattening parameter for the potential;
• 2.  
  Ψ₀—central potential;
• 3.  
  ${v}_{0}=\sqrt{2}\times {\sigma }_{0}$, where σ₀ is the central velocity dispersion;
• 4.  
  R_a —the radius at which the halo rotation curve, if continued at its r = 0 slope, would reach the value $\sqrt{2}{\sigma }_{0}$, a scaling radius for the halo;
• 5.  
  ${{r}_{{ck}}}^{2}=\displaystyle \frac{{{R}_{c}}^{2}}{{{R}_{{\rm{K}}}}^{2}}$—a core smoothing parameter—the ratio of the core radius (R_c ) to the derived King radius (R_K). This is the radius at which the gravitational potential has risen by about $2{{\sigma }_{0}}^{2}$ over its central value, provided that the potential well depth is above $2{{\sigma }_{0}}^{2}$.

For the bulge NEMO model, which stands for the visible part of our galaxy, we need three parameters:

• 1.  
  ρ_b —central density;
• 2.  
  Ψ_c —bulge cutoff potential;
• 3.  
  σ_b —bulge central potential.

The names of the parameters for the mkkd95 code can be found in Table 9 in Appendix E.

First, we relied on the density distributions of the components mentioned in the previous section (the density distribution, Equation (1), for the visible part, and Equation (4) for the DM part of the galaxy) to calculate the parameters v₀, R_a , ${{r}_{{ck}}}^{2}$, ρ_b , Ψ_c , and σ_b . With this approach, we failed to construct the NEMO model using the mkkd95 code. The next idea was to fit the density distributions given by Equations (10) and (12) directly to the density distributions given by Equations (2) and (4). To do this, we need an initial guess on the distribution of the potential derived from our density distributions. We have calculated the combined potential distribution Ψ(R, z = 0) in cylindrical coordinates in the equatorial plane z = 0 (see the Appendices C and D).

For the first model, we need to find the central potential Ψ₀ as the sum of the central prolate DM potential and the central stellar potential. To calculate the central DM potential, we use Equation (D30) for the case r = 0, z = 0. As the μ function for this formula, we used Equation (C9) for the prolate cusped Zhao DM potential. For the Plummer oblate central stellar potential, we used Equation (C5).

After that, we perform a fitting procedure for Equation (10) (with a known Ψ₀ and three free parameters, Ψ_c , σ_b , and ρ_b ) to the Plummer density distribution function, Equation (2). To do this, we need to numerically integrate Equation (D32), the same way as Equation (C3), for the growth of the DM potential over its central value. We also need to get the limit of the formula for the stellar radial potential distribution (Equation (C5)) at the equatorial plane of the galaxy. (We use the value x = 10,000 instead of ∞.)

In Figure 1, the fitted function ρ_bulge(Ψ) from Equation (10), shown by the black line, almost coincides with the data curve ρ_Plummer with the M_p = 20.0 parameter from Equation (2), depicted by the blue dash-line. This value of the parameter makes the stellar density profile the same as the one taken by Hayashi et al. (2016). This profile was also assumed when constructing AGAMA model 3. The red dashed–dotted line represents the King profile taken for AGAMA model 2. This profile lies near the Plummer density profile shown by the blue solid line with the stellar mass ∼20.0 × 10⁶ M_⊙ from McConnachie (2012) and the M_p = 14.0 parameter. The stellar profile for model 2 also lies near the King profile given by Irwin & Hatzidimitriou (1995), shown by the red solid line.

The next step is fitting the function given by Equation (12) with three free parameters, σ₀, B, and C, to the DM density distribution function Equation (4). In Figure 2, the fitted function ρ_halo(Ψ, R) from Equation (12), shown by the black line, almost coincides with the DM density profile ρ(m) from Hayashi et al. (2016; see Equation (4)). This curve is shown by the blue dashed line. Then, we get v₀ from σ₀ and the R_a value from those parameters using the expression in Equation (16). To find the ${{r}_{{ck}}}^{2}$ value, we need to calculate R_c using Equations (14) and (15), and R_K. This value is the radius at which the rise of the combined potential over its central value is equal to $2{{\sigma }_{0}}^{2}$.

Zoom In Zoom Out Reset image size

The fitting procedures were done for the density distributions at the equatorial plane of the galaxy. With these two fits, we have obtained five parameters for the mkkd95 NEMO modeling listed in Table 4.

Table 4. Parameters for the NEMO Models

Model	ρb	ψc	σb	q	Ψ0	v0	Ra	${{r}_{{ck}}}^{2}$	Δ
1	20.75	−277.5	7.796	1.11	−494.6	14.97	642.5	1.3 ∗ 1011	11.8%

Download table as:  ASCIITypeset image

The parameters taken from the fitting procedures (Table 4) are the initial parameters for the NEMO modeling. We can also see a variable Δ in Table 4, which is a measure of the violation of the virial theorem at the initial point of the numerical falcON evolution of the constructed galaxy model, namely, Δ = −(2T/W) − 1, where T is the kinetic energy and W < 0 is the potential energy of the system.

3.2. Orbit-based AGAMA Modeling

The AGAMA framework developed by Vasiliev (2018) is able to form an orbit-based model from some specific and common patterns of the densities and DFs.

For the DM halo, we use the spheroid AGAMA component of the following form:

$$\begin{eqnarray}\begin{array}{rcl}\rho (\bar{r}) & = & {\rho }_{0}{\left(\displaystyle \frac{\bar{r}}{{b}_{\mathrm{halo}}}\right)}^{\alpha }{\left[1+{\left(\displaystyle \frac{\bar{r}}{{b}_{\mathrm{halo}}}\right)}^{-\gamma }\right]}^{\tfrac{-\alpha -\eta }{\gamma }}\\ & & \times \exp \left[-{\left(\displaystyle \frac{\bar{r}}{{r}_{\mathrm{cut}}}\right)}^{\xi }\right],\\ \bar{r} & = & \sqrt{{x}^{2}+{\left(y/p\right)}^{2}+{\left(z/Q\right)}^{2}}.\end{array}\end{eqnarray} \tag{ 17 }$$

Far from the r_cut radius, the above form coincides with the DM profile (Equation (4)) when

$$\begin{eqnarray}&&\gamma =2,\,\eta =3,\,p=1.0.\end{eqnarray} \tag{ 18 }$$

For all models, we set η = 3 and γ = 2.0 (see Table 5). As for the r_cut and ξ parameters, we use the following numbers for all our AGAMA models:

$$\begin{eqnarray}&&{r}_{\mathrm{cut}}=55,\,\xi =2.5.\end{eqnarray} \tag{ 19 }$$

Table 5. Parameters for the AGAMA DM Component

Model	ρ0	bhalo	γ	α	Q
2	87.1	0.617	2.0	0.0	1.11
3	87.1	0.617	2.0	−0.22	1.11

Download table as:  ASCIITypeset image

The visible component is set either by the Plummer or King density profiles. To construct the Plummer component using the orbit-based AGAMA code, we need the parameters q, p, b_p , and M_stars for the total component mass parameter (see Table 6). Taking into account Equation (C7) from Appendix C, we can estimate this mass from the M_p parameter:

$$\begin{eqnarray}&&{M}_{\mathrm{stars}}={{qM}}_{p}.\end{eqnarray} \tag{ 20 }$$

For the King component, we need the M_stars, r_cA , and W₀ = [Ψ(r_t ) − Ψ(0)]/σ² parameters (see Table 6). The latter is the dimensionless potential depth of the generalized King (lowered isothermal) models (see the AGAMA documentation). The identifiers of parameters for the AGAMA Schwarzschild's code can be found in Table 10 in Appendix E.

We have gotten two AGAMA models, one of them having the same cusped DM density profile as Hayashi et al. (2016).

The first AGAMA model 2 assumes the visible density King profile with the velocity anisotropy β = −0.17. The characteristic radius and density for the halo component are from Hayashi et al. (2016), but the profile is cored. The second AGAMA model 3 employs the Plummer visible density profile and more appropriate velocity anisotropy parameters for the model of Hayashi et al. (2016) and the same parameters for the DM component.

Figure 3 demonstrates that the orbit-based AGAMA code produces systems that are in much better equilibrium than the ones from the mkkd95 NEMO code. The falcON code computes the virial ratio for each step of system evolution. The values of this ratio for all our models over the time evolution steps are depicted in Figure 3. We can see that the AGAMA values Δ for the violation of the virial theorem (models 2–3) are negligible in comparison to the ones for the NEMO model.

Zoom In Zoom Out Reset image size

From this figure, all of the models reach equilibrium sooner or later.

4.1. NEMO Modeling

The Dehnen (2002) manipulators for the N-body galaxy snapshot files are able to calculate the profiles of different physical quantities averaged either over the spherical radius or over the cylindrical radius. In the latter case, we use the projection along the line of sight. Figure 4 shows the evolution of the modeled surface density profiles of the visible and DM components and the initial analytical profiles. There are several analytical profiles of the surface density along different axes and the averaged surface density.

Zoom In Zoom Out Reset image size

We also found the best-fit parameters of the King distribution for the time t = 0.4 modeled profile. The fitting dashed–dotted line is hidden by the orange, red and brown symbols standing for t = 0.3, t = 0.4 and t = 10 density profiles. The parameters for the line are included in the legend of Figure 4.

For the model, we can see that the stellar surface density profile (in the left panel of Figure 4) decreases by about 20% at the first evolution steps and then fluctuates about its equilibrium position. This decrease is a factor of 2 faster than the decrease of our other NEMO models not presented here. But this model shows the best agreement with the velocity dispersion data among all our NEMO models.

Earlier than t ∼ 10.0, the virial ratio does not reach its equilibrium value (Figure 3) but we can say that all further fluctuations of the projected profiles lie between the t = 0.1 and t = 0.4 lines.

In the right panels of Figure 4, we can see a smaller relative drop of the surface density during the first evolution steps in comparison with the visible density profile.

We can compare the modeled velocity dispersion projected onto the line of sight and averaged over the cylindrical radius of the sky plane profile and its falcON evolution directly with the observational velocity dispersion points. We need to multiply the projprof Dehnen manipulator velocity data by V_scale = 2.079 in order to convert the ${\sqrt{{10}^{6}{M}_{\odot }{\rm{kpc}}}}^{-1}$ velocity units to km s⁻¹.

We use the χ² test to compare our models:

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i=1}^{N}\displaystyle \frac{{\left({\sigma }_{l}\left({{R}_{\mathrm{obs}}}^{i}\right)-\left({{\sigma }_{\mathrm{obs}}}^{i}\right)\right)}^{2}}{{\left({\delta {\sigma }_{\mathrm{obs}}}^{i}\right)}^{2}},\end{eqnarray} \tag{ 21 }$$

where ${\sigma }_{l}({{R}_{\mathrm{obs}}}^{i})$ are the mean modeled velocity dispersions projected onto the line of sight over the modeled star points spaced ${{R}_{\mathrm{obs}}}^{i}$ apart on the line-of-sight axes, ${{\sigma }_{\mathrm{obs}}}^{i}$ are the observational velocity dispersions at ${{R}_{\mathrm{obs}}}^{i}$ radii, and ${\delta {\sigma }_{\mathrm{obs}}}^{i}$ are observational errors.

We do not perform the fitting procedure for minimizing χ² by varying parameters because each set of parameters stands for one 10⁶ particle N-body simulation so the fitting procedure would take too many computational resources. That is why we use hydrodynamical parameters to start our modeling. The mean value of $\overline{{\chi }^{2}}$ averaged over time for t ≥ 0.35 is 83.5 for model 1. In Appendix A.1, we discuss the time interval for averaging and the equilibrium velocity dispersion profile, which comes so much earlier than the Δ parameter does. Figure 5 represents the velocity dispersion projected onto the line-of-sight profile for several evolution time points of falcON in comparison to the observational data with error bars.

Zoom In Zoom Out Reset image size

4.2. AGAMA Modeling

The AGAMA-modeled visible and DM surface density profiles are depicted in Figure 6. We can see the decrease of the stellar surface density profile for the equilibrium-establishing time interval. Such decrease in the central part of the galaxy is less than that for the NEMO model 1.

Zoom In Zoom Out Reset image size

But in contrast to this model, the changing of the density profiles over time for the AGAMA systems is noticeable only at the central part of the galaxy with r ≲ 0.25 kpc. For the NEMO model, the radius of noticeable change of the surface density profile is not less than 0.5 kpc.

For model 2, we found the Σ₀, r_c , and r_t King profile parameters that fit well the starting surface density profile with the initial W₀, r_c , and M_stars AGAMA parameters. They are shown in the panel for this model in the Figure 6.

For model 3, we depicted the averaged Plummer surface density profile by a black line, taking into account our start parameters.

The second row of Figure 6 shows the DM surface density profiles. All of the initial profiles almost coincide with the analytical profiles taking our initial parameters. We can also see small evolution changes in these profiles. The greatest change for model 2 is less than the smallest one for our NEMO model.

Model 3 has almost the same initial parameters as NEMO model 1, which coincides with the Hayashi et al. (2016) parameters. Model 3 changes less than the corresponding NEMO model, but the ensuing evolution curves of model 1 intersect more data error bars than those of the AGAMA model.

Table 6. Parameters for the AGAMA Visible Component

Model	Mstars	rcA	W0	β	βz	κ
King visible component
2	12.13	0.753	1.786	−0.17	⋯	⋯
Plummer visible component
		bp	q
3	13.19	0.668	0.66	⋯	0.47	1.0

Download table as:  ASCIITypeset image

We can see the average values of χ² over time for t ≥ 0.2 for the AGAMA models and for t ≥ 0.35 for the NEMO model in Table 7. In Appendix A.2, the discussion on this averaging can be found. The mean equilibrium χ² value for NEMO model 1 is less than that of AGAMA model 3.

Table 7. Mean Values of χ²

	1	2	3
$\overline{{\chi }^{2}}$	83.5	36.5	104.7

Download table as:  ASCIITypeset image

The worst χ² is for model 3 with α = −0.22 and β_z = 0.47. It seems that this galaxy shows a better cored DM profile. Among the multitude of our other models, we can say that the exact Hayashi et al. (2016) profile is better with the β_z = 0.47 velocity anisotropy parameter and the α = 0.0 or γ = 1.0 DM profiles are better with the β_z = −0.17 or β = −0.17 parameters.

In Figure 7, we can see that the AGAMA models are compliant with the observational data with a different form of curves. The notation of the points is the same as in Figure 5.

Zoom In Zoom Out Reset image size

Let us compare the best NEMO and AGAMA model stellar component equilibrium profiles. The King parameters for them are listed in Table 8. We can see that for both models, Σ₀ is on the bottom edge of the observed error bar or a bit lower. The structural parameters are closer to the Battaglia et al. (2006) data.

Table 8. King Best-fit Parameters vs. Observations

	Σ0, L⊙ pc−2	rc , kpc	rt , kpc
Observed 1	15.7 ± 5.1	0.75 ± 0.01	2.95 ± 0.02
Observed 2	15.7 ± 5.1	0.58 ± 0.05	3.04 ± 0.17
Model 1	9.7	0.69	1.9
model 2	11.0	0.53	1.66

Note. Fitted to the observation structural parameters from Battaglia et al. (2006)—observed 1; Irwin & Hatzidimitriou (1995)—observed 2.

Download table as:  ASCIITypeset image

Let us compare our results with the investigations of Pascale et al. (2018), Amorisco & Evans (2011), and Strigari et al. (2008). We can estimate the dynamical mass M_dyn(r) of our best model inside spheres of 1.05 kpc and 300 pc radii.

For M_dyn(1.05 kpc), we have 1.29 × 10⁸ M_⊙ for our model 2. The estimate of Pascale et al. (2018) for this mass is 1.38 ± 0.10 × 10⁸ M_⊙ compared with the estimate of Amorisco & Evans (2011) at 1.3 × 10⁸ M_⊙. We can see that our numbers are consistent with previous results.

The mass enclosed in a sphere of 300 pc radius is 0.93 × 10⁷ M_⊙ for model 2. The results of Strigari et al. (2008) for M_dyn(0.3 kpc) is ${1.14}_{-0.12}^{+0.09}\times {10}^{7}{M}_{\odot }$. The smaller value of ${0.44}_{-0.03}^{+0.07}\times {10}^{7}{M}_{\odot }$ is the estimate of Pascale et al. (2018). We can see that our numbers for a smaller radius are closer to the older estimates of Strigari et al. (2008).

We can also compare our results to the DM halo density at 150 pc radius modeled by Read et al. (2019) using their GravSphere code based on the spherical Jeans equations and their sophisticated assumptions about the velocity anisotropy parameter. They obtained the value ${\rho }_{\mathrm{DM}}(150\,{\rm{pc}})\,={0.79}_{-0.19}^{+0.27}\times {10}^{8}{M}_{\odot }\,{{\rm{kpc}}}^{-3}$. The values of this density for all our models are inside this error bar. We can say that our results are also consistent with the estimates of Read et al. (2019).

The main point of this paper is that we first investigate whether the dynamical nonspherical structure of the Fornax dSph galaxy estimated from kinematic analysis by Hayashi et al. (2016) is stable or not. We have constructed N-body models using two source codes for the Fornax galaxy from these hydrodynamical studies and followed their numerical evolution.

Then, we found that the AGAMA model with the inferred dark halo parameters demonstrates that, indeed, such a system may be stable during several dynamical times.

Both methods (NEMO and AGAMA) used in our paper give approximate realizations of stationary models. AGAMA is much more accurate, which is clear from the criterion of the virial test. Moreover, AGAMA can reproduce the nonspherical self-consistent structure of Fornax as the weighted sum of the orbit contributions to the galactic density even though the inferred dark halo parameters come from Jeans analysis, which does not require that any distribution functions should be positive.

We have also changed some assumptions about the stellar and DM density distributions and got the model that fits the kinematical data much better than the models found in previous studies. The better DM profiles are cored profiles. The velocity anisotropy parameter for the stellar component also shows some correlations with the DM profile characteristics. The best stellar profile is the King profile. We choose this model in favor of the best agreement with the velocity dispersion data than with the luminosity data because the velocity dispersion observations are more precise than the luminosity ones.

One should understand that even exactly stationary DF-based models may be dynamically or secularly unstable (there are numerous well-known stationary solutions for collisionless disk-like or ellipsoidal systems, which may be unstable). Thus, N-body evolution is needed to test not only the deviations from the exact stationary solution but also to check the stability. Our results provided by falcON simulations do show that all our models are dynamically stable because they support their shapes with small oscillations for several dynamical timescales.

Using N-body simulation results, one can constrain the prior distributions of the DM halo parameters such as scale radius and scale density. Based on observational results (such as the Jeans analysis by Hayashi et al. 2016), our analysis, applied to other dwarf galaxies, can obtain the limits on the parameter ranges of the DM halo that sustains dynamical equilibrium. Therefore, using the constrained parameter ranges, we may improve the limits on the nature of DM particles through their annihilation (Ando et al. 2020).

The authors are grateful to Eugene Vasiliev for his invaluable help with the AGAMA code and to Marat Potashov for debugging the falcOn code.

We are grateful to the anonymous referee for useful comments, which helped us improve the presentation of our results. We are grateful to Tim Blinnikov for his invaluable help with the English text of the paper.

This work was supported by JSPS KAKENHI grant Nos. 18H04359 and 18J00277 to K.H.

S.B. is grateful to Sternberg Astronomical Institute, MSU, for its generous support.

A.1. NEMO Results

In Figure 8, variations of the χ² values over time are depicted. We can see how the χ² moves toward its equilibrium value for NEMO model 1 for times t < 4 in the right panel of Figure 8. This equilibrium value of χ² can be estimated as the mean values for times about t ≥ 0.35 and is depicted by a horizontal black line in Figure 8. The estimates are included in Table 8. We can see that the χ² quantity starts to oscillate over its equilibrium value much earlier than the Δ quantity (see Figure 3), so we choose this time interval for averaging χ² for the NEMO model.

Zoom In Zoom Out Reset image size

A.2. AGAMA Results

Let us consider in detail the equilibrium characteristics of our AGAMA models. Their virial ratio evolution is depicted in Figure 9. We can see that the absolute magnitude of the Δ start values is not greater than their following evolution values in contrast to our NEMO model. But the sign of the violation for all our N-body systems starts from a positive one tending to a negative value during the first time interval of the evolution, see Figure 9. In terms of equilibrium, all our AGAMA models are equally good.

Zoom In Zoom Out Reset image size

Let us look at the evolution of the χ² values for models 2–3 in Figure 8 in order to compare their conformity with the data. In the right panel of Figure 8, we can see the initial steps of the χ² evolution. The setting in of equilibrium is completed by time t = 0.2 (see the brown vertical line) for model 3 and is not resolved for model 2. So in Figure 8, the lines for the mean values of χ² after t = 0.2 for models 2 and 3 are depicted. We find that the best model with the lowest χ² is model 2 with the core DM density profile and the King best-fit stellar component parameters.

In this appendix, we give a summary of the useful formulas for the King (1962) distribution which are not easily found in the textbook.

B.1. Central Surface Density of the King's Distribution

Now we shall express the characteristic density ρ_b in terms of the observed quantities r_c , t, and Σ₀—the central surface density of the Fornax galaxy. Then, we are to find the expression of Σ₀ for the King's distribution:

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{0}\equiv 2{\int }_{0}^{{r}_{t}}{\rho }_{\mathrm{King}}(r){dr}\equiv 2{\rho }_{b}{r}_{c}{s}_{0}(t),\end{eqnarray} \tag{ B1 }$$

where

$$\begin{eqnarray}&&{s}_{0}(t)=\arctan t-\displaystyle \frac{2\mathrm{log}\left(t+\sqrt{{t}^{2}+1}\right)}{\sqrt{{t}^{2}+1}}+\displaystyle \frac{t}{{t}^{2}+1}.\end{eqnarray} \tag{ B2 }$$

For the characteristic density ρ_b expressed in terms of r_c , t, and Σ₀, we get

$$\begin{eqnarray}&&{\rho }_{b}=\displaystyle \frac{{{\rm{\Sigma }}}_{0}}{2{r}_{c}{s}_{0}(t)}.\end{eqnarray} \tag{ B3 }$$

In this appendix, we cite some formulae from Binney & Tremaine (2008) in order to write down equations for the radial distribution of the potential at the equatorial plane for the Plummer oblate density profile and DM oblate density profiles. The latter formulae are useful for the radial distribution of the potential for the prolate DM density profiles. We need these distributions for fitting the NEMO code input parameters.

Let us consider a gravitational system with equidensity axisymmetric ellipsoids of ellipticity e and axial ratio Q. Then, for our parameter, we take m, so that for cylindrical coordinates of an equidensity ellipsoid (R, z), we have the equation:

$$\begin{eqnarray}&&{m}^{2}={R}^{2}+{z}^{2}/{Q}^{2},\,{Q}^{2}=1-{e}^{2}.\end{eqnarray} \tag{ C1 }$$

Then we can define density as a function of m, ρ(m), and define the function μ(m) as

$$\begin{eqnarray}&&\mu (m)={\int }_{0}^{m}\rho ({\bar{m}}^{2})d({\bar{m}}^{2}).\end{eqnarray} \tag{ C2 }$$

Then, from Binney & Tremaine (2008), we have a formula for the gravitational potential Ψ(R₀, z₀) of such systems. Let us write down the increase of the potential over its central value in terms of Q instead of e:

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Psi }}({R}_{0},{z}_{0})-{\rm{\Psi }}(0,0)=\pi {{GQa}}_{0}{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{\mu (m)d\tau }{(\tau +{{a}_{0}}^{2})\sqrt{\tau +{{b}_{0}}^{2}}},\\ {\rm{where}}\ \tau \ {\rm{is}}\,{\rm{defined}}\,{\rm{by}}\,{\rm{the}}\,{\rm{equation:}}\\ \displaystyle \frac{{{R}_{0}}^{2}}{\tau +{{a}_{0}}^{2}}+\displaystyle \frac{{{z}_{0}}^{2}}{\tau +{{b}_{0}}^{2}}=\displaystyle \frac{{m}^{2}}{{{a}_{0}}^{2}},\,{b}_{0}={{Qa}}_{0}.\end{array}\end{eqnarray} \tag{ C3 }$$

C.1. Plummer Potential Distribution at the Equatorial Plane

Let us find the potential of the prolate Plummer distribution, Equation (2):

$$\begin{eqnarray*}&&\left\{\begin{array}{l}{\rho }_{\mathrm{Plummer}}(R,z)={\rho }_{p}({m}_{\star })=\displaystyle \frac{3{M}_{p}}{4\pi {{b}_{p}}^{3}}{\left[1+\displaystyle \frac{{{m}_{\star }}^{2}}{{{b}_{p}}^{2}}\right]}^{-5/2},\\ {{m}_{\star }}^{2}={R}^{2}+\displaystyle \frac{{z}^{2}}{{q}^{2}},\end{array}\right.\end{eqnarray*}$$

with the oblateness q and ellipticity $e=\sqrt{1-{q}^{2}}$. Then, the function μ from Equation (C2) for this profile will be

$$\begin{eqnarray}&&\mu (m)=\displaystyle \frac{{M}_{p}}{2\pi {b}_{p}}\left(1-\displaystyle \frac{1}{{\left(1+{m}^{2}/{{b}_{p}}^{2}\right)}^{3/2}}\right).\end{eqnarray} \tag{ C4 }$$

Let us find the potential distribution at the equatorial plane for z₀ = 0:

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Psi }}({R}_{0},0.0) & = & -\displaystyle \frac{{M}_{p}{Gq}}{{b}_{p}\sqrt{1-{q}^{2}}}\left(\displaystyle \frac{x}{{\left(k+1\right)}^{3/2}}{F}_{1}\left(\displaystyle \frac{1}{2};-\displaystyle \frac{1}{2};\displaystyle \frac{3}{2};\displaystyle \frac{3}{2};\right.\right.\\ & & -{\left.\left.\left.{x}^{2};-\displaystyle \frac{{x}^{2}}{k+1}\right)\right)\right|}_{q/e}^{\infty },\end{array}\end{eqnarray} \tag{ C5 }$$

where F₁ is a Gauss hypergeometric function.

C.2. Plummer Mass for the Oblate System

In order to calculate the mass of such a system, we shall use formulae for the mass of a thin homeoid between ellipsoids with oblateness q and equatorial radii m and m + dm from Binney & Tremaine (2008):

$$\begin{eqnarray}&&\delta M=4\pi \rho ({m}_{\star }){{m}_{\star }}^{2}{qdm}.\end{eqnarray} \tag{ C6 }$$

Then, the mass inside such an ellipsoid will be equal to

$$\begin{eqnarray}&&M({m}_{\star })=4\pi q{\int }_{0}^{{m}_{\star }}\rho (\tilde{m}){\tilde{m}}^{2}d\tilde{m}={{qM}}_{p}\displaystyle \frac{{m}^{3}}{{\left({{b}_{p}}^{2}+{m}^{2}\right)}^{3/2}}.\end{eqnarray} \tag{ C7 }$$

C.3. Zhao Density Distribution with One Parameter

Let us consider the Hayashi et al. (2016) density distribution of the DM halo. We can find a more common form of such a DM profile in Zhao (1996). We shall rewrite the function ρ(m) from Equation (4) in terms of the variable p = m/b_halo:

$$\begin{eqnarray}&&\rho (p)={\rho }_{0}{p}^{\alpha }{\left(1+{p}^{2}\right)}^{-(\alpha +3)/2},\,p=\displaystyle \frac{m}{{b}_{\mathrm{halo}}},\,\tilde{p}=\displaystyle \frac{\bar{m}}{{b}_{\mathrm{halo}}},\end{eqnarray} \tag{ C8 }$$

and calculate the function μ(m) for our density distribution:

$$\begin{eqnarray}\begin{array}{rcl}\mu (m) & = & \displaystyle \frac{{\rho }_{0}{b}_{\mathrm{halo}}^{2}}{\tfrac{\alpha }{2}+1}{\tilde{p}}^{\alpha +2}2{F}_{1}\left(\displaystyle \frac{\alpha }{2}+1,\displaystyle \frac{\alpha }{2}+1.5,\right.\\ & & {\left.\left.\displaystyle \frac{\alpha }{2}+2,-{\tilde{p}}^{2}\right)\right|}_{0}^{\displaystyle \frac{m}{{b}_{\mathrm{halo}}}},\end{array}\end{eqnarray} \tag{ C9 }$$

where 2F₁(a, b, c, z) is the Gauss hypergeometric function.

Let us derive the analog of formulas (C3) for prolate systems, i.e., equidensity axisymmetric ellipsoids with axial ratio Q > 1. First, we shall examine a thin prolate ellipsoidal shell and then summarize the contribution of such shells to the total potential. By definition, the axial ratio of the shell is Q, and it has a coordinate m, denoting the shell coordinates R_⋆, z_⋆ as in Equation (C1):

$$\begin{eqnarray}&&{{R}_{\star }}^{2}+\displaystyle \frac{{{z}_{\star }}^{2}}{{Q}^{2}}={m}^{2}.\end{eqnarray} \tag{ D1 }$$

Let us define the prolate coordinates as in problem 2.5 of Chapter 2 in Binney & Tremaine (2008),

$$\begin{eqnarray}&&R=a\sinh u\sin v,\,z=a\cosh u\cos v,\,\phi {\rm{\unicode{x02014}azimuthal}}\,{\rm{angle}}.\end{eqnarray} \tag{ D2 }$$

For comparison, the oblate coordinates for Q < 1 are $R\,=a\cosh u\sin v,z=a\sinh u\cos v$, ϕ—the same azimuthal angle. For our coordinates, we can express $\sin v$ and $\cos v$ from Equation (D2) and use the trigonometrical identical equation:

$$\begin{eqnarray}&&\displaystyle \frac{{R}^{2}}{{\sinh }^{2}u}+\displaystyle \frac{{z}^{2}}{{\cosh }^{2}u}={a}^{2}.\end{eqnarray} \tag{ D3 }$$

We are defining our prolate coordinates so, that one of the ellipsoids with u = constant coincides with our thin shell that generates a gravitational potential δΨ. Let us denote the coordinate for this shell as u_⋆. Then, comparing Equations (D2) and (D3) for R_⋆, z_⋆, and u_⋆, we can write for a and u_⋆:

$$\begin{eqnarray}&&\coth {u}_{\star }=Q,\,a=m\,\csc \,{u}_{\star }=m\sqrt{{Q}^{2}-1}.\end{eqnarray} \tag{ D4 }$$

All our ellipsoids with u = constant are confocal with the focus $a=m\sqrt{{Q}^{2}-1}$. We can see that our coordinates are suitable for Q > 1, because $\coth {u}_{\star }\gt 1$ for all positive u_⋆. And for oblate coordinates it will be $\tanh {u}_{\star }=Q\lt 1$ for all positive u_⋆.

In order to find the potential δΨ, we are to solve the equation ∇² δΨ = 0. We shall find the expression for ∇² δΨ in our prolate coordinates. Now we are varying in turn the coordinates u, v, and ϕ and moving along the orthogonal vectors ${\hat{{\boldsymbol{e}}}}_{u},{\hat{{\boldsymbol{e}}}}_{v},{\hat{{\boldsymbol{e}}}}_{\phi }$ by the distances h_u δ u, h_v δ v, h_ϕ δ ϕ:

$$\begin{eqnarray}\begin{array}{rcl}{h}_{u} & = & a\sqrt{{\cosh }^{2}u{\sin }^{2}v+{\sinh }^{2}u{\cos }^{2}v},\\ {h}_{v} & = & a\sqrt{{\cosh }^{2}u{\sin }^{2}v+{\sinh }^{2}u{\cos }^{2}v},\\ {h}_{\phi } & = & R=a\sinh u\sin v.\end{array}\end{eqnarray} \tag{ D5 }$$

After some transformations, we have for h_u , h_v :

$$\begin{eqnarray}&&{h}_{u}={h}_{v}=a\sqrt{{\cosh }^{2}u-{\cos }^{2}v}.\end{eqnarray} \tag{ D6 }$$

For comparison, we can write down the coefficients for oblate systems in oblate coordinates: ${h}_{u}={h}_{v}=a\sqrt{{\sinh }^{2}u+{\cos }^{2}v}$, ${h}_{\phi }=a\cosh u\sin v$. Following the equation for the gradient:

$$\begin{eqnarray}&&{\rm{\nabla }}=\displaystyle \frac{{\hat{{\boldsymbol{e}}}}_{i}}{{h}_{i}}\displaystyle \frac{\partial }{\partial {q}_{i}},\,{q}_{1}=u,\,{q}_{2}=v,\,{q}_{3}=\phi ,\end{eqnarray} \tag{ D7 }$$

we shall write the gradient of the potential δΨ:

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\nabla }}\left(\delta {\rm{\Psi }}\right) & = & \displaystyle \frac{1}{a\sqrt{{\cosh }^{2}u-{\cos }^{2}v}}\left(\displaystyle \frac{\partial \left(\delta {\rm{\Psi }}\right)}{\partial u}{\hat{{\boldsymbol{e}}}}_{u}+\displaystyle \frac{\partial \left(\delta {\rm{\Psi }}\right)}{\partial v}{\hat{{\boldsymbol{e}}}}_{v}\right)\\ & & +\displaystyle \frac{1}{a\sinh u\sin v}\displaystyle \frac{\partial \left(\delta {\rm{\Psi }}\right)}{\partial \phi }{\hat{{\boldsymbol{e}}}}_{\phi }.\end{array}\end{eqnarray} \tag{ D8 }$$

The equation for the Laplacian is

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\nabla }}}^{2}F & = & \displaystyle \frac{1}{{h}_{1}{h}_{2}{h}_{3}}\left(\displaystyle \frac{\partial }{\partial {q}_{1}}\left(\displaystyle \frac{{h}_{2}{h}_{3}}{{h}_{1}}\displaystyle \frac{\partial F}{\partial {q}_{1}}\right)\right.\\ & & \left.+\displaystyle \frac{\partial }{\partial {q}_{2}}\left(\displaystyle \frac{{h}_{3}{h}_{1}}{{h}_{2}}\displaystyle \frac{\partial F}{\partial {q}_{2}}\right)+\displaystyle \frac{\partial }{\partial {q}_{3}}\left(\displaystyle \frac{{h}_{1}{h}_{2}}{{h}_{3}}\displaystyle \frac{\partial F}{\partial {q}_{3}}\right)\right).\end{array}\end{eqnarray} \tag{ D9 }$$

For h₁ = h₂ and h₃, all independent of q₃ we can rewrite this equation in the following way:

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\nabla }}}^{2}F & = & \displaystyle \frac{1}{{{h}_{1}}^{2}{h}_{3}}\left(\displaystyle \frac{\partial }{\partial {q}_{1}}\left({h}_{3}\displaystyle \frac{\partial F}{\partial {q}_{1}}\right)+\displaystyle \frac{\partial }{\partial {q}_{2}}\left({h}_{3}\displaystyle \frac{\partial F}{\partial {q}_{2}}\right)\right)\\ & & +\displaystyle \frac{1}{{{h}_{3}}^{2}}\displaystyle \frac{{\partial }^{2}F}{{\partial }^{2}{q}_{3}}.\end{array}\end{eqnarray} \tag{ D10 }$$

Then, for the Laplacian of the function δΨ we can write down:

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\nabla }}}^{2}\left(\delta {\rm{\Psi }}\right) & = & \displaystyle \frac{1}{{a}^{2}({\sinh }^{2}u-{\cos }^{2}v)}\left(\displaystyle \frac{1}{\sinh u}\displaystyle \frac{\partial }{\partial u}\left(\sinh u\displaystyle \frac{\partial \left(\delta {\rm{\Psi }}\right)}{\partial u}\right)\right.\\ & & \left.+\displaystyle \frac{1}{\sin v}\displaystyle \frac{\partial }{\partial v}\left(\sin v\displaystyle \frac{\partial \left(\delta {\rm{\Psi }}\right)}{\partial v}\right)\right)\\ & & +\displaystyle \frac{1}{{a}^{2}{\sinh }^{2}u{\sin }^{2}v}\displaystyle \frac{{\partial }^{2}\left(\delta {\rm{\Psi }}\right)}{{\partial }^{2}\phi }.\end{array}\end{eqnarray} \tag{ D11 }$$

We are interested in the potentials, constant over the ϕ and v variables, so we can write:

$$\begin{eqnarray}&&\begin{array}{l}\left\{\begin{array}{l}\displaystyle \frac{\partial \left(\delta {\rm{\Psi }}\right)}{\partial v}=0,\\ \displaystyle \frac{\partial \left(\delta {\rm{\Psi }}\right)}{\partial \phi }=0,\\ {{\rm{\nabla }}}^{2}\left(\delta {\rm{\Psi }}\right)=0.\end{array}\right.\Rightarrow \quad \displaystyle \frac{d}{{du}}\left(\sinh u\displaystyle \frac{d\left(\delta {\rm{\Psi }}\right)}{{du}}\right)=0\\ \quad \iff \quad \displaystyle \frac{d\left(\delta {\rm{\Psi }}\right)}{{du}}=\displaystyle \frac{A}{\sinh u},\,\delta {\rm{\Psi }}=A\mathrm{log}\left(\tanh \displaystyle \frac{u}{2}\right)+B.\end{array}\end{eqnarray} \tag{ D12 }$$

In order to find the constants A and B, we are to take the limit of the function δΨ at the infinite value of variables R and z, i.e., infinite $r=\sqrt{{R}^{2}+{z}^{2}}$. The limit of the potential of our ellipsoidal shell should be

$$\begin{eqnarray}&&{\mathrm{lim}}_{r\to \infty }\delta {\rm{\Psi }}=-\displaystyle \frac{G\delta M}{r},\end{eqnarray} \tag{ D13 }$$

where δ M is the mass of the ellipsoidal shell.

Let us express r in terms of the u and v variables:

$$\begin{eqnarray}&&{r}^{2}\equiv {R}^{2}+{z}^{2}={a}^{2}({\sinh }^{2}u+{\cos }^{2}v),\end{eqnarray} \tag{ D14 }$$

and find $\sinh u$ from this expression:

$$\begin{eqnarray}&&{a}^{2}{\sinh }^{2}u={r}^{2}-{a}^{2}{\cos }^{2}v.\end{eqnarray} \tag{ D15 }$$

Now, subtracting Equation (D15) from ${a}^{2}{\cosh }^{2}u$, we can use the identical hypergeometric equation:

$$\begin{eqnarray}&&{a}^{2}{\cosh }^{2}u={r}^{2}+{a}^{2}{\sin }^{2}v.\end{eqnarray} \tag{ D16 }$$

Then we shall rewrite $\tanh \displaystyle \frac{u}{2}$ in terms of r and v, using Equation (D16):

$$\begin{eqnarray}&&\mathrm{log}\left(\tanh \displaystyle \frac{u}{2}\right)=\displaystyle \frac{1}{2}\mathrm{log}\left(\displaystyle \frac{\sqrt{{r}^{2}+{a}^{2}\sin v}-a}{\sqrt{{r}^{2}+{a}^{2}\sin v}+a}\right).\end{eqnarray} \tag{ D17 }$$

Now we shall take a Taylor series expansion of δΨ for small $a/r$:

$$\begin{eqnarray}&&\delta {\rm{\Psi }}=-A\displaystyle \frac{a}{r}+B+o\left(\displaystyle \frac{a}{r}\right).\end{eqnarray} \tag{ D18 }$$

To get the limit of δΨ the same as in Equation (D13), we can write down for constants A and B:

$$\begin{eqnarray}&&A=\displaystyle \frac{G\delta M}{a},B=0.\end{eqnarray} \tag{ D19 }$$

For the points inside our shell, the potential is constant and equal to its value on the shell itself. Let us do some hypertrigonometric transformations to get this value:

$$\begin{eqnarray}&&\mathrm{log}\left(\tanh \displaystyle \frac{{u}_{\star }}{2}\right)=\mathrm{log}\left(Q-\sqrt{{Q}^{2}-1}\right).\end{eqnarray} \tag{ D20 }$$

So, the final result for the potential of a thin prolate ellipsoidal shell is

$$\begin{eqnarray}&&\delta {\rm{\Psi }}=\left\{\begin{array}{l}\displaystyle \frac{G\delta M}{m\sqrt{{Q}^{2}-1}}\mathrm{log}\left(\tanh \displaystyle \frac{u}{2}\right),\quad u\gt {u}_{\star },\mathrm{outside};\\ \displaystyle \frac{G\delta M}{m\sqrt{{Q}^{2}-1}}\mathrm{log}\left(Q-\sqrt{{Q}^{2}-1}\right),\quad u\leqslant {u}_{\star },\mathrm{inside}.\end{array}\right.\end{eqnarray} \tag{ D21 }$$

Let us find the mass δ M of our shell with the density ρ(m). The semiaxes of the ellipsoid are m and Qm, so the total volume inside the ellipsoid is

$$\begin{eqnarray}&&V=\displaystyle \frac{4}{3}\pi {m}^{3}Q,\end{eqnarray} \tag{ D22 }$$

and the volume δ V inside the shell is

$$\begin{eqnarray}&&\delta V2\pi {Qm}\delta \left({m}^{2}\right).\end{eqnarray} \tag{ D23 }$$

Then, the mass δ m of the shell is

$$\begin{eqnarray}&&\delta M=2\pi Q\rho (m)m\delta \left({m}^{2}\right).\end{eqnarray} \tag{ D24 }$$

Now we can rewrite the Equation (D21) for the potential of a thin prolate ellipsoidal shell using the expression for δ M:

$$\begin{eqnarray}&&\begin{array}{l}\delta {\rm{\Psi }}\\ \,\,=\,\left\{\begin{array}{l}2\pi G\rho (m)\displaystyle \frac{Q}{\sqrt{{Q}^{2}-1}}\mathrm{log}\left(\tanh \displaystyle \frac{u}{2}\right)\delta \left({m}^{2}\right),\quad u\gt {u}_{\star };\\ 2\pi G\rho (m)\displaystyle \frac{Q}{\sqrt{{Q}^{2}-1}}\mathrm{log}\left(Q-\sqrt{{Q}^{2}-1}\right)\delta \left({m}^{2}\right),\quad u\leqslant {u}_{\star }\end{array}\right.,\end{array}\end{eqnarray} \tag{ D25 }$$

where u_⋆ is the coordinate of the shell and u is the coordinate of the ellipsoid, which is confocal with the thin shell and passes through the point (R₀, z₀).

D.1. Gravitational Potential for Prolate Systems

Now we are ready to summarize the potentials of prolate ellipsoidal shells at the point (R₀, z₀) to get the total potential Ψ. For each shell, we have the coordinates m, R_⋆, z_⋆, u_⋆ of the shell itself, which satisfy the Equation (D3), and the coordinate of the ellipsoid u_m on which the point R₀, z₀ lies, and this ellipsoid is confocal with the shell. So, the coordinates u_m , R₀, z₀ also satisfy Equation (D3):

$$\begin{eqnarray}&&\displaystyle \frac{{{R}_{0}}^{2}}{{\sinh }^{2}{u}_{m}}+\displaystyle \frac{{{z}_{0}}^{2}}{{\cosh }^{2}{u}_{m}}={m}^{2}\left({Q}^{2}-1\right).\end{eqnarray} \tag{ D26 }$$

Our point (R₀, z₀) lies on the equidensity shell with the coordinate m₀. Then, the point (R₀, z₀) lies outside the shell for coordinate m < m₀ and inside the shell for coordinate m > m₀. Let us denote by δΨ_out(m) the contribution to the total potential Ψ of those shells, for which the point (R₀, z₀) lies outside the shell. Then, for such shells, our coordinate u_m is greater then the shell coordinate u_⋆: u_m > u_⋆. The contribution of the shells, for which the point (R₀, z₀) lies inside the shell, will be denoted by δΨ_in(m). Then, for these shells, our coordinate u_m is less than the shell coordinate u_⋆: u_m < u_⋆. The sum of δΨ_in(m) is the case u ≤ u_⋆ in Equation (D25), and the sum of δΨ_out(m) is the case u > u_⋆ in that formula. We shall write down the sum of δΨ_in(m), taking into account the definition of the function μ(m) (Equation (C2)):

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \sum _{m\gt {m}_{0}}\delta {{\rm{\Psi }}}_{\mathrm{in}}(m) & = & 2\pi G\displaystyle \frac{Q}{\sqrt{{Q}^{2}-1}}\left(\mu (\infty )-\mu ({m}_{0})\right)\\ & & \times \mathrm{log}\left(Q-\sqrt{{Q}^{2}-1}\right).\end{array}\end{eqnarray} \tag{ D27 }$$

The sum of δΨ_out(m) will be

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \sum _{m\lt {m}_{0}}\delta {{\rm{\Psi }}}_{\mathrm{out}}(m) & = & 2\pi G\displaystyle \frac{Q}{\sqrt{{Q}^{2}-1}}\left(\mu (m)\mathrm{log}{\left.\left(\tanh \displaystyle \frac{{u}_{m}}{2}\right)\right|}_{0}^{{m}_{0}}\right.\\ & & -\left.{\displaystyle \int }_{0}^{{m}_{0}}\mu (m)d\left(\mathrm{log}\left(\tanh \displaystyle \frac{{u}_{m}}{2}\right)\right)\right).\end{array}\end{eqnarray} \tag{ D28 }$$

Let us rewrite the differential:

$$\begin{eqnarray}&&d\left(\mathrm{log}\left(\tanh \displaystyle \frac{{u}_{m}}{2}\right)\right)=\displaystyle \frac{d\cosh {u}_{m}}{{\sinh }^{2}{u}_{m}}.\end{eqnarray} \tag{ D29 }$$

For m = m₀, the ellipsoid u_m coincides with the ellipsoid u_⋆, and we can write $\coth {u}_{m}=Q$ as in Equation (D4). So, we can evaluate $\cosh {u}_{m}=\tfrac{Q}{\sqrt{{Q}^{2}-1}}$ for m = m₀. For m = 0 to satisfy Equation (D3), we set $\coth {u}_{m}=\infty$. Summarizing Equations (D27) and (D28) and taking into account the differential Equation (D29), we can derive the gravitational potential at the point Ψ(R₀, z₀):

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Psi }}({R}_{0},{z}_{0}) & = & 2\pi G\displaystyle \frac{Q}{\sqrt{{Q}^{2}-1}}\left(\mu (\infty )\mathrm{log}\left(Q-\sqrt{{Q}^{2}-1}\right)\right.\\ & & +\left.{\displaystyle \int }_{\cosh {u}_{m}=\displaystyle \frac{Q}{\sqrt{{Q}^{2}-1}}}^{\infty }\mu (m)\displaystyle \frac{d\cosh {u}_{m}}{{\sinh }^{2}{u}_{m}}\right).\end{array}\end{eqnarray} \tag{ D30 }$$

Let us rewrite the expressions for $\sinh {u}_{m},\cosh {u}_{m}$ in terms of τ, a₀, and b₀:

$$\begin{eqnarray}\begin{array}{rcl}{b}_{0} & = & {{Qa}}_{0};\,\tau +{{a}_{0}}^{2}={{a}_{0}}^{2}{\sinh }^{2}{u}_{m}\left({Q}^{2}-1\right)\\ & & \Rightarrow \,{{a}_{0}}^{2}{\cosh }^{2}{u}_{m}\left({Q}^{2}-1\right)=\tau +{{b}_{0}}^{2},\\ & \quad & d\tau =2{a}_{0}\sqrt{{Q}^{2}-1}\sqrt{\tau +{{b}_{0}}^{2}}d\cosh {u}_{m}.\end{array}\end{eqnarray} \tag{ D31 }$$

For m = m₀, we can write $\sinh {u}_{m}=\tfrac{1}{\sqrt{{Q}^{2}-1}}$ and then τ = 0. For m = 0 and $\cosh {u}_{m}=\infty$, we can write τ = ∞ . Then, the rise of the potential over its central value could be written as

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Psi }}({R}_{0},{z}_{0})-{\rm{\Psi }}(0,0) & = & \pi {{GQa}}_{0}{\int }_{0}^{\infty }\mu (m)\\ & & \times \displaystyle \frac{d\tau }{(\tau +{{a}_{0}}^{2})\sqrt{\tau +{{b}_{0}}^{2}}}.\end{array}\end{eqnarray} \tag{ D32 }$$

We can see that the definitions of a₀, b₀ and Equations (C3) and (D32) for Ψ(R₀, z₀) − Ψ(0, 0) are identical for oblate and prolate systems!

We can see the vocabulary for the program code parameters in Tables 9 and 10. The first row of the tables are the code parameters in bold font and the second row are the Greek letter notations for the parameters used in that paper. Table 9 shows the NEMO mkkd95 code parameters. Table 10 shows the AGAMA Schwarzschild's modeling code parameters. The notation for the falcON code parameters are as follows:

• 1.  
  kmax—k_max,
• 2.  
  eps—.

Table 9. Code Parameters for the NEMO mkkd95 Runs

q	psi0	v0	ra	rck2	rhob	psicut	sigb
q	Ψ0	v0	Ra	${{r}_{{ck}}}^{2}$	ρb	Ψc	σb

Download table as:  ASCIITypeset image