Lopsidedness of self-consistent galaxies by the external field effect of clusters

TL;DR

This study uses orbit-superposition and N-body simulations to explore how external gravitational fields from galaxy clusters induce lopsided, triaxial shapes in galaxies within the framework of Milgrom's MOND, revealing instability-driven shape evolution.

Contribution

It demonstrates that external fields cause axisymmetric galaxies to become triaxial and unstable, providing new insights into galaxy morphology in MOND environments.

Findings

Galaxies lose symmetry and become triaxial under external fields.

Instability arises from box and low angular momentum orbits.

External fields predict inherent lopsidedness in galaxy shapes.

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. These models represent relatively massive ellipticals with a Hernquist radial profile at various distances from the cluster centre. 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 non-classified orbits with low angular momentum. We also consider a self-consistent isolated system which is then placed in a strong external field and allowed to evolve freely. This model, just as 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 MOND universe generically predicts some lopsidedness of galaxy shapes.

Tables

Table 1. Table 1 Galaxy models and Schwarzschild parameters: The total mass M 𝑀 M of each model is expressed in units of 10 10 ​ M ⊙ superscript 10 10 subscript 𝑀 direct-product 10^{10}M_{\odot} , a , b 𝑎 𝑏 a,~{}b and c 𝑐 c are scale lengths in units of kpc kpc \,{\rm kpc} , and g ext subscript 𝑔 ext g_{\rm ext} denotes the strength of the external field in units of a 0 ≈ 3700 ​ ( km ​ s − 1 ) 2 ​ kpc subscript 𝑎 0 3700 superscript km superscript s 1 2 kpc a_{0}\approx 3700(\,{\rm km\,s}^{-1})^{2}\,{\rm kpc} . As 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 subscript 𝑁 cell N_{\rm cell} is the number of cells to impose self-consistency, N stationary subscript 𝑁 stationary N_{\rm stationary} is the number of stationary starting orbits built in the orbital library, and N ejecting subscript 𝑁 ejecting N_{\rm ejecting} is the number of ejecting orbits starting from the x 𝑥 x - z 𝑧 z plane. The total number of orbits is given by N orbit subscript 𝑁 orbit N_{\rm orbit} .

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

Table 2. Table 2 Model self-consistency and equilibrium: Here δ smooth subscript 𝛿 smooth \delta_{\rm smooth} and δ 𝛿 \delta are the self-consistency parameters with and without smoothing, respectively, N > 0 subscript 𝑁 absent 0 N_{\rm>0} and N > 0 , smooth subscript 𝑁 absent 0 smooth N_{\rm>0,~{}smooth} are the corresponding orbit numbers for non-zero weights (obtained from Eqs.  10 and 12 ), v rms subscript 𝑣 rms v_{\rm rms} is the root-mean-square velocity in units of km ​ s − 1 km superscript s 1 \,{\rm km\,s}^{-1} , and − 2 ​ K / W 2 𝐾 𝑊 -2K/W is the virial ratio of the initial conditions. The radial-to-tangential anisotropic ratio ξ 𝜉 \xi characterizes the radial instability. L c subscript 𝐿 c L_{\rm c} denotes the unit-mass angular momentum with circular velocity v c subscript 𝑣 c v_{\rm c} at radius r c = 1 ​ kpc subscript 𝑟 c 1 kpc r_{\rm c}=1\,{\rm kpc} , where L c = r c ​ v c subscript 𝐿 c subscript 𝑟 c subscript 𝑣 c L_{\rm c}=r_{\rm c}v_{\rm c} is expressed in units of km ​ s − 1 ​ kpc km superscript s 1 kpc \,{\rm km\,s}^{-1}\,{\rm kpc} .

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

Table 3. Table 3 Parameters of the new orbital library after removing the box orbits. The number of total orbits is denoted by N orbit ′ superscript subscript 𝑁 orbit ′ N_{\rm orbit}^{\prime} , δ smooth subscript 𝛿 smooth \delta_{\rm smooth} is the self-consistency parameter, and N > 0 , smooth ′ superscript subscript 𝑁 absent 0 smooth ′ N_{\rm>0,~{}smooth}^{\prime} is the number of orbits with non-zero weights. The values of the virial ratio, − 2 ​ K / W 2 𝐾 𝑊 -2K/W , and of the radial stability parameter, ξ 𝜉 \xi , are also listed.

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

Table 4. Table 4 Parameters of the new orbital libraries after removing box and non-classified orbits. The number of total orbits is denoted by N orbit ′′ superscript subscript 𝑁 orbit ′′ N_{\rm orbit}^{\prime\prime} , δ smooth subscript 𝛿 smooth \delta_{\rm smooth} is the self-consistency parameter, and N > 0 , smooth ′′ superscript subscript 𝑁 absent 0 smooth ′′ N_{\rm>0,~{}smooth}^{\prime\prime} is the number of orbits with non-zero weights. The values of the virial ratio, − 2 ​ K / W 2 𝐾 𝑊 -2K/W , and of the radial stability parameter, ξ 𝜉 \xi , are also listed.

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

Table 5. Table 5 Parameters of the orbital libraries after removing box and long-axis loop orbits. The number of total orbits is denoted by N orbit ∗ N_{\rm orbit}* , δ smooth subscript 𝛿 smooth \delta_{\rm smooth} is the self-consistency parameter, and N > 0 , smooth ∗ N_{\rm>0,~{}smooth}* is the number of orbits with non-zero weights. The values of the virial ratio, − 2 ​ K / W 2 𝐾 𝑊 -2K/W , and of the radial stability parameter, ξ 𝜉 \xi , are also listed.

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

Equations

$$\vec{\nabla }\cdot \left[\mu \left(\frac{|\vec{g}|}{a_0}\right)\vec{g}\right]=4\pi G{\rho }_b,\vec{g}=\vec{\nabla }{\Phi }_{\mathrm{int}}+{\vec{g}}_{\mathrm{ext}}.$$

$$\begin{array}{rl}\mu (x)& \sim xx\ll 1,\\ \mu (x)& \sim 1x\gg 1.\end{array}$$

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

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

$$r=\sqrt{{\left(\frac{x}{a}\right)}^2+{\left(\frac{y}{b}\right)}^2+{\left(\frac{z}{c}\right)}^2},$$

$$\begin{array}{rl}{r}_i& =2\tan \left[(i+0.5)0.5\pi /(n_r+1)\right]\mathrm{kpc},\\ {\theta }_j& =\pi \times (j+0.5)/n_{\theta },\\ {\varphi }_k& =2\pi \times k/n_{\varphi },\end{array}$$

$$\mu \left(\frac{|\vec{g}|}{a_0}\right)\vec{g}={\vec{g}}_N+\vec{\nabla }\times \vec{h},$$

