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

Xufen Wu, Yougang Wang, Martin Feix, Hongsheng Zhao

► To cite this version:

Xufen Wu, Yougang Wang, Martin Feix, Hongsheng Zhao. Lopsidedness of Self-consistent Galaxies Caused by the External Field Effect of Clusters. The Astrophysical Journal, 2017, 844, 10.3847/1538-4357/aa7b8a. insu-03747450

HAL Id: insu-03747450

https://insu.hal.science/insu-03747450v1

Submitted on 27 Aug 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

Xufen Wu<sup>1,2,3</sup>, Yougang Wang<sup>4</sup>, Martin Feix<sup>5,6</sup>, and HongSheng Zhao<sup>3,4</sup>

Received 2016 December 13; revised 2017 June 23; accepted 2017 June 23; published 2017 July 31

## Abstract

Adopting Schwarzschild's orbit-superposition technique, we construct a series of self-consistent galaxy models, embedded in the external field of galaxy clusters in the framework of Milgrom's MOdified Newtonian Dynamics (MOND). These models represent relatively massive ellipticals with a Hernquist radial profile at various distances from the cluster center. Using N-body simulations, we perform a first analysis of these models and their evolution. We find that self-gravitating axisymmetric density models, even under a weak external field, lose their symmetry by instability and generally evolve to triaxial configurations. A kinematic analysis suggests that the instability originates from both box and nonclassified orbits with low angular momentum. We also consider a self-consistent isolated system that is then placed in a strong external field and allowed to evolve freely. This model, just like the corresponding equilibrium model in the same external field, eventually settles to a triaxial equilibrium as well, but has a higher velocity radial anisotropy and is rounder. The presence of an external field in the MOND universe generically predicts some lopsidedness of galaxy shapes.

**Key words:** galaxies: elliptical and lenticular, cD – galaxies: kinematics and dynamics – gravitation – methods: numerical

## 1. Introduction

Elliptical galaxies widely exist in the universe, either isolated or embedded in clusters of galaxies, and they have compact centers. Studies of their dynamical behavior and evolution typically require the building of equilibrium models, which is quite challenging for elongated or triaxial systems. Schwarzschild's orbit-superposition method (Schwarzschild 1979, 1982) is a powerful technique for finding self-consistent solutions for spherical, axisymmetric, and triaxial systems. Adopting Schwarzschild's method, the aim of this work is to construct such equilibrium models for elliptical galaxies in external fields, that is, galaxies that are gravitationally bound within clusters, using the framework of Milgrom's MOdified Newtonian Dynamics (MOND; Milgrom 1983a; Bekenstein & Milgrom 1984). In addition, the kinematic properties and stability of these systems are explored with the help of N-body simulations.

The approach of Schwarzschild can be divided into three basic steps:

1. An analytic density distribution is chosen, and the corresponding gravitational potential is calculated. The whole system is then segmented into many equal-mass cells.

2. A full library of orbits within the previously calculated potential is computed, and the time spent in each cell is recorded.

3. The nonnegative linear superposition of orbits that reproduces the original density profile is determined.

The method of Schwarzschild has been applied to test various density models for self-consistency, including pattern-rotating barred galaxies (Zhao 1996; Wang et al. 2012, 2013). Many early applications of Schwarzschild's approach assumed

constant-density cores. For example, Statler (1987) found self-consistency in the perfect ellipsoid models of de Zeeuw & Lynden-Bell (1985). However, observations showed that almost all elliptical galaxies have central densities that follow a power law $\rho \sim r^{-\gamma}$ (Crane et al. 1993; Ferrarese et al. 1994; Lauer et al. 1995; Moller et al. 1995). Low-luminosity ellipticals have steeper centers, $\gamma \approx 2$, while the most luminous ellipticals have shallower ones, $0 \le \gamma \le 1$. Also, Tremblay & Merritt (1996) found that the intrinsic shapes of ellipticals depend on their luminosity: the short-long axis ratio of the most luminous ellipticals has a peak at 0.75, whereas that of low-luminosity ellipticals peaks at 0.65. Dehnen (1993) proposed a family of models whose density distributions follow $\rho \sim r^{-\gamma}$ in the central region and $\rho \sim r^{-4}$ at distant radii, where $\gamma$ is a free parameter. Various studies (Merritt & Fridman 1996; Rix et al. 1997; Poon & Merritt 2004; Capuzzo-Dolcetta et al. 2007) of Schwarzschild's technique applied to triaxial Dehnen profiles have been conducted, and it has been shown that self-consistent solutions for these models can be constructed in the case of Newtonian gravity.

As an alternative to cold dark matter on galactic scales, the MOND paradigm is built on the tight relation between the distribution of baryons and the gravitational field in spiral galaxies (Famaey et al. 2007b; McGaugh et al. 2007). In fact, its simple formulation leads to excellent predictions of rotation curves for galaxies ranging over five decades in mass (see, e.g., Sanders & McGaugh 2002; Bekenstein 2006), including our own Milky Way (Famaey & Binney 2005; Famaey et al. 2007a). The successes and problems of MOND are extensively discussed in Famaey & McGaugh (2012). Moreover, MOND has recently been used to explain the rotational speed in polar rings (Lüghausen et al. 2013), the formation of the shell

<footnote n="1"><sup>1</sup> CAS Key Laboratory for Research in Galaxies and Cosmology, Department of Astronomy, University of Science and Technology of China, Hefei, 230026, P.R. China; xufenwu@ustc.edu.cn</footnote>

<footnote n="2"><sup>2</sup> School of Astronomy and Space Science, University of Science and Technology of China, Hefei 230026, P.R. China</footnote>

<footnote n="3"><sup>3</sup> School of Physics and Astronomy, University of St Andrews, North Haugh, Fife, KY16 9SS, UK</footnote>

<footnote n="4"><sup>4</sup> Key Laboratory of Computational Astrophysics, National Astronomical Observatories, Chinese Academy of Sciences, Beijing, 100012, P.R. China</footnote>

<footnote n="5"><sup>5</sup> CNRS, UMR 7095 & UPMC, Institut d'Astrophysique de Paris, 98 bis Boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="6"><sup>6</sup> Department of Physics, Israel Institute of Technology—Technion, Haifa 32000, Israel</footnote>
structure in the elliptical galaxy NGC 3923 (Bílek et al. 2013, 2014), the velocity dispersion of Andromeda dwarf galaxies (McGaugh & Milgrom 2013), and the mass discrepancy-acceleration correlation of disk galaxies (Milgrom 1983b; Sanders 1990; McGaugh 2004; Wu & Kroupa 2015) and of pressure-supported systems (Scarpa 2006). It also provides constraints on the mass-to-light ratio derived from the vertical stellar velocity dispersion (Angus et al. 2016).

In contrast to the Newtonian case, the internal dynamics of a gravitating system in MOND is affected by any external fields; that is, even a freely falling system in MOND will exhibit a dynamical evolution different from that of an isolated one. This attribute implies a violation of the strong equivalence principle and is usually referred to as the external field effect (Milgrom 1983b). The impact of external fields has been studied for a variety of different situations, including the motion of probes in the inner solar system (Milgrom 2009), the Roche lobe of binary stars (Zhao & Tian 2006), the kinetics of stars in globular clusters (Milgrom 1983b; Perets et al. 2009), the escape speeds and truncations of galactic rotation curves (Wu et al. 2007), satellites surrounding a host galaxy (Brada & Milgrom 2000; Tiret et al. 2007; Angus 2008), and the phase transition of distant star clusters moving toward the galactic center (Wu & Kroupa 2013).

Schwarzschild's orbit-superposition method has already been applied within Milgromian dynamics: for example, models of ellipsoidal field galaxies were found to be both self-consistent and stable (Wang et al. 2008; Wu et al. 2009). Further, the morphology of elliptical cluster galaxies was discussed by Wu et al. (2010) and exhibits lopsided shapes along the external field direction. In what follows, we will mainly focus on the kinematic aspects of these models and the numerical details of our approach, extending the analysis to triaxial systems in a gravitating environment. Although it yields a less realistic scenario, we only consider external fields that are uniform and constant. Independent of the framework, external tidal fields will influence a system's internal dynamics and may obscure fundamental differences between gravity theories. To maximally distinguish between Newtonian dynamics and MOND, we thus ignore tidal effects in our analysis.<fnref n="7" /> Such an idealized case corresponds to systems either entirely restricted within the tidal radius or moving in a smooth and slowly varying background field, for example, a galaxy circularly orbiting the cluster center. In these situations, an external field is mainly dominated by its uniform part, and tidal effects play only a subordinate role.

Already posing a challenge in Newtonian gravity, the construction of equilibrium models for elongated or triaxial systems in MOND is additionally hampered by the nonlinear modification of Poisson's equation. This is especially true for systems with external fields because their phase-space distribution is determined by both internal density and external field. In their galaxy merger simulations, Nipoti et al. (2007) obtained the distribution function of a Hernquist sphere by Eddington inversion in a MONDian potential. For more complex triaxial models, however, it is generally not possible to obtain analytic solutions to Eddington's equation, and this procedure becomes very difficult. Tiret & Combes

