Evaluating the ability of triaxial Schwarzschild modelling to estimate properties of galaxies from the Illustris simulation

TL;DR

This study assesses the effectiveness of Schwarzschild's orbit-superposition method in accurately estimating galaxy properties from high-resolution simulated data, finding it reliably recovers mass profiles and orbital distributions with some biases.

Contribution

It provides a comprehensive evaluation of Schwarzschild modelling on simulated galaxies, highlighting its strengths and limitations in recovering galaxy properties.

Findings

Total mass within effective radius is recovered within ±15%.

Stellar mass is underestimated by ~24%, dark matter overestimated by ~38%.

Velocity anisotropy and orbital distributions are well reconstructed.

Abstract

We evaluate the capabilities of Schwarzschild's orbit-superposition method by applying it to galaxies from the large scale, high resolution Illustris simulation. Nine early-type galaxies with a range of triaxiality are selected, and we create mock integral field unit data for five line-of-sight projections of each galaxy. Each of the 45 mock data sets is taken as an independent observed galaxy. Using van den Bosch's 2008 triaxial Schwarzschild implementation, we assess model estimates of various galaxy properties, covering mass profiles, intrinsic shapes, stellar orbit distributions and velocity anisotropies. Total mass within $\overline{R_{\rm e}}$ is recovered well with average deviations within $\pm15$ percent. Stellar mass is underestimated by $\sim24$ percent and dark matter overestimated by $\sim38$ percent (assuming an NFW dark matter profile and allowing for degeneracy between stellar mass and dark matter mass). Using a gNFW profile, these values improve to $\sim13$ percent for stellar mass and $\sim18$ percent for dark matter. Axis ratio estimates show a moderate bias of $\Delta (b/a)=0.07$ and $\Delta (c/a)=0.14$ ($a\ge b\ge c$). Distributions of the orbit circularities $\lambda_z$ and $\lambda_x$, representing rotation about the minor and major axes, are well reconstructed. Separating orbits into thermal categories, our models match the average fractions of these categories to within $10$ percent. Velocity anisotropy is well estimated with values matching in the inner regions but becoming slightly radially biased in the outer regions. Overall, the galaxy property estimates we obtained using Schwarzschild modelling are not implausible and are representative of the simulated galaxies we modelled.

Tables

Table 1. Table 1: The catalog of nine simulated galaxies. From top to bottom are: (1) the galaxy morphology; (2) galaxy name (Illustris object number); (3) galaxy short name used in our paper; (4) total stellar mass log ⁡ ( M ∗ / M ⊙ ) subscript 𝑀 subscript 𝑀 direct-product \log(M_{*}/M_{\odot}) ; (5) central black hole mass log ⁡ ( M BH / M ⊙ ) subscript 𝑀 BH subscript 𝑀 direct-product \log(M_{\rm BH}/M_{\odot}) ; (6) number of stellar particles N star subscript 𝑁 star N_{\rm star} ; (7) number of dark matter particles N dark subscript 𝑁 dark N_{\rm dark} ; (8) total gas fraction f gas subscript 𝑓 gas f_{\rm gas} ; (9) average effective radius R e ¯ ¯ subscript 𝑅 e \overline{R_{\rm e}} ; (10) maximum data coverage R m subscript 𝑅 m {R_{\rm m}} ; (11) axis ratios p e subscript 𝑝 e p_{\rm e} , q e subscript 𝑞 e q_{\rm e} and triaxial parameter T e subscript 𝑇 e T_{\rm e} at R e ¯ ¯ subscript 𝑅 e \overline{R_{\rm e}} ; (12) viewing angles θ 𝜃 \theta and φ 𝜑 \varphi of five random projections. There are totally 5 × 9 = 45 5 9 45 5\times 9=45 mock data sets created.

Morphology		Oblates			Triaxials			Prolates
Galaxy name (subhalo)		73666	347122	336920	312924	175435	16942	73664	191220	222715
Short name		O1	O2	O3	T1	T2	T3	P1	P2	P3
$\log \left(M_{\ast }/M_{\odot }\right)$		11.76	11.21	10.74	11.56	11.61	11.54	12.05	11.99	11.72
$\log \left(M_{\mathrm{BH}}/M_{\odot }\right)$		9.02	8.58	7.71	8.73	9.34	8.91	9.84	9.65	9.30
$N_{\mathrm{star}}$ (${10}^5$)		7.01	2.04	0.65	4.52	5.01	4.35	13.98	12.14	6.55
$N_{\mathrm{dark}}$ (${10}^5$)		13.58	6.18	4.82	8.04	35.39	3.99	48.76	56.26	23.12
$f_{\mathrm{gas}}$ ($\%$)		0.19	0.17	1.36	0.40	0.31	0.00	0.27	0.14	0.13
$\overline{R_e}$ (kpc)		10.46	7.06	5.12	7.42	7.06	9.18	8.96	16.09	9.30
$R_m$ (kpc)		32.30	22.36	15.10	22.82	18.46	28.83	24.56	50.56	31.23
Intrinsic shapes	$p_e$	0.98	0.89	0.92	0.80	0.76	0.68	0.71	0.69	0.52
	$q_e$	0.44	0.34	0.59	0.42	0.55	0.45	0.63	0.61	0.50
	$T_e$	0.04	0.23	0.24	0.43	0.59	0.67	0.81	0.83	0.97
Viewing angles	$(\theta ,\phi )$	${(17,87)}^{\circ }$	${(27,78)}^{\circ }$	${(10,58)}^{\circ }$	${(40,119)}^{\circ }$	${(28,107)}^{\circ }$	${(3,77)}^{\circ }$	${(36,71)}^{\circ }$	${(18,19)}^{\circ }$	${(18,103)}^{\circ }$
		${(42,60)}^{\circ }$	${(53,79)}^{\circ }$	${(38,131)}^{\circ }$	${(63,123)}^{\circ }$	${(43,134)}^{\circ }$	${(45,114)}^{\circ }$	${(41,58)}^{\circ }$	${(45,116)}^{\circ }$	${(37,60)}^{\circ }$
		${(55,136)}^{\circ }$	${(61,24)}^{\circ }$	${(55,52)}^{\circ }$	${(70,47)}^{\circ }$	${(57,88)}^{\circ }$	${(49,65)}^{\circ }$	${(52,164)}^{\circ }$	${(52,54)}^{\circ }$	${(66,140)}^{\circ }$
		${(74,54)}^{\circ }$	${(74,151)}^{\circ }$	${(73,29)}^{\circ }$	${(76,3)}^{\circ }$	${(67,110)}^{\circ }$	${(62,165)}^{\circ }$	${(63,42)}^{\circ }$	${(67,107)}^{\circ }$	${(75,146)}^{\circ }$
		${(83,131)}^{\circ }$	${(85,92)}^{\circ }$	${(81,99)}^{\circ }$	${(86,38)}^{\circ }$	${(69,6)}^{\circ }$	${(86,172)}^{\circ }$	${(84,28)}^{\circ }$	${(87,174)}^{\circ }$	${(88,77)}^{\circ }$

Table 2. Table 2: The information of active orbits. From left to right are: (1) the mock data, contain galaxy names and viewing angles; (2) initial conditions n E × n I 2 × n I 3 subscript 𝑛 𝐸 subscript 𝑛 subscript 𝐼 2 subscript 𝑛 subscript 𝐼 3 n_{E}\times n_{I_{2}}\times n_{I_{3}} ; (3) the number of input orbits before dithering; (4) the number of active orbits; (5) the fraction of active orbits.

Mock data	Initial conditions $n_E\times n_{I_2}\times n_{I_3}$	Input orbits	Active orbits	Active fraction
galaxy O1 with ${(42,60)}^{\circ }$	$15\times 7\times 5$	1575	332	$21.1\%$
	$21\times 10\times 7$	4410	479	$10.9\%$
	$31\times 15\times 10$	13950	582	$4.2\%$
galaxy T2 with ${(69,6)}^{\circ }$	$15\times 7\times 5$	1575	343	$21.8\%$
	$21\times 10\times 7$	4410	494	$11.2\%$
	$31\times 15\times 10$	13950	567	$4.1\%$
galaxy P1 with ${(84,28)}^{\circ }$	$15\times 7\times 5$	1575	339	$21.5\%$
	$21\times 10\times 7$	4410	541	$12.3\%$
	$31\times 15\times 10$	13950	696	$5.0\%$

Equations

$$\Sigma (R^{\prime },{\theta }^{\prime })=\sum _{\mathrm{j}=1}^N\frac{L_j}{2\pi {\sigma }_j^{\prime 2}{q}_j^{\prime }}\exp \left[-\frac{1}{2{\sigma }_j^{\prime 2}}(x^{\prime 2}+\frac{y^{\prime 2}}{q_j^{\prime 2}})\right],$$

