Lopsidedness of Self-consistent Galaxies Caused by the External Field Effect of Clusters

The MONDian Poisson's equation including the presence of an external gravitational acceleration ${{\boldsymbol{g}}}_{\mathrm{ext}}$ reads (Wu et al. 2007)

$$\begin{eqnarray}&&{\boldsymbol{\nabla }}\cdot \left[\mu \left(\displaystyle \frac{| {\boldsymbol{g}}| }{{a}_{0}}\right){\boldsymbol{g}}\right]=4\pi G{\rho }_{b},\quad {\boldsymbol{g}}={\boldsymbol{\nabla }}{{\rm{\Phi }}}_{\mathrm{int}}+{{\boldsymbol{g}}}_{\mathrm{ext}}.\end{eqnarray} \tag{ 1 }$$

Here ${{\rm{\Phi }}}_{\mathrm{int}}$ is the internal potential generated by the baryon density ${\rho }_{b}$, and ${a}_{0}\approx 3700\,{\mathrm{km}}^{2}\,{{\rm{s}}}^{-2}\,{\mathrm{kpc}}^{-1}$ is Milgrom's constant. To produce both a Newtonian and a MONDian limit, the interpolating function $\mu (x)$ has to be of the following form $(x={\boldsymbol{g}}/{a}_{0})$:

$$\begin{eqnarray}&&\mu (x)\sim x\qquad x\ll 1,\\ &&\mu (x)\sim 1\qquad x\gg 1.\end{eqnarray} \tag{ 2 }$$

In what follows, we will use the simple form of the μ function adopted by Famaey & Binney (2005), which is

$$\begin{eqnarray}&&\mu (x)=\displaystyle \frac{x}{1+x}.\end{eqnarray} \tag{ 3 }$$

Although observationally excluded in the solar system, this μ function is still a viable choice on galactic scales and is in good agreement with the terminal velocities of the Milky Way and NGC3198. Compared to the "standard" form introduced in Milgrom (1983a), the transition between Newtonian gravity and MOND happens more gradually. For the baryonic density, we adopt a Hernquist profile (Hernquist 1990):

$$\begin{eqnarray}&&\rho (r)=\displaystyle \frac{M}{2\pi {abc}}\displaystyle \frac{1}{r{(1+r)}^{3}},\end{eqnarray} \tag{ 4 }$$

where

$$\begin{eqnarray}&&r=\sqrt{{\left(\displaystyle \frac{x}{a}\right)}^{2}+{\left(\displaystyle \frac{y}{b}\right)}^{2}+{\left(\displaystyle \frac{z}{c}\right)}^{2}},\end{eqnarray} \tag{ 5 }$$

and the constants $a,b,c$ are the typical length scales of the galaxy's major, intermediate, and minor axes, respectively. Here we consider five different galaxy models: four axisymmetric elliptical galaxies without and with external fields, respectively, and a triaxial galaxy model embedded in an external field. The parameters of the galaxy models are listed in Table 1. These models represent medium-sized elliptical galaxies with masses on the order of 10¹⁰–10¹¹ ${M}_{\odot }$, bright enough to observe the outer parts. In this case, the internal accelerations are comparable to several a₀. Because elliptical galaxies should lie on the fundamental plane (Djorgovski & Davis 1987), there is a strong correlation between stellar masses M and effective radii R_e (Figure 13 of Gadotti 2009). The effective radii of galaxies with a total mass of 10¹⁰–10¹¹ ${M}_{\odot }$ range from 0.5 to 3 kpc. For all five models, we set $a=1\,\mathrm{kpc}$, which lies within the R_e dispersion range of observational data points (Gadotti 2009).

Table 1.  Galaxy Models and Schwarzschild Parameters: The Total Mass M of Each Model Is Expressed in Units of ${10}^{10}\,{M}_{\odot }$; a, b, and c are Scale Lengths in Units of kpc; and ${g}_{\mathrm{ext}}$ Denotes the Strength of the External Field in Units of ${a}_{0}\approx 3700{(\mathrm{km}{{\rm{s}}}^{-1})}^{2}\,\mathrm{kpc}$

Model	1	2	3	4	5
$M[{10}^{10}\,{M}_{\odot }]$	5	5	5	5	10
$a:b:c$	1: 1: 0.7	1: 1: 0.7	1: 1: 0.7	1: 1: 0.7	1: 0.86: 0.7
${g}_{\mathrm{ext}}$ [a0]	0	0.01	0.1	1.0	1.0
Direction of ${{\boldsymbol{g}}}_{\mathrm{ext}}$	⋯	positive x–z	positive x–z	negative x–z	negative
		diagonal	diagonal	diagonal	x-axis
Density symmetry	axisym	axisym	axisym	axisym	triaxial
Potential sym. axes	x, y, z	y	y	y	y, z
Starting octants	y	I, II, V, VI	I, II, V, VI	I, II, V, VI	I, II
Reflecting planes	x–y, x–z, y–z	x–z	x–z	x–z	x–y, x–z
${N}_{\mathrm{cell}}$	960	3840	3840	3840	1920
${N}_{\mathrm{stationary}}$	3840	15360	15360	15360	7680
${N}_{\mathrm{ejecting}}$	3000	12000	12000	12000	6000
${N}_{\mathrm{orbit}}$	6840	27360	27360	27360	13680

Note. Because the external field breaks the symmetry of the potential, the symmetry axes of potential and density differ from each other. The starting octants refer to the symmetry plane of the potential. Here ${N}_{\mathrm{cell}}$ is the number of cells to impose self-consistency, ${N}_{\mathrm{stationary}}$ is the number of stationary starting orbits built into the orbital library, and ${N}_{\mathrm{ejecting}}$ is the number of ejecting orbits starting from the x–z plane. The total number of orbits is given by ${N}_{\mathrm{orbit}}$.

Download table as:  ASCIITypeset image

The strengths of the external fields are chosen as ${g}_{\mathrm{ext}}\,=0.01{a}_{0},\,0.1{a}_{0}$, and $1.0{a}_{0}$ for the axisymmetric models and $1.0{a}_{0}$ for the triaxial model, corresponding to weak, intermediate, and strong external fields, respectively. Hence for the strong external field cases, the internal and external accelerations at several R_e (about $10\,\mathrm{kpc}$ to the galactic center) are comparable to each other. For the intermediate case (model 3), the two accelerations become comparable at approximately 70 kpc, whereas for the weak external field case (model 2), this requires moving to very large radii, beyond hundreds of kiloparsecs.

To compute the static potential, we make use of the MONDian N-body solver NMODY (Ciotti et al. 2006; Nipoti et al. 2007). NMODY is a particle-mesh code that assigns particles by cloud-in-cell and solves the modified Poisson's equation on a spherical grid, using a second-order leap-frog scheme for time integration. We adopt a resolution of ${n}_{r}\times {n}_{\theta }\times {n}_{\phi }=512\times 128\times 256$ on a spherical grid, where the grid segments are defined as

$$\begin{eqnarray}{r}_{i} & = & 2\tan [(i+0.5)0.5\pi /({n}_{r}+1)]\,\mathrm{kpc},\\ {\theta }_{j} & = & \pi \times (j+0.5)/{n}_{\theta },\\ {\phi }_{k} & = & 2\pi \times k/{n}_{\phi },\end{eqnarray} \tag{ 6 }$$

where i = 0 .. n_r, j = 0 .. ${n}_{\theta }-1$, and k = 0 .. ${n}_{\phi }-1$.

Wu et al. (2008, 2010) showed that the internal potentials are prolate with respect to the direction of the external field. This effect becomes important when the internal and external fields are roughly of the same order; an even more significant effect occurs at weaker accelerations, in external field dominated regions where $| {{\boldsymbol{g}}}_{\mathrm{int}}| \ll | {{\boldsymbol{g}}}_{\mathrm{ext}}| \ll {a}_{0}$. In this case, one finds that the potentials are not only prolate but also appear distorted. Figure 1 shows isodensity and isopotential contours for model 4 as listed in Table 1. The contours correspond to radii of $2,5,10,15$, and 20 kpc along the major axis. We can easily identify a distortion of the potential along the direction of the external field at large radii, $r\sim 10\,\mathrm{kpc}$, where internal and external fields are comparable. The density contours, however, are still axisymmetric. Consequently, models constructed with external fields may not be necessarily self-consistent. Considering that the radius where the strongest external field (model 4) starts dominating the internal dynamics is at about 10 kpc (enclosing ≈82% of the total mass), the models should not be significantly affected for the most part and should reside in a quasi-equilibrium. Thus we expect Schwarzschild's method to be applicable in these cases.

Zoom In Zoom Out Reset image size

For a simplified view, let us consider a spherically symmetric system embedded into an external field. Integrating the MOND Poisson's equation, we arrive at the following expression (Bekenstein & Milgrom 1984):

$$\begin{eqnarray}&&\mu \left(\displaystyle \frac{| {\boldsymbol{g}}| }{{a}_{0}}\right){\boldsymbol{g}}={{\boldsymbol{g}}}_{{\rm{N}}}+{\boldsymbol{\nabla }}\times {\boldsymbol{h}},\end{eqnarray} \tag{ 7 }$$