$$\sum _{j=1}^{N_{\mathrm{orbit}}}{w}_j{O}_{ij}=m_i.$$

$$O_{ij}=\frac{\Delta {\tau }_j\times {\nu }_{ij}}{\Delta {\tau }_j\times {\nu }_j}=\frac{{\nu }_{ij}}{{\nu }_j},$$

$${\chi }^2=\frac{1}{N_{\mathrm{cell}}}\sum _{i=1}^{N_{\mathrm{cell}}}{\left(m_i-\sum _{j=1}^{N_{\mathrm{orbit}}}{w}_j{O}_{ij}\right)}^2.$$

$$\delta =\sqrt{{\chi }^2}/\overline{m},$$

$${\chi }_{\mathrm{smooth}}^2=\frac{1}{N_{\mathrm{cell}}}\sum _{i=1}^{N_{\mathrm{cell}}}{\left(m_i-\sum _{j=1}^{N_{\mathrm{orbit}}}{w}_j{O}_{ij}\right)}^2+\lambda \sum _{j=1}^{N_{\mathrm{orbit}}}{w}_j^2.$$

$$W=\int \rho \vec{x}\cdot \vec{\nabla }\Phi {\mathrm{d}}^{3}x$$

$$r_i=8\tan \left[(i+0.5)0.5\pi /(n_r+1)\right]\mathrm{kpc},$$

$$\begin{array}{rl}{T}_{\mathrm{simu}}& ={\left(\frac{GM}{a^3}\right)}^{-1/2}\\ & =4.7\times 10^6\mathrm{yr}{\left(\frac{M}{10^{10}{M}_{\odot }}\right)}^{-1/2}{\left(\frac{a}{1\mathrm{kpc}}\right)}^{3/2}.\end{array}$$

$${\tilde{I}}_{\mathrm{xx}}(r)$$

$${\tilde{I}}_{\mathrm{xy}}(r)$$

$$\tilde{a}(r)$$

$$\tilde{b}(r)$$

$$\tilde{c}(r)$$

$$\beta (r)=1-\frac{{\sigma }_{\theta }^2+{\sigma }_{\varphi }^2}{2{\sigma }_r^2},$$

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

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

$$\frac{\max \left(3v_x^2\right)\cdot \max \left(3v_y^2\right)\cdot \max \left(3v_z^2\right)}{\max {\left(v_x^2+v_y^2+v_z^2\right)}^3}<1,$$

$$\max \left(L_x\right)\cdot \min \left(L_x\right)>0.$$

$$\max \left(L_z\right)\cdot \min \left(L_z\right)>0.$$

$$K_{\mathrm{xx}}=\frac{1}{2}\frac{{\sum }_i{m}_i{v}_x^2}{{\sum }_i{m}_i}$$

$$L_{\mathrm{x}}=\frac{1}{M}\sum _{i=1}^{N_p}{m}_i{L}_x^{(i)}$$

Full text

1. Introduction

Elliptical galaxies widely exist in the Universe, either isolated or embedded in clusters of galaxies, and they have compact centres. Studies of their dynamical behaviour 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 to find 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, i.e. 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:

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

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

The non-negative linear superposition of orbits which 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. E.g., Statler (1987) found self-consistency of 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}$ (Moller et al., 1995; Crane et al., 1993; Ferrarese et al., 1994; Lauer et al., 1995). Low-luminosity ellipticals have steeper centres, $\gamma\approx 2$, while the most luminous ellipticals have shallower ones, $0\leq\gamma\leq 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 (McGaugh et al., 2007; Famaey et al., 2007b). 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 shell 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 disc 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, i.e. 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 towards the galactic centre (Wu & Kroupa, 2013).

Schwarzschild’s orbit-superposition method has already been applied within Milgromian dynamics: for example, models of ellipsoidal field galaxies were found 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 which 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.111Generally, it is desirable to have a full treatment of the problem including tidal effects. This would allow to explore other, more complex scenarios such as evolution in fast-varying backgrounds, e.g. a galaxy crossing the centre of cluster, and will be subject to future work. Such an idealised 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 centre. In these situations, an external field is mainly dominated by its uniform part, and tidal effects play only a subordinate role.

Posing already 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 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 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 stabilised by dense centers. Further, De Rijcke & Debattista (2004) found off-centered nuclei in flattened non-rotating 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 counter-rotating 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, there exist no stability studies in external or tidal fields. A numerical study thus provides a first step into this direction.

The paper is organised as follows: In §2, we introduce our basic setup and discuss the effect of external fields on static galaxy potentials. In §3, we use Schwarzschild’s technique to construct quasi-equilibrium models of galaxies, and perform a kinematic analysis of these systems in §4. Finally, we compare our results to the evolution of isolated galaxies, and conclude in §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 $\vec{g}_{\rm ext}$ reads (Wu et al., 2007)

[TABLE]

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

[TABLE]

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

[TABLE]

Although observationally excluded in the Solar System, this $\mu$-function is still a viable choice on galactic scales and 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),

[TABLE]

where