$$\left\{\begin{array}[]{ll}x^{\prime}=R^{\prime}\sin(\theta^{\prime}-\psi),&\\ y^{\prime}=R^{\prime}\cos(\theta^{\prime}-\psi),\end{array}\right.$$

$$\nu (x,y,z)=\sum _{\mathrm{j}=1}^N\frac{L_j}{(\sqrt{2\pi }{\sigma }_j{)}^3{p}_j{q}_j}\exp \left[-\frac{1}{2{\sigma }_j^2}(x^2+\frac{y^2}{p_j^2}+\frac{z^2}{q_j^2})\right],$$

$${\rho }_{\mathrm{star}}(x,y,z)=(M_{\ast }/L)\times \nu (x,y,z).$$

$${\rho }_{\mathrm{DM}}(r)=\frac{{\rho }_0}{\frac{r}{R_{\mathrm{s}}}{\left(1+\frac{r}{R_{\mathrm{s}}}\right)}^2},$$

$${\Phi }_{\mathrm{DM}}(r)=-\frac{4\pi G{\rho }_0{R}_s^3}{r}\ln (1+\frac{r}{R_s}).$$

$$\frac{M_{200}}{R_{200}^3\times 4/3\pi }=200{\rho }_{\mathrm{crit}},$$

$${\rho }_{\mathrm{crit}}=\frac{3H_0^2}{8\pi G},$$

$$M_{200}={\int }_0^{R_{200}}4\pi r^2\rho (r)dr=4\pi {\rho }_0{R}_{\mathrm{s}}^3\left[\ln (\frac{R_s+R_{200}}{R_s})-\frac{R_{200}}{R_s+R_{200}}\right].$$

$$\left\{\begin{array}[]{ll}c=R_{200}/R_{\rm s},&\\ f=M_{200}/M_{*},\end{array}\right.$$

$${\rho }_0=\frac{200}{3}\frac{c^3}{\ln (1+c)-c/(1+c)}\times {\rho }_{\mathrm{crit}},$$

$$R_{\mathrm{s}}={\left[\frac{3}{800\pi }\frac{M_{\ast }f}{{\rho }_{\mathrm{crit}}{c}^3}\right]}^{1/3}.$$

$${\Phi }_{\mathrm{BH}}(x,y,z)=-\frac{G{M}_{\mathrm{BH}}}{\sqrt{r_{\mathrm{soft}}^2+x^2+y^2+z^2}},$$

$$\left[\begin{array}{c}{x}^{\prime \prime }\\ y^{\prime \prime }\\ z^{\prime \prime }\end{array}\right]=\left[\begin{array}{ccc}\sin \psi & -\cos \psi & 0\\ \cos \psi & \sin \psi & 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ z^{\prime }\end{array}\right],$$

$$\left[\begin{array}{c}{x}^{\prime }\\ y^{\prime }\\ z^{\prime }\end{array}\right]=\left[\begin{array}{ccc}-\sin \phi & \cos \phi & 0\\ -\cos \theta \cos \phi & -\cos \theta \sin \phi & \sin \theta \\ \sin \theta \cos \phi & \sin \theta \sin \phi & \cos \theta \end{array}\right]\left[\begin{array}{c}x\\ y\\ z\end{array}\right].$$

$$1-q_j^2=\frac{{\delta }_j^{\prime }[2\cos 2\psi +\sin 2\psi (\sec \theta \cot \phi -\cos \theta \tan \phi )]}{2{\sin }^2\theta [{\delta }_j^{\prime }\cos \psi (\cos \psi +\cot \phi \sec \theta \sin \psi )-1]},$$

$$p_j^2-q_j^2=\frac{{\delta }_j^{\prime }[2\cos 2\psi +\sin 2\psi (\cos \theta \cot \phi -\sec \theta \tan \phi )]}{2{\sin }^2\theta [{\delta }_j^{\prime }\cos \psi (\cos \psi +\cot \phi \sec \theta \sin \psi )-1]},$$

$$u_j^2=\frac{1}{q_j^{\prime 2}}\sqrt{p_j^2{\cos }^2\theta +q_j^2{\sin }^2\theta (p_j^2{\cos }^2\phi +{\sin }^2\phi )},$$

$${\cos }^2\theta =\frac{(u^2-q^2)(q^{\prime 2}{u}^2-q^2)}{(1-q^2)(p^2-q^2)},$$

$${\tan }^2\phi =\frac{(u^2-p^2)(p^2-q^{\prime 2}{u}^2)(1-q^2)}{(1-u^2)(1-q^{\prime 2}{u}^2)(p^2-q^2)},$$

$${\tan }^2\psi =\frac{(1-q^{\prime 2}{u}^2)(p^2-q^{\prime 2}{u}^2)(u^2-q^2)}{(1-u^2)(u^2-p^2)(q^{\prime 2}{u}^2-q^2)},$$

$$T=\frac{1-p^2}{1-q^2}.$$

$${\lambda }_z=\overline{L_z}/(\overline{r}\times \overline{V_{\mathrm{rms}}}),$$

$${\lambda }_x=\overline{L_x}/(\overline{r}\times \overline{V_{\mathrm{rms}}}),$$

$${\beta }_r=1-\frac{\overline{{\sigma }_{\theta }^2}+\overline{{\sigma }_{\phi }^2}}{2\overline{{\sigma }_r^2}}\equiv 1-\frac{\overline{{\sigma }_t^2}}{2\overline{{\sigma }_r^2}},$$

$$f_t=\frac{(\overline{{\sigma }_{\theta }^2}+\overline{{\sigma }_{\phi }^2})/2}{\overline{{\sigma }_r^2}+(\overline{{\sigma }_{\theta }^2}+\overline{{\sigma }_{\phi }^2})/2}\equiv \frac{\overline{{\sigma }_t^2}}{2\overline{{\sigma }_r^2}+\overline{{\sigma }_t^2}}.$$

$${\chi }_{\mathrm{NNLS}}^2={\chi }_{\mathrm{lum}}^2+{\chi }_{\mathrm{kin}}^2.$$

$${\chi }_{\mathrm{lum}}^2=\sum _{\mathrm{n}=1}^N{\left(\frac{S_n^{\ast }-S_n}{0.01S_n}\right)}^2+\sum _{\mathrm{m}=1}^M{\left(\frac{{\rho }_n^{\ast }-{\rho }_n}{0.02{\rho }_n}\right)}^2,$$

$${\chi }_{\mathrm{kin}}^2=\sum _{\mathrm{n}=1}^N\left[{\left(\frac{h_{1n}^{\ast }-h_{1n}}{\Delta h_{1n}}\right)}^2+{\left(\frac{h_{2n}^{\ast }-h_{2n}}{\Delta h_{2n}}\right)}^2\right],$$

$${\chi }^2=\sum _{\mathrm{n}=1}^N\left[{\left(\frac{V_n^{\ast }-V_n}{\Delta V_n}\right)}^2+{\left(\frac{{\sigma }_n^{\ast }-{\sigma }_n}{\Delta {\sigma }_n}\right)}^2\right],$$

Full text

††pubyear: 2018††pagerange: Evaluating the ability of triaxial Schwarzschild modelling to estimate properties of galaxies from the Illustris simulation–LABEL:lastpage

1 Introduction

The challenge we face in galactic astronomy is how to analyse observations of a galaxy and not only understand what the observations themselves are telling us but what can also be understood or implied about the galaxy’s properties that are not directly observable. Such properties include for example the dark matter distribution, orbit structures, variations in velocity dispersion anisotropy and so on.

In general (our own Galaxy being a major exception), we work with surface photometry and integrated line-of-sight spectroscopic data as provided by an Integral Field Unit (IFU). In recent years, good quality data from IFU surveys such as SAURON (Bacon et al., 2001), $\rm ATLAS^{3D}$ (Cappellari et al., 2011), CALIFA (Sánchez et al., 2012), SAMI (Bryant et al., 2015), and MaNGA (Bundy et al., 2015) has become available. This data is characterised by a two dimensional on-sky spatial position with line of sight kinematics or line strength data. Not only do we wish to understand galaxy properties which are not directly observable but we also have to contend with degeneracies generated by projecting or deprojecting our observed data.

Various modelling schemes have been devised to assist investigations such as ours. The Schwarzschild orbit-superposition technique (Schwarzschild, 1979) we use in this paper is one of them. Other schemes include solutions of the Jeans equations like the Jeans-Anisotropic-MGE (JAM) models (Cappellari, 2008), the particle based Made-to-Measure (M2M) method (Syer & Tremaine, 1996), and action-angle distribution functions (e.g.,Binney, 2010) which have mainly been applied to observations of the Milky Way. All these schemes assume that if a model is able to reproduce the observed data from a galaxy then the model can also be used to determine underlying properties of the galaxy.

JAM is computationally quick to run but, having strong assumptions on galaxy properties, may not produce physically realisable solutions, and may not be able to determine the properties we are interested in. Schwarzschild’s method and M2M are more computationally time-consuming but have the advantage that the orbit properties of a galaxy may be assessed and this has contributed to our decision to use Schwarzschild’s method. However the solutions obtained from both these methods are heavily dependent on how well the initial conditions represent the galaxy being modelled.

Ideally before using the modelling schemes with real galaxies, the systematics of the schemes for their intended use should be well understood. This concerns not only whether or not observed data are reproduced but also the accuracy to which model estimates of underlying properties can be made. This will be a key feature of our investigation. Note that in order to assess accuracy we will need test galaxies with known properties. We will use mock (test) galaxies taken from a cosmological simulation for this purpose.

It is not unusual for galaxies to be modelled under various simplifying assumptions. For example, dark matter may be ignored, spherical symmetry may be assumed, and galaxy features ignored (e.g. bars and spiral arms). We will address this in part by including dark matter in our models, and by seeking to model early-type triaxial galaxies. Bars and spiral arms in late-type galaxies add modelling complexity through their lack of dynamical equilibrium and are excluded. Axisymmetric models rather than triaxial ones are often constructed to analyse galaxies. The main reason for this is that axisymmetric models have fewer parameters, and so they are simpler and faster to execute. However, in real observations such as MaNGA, many early-type galaxies have twisted velocity fields which can only be created by triaxial galaxies. Recent observations appear to suggest this (Li et al., 2018). In addition, cosmological simulations such as EAGLE (Schaye et al., 2015; Crain et al., 2015) and Illustris (Vogelsberger et al., 2014a, b; Genel et al., 2014; Nelson et al., 2015) show that many galaxies are in fact triaxial with some being close to oblate or alternatively close to prolate (Velliscig et al., 2015; Li et al., 2016).

Previous attempts to assess the systematic behaviour of Schwarzschild’s method in estimating specific properties of galaxies include van de Ven et al. (2008) (internal orbit sructure) and van den Bosch et al. (2009) (intrinsic shape) who used mock galaxies based on theoretical Abel models (Dejonghe & Laurent, 1991; Mathieu & Dejonghe, 1999). Thomas et al. (2007) considered recovery of mass distributions while Kowalczyk et al. (2017) attempted to recover both mass distribution and orbit anisotropy in mock dwarf galaxies. Zhu et al. (2018b) investigated the recovery of orbit circularity. For other modelling schemes, Li et al. (2016) investigated the behaviour of JAM, and Long & Mao (2012) extended the earlier Schwarzschild and Jeans equation work of Cappellari et al. (2006) by providing a M2M mass-to-light comparison.

We now consolidate the statements so far into a set of objectives for our investigation. Using an existing Schwarzschild implementation capable of modelling triaxial galaxies (van den Bosch et al., 2008) together with mock triaxial galaxy data constructed from a cosmological simulation, our objectives are to understand how well Schwarzschild’s method is able to estimate the underlying properties of our mock galaxies. We consider four properties in particular,

the mass profile including both stellar and dark matter, 2. 2.

galaxy morphology, 3. 3.

orbit circularity, and 4. 4.

velocity dispersion anisotropy.

To our knowledge, this is the first time a quantitative assessment of the estimation of the four properties has been undertaken for early-type triaxial galaxies using Schwarzschild’s method.

The structure of the paper is as follows. In $\S$ 2 we explain the approach we will take to address the objectives. In $\S$ 3 we introduce Schwarzschild’s method while in $\S$ 4 we present lower level details covering the stellar and dark matter potentials, intrinsic shape parameters, and parameters used to characterize orbit distributions. Information about the mock galaxies we select is given in $\S$ 5. In $\S$ 6, we present the results of our parameter estimations based on Schwarzschild modelling and compare our model values with the true galaxy values. In $\S$ 7, we discuss any biases revealed in attempting to recover the properties of our mock galaxies. Finally in $\S$ 8, we summarize our research, indicating where we have met our objectives and where outstanding issues remain to be investigated.

2 Approach

The previous section was concerned with setting out the scientific context and objectives for our investigation. This section describes at a top level the approach we are going to take to address those objectives. To recap, we are using an existing Schwarzschild implementation capable of modelling triaxial galaxies together with mock triaxial galaxy data constructed from a cosmological simulation to try and understand how well Schwarzschild’s method is able to estimate underlying properties of our mock galaxies.

We employ the software implementation of Schwarzschild’s method described in van den Bosch et al. (2008). As indicated in the paper, this implementation draws on design ideas from Rix et al. (1997), van der Marel et al. (1998), Cretton et al. (1999), Verolme et al. (2002) and Cappellari et al. (2006). In particular, it is able to model triaxial galaxies using Gauss-Hermite coefficients from line of sight velocity distributions as constraints, and uses orbit dithering to reduce noise in model orbit-based calculations. The software has not been parallelised. Other than a number of small minor changes targetted at improving primarily parameter handling, no changes have been made to the van den Bosch software. The core software implementing the Schwarzschild mechanisms is unaltered.

The gravitational potential for the stellar matter in a triaxial galaxy is straightforward to calculate using the Multi-Gaussian Expansion (MGE) formalism (Emsellem et al., 1994a; Cappellari, 2002). MGEs are able to handle the multiple viewing angles required by triaxial modelling and their use is an integral feature of the software implementation we have chosen. The shapes of dark matter haloes are hard to determine through observations, so we we use a spherical Navarro-Frenk-White (NFW) profile (Navarro et al., 1996). We treat the central black hole as a point source. The major complication that triaxial galaxies do introduce is that three viewing angles must be taken into account. Axisymmetric galaxies by comparison only have two viewing angles. For the Schwarzschild modelling runs, the viewing angles are treated as free parameters whose values are to be determined.

Our mock galaxies and their observations are taken from the Illustris cosmological simulations (Vogelsberger et al., 2014a, b; Genel et al., 2014; Nelson et al., 2015). The Illustris project is a series of large-scale hydrodynamical cosmological simulations of galaxy formation, using advanced moving mesh code AREPO (Springel, 2010) with comprehensive physical process modelling. Star formation is included (e.g., Springel & Hernquist, 2003; Schaye & Dalla, 2008) as well as primordial and metal-line cooling processes (e.g., Katz et al., 1996; Wiersma et al., 2009a), stellar evolution and chemical enrichment procedures (e.g., Wiersma et al., 2009b; Few et al., 2012), stellar feedback (e.g., Navarro & White, 1993; Springel & Hernquist, 2003) and AGN feedback (e.g., Springel et al., 2005; Di Matteo et al., 2005). With these physical process models, Illustris is able to produce galaxies sufficiently realistic for our purposes. We choose nine early type galaxies for our Schwarzschild evaluation. The morphologies are deliberately selected such that we have three galaxies close to being oblate, three close to prolate, and three definitely triaxial. For each galaxy, we construct IFU-like mock brightness and kinematic data at five different sets of viewing angles.

The cosmological parameters adopted in the Illustris simulations and thus implicit in our Schwarzschild modelling are: $\Omega_{\rm m}=0.2726$, $\Omega_{\rm\Lambda}=0.7274$, $\Omega_{\rm b}=0.0456$, $\sigma_{8}=0.809$, $n_{\rm s}=0.963$, and $H_{0}=\rm 70.4\ km\cdot s^{-1}\cdot Mpc^{-1}$ (Hinshaw et al., 2013).

On a per galaxy basis, the steps we take are

calculate the true values for the galaxy properties we are interested in from the simulation stellar particle data, 2. 2.

construct the mock data, at 5 different sets of viewing angles, to be used as constraints in the Schwarzschild models, 3. 3.

for each set of mock data

1. (a)

execute a series of Schwarzschild models varying the free parameters, 2. (b)

select the best-fitting model matching the constraint data, 3. (c)

from the best-fitting model, calculate the model values for the galaxy properties, 4. (d)

compare the model property values with the true property values, 4. 4.

collate the results for the five viewing angle sets and assess what has been achieved.

The method we use for varying the free parameters is an optimised grid search as described in Zhu et al. (2018a). No amendments to the method or its associated software are required.

3 The Schwarzschild’s method

The Schwarzschild’s method is a flexible orbit-superposition method for building dynamical models of galaxies in dynamical equilibrium and was first proposed by Schwarzschild (1979). Initially, the method was used to create spherical and axisymetric models (e.g., Richstone, 1980, 1982, 1984; Richstone & Tremaine, 1984; Levison & Richstone, 1985). Schwarzschild (1982, 1993) extended the method to more general triaxial shapes. The method has been widely used in dynamical modelling (e.g., Merritt & Fridman, 1996; Zhao, 1996; Rix et al., 1997; Siopis & Kandrup, 2000; Häfner et al., 2000; van de Ven et al., 2006; Capuzzo-Dolcetta et al., 2007; van den Bosch et al., 2008; Zhu et al., 2018a, b), and has been applied to more complicated barred galaxies (e.g., Wang et al., 2012, 2013; Vasiliev & Athanassoula, 2015). The steps we take in modelling a galaxy with Schwarzschild’s method are indicated below, and follow the approach in van den Bosch et al. (2008).

Creating gravitational potentials for our galaxies is straightforward. Our potentials consist of a stellar potential, a dark matter potential and also a potential generated by a central black hole. For the stellar component, we use the Multi-Gaussian Expansion (MGE) formalism (Cappellari, 2002). Within the stellar MGE, we use the triaxial viewing angles to deproject the surface brightness into a three-dimensional luminosity density, and use a constant mass-to-light ratio to turn the density into a mass density from which the potential can be determined.

We create initial conditions for our orbits by sampling from the three integrals of motion defining a triaxial system. Taking the integrals as energy $E$, second integral $I_{2}$, and third integral $I_{3}$, these integral triplets ($E,I_{2},I_{3}$) play the same role as ($E,L_{z},I_{3}$) in an axisymmetric system (Binney & Tremaine, 2008). In this paper, we take $21\times 10\times 7$ combinations of ($E,I_{2},I_{3}$) to form initial conditions for our orbits.

Box orbits are crucial for supporting the triaxial shape, so we add further initial conditions specifically for box orbits (Schwarzschild, 1993; van den Bosch et al., 2008). As box orbits always touch equipotentials (Schwarzschild, 1979), we construct our initial conditions on equipotential surfaces using energy $E$ and two spherical angles ($\theta,\varphi$) giving a further $21\times 10\times 7$ orbits. Since retrograde stars may exist in elliptical (Bender, 1988) and lenticular galaxies (Kuijken et al., 1996), we also produce initial conditions which will result in counter-rotating orbits.

In total we have 3 sets of $21\times 10\times 7$ orbits giving 4410 orbits in our orbit library. The library comprises a typical set of ($E,I_{2},I_{3}$) orbits, a set of counter-rotating orbits ($E,-I_{2},I_{3}$), and the box orbits ($E,\theta,\varphi$).

As in Cappellari et al. (2006) and van den Bosch et al. (2008), we dither every orbit to give $5^{3}$ orbits by perturbing the initial conditions slightly. Once the orbit trajectories have been created the trajectories will be co-added to form a single orbit trajectory in our orbit library. The rationale behind dithering orbits is that the mechanism should reduce the Poisson noise in calculating model observables. In creating the trajectories, we integrate along the orbits in our gravitational potential for 200 orbital periods and aim to achieve a relative accuracy of $10^{-5}$ in energy conservation.

In arriving at the orbit weights, we use surface brightness and luminosity density from a galaxy’s MGE together with Gauss-Hermite coefficients (van der Marel & Franx, 1993; Gerhard, 1993) of the line-of-sight velocity distribution as constraints. Schwarzschild’s method weights orbit contributions to model observables in an attempt to match the observed values. In so doing the model and observed values are divided by the observational error so that a $\chi^{2}$ comparison is achieved. The weights themselves are determined by the van den Bosch et al. (2008) implementation using the Lawson & Hanson (1974) non-negative least squares (NNLS) implementation.

4 Formulas and Parameters

In this section, we introduce the key formulas and parameters used in our modelling. $\S$ 4.1 to $\S$ 4.4 deal with the gravitational potential; $\S$ 4.5, orbits; and $\S$ 4.6, error analysis.

4.1 Stellar potential

We construct our stellar potential by using the Multi-Gaussian Expansion (MGE) formalism, the starting point for which is to model a galaxy’s surface brightness by a sum of two-dimensional elliptical Gaussians. Bendinelli (1991) used MGEs for the reconstruction of galaxy images, but only for spherical systems. Monnet et al. (1992) generalized MGEs to triaxial systems and made their application to real galaxies possible. Further development of the formalism and its application was undertaken by Emsellem et al. (1994a, b). Cappellari (2002) develop a sectors fitting code to perform MGE fits to galaxy images, and this code made MGEs more readily accessible.

The projected surface brightness can be written as a sum of two-dimensional Gaussians

[TABLE]

with

[TABLE]

where $N$ is the number of Gaussians, the subscript $j$ means the value for $j$-th Gaussian, $L_{j}$ is the total luminosity, $q_{j}^{\prime}=b_{j}^{\prime}/a_{j}^{\prime}$ is the axis ratio and $\sigma_{j}^{\prime}$ is the dispersion along major axis. $(R^{\prime},\theta^{\prime})$ are the polar coordinates on the plane of the sky $(x^{\prime},y^{\prime})$. $\psi$ is the position angle (PA), which is the angle between the Gaussian major axis and $y^{\prime}$ axis, measured counterclockwise. Here we note that we do not have subscript $j$ in $x$ and $y$ as Equation (1) in Cappellari (2002) and Equation (4) in van den Bosch et al. (2008), because we assume the Gaussians are concentric with their major axes aligned. By using the deprojection method given in Emsellem et al. (1994a), we can obtain the three-dimensional luminosity distribution (see Equation 8 in Emsellem et al., 1994a)

[TABLE]

where $p_{j}=b_{j}/a_{j}$ and $q_{j}=c_{j}/a_{j}$ are the axis ratios for the ellipsoid ($a_{j}\geq b_{j}\geq c_{j}$), while other terms have the same meaning as Equation (1). The combination between Equation (1) and Equation (3) are described in $\S$ 4.4. We use a constant stellar mass to light ratio $M_{*}/L$ so that the stellar mass density is given by

[TABLE]

The ratio $M_{*}/L$ is a free parameter in our modelling. The corresponding gravitational potential $\Phi_{\rm star}$ is calculated based on the classical Chandrasekhar (1969) formula (see $\S 3.8$ in van den Bosch et al., 2008 for details).

4.2 NFW dark matter halo

We use the spherical Navarro-Frenk-White (NFW) profile described in Navarro et al. (1996) for creating our dark matter haloes. Using the profile, the density of dark matter can be written as

[TABLE]

where $\rho_{0}$ and the “scale radius” $R_{\rm s}$ are two free parameters. By solving Poisson’s equation, the potential of NFW profile is

[TABLE]

The virial radius $R_{200}$ is defined as the radius within which the average density is 200 times the critical density $\rho_{\rm crit}$. As the total enclosed mass within $R_{200}$ is $M_{200}$, we have

[TABLE]

with

[TABLE]

where $H_{0}$ is the Hubble constant and $G$ is the gravitational constant.

From the definition of density profile in Equation 7, we obtain

[TABLE]

Instead of $\rho_{0}$ and $R_{\rm s}$, we use the “concentration parameter” $c$ and “the fraction of dark matter within $R_{\rm 200}$” $f$ to represent a NFW halo

[TABLE]

where $M_{*}$ is the total stellar mass.

Based on equations  7 to  10, the parameters in the NFW profile can be rewritten as

[TABLE]

[TABLE]

The two parameters $c$ and $f$ of the NFW profile are free parameters in our modelling.

4.3 Black hole

We assume a central black hole can be represented by a Plummer potential (van den Bosch et al., 2008):

[TABLE]

where $r_{\rm soft}$ is the softening length, which is introduced to prevent the central potential to be infinite, and equals to 0.001 arcsec in our modelling runs. All our galaxies have central black holes but, given the nature of our research, they will not strongly influence our results. We choose therefore to only use a token black hole which also helps in reducing our computing costs. We therefore fix the black hole mass $M_{\rm BH}=10^{6}M_{\odot}$ for all mock data sets.

4.4 Intrinsic shapes and viewing angles

A two-dimensional Gaussian has two further parameters in addition to luminosity: the deviation along major axis $\sigma_{j}^{\prime}$ and the axis ratio $q_{j}^{\prime}$. For a three-dimensional Gaussian, it becomes $\sigma_{j}$ and axis ratios $(p_{j},q_{j})$. Here we use the compression factor $u_{j}=\sigma_{j}^{\prime}/\sigma_{j}$ to present the difference between the length of the intrinsic major axis and the major axis on the projected plane. For a triaxial galaxy, we set the major axis to be the $x-$axis in Cartesian coordinates. There are three viewing angles to be considered, angles in spherical coordinates $(\theta,\varphi)$ and the rotation $\psi$ of the projected major axis on the sky-plane, measured counterclockwise from the $y^{\prime}-$axis. The original coordinates $(x,y,z)$, the line of sight coordinates $(x^{\prime},y^{\prime},z^{\prime})$ and the projected coordinates after rotation $(x^{\prime\prime},y^{\prime\prime},z^{\prime\prime})$ are related by

[TABLE]

[TABLE]

After assuming the viewing angles $(\theta,\varphi,\psi)$ in models, we can covert observed quantities $(\sigma_{j}^{\prime},q_{j}^{\prime})$ into intrinsic parameters $(\sigma_{j},p_{j},q_{j})$, based on formulas below (van den Bosch et al., 2008, Cappellari, 2002)

[TABLE]

[TABLE]

[TABLE]

where $\delta_{j}^{\prime}=1-q_{j}^{\prime 2}$.

In reverse, if we know the intrinsic shapes $(p,q)$ combined with the compression factor $u$ of any deprojected Gaussian along with the projected quantity $q^{\prime}$, we can obtain the space orientation parameters $(\theta,\varphi,\psi)$

[TABLE]

[TABLE]

[TABLE]

where the parameters should satisfy $q\leq p\leq 1$, $q\leq q^{\prime}$ and $\max(q/q^{\prime},p)\leq u\leq\min(p/q^{\prime},1)$ for validation of the equation.

In the modelling, we take $(p,q,u)$ as free parameters instead of $(\theta,\varphi,\psi)$ as they directly represent the intrinsic morphologies of galaxies. We use $p$ and $q$ in our definition of triaxiality (Binney & Tremaine, 2008) and take

[TABLE]

4.5 Orbit distribution parameters

We introduce three parameters to characterize orbits, time-averaged radius $r$, circularity $\lambda_{z}$ and $\lambda_{x}$. We define orbit circularity $\lambda_{z}$ as in Zhu et al. (2018a, b) as a ratio of time-averaged quantities

[TABLE]

where $\overline{L_{z}}=\overline{xv_{y}-yv_{x}}$, $\overline{r}=\overline{\sqrt{x^{2}+y^{2}+z^{2}}}$ and $\overline{V_{\rm rms}}=\sqrt{\overline{v_{x}^{2}+v_{y}^{2}+v_{z}^{2}+2v_{x}v_{y}+2v_{x}v_{z}+2v_{y}v_{z}}}$. The numerator is the the time averaged z-component of the the orbit’s angular momentum, and the denominator is the product of two time-averaged root mean square (rms) quantities calculated from the orbit, the average radial position and velocity. If we treat the rms velocity as a circular velocity in the equatorial plane, then the product represents the angular momentum of a typical circular orbit associated with the original orbit. The ratio of the two angular momentum terms gives us our measure of orbit circularity. Note that it is the circularity distribution as a function of radius which is important for comparison purposes. Whilst there is a physical explanation for our orbit circularity definition, it can not be justified precisely mathematically, and should be treated as empirical. Thus $\lambda_{z}\sim 1$ represents highly rotating short-axis tube orbits, while $\lambda_{z}\sim 0$ is mostly long-axis tube and box orbits. Since there are prolate-like galaxies in our sample and $\lambda_{z}$ can not distinguish long-axis tube orbits from box orbits, we also define $\lambda_{x}$, which represents the normalized angular momentum around the major axis ($x$-axis)

[TABLE]

where $\overline{L_{x}}=\overline{yv_{z}-zv_{y}}$. The orbit circularity distribution can be thought of the probability densities of orbits on the coordinate space $\lambda_{z}$ versus $r$ or $\lambda_{x}$ versus $r$.

Velocity dispersion anisotropy parameters are widely used as indicators of the underlying orbit distribution of a galaxy. The velocity anisotropy $\beta_{r}$ (Binney & Tremaine, 2008) is defined as

[TABLE]

where $\sigma_{r}$, $\sigma_{\theta}$ and $\sigma_{\varphi}$ are the three components of the velocity dispersion in spherical coordinates, and $\sigma_{t}^{2}=\sigma_{\theta}^{2}+\sigma_{\varphi}^{2}$. $\beta_{r}>0$ indicates radial anisotropy and $\beta_{r}<0$ indicates tangential anisotropy. However, $\beta_{r}$ ranges from negative infinity to 1, and is not uniformly distributed in the meaning. In order to overcome this shortcoming of $\beta_{r}$ for quantitative analysis, we define a new anisotropy parameter the “tangential fraction” ($0\leq f_{t}\leq 1$), as

[TABLE]

$f_{t}>0.5$ represents tangential anisotropy while $f_{t}<0.5$ represents radial anisotropy, and $f_{t}=0.5$ is for isotropy.

4.6 $\chi^{2}$ and the best-fitting model

We have a total of six free parameters when we run models. They are intrinsic shape parameters $p$ , $q$ and $u$ ($\S$ 4.4), NFW dark matter halo parameters $c$ and $f$ ($\S$ 4.2) and the stellar mass to light ratio $M_{*}/L$ ($\S$ 4.1). We use both kinematic data and luminosity distributions as constraints. The $\chi_{\rm NNLS}^{2}$ function that needs to be minimized by NNLS contains two components:

[TABLE]

For the luminous component, we use both the two-dimensional surface brightness $S_{n}$ ($n$-th Voronoi bin in the observing plane) and the three-dimensional luminosity distribution $\rho_{m}$ ($m$-th Schwarzschild cell in 3D space) as constraints. Following (Zhu et al., 2018a), we set relative errors of $S_{n}$ to be $1\%$ and $\rho_{n}$ to be $2\%$, thus

[TABLE]

where the variables marked with a “$*$” indicates model predictions and those without a “$*$” are observational data which are constructed from the MGE formulae.

For kinematic components, we use Gauss-Hermite series (van der Marel & Franx, 1993; Gerhard, 1993; Rix et al., 1997) to fit the velocity distribution profile in each bin. Considering Gauss-Hermite coefficients $h_{1}$ and $h_{2}$, we have

[TABLE]

where the markers “$\Delta$” represent observational errors. The kinematics from the model are luminosity-weighted in the same way as the observations. After modelling, $h_{1}$, $h_{2}$, $\Delta h_{1}$ and $\Delta h_{2}$ will be converted to $V$, $\sigma$, $\Delta V$ and $\Delta\sigma$ for comparison purposes (Zhu et al., 2018a).

We explore the parameter space by changing the values of our free parameters and reconstructing our gravitational potential accordingly. The modelling begins from assumed start parameters and iterates towards models with smaller chi-square values in fixed steps. The differences between model values and observed values are calculated as

[TABLE]

where $V^{*}_{n}$, $\sigma^{*}_{n}$ are model values for $n$-th bin and $V_{n}$, $\sigma_{n}$ mean observed values, while $\Delta V_{n}$, $\Delta\sigma_{n}$ represent observational errors. We use $\chi^{2}$ (Equation 30) to search for the best-fitting model. Considering Equation 27, $\chi^{2}_{\rm lum}$ is usually small, $\chi^{2}$ is in fact strongly correlated with $\chi^{2}_{\rm NNLS}$. Typically, a completed model with six free parameters require $1000\sim 2000$ Schwarzschild modelling runs depending on how close the start parameter values are to the best-fitting model. We normalize $\chi^{2}$ by forcing $\chi^{2}_{\rm min}=N_{\rm obs}$, thus the normalized $\chi^{\prime 2}=\chi^{2}/\chi^{2}_{\rm min}\times N_{\rm obs}$.

In van den Bosch et al. (2008), the $1\sigma$ confidence level is defined as

[TABLE]

where $N_{\rm obs}=2\times N_{\rm kin}$ when we only use $V$ and $\sigma$ as model constraints. $N_{\rm par}$ means the number of free parameters and equals to 6 here, while $N_{\rm kin}$ is usually around several hundreds or one thousand in our mock data sets.

We can see with the increase of $N_{\rm obs}$, $\Delta\chi^{\prime 2}/\chi^{\prime 2}$ goes down, which means $\sqrt{2(N_{\rm obs}-N_{\rm par})}$ could be smaller than the model $\chi^{\prime 2}$ fluctuation thus underestimate the $1\sigma$ region when $N_{\rm obs}$ is large. In order to solve this problem, we try to re-scale $1\sigma$ confidence level by

[TABLE]

where $A$ is a constant parameter to be determined. For each mock galaxy, we perturb the kinematic maps $(V_{i},\sigma_{i})$ with their error maps $(\Delta V_{i},\Delta\sigma_{i})$: $V^{\prime}_{i}=V_{i}+\Delta V_{i}\times a$, $\sigma^{\prime}_{i}=\sigma_{i}+\Delta\sigma_{i}\times b$, where $a$, $b$ are random values with standard deviation equals unity. By fixing the gravitational potential, we re-solve the orbits weight with NNLS, and its result $\chi^{\prime 2}_{\rm perturb}$ is slightly different with $\chi^{\prime 2}$. We do the perturbations and find orbit weight solutions for 100 times. The standard deviation of these 100 $\chi^{\prime 2}_{\rm perturb}$ is taken as the model $\chi^{\prime 2}$ fluctuation. We find the parameter $A$ ranges from $1.03$ to $3.19$ for different mock galaxies, with an average of $2.01$ (this is consistent with the CALIFA model fluctuations in Zhu et al., 2018a). Thus we roughly take $A=2$ and Equation 32 become

[TABLE]

which will be used as the $1\sigma$ confidence level in our models constrained by MaNGA/CALIFA quality data.

The $1\sigma$ confidence level we define here represents the model fluctuations, and we take it as the modelling error in our results analysis. A standard $\chi^{2}$ analysis is not valid here for two reasons. First, for our models, the number of active orbits could be larger than the number of data points used as model constraints, so that the degrees of freedom $N_{\rm kin}-N_{\rm orb}$ in a $\chi^{2}$ distribution will be negative (Press et al., 1992). Second, the typical $\Delta\chi^{2}$ of the standard $\chi^{2}$ analysis is equal to 7.04 (joint constraints on 6 free parameters simultaneously) for 1 $\sigma$ confidence level, while the model fluctuation is relatively large ($\Delta\chi^{2}\sim 130$ for $N_{\rm obs}\sim 2000$). That means if we use the standard $\chi^{2}$ analysis, $\Delta\chi^{2}$ (errors of the model) will be significantly under-estimated. Using model fluctuations as confidence levels is an empirical way of calculating the error bars. We refer the readers to Morganti et al. (2013), who applied the M2M method to NGC 4494, for detailed discussions on the approach and justifications.

5 Simulated galaxies and Mock data

5.1 Simulated galaxies

We use the highest resolution simulation, Illustris-1, in the Illustris project (Vogelsberger et al., 2014a, b; Genel et al., 2014; Nelson et al., 2015), which contains $2\times 1820^{3}$ resolution elements in a volume of $(106.5\ \rm Mpc)^{3}$, with $6.26\times 10^{6}M_{\odot}$ dark matter mass resolution and an initial baryonic mass resolution of $1.26\times 10^{6}M_{\odot}$.

We select nine early-type Illustris galaxies with stellar mass ranging from $10^{10.74}$ to $10^{12.05}M_{\odot}$. Since effective radius $R_{\rm e}$ is usually calculated from a galaxy’s surface brightness, its value may change with viewing angles. We use the average effective radius of five different projections $\overline{R_{\rm e}}$ as the typical size of a galaxy. We calculate the intrinsic shapes of these galaxies by using the reduced inertia tensor method (Allgood et al., 2006). For the stellar components, these nine galaxies have different morphologies and comprise three oblate galaxies subhalo73666, subhalo347122, subhalo336920 (galaxy O1, O2, O3 for short) with their triaxial parameter values at $\overline{R_{\rm e}}$, $T_{\rm e}$, typically $\leq 0.3$; three triaxial galaxies subhalo312924, subhalo175435, subhalo16942 (galaxy T1, T2, T3 for short) with $0.3<T_{\rm e}<0.7$; and three prolate galaxies subhalo73664, subhalo191220, subhalo222715 (galaxy P1, P2, P3 for short) with $T_{\rm e}\geq 0.7$. We refer to them as oblates, triaxials and prolates. The dark matter haloes in these galaxies are triaxial with short to long axis ratios between 0.63 and 0.92, and medium to long axis ratios between 0.70 and 0.95 at 30 kpc. Detailed information on the sample galaxies is shown in Table 1.

5.2 Mock data

The Illustris galaxies contain both stellar particles and dark matter particles with a resolution of $\sim 10^{6}M_{\odot}$, as mentioned in $\S$ 5.1. We project the stellar particles in each galaxy onto five different observational planes by setting specific viewing angles. The angles $(\theta,\varphi)$ are randomly chosen in spherical coordinates, with probability $dP=d\Omega/4\pi=\sin\theta d\theta d\varphi/4\pi$, where $\Omega$ represents the solid angle. Due to symmetry, we restrict $0\leq\theta\leq\pi/2$ and $0\leq\varphi<\pi$. The angles $(\theta,\varphi)$ we have chosen for the mock data are shown in Table 1. We create surface brightness and kinematic maps, which are used to generate model constraints, for each simulated galaxy with the chosen viewing angles $(\theta,\varphi)$. Making allowances for spatial differences, the way we generate luminosity input and kinematic input follows the method in Li et al. (2016). We create $5\times 9=45$ mock data sets in total, and treat each as an independent observed galaxy in the modelling.

The surface mass density $\Sigma^{*}(n)$ is calculated with a pixel size of 0.5kpc$\times$0.5kpc. Instead of using the luminosity distribution given in the simulation directly, we use the stellar mass distribution to construct mock galaxies with a constant stellar mass-to-light ratio $M_{*}/L=5$. We can obtain the surface brightness $\Sigma(n)=\Sigma^{*}(n)/5$ and fit it with the MGE method to generate the luminosity input (see $\S$ 4.1). Finally, we rotate the major axis of surface brightness to be horizontal in the observational plane, ensuring the other viewing angle $\psi\approx 90^{\circ}$ (see $\S$ 4.4).

Kinematic data can be obtained from the velocity distribution of stellar particles, also in a 0.5 kpc$\times$0.5 kpc pixels. In this paper, we use the luminosity-weighted mean velocity $V$ and velocity dispersion $\sigma$ (standard deviation of line-of-sight velocity) together with their observational errors, converted to Gauss-Hermite coefficients $h_{1}$ and $h_{2}$, as constraints. In order to achieve consistency in approach between the simulation tests and real observations, we do not model with skewness $h_{3}$ and kurtosis $h_{4}$ as they are not always able to be determined well in real IFU surveys like MaNGA. We bin the pixels to nearly constant signal-to-noise ratio $S/N=20$ per bin using a Voronoi 2D-binning method (Cappellari & Copin, 2003). By projecting the stellar particles onto the observing plane and obtaining the line of sight velocities of stellar particles located in each bin, we calculate the luminosity-weighted mean velocity $V$ and velocity dispersion $\sigma$ by fitting the binned data with a Gaussian. Observational errors are estimated by using the bootstrap method (random sampling with replacement). We randomly select half the number of the stellar particles in each bin with replacement 500 times, and calculate $V_{1},V_{2},...V_{500},\sigma_{1},\sigma_{2},...,\sigma_{500}$ for our selections. The standard deviation of these 500 mean velocities and velocity dispersions are taken as observational errors $\Delta V$ and $\Delta\sigma$. As triaxial models are point-symmetric, we lastly point-symmetrize the kinematic data in the Voronoi bins.

We take galaxy O1 with viewing angles $(\theta,\varphi)=(83,131)^{\circ}$ as an example in Figure 1 and Figure 2. Figure 1 shows the surface brightness distribution and the contours obtained from MGE fitting. Figure 2 presents the kinematic inputs and show the observational mean velocity $V_{\rm obs}$, and velocity dispersion $\sigma_{\rm obs}$ together with their errors $V_{\rm obserr}$, $\sigma_{\rm obserr}$. For ease of comparsion we also show the model estimated $V_{\rm model}$, $\sigma_{\rm model}$ together with the standardized residuals between the inputs and outputs $V_{\rm stdres}=(V_{\rm model}-V_{\rm obs})/V_{\rm obserr}$, $\sigma_{\rm stdres}=(\sigma_{\rm model}-\sigma_{\rm obs})/\sigma_{\rm obserr}$. Note that we use units of arcsec in Figure 1 and Figure 2 by assuming the galaxy is located at the redshift of $z=0.03$.

6 Results

In this section, we show model estimates of mass distributions and orbit distributions, including enclosed mass profiles in $\S$ 6.1, galaxy morphologies in $\S$ 6.2, circularity distributions $f(r,\lambda_{z})$, $f(r,\lambda_{x})$ in $\S$ 6.3 and anisotropy parameters $(\beta_{r},f_{t})$ in $\S$ 6.4. In addition, we report on a simple test concerning the influence of the initial conditions $(E,I_{2},I_{3})$ in $\S$ 6.5.

6.1 Mass profile

Mass estimation is the first capability, in line with our objectives, that we consider. We estimate the mass variation with spherical radius $r$ as it is independent of the morphology assumption. Here we compare the total mass $M_{\rm tot}(\textless r)$, the stellar mass $M_{*}(\textless r)$ and the dark matter mass $M_{\rm dark}(\textless r)$.

Figure 3 shows the true mass profile $M(\textless r)$ calculated from the simulation and the profile obtained from modelling. The blue lines and symbols represent the total mass, the red lines and symbols the stellar mass, while the green ones are for dark matter. We plot the true profiles with solid lines, and use dashed lines with different markers to indicate the model results from mock data with different viewing angles $(\theta,\varphi)$. The blue triangles with error bars in the top of each panel represent the typical error bars ($1\sigma$ confidence levels, see Equation 33) of the model estimated total mass. The dots in the bottom of each panel show the relative average difference $1-\langle M_{\rm true}(\textless r)/M_{\rm model}(\textless r)\rangle$ between true values and model estimated values from our mock data with five different projections, and the corresponding error bars represent average confidence levels within $1\sigma$. The vertical dashed line represents $\overline{R_{\rm e}}$ of each galaxy and the value of maximum data coverage $R_{\rm m}$, which means the points outside this value are not reliable, is shown as text. As seen in Figure 3, the estimation of total mass within $\overline{R_{\rm e}}$ is good for all galaxies (average bias $\textless 15$ percent) except P2. However the stellar mass is systematically underestimated and the dark matter mass overestimated, especially in the inner parts of galaxies.

In order to quantify the difference between true profiles and model results, we compare in Figure 4 the enclosed mass within $\overline{R_{\rm e}}$ with the true values of all 45 mock observations. The blue markers represent different viewing angles, from face-on to edge-on are circles, triangles, squares, diamonds, pentagrams. The red triangles with error bars have the same meaning as the error bars in Figure 3. For each galaxy except P2, the estimation of total mass $M_{\rm tot}(\textless\overline{R_{\rm e}})$ is satisfactory, with average relative deviations being inside $\pm 15$ percent region. The stellar mass $M_{*}(\textless\overline{R_{\rm e}})$ is $\sim 24$ percent lower than the true value on average and dark matter mass $M_{\rm dark}(\textless\overline{R_{\rm e}})$ is $\sim 38$ percent higher.

The potential in our model is generated by both the stellar mass and dark matter mass. The stellar kinematics, in principle, can only constrain the total mass profile $M_{\rm tot}(r)$. The only way to separate stellar mass and dark matter is based on the difference in their density slopes $\gamma$ (see Figure 5) . Since we use the NFW profile to construct the dark matter potential, the inner density slopes of our model results in $\gamma$ approximate to $-1$, while the true profiles are more likely to be “cored” ($\gamma=0$). Thus we have good recovery for total mass, but the stellar mass and dark matter mass estimates have larger uncertainties due to our choice of dark matter model and to the degeneracy between them. By changing the NFW dark matter profile to the generalized NFW (gNFW, Cappellari et al., 2013), whose inner density slopes are not fixed, we re-run 9 mock data sets (one for each galaxy). The average underestimation of stellar mass decrease from $\sim 24$ to $\sim 13$ percent and the average overestimation of dark matter goes down from $\sim 38$ to $\sim 18$ percent.

In Figure 4, we show comparisons between our Schwarzschild modelling and an oblate JAM technique (Li et al., 2016). We use the same mock data as input to the two different methods. The green and dark symbols represent JAM results, corresponding to the red and blue markers for the results from Schwarzschild modelling. The recovery of total mass $M_{\rm tot}(\textless\overline{R_{\rm e}})$ is also good for JAM modelling. The scatter of relative deviations of all 45 mock data sets are $18$ percent for Schwarzschild modelling and $16$ percent for JAM modelling, which means these two different methods have similar abilities to recover the total mass. However, JAM modelling gives very small error bars, which seems to be optimistic. JAM modelling can yield more narrowly constrained parameters purely because of more restrictive and sometimes ad-hoc assumptions, for example on the velocity anisotropy which is left completely free in Schwarzschild modelling. For the stellar mass and dark matter mass, both techniques give large uncertainties. Using a gNFW profile for dark matter, JAM modelling seems to underestimate dark matter mass and overestimate stellar mass slightly.

Our mass accuracy is consistent with Thomas et al. (2007). They construct axisymmetric Schwarzschild models for N-body merger remnants using a different orbit sampling method (Thomas et al., 2004). They find the total mass has a bias of between 3 to 20 percent for edge-on views and up to 50 percent for face-on views around $1\sim 2R_{\rm e}$. They also adopt the NFW profile for dark matter and find a systematic underestimation for stellar mass: $(M_{*}/L)_{\rm model}/(M_{*}/L)_{\rm true}$ is between 0.5 and 0.9.

6.2 Morphology

The different intrinsic shapes of the nine simulated galaxies make them a good sample with enough variations to assess estimation of triaxiality. Here we use the axis ratios $p=b/a$, $q=c/a$ and corresponding triaxial parameter $T=(1-p^{2})/(1-q^{2})$ to describe galaxy morphology, where $a$, $b$ and $c$ are the major, intermediate and minor axes of the three-dimensional luminosity distribution of the galaxy. Since the stellar potential is generated by a superposition of three-dimensional Gaussians (see $\S$ 4.1), the equipotential planes are usually not regular ellipsoids. Since we choose to ignore any misalignment between different Gaussians, we can simply calculate the model’s intrinsic shapes at different radii by applying the least-square ellipse fitting method to luminous contours on projected planes. We thus obtain $p$ as a function of radius $r$, $p(r)$, from the $x$-$y$ plane and obtain $q(r)$ from the $x$-$z$ plane.

The variation of intrinsic shapes with radius for all nine galaxies is shown in Figure 6. We show variations of intrinsic shapes from $0.5\overline{R_{\rm e}}$ to $2\overline{R_{\rm e}}$. The black dots are for true shapes, while the blue, magenta, green, orange, red dots represent the shapes recovered from mock data of face-on to edge-on views. The size of the dots indicate the radius: larger symbols represent larger radii. For oblate and triaxial galaxies, the model estimates improve with increasing inclination $\theta$. The model shapes constrained by edge-on mock data (red) are close to the real morphology (black). This is because we know the real kinematic information well if we observe galaxies dominated by minor-axis rotation in edge-on views, and the modelling thus gives a better restriction for the axis ratio $q$ (see Equation 21). For the prolate galaxy P3, we find that the model estimate with an edge-on view is quite different. For the viewing angles $(\theta,\varphi)=(88,77)^{\circ}$, we are nearly observing this galaxy along the major-axis. We therefore are unable to obtain enough information about the major-axis rotation, which dominates prolate galaxies. For all the galaxies, we find that some models tend to overestimate $p$ and $q$, which means they are slightly biased towards sphericals.

We show a quantitative comparison of the true and model values for intrinsic shapes and the triaxial parameter at $\overline{R_{\rm e}}$ in Figure 7. The green dots represent the true values calculated from simulation directly, and the circles, triangles, squares, diamonds and pentagrams blue markers are for those estimated from mock data face-on to edge-on views. The red triangles with error bars show the mean values of results from mock data with five different projections and average confidence levels within $1\sigma$. For the axis ratio $p$, four of nine galaxies (O2, O3, T1, P2) give relatively good recoveries, while four galaxies overestimate $p$ by $0.13\sim 0.18$ and one galaxy underestimate $p$ by $0.10$. For the axis ratio $q$, only the galaxy P2 shows a good recovery, while other galaxies have average overestimations between $0.10$ and $0.26$. The mean bias of these two parameters are $\Delta p=0.07$ and $\Delta q=0.14$. Although there are some systematic biases of $p$ and $q$, we can still see a trend in the distribution of triaxial parameter $T$: the model recovered $T_{\rm e}$ rises as the true $T_{\rm e}$ rises but with a shallower slope. For galaxies O1, O2, O3 and T1, $T_{\rm e}$ is overestimated, while it is underestimated for T2, T3, P1, P2, P3. This means that on average the modelling tends to make galaxies more triaxial.

6.3 Circularity

As mentioned in $\S$ 4.5, the $\lambda_{z}$ circularity parameter is used to characterize the orbits. In the modelling, we calculate the circularity parameter for each orbit. Based on the orbit weights from the best-fitting model, we obtain the probability density distribution on the 2D coordinate space $f(r,\lambda_{z})$. In the simulations, the orbital circularity of each stellar particle is obtained through running 20-orbit revolutions and recording time-averaged quantities of $L_{z}$, $r$ and $V_{\rm rms}$. In order to do so, we freeze the gravitational potential of the galaxy at the current snapshot and use the phase space information of given stellar particles as their initial condition. A second-order leapfrog integrator is adopted to update the positions and velocities during orbital evolution (Dehnen & Read, 2011). Thus we can also obtain the distribution $f(r,\lambda_{z})$ from the simulations.

We show the comparison of this distribution for an example galaxy T1 in Figure 8. The left panel shows the true circularity distribution while the right panel shows the model estimate for the mock data set with viewing angles $(70,47)^{\circ}$. Each map is equally divided into $n_{r}\times n_{\rm\lambda}=10\times 21$ rectangle bins. Each colored pixel represents the probability density of orbits in the bins, from blue (low density) to red (high density). The maps are smoothed and the color bar is linear. In order to further quantify the model estimate, we divide orbits into different components based on $f(r,\lambda_{z})$. Following Zhu et al. (2018b), we classify orbits with different $\lambda_{z}$ to be: dynamically cold ($\lambda_{z}\geq 0.8$), dynamically warm ($0.25<\lambda_{z}<0.8$), dynamically hot ($-0.25\leq\lambda_{z}\leq 0.25$) and counter-rotating ($\lambda_{z}<-0.25$) components. From our models, the hot orbits mainly consist of box orbits and long-axis tube orbits, while the warm and cold orbits are short-axis tube orbits. These four different components are indicated by the horizontal dashed lines in the maps. We can see the model orbit distribution $f(r,\lambda_{z})$ is a reasonable match to the true distribution.

In Figure 8, there is a clear peak around $\lambda_{z}\sim 0$, which comprises mainly box orbits and long-axis tube orbits that can not be separated in the distribution of $\lambda_{z}$. For the purpose of distinguishing them, we also need to examine the distribution of $\lambda_{x}$, which quantifies rotation about the major axis. We calculate both $\lambda_{z}$ and $\lambda_{x}$ for all active orbits within $\overline{R_{\rm e}}$ for the same example galaxy and plot the probability density distributions in Figure 9. The top left panel shows the true distribution while the top right panel gives the model estimate. The bottom panel indicate how we separate orbits with different circularities. In addition to separating cold, warm, hot and counter-rotating components based on $\lambda_{z}$, we also separate the hot orbits based on $\lambda_{x}$: prograde long-axis tube components ($\lambda_{x}>0.25$), counter-rotating long-axis tube components ($\lambda_{x}<-0.25$), box components ($\left|\lambda_{x}\right|\leq 0.25$) and slowly-rotating components ($\left|\lambda_{x}\right|,\left|\lambda_{z}\right|\leq 0.25,\left|\lambda_{x}\right|$ or $\left|\lambda_{z}\right|>0.05$).

We show a comparison of cold, warm, hot, and counter-rotating orbit fractions within $\overline{R_{\rm e}}$ in Figure 10. The magenta circles, orange squares and red triangles with error bars represent the mean values and $1\sigma$ uncertainties for oblates, triaxials and prolates. The blue symbols with different markers represent model values obtained from mock data with different viewing angles $(\theta,\varphi)$. The solid lines mean equal values, while the dashed lines are $\pm 0.1$ away from the solid lines. Note that the coordinate ranges are different in each panel. For $\lambda_{z}$, the comparison of cold components is quite good, with an average bias $\approx 0$. The warm fractions are $0.07$ overestimated and the counter-rotating fractions are $0.05$ overestimated on average, while the hot fractions have an average underestimation equal to $0.12$. Ignoring galaxy morphologies, these results are consistent with Zhu et al. (2018b), whose sample mainly consist of spiral galaxies, together with some oblate and triaxial early-type galaxies.

Comparison of the luminosity fractions of long-axis tubes, box orbits and slowly-rotating orbits is separately shown in Figure 11 (symbols are the same as in Figure 10). For oblates and triaxials with face-on views, modelling tends to overestimate both prograde and counter-rotating long-axis tubes, while correspondingly box orbits are underestimated. For prolate galaxies, the long-axis tube fractions for most viewing angles are lower than the true values, while the box orbit fraction is higher. These results indicate a clear modelling degeneracy between long-axis tube orbits and box orbits. Slowly rotating orbit fractions are underestimated for almost all mock data sets, and is the main reason behind the bias of hot orbit fractions.

6.4 Velocity anisotropy and tangential fraction

We calculate the luminosity-weighted velocity anisotropy $\beta_{r}(r)$ along spherical radius $r$ in the same way from the simulations and from our models. The three-dimensional space is divided into different cells. For the cell $i$ located at radius $r_{i}$, we can obtain the stellar mass $M_{*}(r_{i})$ and the velocity anisotropy $\beta_{r}(r_{i})$. Thus we have $\beta_{r}(r)=\Sigma_{i}[M_{*}(r_{i})\times\beta_{r}(r_{i})]/M_{*}(r)$

In Figure 12, we compare $\beta_{r}(r)$ from the simulated galaxies with $\beta_{r}(r)$ from the models. The green lines represent the true velocity anisotropy profiles and blue dashed lines with different markers are model values obtained from the mock data at different viewing angles. The red dots with error bars show the mean values of results from mock data with the five different projections and the average confidence levels within $1\sigma$. Due to the limited resolution of the Illustris simulation, the green lines do not extend to the most inner regions ($\sim 0.5$ kpc) of some galaxies. The oblate galaxies tend to be tangential at $\sim 10$ kpc, indicating that they are dominated by short-axis tube orbits in these regions. Triaxials and prolates become more and more radial with increasing radius. Our models roughly match the true $\beta_{r}(r)$ profiles for our different types of galaxies, although not without some bias.

In order to analyse the model biases, we plot the tangential fraction $f_{t}$ (see definition in Equation 26) in Figure 13. We change the linear radius to log radius in this figure. The symbols are the same as Figure 12. The model profiles track the true profiles, but with some underestimation for most galaxies, especially in the outer parts. This means that our models are generally more radial than our simulated galaxies. Long & Mao (2012) also find a similar bias towards radial anisotropy, using a different dynamical method Made-to-Measure (M2M).

We calculate the differences in tangential fraction $\Delta f_{t}(r)=f_{t}(r)_{\rm model}-f_{t}(r)_{\rm true}$ between the model and true profiles and show them in Figure 14. By performing a simple straight line fit for $\Delta f_{t}(r)$, we find a clear systematic trend. The tangential fraction does not have obvious bias in the inner regions. However, with increasing radius, the underestimation becomes larger. The average bias $\Delta f_{t}(r)$ at $\overline{R_{\rm e}}$ is around 0.064 and reaches 0.1 at large radii.

6.5 Testing the effect of initial conditions

How to select orbits is always important in Schwarzschild modelling. We do not attempt to vary the method by which the orbit initial conditions are created, but examine varying the number of orbits and assessing the impact.

Our results in $\S$ 6.1 to $\S$ 6.4 are obtained from empirical settings of the initial conditions: $n_{E}\times n_{I_{2}}\times n_{I_{3}}=21\times 10\times 7$. In order to check whether some systematic biases are caused by the settings, we performed some tests on the effect of orbit sampling. We select three observations: galaxy O1 with viewing angles $(\theta,\varphi)=(42,60)^{\circ}$, galaxy T2 with viewing angles $(\theta,\varphi)=(69,6)^{\circ}$ and galaxy P1 with viewing angles $(\theta,\varphi)=(84,28)^{\circ}$. We change the values of $n_{E}\times n_{I_{2}}\times n_{I_{3}}$ to be $15\times 7\times 5$ (fewer input orbits) and $31\times 15\times 10$ (more input orbits), and then rerun the models. We check the estimates of mass distributions, morphologies and orbit distributions for these new models, and do not find any significant differences. The only noticeable difference is that the $1\sigma$ confidence level given by the modelling becomes larger with increasing numbers of input orbits. This could be caused by the NNLS program — the numbers of active orbits (output orbits with nonzero weights) is not directly proportional to the number of input orbits. The estimates do not have significant changes, but more freedom is allowed in the modelling with more orbits sampled. The detailed results concerning active orbits are shown in Table 2.

7 Discussion

From the results shown in $\S$ 6.1, Schwarzschild modelling and JAM modelling have similar abilities in being able to recover total mass. Schwarzschild modelling has more freedom, and usually gives a number of reasonable models around the best-fitting model, while JAM modelling tends to underestimate errors due to its strong assumptions. For stellar mass and dark matter mass, both methods give results with large biases. The computer time required in using each scheme is quite different. As Schwarzschild modelling is an orbit-superposition method, it is time-consuming to construct the orbits particularly when dithering is used. For each combination of free parameters, we need to integrate $4410\times 5^{3}\approx 5\times 10^{5}$ orbits. A typical time scale for arriving at a best-fitting triaxial Schwarzschild model for our mock galaxies is $2\sim 3$ days on a 256 cores cluster, while JAM is around 100 times faster. For investigators who only want to determine the masses of galaxies, using JAM is more efficient, and little mass information will be lost by comparison with using Schwarzschild’s method. The advantage of Schwarzschild’s method is orbit information that JAM can not provide.

Since only the total mass is constrained by the kinematic data, it is difficult to separate the stellar and dark matter components directly from kinematics. The only way to distinguish them is through the inner mass density slope $\gamma$. As our modelling takes the projected surface brightness as input, we are able to give a good estimate of the stellar density slope, but the dark matter density slope can only come from hypothesis. We have assumed that the dark matter mass follows the NFW profile, so the inner dark matter density slope equals $-1$, while the stellar density slope is usually steeper than $-1$. However, the dark matter density slope from simulations is close to 0 (see Figure 5). This means that the modelling must control the fraction of stellar mass to ensure that the total mass fits the true value. This directly relates to why the stellar mass has a systematic underestimation. Clearly, a good assumption of the dark matter profile is crucial for estimating these two mass components more accurately. By changing the dark matter profile from NFW to gNFW (generalized NFW), we find the biases in stellar mass and dark matter reduce appreciably. However, there is still a large uncertainty for the inner density slope of dark matter and it is not recovered well. Due to computer time limitations, we only model a few galaxies with the gNFW profile as it has more free parameters. We take galaxy O1 with viewing angles $(42,60)^{\circ}$ as an example to show the difference between density profiles in our modelling by using NFW and gNFW haloes in Figure 15.

In $\S$ 6.2, we find intrinsic shapes have large uncertainties, especially for oblate and triaxial galaxies with face-on viewing angles. We suggest there are two reasons for this. The main reason is that the MGE method can only restrict the minimum value of the axis ratio $q_{\rm min}$ from the projected surface brightness. For face-on views, it is very possible for modelling to overestimate $q$. The other axis ratio $p$ is related to $q$ $(p>q)$, so that model estimates are more nearly spherical. The second reason is that, in our models, dark matter is taken to be spherical while the stellar component is triaxial. This means that the degeneracy between stellar mass and dark matter mass may influence the intrinsic shapes. Overall, we believe that our existing information is not sufficient to constrain the intrinsic shapes very well.

By combining the estimation of intrinsic shapes and circularity $\lambda_{x}$ distributions, we find the bias of axis ratios have a clear relationship with the degeneracy between long-axis tubes and box orbits: both the overestimation of long-axis tube fractions for face-on views in oblates and triaxials and the underestimation of long-axis tube fractions in prolates requires the intrinsic shapes to become more triaxial, which is consistent with our results in Figure 7.

In $\S$ 6.3, there is a systematic bias identified in the hot orbit fractions. We separate hot orbits into four parts and find the bias is mainly caused the slowly-rotating orbits. These orbits act like “bridges” that connect box orbits and tube orbits. In galaxy formation, the orbit distribution is usually continuous so the circularity ranges of box orbits and tube orbits could be broader. However, in our modelling, we sample many box orbits on equipotential planes, which means these box orbits will concentrate on the peaks at $\lambda_{z}=0$ and $\lambda_{x}=0$. Our models lack these slowly rotating orbits and in that respect our orbit initial conditions are deficient.

We analyse the orbit distributions in $\S$ 6.3 and $\S$ 6.4. From the anisotropy parameters $\beta_{r}$ and $f_{t}$, we find a systematic trend that the model values are more radial than our simulated galaxies. This trend does not appear in circularity distributions. Although both velocity anisotropy and circularity are orbit properties, what they tell us is different. A key point is that circularity could not distinguish different kinds of box orbits, which are all located on the central peak $\lambda_{z}\sim 0$ and $\lambda_{x}\sim 0$. Vietri & Schwarzschild (1983) studied box orbits in triaxial galaxies and introduced a width parameter to classify box orbits. From their paper, some box orbits are “wide”, like peanuts; while some are “narrow”, like pens. These “narrow” box orbits could be quite radial, while the “wide” box orbits may not be — they can also have large fractions of tangential velocity components. A possible explanation of our results is that modelling overestimates radial box orbits but underestimate other box orbits. A detailed investigation of the properties of these different orbit families may help us to understand this problem. Since our results are almost independent of sampling orbits (see $\S$ 6.5), the way to reduce the biases may be by finding some useful restriction to the orbit families.

In assessing our work, we are mindful that any results where we are comparing model 6-dimensional data (three positions and three velocities) with our simulated galaxy data are subject to deprojection effects because of the 3-dimensional data (two positions and one velocity) we use as constraints. This means that orbit circularity and velocity dispersion model results, for example, at best can only be illustrative of what the true galaxy might be like, and can not be regarded as being accurate. This deprojection issue is compounded by having to make assumptions as to what orbit initial conditions will be required to model any given galaxy. As we have seen above, these initial conditions do influence the results we are able to achieve.

8 Summary

Our objectives were to understand how well Schwarzschild’s method is able to estimate the underlying properties of our test galaxies. We considered four properties in particular of our simulated galaxies,

the mass profile including both stellar and dark matter, 2. 2.

galaxy morphology, 3. 3.

the orbit circularity distribution of each galaxy, and 4. 4.

its velocity dispersion anisotropy.

We have successfully met our objectives and are now in a position where we can make a recovery or estimation assessment for each property.

• •

For the total mass within one $\overline{R_{\rm e}}$, the estimate is quite good for most galaxies (eight out of nine) with absolute average relative deviations within 15 percent. The one exception is a prolate galaxy where the deviation is 36 percent. Separating the mass types, the stellar mass within one $\overline{R_{\rm e}}$ is on average $\sim 24$ percent lower than the true values while the dark matter mass is $\sim 38$ percent higher. These deviations are comparable with those elsewhere in the literature, notably Thomas et al. (2007). Using a reduced galaxy sample (nine mock data sets) with a gNFW profile, these values improve to $\sim 13$ percent underestimation for stellar mass and $\sim 18$ percent overestimation for dark matter. Perhaps the anisotropy increase noted below is partially responsible.

• •

For galaxy morphology, both the intrinsic shape parameters ($p$ and $q$) tend to be overestimated, with, with $\Delta p=0.07$ and $\Delta q=0.14$. The triaxiality parameter is generally recovered well but with model oblate and prolate galaxies tending to be slightly more triaxial than the actual test galaxies.

• •

Estimates of the $\lambda_{z}$ circularity distribution are plausible, with the fractions of warm and counter-rotating categories being 0.05 and 0.07 overestimated, while the hot category is under-estimated by 0.12. The estimate of $\lambda_{x}$ is also plausible for three categories except the slowly-rotating orbits, which is on average 0.13 underestimated and dominates the bias in the hot category.

• •

Comparisons of the velocity anisotropy and tangential fraction profiles with their galaxy profiles show that the galaxy models tend to be more radial in their outer regions, with the average bias in tangential fraction $f_{t}$ being $\sim 0.064$ at one $\overline{R_{\rm e}}$ and $\sim 0.1$ at higher radii. The match in the inner regions is quite good.

For the future, we are in the process of applying triaxial Schwarzschild modelling to observations of real galaxies such as in the MaNGA survey. We expect to take the findings of this paper into account in interpreting the results we achieve.

Acknowledgements

We thank R. C. E. van den Bosch for providing us his triaxial Schwarzschild software and M. Cappellari for making his MGE software publicly available. We are also grateful to the Illustris team for the their simulated galaxy data. The modelling was accomplished on the “Zen” cluster at National Astronomical Observatories, Chinese Academy of Sciences (NAOC) and on “Venus” at Tsinghua University. This work is partly supported by the National Key Basic Research and Development Program of China (No. 2018YFA0404501 to SM), by the National Science Foundation of China (Grant No. 11821303, 11333003, 11390372 and 11761131004 to SM). LZ acknowledges support from Shanghai Astronomical Observatory, Chinese Academy of Sciences under grant NO.Y895201009. GvdV acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 724857 (Consolidator Grant ArcheoDyn).