where g_N is the Newtonian gravitational field, and ${\boldsymbol{\nabla }}\times {\boldsymbol{h}}$ is a solenoidal vector field determined by the condition that ${\boldsymbol{g}}$ can be expressed as the gradient of a scalar potential. Restricting ourselves to the axis parallel to the external field's direction, ${\boldsymbol{g}}$ and ${{\boldsymbol{g}}}_{{\rm{N}}}$ must be either parallel or antiparallel, assuming that the symmetry center coincides with the coordinate origin. Hence the curl term ${\boldsymbol{\nabla }}\times {\boldsymbol{h}}$ vanishes. The strength of the total gravitational acceleration $g=| {\boldsymbol{g}}|$ along the negative semi-axis in the external field's direction, that is, where the external field cancels part of the internal one, is ${g}^{-}=| {g}_{\mathrm{ext}}-{g}_{\mathrm{int}}|$, while on the positive semi-axis it is ${g}^{+}={g}_{\mathrm{int}}+{g}_{\mathrm{ext}}$. The two different sides have different values of the μ function at the same radii, leading to a larger MONDian enhancement of gravitation along the negative semi-axis. Clearly, for the underlying spherically symmetric density distribution, the potential and its derivatives are axisymmetric. Applied to a typical triaxial system, however, the result is approximately the same. For an external field pointing in an arbitrary direction, such a system has no symmetries anymore, but the curl term in Equation (7) only accounts for corrections on the level of 10% (Brada & Milgrom 1995).

Returning to the axisymmetric model 4, Figure 1 confirms the above considerations. On the left panel, the isopotential contours are denser in the first octant than in the third one. The internal potential is shallower in the first octant where ${g}^{+}=| {g}_{\mathrm{ext}}+{g}_{\mathrm{int}}|$ and steeper in the third octant where ${g}^{-}={g}_{\mathrm{int}}-{g}_{\mathrm{ext}}$. As can be seen from the right panel of Figure 1, the lopsidedness of the potential reaches its maximum at roughly 10 kpc, where the external and internal fields are comparable to each other. The semi x-axis ratio ${r}^{-}:{r}^{+}$ of the isopotential contours reaches about 1.14 between 10 and 11 kpc. At small radii ($r\lt 3\,\mathrm{kpc}$), the internal gravitational field dominates, and thus ${r}^{-}:{r}^{+}$ is basically 1. At much larger radii ($r\gg 10\,\mathrm{kpc}$), the relative contribution of the external field increases, and ${r}^{-}:{r}^{+}$ falls down to 1 again because the μ function approaches a constant, $\mu =\mu (| {{\boldsymbol{g}}}_{\mathrm{ext}}| /{a}_{0})$.

As stated in the previous section, all orbits are integrated for $100\,{T}_{\mathrm{cir}}$. In Figure 2, we plot the circular orbital time against the radius for all five models, using a logarithmic scaling. The circular orbital times of models 1 and 4 start to differ from each other at about 10 kpc. While the slope of model 1 is approaching unity at large radii, those of models 2–5 approach a value of 3/2 at infinity. For the isolated MOND model 1, the circular velocity turns constant at large radii, ${v}_{{\rm{c}}}={({{GMa}}_{0})}^{1/4}$, and hence ${T}_{\mathrm{cir}}\propto r$. In the case of strong external fields, however, the interpolating function μ becomes a constant far away from the system, leading to a Newtonian-like behavior ${T}_{\mathrm{cir}}\propto {r}^{3/2}$.

Zoom In Zoom Out Reset image size

Once the orbit libraries are built and the O_ij array is recorded, we can calculate the orbital weights in Equation (8). The right-hand side of Equation (8) denotes the mass of the $i\mathrm{th}$ cell, which we obtain from Monte Carlo simulations using the analytic density distribution. The linear system is then solved by applying a nonnegative least-squares (NNLS) method, which minimizes the following quantity:

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \frac{1}{{N}_{\mathrm{cell}}}\displaystyle \sum _{i=1}^{{N}_{\mathrm{cell}}}{\left({m}_{i}-\displaystyle \sum _{j=1}^{{N}_{\mathrm{orbit}}}{w}_{j}{O}_{{ij}}\right)}^{2}.\end{eqnarray} \tag{ 10 }$$

Furthermore, we introduce the self-consistency parameter δ (Merritt & Fridman 1996) as

$$\begin{eqnarray}&&\delta =\sqrt{{\chi }^{2}}/\overline{m},\end{eqnarray} \tag{ 11 }$$

where $\overline{m}$ is the mean Monte Carlo mass in cells. For self-consistent models, the value of δ is expected to quickly decrease with an increasing number of orbits, and it should be very close to zero if a large number of orbits is adopted. In Table 2, we list the self-consistency parameters of all five models. We find that the isolated model and the triaxial one in a strong external field are the most self-consistent. As a result of the broken symmetries, models in external fields should generally exhibit a lower level of self-consistency. For the axisymmetric models 2–4, the δ values are on the order of 10⁻². Compared to model 5, these systems feature only a single symmetry axis of the potential. The mass distribution reconstructed from the orbits in the distorted potential becomes lopsided in the outer parts and does not accurately reproduce the analytic density profile. This is in accordance with the observation that the estimated δ values grow with increasing external field strength.

Table 2.  Model Self-consistency and Equilibrium: Here ${\delta }_{\mathrm{smooth}}$ and δ Are the Self-consistency Parameters with and without Smoothing, Respectively, ${N}_{\gt 0}$ and ${N}_{\gt 0,\mathrm{smooth}}$ Are the Corresponding Orbit Numbers for Nonzero Weights (Obtained from Equations (10) and (12)), ${v}_{\mathrm{rms}}$ is the Root Mean Square Velocity in Units of $\,\mathrm{km}\,{{\rm{s}}}^{-1}$, and $-2K/W$ Is the Virial Ratio of the Initial Conditions

Model	1	2	3	4	5
${N}_{\mathrm{orbit}}$	6840	27360	27360	27360	13680
δ	$1.5\times {10}^{-15}$	$1.37178\times {10}^{-2}$	$2.44238\times {10}^{-2}$	$4.41140\times {10}^{-2}$	$6.2\times {10}^{-5}$
${\delta }_{\mathrm{smooth}}$	$6.7\times {10}^{-5}$	$1.37180\times {10}^{-2}$	$2.44239\times {10}^{-2}$	$4.41141\times {10}^{-2}$	$6.3\times {10}^{-5}$
${N}_{\gt 0}$	960	2378	2102	1677	1905
${N}_{\gt 0,\mathrm{smooth}}$	5995	2379	2103	1677	1949
${v}_{\mathrm{rms}}[\,\mathrm{km}\,{{\rm{s}}}^{-1}]$	225.55	225.20	224.67	220.28	309.58
$-2K/W$	1.01	1.02	1.01	1.00	1.01
ξ	1.26	1.74	1.57	1.65	1.60
${L}_{{\rm{c}}}[\,\mathrm{km}\,{{\rm{s}}}^{-1}\,\mathrm{kpc}]$	253.25	251.59	252.16	253.99	352.10

Note. The radial-to-tangential anisotropic ratio ξ characterizes the radial instability. Here, ${L}_{{\rm{c}}}$ denotes the unit-mass angular momentum with circular velocity ${v}_{{\rm{c}}}$ at radius ${r}_{{\rm{c}}}=1\,\mathrm{kpc}$, where ${L}_{{\rm{c}}}={r}_{{\rm{c}}}{v}_{{\rm{c}}}$ is expressed in units of $\,\mathrm{km}\,{{\rm{s}}}^{-1}\,\mathrm{kpc}$.

Download table as:  ASCIITypeset image

To construct sufficiently self-consistent models, one needs to use a large number of orbits, often far more than the number of cells ${N}_{\mathrm{cell}}$. The best solutions are typically nonunique with a very noisy phase-space distribution characterized by ${N}_{\mathrm{orbit}}-{N}_{\mathrm{cell}}$ and ${N}_{\mathrm{cell}}$ zero-weight and non-zero-weight orbits, respectively (Merritt & Fridman 1996; Zhao 1996; Rix et al. 1997). As the mass distributions given by Equation (4) are smooth, however, one would desire to select orbits in a less noisy way, that is, with little oscillation in the weights of neighboring orbits in phase space. Introducing a regularization mechanism allows one to construct a physically more plausible model. For instance, Zhao (1996) smoothed the orbits by averaging the weights of the nearest 26 neighboring orbits when solving the NNLS.

Here we apply a simpler method of regularization: we minimize the scatter of orbital weights by introducing a smoothing parameter λ, where $\lambda ={N}_{\mathrm{orbit}}^{-2}$ is chosen as in Zhao (1996). The regularization method used here is very similar to that of Merritt & Fridman (1996) and ensures the least number of orbits with zero weights. Hence, the fluctuations of weights become smaller, and the contribution of orbits to the mass distribution becomes smoother. To this end, Equation (10) is modified as

$$\begin{eqnarray}&&{\chi }_{\mathrm{smooth}}^{2}=\displaystyle \frac{1}{{N}_{\mathrm{cell}}}\displaystyle \sum _{i=1}^{{N}_{\mathrm{cell}}}{\left({m}_{i}-\displaystyle \sum _{j=1}^{{N}_{\mathrm{orbit}}}{w}_{j}{O}_{{ij}}\right)}^{2}+\lambda \displaystyle \sum _{j=1}^{{N}_{\mathrm{orbit}}}{w}_{j}^{2}.\end{eqnarray} \tag{ 12 }$$