(2007, 2008) employed Newtonian equilibrium models for spiral galaxies embedded into a Plummer-type dark halo, and they replaced gravity by MONDian dynamics afterwards. Similarly, the simulations by Haghi et al. (2009) made use of Newtonian models initially, but then particle velocities were increased to avoid gravitational collapse of globular cluster models in the Milky Way. If set up in this way, such initial conditions (ICs) will immediately cause the system in question to relax until it reaches a new state of equilibrium. This will remain true when using Schwarzschild's approach, which was designed for genuine equilibrium systems. For example, initially axisymmetric models will start developing asymmetric shapes (Wu et al. 2010). Below we shall investigate this relaxation process in more detail, including both axisymmetric and triaxial configurations.

The stability of disk galaxies hosted by dark matter halos in Newtonian gravity has been studied by Sellwood & Evans (2001). It has been shown that the lopsided instability ($m=1$ mode) can be avoided when their massive outer disks are tapered and the galaxies are stabilized by dense centers. Further, De Rijcke & Debattista (2004) found off-center nuclei in flattened nonrotating systems, and a promising mechanism is the destruction of box orbits. The growth of two different $m=1$ modes, associated with a Jeans-type instability for counterrotating disk models and swing amplification for fully rotating disk models, was studied by Dury et al. (2008). As most observed galaxies appear lopsided ($m=1$ perturbations; Rix & Zaritsky 1995; Haynes et al. 1998; de Zeeuw et al. 2002), this motivates investigating such mechanisms also in the context of MOND. From an analytic point of view, there is still little known about the stability of galaxies in MOND (Brada & Milgrom 1999; Wu et al. 2009; Nipoti et al. 2011), and, in particular, no stability studies exist in external or tidal fields. A numerical study thus provides a first step in this direction.

The paper is organized as follows. In Section 2, we introduce our basic setup and discuss the effect of external fields on static galaxy potentials. In Section 3, we use Schwarzschild's technique to construct quasi-equilibrium models of galaxies, and we perform a kinematic analysis of these systems in Section 4. Finally, we compare our results to the evolution of isolated galaxies and conclude in Section 8.

## 2. Mass Models and Static Potentials

### 2.1. Density Profiles and External Fields

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

$$ \nabla \cdot \left[ \mu \left( \frac{|g|}{a_0} \right) g \right] = 4\pi G\rho_b, \quad g = \nabla \Phi_{\text{int}} + g_{\text{ext}}. \quad (1) $$

Here $\Phi_{\text{int}}$ is the internal potential generated by the baryon density $\rho_b$, and $a_0 \approx 3700 \text{ km}^2 \text{ s}^{-2} \text{ 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 = g/a_0$):

$$ \begin{aligned} \mu(x) &\sim x & x &\ll 1, \\ \mu(x) &\sim 1 & x &\gg 1. \end{aligned} \quad (2) $$

<footnote n="7">⁷ Generally, it is desirable to have a full treatment of the problem including tidal effects. This would allow us to explore other more complex scenarios such as evolution in fast-varying backgrounds, e.g., a galaxy crossing the center of cluster, and will be subject to future work.</footnote>
<table><caption>**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*<sub>ext</sub> Denotes the Strength of the External Field in Units of $a_0 \approx 3700(\text{km s}^{-1})^2$ kpc</caption><thead><tr><th>Model</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td><i>M</i> [10<sup>10</sup> <i>M</i><sub>☉</sub>]</td><td>5</td><td>5</td><td>5</td><td>5</td><td>10</td></tr><tr><td><i>a</i>: <i>b</i>: <i>c</i></td><td>1: 1: 0.7</td><td>1: 1: 0.7</td><td>1: 1: 0.7</td><td>1: 1: 0.7</td><td>1: 0.86: 0.7</td></tr><tr><td><i>g</i><sub>ext</sub> [a<sub>0</sub>]</td><td>0</td><td>0.01</td><td>0.1</td><td>1.0</td><td>1.0</td></tr><tr><td>Direction of <i>g</i><sub>ext</sub></td><td>...</td><td>positive x-z<br>diagonal</td><td>positive x-z<br>diagonal</td><td>negative x-z<br>diagonal</td><td>negative<br>x-axis</td></tr><tr><td>Density symmetry</td><td>axisym</td><td>axisym</td><td>axisym</td><td>axisym</td><td>triaxial</td></tr><tr><td>Potential sym. axes</td><td><i>x</i>, <i>y</i>, <i>z</i></td><td><i>y</i></td><td><i>y</i></td><td><i>y</i></td><td><i>y</i>, <i>z</i></td></tr><tr><td>Starting octants</td><td><i>y</i></td><td>I, II, V, VI</td><td>I, II, V, VI</td><td>I, II, V, VI</td><td>I, II</td></tr><tr><td>Reflecting planes</td><td><i>x</i>-<i>y</i>, <i>x</i>-<i>z</i>, <i>y</i>-<i>z</i></td><td><i>x</i>-<i>z</i></td><td><i>x</i>-<i>z</i></td><td><i>x</i>-<i>z</i></td><td><i>x</i>-<i>y</i>, <i>x</i>-<i>z</i></td></tr><tr><td><i>N</i><sub>cell</sub></td><td>960</td><td>3840</td><td>3840</td><td>3840</td><td>1920</td></tr><tr><td><i>N</i><sub>stationary</sub></td><td>3840</td><td>15360</td><td>15360</td><td>15360</td><td>7680</td></tr><tr><td><i>N</i><sub>ejecting</sub></td><td>3000</td><td>12000</td><td>12000</td><td>12000</td><td>6000</td></tr><tr><td><i>N</i><sub>orbit</sub></td><td>6840</td><td>27360</td><td>27360</td><td>27360</td><td>13680</td></tr></tbody></table>

**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_{cell}$ is the number of cells to impose self-consistency, $N_{stationary}$ is the number of stationary starting orbits built into the orbital library, and $N_{ejecting}$ is the number of ejecting orbits starting from the x-z plane. The total number of orbits is given by $N_{orbit}$.

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

$$ \mu(x) = \frac{x}{1 + x}. \qquad (3) $$

Although observationally excluded in the solar system, this $\mu$ 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):

$$ \rho(r) = \frac{M}{2\pi abc} \frac{1}{r(1 + r)^3}, \qquad (4) $$

where

$$ r = \sqrt{\left(\frac{x}{a}\right)^2 + \left(\frac{y}{b}\right)^2 + \left(\frac{z}{c}\right)^2}, \qquad (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}-10^{11} M_{\odot}$, bright enough to observe the outer parts. In this case, the internal accelerations are comparable to several $a_0$. 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}-10^{11} M_{\odot}$ range from 0.5 to 3 kpc. For all five models, we set *a* = 1 kpc, which lies within the $R_e$ dispersion range of observational data points (Gadotti 2009).

The strengths of the external fields are chosen as $g_{ext} = 0.01a_0$, $0.1a_0$, and $1.0a_0$ for the axisymmetric models and

$1.0a_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 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.