[TABLE]

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$–$3\,{\rm kpc}$. For all five models, we set $a=1\,{\rm 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_{\rm 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\,{\rm kpc}$ to the galactic centre) are comparable to each other. For the intermediate case (model $3$), the two accelerations become comparable at approximately $70\,{\rm kpc}$, whereas for the weak external field case (model $2$), this requires moving to very large radii, beyond hundreds of kpc.

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

[TABLE]

where $i=0$ .. $n_{r}$, $j=0$ .. $n_{\theta}-1$, and $k=0$ .. $n_{\phi}-1$.

Wu et al. (2008) and Wu et al. (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, i.e. in external field dominated regions where $|\vec{g}_{\rm int}|<<|\vec{g}_{\rm ext}|<<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\,{\rm 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\,{\rm kpc}$, where internal and external field are comparable. The density contours, however, are still axisymmetric. Consequently, models constructed with external fields may not be necessarily self-consistent. Considering that the radius where the strongest external field (model $4$) starts dominating the internal dynamics is at about $10~{}\,{\rm kpc}$ (enlcosing $\approx 82\%$ of the total mass), the models should not be significantly affected for the most part and 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):

[TABLE]

where $g_{N}$ is the Newtonian gravitational field, and $\vec{\nabla}\times\vec{h}$ is a solenoidal vector field determined by the condition that $\vec{g}$ can be expressed as the gradient of a scalar potential. Restricting ourselves to the axis parallel to the external field’s direction, $\vec{g}$ and $\vec{g}_{N}$ must be either parallel or anti-parallel, assuming that the symmetry centre coincides with the coordinate origin. Hence the curl term $\vec{\nabla}\times\vec{h}$ vanishes. The strength of the total gravitational acceleration $g=|\vec{g}|$ along the negative semi-axis in the external field’s direction, i.e. where the external field cancels part of the internal one, is $g^{-}=|g_{\rm ext}-g_{\rm int}|$ while on the positive semi-axis it is $g^{+}=g_{\rm int}+g_{\rm 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 into an arbitrary direction, such a system has no symmetries anymore, but the curl term in Eq. 7 only accounts for corrections on the level of $10\%$ (Brada & Milgrom, 1995).

Returning to the axisymmetric model $4$, Fig. 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_{\rm ext}+g_{\rm int}|$, and steeper in the third octant where $g^{-}=g_{\rm int}-g_{\rm ext}$. As can be seen from the right panel of Fig. 1, the lopsidedness of the potential reaches its maximum at roughly $10\,{\rm 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\,{\rm kpc}$. At small radii ($r<3\,{\rm kpc}$), the internal gravitational field dominates, and thus $r^{-}:r^{+}$ is basically $1$. At much larger radii ($r>>10\,{\rm 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(|\vec{g}_{\rm 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_{\rm orbit}$ be the number of the orbits in the library ($j\in N_{\rm orbit}$) and $N_{\rm cell}$ the total number of grid cells segmenting space ($i\in N_{\rm cell}$). Further, let $O_{ij}$ denote the fraction of time spent by the $j{\rm th}$ orbit in the $i{\rm th}$ cell. The weight and mass of the $j{\rm th}$ orbit are defined by $w_{j}$ and $m_{i}$, respectively, and they are related by the following set of linear equations:

[TABLE]

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{\rm th}$ orbit in equal time intervals $\Delta\tau_{j}$, and counting the numbers of the output points $\nu_{ij}$ in the $i{\rm th}$ cell. After that the elements are determined according to

[TABLE]

where $\nu_{j}$ is the total output number of the $j{\rm th}$ orbit. The previous analyses in Wang et al. (2008) and Wu et al. (2009) used non-equal time interval outputs, given by the variation of the gravitational field strength, $\Delta\tau_{j}^{\prime}\sim 1/|\vec{\nabla}\cdot\vec{g}|^{1/2}$. In this case, the real time intervals $\Delta\tau_{j}^{\prime}$ of the $j{\rm th}$ orbit are not constant anymore. However, in these previous studies the unevenness of the time intervals between the outputed 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 centres. Our present analysis does not make this approximation. Each orbit is integrated for $100$ times its circular orbital time, $100~{}T_{\rm cir}$, hence the equal time interval $O_{ij}$ does not give rise to additional inaccuracy. We will discuss the time integration of the orbits in Sec. 3.1. Further details related to grid segmentation, initial conditions, 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_{\rm cir}$. In Fig. 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~{}\,{\rm kpc}$. While the slope of model $1$ is approaching unity at large radii, those of models $2$–$5$ approach a value of $3/2$ at infinity. For the isolated MOND model $1$, the circular velocity turns constant at large radii, $v_{\rm c}=(GMa_{0})^{1/4}$, hence $T_{\rm 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 Newton-like behavior $T_{\rm 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 Eq. 8. The right-hand side of Eq. 8 denotes the mass of the $i{\rm th}$ cell which we obtain from Monte-Carlo simulations using the analytic density distribution. The linear system is then solved by applying a non-negative least-squares (NNLS) method which minimises the following quantity:

[TABLE]

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

[TABLE]

where $\overline{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 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_{\rm cell}$. The best solutions are typically non-unique with a very noisy phase-space distribution characterised by $N_{\rm orbit}-N_{\rm cell}$ and $N_{\rm cell}$ zero-weight and non-zero-weight orbits, respectively (Merritt & Fridman, 1996; Zhao, 1996; Rix et al., 1997). As the mass distributions given by Eq. 4 are smooth, however, one would desire to select orbits in a less noisy way, i.e. with little oscillations in the weights of neighbouring orbits in phase space. Introducing a regularisation 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$ neighbouring orbits when solving the NNLS.

Here we apply a simpler method of regularisation: We minimise the scatter of orbital weights by introducing a smoothing parameter $\lambda$, where $\lambda=N_{\rm orbit}^{-2}$ is chosen as in Zhao (1996). The regularisation 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, Eq. 10 is modified as

[TABLE]

Due to the regularisation, 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_{\rm smooth}=\sqrt{\chi^{2}_{\rm smooth}}/\overline{m}$ (where $\overline{m}$ is the mean Monte-Carlo mass in cells), for the smoothed models. Comparing with second line, we find that the regularisation leads to an increase of $\chi^{2}$ on the order of $10^{-7}-10^{-5}$. Since the models $2$–$5$ feature symmetrised orbits, i.e., additional orbits starting from the other three octants (the second, fifth, and sixth octants) with the same initial conditions (see Appendix A), the increment of orbits is equivalent to smoothing the orbital structure. The regularisation in Eq. 12 only slightly changes the number of orbits with non-zero weights (see Table 2), and thus does not considerably change the accumulative fraction of individual orbit families.

Fig.  3 shows the integrated contributions of orbits (for energies $<E$) to the system’s mass with (right panels) and without (left panels) regularisation. The individual contributions of long-axis loop (green), short-axis loop (bright blue), box (dark blue lines), and non-classified (purple lines) orbits are plotted against the energy $E$. The fraction of box orbits in model $1$ is clearly increases when the regularisation 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 regularisation. Reducing the number of cells and orbits to that of model $1$, we further recomputed orbits for model $2$ using Eq. 12. The resulting orbital structure at low energy turns out very similar to the non-smoothed solution of model $1$, indicating that the use of reflecting (symmetry) planes changes the orbital structure for models $2$–$5$. The total amount of orbits with low angular momentum, i.e. box and non-classified orbits, are quite large for these models. These orbits might introduce instability which will be studied in the later sections.

For the 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 regularisation. 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 Fig. 1), this simply reflects that stronger external fields impact a larger fraction of orbits. The fraction of non-classified orbits follows a similar behaviour. 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 make 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 non-classified orbits, making up almost 70% of the total mass. In Newtonian gravity, Merritt & Fridman (1996) have demonstrated that a galaxy constructed completely by 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 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 non-classified orbits become important for systems with lower symmetry.

As solutions obtained from Schwarzschild’s method are not unique, the orbital structure could change when adopting different regularisation 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 into this direction by performing a kinematic analysis of the previously introduced models, starting with $N$-body initial conditions (ICs) given by Schwarzschild’s approach.

4.1. Initial conditions 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 $j{\rm th}$ orbit is $n_{j}=w_{j}N$, where $w_{j}$ is the weight of the $j{\rm th}$ orbit. Particles are placed on the $j{\rm th}$ 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 $j{\rm th}$ orbit (see Fig. 2 ).222One can also randomly sample particles on the $j_{th}$ orbit from a uniform distribution. Since most of the orbits have small positive weights in our simulations, the number of particles on $j_{th}$ are quite small. A random sampling might introduce numerical noise $\propto 1/\sqrt{n_{j}}$, and could, therefore, have problems to reflect the real shape of the orbit if the weight is small. To avoid such problems, we choose an isochronous sampling.

Our galaxy models exhibit special symmetries which significantly decrease the amount of computing work (see Table 1 and Appendix A). Also, as the NNLS selects hundreds of non-classified orbits which 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, i.e. $W+2K=0$, where

[TABLE]

is the Clausius integral and $K$ is the kinetic energy of the system (Binney & Tremaine 1987). The virial ratios $-2K/W$ and the root-mean-square velocities of the ICs are listed in Table 2. We find that all five models satisfy the scalar virial theorem very well, with $-2K/W=1\pm 0.02$. The $v_{\rm rms}$-values in models $2$–$4$ are slightly smaller than in model $1$ due to the presence of external fields. Since the potential is shallower in an external field, $v_{\rm rms}\approx\sigma$ for pressure-supported systems becomes smaller. For the strongest external field (model $4$), the $v_{\rm rms}$-value is smallest, $\approx 5~{}\,{\rm km\,s}^{-1}$ less than in model $1$. Even in this case, the decrease of $v_{\rm 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, $v_{\rm 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~{}\,{\rm kpc}$. The predicted value of $\sigma$ in this galaxy is only $2.1~{}\,{\rm 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 subject to 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 non-classified orbits are characterised by low angular momentum, and thus 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_{\rm ip=1}^{\rm N}\frac{1}{2}m_{\rm ip}\sigma_{\rm r,ip}^{2}$, where the index ${\rm 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_{\rm crit}$. Various studies have found different values of $\xi_{\rm crit}$ (e.g., May & Binney 1986, $\xi_{\rm crit}\approx 2.2$; Saha 1991, $\xi_{\rm crit}\approx 1.4$; Saha 1992, $\xi_{\rm crit}\approx 2.3$; Bertin & Stiavelli 1989, $\xi_{\rm crit}\approx 1.9$; Bertin et al. 1994, $\xi_{\rm 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 centre 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 non-classified ones, nor 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_{\rm 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_{\rm crit}$. Since the external field breaks the potential symmetries associated with these models, however, one cannot make a conclusive 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

[TABLE]

and the other two remaining grid segments are the same as in § 2. To reduce the computational workload, the chosen angular resolution is lower than that adopted for the Schwarzschild modeling. The upper panel of Fig. 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~{}\,{\rm 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 which has previously been studied in Wu et al. (2008). The time unit used in the NMODY code is given by (Wu et al., 2009)

[TABLE]

Here the quantity $2\pi T_{\rm 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 Mpc. However, their central regions where there exist strong and nearly uniform gravitational backgrounds 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 time scale where this approximation holds, this gives around $60~{}T_{\rm 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 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_{\rm simu}$ (twice the value of $60~{}T_{\rm simu}$) to examine the systems’ behaviour beyond the actual simulation time interval. In what follows, we will restrict the discussion to within $60~{}T_{\rm simu}$; only virial ratios are presented for the fully simulated range of $120~{}T_{\rm 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 ($\pm 2\%$), and we expect these quasi-equilibrium ICs to quickly relax to dynamically virialsed systems at the beginning of the simulations. Figure 5 shows the evolution of the virial ratio $2K/|W|$ within $120~{}T_{\rm simu}$. As can be seen from the figure, the models revirialise within a few simulations times and then the virial ratios oscillate around $1$ for all five models, with an oscillation amplitude of roughly $0.05$ within $120~{}T_{\rm simu}$. Since the overall residuals of virial ratios at $T=0$ are at the level of few percent, the deviation from the exact equilibrium state is an 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 Fig. 6 illustrate the spherical radii $r_{0}(M)$ enclosing different fractions of the total mass $M$, i.e., 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 less symmetries for models $2$–$5$, and thus there are less mirror particles in phase space. Such a situation could result in a self-rotation which cannot be canceled due to 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 §B.1, we will further comment 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

[TABLE]

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 diagonalised, yielding eigenvalues ${\tilde{I}_{\rm x^{\prime}x^{\prime}}}(r),~{}{\tilde{I}_{\rm y^{\prime}y^{\prime}}}(r)$ and ${\tilde{I}_{\rm z^{\prime}z^{\prime}}}(r)$, where the primed coordinate system refers to the corresponding eigenframe. The associated principle axes,

[TABLE]

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 Fig. 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, i.e. $r_{90\%}$.333The oscillations of the radial densities emerge from numerical noise. The radial densities along the major axes are computed from the $N$-body simulation grids which are closest to the major axes. Approximately, there are $N_{P}/(n_{\theta}\times n_{\phi})\approx 250$ particles along the major axes, and Eq. 14 is used for radial binning. The $m_{th}$ radial bin contains $n_{p,m}$ particles, and the associated numerical noise is $\approx 1/\sqrt{n_{p,m}}$. The initial density distributions (black dotted curves) along 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_{\rm simu}$ and reach stable states afterwards. For model $2$, the density along major axes within $1$–$5~{}\,{\rm 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]~{}\,{\rm 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~{}\,{\rm 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$–$5$ are generally due to 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 Fig. 7 shows the axial ratios of all five models defined for $r_{90\%}$, ${\tilde{p}}_{\rm 90\%}={\tilde{b}}_{\rm 90\%}/{\tilde{a}}_{\rm 90\%}$ and ${\tilde{q}}_{\rm 90\%}={\tilde{c}}_{\rm 90\%}/{\tilde{a}}_{\rm 90\%}$. The ratios ${\tilde{p}}_{\rm 90\%}$ and ${\tilde{q}}_{\rm 90\%}$ significantly evolve for all models during $60~{}T_{\rm simu}$, where the evolution of ${\tilde{p}}_{\rm 90\%}$ is generally stronger than that of ${\tilde{q}}_{\rm 90\%}$. For all models, the values of ${\tilde{p}}_{\rm 90\%}$ become considerably smaller. We observe that the axisymmetric models evolve into a triaxial configuration, and the initially triaxial model turns slightly more prolate.

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

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

The models $2$–$4$ eventually evolve into triaxial systems, and the final state is reached after $\sim 55~{}T_{\rm simu}$ for the models $2,3$ and $\sim 25~{}T_{\rm simu}$ for model $4$. Note that the model analysis is made after rotating the systems into the reference frame defined by their principle axes. Hence, any effects of pattern rotation (see the evolution of angular momentum in Fig. 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}}_{\rm 90\%}:~{}{\tilde{b}}_{\rm 90\%}:~{}{\tilde{c}}_{\rm 90\%}~{}=~{}1.00:~{}0.86:~{}0.71$, which is close to the isolated case discussed in Wu et al. (2009). Both ${\tilde{p}}_{\rm 90\%}$ and ${\tilde{q}}_{\rm 90\%}$ start to increase within $10~{}T_{\rm simu}$, but then decrease until $25~{}T_{\rm simu}$, where ${\tilde{a}}_{\rm 90\%}:~{}{\tilde{b}}_{\rm 90\%}:~{}{\tilde{c}}_{\rm 90\%}~{}=~{}1.00:~{}0.87:~{}0.69$. Again, there is no significant evolution beyond this point, and one finds ${\tilde{a}}_{\rm 90\%}:~{}{\tilde{b}}_{\rm 90\%}:~{}{\tilde{c}}_{\rm 90\%}~{}=~{}1:00:~{}0.85:~{}0.69$ at $60~{}T_{\rm simu}$. Since the external field for this model points into the $x$-direction, i.e. perpendicular to $y-z$ plane, the curves for the components along $y$- and $z$-axis display a very similar running behaviour. The instability appears within the first $25~{}T_{\rm simu}$ for the models $4$ and $5$, after which they become stable. We also note that the evolution for the models $1$–$3$ is not too different, since the external field effect 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 Fig. 7. For the initial models (solid curves), the axial ratio agrees well with the analytic density distribution and does not evolve considerably with increasing mass for the isolated model (model $1$). For the models in external fields (models $2$–$4$), there are deviations between the analytical and model’s principal axes in the inner regions. There are small oscillations or spikes for the resulting $p={\tilde{b}}/{\tilde{a}}$ and $q={\tilde{c}}/{\tilde{a}}$ in model $5$, with amplitudes $<\pm 0.1$. These spikes are purely numerical features and emerge from the iterative procedure, used to determine the principal axes, which is quite sensitive to the chosen initial guess in this case.

The initial axial ratios of the 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 centres 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 centre has to be placed at the computational grid’s centre. In our simulations, we transformed the coordinates at every time step, moving the centre of mass (CoM) to the grid centre and changing to the frame where its velocity is zero. The CoM does not need to coincide with the galactic centre, i.e. the centre 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 an 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 to the external field, i.e. 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 the 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_{\rm simu}$ (dotted curves), the values of $p$ for the models $1$–$3$ drop to $\approx 0.94$ in the inner region and to $\approx 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 Fig. 7. The 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 non-classified orbits (see the orbital fraction in the right panels of Fig. 3) and changes the shape of the model. Model $4$ turns out more prolate when compared to the 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 summarise, our results imply that

an axisymmetric model with and without the external field effect is not stable in MOND, and that 2. 2.

the principle axes, describing the shape of a galaxy, generally evolve more substantially in the presence of a strong external field, and that 3. 3.

radial instability caused by both box orbits and non-classified orbits with low angular momentum can change the axial ratios, i.e. the shape of a galaxy. The origin of the instability will be studied later in Sec. 5.

4.4. Kinetics

From §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

[TABLE]

where $r$ is defined in Eq. 4 and $\sigma_{\theta}$, $\sigma_{\phi}$ are the tangential and azimuthal velocity dispersions, respectively.

The left panels of Fig. 8 show the time evolution of the radial velocity dispersion $\sigma_{r}(r)$.444The radial velocity dispersions are plotted in 21 radial bins of equal mass (see Sec. A). Due to the increased particle number in these bins, the numerical noise is much smaller than in the right panels of Fig. 6. The same binning procedure is applied for other quantities that are discussed hereafter. 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-30T_{\rm 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\,{\rm kpc})\approx 170\,{\rm km\,s}^{-1}$ for the 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_{\rm simu}$. At larger radii, dispersions $\sigma_{r}(r>1\,{\rm kpc})$ are about $10-20~{}\,{\rm km\,s}^{-1}$ higher than in the initial state. The inner region of $\sigma_{r}(r)$ for model $5$ decreases from $300~{}\,{\rm km\,s}^{-1}$ to $250~{}\,{\rm km\,s}^{-1}$ within $30T_{\rm simu}$, and the profile also appears flat within $1~{}\,{\rm kpc}$. The evolution from the ICs to the final stable state indicates a radial instability which might originate from box orbits or from non-classified orbits with low angular momentum.

For the axisymmetric models $1$–$4$, the profiles of the anisotropy parameter (right panels of Fig. 8), $\beta(r)$, are quite different for the isolated and non-isolated models. In the case of an isolated galaxy, model $1$, $\beta(r)$ is initially almost constant, $\approx 0.2$, in accordance with the results of Wu et al. (2009). When there is an external field (models $2$–$4$), however, the initial $\beta(r)$-profiles are no longer constant. In general, $\beta(r)$ declines in the inner region and start to increase from $\simeq 0.5~{}\,{\rm kpc}$ to $\simeq 4~{}\,{\rm kpc}$. The model with weakest external field has the highest radial anisotropy at $4\,{\rm kpc}$, $\beta(r\geq 5\,{\rm kpc})\approx 0.8$. The $\beta(r)$-profile stays almost constant in the outer region where $r>5~{}\,{\rm kpc}$. Given that the external field is very weak in model $2$, it seems unexpected that there is such a large difference between the 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 the models $2$–$4$ (see Table 1). The increment 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, labelled $\tilde{1}$, by launching orbits from the same octants as in model $2$ and increasing the number of orbits to $27360$. Generating the corresponding $N$-body ICs as before, we then computed the anisotropy profile of model $\tilde{1}$ and the result is shown in Fig. 9. We find that $\beta(r)$ for model $\tilde{1}$ is very similar to that of model $2$, confirming that the different numbers of cells and orbits yield the observed discrepancy between the $\beta(r)$-profiles of models $1$ and $2$ at $T=0$.

When the external fields become stronger (models $3$ and $4$), the maximal values of $\beta(r)$ are about $0.6$ at $4~{}\,{\rm kpc}$ and the $\beta(r)$-profiles fall off again at larger radii. The $\beta(r)$-profiles for ICs with the strongest external fields, e.g. 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 regularisation method is the same. It seems safe to conclude that for models of the same mass derived from the regularisation in Eq. 12, stronger external fields will generally lead to less radially anisotropic models at larger radii.

The coloured 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~{}\,{\rm kpc}$, but decreases in the inner region within $15-30~{}T_{\rm simu}$ after which it remains stable. The models with the strongest external field show the most prominent evolution within $1~{}\,{\rm kpc}$. At $60~{}T_{\rm simu}$, $\beta(r)$ is nearly zero in the central region, i.e. the centre is basically isotropic. This is consistent with the left panels of Fig. 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 non-classified orbits might be the origin of radial instability for the triaxial model $5$.

Our results are consistent with the shape evolution illustrated in Fig. 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 of Sec. 4 showed that axisymmetric models are unstable and that the models exhibit large fractions of box and non-classified 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 the models $1$–$4$. Thus box orbits no longer contribute to the mass of these systems, and the corresponding models are labelled through a prime, i.e. $1^{\prime}$, $2^{\prime}$, $3^{\prime}$ and $4^{\prime}$. The numbers of orbits, $N_{\rm orbit}^{\prime}$, in the new libraries are listed in Table 3. There are at least four thousand orbits for the axisymmetric models. The best-fit solutions of Eq.  12 are then computed based on the new orbital library. The new $\delta_{\rm smooth}$-parameters are on the order of $10^{-2}$ and all larger than for models which 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 non-zero weights in the new models is smaller than in the original ones. Fig. 10 shows the orbital structure of the four new axisymmetric models. In these models, the fractions of short-axis loop orbits generally increases. The isolated system (model $1^{\prime}$) is almost entirely comprised of short-axis loop orbits. The fraction of non-classified and long-axis orbits increases with growing external field strength. Model $2^{\prime}$ is dominated by short-axis ($\approx 60\%$) and non-classified orbits ($\approx 40\%$). The fraction of non-classified orbits in models $3^{\prime}$ and $4^{\prime}$ (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^{\prime}$.

Generating $N$-body ICs according to the scheme presented in Sec. 4.1.1, we have conducted the same stability tests as before. The virial ratios of the new ICs are listed in Tab. 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 Tab. 3). We find that $\xi<<\xi_{\rm crit}$ for model $1^{\prime}$ 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^{\prime}$ is expected to be stable, the situation for the other models is unclear and needs to be analysed in more detail.

The shape evolution of the new models at $r_{90\%}$ for $T=0$–$60~{}T_{\rm simu}$ is presented in the left panel of Fig. 11. In contrast to the original isolated model with box orbits, the global axial ratios for model $1^{\prime}$ does not evolve with time. The axial ratios as a function of enclosed mass are shown in the right panel of Fig. 11. The black dotted curves confirm that the shape of model $1^{\prime}$ does has not changed after $60~{}T_{\rm simu}$. Therefore, model $1^{\prime}$, which is mostly comprised of short-axis loop orbits and has only a tiny fraction of non-classified orbits, is stable. The global axial ratios, ${\tilde{p}}_{\rm 90\%}$ and ${\tilde{q}}_{\rm 90\%}$, of the models $2^{\prime}$–$4^{\prime}$ 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 non-classified orbits which contain little angular momentum. The right panel of Fig. 11 shows that the shape of model $2^{\prime}$ at $r_{40\%}$ has not significantly change after $60~{}T_{\rm simu}$. The values of $p$ and $q$ slowly decline with increasing enclosed mass, but rise again at $r_{80\%}$. At $r_{90\%}$, the axial ratios are ${\tilde{a}}~{}:~{}{\tilde{b}}~{}:~{}{\tilde{c}}=1.00~{}:~{}0.96~{}:~{}0.66$. The instability of model $2^{\prime}$ is driven by non-classified orbits which redistribute to higher energies. The axial ratio of model $3^{\prime}$ at $r_{\rm 90\%}$ declines to $1.0~{}:~{}0.64~{}:~{}0.52$ at $60~{}T_{\rm simu}$. Compared to the results for model $3$, this indicates that the shape of model $3^{\prime}$ becomes more prolate. For model $4^{\prime}$, 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 non-classified orbits play an important role for the instability in this case.

5.2. Box orbits and non-classified orbits removed

Non-classified orbits are defined as low angular momentum orbits in which the velocity cannot be restricted to a rectangular velocity box (Eq. A1). These orbits behave stochastically in grid space. In order to characterise the origin of the instability in the models $2^{\prime}$–$4^{\prime}$, we perform a further test by removing all non-classified orbits from their orbital library. The such obtained models are mostly comprised of short-axis loop orbits and labelled as models $2^{\prime\prime},~{}3^{\prime\prime}$ and $4^{\prime\prime}$. Their orbital structure is illustrated in the upper panels of Fig. 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^{\prime}$–$4^{\prime}$. This implies that the $N$-body ICs sampled from these Schwarzschild models are not in equilibrium. Therefore, non-classified orbits are a necessary ingredient for the self-consistency of the models $2$–$4$. The stability of these models is again examined using $N$-body runs as described in Sec. 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_{\rm simu}$ are illustrated in the upper panels of Fig. 13. We find that the models are still axisymmetric after $60~{}T_{\rm 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^{\prime\prime}$ at $60~{}T_{\rm simu}$ turn out very similar to that of model $1^{\prime}$, indicating that by removing box and the non-classified orbits, axisymmetric models embedded in weak external fields evolve almost like isolated ones. Nevertheless, the models become eventually 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 non-classified orbits with low angular momentum likely play an important role for the instability in the original models $3$ and $4$ discussed in Sec. 4.

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^{\prime}$–$4^{\prime}$ to another set of models in which these illegal orbits are removed. The latter are labelled as models $2*,~{}3*$ and $4*$, and their orbital structure and evolution is analysed following the same procedure as in Sec. 5.1. The orbital structure of the models $2*$–$4*$ is shown in the lower panel of Fig. 12. These additional models are comprised of only short-axis loop and non-classified orbits. As can read from Table 5, all models yield self-consistency parameters $\delta\approx 0.05$, which is close to the corresponding values for models $2^{\prime}$–$4^{\prime}$. If long-axis loop orbits introduced another source of instability, the inferred axial ratios of models $2*$–$4*$ should differ from that of the models $2^{\prime}$–$4^{\prime}$ where only box orbits have been removed. The axial ratios of models $2*$–$4*$ are illustrated in the lower panel of Fig. 13. The time evolution of axial ratios at $r_{90\%}$ (lower left panel) is quite similar to that found in models $2^{\prime}$–$4^{\prime}$. Moreover, the axial ratios as a function of enclosed mass at $T=0$ and $T=60~{}T_{\rm simu}$ (lower right panel) are not significantly altered relative to that found for the models $2^{\prime}$–$4^{\prime}$. 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 the models $2^{\prime}$–$4^{\prime}$ 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^{\prime}$ drops more quickly from $r_{40\%}$ and suddenly increases at $r_{80\%}$. Although the orbital structures of the models $2^{\prime}$ and $2*$ are nearly identical, there is an $\approx 1\%$ fraction of long-axis loop orbits in model $2^{\prime}$ and, accordingly, the fraction of non-classified 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\%}$, i.e., a better agreement with the axial ratios of the original profiles. The slight changes in the orbital structure and non-classified orbit fractions are probably the main reasons for the observed differences between the two models. For the models $4*$ and $4^{\prime}$, the changes in the orbital structure are more pronounced, with a $10\%$ fraction of long-axis loop orbits in model $4^{\prime}$ and an increased $\approx~{}10\%$ fraction of non-classified orbits in model $4*$. While the shape of model $4^{\prime}$ at $r_{50\%}$ is more prolate at $T=60~{}T_{\rm simu}$, that of model $4*$ turns out more triaxial. The difference between the models $4^{\prime}$ and $4*$ is due to the different amounts of non-classified orbits.

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

6. The effects of an external field shock

The 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 which 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 time scale for the model,

[TABLE]

where $r_{50\%}$ is the half-mass radius. A MOND system revirialising in a new environment has been studied by Wu & Kroupa (2013), where the system is initially embedded in a strong external field and the considered in isolation. The system fully revirialises within $5~{}t_{\rm dyn}$ defined at 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 Fig. 14. The system appears hotter than the 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_{\rm simu}\approx 20~{}t_{\rm 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\,{\rm kpc}$ (see the deviation of the circular orbital time between models $1$ and $4$ in Fig. 2). The circular orbital time for model $1$ at $8\,{\rm kpc}$, where there is over $80\%$ of the total mass enclosed, is $\approx 0.3~{}Gyr\approx 60~{}T_{\rm simu}\approx 20~{}t_{\rm dyn}$. The crossing time for the furthest particles in model $1$ is approximately $50~{}T_{\rm simu}$, which is comparable to the revirialisation time. It explains why the revirialisation 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_{\rm simu}$.

The axial ratios for particles within $r_{90\%}$ for the shocked model are shown in the left panel of Fig. 15. The values of ${\tilde{p}}_{\rm 90\%}$ and ${\tilde{q}}_{\rm 90\%}$ evolve within the first $40~{}T_{\rm simu}$, where ${\tilde{p}}_{\rm 90\%}$ decreases from $1.00$ to $0.94$, ${\tilde{q}}_{\rm 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_{\rm simu}$. At $T=80~{}T_{\rm simu}$, there is an inflection for both ${\tilde{p}}_{\rm 90\%}$ and ${\tilde{q}}_{\rm 90\%}$, and the model becomes even more prolate, with an axial ratio as low as $1.00~{}:~{}0.88~{}:~{}0.80$ at $T\approx 90T_{\rm simu}$. The ratios return to $1.00~{}:~{}0.92~{}:~{}0.87$ at $T\approx 105T_{\rm simu}$, and the system evolves into a triaxial configuration. The behaviour within $80$–$105~{}T_{\rm simu}$ implies an instability in the model’s evolution.

The right panel of Fig. 15 shows the model’s axial ratios at $T=0$ (solid curves), after revirialisation at $T=60~{}T_{\rm simu}$ (dotted curves), and at $T=120~{}T_{\rm simu}$ (dashed curves). The oblate axisymmetric model turns prolate into a prolate one after revirialisation 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_{\rm 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 which 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, 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 revirialisation and further evolution. This is confirmed in Fig. 16. The anisotropy parameter $\beta(r)$ grows up to $0.4$ at $T=60~{}T_{\rm simu}$ and to $0.5$ at the end of the simulation. Within $1~{}\,{\rm kpc}$, $\beta(r)$ decreases and eventually reaches values close to zero. The final $\beta(r)$-profile within $1~{}\,{\rm 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, i.e. a system embedded in a strong external field that is suddenly isolated, Wu & Kroupa (2013) have not found any distinct features between the behaviour 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 Fig. 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 characterised by an offset between the density peak and the centre of mass, 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

[TABLE]

where $\pm$ denotes the upper/lower half space with respect to the $x$-axis, we can study the the shape of mass distribution at different times during the evolution. In both panels of Fig. 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 within $30~{}T_{\rm 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 centres 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 are drastically lopsided. The associated typical centroid offset of these galaxies is about $100~{}\,{\rm pc}$ (assuming that Virgo is at a distance of $20~{}\,{\rm Mpc}$). The other dwarfs in the sample appear also lopsided, but less dramatic (see Fig. 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 centre, the gravitational environment becomes weaker and one expects smaller offsets. Figure 18 shows the absolute offsets of nuclei (left panel, in units of arc-second) and the relative offsets ${\delta_{r}/r_{\rm eff}}$ (right panel). Most of these dwarf galaxies exhibit offsets that are larger than about 5% of their effective radii $r_{\rm 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 centre 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~{}\,{\rm pc}$ (Wu et al., 2010), we created several random realizations of $\sim 100$ galaxies within a Virgo-like potential, adopting the NFW 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 into uniform external fields within the framework of MOND. For these models, we performed instability tests and kinematic analyses by means of $N$-body simulations which 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 non-equilibrium initial conditions which relax to a dynamic equilibrium within a few simulation times (see Fig. 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 $\sim 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_{\rm 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 non-classified orbits with low angular momentum. It evolves towards 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 non-classified 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 non-classified 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 which is perturbed by a strong external field. Shocked by the external field, this model revirialised to a new equilibrium state after $60~{}T_{\rm simu}\approx 20~{}t_{\rm dyn}$, and then evolved through instability during the following $45~{}T_{\rm 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 datasets. Since the MOND lopsidedness appears on the system’s 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 Eq. 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~{}\,{\rm kpc}$ for models $1-4$ and $\approx~{}11~{}\,{\rm kpc}$ for model $5$; see Fig. 2). The use of different interpolating function in the Schwarzschild modelling 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 non-classified 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 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 fly-by within 30 kpc of the Milky Way at about 7 Gyrs ago (cf. 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 (Banik, Renaud, private communications, 2017). From the view of the external field problem, one would expect that at the time of fly-by, both galaxies suffered an external field shock with a peak value of $g_{\rm ext}\approx(200\,{\rm km\,s}^{-1})^{2}/30\,{\rm kpc}\approx 0.3a_{0}$. This external field might have had a stronger effect on the smaller system, i.e. our Galaxy, and could have triggered the formation of the triaxial bar in the centre.

It is also worth cautioning that the external field effect is theory-dependent. For instance, the superfluid theory (Khoury, 2015, 2016; Berezhiani & 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.

9. Acknowledgments

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, 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” Programme. XW thanks Simon Portegies Zwart for his warm hospitality at Leiden University. XW is grateful to James Binney and Keith Horne for their valuable comments and suggestions at an early stage of the project. XW thanks the support from “Hundred Talents Project of Anhui Province”. Part of this project was carried out during the PhD programme of XW. XW, MF, and HZ acknowledge partial support from the Scottish Universities Physics Alliance (SUPA). YGW acknowledges support of the 973 Program (No. 2014CB845703), and the NSFC grants 11390372, 1163004. The research of MF 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. MF acknowledges support through a fellowship from the Minerva Foundation.

Appendix A A A. Integration of orbits

A.1. A.1. Symmetry and grid segmentation of the models

We use an orbital integration code which 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.555The work of Wang et al. (2008) ignored the two outermost sectors, thus 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.

In a previous study, Wang et al. (2008) and Wu et al. (2009) considered isolated triaxial models. Due to the three-folded symmetry, they took only the first octant into account, and then symmetrised the $O_{ij}$ by reflecting all orbits at the octant’s boundaries, i.e. 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 into negative $x$-direction, and thus orbits starting from the positive $x$-semi-space are not the mirror orbits of those starting from the negative $x$-semi-space, i.e. orbits on different sides of the $y$-$z$ plane behave differently. However, the orbits will still have a two-folded symmetry with respect to the $x$-$z$ and $x$-$y$ planes. Hence we used the first ($|x|$,$|y|$,$|z|$) and second ($-|x|$,$|y|$,$|z|$) octants to calculate orbits for model $5$. The total number of cells in this case is set to $960\times 2=1920$, where the outermost $96$ cells are again excluded.

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

A.1.1 A.2. Initial conditions 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); Wu et al. (2009). Here we only give a brief description. There are two sets of ICs for the starting points:

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 centre of the system and change the sign of their angular momenta. It is known that stationary initial conditions can produce box orbits. Defining a velocity box through

[TABLE]

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

Ejecting orbits starting from the $x$-$z$ plane with initial velocities $(v_{\rm x},v_{\rm y},v_{\rm z})=(0,\sqrt{2(E-\Phi_{\rm int})},0)$. Here $E$ corresponds to the total potential energy of the 20 equal mass sectors along the $x$-axis, and $\Phi_{\rm 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 which 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

[TABLE]

Similarly, we define short-axis loop orbits as

[TABLE]

These orbits cannot move across the centre 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 non-classified 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 the models $2$–$4$, we use four mirror octants as mentioned above. In Table 1, we summarise the ICs for the five models, where $N_{\rm stationary}$ and $N_{\rm 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 (Eq. 8) increase, and building the orbit library becomes computationally more expensive as has been specified in Sec 3.1.

Appendix B B B. Further dynamical studies on the systems’ stability

B.1. B.1. Kinetic energy

Since the inner and outermost parts of non-isolated 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 Fig. 6 show that the $90\%$ Lagrangian radii for all the models stay constant within $60~{}T_{\rm simu}$, we shall use the particles within $r_{\rm 90\%}$.

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

[TABLE]

for $r<r_{90\%}$ and analogously for the other directions. We highlight that $K_{\rm xx}$, $K_{\rm yy}$ and $K_{\rm zz}$ are computed with respect to the original $(x,~{}y,~{}z)$ coordinate frame instead of the eigenframe considered in Sec. 4.3. Since all particles have equal mass $m_{i}=M/N$, Eq.  B1 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_{\rm simu}$. The components $K_{\rm xx}$, $K_{\rm yy}$, and $K_{\rm zz}$ of model $1$ evolve almost identically. Compared to the initial values, oscillations in the individual components are smaller than $10\%$ within $60~{}T_{\rm simu}$. At the beginning, the kinetic energy components $K_{\rm xx}$ (solid lines) and $K_{\rm yy}$ (dotted line) of model $2$ (magenta), model $3$ (cyan) and model $4$ (yellow) coincide since the ICs are axisymmetric. The component $K_{\rm 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_{\rm yy}$ of the models $2$ and $3$ evolve quite different from that of model $1$. The component $K_{\rm yy}$ of the models $2$ and $3$ slightly drops from about $20~{}T_{\rm simu}$. Compared to their initial values, $K_{\rm yy}$ decreases by about $5\%$ and $13\%$ for model $2$ and $3$, respectively, after $60~{}T_{\rm simu}$. For model $4$, where the external field is the strongest, $K_{\rm yy}$ falls from about $50~{}T_{\rm simu}$ and becomes around $8\%$ smaller than at $T=0$. The component $K_{\rm zz}$ of model $2$ starts to decrease from $10_{\rm simu}$ and is reduced by about $8\%$ at the end of the simulation; for the 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 Sec. A, we have seen that the fraction of box orbits in the models $2$–$4$ decreases while the fraction of non-classified 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_{\rm xx}>K_{\rm yy}>K_{\rm zz}$ since the model is triaxial. The component $K_{\rm xx}$ starts to drop whereas $K_{\rm yy}$ rises, and they intersect at $\approx 14~{}T_{\rm simu}$. The component $K_{\rm yy}$ reaches its maximum value, (up to $\approx 10\%$ relative to its original value) at $\approx 25~{}T_{\rm simu}$, and $K_{\rm xx}$ becomes minimal, (reduced by $\approx 13\%$ relative to the original value, at $\approx 30~{}T_{\rm simu}$. The two components meet again at $\approx 42~{}T_{\rm simu}$. The component $K_{\rm zz}$ does not considerably evolve, and oscillates around the initial value with an amplitude less than $8\%$.

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

B.2. 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 momenta for all models. The results are shown in Fig. 20 which illustrates the unit-mass angular momentum components $L_{\rm x}$, $L_{\rm y}$ and $L_{\rm z}$ in units of $L_{\rm c}$,

[TABLE]

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

A non-zero 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 Fig. 20 confirm this. The angular momentum of model $2$ is also conserved during the simulation. Although there is a tiny non-zero value of $L_{\rm y}$, the values do not change noticeably within $60~{}T_{\rm simu}$. For the models $3$–$5$, there are interesting evolutions of angular momentum. The component $L_{\rm y}$ of model $3$ starts increasing from $30~{}T_{\rm simu}$ and has a small value of $0.007~{}L_{\rm c}$, which is comparable to the constant $L_{\rm 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_{\rm y}$ changes from the beginning of the simulation and grows up to $0.04~{}L_{\rm c}$ at $60~{}T_{\rm simu}$. The components $L_{\rm x}$ and $L_{\rm z}$ begin to increase from $\approx 20~{}T_{\rm 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_{\rm x}\approx 0.025~{}L_{\rm c}$ and $L_{\rm z}\approx 0.022~{}L_{\rm c}$ at the end of the simulation. The components $L_{\rm x}$ and $L_{\rm 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_{\rm z}$ of model $5$ grows in the opposite direction from $\approx 20~{}T_{\rm simu}$ and $L_{\rm z}\approx-0.02$ at $60~{}T_{\rm simu}$. The growth of angular momentum at around $20-30~{}T_{\rm 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 of angular momentum components are relatively small, corresponding to a few $\,{\rm kpc}~{}\,{\rm km\,s}^{-1}$. Consequently, any self-rotation in the models caused by the external fields must be tiny.