Due to the regularization, the models acquire larger ${\chi }^{2}$ values, and therefore, the solution loses part of the self-consistency. The third line of Table 2 shows the self-consistency parameters, ${\delta }_{\mathrm{smooth}}=\sqrt{{\chi }_{\mathrm{smooth}}^{2}}/\overline{m}$ (where $\overline{m}$ is the mean Monte Carlo mass in cells), for the smoothed models. Comparing with the second line, we find that the regularization leads to an increase of ${\chi }^{2}$ on the order of ${10}^{-7}\mbox{--}{10}^{-5}$. Since the models 2–5 feature symmetrized orbits, that is, additional orbits starting from the other three octants (the second, fifth, and sixth octants) with the same ICs (see Appendix A), the increase of orbits is equivalent to smoothing the orbital structure. The regularization in Equation (12) only slightly changes the number of orbits with nonzero weights (see Table 2) and thus does not considerably change the accumulative fraction of individual orbit families.

Figure 3 shows the integrated contributions of orbits (for energies $\lt E$) to the system's mass with (right panels) and without (left panels) regularization. The individual contributions of long-axis loop (green), short-axis loop (bright blue), box (dark blue lines), and nonclassified (purple lines) orbits are plotted against the energy E. The fraction of box orbits in model 1 clearly increases when the regularization is applied. Such a large amount of box orbits might result in radial instability of the model. For models 2–5, the smoothing procedure does not change the fractions of orbital families noticeably; the numbers of orbits are large, and thus the solutions are smooth enough before applying the regularization. Reducing the number of cells and orbits to that of model 1, we further recomputed orbits for model 2 using Equation (12). The resulting orbital structure at low energy turns out very similar to the nonsmoothed solution of model 1, indicating that the use of reflecting (symmetry) planes changes the orbital structure for models 2–5. The total amount of orbits with low angular momentum, that is, box and nonclassified orbits, are quite large for these models. These orbits might introduce instability, which will be studied in later sections.

Zoom In Zoom Out Reset image size

For models 1–4, we find that short-axis loop orbits provide a large mass fraction, comprising over 40% of the total mass, even after the regularization. Long-axis loop orbits typically appear if an external field is applied (models 2–4). The stronger the external field, the more long-axis loop orbits emerge. Since the potential symmetry is broken at smaller radii (see Figure 1), this simply reflects that stronger external fields impact a larger fraction of orbits. The fraction of nonclassified orbits follows a similar behavior. While their fraction at fixed energies increases, the number of box orbits is simultaneously reduced. This might imply that the external field destroys the well-defined box orbits and makes them appear stochastic in phase space.

For the triaxial model 5, however, the orbital contributions significantly change. The smoothed result shows that the model is dominated by nonclassified orbits, making up almost 70% of the total mass. In Newtonian gravity, Merritt & Fridman (1996) have demonstrated that a galaxy constructed completely of regular orbits is not self-consistent, and most orbits in their galaxy models are irregular, especially for the cases with strong cusps. Our result is consistent with the Newtonian models of Merritt & Fridman (1996). The two families of loop orbits are the least important components, with mass fractions close to zero at all energy levels. Considering their total fraction in all five models, we conclude that nonclassified orbits become important for systems with lower symmetry.

As the solutions obtained from Schwarzschild's method are not unique, the orbital structure could change when adopting different regularization methods. The general conclusions, however, should remain the same.

4.1.1. Generating ICs from Orbital Libraries

We follow Zhao (1996) and Wu et al. (2009) to generate the ICs for N-body simulations. In brief, for an N-particle system, the number of particles on the $j\mathrm{th}$ orbit is ${n}_{j}={w}_{j}N$, where w_j is the weight of the $j\mathrm{th}$ orbit. Particles are placed on the $j\mathrm{th}$ orbit on equidistant times given by the interval ${\rm{\Delta }}{t}_{j}={T}_{j}/{n}_{j}$, where T_j is the total integrated time for the $j\mathrm{th}$ orbit (see Figure 2).⁸

Our galaxy models exhibit special symmetries that significantly decrease the amount of computing work (see Table 1 and Appendix A). Also, as the NNLS selects hundreds of nonclassified orbits that have not completely relaxed within 100 circular orbital times, the systems' phase-space symmetries are broken when placing particles on these orbits. Therefore we need to consider additional mirror particles in phase space, where the "mirrors" are the corresponding reflecting planes specified in Table 1. In the simulations, we use 10⁶ particles for each model after taking into account these symmetry considerations. These particles represent the inner 20 mass sectors (21 mass sectors in total) of a Hernquist model. The details of segmentation of the models are shown in Appendix A and in Wu et al. (2009).

If the ICs generated by Schwarzschild's technique are in quasi-equilibrium, the scalar virial theorem should be approximately valid, that is, $W+2K=0$, where

$$\begin{eqnarray}&&W=\int \rho {\boldsymbol{x}}\cdot {\boldsymbol{\nabla }}{\rm{\Phi }}{d}^{3}x\end{eqnarray} \tag{ 13 }$$

is the Clausius integral and K is the kinetic energy of the system (Binney & Tremaine 1987). The virial ratios $-2K/W$ and the root mean square velocities of the ICs are listed in Table 2. We find that all five models satisfy the scalar virial theorem very well, with $-2K/W=1\pm 0.02$. The ${v}_{\mathrm{rms}}$ values in models 2–4 are slightly smaller than in model 1, due to the presence of external fields. Since the potential is shallower in an external field, ${v}_{\mathrm{rms}}\approx \sigma$ for pressure-supported systems becomes smaller. For the strongest external field (model 4), the ${v}_{\mathrm{rms}}$ value is smallest, $\approx 5\,\mathrm{km}\,{{\rm{s}}}^{-1}$ less than in model 1. Even in this case, the decrease of ${v}_{\mathrm{rms}}$ is only 2%, implying that the dynamics is dominated by its self-gravity for the most part.

Consider a sizable low-mass galaxy dominated by an external field. Compared to the case of an isolated MOND galaxy or a CDM-dominated dwarf galaxy, the ${v}_{\mathrm{rms}}$ is expected to be much smaller. Crater II in the Local Group is such a galaxy and has recently been studied by McGaugh (2016). We know that Crater II is a very diffuse dwarf galaxy that is dominated by the weak external field of the Milky Way at a Galactic distance of 120 kpc. The predicted value of σ in this galaxy is only 2.1 $\,\mathrm{km}\,{{\rm{s}}}^{-1}$ if one accounts for this external field, but approximately twice as large for an isolated model. Given the magnitude of this effect, an analysis of such systems could be very rewarding. A detailed study of very diffuse systems that are entirely dominated by external fields is beyond the scope of this paper and will be the subject of a follow-up project.

4.1.2. Radial Instability of the Models

4.1.3. Numerical Setup and the Virial Ratio

Since the inclusion of external fields can be achieved by means of suitable boundary conditions, the NMODY Poisson solver does not need to be substantially altered and can be easily adapted to our purposes. For the simulations, we use a grid resolution of ${n}_{r}\times {n}_{\theta }\times {n}_{\phi }=1500\times 64\times 64$ in spherical coordinates ($r,\theta ,\phi$), where the radial grid segments are defined as

$$\begin{eqnarray}&&{r}_{i}=8\tan [(i+0.5)0.5\pi /({n}_{r}+1)]\,\mathrm{kpc},\end{eqnarray} \tag{ 14 }$$

and the other two remaining grid segments are the same as in Section 2. To reduce the computational workload, the chosen angular resolution is lower than that adopted for the Schwarzschild modeling. The upper panel of Figure 4 shows the static potential of model 3 along the $x,\,y$, and z axes for the high (black curves) and low resolutions (magenta curves). The two potentials agree well, with a relative difference of less than 3% within 40 kpc (where over 95% of the total mass is enclosed). Therefore, we do not expect any significant relaxation or evolution due to the reduction of angular resolution in the N-body simulations. We also note that the potential is very similar along the $x,\,y$, and z axes. This indicates that MOND potentials turn out rounder than their underlying density distributions in the regions where the dynamics is not dominated by the external field, an effect that has previously been studied in Wu et al. (2008). The time unit used in the NMODY code is given by (Wu et al. 2009)

$$\begin{eqnarray}{T}_{\mathrm{simu}} & = & {\left(\displaystyle \frac{{GM}}{{a}^{3}}\right)}^{-1/2}\\ & = & 4.7\times {10}^{6}\,\mathrm{year}{\left(\displaystyle \frac{M}{{10}^{10}{M}_{\odot }}\right)}^{-1/2}{\left(\displaystyle \frac{a}{1\mathrm{kpc}}\right)}^{3/2}.\end{eqnarray} \tag{ 15 }$$

Here the quantity $2\pi {T}_{\mathrm{simu}}$ has the physical meaning of one Keplerian time at a radius of r = a, and the system's total mass M is expressed in units of $1.0\times {10}^{10}\,{M}_{\odot }$.

Zoom In Zoom Out Reset image size

It is well known that the typical size of galaxy clusters is on the scale of several megaparsecs. However, their central regions, where strong and nearly uniform gravitational backgrounds exist, are much smaller. To give a rough estimate, the size is typically one order of magnitude smaller than the size of the cluster, which is around 0.1 Mpc. Galaxies are accelerated in an almost constant field at this scale. Converted into a physical timescale where this approximation holds, this gives around $60\,{T}_{\mathrm{simu}}\sim 0.3$ Gyrs. Of course, real galaxies are accelerated within inhomogeneous fields, but this general case is still too complex to be modeled at present, and most of the physics we are interested in at the moment can be explored in a constant background. More details about the time steps used in the code are discussed in Wu et al. (2009). We have simulated our models up to $120\,{T}_{\mathrm{simu}}$ (twice the value of $60\,{T}_{\mathrm{simu}}$) to examine the systems' behavior beyond the actual simulation time interval. In what follows, we will restrict the discussion to within $60\,{T}_{\mathrm{simu}};$ only virial ratios are presented for the fully simulated range of $120\,{T}_{\mathrm{simu}}$.