## 2.2. Distortion of Static Potentials in External Fields

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{aligned} r_i &= 2 \tan[(i + 0.5)0.5\pi/(n_r + 1)] \text{kpc}, \\ \theta_j &= \pi \times (j + 0.5)/n_{\theta}, \\ \phi_k &= 2\pi \times k/n_{\phi}, \end{aligned} \qquad (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 $|g_{int}| \ll |g_{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$ 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.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Left panel: isodensity (dotted) and isopotential (solid) contours of model 4. The isodensity contours correspond to ellipsoidal radii of 2, 5, 10, 15, and 20 kpc. The isopotential contours are located at the same radii along the major axis. The direction of the external field is along the negative x-z diagonal. Right panel: the ratio of isopotentials, ($r_{\perp}: r_{\parallel})_{\Phi(r_{\perp})}$, evaluated at the same radius in antiparallel and parallel directions of the external field for model 4.</figcaption></figure>

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.

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):

$$ \mu \left( \frac{|g|}{a_0} \right) g = g_N + \nabla \times h, \quad (7) $$

where $g_N$ is the Newtonian gravitational field, and $\nabla \times h$ is a solenoidal vector field determined by the condition that $g$ can be expressed as the gradient of a scalar potential. Restricting ourselves to the axis parallel to the external field's direction, $g$ and $g_N$ must be either parallel or antiparallel, assuming that the symmetry center coincides with the coordinate origin. Hence the curl term $\nabla \times h$ vanishes. The strength of the total gravitational acceleration $g = |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_{\text{ext}} - g_{\text{int}}|$, while on the positive semi-axis it is $g^+ = g_{\text{int}} + g_{\text{ext}}$. The two different sides have different values of the $\mu$ 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_{\text{ext}} + g_{\text{int}}|$ and steeper in the third octant where $g^- = g_{\text{int}} - g_{\text{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 < 3$ kpc), the internal gravitational field dominates, and thus $r^-: r^+$ is basically 1. At much larger radii ($r \gg 10$ kpc), the relative contribution of the external field increases, and $r^-: r^+$ falls down to 1 again because the $\mu$ function approaches a constant, $\mu = \mu(|g_{\text{ext}}|/a_0)$.

### 3. Schwarzschild Technique and Model Self-consistency

Schwarzschild (1979, 1982) proposed the orbit-superposition method to reproduce the density distribution of galaxies and build triaxial galaxy models. The basic idea is to compute a large library of orbits in a given potential and determine the superposition of orbits that provides the best fit to the observational density distribution or the underlying density model. Let $N_{\text{orbit}}$ be the number of the orbits in the library ($j \in N_{\text{orbit}}$) and $N_{\text{cell}}$ the total number of grid cells segmenting space ($i \in N_{\text{cell}}$). Further, let $O_{ij}$ denote the fraction of time spent by the $j$th orbit in the $i$th cell. The weight and mass of the $j$th orbit are defined by $w_j$ and $m_i$, respectively, and they are related by the following set of linear equations:

$$ \sum_{j=1}^{N_{\text{orbit}}} w_j O_{ij} = m_i. \quad (8) $$

Schwarzschild's method is widely used to build spherical, axisymmetric, and triaxial models for galaxies (Richstone 1980, 1984; Pfenniger 1984; Richstone & Tremaine 1984; Zhao 1996; Rix et al. 1997; van der Marel et al. 1998; Binney 2005; Capuzzo-Dolcetta et al. 2007; Wu et al. 2009).

The $O_{ij}$ array is obtained by computing the superpositions of the $j$th orbit in equal time intervals $\Delta\tau_j$, and counting the number of the output points $\nu_{ij}$ in the $i$th cell. After that, the elements are determined according to

$$ O_{ij} = \frac{\Delta\tau_j \times \nu_{ij}}{\Delta\tau_j \times \nu_j} = \frac{\nu_{ij}}{\nu_j}, \quad (9) $$

where $\nu_j$ is the total output number of the $j$th orbit. The previous analyses in Wang et al. (2008) and Wu et al. (2009) used nonequal time interval outputs, given by the variation of the gravitational field strength, $\Delta\tau'_j \sim 1/|\nabla \cdot g|^2$. In this case, the real time intervals $\Delta\tau'_j$ of the $j$th orbit are not constant anymore. However, in these previous studies, the unevenness of the time intervals between the outputted points along an orbit was neglected, and the $O_{ij}$ were calculated using $O_{ij} \approx \nu_{ij}/\nu_j$, which systematically increases the number of output points in cuspy centers. Our present analysis does not make this approximation. Each orbit is integrated for 100 times its circular orbital time, 100 $T_{\text{cir}}$, so the equal time interval $O_{ij}$ does not give rise to additional inaccuracy. We will discuss the time integration of the orbits in Section 3.1. Further details related to grid segmentation, ICs, and orbital classification are separately given in Appendix A.

#### 3.1. Integrating the Orbits

As stated in the previous section, all orbits are integrated for 100 $T_{\text{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
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The circular orbital time in units of Gyr plotted against radius in units of kpc, using a logarithmic scaling. The black, magenta, cyan, yellow, and green lines represent models 1–5 listed in Table 1, respectively.</figcaption></figure>

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_c = (GMa_0)^{1/4}$, and hence $T_{cir} \propto r$. In the case of strong external fields, however, the interpolating function $\mu$ becomes a constant far away from the system, leading to a Newtonian-like behavior $T_{cir} \propto r^{3/2}$.

## 3.2. Smoothed Solutions

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 ith 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:

$$ \chi^2 = \frac{1}{N_{\text{cell}}} \sum_{i=1}^{N_{\text{cell}}} \left( m_i - \sum_{j=1}^{N_{\text{orbit}}} w_j O_{ij} \right)^2 . \quad (10) $$

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

$$ \delta = \sqrt{\chi^2 / \bar{m}}, \quad (11) $$

where $\bar{m}$ is the mean Monte Carlo mass in cells. For self-consistent models, the value of $\delta$ 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 $\delta$ values are on the order of $10^{-2}$. 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 $\delta$ values grow with increasing external field strength.

To construct sufficiently self-consistent models, one needs to use a large number of orbits, often far more than the number of cells $N_{cell}$. The best solutions are typically *nonunique* with a very noisy phase-space distribution characterized by $N_{orbit} - N_{cell}$ and $N_{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 $\lambda$, where $\lambda = N_{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

$$ \chi^2_{\text{smooth}} = \frac{1}{N_{\text{cell}}} \sum_{i=1}^{N_{\text{cell}}} \left( m_i - \sum_{j=1}^{N_{\text{orbit}}} w_j O_{ij} \right)^2 + \lambda \sum_{j=1}^{N_{\text{orbit}}} w_j^2. \quad (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_{smooth} = \sqrt{\chi^2_{smooth}/\bar{m}}$ (where $\bar{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}-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 $<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
<caption>Table 2</caption>

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

<table><thead><tr><th>Model</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>$N_{\text{orbit}}$</td><td>6840</td><td>27360</td><td>27360</td><td>27360</td><td>13680</td></tr><tr><td>$\delta$</td><td>1.5 × 10<sup>-15</sup></td><td>1.37178 × 10<sup>-2</sup></td><td>2.44238 × 10<sup>-2</sup></td><td>4.41140 × 10<sup>-2</sup></td><td>6.2 × 10<sup>-5</sup></td></tr><tr><td>$\delta_{\text{smooth}}$</td><td>6.7 × 10<sup>-5</sup></td><td>1.37180 × 10<sup>-2</sup></td><td>2.44239 × 10<sup>-2</sup></td><td>4.41141 × 10<sup>-2</sup></td><td>6.3 × 10<sup>-5</sup></td></tr><tr><td>$N_{>0}$</td><td>960</td><td>2378</td><td>2102</td><td>1677</td><td>1905</td></tr><tr><td>$N_{>0,\text{smooth}}$</td><td>5995</td><td>2379</td><td>2103</td><td>1677</td><td>1949</td></tr><tr><td>$v_{\text{rms}}$ [km s<sup>-1</sup>]</td><td>225.55</td><td>225.20</td><td>224.67</td><td>220.28</td><td>309.58</td></tr><tr><td>$-2K/W$</td><td>1.01</td><td>1.02</td><td>1.01</td><td>1.00</td><td>1.01</td></tr><tr><td>$\xi$</td><td>1.26</td><td>1.74</td><td>1.57</td><td>1.65</td><td>1.60</td></tr><tr><td>$L_c$ [km s<sup>-1</sup> kpc]</td><td>253.25</td><td>251.59</td><td>252.16</td><td>253.99</td><td>352.10</td></tr></tbody></table>

Note. The radial-to-tangential anisotropic ratio $\xi$ characterizes the radial instability. Here, $L_c$ denotes the unit-mass angular momentum with circular velocity $v_c$ at radius $r_c = 1$ kpc, where $L_c = r_c v_c$ is expressed in units of km s$^{-1}$ kpc.

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.

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. Instability of Cluster Galaxies

It is unknown whether quasi-equilibrium models constructed with Schwarzschild's approach are stable. The direct way to test the stability and evolution is to use N-body tools. Due to the external field, the potentials of axisymmetric density profiles are lopsided, and orbits running in these potentials also become lopsided. For an arbitrary orbit integrated in a given potential for a long enough time, the mass reproduced by this orbit will also be lopsided. Thus the uncertainty on the

model's stability increases in this case. It is an important issue to investigate the stability of MOND models in external fields since there are many elliptical galaxies observed in clusters. In what follows, we want to take a first step in this direction by performing a kinematic analysis of the previously introduced models, starting with N-body ICs given by Schwarzschild's approach.

### 4.1. ICs and Numerical Setup

#### 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 jth orbit is $n_j = w_j N$, where $w_j$ is the weight of the jth orbit. Particles are placed on the jth orbit on equidistant times given by the interval $\Delta t_j = T_j/n_j$, where $T_j$ is the total integrated time for the jth orbit (see Figure 2).<fnref n="8" />

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^6$ 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

$$W = \int \rho x \cdot \nabla \Phi d^3 x \quad (13)$$

<footnote n="8">⁸ One can also randomly sample particles on the jth orbit from a uniform distribution. Since most of the orbits have small positive weights in our simulations, the number of particles on the jth are quite small. A random sampling might introduce numerical noise $\propto 1/\sqrt{n_j}$ and could, therefore, have problems reflecting the real shape of the orbit if the weight is small. To avoid such problems, we choose an isochronous sampling.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Left panels: the integrated contributions of different orbit families (for energies <E) to the mass as a function of energy for models 1–5 (from top to bottom panels), assuming Equation (10), that is, no smoothing. Right panels: the same as the left panels, but now assuming the regularization given by Equation (12).</figcaption></figure>

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 ± 0.02. The vrms 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_{\text{rms}} \approx \sigma$ for pressure-supported systems becomes smaller. For the strongest external field (model 4), the $v_{\text{rms}}$ value is smallest, $\approx 5$ km s$^{-1}$ less than in model 1. Even in this case, the decrease of $v_{\text{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_{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 $\sigma$ in this galaxy is only 2.1 km 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

As discussed in Section 3, large populations of box orbits are selected to fit the underlying density distribution. In addition, the nonclassified orbits are characterized by low angular momentum and thus are highly radially anisotropic. It is therefore important to examine the radial instability of the model ICs. To this end, we consider the anisotropy parameter $\xi = \frac{2K_r}{K_t}$, which has been introduced for spherical Newtonian systems (Polyachenko & Shukhman 1981; Saha 1992; Bertin et al. 1994; Trenti & Bertin 2006). Here $K_r$ and $K_t = K_\theta + K_\phi$ are the radial and tangential kinetic energy components, $K_r = \sum_{ip=1}^N \frac{1}{2} m_{ip} \sigma_{r,ip}^2$, where the index ip runs over all particles, and $K_\theta$ and $K_\phi$ can be defined in the same manner. Spherical Newtonian systems are unstable if $\xi$ lies above a critical value, $\xi_{crit}$. Various studies have found different values of $\xi_{crit}$ (e.g., May & Binney 1986, $\xi_{crit} \approx 2.2$; Saha 1991, $\xi_{crit} \approx 1.4$; Saha 1992, $\xi_{crit} \approx 2.3$; Bertin & Stiavelli 1989, $\xi_{crit} \approx 1.9$; Bertin et al. 1994, $\xi_{crit} \approx 1.6$). Spherical models are generally unstable when $\xi \gtrsim 2.3$, but the models may also be unstable for smaller values of $\xi$ if their distribution functions increase rapidly with low angular momentum. The radial instability transforms an originally spherical system into a triaxial one and also alters the spatial distribution of the velocity dispersion. The $\sigma(r)$ profiles become more isotropic in the center and radially anisotropic in the outer regions (Barnes et al. 2005; Bellovary et al. 2008). Moreover, axisymmetric models are unstable within a major part of the parameter space (Levison et al. 1990). For triaxial models, a collective radial instability has been studied by Antonini et al. (2008). Such instabilities are caused by the box-like orbits, rather than nonclassified ones or by any deviations in the model self-consistency. The radial instability causes triaxial models to become more prolate (Antonini et al. 2008, 2009).

In the context of MOND, the parameter $\xi$ has been studied for spherical Osipkov-Merritt radially anisotropic $\gamma$ models (with $\gamma = 0$ and 1, where the latter recovers the Hernquist model; Nipoti et al. 2011). It was found that a MOND system with radial anisotropy is more stable than a pure Newtonian model with exactly the same density distribution. As the anisotropic radius in MOND is larger than that in a pure Newtonian model, a larger fraction of radial orbits can exist in the outer region of a MOND system. The inferred $\xi_{crit}$ values for MOND systems appear within the range [2.3, 2.6]. Considering our galaxy models, we estimate $1.2 < \xi < 1.8$, which is well below the corresponding values of $\xi_{crit}$. Since the external field breaks the potential symmetries associated with these models, however, one cannot make a conclusive

<figure><img src="image_6.png" /><figcaption>Figure 4. Upper panel: the static potential of model 3 along the x (solid curves), y (dotted), and z (dashed) axes for the Schwarzschild modeling ($\Phi_{Schwarzschild}$, black curves) and for the N-body simulations ($\Phi_{N-body}$, magenta curves). Lower panel: relative deviation between the differently resolved potentials, $\Delta\Phi = (\Phi_{Schwarzschild} - \Phi_{N-body})/\Phi_{Schwarzschild}$.</figcaption></figure>

statement about the presence of radial instabilities in these cases. While a stability study for isolated triaxial MOND models has been carried out by Wu et al. (2009), an analysis of MOND models in external fields is still missing. In the following sections, we will present a first investigation on the stability of such 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

$$r_i = 8 \tan[(i + 0.5)0.5\pi/(n_r + 1)] \text{kpc}, \quad (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
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Evolution of the virial ratio: the black, magenta, cyan, yellow, and green lines represent models 1–5, respectively. Here $2\pi(GM/a^3)^{-1/2}$ is the Keplerian rotation time at the length scale $a = 1$ kpc for a total mass $M$.</figcaption></figure>

has previously been studied in Wu et al. (2008). The time unit used in the NMODY code is given by (Wu et al. 2009)

$$ T_{\text{simu}} = \left(\frac{GM}{a^3}\right)^{-1/2} = 4.7 \times 10^6 \text{ year} \left(\frac{M}{10^{10} M_{\odot}}\right)^{-1/2} \left(\frac{a}{1 \text{ kpc}}\right)^{3/2} . \quad (15) $$

Here the quantity $2\pi T_{\text{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}$.

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_{\text{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_{\text{simu}}$ (twice the value of $60 T_{\text{simu}}$) to examine the systems' behavior beyond the actual simulation time interval. In what follows, we will restrict the discussion to within $60 T_{\text{simu}}$; only virial ratios are presented for the fully simulated range of $120 T_{\text{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_{\text{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_{\text{simu}}$.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Left panels: the evolution of Lagrangian radii for models 1–5 (from top to bottom, respectively). Right panels: density distribution along the galaxy major axes. The analytic initial density and the density constructed from $N$ particles correspond to the black solid and dashed lines, respectively. The violet, blue, yellow, and green lines correspond to simulation times of 15, 30, 45, and $60 T_{\text{simu}}$, respectively.</figcaption></figure>

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.

## 4.2. Mass Distribution

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.

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{aligned}
\tilde{I}_{xx}(r) &= \frac{\sum_i m_i (y_i^2 + z_i^2) / r_i^2}{\sum_i m_i}, \\
\tilde{I}_{xy}(r) &= \frac{\sum_i -m_i x_i y_i / r_i^2}{\sum_i m_i},
\end{aligned}
\quad (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'x'}(r)$, $\tilde{I}_{y'y'}(r)$, and $\tilde{I}_{z'z'}(r)$, where the primed coordinate system refers to the corresponding eigenframe. The associated principal axes,

$$ 
\begin{aligned}
\tilde{a}(r) &= \sqrt{[\tilde{I}_{y'y'}(r) + \tilde{I}_{z'z'}(r) - \tilde{I}_{x'x'}(r)]/2}, \\
\tilde{b}(r) &= \sqrt{[\tilde{I}_{x'x'}(r) + \tilde{I}_{z'z'}(r) - \tilde{I}_{y'y'}(r)]/2}, \\
\tilde{c}(r) &= \sqrt{[\tilde{I}_{x'x'}(r) + \tilde{I}_{y'y'}(r) - \tilde{I}_{z'z'}(r)]/2},
\end{aligned}
\quad (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^{-3}$.

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}$.<fnref n="9" /> 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–30 $T_{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]$ 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 < 5.0$ 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.

### 4.3. Axial Ratios

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_{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.

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_{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 $\pm 0.05$ within 60 $T_{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_{simu}$. For model 4, the values of $\tilde{p}_{90}$ decrease to around 0.9 within the first 20 $T_{simu}$ and stay constant afterwards. The ratio $\tilde{q}_{90}$ oscillates around 0.75 within 20 $T_{simu}$, and then jumps to 0.82 at 25 $T_{simu}$. There is no significant evolution of $\tilde{q}_{90}$ at later times. The different evolution of axial ratios between models 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_{simu}$ for models 2 and 3 and $\sim 25 T_{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_{simu}$, but then decrease until 25 $T_{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_{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_{simu}$ for models 4 and 5, after which they become stable. We also note that the

<footnote n="9">⁹ The oscillations of the radial densities emerge from numerical noise. The radial densities along the major axes are computed from the N-body simulation grids that are closest to the major axes. Approximately, there are $N_P/(n_\theta \times n_\phi) \approx 250$ particles along the major axes, and Equation (14) is used for radial binning. The mth radial bin contains $n_{p,m}$ particles, and the associated numerical noise is $\approx 1/\sqrt{n_{p,m}}$.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Left panel: time evolution of axial ratios at radii enclosing 90% of the total mass. The colors black, magenta, cyan, yellow, and green represent the models 1–5, respectively, and illustrate the axial ratios $\bar{p}_{90\%} = \bar{b}_{90\%}/\bar{a}_{90\%}$ (solid lines) and $q_{90\%} = \bar{c}_{90\%}/\bar{a}_{90\%}$ (dashed lines). Right panel: axial ratios of the five models as a function of enclosed masses. The upper set of curves refers to $p = \bar{b}/\bar{a}$, and the lower set to $q = \bar{c}/\bar{a}$. The solid and dotted curves represent initial ($T=0$) and final ($T=60 T_{\text{simu}}$) ratios, respectively. The colors are the same as in the left panel.</figcaption></figure>

evolution from models 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 = \bar{b}/\bar{a}$ and $q = \bar{c}/\bar{a}$ in model 5, with amplitudes <±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 $\delta$ 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_{\text{simu}}$ (dotted curves), the values of $p$ for models 1–3 drop to ≈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 <+0.02/-0.05. 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.

### 4.4. Kinetics

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

$$ \beta(r) = 1 - \frac{\sigma_{\theta}^{2} + \sigma_{\phi}^{2}}{2\sigma_{r}^{2}}, \qquad (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)$.<fnref n="10" /> 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–30$T_{simu}$. In addition, we find that $\sigma_r(r)$ is flattened in the inner parts (<1.0 kpc) of all systems, with $\sigma_r(r < 1.0 \text{ kpc}) \approx 170 \text{ km 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_{simu}$. At larger radii, dispersions $\sigma_r(r > 1 \text{ kpc})$ are about

<figure><img src="image_10.png" /><figcaption>Figure 8. Left panels: time evolution of the radial velocity dispersion. The panels (from top to bottom) show the results for models 1–5. Right panels: the anisotropy parameter $\beta$. The black lines represent the ICs, and the different colors are defined as in Figure 6.</figcaption></figure>

10–20 km s$^{-1}$ higher than in the initial state. The inner region of $\sigma_r(r)$ for model 5 decreases from 300 to 250 km s$^{-1}$ within 30$T_{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.

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)$

<footnote n="10">¹⁰ The radial velocity dispersions are plotted in 21 radial bins of equal mass (see Appendix A). Due to the increased particle number in these bins, the numerical noise is much smaller than in the right panels of Figure 6. The same binning procedure is applied for other quantities that are discussed hereafter.</footnote>
<figure><img src="image_11.png" /><figcaption>Figure 9. The anisotropy profile of model 1.</figcaption></figure>

profiles are no longer constant. In general, $\beta(r)$ declines in the inner region and starts to increase from $\approx 0.5$ kpc to $\approx 4$ kpc. The model with the weakest external field has the highest radial anisotropy at 4 kpc, $\beta(r) \ge 5$ kpc $\approx 0.8$. The $\beta(r)$ profile stays almost constant in the outer region, where $r > 5$ 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 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 1, and the result is shown in Figure 9. We find that $\beta(r)$ for model 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$.

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

<table><caption>Table 3
Parameters of the New Orbital Library after Removing the Box Orbits</caption><thead><tr><th>Model</th><th>1'</th><th>2'</th><th>3'</th><th>4'</th></tr></thead><tbody><tr><td>$N'_{\text{orbit}}$</td><td>4178</td><td>19361</td><td>21424</td><td>23372</td></tr><tr><td>$\delta_{\text{smooth}}$</td><td>$1.52 \times 10^{-2}$</td><td>$5.33 \times 10^{-2}$</td><td>$4.32 \times 10^{-2}$</td><td>$4.90 \times 10^{-2}$</td></tr><tr><td>$N'_{\ge 0, \text{smooth}}$</td><td>617</td><td>1517</td><td>1601</td><td>1465</td></tr><tr><td>$-2K/W$</td><td>1.01</td><td>1.01</td><td>1.01</td><td>1.00</td></tr><tr><td>$\xi$</td><td>0.98</td><td>1.58</td><td>1.79</td><td>1.65</td></tr></tbody></table>

Note. The number of total orbits is denoted by $N'_{\text{orbit}}$, $\delta_{\text{smooth}}$ is the self-consistency parameter, and $N'_{\ge 0, \text{smooth}}$ 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.

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 > 1.0$ kpc, but decreases in the inner region within 15–30 $T_{\text{simu}}$, after which it remains stable. The models with the strongest external field show the most prominent evolution within 1 kpc. At 60 $T_{\text{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.

# 5. The Origin of the Instability

## 5.1. Box Orbits Removed

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'_{\text{orbit}}$, 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_{\text{smooth}}$ parameters are on the order of $10^{-2}$ 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
<figure><img src="image_12.png" /><figcaption>Figure 10. The orbital structure of the new models 1′–4′ selected by Equation (12) after removing box orbits from their orbital library.</figcaption></figure>

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'.

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 $\xi$ values of the new models (see Table 3). We find that $\xi \ll \xi_{crit}$ for model 1' while $\xi \in [1.58, 1.79]$ for the other ones. These values lie within the $\xi$ 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–60$ $T_{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_{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_{simu}$. The values of *p* and *q* slowly decline with increasing

<figure><img src="image_13.png" /><figcaption>Figure 11. Left panel: time evolution of axial ratios at $r_{90\%}$ for models 1' (black), 2' (magenta), 3' (cyan), and 4' (yellow). The solid, dotted, and dashed lines denote the different components $\tilde{a}$, $\tilde{b}$, and $\tilde{c}$, respectively. Right panel: axial ratios of the four models as a function of enclosed mass. The upper set of curves refers to $p = \tilde{b}/\tilde{a}$ and the lower set to $q = \tilde{c}/\tilde{a}$. The solid, dotted, and dashed curves represent initial ($T=0$) and final ($T = 60$ $T_{simu}$) axial ratios.</figcaption></figure>

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_{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.

## 5.2. Box Orbits and Nonclassified Orbits Removed

Nonclassified orbits are defined as low angular momentum orbits in which the velocity cannot be restricted to a rectangular
<figure><img src="image_14.png" /><figcaption>**Figure 12.** The upper panels illustrate the orbital structure of models 2''–4'' computed with Equation (12) after removing box and nonclassified orbits from the orbital library. The lower panels show the orbital structures of models 2*–4*, which have been obtained in the same way after removing the box and long-axis loop orbits.</figcaption></figure>

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 $\delta$ 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_{simu}$ are illustrated in the upper panels of Figure 13. We find that the models are still axisymmetric after 60 $T_{simu}$. No triaxial configurations emerge during the N-body experiments, which agrees with the small $\xi$ values estimated for these models. Ignoring the innermost particles within $r_{20\%}$, the axial ratios of model 2'' at 60 $T_{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.

<table><caption>**Table 4**
Parameters of the New Orbital Libraries after Removing Box and Nonclassified Orbits</caption><thead><tr><th>Model</th><th>2"</th><th>3"</th><th>4"</th></tr></thead><tbody><tr><td>$N''_{orbit}$</td><td>12000</td><td>11979</td><td>11297</td></tr><tr><td>$\delta''_{smooth}$</td><td>0.087</td><td>0.087</td><td>0.108</td></tr><tr><td>$N''_{>0,smooth}$</td><td>6077</td><td>5274</td><td>820</td></tr><tr><td>$-2K/W$</td><td>1.01</td><td>1.01</td><td>1.01</td></tr><tr><td>$\xi$</td><td>0.98</td><td>0.83</td><td>0.98</td></tr></tbody></table>

**Note.** The number of total orbits is denoted by $N''_{orbit}$, $\delta''_{smooth}$ is the self-consistency parameter, and $N''_{>0,smooth}$ 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.

### 5.3. Box Orbits and Long-axis Loop Orbits Removed

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
<figure><img src="image_15.png" /><figcaption>**Figure 13.** The upper panels show the time evolution of axial ratios at $r_{90\%}$ (left) and axial ratios as a function of enclosed masses (right) for models 2'' (magenta), 3'' (cyan), and 4'' (yellow). The line types are defined as in Figure 11. The lower panels illustrate the same results as the upper ones, but not for models 2* (magenta), 3* (cyan), and 4* (yellow).</figcaption></figure>

<table><caption>**Table 5**
Parameters of the Orbital Libraries after Removing Box and Long-axis Loop Orbits</caption><thead><tr><th>Model</th><th>2*</th><th>3*</th><th>4*</th></tr></thead><tbody><tr><td>$N_{\text{orbit}*}$</td><td>19097</td><td>20381</td><td>21079</td></tr><tr><td>$\delta_{\text{smooth}*}$</td><td>0.053</td><td>0.046</td><td>0.059</td></tr><tr><td>$N_{>0,\text{smooth}*}$</td><td>1502</td><td>1570</td><td>1385</td></tr><tr><td>$-2K/W$</td><td>1.01</td><td>1.01</td><td>1.00</td></tr><tr><td>$\xi$</td><td>1.58</td><td>1.81</td><td>1.80</td></tr></tbody></table>

**Note.** The number of total orbits is denoted by $N_{\text{orbit}*}$, $\delta_{\text{smooth}}$ is the self-consistency parameter, and $N_{>0,\text{smooth}*}$ 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.

$T = 60 T_{\text{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 $\approx 1\%$ fraction of long-axis loop orbits in model 2' and, accordingly, the fraction of nonclassified orbits is slightly larger (by about $\approx 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

<figure><img src="image_16.png" /><figcaption>**Figure 14.** The virial ratio of model 1 perturbed by an external field of 1 $a_0$.</figcaption></figure>

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 $\approx 10\%$ fraction of nonclassified orbits in model 4*. While the shape of model 4' at $r_{50\%}$ is more prolate at $T = 60 T_{\text{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.

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.

## 6. The Effects of an External Field Shock

Models 1–5 have been constructed in (quasi-)equilibrium from the Schwarzschild models. This is the first attempt to adopt Schwarzschild’s technique for studying models embedded in external fields. It is therefore interesting to compare an equilibrium model embedded in an external field to an isolated model that suffers from a sudden perturbation due to an external field. The latter model is expected to relax to a new equilibrium state within the applied external field. Let us define a dynamical timescale for the model,

$$t_{\text{dyn}} \equiv \frac{r_{50\%}}{v_{\text{circ},r50\%}} = r_{50\%} (0.5 \times GM a_0)^{-1/4} \approx 3 T_{\text{simu}}, \quad (19)$$

where $r_{50\%}$ is the half-mass radius. A MOND system revirializing in a new environment has been studied by Wu & Kroupa (2013), where the system is initially embedded in a strong external field and then considered in isolation. The system fully revirializes within 5 $t_{\text{dyn}}$, defined as the Plummer radius for a Plummer sphere. Here we are interested in the inverse process.

We apply an external field with an amplitude of 1 $a_0$ to model 1, and the direction of the external field is assumed to be the same as for model 4. The new system is then evolved using the NMODY code. The time evolution of the virial ratio is presented in Figure 14. The system appears hotter than the
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Left panel: time evolution of the principal axes at $r_{90\%}$ for model 1 perturbed by a strong external field. The solid, dotted, and dashed lines denote the different components $\tilde{a}$, $\tilde{b}$, and $\tilde{c}$, respectively. Right panel: axial ratios at three evolutionary times as a function of enclosed mass. The upper set of curves refers to $p = \tilde{b}/\tilde{a}$ and the lower set to $q = \tilde{c}/\tilde{a}$. The solid, dotted, and dashed curves represent the initial ($T=0$), intermediate ($T=60 T_{\text{simu}}$), and final ($T=120 T_{\text{simu}}$) axial ratios.</figcaption></figure>

equilibrium state at $T=0$ since the potential becomes shallower in the strong external field, and the virial ratio is about 1.05 initially. The system relaxes to a new equilibrium state within about 60 $T_{\text{simu}} \approx 20 t_{\text{dyn}}$, after which the virial ratio oscillates around unity, with an amplitude less than 4%. Although the external field is as strong as $1 a_0$, the most impacted region is at larger radii where $r > 8$ kpc (see the deviation of the circular orbital time between models 1 and 4 in Figure 2). The circular orbital time for model 1 at 8 kpc, where there is over 80% of the total mass enclosed, is $\approx 0.3$ Gyr $\approx 60 T_{\text{simu}} \approx 20 t_{\text{dyn}}$. The crossing time for the farthest particles in model 1 is approximately 50 $T_{\text{simu}}$, which is comparable to the revirialization time. It explains why the revirialization time is longer than that observed in Wu & Kroupa (2013). Further, to be fully relaxed to and further test the stability of the new equilibrium state, the system has been evolved for another 60 $T_{\text{simu}}$.

The axial ratios for particles within $r_{90\%}$ for the shocked model are shown in the left panel of Figure 15. The values of $\tilde{p}_{90\%}$ and $\tilde{q}_{90\%}$ evolve within the first 40 $T_{\text{simu}}$, where $\tilde{p}_{90\%}$ decreases from 1.00 to 0.94, $\tilde{q}_{90\%}$ increases from 0.70 to 0.84, and the model becomes prolate, with axial ratios 1.00:0.94:0.84. The global shape then remains stable for the next 40 $T_{\text{simu}}$. At $T=80 T_{\text{simu}}$, there is an inflection for both $\tilde{p}_{90\%}$ and $\tilde{q}_{90\%}$, and the model becomes even more prolate, with an axial ratio as low as 1.00:0.88:0.80 at $T \approx 90 T_{\text{simu}}$. The ratios return to 1.00:0.92:0.87 at $T \approx 105 T_{\text{simu}}$, and the system evolves into a triaxial configuration. The behavior within 80–105 $T_{\text{simu}}$ implies an instability in the model’s evolution.

The right panel of Figure 15 shows the model’s axial ratios at $T=0$ (solid curves), after revirialization at $T=60 T_{\text{simu}}$ (dotted curves), and at $T=120 T_{\text{simu}}$ (dashed curves). The oblate axisymmetric model turns into a prolate one after revirialization in the external field at all Lagrangian radii. After reaching a new equilibrium state, the model further exhibits instability and becomes triaxial after another 60 $T_{\text{simu}}$ with axial ratios 1.00:0.92:0.87 at $r_{90\%}$. The resulting density configuration is less prolate than in model 4, suggesting that such shocked external field models are rounder than their counterparts that are initially constructed in (quasi-)equilibrium.

As the external field leads to a shallower potential, the model particles are less bound than in the unperturbed case. Hence,

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Anisotropy parameter $\beta(r)$ for the external field shocked model. The solid, dotted, and dashed curves represent the initial ($T=0$), intermediate ($T=60 T_{\text{simu}}$), and final ($T=120 T_{\text{simu}}$) states.</figcaption></figure>

the particles can move along more radial orbits, which is especially true for particles in the outer region. We thus expect an increase of $\beta(r)$ in the model’s outer region after revirialization and further evolution. This is confirmed in Figure 16. The anisotropy parameter $\beta(r)$ grows to 0.4 at $T=60 T_{\text{simu}}$ and to 0.5 at the end of the simulation. Within 1 kpc, $\beta(r)$ decreases and eventually reaches values close to zero. The final $\beta(r)$ profile within 1 kpc is very similar to that of model 1 because the external field does not significantly affect the orbits of the innermost particles.

For an adiabatic scenario in which the external field is gradually increased, we do not expect any prominent differences from the shocked model. Considering the opposite setup, that is, a system embedded in a strong external field that is suddenly isolated, Wu & Kroupa (2013) have not found any distinct features between the behaviors in the shocked and adiabatic situations, although the evolution of $\beta(r)$ is milder in the adiabatic case. Investigating this issue in more detail, however, is beyond the scope of this paper.

## 7. Lopsidedness

In Figure 1, we have illustrated the lopsided shape of the potential for an axisymmetric input density in an external field. The difference between the potential and the density distribution will lead to an additional relaxation of the model, and the resulting lopsided morphologies, which are further characterized by an offset between the density peak and the CoM, have previously been studied in Wu et al. (2010). Nevertheless, it is unknown whether the found lopsided shape is a stable feature. By defining the axis ratio

$$r^{-}: r^{+} = \sqrt{\sum_{i} m_{i} x_{i,-}^{2} / \sum_{i} m_{i}}: \sqrt{\sum_{j} m_{j} x_{j,+}^{2} / \sum_{j} m_{j}}, \quad (20)$$

where $\pm$ denotes the upper/lower half space with respect to the x axis, we can study the shape of the mass distribution at different times during the evolution. In both panels of Figure 17, we find that the axis ratios of the ICs are symmetric within $r \sim r_{90\%}$. After the evolution, the models become lopsided at large radii. The lopsided feature becomes stable
<figure><img src="image_19.png" /><figcaption>**Figure 17.** The evolution of lopsided shapes in terms of $r^{-}: r^{+}$, assuming model 4 (left panel) and model 5 (right panel). Different line types refer to different simulation snapshots.</figcaption></figure>

within 30 $T_{\text{simu}}$ for both models 4 and 5. The resulting values of $r^{-}: r^{+}$ for model 4 and model 5 are around 1.05 and 1.07, respectively. Both models 4 and 5 correspond to very bright ellipticals, and such lopsided features are expected to be more prominent for low-luminosity ellipticals within similar external fields.

This suggests that one cannot find perfectly symmetric galaxies in the centers of clusters within the framework of MOND. External fields will distort the original symmetry when the gravitational field strength starts to become comparable to that of the internal one, which is especially true for the outer parts of a galaxy.

The observations of 78 “nucleated” dwarf elliptical galaxies in Virgo (Binggeli et al. 2000) show that 20% of the sample galaxies are drastically lopsided. The associated typical centroid offset of these galaxies is about 100 pc (assuming that Virgo is at a distance of 20 Mpc). The other dwarfs in the sample also appear lopsided, but less dramatically so (see Figure 18). In the context of MOND, the magnitude of the offset associated with the lopsided shape should correlate with the strength of the external field. For galaxies far from the cluster center, the gravitational environment becomes weaker, and one expects smaller offsets. Figure 18 shows the absolute offsets of nuclei (left panel, in units of arcsec) and the relative offsets $\delta_r/r_{\text{eff}}$ (right panel). Most of these dwarf galaxies exhibit offsets that are larger than about 5% of their effective radii $r_{\text{eff}}$.

To assess whether observations of this kind could be used to constrain the MOND paradigm, we need to investigate how a three-dimensional correlation between offsets and distances from the cluster center would appear when looking at the corresponding projected quantities. For this purpose, we have considered a simple numerical experiment. Assuming a linear relation between $\delta_r$ and the external field such that $\delta_r(a_0) = 200$ pc (Wu et al. 2010), we created several random realizations of ~100 galaxies within a Virgo-like potential, adopting the Navarro–Frenk–White (NFW; Navarro et al. 1997) profile of McLaughlin (1999). Ignoring any further errors or physical effects, we then simply compare estimates of the linear correlation coefficient $Q$ between offsets and radii before and after projecting along the line of sight. As a result, we find that the correlation is significantly reduced, with $-Q$ dropping from values of 0.9–0.95 down to 0.4–0.45. Given additional uncertainties such as the presence of tidal fields,

selection biases, and other degenerate effects (Binggeli et al. 2000), we conclude that the current data are not enough to falsify or support MOND. Future high-resolution observations for nearby clusters, however, might improve this situation.

## 8. Conclusions and Discussions

Following Schwarzschild's approach, we constructed (quasi-) equilibrium models for galaxies with a central cusp embedded in uniform external fields within the framework of MOND. For these models, we performed instability tests and kinematic analyses by means of N-body simulations that operate on a spherical grid (Londrillo & Nipoti 2009).

When applying Schwarzschild's method to galaxies in external fields, the internal potential is distorted, and the models are not exactly self-consistent with respect to the original analytic density profile. This leads to nonequilibrium ICs that relax to a dynamic equilibrium within a few simulation times (see Figure 5). Since the overall residual is only at the level of a few percent, the deviation from the equilibrium state is rather minor. For comparison, the isolated models discussed in Wu et al. (2009) are found in a perfect equilibrium.

Interpreting the results of our N-body simulations, we conclude that galaxy models within external fields appear unstable over at most ~30 dynamical times. The shapes of these systems evolve due to instability. The density profiles along the major axes of the strong external field models clearly change after 30 $T_{\text{simu}}$, while those of other models remain close to the initial Hernquist profiles. In contrast to the stable, isolated triaxial systems considered in Wu et al. (2009), the isolated axisymmetric model is also unstable due to illegal box and nonclassified orbits with low angular momentum. It evolves toward a triaxial model with axial ratio close to 1.00:0.90:0.71 within $r_{90\%}$. If box orbits are removed from the orbital library at the cost of self-consistency, the orbital selection procedure suppresses nonclassified orbits, and the isolated axisymmetric model becomes stable.

For equilibrium systems embedded in weak and intermediate external fields, the final evolved shapes are similar to that of the isolated axisymmetric model. The presence of a strong external field yields a more prolate shape for axisymmetric models based on the same initial density profile. Similarly, the final shape of the triaxial model we considered within a strong external field turns out slightly more prolate. For these models, both box and nonclassified orbits contribute to the instability. Long-axis loop orbits, which appear in the asymmetric potential due to the external field, however, do not. Further evidence for the instability of all models (including the isolated case) is provided by the temporal evolution of $\sigma_r(r)$ and $\beta(r)$, especially in their inner regions.

We have also studied the case of an isolated axisymmetric model that is perturbed by a strong external field. Shocked by the external field, this model revirialized to a new equilibrium state after 60 $T_{\text{simu}} \approx 20$ $t_{\text{dyn}}$, and then evolved through instability during the following 45 $T_{\text{simu}}$. The final shape of the model is also prolate, but rounder than that of the corresponding self-consistent model in the same external field. The shocked case further shows an increase of radial anisotropy in the outer region.

While there is observational evidence for lopsided dwarf galaxies in Virgo, it is inconclusive whether this could be clearly linked to external field effects in MOND for future data sets. Since the MOND lopsidedness appears on the system's
<figure><img src="image_20.png" /><figcaption>**Figure 18.** The offsets of nuclei in dwarf elliptical galaxies inside the Virgo cluster: the colored open circles are the 78 dwarf ellipticals observed on the sky (Binggeli et al. 2000), and the red cross denotes the position of Virgo's center. The colors on the left panel illustrate the centroids' offsets $\delta_r$ in units of arcsec ($1'' \sim 100$ pc), and those on the right panel relate to the relative offsets $\delta_r/r_{eff}$, where $r_{eff}$ is the effective radius of a dwarf galaxy.</figcaption></figure>

outskirts, accounting just for a small fraction of the total mass, the effect is expected to be rather small.

Here we used the simple $\mu$ function defined by Equation (3) for all considered models. Compared to the standard form used in Milgrom (1983a), it leads to a more gradual transition from the deep MOND regime to Newtonian gravity. If the standard $\mu$ function was adopted, the models would be subject to weaker MOND effects in the intermediate field regions where gravity is comparable to $a_0$ (i.e., at the radius $\approx 8$ kpc for models 1–4 and $\approx 11$ kpc for model 5; see Figure 2). The use of a different interpolating function in the Schwarzschild modeling will result in a change of the orbital structure. Since a MOND spherical system with radial anisotropy is more stable than a pure Newtonian model with exactly the same density distribution (Nipoti et al. 2011), one may expect more severe instability in MOND models based on the standard $\mu$ function if the fraction of box and nonclassified orbits with low angular momentum is comparable to what we have found in our analysis.

It is interesting to consider the consequences of an evolving external field associated with the formation of galaxy clusters or galaxy groups. As galaxy clusters grow in density, the external field amplitude grows as well, which suggests an increase of radially anisotropic triaxial cluster galaxies as a function of decreasing redshift. Such a trend might be testable through detailed statistics of substructure shapes from lensing data of galaxy clusters. Considering the local universe, one surprising prediction of MOND is that the Andromeda galaxy should have had a flyby within 30 kpc of the Milky Way at about 7 Gyr ago (see Zhao et al. 2013; Banik & Zhao 2016, 2017). In this scenario, the two galaxies were drawn together by the greatly enhanced mutual attraction in MOND and flung out efficiently with very little dynamical friction, as shown in preliminary N-body and hydro simulations (I. Banik & F. Renaud 2017, private communications). From the view of the external field problem, one would expect that at the time of flyby, both galaxies suffered an external field shock

with a peak value of $g_{ext} \approx (200 \text{ km s}^{-1})^2/30 \text{ kpc} \approx 0.3a_0$. This external field might have had a stronger effect on the smaller system, our Galaxy, and could have triggered the formation of the triaxial bar in the center.

It is also worth cautioning that the external field effect is theory-dependent. For instance, the superfluid theory (Berezhiani & Khoury 2015, 2016; Khoury 2015, 2016) is one of the recent attempts to conciliate MOND with the cold dark matter phenomenology in galaxy clusters. In such theories, the MOND effect is achieved only on scales of 50 kpc inside a small superfluid core of the cluster, and our predictions would *not* apply to galaxies located in the bulk or the outer parts of clusters.

We thank the anonymous referee for careful and helpful suggestions and comments to the earlier version of the manuscript. The authors thank Luca Ciotti, Pasquale Londrillo, and Carlo Nipoti for generously sharing their code. The main body of the work has been accomplished at the University of Science and Technology of China and is supported by the NSFC grants 11503025, 11421303, and by “the Fundamental Research Funds for the Central Universities.” An early stage of this work has been performed under the Project HPC-EUROPA (211437), with the support of the European Community—Research Infrastructure Action under the FP8 “Structuring the European Research Area” program. X.W. thanks Simon Portegies Zwart for his warm hospitality at Leiden University. X.W. is grateful to James Binney and Keith Horne for their valuable comments and suggestions at an early stage of the project. X.W. thanks the support from “Hundred Talents Project of Anhui Province.” Part of this project was carried out during the Ph.D. program of X.W. X.W., M.F., and H.Z. acknowledge partial support from the Scottish Universities Physics Alliance (SUPA). Y.G.W. acknowledges support of the 973 Program (No. 2014CB845703) and the NSFC grants 11390372 and 1163004. The research of M.F. was supported by the I-CORE Program of the Planning and Budgeting Committee, the Israel Science Foundation (grants No. 1829/12
and No. 203/09), and the Asher Space Research Institute. M.F. acknowledges support through a fellowship from the Minerva Foundation.

# Appendix A
## Integration of Orbits

### A.1. Symmetry and Grid Segmentation of the Models

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.<fnref n="11" />

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 × 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 × 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 × 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.

### A.2. ICs for Orbits and Orbital Classification

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

$$ \frac{\max(3v_x^2) \cdot \max(3v_y^2) \cdot \max(3v_z^2)}{\max(v_x^2 + v_y^2 + v_z^2)^3} < 1, \quad (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 - \Phi_{int})}$, 0). Here $E$ corresponds to the total potential energy of the 20 equal-mass sectors along the x axis, and $\Phi_{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

$$ \max(L_x) \cdot \min(L_x) > 0. \quad (22) $$

Similarly, we define short-axis loop orbits as

$$ \max(L_z) \cdot \min(L_z) > 0. \quad (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_{stationary}$ and $N_{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.

<footnote n="11">¹¹ The work of Wang et al. (2008) ignored the two outermost sectors, so there were only 912 cells in the models. Here we discard only the last sector because the galactic outskirts might extend to larger radii in an external field. In addition, our available hardware has significantly improved, allowing us to consider these very diffuse outskirts.</footnote>
<figure><img src="image_21.png" /><figcaption>**Figure 19.** Left panel: time evolution of the kinetic energy tensor components $K_{xx}$, $K_{yy}$, and $K_{zz}$. The colors black, magenta, cyan, yellow, and green represent models 1–5, respectively. The solid, dotted, and dashed lines denote the different components. Right panel: time evolution of the root mean square velocities. The colors are defined as in the left panel.</figcaption></figure>

## Appendix B

### Further Dynamical Studies on the Systems' Stability

#### B.1. Kinetic Energy

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_{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

$$K_{xx} = \frac{1}{2} \frac{\sum_{i} m_i v_x^2}{\sum_{i} m_i} \qquad (24)$$

for $r < 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_{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_{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_{simu}$. Compared to their initial

values, $K_{yy}$ decreases by about 5% and 13% for models 2 and 3, respectively, after 60 $T_{simu}$. For model 4, where the external field is the strongest, $K_{yy}$ falls from about 50 $T_{simu}$ and becomes around 8% smaller than at $T = 0$. The component $K_{zz}$ of model 2 starts to decrease from 10 $T_{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.

The situation for model 5 is more complex. At the beginning of the simulation, we have $K_{xx} > K_{yy} > 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_{simu}$. The component $K_{yy}$ reaches its maximum value (up to $\approx 10\%$ relative to its original value) at $\approx 25 T_{simu}$, and $K_{xx}$ becomes minimal (reduced by $\approx 13\%$ relative to the original value, at $\approx 30 T_{simu}$). The two components meet again at $\approx 42 T_{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_{rms}$. For all models, the $v_{rms}$ profiles remain stable with tiny ($\pm 4\%$) and smooth oscillations. Therefore, the overall kinetic energy is conserved during the simulations.

#### B.2. Angular Momentum

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
<figure><img src="image_22.png" /><figcaption>**Figure 20.** Evolution of angular momenta $L_x$, $L_y$, and $L_z$ of models with external fields (models 1 and 3). The colors black, magenta, cyan, yellow, and green represent models 1–5, respectively. The solid, dotted, and dashed lines denote the different components.</figcaption></figure>

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_c$,

$$L_x = \frac{1}{M} \sum_{i=1}^{N_p} m_i L_x^{(i)}, \qquad (25)$$

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

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_{simu}$. For models 3–5, there are interesting evolutions of angular momentum. The component $L_y$ of model 3 starts increasing from 30 $T_{simu}$ and has a small value of 0.007 $L_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_c$ at 60 $T_{simu}$. The components $L_x$ and $L_z$ begin to increase from ≈20 $T_{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_c$ and $L_z \approx 0.022 L_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 ≈20 $T_{simu}$ and $L_z \approx -0.02$ at 60 $T_{simu}$. The growth of angular momentum at around 20–30 $T_{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 kpc km s⁻¹. Consequently, any self-rotation in the models caused by the external fields must be tiny.

## References

* Angus, G. W. 2008, MNRAS, 387, 1481

* Angus, G. W., Gentile, G., & Famaey, B. 2016, A&A, 585, A17

* Antonini, F., Capuzzo-Dolcetta, R., & Merritt, D. 2008, AN, 329, 900

* Antonini, F., Capuzzo-Dolcetta, R., & Merritt, D. 2009, MNRAS, 399, 671

* Banik, I., & Zhao, H. 2016, MNRAS, 459, 2237

* Banik, I., & Zhao, H. 2017, MNRAS, 467, 2180

* Barnes, E. I., Williams, L. L. R., Babul, A., & Dalcanton, J. J. 2005, ApJ, 634, 775

* Bekenstein, J., & Milgrom, M. 1984, ApJ, 286, 7

* Bekenstein, J. D. 2006, ConPh, 47, 387

* Bellovary, J. M., Dalcanton, J. J., Babul, A., et al. 2008, ApJ, 685, 739

* Berezhiani, L., & Khoury, J. 2015, PhRvD, 92, 103510

* Berezhiani, L., & Khoury, J. 2016, PhLB, 753, 639

* Bertin, G., Pegoraro, F., Rubini, F., & Vesperini, E. 1994, ApJ, 434, 94

* Bertin, G., & Stiavelli, M. 1989, ApJ, 338, 723

* Bílek, M., Bartošková, K., Ebrová, I., & Jungwiert, B. 2014, A&A, 566, A151

* Bílek, M., Jungwiert, B., Jílková, L., et al. 2013, A&A, 559, A110

* Binggeli, B., Barazza, F., & Jerjen, H. 2000, A&A, 359, 447

* Binney, J. 2005, MNRAS, 363, 937

* Binney, J., & Tremaine, S. 1987, Galactic Dynamics (Princeton, NJ: Princeton Univ. Press)
Brada, R., & Milgrom, M. 1995, MNRAS, 276, 453

Brada, R., & Milgrom, M. 1999, ApJ, 519, 590

Brada, R., & Milgrom, M. 2000, ApJ, 541, 556

Capuzzo-Dolcetta, R., Leccese, L., Merritt, D., & Vicari, A. 2007, ApJ, 666, 165

Ciotti, L., Londrillo, P., & Nipoti, C. 2006, ApJ, 640, 741

Crane, P., Stiavelli, M., King, I. R., et al. 1993, AJ, 106, 1371

De Rijcke, S., & Debattista, V. P. 2004, ApJL, 603, L25

de Zeeuw, P. T., Bureau, M., Emsellem, E., et al. 2002, MNRAS, 329, 513

de Zeeuw, P. T., & Lynden-Bell, D. 1985, MNRAS, 215, 713

Dehnen, W. 1993, MNRAS, 265, 250

Djorgovski, S., & Davis, M. 1987, ApJ, 313, 59

Dury, V., de Rijcke, S., Debattista, V. P., & Dejonghe, H. 2008, MNRAS, 387, 2

Famaey, B., & Binney, J. 2005, MNRAS, 363, 603

Famaey, B., Bruneton, J.-P., & Zhao, H. 2007a, MNRAS, 377, L79

Famaey, B., Gentile, G., Bruneton, J.-P., & Zhao, H. 2007b, PhRvD, 75, 063002

Famaey, B., & McGaugh, S. S. 2012, LRR, 15, 10

Fehlberg, E. 1968, Classical Fifth-, Sixth-, Seventh- and Eighth-Order Rung-Kutta Formullas with Stepsize Control, NASA Technical Report R-287c (Washington, DC: NASA)

Ferrarese, L., van den Bosch, F. C., Ford, H. C., Jaffe, W., & O'Connell, R. W. 1994, AJ, 108, 1598

Gadotti, D. A. 2009, MNRAS, 393, 1531

Haghi, H., Baumgardt, H., Kroupa, P., et al. 2009, MNRAS, 395, 1549

Haynes, M. P., van Zee, L., Hogg, D. E., Roberts, M. S., & Maddalena, R. J. 1998, AJ, 115, 62

Hernquist, L. 1990, ApJ, 356, 359

Khoury, J. 2015, PhRvD, 91, 024022

Khoury, J. 2016, PhRvD, 93, 103533

Lauer, T. R., Ajhar, E. A., Byun, Y.-L., et al. 1995, AJ, 110, 2622

Levison, H. F., Duncan, M. J., & Smith, B. F. 1990, ApJ, 363, 66

Londrillo, P., & Nipoti, C. 2009, MSAIS, 13, 89

Lüghausen, F., Famaey, B., Kroupa, P., et al. 2013, MNRAS, 432, 2846

May, A., & Binney, J. 1986, MNRAS, 221, 13P

McGaugh, S., & Milgrom, M. 2013, ApJ, 775, 139

McGaugh, S. S. 2004, ApJ, 609, 652

McGaugh, S. S. 2016, ApJL, 832, L8

McGaugh, S. S., de Blok, W. J. G., Schombert, J. M., Kuzio de Naray, R., & Kim, J. H. 2007, ApJ, 659, 149

McLaughlin, D. E. 1999, ApJL, 512, L9

Merritt, D., & Fridman, T. 1996, ApJ, 460, 136

Milgrom, M. 1983a, ApJ, 270, 384

Milgrom, M. 1983b, ApJ, 270, 365

Milgrom, M. 2009, MNRAS, 399, 474

Moller, P., Stiavelli, M., & Zeilinger, W. W. 1995, MNRAS, 276, 979

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Nipoti, C., Ciotti, L., & Londrillo, P. 2011, MNRAS, 414, 3298

Nipoti, C., Londrillo, P., & Ciotti, L. 2007, MNRAS, 381, L104

Perets, H. B., Wu, X., Zhao, H. S., et al. 2009, ApJ, 697, 2096

Pfenniger, D. 1984, A&A, 134, 373

Polyachenko, V. L., & Shukhman, I. G. 1981, SvA, 25, 533

Poon, M. Y., & Merritt, D. 2004, ApJ, 606, 774

Richstone, D. O. 1980, ApJ, 238, 103

Richstone, D. O. 1984, ApJ, 281, 100

Richstone, D. O., & Tremaine, S. 1984, ApJ, 286, 27

Rix, H.-W., de Zeeuw, P. T., Cretton, N., van der Marel, R. P., & Carollo, C. M. 1997, ApJ, 488, 702

Rix, H.-W., & Zaritsky, D. 1995, ApJ, 447, 82

Saha, P. 1991, MNRAS, 248, 494

Saha, P. 1992, MNRAS, 254, 132

Sanders, R. H. 1990, A&ARv, 2, 1

Sanders, R. H., & McGaugh, S. S. 2002, ARA&A, 40, 263

Scarpa, R. 2006, in AIP Conf. Ser. 822, First Crisis in Cosmology Conf., ed. E. J. Lerner & J. B. Almeida (San Francisco, CA: ASP), 253

Schwarzschild, M. 1979, ApJ, 232, 236

Schwarzschild, M. 1982, ApJ, 263, 599

Schwarzschild, M. 1993, ApJ, 409, 563

Sellwood, J. A., & Evans, N. W. 2001, ApJ, 546, 176

Statler, T. S. 1987, ApJ, 321, 113

Tiret, O., & Combes, F. 2007, A&A, 464, 517

Tiret, O., & Combes, F. 2008, A&A, 483, 719

Tiret, O., Combes, F., Angus, G. W., Famaey, B., & Zhao, H. S. 2007, A&A, 476, L1

Tremblay, B., & Merritt, D. 1996, AJ, 111, 2243

Trenti, M., & Bertin, G. 2006, ApJ, 637, 717

van der Marel, R. P., Cretton, N., de Zeeuw, P. T., & Rix, H.-W. 1998, ApJ, 493, 613

Wang, Y., Mao, S., Long, R. J., & Shen, J. 2013, MNRAS, 435, 3437

Wang, Y., Zhao, H., Mao, S., & Rich, R. M. 2012, MNRAS, 427, 1429

Wang, Y.-G., Wu, X., & Zhao, H. 2008, ApJ, 677, 1033

Wu, X., Famaey, B., Gentile, G., Perets, H., & Zhao, H. 2008, MNRAS, 386, 2199

Wu, X., & Kroupa, P. 2013, MNRAS, 435, 728

Wu, X., & Kroupa, P. 2015, MNRAS, 446, 330

Wu, X., Zhao, H., Famaey, B., et al. 2007, ApJL, 665, L101

Wu, X., Zhao, H., & Famaey, B. 2010, JCAP, 6, 10

Wu, X., Zhao, H., Wang, Y., Linares, C., & Knebe, A. 2009, MNRAS, 396, 109

Zhao, H., Famaey, B., Lüghausen, F., & Kroupa, P. 2013, A&A, 557, L3

Zhao, H., & Tian, L. 2006, A&A, 450, 1005

Zhao, H. S. 1996, MNRAS, 283, 149