As mentioned above, the external field models are not exactly self-consistent with respect to the original analytic density profile. The virial ratios of the ICs slightly deviate from unity (±2%), and we expect these quasi-equilibrium ICs to quickly relax to dynamically virialized systems at the beginning of the simulations. Figure 5 shows the evolution of the virial ratio $2K/| W|$ within $120\,{T}_{\mathrm{simu}}$. As can be seen from the figure, the models revirialize within a few simulation times, and then the virial ratios oscillate around 1 for all five models, with an oscillation amplitude of roughly 0.05 within $120\,{T}_{\mathrm{simu}}$. Since the overall residuals of virial ratios at T = 0 are at the level of a few percent, the deviation from the exact equilibrium state is a minor effect. Hence we conclude that the virial theorem is valid for all considered models.

Zoom In Zoom Out Reset image size

The presence of an external field gives rise to a lopsided potential. Thus the mass density will redistribute inside the total potential until the density with its associated (internal) potential reaches an equilibrium configuration. In the following, we want to address how far-reaching this evolution in an external field is.

The left panels of Figure 6 illustrate the spherical radii ${r}_{0}(M)$ enclosing different fractions of the total mass M, that is, Lagrangian radii, increasing from 10% to 90% in steps of 10%. For models 1–3, these radii are very stable, showing only tiny oscillations within 60 simulation times. Only the innermost 10% mass radii of models 4 and 5 slightly increase by about 15%. This implies that the global radial mass distributions do not significantly evolve with time. For models in strong external fields, the radial mass distributions slightly decrease in the innermost regions and stay almost constant elsewhere.

Zoom In Zoom Out Reset image size

Note that there are fewer symmetries for models 2–5, and thus there are fewer mirror particles in phase space. Such a situation could result in a self-rotation that cannot be canceled out by the lack of counteracting mirror particles. This is especially true for the outer regions, where the impact of the external field starts to become important. In Appendix B.1, we will further comment on this issue. As a consequence, the major axes of models 2–5 may not coincide with the x axis anymore.

To this end, we determine the system's principal axes according to the following approach. Starting from an initial guess $p=q=1$, we consider all particles within an ellipsoid defined by ${x}^{2}+{(y/p)}^{2}+{(z/q)}^{2}={r}^{2}$. These particles are then used to compute the components of the weighted moments of inertia tensor given by

$$\begin{eqnarray}{\tilde{I}}_{{xx}}(r) & = & \displaystyle \frac{{\displaystyle \sum }_{i}\,{m}_{i}({y}_{i}^{2}+{z}_{i}^{2}){/r}_{i}^{2}}{{\displaystyle \sum }_{i}\,{m}_{i}},\\ {\tilde{I}}_{{xy}}(r) & = & \displaystyle \frac{{\displaystyle \sum }_{i}-{m}_{i}{x}_{i}{y}_{i}{/r}_{i}^{2}}{{\displaystyle \sum }_{i}\,{m}_{i}},\end{eqnarray} \tag{ 16 }$$

and similar expressions for the other components. Here we adopt the weighted moments of inertia tensor to mitigate the noisy contribution of particles at larger radii. The resulting inertia tensor is diagonalized, yielding eigenvalues ${\tilde{I}}_{x^{\prime} x^{\prime} }(r),\,{\tilde{I}}_{y^{\prime} y^{\prime} }(r)$, and ${\tilde{I}}_{z^{\prime} z^{\prime} }(r)$, where the primed coordinate system refers to the corresponding eigenframe. The associated principal axes,

$$\begin{eqnarray}\tilde{a}(r) & = & \sqrt{[{\tilde{I}}_{y^{\prime} y^{\prime} }(r)+{\tilde{I}}_{z^{\prime} z^{\prime} }(r)-{\tilde{I}}_{x^{\prime} x^{\prime} }(r)]/2},\\ \tilde{b}(r) & = & \sqrt{[{\tilde{I}}_{x^{\prime} x^{\prime} }(r)+{\tilde{I}}_{z^{\prime} z^{\prime} }(r)-{\tilde{I}}_{y^{\prime} y^{\prime} }(r)]/2},\\ \tilde{c}(r) & = & \sqrt{[{\tilde{I}}_{x^{\prime} x^{\prime} }(r)+{\tilde{I}}_{y^{\prime} y^{\prime} }(r)-{\tilde{I}}_{z^{\prime} z^{\prime} }(r)]/2},\end{eqnarray} \tag{ 17 }$$

are used to determine new values of $p\equiv \tilde{b}/\tilde{a}$ and $q\equiv \tilde{c}/\tilde{a}$, and the particle coordinates are rotated into the inertia tensor's eigenframe. The procedure is repeated iteratively until both axial ratios p and q converge to a relative deviation of less than 10⁻³.

The right panels of Figure 6 show the models' density distributions along their major axes obtained from the principal axes at a spheroidal radius enclosing 90% of the total mass, ${r}_{90 \% }$.⁹ The initial density distributions (black dotted curves) along the major axes agree well with the analytical densities of Hernquist profiles (black solid curves). The density profile of model 1 does not change substantially during the simulation. For models embedded in external fields, the density profiles quickly evolve to new profiles within $15\mbox{--}30\,{T}_{\mathrm{simu}}$ and reach stable states afterward. For model 2, the density along major axes within 1–5 kpc increases about 40% relative to the initial configuration. For models 3 and 4, the densities along the major axes decrease over the range $r\in [1.0,\,10.0]\,\mathrm{kpc}$ by approximately 50% and 60%, respectively. The density evolution becomes more pronounced for stronger external fields (model 4). The density profile of model 5 increases by about 60% in the inner region ($r\lt 5.0\,\mathrm{kpc}$), but there is no noticeable evolution in the outer parts. The differences in the evolution of density profiles along the major axes between models 2 and 5 are generally due to the varied strengths and directions of the external field.

Using the iterative procedure introduced in the previous section, we are able to obtain the shapes of the systems. The left panel of Figure 7 shows the axial ratios of all five models defined for ${r}_{90 \% }$, ${\tilde{p}}_{90 \% }={\tilde{b}}_{90 \% }/{\tilde{a}}_{90 \% }$, and ${\tilde{q}}_{90 \% }={\tilde{c}}_{90 \% }/{\tilde{a}}_{90 \% }$. The ratios ${\tilde{p}}_{90 \% }$ and ${\tilde{q}}_{90 \% }$ significantly evolve for all models during $60\,{T}_{\mathrm{simu}}$, where the evolution of ${\tilde{p}}_{90 \% }$ is generally stronger than that of ${\tilde{q}}_{90 \% }$. For all models, the values of ${\tilde{p}}_{90 \% }$ become considerably smaller. We observe that the axisymmetric models evolve into a triaxial configuration, and the initially triaxial model becomes slightly more prolate.

Zoom In Zoom Out Reset image size

Unlike the isolated triaxial models with axial ratios ${\tilde{a}}_{90 \% }:{\tilde{b}}_{90 \% }:{\tilde{c}}_{90 \% }=1.00:0.86:0.70$ studied in Wu et al. (2009), an instability appears for model 1 (black curves) within $60\,{T}_{\mathrm{simu}}$, yielding final axial ratios ${\tilde{a}}_{90 \% }:{\tilde{b}}_{90 \% }:{\tilde{c}}_{90 \% }\,=1.00:0.90:0.71$. There are oscillations of ${\tilde{q}}_{90 \% }$ around 0.7 with an amplitude of ±0.05 within $60\,{T}_{\mathrm{simu}}$. This suggests that the triaxial model in MOND is more stable than the axisymmetric one when the system is isolated.

For models 2–4 (axisymmetric models in external fields), however, the axial ratios have initially the same value as in model 1, ${\tilde{a}}_{90 \% }:{\tilde{b}}_{90 \% }:{\tilde{c}}_{90 \% }=1.0:1.0:0.7$. With external fields pointing in the (diagonal) positive or negative x–z direction, ${\tilde{a}}_{90 \% }$ and ${\tilde{b}}_{90 \% }$ start to separate from each other at the beginning of the simulations. For models 2 and 3, ${\tilde{p}}_{90 \% }$ decreases gradually to 0.85, and ${\tilde{q}}_{90 \% }$ drops slowly to around 0.65 at $60\,{T}_{\mathrm{simu}}$. For model 4, the values of ${\tilde{p}}_{90 \% }$ decrease to around 0.9 within the first $20\,{T}_{\mathrm{simu}}$ and stay constant afterwards. The ratio ${\tilde{q}}_{90 \% }$ oscillates around 0.75 within $20\,{T}_{\mathrm{simu}}$, and then jumps to 0.82 at $25\,{T}_{\mathrm{simu}}$. There is no significant evolution of ${\tilde{q}}_{90 \% }$ at later times. The different evolution of axial ratios betweenmodels $2,3$, and 4 is driven by the strength of the external field, which is strongest for model 4.

Models 2–4 eventually evolve into triaxial systems, and the final state is reached after $\sim 55\,{T}_{\mathrm{simu}}$ for models 2 and 3 and $\sim 25\,{T}_{\mathrm{simu}}$ for model 4. Note that the model analysis is made after rotating the systems into the reference frame defined by their principal axes. Hence, any effects of pattern rotation (see the evolution of angular momentum in Figure 20) due to the external field are eliminated.

For model 5 (green lines), the situation appears simpler. Initially, the axial ratios are approximately ${\tilde{a}}_{90 \% }:{\tilde{b}}_{90 \% }:{\tilde{c}}_{90 \% }\,=1.00:0.86:0.71$, which is close to the isolated case discussed in Wu et al. (2009). Both ${\tilde{p}}_{90 \% }$ and ${\tilde{q}}_{90 \% }$ start to increase within $10\,{T}_{\mathrm{simu}}$, but then decrease until $25\,{T}_{\mathrm{simu}}$, where ${\tilde{a}}_{90 \% }:{\tilde{b}}_{90 \% }:{\tilde{c}}_{90 \% }=1.00:0.87:0.69$. Again, there is no significant evolution beyond this point, and one finds ${\tilde{a}}_{90 \% }:{\tilde{b}}_{90 \% }:{\tilde{c}}_{90 \% }=1:00:0.85:0.69$ at $60\,{T}_{\mathrm{simu}}$. Since the external field for this model points in the x direction, that is, perpendicular to the y–z plane, the curves for the components along the y and z axes display a very similar running behavior. The instability appears within the first $25\,{T}_{\mathrm{simu}}$ for models 4 and 5, after which they become stable. We also note that the evolution formodels 1–3 is not too different, since the external field effects in these cases (models 2 and 3) are mild. Again, this shows that the evolution of axial ratios is closely related to the external field strength.

The axial ratios as a function of enclosed masses are presented for all five models in the right panel of Figure 7. For the initial models (solid curves), the axial ratio agrees well with the analytic density distribution and does not evolve considerably with increasing mass for the isolated model (model 1). For the models in external fields (models 2–4), there are deviations between the analytical and model's principal axes in the inner regions. There are small oscillations or spikes for the resulting $p=\tilde{b}/\tilde{a}$ and $q=\tilde{c}/\tilde{a}$ in model 5, with amplitudes $\lt \pm 0.1$. These spikes are purely numerical features and emerge from the iterative procedure, used to determine the principal axes, which is quite sensitive to the chosen initial guess in this case.

The initial axial ratios of models 2–5 deviate from the original analytic models, especially within ${r}_{20 \% }$, the Lagrangian radius enclosing 20% of the total mass. These effects are caused by the offsets between the cusp centers and the center of masses (CoMs) of the models (Wu et al. 2010). Considering models with an external field applied, the systems are accelerated and moving due to the constant background field. To keep the angular resolution of the simulation at a reasonable size, the galactic center has to be placed at the computational grid's center. In our simulations, we transformed the coordinates at every time step, moving the CoM to the grid center and changing to the frame where its velocity is zero. The CoM does not need to coincide with the galactic center, the center of the cusp. Due to the lopsided potential, such an offset within the density distribution can develop within MOND. However, this will happen neither in Newtonian gravity nor for isolated models in MOND. To better understand this effect, Wu et al. (2010) have designed a simple experiment: a spherically symmetric Plummer model is placed into a homogeneous external field. In the case of Newtonian gravity, the superposition principle applies, and the whole system is equally accelerated along the external field's direction. Hence the position of the CoM does not move in the internal system. For MOND, however, the internal gravity is determined by both the external field and the internal matter distribution. Hence, the outer parts of the galaxy will become lopsided, generating an additional external field itself that will act on the inner part. The additional field will cause the CoM to slightly shift away from the point where gravity equals the external field, that is, the point where internal gravitational forces cancel. Besides discrepancies in the innermost regions, the axial ratios as a function of enclosed masses differ from the analytic axis-symmetric models in external fields also at larger Lagrangian radii. For instance, $p\approx 0.96$ and $q\approx 0.74$ in model 2. Such deviations show that our ICs for the N-body simulations do not perfectly describe the shapes of the original analytic axis-symmetric models embedded in external fields. This is likely due to departures from self-consistency as the corresponding δ parameters of models 2–4 are around a few percent and the ICs of the models are in quasi-equilibrium. The shapes of models 1 and 5 (again, in the outer regions) agree very well with that of the analytic models, and their self-consistency parameters are about three orders of magnitudes smaller.

At $T=60\,{T}_{\mathrm{simu}}$ (dotted curves), the values of p for models 1–3 drop to $\approx 0.94$ in the inner region and to ≈0.86–0.90 at ${r}_{90 \% }$. The values of q do not change as much as that of p, with amplitudes ${\lt }_{-0.05}^{+0.02}$. This agrees with the results in the left panel of Figure 7. Models 1–3 evolve into triaxial configurations due to the radial instability. The evolution of model 4 is much more striking. The value of p for the innermost particles within ${r}_{5 \% }$ is about 0.98. It then decreases to 0.88 at ${r}_{40 \% }$ and jumps up to 0.92 for an enclosed mass of 45%. The values of p stay almost constant in the outer region of the model. The sudden increase of p is an effect of the offset between the CoM and the cusp center. The values of q are larger than that of the ICs. At ${r}_{5 \% }$, $q\approx 0.90$ and decreases to around 0.82 at ${r}_{90 \% }$. The instability is mainly caused by nonclassified orbits (see the orbital fraction in the right panels of Figure 3) and changes the shape of the model. Model 4 turns out more prolate when compared to models 1–3. For model 5, the axial ratio p does not evolve substantially in the strong external field. The values of q are generally about 0.05 lower than in the ICs. Compared to model 4, the shape evolution of model 5 is considerably suppressed. There are two possible reasons: either a triaxial model is more stable than an axisymmetric model in MOND, or the symmetry of model 4 is broken in a more complex way since the external field direction is along the negative x–z diagonal.

To summarize, our results imply that

• 1.  
  an axisymmetric model with and without the external field effect is not stable in MOND, and that
• 2.  
  the principal axes, describing the shape of a galaxy, generally evolve more substantially in the presence of a strong external field, and that
• 3.  
  radial instability caused by both box orbits and nonclassified orbits with low angular momentum can change the axial ratios, that is, the shape of a galaxy. The origin of the instability will be studied later in Section 5.

From Section 4.2, we already know that the major axes of inner regions of the galaxy models in an external field are slightly expanded compared to the isolated case. Additional information on these models can be obtained by exploring their properties in the full phase space. To study the dynamical evolution of our systems, we calculate the radial velocity dispersion ${\sigma }_{r}(r)$ and the anisotropy parameter

$$\begin{eqnarray}&&\beta (r)=1-\displaystyle \frac{{\sigma }_{\theta }^{2}+{\sigma }_{\phi }^{2}}{2{\sigma }_{r}^{2}},\end{eqnarray} \tag{ 18 }$$

where r is defined in Equation (4) and ${\sigma }_{\theta }$ and ${\sigma }_{\phi }$ are the tangential and azimuthal velocity dispersions, respectively.

The left panels of Figure 8 show the time evolution of the radial velocity dispersion ${\sigma }_{r}(r)$.¹⁰ Obviously, all models start off with an instability and seem to reach a stable state after a short evolutionary phase. The dispersions ${\sigma }_{r}(r)$ of all models converge to stable states within $15\mbox{--}30{T}_{\mathrm{simu}}$. In addition, we find that ${\sigma }_{r}(r)$ is flattened in the inner parts ($\lt 1.0$ kpc) of all systems, with ${\sigma }_{r}(r\lt 1.0\,\mathrm{kpc})\approx 170\,\mathrm{km}\,{{\rm{s}}}^{-1}$ for models 1–4. As expected, the models in external fields (i.e., models 2–5) exhibit more evolution than the isolated model 1. For the same density model (models 1–4), a stronger external field also leads to more evolution in the ${\sigma }_{r}(r)$ profiles. Velocity dispersions in model 5 are much larger, since the mass of model 5 is a factor of 2 larger than other models. The ${\sigma }_{r}(r)$ profile of the triaxial model 5 slightly increases after the relaxation within $30\,{T}_{\mathrm{simu}}$. At larger radii, dispersions ${\sigma }_{r}(r\gt 1\,\mathrm{kpc})$ are about $10\mbox{--}20\,\mathrm{km}\,{{\rm{s}}}^{-1}$ higher than in the initial state. The inner region of ${\sigma }_{r}(r)$ for model 5 decreases from 300 to 250 $\,\mathrm{km}\,{{\rm{s}}}^{-1}$ within $30{T}_{\mathrm{simu}}$, and the profile also appears flat within 1 kpc. The evolution from the ICs to the final stable state indicates a radial instability, which might originate from box orbits or from nonclassified orbits with low angular momentum.

Zoom In Zoom Out Reset image size

For the axisymmetric models 1–4, the profiles of the anisotropy parameter (right panels of Figure 8), $\beta (r)$, are quite different for the isolated and nonisolated models. In the case of an isolated galaxy, model 1, $\beta (r)$ is initially almost constant, $\approx 0.2$, in accordance with the results of Wu et al. (2009). When there is an external field (models 2–4), however, the initial $\beta (r)$ profiles are no longer constant. In general, $\beta (r)$ declines in the inner region and starts to increase from $\simeq 0.5\,\mathrm{kpc}$ to $\simeq 4\,\mathrm{kpc}$. The model with the weakest external field has the highest radial anisotropy at 4 kpc, $\beta (r\geqslant 5\,\mathrm{kpc})\approx 0.8$. The $\beta (r)$ profile stays almost constant in the outer region, where $r\gt 5\,\mathrm{kpc}$. Given that the external field is very weak in model 2, it seems unexpected that there is such a large difference between models 1 and 2. The main reason is that the numbers of cells in the grids and the numbers of orbits in the orbital libraries are different. As previously mentioned, there are three times more cells and orbits in models 2–4 (see Table 1). The increase of cells and orbits is equivalent to additional smoothing for the orbital structure. To examine the impact of increasing the number of cells and orbits, we constructed a new model, labeled $\tilde{1}$, by launching orbits from the same octants as in model 2 and increasing the number of orbits to 27,360. Generating the corresponding N-body ICs as before, we then computed the anisotropy profile of model $\tilde{1}$, and the result is shown in Figure 9. We find that $\beta (r)$ for model $\tilde{1}$ is very similar to that of model 2, confirming that the different numbers of cells and orbits yield the observed discrepancy between the $\beta (r)$ profiles of models 1 and 2 at T = 0.

Zoom In Zoom Out Reset image size

When the external fields become stronger (models 3 and 4), the maximal values of $\beta (r)$ are about 0.6 at 4 kpc, and the $\beta (r)$ profiles fall off again at larger radii. The $\beta (r)$ profiles for ICs with the strongest external fields, for example, model 4, have the smallest values at large radii among the external field models. As the external field amplitude increases, the radial anisotropy is reduced in the outer region of the models.

It is important to keep in mind that the solution of Schwarzschild's method is not unique. Therefore, possible changes in the orbital structure could lead to different anisotropy profiles of the velocity dispersion. For all the models in this work, however, the regularization method is the same. It seems safe to conclude that for models of the same mass derived from the regularization in Equation (12), stronger external fields will generally lead to less radially anisotropic models at larger radii.

The colored curves show the $\beta (r)$ profiles at different simulation times. For the isolated model, the $\beta (r)$ profile does not evolve very much, which is again similar to the results presented in Wu et al. (2009). For the axisymmetric models embedded in external fields along the diagonal x–z direction, $\beta (r)$ does not substantially evolve in the outer regions where $r\gt 1.0\,\mathrm{kpc}$, but decreases in the inner region within $15\mbox{--}30\,{T}_{\mathrm{simu}}$, after which it remains stable. The models with the strongest external field show the most prominent evolution within 1 kpc. At $60\,{T}_{\mathrm{simu}}$, $\beta (r)$ is nearly zero in the central region; that is, the center is basically isotropic. This is consistent with the left panels of Figure 8. Model 5 eventually becomes less radially anisotropic in the inner region and more anisotropic at larger radii. Since the fraction of box orbits is quite small, the numerous nonclassified orbits might be the origin of the radial instability for triaxial model 5.

Our results are consistent with the shape evolution illustrated in Figure 7. Systems embedded in external fields appear stable after an early evolution caused by radial instability. Nevertheless, more dynamical quantities need to be investigated before making any such statements. Further dynamical studies related to kinetic energy and angular momentum are presented in Appendix B.

The analysis in Section 4 showed that axisymmetric models are unstable and that the models exhibit large fractions of box and nonclassified orbits with low angular momentum. In order to examine the origin of the instability, we shall perform a further test here. In our test, the box orbits are removed from the global orbital library of models 1–4. Thus box orbits no longer contribute to the mass of these systems, and the corresponding models are labeled with a prime: 1', 2', 3', and 4'. The numbers of orbits, ${N}_{\mathrm{orbit}}^{\prime }$, in the new libraries are listed in Table 3. There are at least 4000 orbits for the axisymmetric models. The best-fit solutions of Equation (12) are then computed based on the new orbital library. The new ${\delta }_{\mathrm{smooth}}$ parameters are on the order of 10⁻² and are all larger than for models that include box orbits. This is reasonable because removing these orbits is equivalent to introducing additional constraints on the weights for the original orbital library. The number of nonzero weights in the new models is smaller than in the original ones. Figure 10 shows the orbital structure of the four new axisymmetric models. In these models, the fraction of short-axis loop orbits generally increases. The isolated system (model 1') is almost entirely composed of short-axis loop orbits. The fraction of nonclassified and long-axis orbits increases with growing external field strength. Model 2' is dominated by short-axis (≈60%) and nonclassified orbits (≈40%). The fractions of nonclassified orbits in models 3' and 4' (contributing over 50% of the total mass) are larger than that of short-axis orbits. The fraction of long-axis loop orbits reaches up to 12% in model 4'.

Zoom In Zoom Out Reset image size

Table 3.  Parameters of the New Orbital Library after Removing the Box Orbits

Model	1'	2'	3'	4'
${N}_{\mathrm{orbit}}^{\prime }$	4178	19361	21424	23372
${\delta }_{\mathrm{smooth}}$	$1.52\times {10}^{-2}$	$5.33\times {10}^{-2}$	$4.32\times {10}^{-2}$	$4.90\times {10}^{-2}$
${N}_{\gt 0,\mathrm{smooth}}^{\prime }$	617	1517	1601	1465
$-2K/W$	1.01	1.01	1.01	1.00
ξ	0.98	1.58	1.79	1.65

Note. The number of total orbits is denoted by ${N}_{\mathrm{orbit}}^{\prime }$, ${\delta }_{\mathrm{smooth}}$ is the self-consistency parameter, and ${N}_{\gt 0,\mathrm{smooth}}^{\prime }$ is the number of orbits with nonzero weights. The values of the virial ratio, $-2K/W$, and of the radial stability parameter, $\xi$, are also listed.

Download table as:  ASCIITypeset image

Generating N-body ICs according to the scheme presented in Section 4.1.1, we have conducted the same stability tests as before. The virial ratios of the new ICs are listed in Table 3. Again, all models satisfy the scalar virial theorem very well, which implies that the new ICs are close to equilibrium. The first step of the stability test is to examine the ξ values of the new models (see Table 3). We find that $\xi \ll {\xi }_{\mathrm{crit}}$ for model 1' while $\xi \in [1.58,\,1.79]$ for the other ones. These values lie within the ξ range estimated for the original models 2–4. While model 1' is expected to be stable, the situation for the other models is unclear and needs to be analyzed in more detail.

The shape evolution of the new models at ${r}_{90 \% }$ for $T=0\mbox{--}60\,{T}_{\mathrm{simu}}$ is presented in the left panel of Figure 11. In contrast to the original isolated model with box orbits, the global axial ratios for model 1' do not evolve with time. The axial ratios as a function of enclosed mass are shown in the right panel of Figure 11. The black dotted curves confirm that the shape of model 1' has not changed after $60\,{T}_{\mathrm{simu}}$. Therefore, model 1', which is mostly composed of short-axis loop orbits and has only a tiny fraction of nonclassified orbits, is stable. The global axial ratios, ${\tilde{p}}_{90 \% }$ and ${\tilde{q}}_{90 \% }$, of models 2'–4' evolve similarly to that of the original models 2–4, indicating that these models are not stable. The instability might arise from the large population of nonclassified orbits that contain little angular momentum. The right panel of Figure 11 shows that the shape of model 2' at ${r}_{40 \% }$ has not significantly changed after $60\,{T}_{\mathrm{simu}}$. The values of p and q slowly decline with increasing enclosed mass, but rise again at ${r}_{80 \% }$. At ${r}_{90 \% }$, the axial ratios are $\tilde{a}:\tilde{b}:\tilde{c}=1.00:0.96:0.66$. The instability of model 2' is driven by nonclassified orbits that are redistributed to higher energies. The axial ratio of model 3' at ${r}_{90 \% }$ declines to 1.0:0.64:0.52 at $60\,{T}_{\mathrm{simu}}$. Compared to the results for model 3, this indicates that the shape of model 3' becomes more prolate. For model 4', the axial ratios as a function of enclosed mass are very similar to that of model 4. Since the fraction of box orbits in model 4 is less than 10%, removing them does not substantially affect the evolution of the model. Therefore, we can infer that nonclassified orbits play an important role for the instability in this case.

Zoom In Zoom Out Reset image size

Nonclassified orbits are defined as low angular momentum orbits in which the velocity cannot be restricted to a rectangular velocity box (Equation (21)). These orbits behave stochastically in grid space. In order to characterize the origin of the instability in models 2'–4', we perform a further test by removing all nonclassified orbits from their orbital library. The such-obtained models are mostly composed of short-axis loop orbits and are labeled as models 2'', 3'', and 4''. Their orbital structures are illustrated in the upper panels of Figure 12. The models exhibit a small fraction of long-axis loop orbits, and their amount grows with increasing external field strength. It is known that long-axis loop orbits are illegal orbits in both Newtonian axisymmetric models and isolated axisymmetric MOND models. Their existence might cause additional instability, which we will investigate further below. The resulting self-consistency parameters δ are listed in Table 4. In all cases, $\delta \sim 0.1$, which is significantly larger than for models 2–4 and 2'–4'. This implies that the N-body ICs sampled from these Schwarzschild models are not in equilibrium. Therefore, nonclassified orbits are a necessary ingredient for the self-consistency of models 2–4. The stability of these models is again examined using N-body runs, as described in Section 5.1. The estimated axial ratios at ${r}_{90 \% }$ as a function of time and as a function of the enclosed mass at T = 0 and $T=60\,{T}_{\mathrm{simu}}$ are illustrated in the upper panels of Figure 13. We find that the models are still axisymmetric after $60\,{T}_{\mathrm{simu}}$. No triaxial configurations emerge during the N-body experiments, which agrees with the small ξ values estimated for these models. Ignoring the innermost particles within ${r}_{20 \% }$, the axial ratios of model 2'' at $60\,{T}_{\mathrm{simu}}$ turn out very similar to that of model 1', indicating that by removing box and nonclassified orbits, axisymmetric models embedded in weak external fields evolve almost like isolated ones. Nevertheless, the models eventually become rounder, where the effect is more prominent for stronger external fields. This might be a result of the model's considerable departure from self-consistency and should be further explored in future work. Our results suggest that nonclassified orbits with low angular momentum likely play an important role in the instability in the original models 3 and 4 discussed in Section 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4.  Parameters of the New Orbital Libraries after Removing Box and Nonclassified Orbits

Model	2''	3''	4''
${N}_{\mathrm{orbit}}^{\prime\prime }$	12000	11979	11297
${\delta }_{\mathrm{smooth}}^{\prime\prime }$	0.087	0.087	0.108
${N}_{\gt 0,\mathrm{smooth}}^{\prime\prime }$	6077	5274	820
$-2K/W$	1.01	1.01	1.01
ξ	0.98	0.83	0.98

Note. The number of total orbits is denoted by ${N}_{\mathrm{orbit}}^{\prime\prime }$, ${\delta }_{\mathrm{smooth}}$ is the self-consistency parameter, and ${N}_{\gt 0,\mathrm{smooth}}^{\prime\prime }$ is the number of orbits with nonzero weights. The values of the virial ratio, $-2K/W$, and of the radial stability parameter, ξ, are also listed.

Download table as:  ASCIITypeset image

As mentioned before, external fields break the symmetry of the potential and introduce long-axis loop orbits. However, the density distribution is still axisymmetric, which conflicts with the presence of long-axis loop orbits. One question that naturally arises is the following: will these long-axis loop orbits import instability? To address this question, we compare models 2'–4' to another set of models in which these illegal orbits are removed. The latter are labeled as models 2∗, 3∗, and 4∗, and their orbital structure and evolution are analyzed following the same procedure as in Section 5.1. The orbital structure of models 2∗–4∗ is shown in the lower panel of Figure 12. These additional models are composed of only short-axis loop and nonclassified orbits. As can be read from Table 5, all models yield self-consistency parameters $\delta \approx 0.05$, which is close to the corresponding values for models 2'–4'. If long-axis loop orbits introduced another source of instability, the inferred axial ratios of models 2∗–4∗ should differ from that of models 2'–4' where only box orbits have been removed. The axial ratios of models 2∗–4∗ are illustrated in the lower panel of Figure 13. The time evolution of axial ratios at ${r}_{90 \% }$ (lower left panel) is quite similar to that found in models 2'–4'. Moreover, the axial ratios as a function of enclosed mass at T = 0 and $T=60\,{T}_{\mathrm{simu}}$ (lower right panel) are not significantly altered relative to that found for models 2'–4'. Consequently, long-axis loop orbits do not contribute additional instability to the models embedded in external fields. There are, however, small differences in the final configurations between models 2'–4' and 2∗–4∗. For instance, the ratio p for model 2∗ declines mildly (and smoothly) from ${r}_{50 \% }$ to $p\approx 0.96$ at ${r}_{90 \% }$, whereas that of model 2' drops more quickly from ${r}_{40 \% }$ and suddenly increases at ${r}_{80 \% }$. Although the orbital structures of models 2' and 2∗ are nearly identical, there is an ≈1% fraction of long-axis loop orbits in model 2' and, accordingly, the fraction of nonclassified orbits is slightly larger (by about ≈0.26%). The long-axis loop orbits are mainly distributed at larger radii where about 80% and more of the total mass is enclosed. This is likely responsible for the sudden rise at ${r}_{80 \% }$, that is, a better agreement with the axial ratios of the original profiles. The slight changes in the orbital structure and nonclassified orbit fractions are probably the main reasons for the observed differences between the two models. For models 4∗ and 4', the changes in the orbital structure are more pronounced, with a 10% fraction of long-axis loop orbits in model 4' and an increased ≈10% fraction of nonclassified orbits in model 4∗. While the shape of model 4' at ${r}_{50 \% }$ is more prolate at $T=60\,{T}_{\mathrm{simu}}$, that of model 4∗ turns out to be more triaxial. The difference between models 4' and 4∗ is due to the different amounts of nonclassified orbits.

Table 5.  Parameters of the Orbital Libraries after Removing Box and Long-axis Loop Orbits

Model	2∗	3∗	4∗
${N}_{\mathrm{orbit}}\ast$	19097	20381	21079
${\delta }_{\mathrm{smooth}}\ast$	0.053	0.046	0.059
${N}_{\gt 0,\mathrm{smooth}}\ast$	1502	1570	1385
$-2K/W$	1.01	1.01	1.00
ξ	1.58	1.81	1.80

Note. The number of total orbits is denoted by ${N}_{\mathrm{orbit}}\ast$, ${\delta }_{\mathrm{smooth}}$ is the self-consistency parameter, and ${N}_{\gt 0,\mathrm{smooth}}\ast$ is the number of orbits with nonzero weights. The values of the virial ratio, $-2K/W$, and of the radial stability parameter, ξ, are also listed.

Download table as:  ASCIITypeset image

To summarize the results of this section, we conclude that both box and nonclassified orbits result in instability for isolated models and for models that are embedded in external fields.

We use an orbital integration code that adopts a 7/8-order Runge–Kutta method (Fehlberg 1968) to ensure sufficient accuracy of the orbits (Wang et al. 2008). The way we divide the spatial grid is exactly the same as in Wang et al. (2008; also see Wu et al. 2009 for numerical details). Here we consider a total of 1008 equal-mass cells. The outermost 48 cells are not taken into account because their boundaries extend to infinity. Since the density evolves as $\rho \sim {r}^{-4}$ at large radii, the orbits in this sector should contribute much less than all other orbits, and thus they are negligible.¹¹

In a previous study, Wang et al. (2008) and Wu et al. (2009) considered isolated triaxial models. Due to the three-fold symmetry, they took only the first octant into account and then symmetrized the O_ij by reflecting all orbits at the octant's boundaries, the x–y, x–z, and y–z planes. We use the same approach for the isolated model 1.

As the presence of an external field breaks the symmetry of the resulting potentials, we cannot simply reflect the orbits at these planes anymore. For model 5, the external field points in the negative x direction, and thus orbits starting from the positive x semispace are not the mirror orbits of those starting from the negative x semispace; that is, orbits on different sides of the y–z plane behave differently. However, the orbits will still have a two-fold symmetry with respect to the x–z and x–y planes. Hence we used the first ($| x| ,| y| ,| z|$) and second ($-| x| ,| y| ,| z|$) octants to calculate orbits for model 5. The total number of cells in this case is set to $960\times 2=1920$, where the outermost 96 cells are again excluded.

When the external field is pointing in an arbitrary direction, the symmetry of the model is further reduced. Switching to an axisymmetric model and forcing the external field to point in a direction perpendicular to the symmetry axis (see any of the models 2–4 in Table 1), however, we can always find a coordinate frame such that the external field is parallel to the x axis. As the potential will be distorted along the external field direction, we lose the symmetry along the x axis in addition to the one along the z axis. However, the potential will still be symmetric with respect to the x–z plane, which can be exploited to simplify the numerical computation. We fold the system at the x–z plane and consider four octants: the first ($| x| ,| y| ,| z|$), second ($-| x| ,| y| ,| z|$), fifth ($| x| ,| y| ,-| z|$), and sixth ($-| x| ,| y| ,-| z|$). We calculate orbits by reflecting them at the x–z plane. In this case, we have $960\times 4=3840$ mass cells (the outermost 192 cells are excluded) dividing the semispace. Note that if the system is triaxial and the external field direction does not coincide with one of the axes, we need to consider all eight octants, using $960\times 8=7680$ cells. Since this is computationally expensive, we limit ourselves to axisymmetric systems (models 2–4) and a single triaxial system where the external field is antiparallel to the x direction.

Following Schwarzschild (1993) and Merritt & Fridman (1996), we set up the ICs for the orbit library. The numerical details of this procedure can be found in Wang et al. (2008) and Wu et al. (2009). Here we only give a brief description. There are two sets of ICs for the starting points:

• 1.  
  Stationary starting orbits from 20 equal-potential surfaces, corresponding to the 20 nodes of equal-mass sectors on the x axis. These are freely falling into the system's internal potential with their initial velocities set to zero. Since these orbits do not carry any angular momentum initially, they may cross the center of the system and change the sign of their angular momenta. It is known that stationary ICs can produce box orbits. Defining a velocity box through

  $$\begin{eqnarray}&&\displaystyle \frac{\max (3{v}_{x}^{2})\cdot \max (3{v}_{y}^{2})\cdot \max (3{v}_{z}^{2})}{\max {({v}_{x}^{2}+{v}_{y}^{2}+{v}_{z}^{2})}^{3}}\lt 1,\end{eqnarray} \tag{ 21 }$$

  we classify orbits as box orbits if they satisfy the above inequality, which describes an almost rectangular box inside a maximum-energy sphere in velocity space. Note that the time-averaged angular momentum is zero.
• 2.  
  Ejecting orbits starting from the x–z plane with initial velocities $({v}_{x},{v}_{y},{v}_{z})=(0,\sqrt{2(E-{{\rm{\Phi }}}_{\mathrm{int}})},0)$. Here E corresponds to the total potential energy of the 20 equal-mass sectors along the x axis, and ${{\rm{\Phi }}}_{\mathrm{int}}$ is the (smooth) internal potential. Since most of the ejecting orbits are loop orbits, they cannot enter the system's central part. There are two families of loop orbits since the initial components L_x and L_z are not zero. We classify orbits that conserve the sign of their angular momentum around the long axis during the whole simulation as the long-axis loop orbits. Alternatively, this is expressed in terms of

  $$\begin{eqnarray}&&\max ({L}_{x})\cdot \min ({L}_{x})\gt 0.\end{eqnarray} \tag{ 22 }$$

  Similarly, we define short-axis loop orbits as

  $$\begin{eqnarray}&&\max ({L}_{z})\cdot \min ({L}_{z})\gt 0.\end{eqnarray} \tag{ 23 }$$

  These orbits cannot move across the center of the system because their angular momentum at this point is zero. Other types of orbits, which are beyond the definition of boxes and loops, are referred to as nonclassified orbits.

With the above ICs, one can generate most of the orbits in full phase space for a given potential (Schwarzschild 1993). The starting points of the orbits for model 1 are the same as in Wu et al. (2009). For model 5, we add mirror positions in the second octant, and for models 2–4, we use four mirror octants as mentioned above. In Table 1, we summarize the ICs for the five models, where ${N}_{\mathrm{stationary}}$ and ${N}_{\mathrm{ejecting}}$ are the numbers of the starting points. Note that with a higher number of cells and orbits, both the size of the array O_ij and the linear system (Equation (8)) increase, and building the orbit library becomes computationally more expensive, as has been specified in Section 3.1.

Since the inner and outermost parts of nonisolated models exhibit some radial instability, probably causing a few particles to leave the system during this phase, we only consider the remaining (more stable) fraction for further study. Since the left panels of Figure 6 show that the 90% Lagrangian radii for all of the models stay constant within $60\,{T}_{\mathrm{simu}}$, we shall use the particles within ${r}_{90 \% }$.

As the velocities redistribute substantially in models 2–5, the kinetic energy of the systems could also change by a significant amount. The left panel of Figure 19 illustrates the kinetic energy tensor components ${K}_{{xx}}$, ${K}_{{yy}}$, and ${K}_{{zz}}$, which are given by

$$\begin{eqnarray}&&{K}_{{xx}}=\displaystyle \frac{1}{2}\displaystyle \frac{{\displaystyle \sum }_{i}\,{m}_{i}{v}_{x}^{2}}{{\displaystyle \sum }_{i}\,{m}_{i}}\end{eqnarray} \tag{ 24 }$$

for $r\lt {r}_{90 \% }$ and analogously for the other directions. We highlight that ${K}_{{xx}}$, ${K}_{{yy}}$, and ${K}_{{zz}}$ are computed with respect to the original $(x,\,y,\,z)$ coordinate frame instead of the eigenframe considered in Section 4.3. Since all particles have equal mass ${m}_{i}=M/N$, Equation (24) simplifies to ${K}_{{xx}}=0.5\langle {v}_{x}^{2}\rangle$. We find that model 1 (black lines) does not exhibit any prominent alterations during $60\,{T}_{\mathrm{simu}}$. The components ${K}_{{xx}}$, ${K}_{{yy}}$, and ${K}_{{zz}}$ of model 1 evolve almost identically. Compared to the initial values, oscillations in the individual components are smaller than 10% within $60\,{T}_{\mathrm{simu}}$. At the beginning, the kinetic energy components ${K}_{{xx}}$ (solid lines) and ${K}_{{yy}}$ (dotted line) of model 2 (magenta), model 3 (cyan), and model 4 (yellow) coincide since the ICs are axisymmetric. The component ${K}_{{xx}}$ does not particularly change with time for these models. The amplitudes of oscillations around the original values are less than 6%. Although initially having the same density profile, the components ${K}_{{yy}}$ of models 2 and 3 evolve quite differently from that of model 1. The component ${K}_{{yy}}$ of models 2 and 3 slightly drops from about $20\,{T}_{\mathrm{simu}}$. Compared to their initial values, ${K}_{{yy}}$ decreases by about 5% and 13% for models 2 and 3, respectively, after $60\,{T}_{\mathrm{simu}}$. For model 4, where the external field is the strongest, ${K}_{{yy}}$ falls from about $50\,{T}_{\mathrm{simu}}$ and becomes around 8% smaller than at T = 0. The component ${K}_{{zz}}$ of model 2 starts to decrease from ${10}_{\mathrm{simu}}$ and is reduced by about 8% at the end of the simulation; for models 3 and 4, it oscillates around the initial values with amplitudes smaller than 10%. Variations of the kinetic energy components signal model instability. In Appendix A, we have seen that the fraction of box orbits in models 2–4 decreases while the fraction of nonclassified orbits increases. The redistribution of kinetic energy for these axisymmetric models might mainly be influenced by instability due to box orbits.

Zoom In Zoom Out Reset image size

The situation for model 5 is more complex. At the beginning of the simulation, we have ${K}_{{xx}}\gt {K}_{{yy}}\gt {K}_{{zz}}$ since the model is triaxial. The component ${K}_{{xx}}$ starts to drop whereas ${K}_{{yy}}$ rises, and they intersect at $\approx 14\,{T}_{\mathrm{simu}}$. The component ${K}_{{yy}}$ reaches its maximum value (up to $\approx 10 \%$ relative to its original value) at $\approx 25\,{T}_{\mathrm{simu}}$, and ${K}_{{xx}}$ becomes minimal (reduced by $\approx 13 \%$ relative to the original value, at $\approx 30\,{T}_{\mathrm{simu}}$). The two components meet again at $\approx 42\,{T}_{\mathrm{simu}}$. The component ${K}_{{zz}}$ does not considerably evolve and oscillates around the initial value with an amplitude less than 8%.

The right panel of Figure 19 shows the time evolution of the root mean square velocities ${v}_{\mathrm{rms}}$. For all models, the ${v}_{\mathrm{rms}}$ profiles remain stable with tiny (±4%) and smooth oscillations. Therefore, the overall kinetic energy is conserved during the simulations.

The external field generally reduces the symmetry of a given model, which could give rise to a self-rotation. Another consequence of the external field is that—unlike the case of an isolated system—the system's (local) angular momentum is typically no longer conserved. The external field introduces additional torque into the system, and hence we already expect to encounter a variation of this quantity in our simulations. To further investigate these effects, we have calculated the angular momenta for all models. The results are shown in Figure 20, which illustrates the unit-mass angular momentum components ${L}_{x}$, ${L}_{y}$, and ${L}_{z}$ in units of ${L}_{{\rm{c}}}$,

$$\begin{eqnarray}&&{L}_{x}=\displaystyle \frac{1}{M}\displaystyle \sum _{i=1}^{{N}_{p}}{m}_{i}{L}_{x}^{(i)},\end{eqnarray} \tag{ 25 }$$

and analogously for the other directions, where N_p is the number of particles inside ${r}_{90 \% }$ and ${L}_{{\rm{c}}}$ is defined as the unit-mass angular momentum with circular velocity ${v}_{{\rm{c}}}$ at a radius ${r}_{{\rm{c}}}=1\,\mathrm{kpc}$, or ${L}_{{\rm{c}}}={r}_{{\rm{c}}}{v}_{{\rm{c}}}$. The values of ${L}_{{\rm{c}}}$ are listed in Table 2. We emphasize that the angular momentum components are computed with respect to the original $(x,\,y,\,z)$ frame.

Zoom In Zoom Out Reset image size

A nonzero angular momentum is not expected for model 1 because it is an axisymmetric model and the system is initially not globally rotating. The black curves in Figure 20 confirm this. The angular momentum of model 2 is also conserved during the simulation. Although there is a tiny nonzero value of ${L}_{y}$, the values do not change noticeably within $60\,{T}_{\mathrm{simu}}$. For models 3–5, there are interesting evolutions of angular momentum. The component ${L}_{y}$ of model 3 starts increasing from $30\,{T}_{\mathrm{simu}}$ and has a small value of $0.007\,{L}_{{\rm{c}}}$, which is comparable to the constant ${L}_{y}$ of model 2. The increase of angular momentum is caused by the external field. For model 4, the evolution of angular momentum is more significant. In this case, ${L}_{y}$ changes from the beginning of the simulation and grows to $0.04\,{L}_{{\rm{c}}}$ at $60\,{T}_{\mathrm{simu}}$. The components ${L}_{x}$ and ${L}_{z}$ begin to increase from $\approx 20\,{T}_{\mathrm{simu}}$ and grow in nearly the same fashion, again since the direction of the strong external field is along the diagonal negative x–z direction. We find ${L}_{x}\approx 0.025\,{L}_{{\rm{c}}}$ and ${L}_{z}\approx 0.022\,{L}_{{\rm{c}}}$ at the end of the simulation. The components ${L}_{x}$ and ${L}_{y}$ of model 5 do not evolve since the external field lies along the x axis and there is no additional torque introduced by the external field. Nevertheless, ${L}_{z}$ of model 5 grows in the opposite direction from $\approx 20\,{T}_{\mathrm{simu}}$ and ${L}_{z}\approx -0.02$ at $60\,{T}_{\mathrm{simu}}$. The growth of angular momentum at around $20\mbox{--}30\,{T}_{\mathrm{simu}}$ for models 3–5 might be a result of model instability because of the time coincidence. However, the evolution of angular momentum is quite sophisticated and less intuitive as it depends on the assumed direction of the external field, the actual orbit composition, and the level of numerical noise, which increases with simulation time.

Finally, we point out that the observed changes in the angular momentum components are relatively small, corresponding to a few $\,\mathrm{kpc}\,\mathrm{km}\,{{\rm{s}}}^{-1}$. Consequently, any self-rotation in the models caused by the external fields must be tiny.