{'fraction_non_alphanumeric': 0.03308705124561902, 'fraction_numerical': 0.014900066306715923, 'mean_word_length': 4.491053781337771, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
SMART: A new implementation of Schwarzschild's Orbit Superposition technique for triaxial galaxies and its application to an N-body merger simulation

We present SMART, a new 3D implementation of the Schwarzschild Method and its application to a triaxial N-body merger simulation. SMART fits full line-of-sight velocity distributions (LOSVDs) to determine the viewing angles, black hole, stellar and dark matter (DM) masses and the stellar orbit distribution of galaxies. Our model uses a 5D orbital starting space to ensure a representative set of stellar trajectories adaptable to the integrals-of-motion space and it is designed to deal with non-parametric stellar and DM densities. SMART's efficiency is demonstrated by application to a realistic N-body merger simulation including supermassive black holes which we model from five different projections. When providing the true viewing angles, 3D stellar luminosity profile and normalized DM halo, we can (i) reproduce the intrinsic velocity moments and anisotropy profile with a precision of ~1% and (ii) recover the black hole mass, stellar mass-to-light ratio and DM normalization to better than a few percent accuracy. This precision is smaller than the currently discussed differences between initial-stellar-mass functions and scatter in black hole scaling relations. Further tests with toy models suggest that the recovery of the anisotropy in triaxial galaxies is almost unique when the potential is known and full LOSVDs are fitted. We show that orbit models even allow the reconstruction of full intrinsic velocity distributions, which contain more information than the classical anisotropy parameter. Surprisingly, the orbit library for the analysed N-body simulation's gravitational potential contains orbits with net rotation around the intermediate axis that is stable over some Gyrs.


INTRODUCTION
Early-type galaxies (ETGs) have long been believed to emerge from collisions between other smaller progenitor galaxies (first proposed by Toomre & Toomre 1972) but nowadays it is clear that their formation history is more complex (e.g. Oser et al. 2010). Their structural and kinematic properties divide them into (1) fainter (absolute magnitude M B > −20.5) and coreless fast rotators which are nearly axisymmetric and have discy-distorted isophotes and (2) brighter and more massive slow rotators with flat cores, which are moderately triaxial and have boxy-distorted isophotes (Faber et al. E-mail:bneu@mpe.mpg.de 1987;Bender 1988a;Bender et al. 1989;Kormendy & Bender 1996;Cappellari et al. 2007;Emsellem et al. 2007). For the formation of fainter elliptical galaxies dissipational processes are believed to be important (e.g. Bender et al. 1992; Barnes & Hernquist 1996;Genzel et al. 2001;Tacconi et al. 2005;Cappellari et al. 2007;Hopkins et al. 2008;Johansson et al. 2009), whereas the latest evolutionary phases in the formation of massive ellipticals are dominated by collisionless processes (e.g. Naab et al. 2006;Cappellari 2016;Naab & Ostriker 2017;Moster et al. 2019). In general, these merging events modify the potential structure and populate a rich diversity of stellar orbits (Röttgers et al. 2014). The intrinsic shape and orbital structure of such galaxies are not directly observable. Instead, sophisticated dynamical models are needed to process kinematic and photometric observational data to extract all the information about the orbital structure and internal composition of the galaxy. Dynamical models are based on the collisionless Boltzmann equation which governs the motion of stars in elliptical galaxies. Dynamical models that go beyond the recovery of velocity moments and aim at reconstructing the entire galaxy structure, additionally take advantage of the Jeans theorem (e.g. Binney & Tremaine 2008). This implies that the distribution function, which is the most general description of a system of stars, is constant along individual trajectories in phase-space. In this regard, Schwarzschild (1979) pioneered an orbit superposition technique, where the equations of motion are numerically integrated for a finite number of stellar trajectories embedded in an assumed gravitational potential with contributions from the stars and possibly dark components. The weighted superposition of the orbits is determined for which the surface brightness and projected velocity distributions of the model match the observed ones in a least squares sense (e.g. Richstone & Tremaine 1984). Besides the orbital weights, all unknown quantities like the central black hole mass and dark matter distribution are varied between different models. The model producing the best fit to the projected velocity distributions is then associated with the correct model parameters. Any galaxy in a steady state can be modeled by SchwarzschildâĂŹs orbit superposition technique. In order to determine both, the mass and internal motions of the stars, and solve an underlying mass-anisotropy-entanglement one needs to describe the deviation of the observed absorption lines from a Gaussian profile by additional Gauss-Hermite functions of at least third and fourth order (Gerhard 1993;van der Marel & Franx 1993;Bender et al. 1994), or, preferably if the signal-to-noise ration permits, measure the line shape non-parametrically (e.g. Mehrgan et al. 2019). Early applications of Schwarzschild's orbit superposition technique concentrated on spherical models (e.g. Richstone & Tremaine 1985;Rix et al. 1997). Since the simplified assumption of spherical symmetry is not true for most galaxies, later applications of Schwarzschild's orbit superposition technique assumed axisymmetry (e.g. van der Marel et al. 1998;Cretton et al. 1999;Gebhardt et al. 2000;Thomas et al. 2004;Valluri et al. 2004). However, it is nowadays known that the most massive galaxies are neither spherical nor axisymmetric but triaxial objects. Observational indications are provided by isophotal twists in the surface brightness distribution, velocity anisotropy, minor axis rotation, kinematically decoupled cores and the statistical distribution of the ellipticity of the isophotes (Illingworth 1977;Bertola & Galletta 1978, 1979Williams & Schwarzschild 1979;Schechter & Gunn 1978;Bender 1988b;Franx & Illingworth 1988;Vincent & Ryden 2005). Schwarzschild (1979) proved the existence of self-consistent triaxial stellar systems in dynamical equilibrium with numerical orbit superposition models. Also, N-body simulations supported the idea of triaxial ellipsoidal stellar bulges and dark matter halos (e.g. Aarseth & Binney 1978;Hohl & Zang 1979;Miller & Smith 1979;Barnes 1992;Jing & Suto 2002;Naab & Burkert 2003;Bailin & Steinmetz 2005). High mass galaxies are of particular interest in several astrophysical aspects, e.g., the stellar initial mass function (IMF) is discussed to vary among galaxies, with the largest excess stellar mass compared to the locally measured Kroupa (Kroupa 2001) or Chabrier (Chabrier 2003) IMF occuring in the most massive galaxies (e.g. van Dokkum & Conroy 2010;Cappellari et al. 2012;Vazdekis et al. 2015;Parikh et al. 2018; Thomas et al. 2011;Posacki et al. 2015;Treu et al. 2010;Smith et al. 2015). Moreover, different growth mod-els for supermassive black holes (SMBH) predict different amounts of intrinsic scatter at the high-mass end of SMBH scaling relations (e.g. Peng 2007;Jahnke & Macciò 2011;Somerville & Davé 2015;Naab & Ostriker 2017). In order to address these questions precision dynamical mass measurements of the stars and SMBHs are required and these are directly linked to triaxial modeling to avoid artificial scatter introduced by wrong symmetry assumptions. Thomas et al. (2007) showed that the stellar mass-to-light-ratio (and, thus, indirect inferences about the IMF) can be biased by up to 50% in extreme cases when using axisymmetric models for a maximally triaxial galaxy. Moreover, a wrongly assumed mass-to-light ratio influences the determination of the mass of the central black hole in the model. The work by van den Bosch & de Zeeuw (2010) suggests that the assumption of axisymmetry may bias black-hole measurements in massive ellipticals. They find that the best-fitting black hole mass estimate doubles when modeling NGC 3379 with their triaxial code (van den Bosch et al. 2008) in comparison to axisymmetric models. Triaxial dynamical modeling routines are therefore required to recover unbiased stellar mass-to-light ratios and black hole masses with the best possible accuracy. To understand the uncertainties and ambiguities of triaxial modeling one has to understand the following three essentially different effects: 1) the intrinsic uncertainty of the applied dynamical modeling algorithm, which can only be tested under circumstances where the solution is designed to be unique. This is one of the aspects covered in this paper. 2) the uncertainty in the reconstruction of orbital and mass parameters from typical observational data given the right deprojection, which is also addressed in this paper, 3) the uncertainty of the deprojection routine. This topic is covered in de Nicola et al. 2020.
All previously described effects need to be combined to evaluate the uncertainties in the whole modeling process. This will be investigated in a future paper.
So far, there exist two dynamical modeling codes using Schwarzschild's orbit superposition technique dealing with triaxiality by van den Bosch et al. (2008) and Vasiliev & Valluri (2019). Their estimated precision and efficiency will be later mentioned in the discussion. In this paper, we present our newly developed triaxial Schwarzschild code called SMART ("Structure and MAss Recovery of Triaxial galaxies") and test the code on a realistic high-resolution numerical merger simulation including supermassive black holes by Rantala et al. (2018). This paper is structured as follows: We introduce SMART and describe its specific benefits in Section 2. We then discuss the most important aspects for choosing this particular simulation in Section 3. In this section we furthermore explain all relevant steps to extract the data needed for modeling the simulation. In Section 4 we will show the results of these models. Section 5 deals with the quasi-uniqueness of the anisotropy recovery when fitting full lineof-sight velocity distributions. This is followed by a short discussion about remaining sources of systematics and comparison to other triaxial modeling codes in Section 7. We summarize our results and conclusions in Section 8. but it is unique in calculating the potential by expansion into spherical harmonics, setting up an adaptive orbit library and computing the orbit superposition by maximizing an entropy-like quantity (green panels). This results in specific advantages (blue panels), e.g. that our code is able to deal with realistic changes in the gravitational potential, e.g. when the SMBH causes a more spherical potential in the center. Our orbit library is adaptive and responds to changes in the integrals of motion. SMART can process density output from any deprojection routine, e.g. a non-parametric deprojection dealing with degeneracies.

TRIAXIAL SCHWARZSCHILD CODE SMART
SMART ("Structure and MAss Recovery of Triaxial galaxies") is a fully 3D orbit superposition code based on the axisymmetric code of Thomas et al. (2004) and its original extension to three dimensions and non-axisymmetric densities by Finozzi (2018). It is written in FORTRAN 90/95 (Brainerd et al. 1996). SMART follows the classical Schwarzschild method consisting of the computation of the potential and forces for a given density, the setting up of an orbit library and the subsequent superposition of the orbits. In the following sections we will explain in more detail how SMART creates self-consistent density-, potential-and orbit-configurations and how we weight the orbits in order to fit the input density and velocity structure. One main feature of our code is that it uses a five dimensional starting space for the orbit library to adapt to potentials with a radially varying structure of the integrals-of-motion space. Fig. 1 gives a schematic overview of the code's main modules and the benefits of their specific implementation.

Coordinate systems and binning
We use two different coordinate systems to describe the intrinsic and projected properties of a galaxy. To transform between the intrinsic coordinates (x, y, z), adapted to the symmetry of the object, and the coordinates (x , y , z ) adapted to the sky projection, two matrices P and R are used: and P = − sin ϕ cos ϕ 0 − cos ϑ cos ϕ − cos ϑ sin ϕ sin ϑ sin ϑ cos ϕ sin ϑ sin ϕ cos ϑ . (3) P and its corresponding viewing angles ϑ and ϕ project to the plane of the sky with z being the line of sight. R and its corresponding rotation angle ψ rotate the coordinates x and y in the plane of the sky along z . If not stated otherwise, the intrinsic long axis is hereafter assumed to coincide with x, the intrinsic intermediate axis with y and the intrinsic short axis with z.
SMART works with a cell structure based on spherical coordinates. Intrinsic properties like the stellar or dark matter distribution or individual orbital properties are integrated over small cells in configuration and/or velocity space. We use a linear sampling for the longitude θ ∈ [−90 • , 90 • ] and the azimuth φ ∈ [0 • , 360 • ]. Radial bins are spaced in even intervals of the radial binning index The constant c allows to adapt the central binning scheme from logarithmic (c = 0) to linear (c = 1). The constants a and b are determined once the radial extent of the library r min r r max and the number of radial bins N r are set and they are chosen so that the minimum radius r min lies within the first radial bin and the maximum radius r max in the last one (see also Siopis et al. 2009). Similar to the spatial properties, the line-of-sight velocity distributions (LOSVDs) are integrated over small cells in phase space given by the spatial pattern of the observations (e.g. Voronoi bins; Cappellari & Copin 2003) and the velocity resolution of the LOSVD data. Like its axisymmetric predecessor (Thomas et al. 2004), the code uses the entire information contained in the full LOSVDs. See Sections 3.1.3 and 3.1.4 for more details.

Density and potential
The total gravitational potential is composed as the sum of a Keplerian contribution from a supermassive black hole (Φ SMBH ) and the contributions from the stars (Φ * ) and dark matter (Φ DM ). SMART allows to use non-parametric densities for both the stars and the dark matter. The stellar density is generally assumed to be provided in 3D tabulated form (as for example returned from a non-parametric deprojection, e.g. de Nicola et al. 2020). The same holds for the dark matter halo. However, the code can also run with quasi-parametric deprojections (e.g. Multi Gaussian Expansion or MGE models, Monnet et al. 1992;Emsellem et al. 1994;Cappellari 2002). It can also run with parametric dark matter halos (e.g. NFW profiles; Navarro et al. 1996). The solution to the Poisson equation is obtained with the help of an expansion in spherical harmonics (Binney & Tremaine 2008). For this, the stellar and dark matter density are individually interpolated by first performing a bi-linear interpolation among the elevation and azimuthal angle bins and afterwards a linear interpolation among the logarithm of the radial bins. For integrating these interpolated densities we use a 10-point Gaussian quadrature algorithm from Press et al. (1996). The advantage of calculating the potential by expansion into spherical harmonics in comparison to other techniques, e.g. by using the MGE method, is its ability to deal with non-parametric densities and halos.

Orbit library
Every orbit in a gravitational potential is uniquely defined by its integrals of motion (e.g. Binney & Tremaine (2008)). The number of (isolating) integrals of motion depends on the given potential. Furthermore, every integral of motion reduces the dimensionality of the trajectories of the stars in the galaxy. Regular potentials admit in general three integrals of motion, one of which being the energy E. In the axisymmetric case another classical integral is known explicitly: the z-component of the angular momentum, L z . The third integral, I 3 , is usually not known explicitly. In the general triaxial case, only the energy is given explicitly. Near the central SMBH, however, the potential becomes more and more Keplerian and the number of isolating integrals of motion in a Keplerian potential is five. An example for a system in an almost Keplerian potential that is described by a 5-integral distribution function is the asymmetric disc in the center of M31 Brown & Magorrian 2013). Since our code is not restricted to axisymmetric or triaxial symmetries, we aim for a 5D starting space for the stellar orbits, which we gain by systematically sampling E, L z , v r , r and φ. The details of the initial conditions sampling technique are described in Section C. In total, SMART sets up and integrates ∼ 50000 orbits for 100 surfaces of section (SOS) crossings, i.e. for 100 crossings of the equatorial plane in upwards direction (for a more detailed description of SOS see App. C). This number of orbits was intentionally chosen to be higher than in the axisymmetric predecessor code of Thomas et al. (2004) to address all necessary complexities given by, e.g., a radially changing structure in the integrals-of-motion space. Moreover, this amount of orbits proves to be sufficient to directly recover the (phase-space) density without any dithering of orbits in the starting space.
If the potential at hand is, e.g., axisymmetric, then all orbits conserve L z and precess around the rotation axis. In this case, the orbits will fill the φ dimension automatically. Likewise, orbits will be represented by invariant curves in the (r,v θ )-plane, due to the conservation of I 3 . Hence, when sequentially sampling the orbital launch conditions, the dimension of the submanifold containing all the orbital initial conditions that are not yet represented shrinks automatically, according to the number of integrals of motion provided by the gravitational potential under study. Since, in general, triaxial potentials have three integrals of motion, a 2D starting space at a given energy (de Zeeuw 1985;Schwarzschild 1993) would provide a sufficient orbit sampling: One could sample initial conditions from the (x,z)-plane producing mainly tube orbits and compensate this with launching additional box orbits from the equipotential surfaces (van den Bosch et al. 2008). However, it is not clear whether the distribution function of realistic triaxial galaxies requires a 5D starting space near the SMBH in the center. With our choice of a 5D starting space we guarantee that our set of orbits adapts to the actual complexity of the integrals-of-motion space. In a realistic triaxial galaxy, like in the studied simulation, it changes from a more spherical center (requiring at least four integrals of motion) into nearly prolate outskirts. Furthermore, it allows us to model systems like eccentric discs with distribution functions that obviously depend on more than three integrals of motion. In Fig. C1 we show that our implemented orbit sampling and integration routine (cf. Section 2.3.1) yields a homogeneous and dense coverage of phase space.

Orbit integration and classification
SMART integrates the orbital equations of motions d ì where i denotes the orbit index, in cartesian coordinates by means of the Cash-Karp algorithm (Cash & Karp 1990). The 5th order Runge-Kutta method is implemented by using an adaptive integration step-size (see Press et al. 1996). The default integration time of the individual orbits corresponds to 100 SOS-crossings. At each integrated time-step the contribution of orbit i to the luminosity, internal velocity moments and projected LOSVDs is calculated as the fraction of time the orbit spends in the corresponding bins. Projected quantities are convolved with the relevant PSF (point spread function) in every time step and before binning. The PSF can either be provided as a parametrised two dimensional Gaussian or in terms of a PSF image. The convolution is performed via Monte-Carlo method by randomly perturbing the coordinates x (t), y (t) (cf. section 2.1) in dependence of the respective PSF. Modeling a galaxy with SMART does not require an exact orbit classification analysis. However, we built in an approximate classification method. For this, SMART checks the sign conservation of the angular momentum in x-, y-and z-direction for every SOS-crossing event. If the sign of L x is conserved for the whole integration time and if this is not true for L y and L z , the orbit gets classified as x-tube. The same applies to the other directions (cf. Barnes 1992). If the 100 SOS crossings do not hold an angular momentum sign conservation along any direction, the orbit gets classified as box/chaotic orbit. If the sign conservation is true for every direction or if the orbit shows no radial and azimuthal change during the integration time, the orbit is classified as spherical/Kepler orbit.

Orbit superposition
The orbital weights w i , which are decisive for the consistency between the observed and modeled luminosity as well as for the projected velocity profiles, are iteratively changed until the difference χ 2 between the observed LOSVDs L data and modeled LOSVDs L mod is minimal: Here, j describes the spatial bin index of the N losvd data cells and k describes the velocity bin index of the N vlos velocity bins. ∆L data j k is the error of the data in the specific bin. An advantage of SMART is that it uses the full information contained in the LOSVDs and not only the Gauss-Hermite parameters alone (cf., e.g., Mehrgan et al. 2019 for a discussion of the benefits of using non-parametric LOSVDs in measuring galaxy masses). The luminosity density serves as a boundary condition for the choice of the orbital weights.
The problem of solving for the weights w i is usually underdetermined because the number of orbits is much larger than the number of data points. We therefore regularize our models by maximizing an entropy-like quantitŷ where In the absence of any other constraints, the entropy maximisation yields w i ∝ ω i (cf. Section 5). Thus, the ω i are bias factors for the orbital weights w i and can be used to smooth the orbit model. Moreover, they can be used to construct orbit models with specific properties, e.g. orbit models dominated by certain families of orbits (cf. Section 5 for examples).
The particular form of the entropy in Equation 8 guarantees the positivity of the orbital weights. The maximum-entropy technique is flexible, however. Other choices for the entropy allow for negative weights as well (Richstone & Tremaine (1988)).
Technically, for each regularization value α, SMART maximiseŝ S by computing the relevant Lagrange multipliers. The iterative adjustment of the w i 's is performed by using Newton's method. The implemented method is based on Richstone & Tremaine (1988). A detailed description of the algorithm as well as tests demonstrating the high accuracy performance of SMART can be found in Appendix D.
As we will describe in full detail in Section 5, the entropy term in equation 7 makes the solution of the orbital weights unique. While this might be advantageous from an algorithmic point of view it comes in principle with the danger of a potential bias. As we will show, by varying the orbital bias factors ω i the maximum-entropy technique allows in principle to reconstruct any of the potentially degenerate solutions for the orbital weights (cf. Section 5). However, the results of the following Sections imply that when fitting all the information contained in the full LOSVDs, the remaining degeneracies in the weight reconstruction have only little impact on the 'macroscopic' galaxy parameters of interest, like the anisotropy for example.
One natural choice for the orbital bias factors is ω i = V i , where V i is the phase-space volume represented by orbit i (cf. e.g. Richstone & Tremaine 1988;Thomas et al. 2004). With this choice equals the Boltzmann entropy. Since the Boltzmann entropy increases during dissipationless evolutionary processes due to phase mixing and violent relaxation, galaxy models with large Boltzmann entropy are more likely than those with a small one (Richstone & Tremaine 1988;Thomas et al. 2007). However, in collisionless selfgravitating systems, every entropy-like functional is assumed to increase in phase-space, such as the generalized H-function (Tremaine et al. 1986), the entropy of the ideal gas (White & Narayan 1987) or the Tsallis entropy (Tsallis 1988). The choice of ω i is thus arbitrary to some degree. This and the fact that in case of a triaxial potential it is computationally expensive to calculate the correct phase-space volume V i for every orbit motivated us to set This functional form was also tested by de Lorenzi et al. (2007) in a slightly different context of a made-to-measure (M2M; Syer & Tremaine 1996;Bissantz et al. 2004) algorithm for N-body particle models. Compared to the Boltzmann entropy, ω i = const. leads to relative preference of orbits with actually small V i while orbits with large V i are relatively suppressed. With this choice of constant orbital bias factors ω i , the entropy equation (cf. eq. 8) resembles the Shannon entropy and yields the least "informed" set of orbital weights.

Mass Optimisation with SMART
SMART is conceived to determine the viewing angles and the mass components, like the dark matter halo, the stellar mass-to-light ratio and black hole mass by looking for the model with the smallest χ 2 . To deal with this multi-dimensional parameter space, SMART uses NOMAD (Nonlinear Optimisation by Mesh Adaptive Direct search), a software optimised for time-consuming constrained blackbox optimisations (Audet & Dennis, Jr. 2006;Le Digabel 2011). NOMAD is able to optimise a noisy function with unknown derivatives to converge to the best fitting model by using a direct-search scheme. SMART runs on multiple computer cores. The orbit processing (including the setup of the initial conditions, orbit integration, orbit classification and computation of the internal and projected velocity distributions) as well as the relevant linear algebra operations applied for their superposition are parallelised.

THE N-BODY SIMULATION
In order to test SMART on a realistic mock galaxy, we use a highresolution collisionless numerical merger simulation by Rantala et al. (2018). The simulation is a single generation binary galaxy merger of two equal-mass elliptical galaxies with an effective radius of 7 kpc hosting a supermassive black hole of 8.5×10 9 M each and corresponds to the so-called γ-1.5-BH-6 simulation in Rantala et al. (2018). The two initial galaxies are set up by using a spherically symmetric Dehnen density-potential (Dehnen 1993) with an initial inner stellar density slope of ρ ∝ r −3/2 for both progenitor galaxies. The merger results in a remnant triaxial galaxy with a supermassive black hole of 1.7 × 10 10 M with a sphere of influence 1 of r SOI ∼ 1 kpc and an effective radius of r e ∼ 14 kpc. With that, the remnant resembles NGC 1600, a galaxy showing a very large core with a tangentially biased central stellar orbit distribution (Thomas et al. 2016). The simulation is based on the hybrid treeâĂŞN-body code KETJU (Rantala et al. 2017;Karl et al. 2015) which is able to accurately compute the dynamics close to the black hole due to the algorithmic chain regularization method AR-CHAIN (Mikkola & Merritt 2006. The computation of the global galactic dynamics is based on the tree code GADGET-3 (Springel 2005). We analyse a snapshot of the simulation ∼1.4 Gyr after the galaxy centers have merged, such that the remnant can be assumed to be in a steady state. At this stage, the actual distance of the two merging black holes in the simulation is 5 pc. The merger remnant shows a radially varying triaxiality parameter, being more oblate at smaller radii and increasingly more round in the centre with a maximum of triaxiality, i.e., T = 0.5, at about 3 kpc and more prolate outskirts. The simulation contains 8.3 × 10 6 stellar particles with masses of 10 5 M each, leading to a total stellar mass of 8.3 × 10 11 M . The mass ratio of one supermassive black hole and one stellar particle is M BH /M * = 8.5 × 10 4 and therefore sufficiently large to investigate a realistic interaction of the SMBH binary with the stars (Mikkola & Valtonen 1992). The number of dark matter particles is 2 × 10 7 with masses of 7.5 × 10 6 M each, leading to a total dark matter mass of 1.5 × 10 14 M . The simulation is particularly suitable to test SMART because of its (1) very high resolution (including properly resolved black-hole dynamics); (2) realistic orbital structure and shape; (3) realistic mass composition (black hole, stars, and dark matter) with a realistic stellar density core (e.g. Thomas et al. 2009;Rantala et al. 2019). Finally, the number of stellar particles is large enough to measure fully resolved LOSVDs from the central sphere of influence of the black holes out to dark matter dominated regions.

Orientation of the Simulation
We aim at orienting our intrinsic coordinate system as closely as possible to the intrinsic symmetry axes of the merger remnant. However, the stellar and dark matter principal axes of the simulated remnant are not aligned. Hence, the orientation of the main axes depends on the radius and on the mass component for which the reduced inertia tensor (see e.g. Bailin & Steinmetz 2005) is calculated. Such a shift between the stellar and dark matter halo axes is not unexpected for collisionless merger simulations (see e.g. Novak et al. 2006). We decided to center the remnant on the stars and black holes and afterwards orientate it by using the reduced inertia tensor for stars and dark matter within 30 kpc. With this, the stellar elliptical isophotes for the three different projections are well aligned with the projected principal axes within the field of view of 15 kpc x 15 kpc. There is a negligible residual misalignment which is strongest for the major axis projection but nowhere larger than ∼ 5 • (see Fig. 3 and A1).

Density
Due to the good alignment and taking advantage of the nearly triaxial intrinsic symmetry of the merger remnant, we increase the resolution of the simulation for computing the density by a factor of 8 by folding all stellar and dark matter particles into one octant. The stars are then binned into concentric radial shells with 1000 stars in each shell and the dark matter particles are binned into radial shells with 5000 dark matter particles each. The single shells are subdivided into angular bins, such that the elevation angle θ ∈ [0 • , 90 • ] increases in constant sin(θ)-steps of 0.1 and that the azimuthal angle φ ∈ [0 • , 90 • ] increases in constant φ-steps of 10 • . Within r < 0.28 kpc for the stellar and r < 10 kpc for the dark matter particles, the resolution is too low to extract a smooth density. Here, we extrapolate the logarithmic densities from the outer parts by a first order polynomial fit in the logarithm of the radius. For the stellar densities we use the slope and vertical intercept averaged over all angular bins. The dark matter density is extrapolated by using the slopes and vertical intercept values of the individual angular bins. We ensure that the total enclosed mass is well covered by the extrapolated density. To smooth the radial density profiles, SciPy's gaussian_filter()-function (Virtanen et al. 2019) is used.

Spatial binning
For spatially binning the kinematic input data we use the Voronoi tessellation method of Cappellari & Copin (2003). For each tested projection (see Section 3.2) we construct a separate set of Voronoi bins. To end up with a roughly constant number of stellar simulation particles N * in each bin, we define the signal-to-noise ratio as When calculating the average over the five different projections, this results in a total number of N losvd = N voronoi = 227 Voronoi bins within the whole field of view and 54 Voronoi bins within r SOI . This resolution is chosen to conform with realistic observational data and does not exceed high-resolution wide-field spectral observations by, e.g., MUSE (cf. Mehrgan et al. 2019) but still proves to be sufficiently high for this analysis.

Kinematic data and velocity binning
For each spatial bin we calculate the line-of-sight velocity distributions for N vlos = 45 equally-sized velocity bins with v max min = ±1600 km s −1 , which is chosen so that it covers about 10 times the velocity dispersion. This results in a velocity resolution of ∆v vlos = 71.11 km s −1 . It is not trivial which "error" for the kinematic input data of the simulation should be used since we do not have an error in the simulation in the sense that repeated measurements give the same results. However, the χ 2 -minimisation formally requires information about an "error". We tested several assumed "error-bars", such as -the difference between two kinematic datasets, each determined by using half of the simulation particles, -the Poisson noise and -a constant absolute error for each LOSVD as 10 percent of the maximum per LOSVD. In order to prevent an underestimation of the relative error for the major axis projection holding more particles along the line of sight than the other projections, the constant absolute error proves to be the most suitable method and is used in the present analysis. The choice of setting the error value to 10% of the maximum value of each LOSVD corresponds to a velocity uncertainty of ∆v = 13 km s −1 , dispersion error of ∆σ = 13 km s −1 (or 5% relative error) and ∆h n = 0.03 for the higher-order Gauss-Hermite moments. This is a reasonable choice both in terms of real observational errors and of the scatter in the kinematic maps of the merger simulation (see Fig. 3 and B1).

Applying SMART to the simulation
For the purpose of testing our code, SMART is provided with the correct viewing angles and with the 3D normalized stellar (i.e. luminosity) density ρ * from the simulation as well as with the normalized 3D dark matter density ρ DM . The DM scaling parameter s DM is to be determined by SMART. We skipped any surface brightness deprojection, since degeneracies in the deprojection (see e.g. de Nicola et al. 2020) would only hamper a correct evaluation of our code. We parameterise the density as where our fit parameters are the black hole mass M BH , stellar mass-to-light ratio Υ and the multiplication factor s DM defining the magnitude of the dark matter density profile favored by SMART. We determine these parameters by finding the minimum in χ 2 . We model and analyse five different projections with (1) ϑ = 90 • , ϕ = 0 • , i.e., the major axis projection, (2) ϑ = 90 • , ϕ = 90 • , i.e., the intermediate axis projection, (3) ϑ = 0 • , ϕ = 90 • , i.e. the minor axis projection, (4) ϑ = 90 • , ϕ = 10 • , i.e. a projection 10 • off the major axis in azimuthal direction and (5) ϑ = 90 • , ϕ = 45 • , i.e. a projection in between the major and intermediate axis projection. Without loss of generality, ψ was set to 90 • for all three viewing directions, i.e., R (see eq. 2) equals the unit matrix. The field of view is chosen to be 15 kpc x 15 kpc. The minimum sampled starting radius is set to r min = 0.05 kpc and the maximum sampled starting radius is set to r max = 80 kpc. We find optimal results for a central binning with c = 0.5 in Equation 4 because this guarantees that in case of the simulated merger remnant the difference of the circular velocity within one radial bin equals the model's velocity resolution ∆v vlos = 71.11 km s −1 at a radius of r = 0.16 kpc = 0.16 · r SOI . For c = 1 this would be only reached at r = 0.38 kpc = 0.38 · r SOI resulting in a deteriorated black hole mass recovery by ∼ 10%.
For each tested projection we model two halfs of the LOS kinematic data: After correct projection onto the plane of the sky (see Section 2.1), we separately model the half of the kinematic data with positive x -coordinates (hereafter called 'right half' of the galaxy) and the half of the kinematic data with negative x -coordinates (hereafter called 'left half' of the galaxy).

Choice of Regularization
When applying the code to realistic noisy measurements, regularization becomes important to prevent the orbital weights to fit the noise in the data. The optimal regularization parameter α for a specific observational data set can be determined by running Monte-Carlo simulations on kinematic mock data (Thomas et al. 2005) and is given as the one providing the minimum deviation of the intrinsic properties (like the distribution function, or velocity moments, or mass parameters) in comparison to the default model. When fitting noiseless ideal data, one would expect best results for χ 2 → 0, or α → ∞ (neglecting recovery degeneracies and assuming an "error" can be defined). Even in that case, however, due to residual systematics (like finite resolution of the orbit library etc.) and due to the intrinsic noise in the N-body simulation, we still expect that the best result may not necessarily be achieved asymptotically for very large α, but already for some finite value of the regularization parameter. To take this into account, we split the code test into two phases: (1) We fit the orbit model with the correct mass parameters and determine that value of α for which the internal structure of the simulation is best recovered. Specifically, we use the 2nd order velocity moments for this comparison.
(2) We then also vary the mass parameters and test how well they can be recovered. Our benchmark is the optimised α from the comparison of the moments, but we will discuss the results for all α to demonstrate their robustness. For determining the deviation between the model's velocity dispersions σ r , σ θ and σ φ and the real ones from the simulation we define where the index i denotes the three coordinates r, θ, φ.
These rms σ -profiles of the five tested projections and their respective modeled halfs in dependence of the regularization parameter α are plotted as thick lines in the top row of Fig. 2. The x-axis ticks thereby mark all α-values which were tested in this analysis. The second row in Fig. 2 shows the quality of the fit as χ 2 /N dataprofile (for definition of χ 2 see formula 6) again plotted against the regularization parameter. χ 2 as deviation from the kinematic input data with the modeled fit is here normalized over the number of input data N data , which is composed of the number of Voronoi bins N voronoi times the number of kinematic bins N vlos . As expected, the fit to the data is poor when α is low (high χ 2 /N data ). In this regime, it is the entropy term which is essentially maximised and the data (via α · χ 2 , cf. eq. 7) have little influence on the fit. With this, the anisotropy strongly depends on the ω i and, in our case, happens to be a poor representation of the internal moments of the merger (high rms values, see Fig. 2, first row). With increasing α, both the fit quality and the agreement with the merger structure improve. However, at very high α-values further improvements of the fit do not make the internal moments better since we are dominated by the noise of the N-body simulation. The most suitable choices of regularization can be read off from the minima min(rms σ ) of the rms σ -profiles and are marked as thick vertical lines in Fig. 2. Their average value is α(min(rms σ ))=0.41. This value, however, depends on the specific implemented setup of SMART as well as of the input data. Evaluated at the individual most suitable regularization values and afterwards averaged over all ten models (5 projections and two halfs each) we get a minimum value of only min(rms σ ) = 0.008. This extremely good agreement demonstrates that our orbit sampling represents the phase space very well. We also test the influence of rotation on the comparison. Due to the overall small angular momentum and thus small absolute velocity amplitude, however, the relative errors in v are sometimes large. The absolute error of the first order velocity moments, averaged over all angular and radial bins, is only ∆v = 3.9 km s −1 , but the maximum velocity over these bins is likewise only 23.4 km s −1 . We therefore here define rms v,σ as normalized deviation between the first internal moments and the velocity dispersions together as where the index i again denotes the three coordinates r, θ, φ. The rms v,σ -profiles are plotted against the regularization parameter as thin lines in the top row of Fig. 2. As exptected, their minima appear at regularization values similar to the ones of the rms σ -profiles. They are marked as thin vertical lines and their average value is α(min(rms v,σ ))=0.012. Again, this value depends on the specific setup. Evaluated at the individual α(min(rms v,σ ))-values and afterwards averaged over the five different projections and their respective modeled halfs we gain a value of only min(rms v,σ ) = 0.012. All values in the vicinity of α(min(rms σ )) and α(min(rms v,σ )) are good regularization choices. Within this regularization region SMART is able to well fit the kinematic input data for each tested projection (see Fig. 3 and B1). Averaged over the five different projections and their respective modeled halfs we receive The y-axis-ticks for the third and fourth row symbolise the concrete masses which were tested and used as input values for the models. The black dotted line marks unity which is achieved when the model correctly recovers the mass. The stellar mass-to-light ratio and black hole mass were recovered with an accuracy better than 5%. Such an intrinsic precision under similar conditions has not yet been demonstrated with other Schwarzschild modeling codes.
mean values and deviations ofv = (0.11 ± 2.09) km s −1 ,σ = (309.75 ± 2.54) km s −1 ,h 3 = 0.00 ± 0.01 andh 4 = 0.01 ± 0.01, again evaluated at α(min(rms σ )). The χ 2 -maps in Fig. 3 and B1 show that the models for each projection are able to fit the kinematic input data homogeneously well over the field of view with slightly larger deviations in the center. If not specifically annotated, the results shown in the further analysis are evaluated at α(min(rms σ )). However, the quality of the models does not strongly depend on the exact regularization value because a broader range of regularization values around these determined α-values is sufficiently appropriate and results in equally good mass parameter reproductions within the overall scatter (see also third and fourth line in Fig. 2 which will be explained in Section 4.3).

Reproduction of internal moments and orbit structure
The previously shown small rms σ -and rms v,σ -values already demonstrate the very good recovery of the internal moments by the model when providing the correct mass parameters and viewing angles. The high level of agreement between model and simulation is further illustrated in Fig. 4. The anisotropy parameter The contour lines correspond to isodensity surfaces. The fifth panel in the bottom row shows χ 2 losvd /N vlos as deviation from the kinematic input data with the modeled fit. We show the result for the model with the correct stellar mass-to-light ratio, black hole mass and dark matter scale factor as input parameters evaluated at the most suitable regularization parameter of α(min(rms σ )). The maps in the second row consist of the results of the modeled left and right half of the galaxy. Therefore, the χ 2 -maps show two colorbars for the two different halfs.
one (see Fig. 4, middle panel). The model is also able to reproduce the negative β within the core radius r b , which equals the black hole sphere of influence (Thomas et al. 2016), reflecting the tangential orbit distribution due to black hole 'core scouring': Within the sphere of influence ETGs at the high-mass end exhibit central regions which are fainter than an extrapolation of a SÃľrsic function (Sérsic 1963) as fit to the outer surface brightness profile would suggest. The commonly accepted theory for the formation of these 'cores' is a gravitational slingshot process of stars on radial orbits caused by SMBH binaries which were arised by galaxy mergers (e.g. Begelman et al. 1980;Hills & Fullerton 1980;Ebisuzaki et al. 1991;Milosavljević & Merritt 2001;Merritt 2006;Rantala et al. 2018).
The classification of the integrated orbits is also plotted in Fig. 4 (bottom panel). Orbits in the immediate vicinity of the black hole are spherical and Keplerian orbits as expected due to the SMBH. Z-tubes are the predominant orbits within 0.3 kpc < r < 3 kpc causing the more oblate shape of the galaxy in this range. For r > 3 kpc the majority of orbits are classified as x-tubes corresponding to the prolate shape in the outskirts of the simulation as described in Section 3. Our orbit classification analysis well matches the one by Frigo et al. in prep., which is done via an orbit frequency analysis of the simulation.

Mass recovery
So far, we have provided SMART with the correct mass parameters of the stellar mass-to-light ratio Υ sim , black hole mass M BH,sim and dark matter multiplication factor s DM,sim of the simulation. Even though SMART can be provided with any type of dark matter profile, we here transfer the dark matter density profile shape of the simulation and concentrate on finding the correct mass multiplication scale factor s DM of this predetermined halo shape (cf. eq. 12).
The following sections will show the one-dimensional mass reproduction results of individually determining the favored stellar mass-to-light ratio Υ or black hole mass M BH (Section 4.3.1) as well as the two-dimensional mass recovery results of simultaneously determining the favored Υ and M BH (Section 4.3.2) or Υ and s DM (Section 4.3.3). We skip any 3-dimensional mass parameter recovery since this would not provide more information than the combined 2-dimensional recoveries in the context of testing the orbit library. Fig. 2 shows the 1-dimensional mass recovery results of Υ (third row) and M BH (fourth row) for the five different projections and their respective modeled halfs. For testing the recovery of the black hole mass, we provide the model with the correct Υ simand s DM,sim -values and run nine models with different black hole masses within M BH ∈ [0.79M BH,sim , 1.21M BH,sim ] including the correct one. The tested mass grid has a smaller grid size close to M BH /M BH,sim = 1 and the exact tested values can be read off from the ordinate ticks in Fig. 2. For testing the 1-dimensional mass recovery of the stellar massto-light ratio we provide M BH,sim and s DM,sim and test the same Υ/Υ sim -values as for the black hole mass analysis. Fig. 2 shows the favored mass parameters, i.e., the mass parameters where χ 2 /N data is smallest, as a function of α. As one can see, for our fiducial choice of α (i.e. α = α(min(rms σ )) as the best recovery of the velocity dispersions), the average mass recovery performs excellently, with ∆M BH = 5% and ∆Υ = 2%. In fact, above log α −3, the results are very robust with little dependency on α. Within the overall minor scatter, all models, independent from the chosen half of the galaxy or projection, show equally good fits and reproductions of the internal moments and mass parameters. and azimuthal direction σ φ (dashed line) from the model (black lines) averaged over the five different projections accurately follow the real ones from the simulation (blue lines) out to 15 kpc as field of view. The grey shaded lines mark the deviations between the different projections. Middle panel: Also the anisotropy parameter β is well reproduced and represents the tangentially anisotropic orbit distribution (β < 0) within the core radius r b = r SOI ∼ 1 kpc as well as the radially anisotropic orbit distribution (β > 0) outside r b . Bottom panel: Radial distribution of the orbit fractions f orbitclass classified by SMART.

2-dimensional mass recovery of Υ and M BH
For simultaneously recovering Υ and M BH by SMART we sample a two-dimensional grid of input masses for 49 models per projection with Υ ∈ [0.79Υ sim , 1.21Υ sim ] and M BH ∈ [0.79M BH,sim , 1.21M BH,sim ]. We again model all five different projections but fit only the right half of the galaxy since the 1dimensional mass recovery showed no significant difference between the respective halfs of each projection. The results are plotted in Fig. 5. Evaluated at the regularization value α(min(rms σ )) the stellar mass-to-light ratio is in every case correctly recovered and the black hole mass is reproduced with an accuracy of 6% averaged over the different projections.

2-dimensional mass recovery of Υ and s DM
For recovering Υ and s DM we sample a two-dimensional grid of input masses for 49 models along the major axis projection with Υ ∈ [0.79Υ sim , 1.21Υ sim ] and s DM ∈ [0.7s DM,sim , 1.3s DM,sim ]. We here model the right half of the major axis projection. The result is shown in Fig. 6. Evaluated at α(min(rms σ )), the dark matter scale factor is slightly overestimated by 10% and the stellar mass-to-light ratio is slightly underestimated by 7% confirming the accurate mass reconstruction from the previous tests.
In conclusion, these mass recovery results demonstrate that our orbit sampling and superposition algorithms allow for a very accurate reconstruction of the mass composition and orbital structure of triaxial systems with known density shapes.

Beyond 2nd order velocity moments
So far, we have focused on the reproduction of the first and second order internal velocity moments. The previous sections have shown that the full shape of the LOSVDs contains enough information to accurately reconstruct the mass and the anisotropy structure of the orbit distribution. In Fig. 7 we show the full mass-weighted stellar velocity distributions per velocity bin against the velocity in km s −1 in radial, longitudinal and azimuthal direction when integrating over the other velocity components and respective spatial bins of the model. We find that the central azimuthal and longitudinal velocity distributions within r < r SOI have two maxima which become more pronounced closer to the center (see Fig. 7). This likely reflects the strong tangential anisotropy produced during the formation of the core and is probably linked to the negative h 4 parameter at the center of the merger remnant (see Fig. 3 and B1), which will be investigated in more detail in a separate paper. The whole internal velocity distributions contain more information about the formation process than the velocity moments alone can do. So, we extended SMART to calculate the internal velocity distributions for 31 velocity bins within the positive and negative escape velocity, i.e., v max min = ±v esc (r, θ, φ), evaluated at each radial and angular bin of the SMART-specific grid (see Section 2.1). The internal velocity distributions averaged within spherical shells reproduce the ones from the simulation sufficiently well, with a deviation of rms = 0.07 averaged over all velocity bins with v < 1000 km s −1 and radial bins within r ∈ [0.25 kpc, 15 kpc] (see Fig. 7). Including the outer wings of the internal velocity distributions with v > 1000 km s −1 , the rms increases, however, the number of simulation particles in these bins is very small. We will apply this ability of SMART to model the whole internal velocity distribution in future studies of real observational data. We checked the behavior of the rms-profile in dependence of the regularization as deviation of the whole internal velocity distributions and compared it with the rms σ -and rms v,σ -profiles. It thereby showed the same form and its minimum appeared in the same regularization region and therefore does not provide additional information when determining the most suitable regularization value.

Robustness and Uniqueness checks
In order to test the robustness of the results against modifications of our fiducial setup and in order to test the uniqueness of the results we make the following checks and considerations:

PSF convolution and noise
While Valluri et al. (2004) describe that three-integral, axisymmetric, orbit-based modeling algorithms in general show a flatbottomed χ 2 distribution, being unable to determine the black hole   Figure 6: 2-dimensional mass recovery results of Υ and s DM for the positive half of the major axis projection. s DM thereby is the mass multiplication scale factor of the predetermined halo shape of the simulation. We evaluate 49 models with different Υ-and s DM -input-masses covering a 2-dimensional grid with a step size of 7% for the stellar mass-to-light ratio and 10% for the dark matter scaling multiplication factor. The favored model is marked with a grey cross. The model slightly overestimates the dark matter halo scale factor by 10% and slightly underestimates the stellar mass-tolight ratio by 7%.
mass to better than a factor of ∼ 3.3, Magorrian (2006) demonstrates that this is only true for noiseless data. According to this, our model should not be able to precisely determine the correct black hole mass of the simulation due to the lack of an "error". However, the previous results have already shown that SMART achieves a well defined black hole mass due to a well defined minimum in the χ 2 -profile, which was probably supported by the intrinsic noise of the N−body simulation. Nevertheless, we check whether the minimum in the χ 2 -curve changes when simulating an "error" in the kinematic input data. For this, we model the positive half of the intermediate axis projection by adding Gaussian noise to the simulation chosen so that the velocity dispersion of the noisy kinematic simulation data results in an observationally realistic error of ∼ 3% (v = (0.14 ± 7.64) km s −1 ,σ = (288.86 ± 7.80) km s −1 ,h 3 = 0.00 ± 0.02,h 4 = 0.02 ± 0.02). To achieve even more realistic conditions we furthermore smooth the data by simulating a psf convolution with a FWHM of 2.43 arcsec which corresponds to 0.24 kpc, i.e., about a fourth of the sphere of influence, when assuming the galaxy to be at a distance of 20 Mpc. The so constructed velocity maps can be seen in Fig. B2. We provide SMART with the information about the used FWHM-value for the two dimensional psf convolution and test the same Υ and M BH input masses as in Section 4.3. The corresponding results for the 1-dimensional mass recoveries when modeling this modified input data are plotted as turquoise lines in Fig. 8, and do not differ decisively from the ones without psf convolution and noise. The stellar mass-tolight ratio is slightly underestimated by ∆Υ(α(min(rms σ )) = 3.5% or ∆Υ(α(min(rms v,σ )) = 7% and the black hole mass is overestimated by ∆M BH (α(min(rms σ )) = 3.5% or underestimated by ∆M BH (α(min(rms v,σ )) = 7%. The rms σ −, rms v,σ − and χ 2 /N data -profiles are of course shifted upwards but follow the same form than the ones without noise and psf convolution (cf. Fig. 2). Fig. 9 shows the 2-dimensional mass recovery result of Υ and M BH for this model. The stellar mass-to-light ratio is again correctly reproduced and the black hole mass is underestimated by only 7%. With this, SMART is able to model the noisy and psf convolved kinematics as well as the default kinematics without any "error" equally well within the overall scatter.

Changing the orbital bias factors ω i
As discussed in Section 2.4, the orbital bias factors ω i of Equation 8 can be used to control the orbital weights w i . In the absence of other constraints w i ∼ ω i . Our choice of ω i = 1 is somewhat arbitrary. In fact, it biases the w i s strongly away from the true solution. Thus, we want to test whether this affects our fits. Specifically, we test the opposite extreme: We remodel the positive half of the major axis, abbreviated below as X + , by setting ω i = w i (X + ), where w i (X + ) are the orbital weights of the best-fit model for X + . We also remodel X + by setting the bias factors ω i = w i (X − ) to the orbital weights of the best-fit model for the negative half of the major axis X − . Fig. 10 shows the results for these completely independent model fits. As expected, the rms-and χ 2 /N data -profiles start with smaller values, since at α = 0, the ω i s bias the orbital weights the strongest. As motivated above, in the specific case here, the weights are biased towards a previous fit, which explains the better initial χ 2 and rms. The internal velocity distributions are calculated for 31 velocity bins within the positive and negative escape velocity evaluated at each radial and angular bin of the SMART-specific grid. We here show the mass-weighted stellar internal velocity distributions per velocity bin averaged within spherical shells. The non-shaded lines correspond to the simulation data and the shaded lines show the modeled results. In closer vicinity of the black hole, the elevation and azimuthal distribution show two maxima corresponding to the tangentially anisotropic orbit distribution. SMART is able to reproduce the internal velocity distributions in general and is also able to follow this specific behaviour in the central bins. The rms-and χ 2 /N data -profile show the same shape and suitable regularization region than for the noiseless and non-psf-convolved case. Also, the Υ-and M BHreproduction shows no remarkable change within the overall scatter.
However, at our fiducial α range, the reproduction of the internal moments and the quality of the fit is the same as for the case with identical ω i s. As a consequence, this implies that the simplifying assumption of constant ω i does not change the modeling results significantly.

Degeneracy
When redoing the just described analysis (Subsection 4.5.2) but using the orbital bias factors of the right half of the intermediate axis projection ω i (Y + ) (brown line in Fig. 11) as initial values for remodeling X + we gain a minor improvement in the reconstruction of the internal moments since the rms σ -and rms v,σ -values are a bit  Figure 9: 2-dimensional mass recovery results of Υ and M BH for the positive half of a psf convolved noisy version of the intermediate axis projection. We again evaluate 49 models with different Υ-and M BH -input-masses covering a 2-dimensional grid with a step size of 7% around the correct parameters. The favored model is marked with a grey cross. The model finds the correct stellar mass-to-light ratio and slightly underestimates the black hole mass by 7%. Within the overall scatter this resembles the results of the models without noise and psf convolution.
smaller compared to the default model (black line), though the mass recovery shows equal results. The fact that, for the same quality of fit (i.e. same χ 2 ), the rms of the internal moments is smaller for the model with ω i assumed to be equal to orbital weights from another modeled projection direction (in this case from the intermediate axis projection Y ) suggests that these ω i contain some information in addition to the kinematics of the given projection (in this case the major axis projection X). The fact that information from another line of sight can improve the model is not surprising. In fact, it shows that some of the already very small residual rms in the internal velocity moments is due to remaining degeneracies in the recovery of the orbital weights. These degeneracies seem to be  Figure 11: Choice of regularization and 1-dimensional mass recovery results for the positive half of the major axis projection when using the default constant initial orbital bias factors ω i = const. (black line) and when using the best-fit orbital weights w i (Y + ) of the positive half of the intermediate axis as orbital bias factors, i.e., ω i = w i (Y + ) (brown line). The additional information from the second projection axis appears to reduce the degeneracy leading to a minor improvement in the internal moments reproduction as seen in the lower rms σ -and rms v,σ -values. surprisingly small. However, our results imply that as long as the deprojected light profile and normalized DM halo are known and the orbit sampling is dense, masses and anisotropies can be recovered with very high accuracy, independent of the viewing angle.

Change in orbit sampling technique
Even though our results already prove the efficient ability of our orbit library to completely reproduce all necessary orbits in a triaxial potential, we want to check the robustness of our orbit sampling method (cf. App. C). We therefore change the angular momentum sampling direction from the minor axis to the major axis. This corresponds to a change of v φ -to v θ -sampling. The basic idea is that a more homogeneous angular momentum sampling along the major axis instead of the minor axis might improve the mass recovery. Thus, we change the major and minor axis coordinates within the simulation and rerun the models for the new major-(left column in Fig. 12) and intermediate-axis-kinematics (right column in Fig. 12). With this, the angular momentum orbit sampling is proceeded along the major axis of the simulated remnant. The results are shown in Fig. 12. One can see that both orbit libraries produce the same results. In fact, the modified orbit sampling procedure reveals a slightly better mass recovery but slightly worse internal moments reproduction (higher rms σand rms v,σ -values) for equally good kinematic fits. Nevertheless, regardless of the chosen angular momentum sampling axis, our orbit sampling technique (cf. App. C) of creating initial conditions, which belong to certain energy shells and angular momentum sequences, has proven to be a highly efficient technique. It produces a general and complete set of orbits for triaxial potentials being able to deal with changing structure in the integrals-of-motion space since it manages to reproduce all relevant properties and to radially adapt itself to the more spherical center as well as the more prolate outskirts of the simulated galaxy.

Summary of Robustness and Uniqueness Checks.
In conclusion, these checks prove that SMART is robust against minor internal modifications as well as input data changes. SMART has proved its ability to handle with noisy and psf convolved data. Furthermore we have shown that constant orbital bias factors are a good approximation. Even though it is impossible to get observations from two viewing points, we have demonstrated that this would allow to reduce the minor degeneracies allowed by the kinematic data even further. We have verified this by using the orbital bias factors from a second projection direction. Moreover, the results are not affected by changing the orbit sampling technique from setting up L z -sequences to setting up L x -sequences. This shows that, as expected, the choice of sampling axis is not decisive and that the orbit sampling routine is universal.

THE QUASI-UNIQUENESS OF THE ANISOTROPY RECONSTRUCTION WHEN FITTING FULL LOSVDS
The results of the previous Section have shown that the anisotropy of the N-body merger remnant can be reconstructed with very high accuracy from the Schwarzschild models that we fitted to the full LOSVDs. This not only demonstrates the high accuracy of our orbit superposition model, but it also implies that when models can exploit the full information contained in the entire LOSVDs, the remaining degeneracy in the recovery of the distribution function can not affect the anisotropy or mass recovery significantly. In this Section we want to use the maximum entropy technique to explore this in more depth.

The Maximum Entropy Technique and the Mathematical Structure of the Solution Space
As already described in Section 2.4, we solve for the orbital weights w i by using a maximum entropy technique (cf. eqs. 7 and 8). The χ 2 term in Equation 7 contains the kinematical constraints and is the deviation between observed LOSVDs L data and the model prediction L mod , i.e. the weighted sum over the contributions of all orbits to LOSVD j and line-of-sight velocity bin k, It is convenient to think of the observed ì L data and the model predictions ì L mod as vectors with N data = N losvd × N vlos elements. Then, where L orb is a matrix with N data rows and N orbit columns. In addition to the kinematical observations ì L data , the orbital weights ì w are subject to photometric constraints. In analogy to Equation 16: where ì p mod is a vector with the model predictions for the 3D luminosity density at spatial position j 3 in the galaxy ( j 3 = 1, . . . , N phot ). To guarantee the self-consistency of our model, the respective observed ì p data are not included via a χ 2 term. Instead, we treat them as boundary conditions for the fit: Hence, we seek for the maximum of Equation 7 subject to the linear equality constraints For convenience, we normalise the ì p data such that ì p data = 1. Since we are only interested in positive orbital weights (see below), the orbital weights obey 0 w i 1 and we can restrict the maximization to the respective N orbit -dimensional convex quader. To show that Equation 7 has a unique global maximum that can be controlled through the bias factors ω i it is convenient to consider the equivalent minimisation problem for f ≡ −Ŝ given by multiplying Equation 7 with −1. Let us first consider Equation 7 without the entropy term S. The χ 2 term can be written as where C v is the covariance matrix of the observed LOSVDs and is positive definite. The Hesse matrix of χ 2 reads Because L orb is positive by construction, the symmetric matrix ∇ 2 χ 2 is at least positive semi-definite and χ 2 is convex. The minimisation of χ 2 alone, subject to the linear equality constraints Equation 19, is therefore a convex optimisation problem with affine equality constraints in standard form. As such, it only has a global minimum (e.g. Boyd & Vandenberghe 2004). In general we cannot assume that χ 2 is strictly convex (i.e. that ∇ 2 χ 2 is positive definite). The set of orbital weights that solve χ 2 ( ì w) = χ 2 min may therefore be non-unique. This is not surprising given that the linear equation will in general be underconstrained if N orbit > N data . As already mentioned above, the reconstruction of the distribution function is not unique even if the deprojected density and the potential are known. We give an example in Appendix D.
Because the Hesse matrix of −S is diagonal with the i-th element equal to 1/w i , the entropy −S is strictly convex. In contrast to the case of χ 2 , the set of orbital weights that minimise −S (or, equivalently, that maximise S) is unique. For the entropy alone this is easy to see since can be solved analytically: w i = exp(−1) · ω i . The constraints from Equation 19 will shift the solution, but the strict convexity still guarantees it to remain unique (e.g. Boyd & Vandenberghe 2004). For general α > 0, the Hesse matrix of f is the sum ∇ 2 (−S)+α∇ 2 χ 2 and f is always strictly convex. Equation 7 subject to the constraints (19) hence always has a unique solution. As already mentioned above, from an algorithmic point of view it is advantageous to have a unique solution. From the physical point of view this is not desirable because any algorithm that picks up only one of the potentially many solutions that minimise χ 2 in Equation 6 may lead to a bias. However, suppose ì s 1 and ì s 2 are two solutions which lead to the same χ 2 min . By setting ω i = s 1,i · exp(1) we can make ì s 1 the global solution of Equation 7. Likewise, by setting ω i = s 2,i · exp(1) we can make ì s 2 the global solution of Equation 7. This shows that the maximum-entropy formulation of the problem allows in principle the reconstruction of the entire solution space (via variation of the ω i ).

Testing the Uniqueness of the Anisotropy Recovery
The previous Section has shown that the recovery of the distribution function (or, equivalently, of the orbital weights) is in general not unique, even when the deprojection and the potential are known. The maximum entropy technique recovers one of the many possible solutions, which is unique for every given set of ω i . Variation of the ω i allows to sample the full solution space. On the other hand, the fits of the N-body merger remnant have shown that the anisotropy recovery is very accurate and stable even to variations of the ω i . Moreover, in Appendix D we explicitly construct two different phase-space distribution functions that fit a given set of kinematics equally well. Even though they have different orbital weights, they reveal very similar anisotropies in the second-order velocity moments. This suggests that while the recovery of the full distribution function is non-unique, the anisotropy of the secondorder moments is actually very well constrained by the information contained in the full LOSVDs. To investigate this further, we create kinematic data of several toy models with different intrinsic properties: with randomised orbit weights (called RANDOM in the following), with an overpopulation of box/chaotic orbits (called BOX), with an overpopulation of z-tubes (called ZTUBE), with an overpopulation of only the prograde z-tubes (called ZROT) and with an overpopulation of x-tubes (called XTUBE). We then fit the mock kinematic data of these toy models under different choices for the ω i , trying to push the fitted orbit model towards extreme shapes. The goal is to test how accurate and stable the recovery of the internal velocity anisotropy is, when we fit the entire information contained in the LOSVDs. The toy models are constructed as maximum-entropy models (i.e. through maximisation of eq. 7 with α = 0; cf. Thomas et al. 2007). For the BOX model, we increase the ω i of all box orbits by a factor of 1000, for the XTUBE and ZTUBE models we increase the ω i of the respective tube orbits by the same amount. For the ZROT toy model we increase the ω i only for prograde z-tubes and for the RANDOM model we use randomised ω i (cf. Appendix D). The weights are then still forced to satisfy the density constraints of the N-body simulation. For each toy model, we create kinematic mock data for four different projections (major-, intermediate-, minor-axis projection and the θ = 45 • , φ = 45 • -projection). We then model every projection of all toy models (in total 20 different input data) and test four different methods for the fits: We use (i) our default constant orbital bias factors, i.e. ω i = 1 corresponding to the Shannon-entropy (abbreviated as 'shannon' in Fig. 13), (ii) increased orbital bias factors by a factor of 10 for the box/chaotic orbits (box10), (iii) increased orbital bias factors by a factor of 10 for the z-tubes (ztube10) and (iv) increased orbital bias factors by a factor of 10 for the x-tubes (xtube10). Figure 13 shows the resulting anisotropy-profiles (when averaging over shells) for these models. The different rows correspond to the different input toy models and the different columns show the individual projection directions. The dotted data points symbolise the anisotropy profiles of the input toy models (which are of course the same for the different projections). The colored lines (green, blue, orange, red) show the recovered anisotropies of SMART fits using entropy functions with enhanced bias factors for specific orbit types as described under (i)-(iv) above. The grey lines, for comparison, show the anisotropies that result when we maximise the above entropy functions (i)-(iv) without fitting the mock LOSVDs of the toy models. The grey lines therefore illustrate the variety of different anisotropy profiles that can be constructed by varying the orbital bias factors ω i . They also indicate the range of different anisotropy profiles that are consistent with the given density distribution. As one can see, even though the entropy functions tend to push the fits into extreme directions, the range of anisotropy profiles recovered after the fit to the LOVSDs is very narrow. When averaging over all radii and toy models fitted with different entropies, the mean deviation to the input models (dots) is |∆β| = 0.05. The average spread in beta inside the sphere of influence is slightly larger |∆β(r < r SOI )| = 0.09 than |∆β(r SOI < r < r FOV )| = 0.03 outside r SOI . One possible explanation for this might be the increase of degrees of freedom of spherical orbits near the center. In addition, because β involves the ratio of the intrinsic dispersions, the same fractional error in the intrinsic dispersions results in a 4 times larger |∆β| when the anisotropy is as tangential as β = −1.5 compared to the isotropic case. Overall, the Shannon entropy (red line) is able to recover the beta anisotropy best. This is our default entropy used in SMART.
These results together with Section 4 strongly suggest that the information contained in the full LOSVDs constrains the anisotropy in the second-order velocity moments very well. In turn, this is the reason why our models can reproduce the mass of the black hole and of the stars in the N-body simulation very well. At larger radii, solely the reconstruction of the intrinsic anisotropy of the toy model with enhanced prograde z-tubes (i.e. ZROT) turns out to be difficult when viewed along the minor axis. This, however, is expected since any rotation around the minor-axis can not be observed and, thus, not be reconstructed from this viewing direction. Since we use equal ω i for prograde and retrograde orbits in our fitted models, these models do not have intrinsic rotation in the z-tubes for this projection. Consequently, the tangential velocity dispersion is larger than in the toy model and the fitted β becomes too negative. We checked that if we use the true second-order velocity moment rather than the velocity dispersion in the tangential direction, then the differences between the outer profiles of the ZROT model and the fits along the minor axis disappear.
For dynamical models which aim for a full phase-space reconstruction (like Schwarzschild models) and which use the full information encoded in the LOSVDs (see also Vasiliev & Valluri 2019) the anisotropy should be recoverable with a typical error of |∆β| ≈ 0.05. We solely found larger anisotropy discrepancies (up to |∆β| = 0.5) in extremely tangentially biased regions inside the sphere of influence. As an example of the corresponding accordance of the reconstructed orbit fractions, Figure 14 shows the case of z-tubes. The intended overpopulation of the z-tubes in the ZTUBE toy model (third row) and the ZROT toy model (fourth row) can be clearly seen in comparison to the other toy models. Independent of the chosen line of sight and entropy method (i.e. the bias factors ω i ), the fraction of z-tubes is qualitatively recovered and follows the enhancement tendencies. The orbit fractions however are less well determined than the anistropy by the data and show a stronger dependence on the entropy. The same is true for the x-tubes, box/chaotic and spherical/Kepler orbits shown in Figures E1 to E3.

INTERMEDIATE AXIS ROTATION
One of the first investigations to check whether tube orbits with net rotation around the intermediate axis are stable in a triaxial ellipsoid was done by Heiligman & Schwarzschild (1979), who studied a triaxial model with fixed axis ratios of 1:1.25:2 by numerical methods and stated that âĂĲY-tube orbits nearly certainly do not exist in the adopted modelâĂİ. Also Binney & Tremaine (2008, p. 263) assert that âĂĲtube orbits around the intermediate axis are unstableâĂİ in a triaxial potential. Adams et al. (2007) analysed the orbit instability of orbits in a triaxial cusp potential, which are initially confined to one of the three principal axes, under a perturbation along the perpendicular direction. They found that orbits around any of the principal axes are unstable to perpendicular motions. However, according to previous results, they again state that orbits around the intermediate axis are more likely to be unstable. This instability is strongest for original box orbits lying in the x-z plane when the axis ratio of these two axes in the original plane is largest. Our orbit classification routine in SMART finds no y-tubes in the sense that there is no sign conservation of the angular momentum along the intermediate axis over our default integration period of 100 SOS-crossings, agreeing with the aforementioned works done by other groups. However, some of the orbits integrated for the N-body models do show y-rotation for a limited time-span. When providing SMART with artificial projected input line-of-sight velocity distributions that mimick a net rotation along the intermediate axis, the model is able to produce a y-rotation signal of the order of 10 km s −1 in the fit. A rotation of this magnitude is small but in principle detectable with today's telescopes' resolution. Fig. 15 shows the SMART fit to the major-axis projection of the simulation (cf. Fig. 3) when assuming the viewing angles of the minor-axis projection. With this, we simulate a hypothetical rotation along the intermediate axis. The top row shows the velocity-and h 3 -map when stopping the orbit integration after 2 Gyrs if this is shorter than the time needed for 100 SOS-crossings. Indeed, the model reproduces a y-rotation signal. The amplitude of this residual y-rotation becomes smaller and smaller when the orbital integration time is increased, and vanishes when all orbits are integrated for 1000 SOS-crossings (bottom row). This analysis indicates that the model's triaxial potential, which is constructed based on the 3D density from the realistic N-body merger simulation, -contrary to expectations-contains orbits with y-rotation for a physically relevant time span. What remains unclear at this moment is whether such a y-rotation indeed appears in real elliptical galaxies. If so, then the connection between kinematic misalignments and photometric twists is less constrained than often assumed when only rotation around the intrinsic long and short axes is considered.

Remaining sources of systematics
All relevant properties of the simulated merger remnant galaxy were proven to be recovered with a convincing precision and the deviations we found are almost negligible. The remaining deviations (in the ∼ 5 − 10% level) can be either originated by SMART or the simulation. One remaining contribution to the scatter in the final mass recovery certainly comes from the finite binning resolution of the simulation data and of the SMART models (see Sections 2 and 3). Especially the need of extrapolating the density towards the center due to limited resolution of the simulation holds uncertainties. One more inaccuracy is potentially induced by the softening length in the simulation used to avoid unrealistic 2-body encounters between massive particles. Force calculations for radii smaller than the softening length are consequently modified by the softening. In the close vicinity of the black hole, the stellar particles in the simulation are modeled using a nonsoftened algorithmic chain regularization technique (ARCHAIN; Mikkola & Merritt 2006 including post-Newtonian corrections (e.g., Will 2006). The particles outside this chain radius are treated by using softened gravitational force calculations based on the GADGET-3 (Springel 2005) leapfrog integrator. The chain radius r chain is chosen to be at least 2.8 times larger than the GADGET-3 softening length (r chain > 2.8 ) to ensure that the particles within the chain remain nonsoftened. The softening length in the simulation for the stellar particles is * = 3.5 pc and the softening length for the dark matter particles is DM = 100 pc. The mass recovery of the global stellar mass-to-light ratio and dark matter scale factor will be unaffected by these relatively small values. However, there might be a remaining influence on the black hole mass recovery within r SOI . Nevertheless, this effect is expected to be small.

Comparison to other triaxial Schwarzschild Models
We found two other dynamical modeling codes in the literature using Schwarzschild's orbit superposition technique dealing with triaxiality by van den Bosch et al. (2008) and Vasiliev & Valluri (2019). The code FORSTAND by Vasiliev & Valluri (2019) is applicable to galaxies of all morphological types. When assuming that the deprojection and dark matter halo is known and provided to the code, the models of noise-free axisymmetric disc mock datasets taken from N-body simulations showed very weak constraints on M BH : any value between zero and 5-10 times the true black hole mass was equally consistent with the data. Read from Fig. 2 in their paper, the stellar mass-to-light ratio showed a variation of around 20%. Jin et al. (2019) tested the triaxial Schwarzschild code by van den Bosch et al. (2008) by applying it to nine triaxial galaxies from the large scale, high resolution Illustris-1 simulation (Vogelsberger et al. 2014), which provides a stellar and dark matter resolution of ∼ 10 6 M . When fixing the black hole mass and allowing the model to deproject the mock dataset, the stellar mass within an average effective radius is underestimated by ∼24% and the dark matter is overestimated by ∼38%. Their averaged model results obtained from mock data with different viewing angles tend to be too radial in the outer regions with better anisotropy matches in the inner region. Of course, these results cannot be used for direct comparison due to a widely varying resolution and in case of the analysis by Jin et al. (2019) the deprojection probably causes the major deviation. However, SMART for sure is able to add further progress in modeling triaxial galaxies and convinces with its proved precision.

SUMMARY AND CONCLUSION
We have developed a new triaxial dynamical Schwarzschild code called SMART and tested its efficiency and reliability by applying it to an N-body merger simulation including supermassive black holes.  Figure 13: Recovery of the anisotropy profiles of different input toy models by SMART fits with different entropy functions. The input models (black dots) are constructed to have specific orbit classes overrated, resulting in different anisotropy profiles (the overrated orbit type is labelled on the y-axis of each individual row, details in the text). All the five toy models can be well recovered. This is independent of the choice of the projection axis (different columns; from left to right major, intermediate, minor and a diagonal axis) and from the assumed entropy function in the fit (colored lines). For example, the anisotropy profile of the toy model with an overpopulation of box/chaotic orbits (called BOX; second row) is well reproduced by models maximising the Shannon-entropy (red lines, labelled shannon) but also with other entropy functions that use a ten times higher bias factor for box-orbits (green, labelled box10), for x-tubes (blue, xtube10) or for z-tubes (orange, ztube10). The grey lines correspond to the anisotropies implied by maximising these four different entropy functions without fitting the kinematic data. They symbolise the variety of anisotropy profiles which are in principle possible for different choices of ω i . After fitting the kinematic data, the average deviation (averaged over all radii and toy models fitted with different entropy techniques) between recovered and input anisotropy is very small, |∆β| = 0.05.
The simulation was deliberately selected due to its high accuracy, reasonable formation process, realistic internal structure and ability to precisely calculate the dynamics close to the central black hole. This ensured the possibility to check whether SMART is able to recover all relevant properties including the mass of a supermassive black hole of a realistic triaxial galaxy when providing the deprojected light profile and normalized DM halo. SMART is assembled with the feature to compute the potential and force by expansion into spherical harmonics allowing to deal with non-parametric densities and halos. Its orbit library contains 50000 integrated orbits which are set up by creating random initial radial and velocity values within given energy shells and angular momentum sequences and by filling the surfaces of section. This ensures the ability to adapt itself to a radially changing number of integrals of motion. The orbit superposition is executed by maximizing an entropy-like quantity and by using the full line-of-sight velocity distributions instead of only Gauss-Hermite parameters alone. These benefits enable SMART to reconstruct all relevant properties and features of the merger remnant with an excellent precision. We will now recap the requirements which were set on SMART and proved to be fulfilled in this analysis: -SMART is able to reproduce the anisotropy profile and internal velocity dispersions with an rms σ of only 1.2%.
-SMART reproduces the stellar mass-to-light ratio, -black hole mass, -and mass scale factor of the dark matter density profile with a precision on the 5-10% level. To our knowledge this is the first time that the intrinsic precision at given deprojected light profile has been quantified to be so high. This sets the basis for further investigations of the whole modeling procedure.  Figure 14: Recovery of the z-loop orbit fractions of different input toy models by models with different entropy methods. We here show the same analysis as in Figure 13 but now for the reconstruction of the fraction of orbits classified as z-tubes. The color coding is adapted to Fig. 4 and C1. Independent of the tested projection, the z-loop fractions of the individual input toy models are well recovered by the models using different entropy methods. The same is true for the other orbit class fractions (see Fig. E1-E3).
precision can be achieved for noiseless as well as noisy and psf convolved kinematic input data.
We extensively discuss that the maximum-entropy technique provides an elegant technique to study the range of possible orbit distributions consistent with a given set of data. Our tests with the N-body data and with additional toy models strongly suggest that when the full information contained in the entire LOSVDs is used to constrain the model, then the remaining degeneracies in the recovery of the exact phase-space distribution function do not affect 'macroscopic' properties of the galaxy models, like the anisotropy in the second-order velocity moments. This is the basis for the very good reconstruction of the orbital structure and mass of the black hole and stars with SMART.
It was shown that the orbit library is robust against axis changes and generates a complete set of well superpositioned orbits necessary to model a triaxial galaxy with all corresponding internal structures. Also the accurate mass parameter recovery accomplished by SMART suggests only minor degeneracies contained in the projected kinematic data, provided that the deprojection is known. We showed that these remaining minor degeneracies could in principle be narrowed even more if information about the orbital bias factors from a second projection direction were provided.
When analysing the elevation and azimuthal internal velocity distributions of the simulation we find that the central radial bins show two maxima. This corresponds to the negative h 4 -parameter in the center and the strong tangential anisotropy produced during the core formation. SMART is able to reconstruct this phenomenon with an accuracy of ∼ 7%.
One more discovery of scientific interest is intermediate axis rotation which is produced by orbits contained in the model's orbit library representing the simulation's triaxial potential. Independent of the question whether such intermediate axis rotation really appears in the real universe, it was shown that our model contains orbits with y-rotation whose stability was empirically found to be maintained up to at least 2 Gyrs.

APPENDIX A: SURFACE BRIGHTNESS-, ELLIPTICITY-AND POSITION-ANGLE-PROFILE
When fitting elliptical isophotes (Bender & Moellenhoff 1987) to the surface brightness maps for the principal axis directions of the simulation in units of stellar simulation particles one can see that all three projections differ in their surface-brightness-(SB) and ellipticity-profile ( = 1 − b/a with a being the semi-major and b being the semi-minor axis of the elliptical isophotes). The positionangle-profile (PA) shows some isophotal twists. Overall, the simulation appears to show typical triaxial behavior (see Fig. A1). Fig. B1 and B2 show the velocity maps of the kinematic input data and their fits by SMART for four different projection axes and an additional psf convolved version of the intermediate axis. The plots also contain surface brightness maps in form of stellar simulation particles and χ 2 -maps as deviation between the simulation data and modeled fit. As already described in Section 4.1 one can see that SMART manages to fit the kinematic input data independent of the chosen projection.

APPENDIX C: ORBITAL REPRESENTATION OF THE PHASE SPACE
For generating our orbital initial conditions we set up energy shells by calculating the maximum allowed potential energy for N E radii between the minimum allowed pericenter chosen to be r peri,min = r min and the maximum allowed apocenter of r apo,max = r max . The radial bins are sampled by using formula 4. Each energy shell i E ∈ [1, N E ] is then subdivided in i E (L z ) sequences. We empirically find i E (L z ) = i 2 E as appropriate choice for the number of sequences per energy shell. The orbits for a specific sequence share the same initial z-component of the angular momentum L z,i E . These L z,i Esequences are logarithmically sampled between the minimum and maximum allowed L z -value in dependence of the radius r(i E ) of the specific bin. The final orbital initial spatial conditions (r i , θ i , φ i ) are then given -by randomly sampling azimuthal angles φ i ∈ [0, 2π], -by setting the elevation angles θ i = 0 because every orbit crosses the equatorial plane and -by randomly sampling possible radial values r i between the minimum and maximum allowed apsis, i.e., r i ∈ r peri (L z,i E ), r apo (L z,i E ) . While sampling the energy as an integral of motion is straightforward, the sampling of L z is an arbitrary way to qualify the azimuthal velocity v φ as the parameter we actually want to sample. In analogy, one could also sample the elevation velocity v θ in terms of the x-component of the angular momentum L x , which is discussed in more detail in Section 4.5.4. The computation of r peri (L z,i E ) and r apo (L z,i E ) is thereby done under the wrong assumption that L z,i E was conserved. Strictly speaking, this is only true for an axisymmetric potential as limiting case. However, our tests above show that this does not bias the results obtained with our orbit library. The remaining initial velocity conditions (v r,i , v θ,i , v φ,i ) can then be subsequently computed as where E is the total energy and Φ the potential and r,i . More initial conditions are generated by filling the surfaces of section (SOS; Henon & Heiles 1964;Richstone 1982). The surface of section, also called PoincarÃľ section, consists of the r i -and v r,i -values of the orbits crossing the equatorial plane, i.e., θ i = 0, in upwards direction. The topology of the SOS depends on the integrals of motion. The orbit library representatively fills the phase-space if the SOS is properly filled. For this, we randomly generate azimuthal angles φ i which get divided in N sector azimuthal sectors. For each azimuthal sector we then sample (r i , v r,i )-tuples which are chosen so that their distance to the nearest imprint-tuple is maximised. The remaining initial conditions v φ,i and v θ,i are generated as before.
To create a homogeneous distribution of pro-and retrograde orbits, the orbits contained in the library which are later on classified as regular orbits (see Section 2.3.1) get duplicated by changing the sign of their direction of rotation. For an illustration of our 5D starting space of the stellar orbits in Fig. C1 we show a series of 2D plots with velocity spheres on annulus sectors (chosen radial and azimuthal sections) in the equatorial plane (i.e. θ = 0). Since the magnitude of the velocity vector v 2 x + v 2 y + v 2 z is determined by the energy E, it is sufficient to plot v x , v y -imprints for a certain energy shell (the remaining v z -component is then given as velocity shell in dependence of the value of E). We show both, the velocity components v x , v y of the orbital initial starting point (marked as thick dots) as well as the orbital imprints for every crossing event in the subsequent time evolution (thin dots). The velocity components are normalized by v max (r, θ, φ) = 2(E − Φ(r, θ, φ)), with Φ being the gravitational potential. The different colors correspond to the different orbit types. Fig. C1 shows that the orbits of our orbit library homogeneously represent the accessible phase space. A sufficient coverage of the phase space is crucial for a reliable model. In addition, one can see that the actual starting points of the orbits can be chosen in very different ways. Our random choice is only one possibility that leads to an appropriate phase-space representation. Most of the space gets automatically filled by the orbital imprints over time.

APPENDIX D: OPTIMIZATION ALGORITHM AND TEST
To find the unique set of orbital weights that maximisesŜ our code conceptionally seeks for a solution subject to the constraints (19). Technically, we follow Richstone & Tremaine (1988) and treat the orbital weights ì w and the kinematical model predictions as independent variables. The actual dependency ì y = ì y( ì w) (Equation 16) is added to the constraints.
Thus, let be the N orbit + N data new variables, then the model solves for the x i and for N data + N phot Lagrange multipliers λ j (combined in the vector ì λ). The matrix C ≡ P orb 0 L orb 1 (D5) contains the constraint equations Equation D4 is nonlinear and solved with a Newton method: or, At each iteration, the Lagrange multipliers are updated via (D6), which leads to The Jacobian dg dx is diagonal and its inverse is easy to compute. The only matrix that needs to be inverted numerically is A and its dimension is N data + N phot < N orbit .

Modulo the scaling by
The required linear algebra operations are implemented through the LAPACK library (Anderson et al. 1999).
In practice, to find the best-fitting orbit model, we start with a very small value α 0 and then iteratively increase α until χ 2 does not change anymore. Everytime we change α → α + dα the unique solution to Equation 7 will shift slightly from ì w to ì w + d ì w in the space of the orbital weights. Because f ≡ −Ŝ is differentiable with respect to α we can keep the new solution ì w + d ì w in an arbitrarily small neighborhood around ì w by limiting dα accordingly. For strictly convex f , the Newton method converges locally (and the convergence is quadratic). This means that once we have found a solution for some value of α, the convergence of the Newton steps for neighbouring α + ∆α is guaranteed (and fast) if we use the solution at α as starting point for the search at α + ∆α: it is only a question of limiting ∆α, if necessary. In this way, we can iteratively solve for any value of α starting from α 0 . In practice, we increase α by a factor of 2 in each iteration. Intermediate steps with a smaller α increment are inserted in case an iteration has not converged (which almost never happens).
At the initial α 0 we start from a homogeneous distribution of orbital weights. Methods with guaranteed convergence even from starting points outside of the constraint conditions have been developed. We do not use such methods, because in practice we always find the solution of Equation D4 even from an initially homogeneous distribution of orbital weights. Since the solution is global and unique for each α, convergence at any step of the α sequence can be tested based on Equation D4 itself.
To test our implementation we reconstruct the kinematical data of a toy model with randomised orbital weights. The toy model is constructed as a maximum-entropy model (i.e. through maximisation ofŜ in eq. 7 with α = 0). The bias factors ω i for the orbital weights are set to ω i = 10 r i , where r i ∈ [0, 3] are random numbers. After maximisation ofŜ subject to the density constraints, the resulting orbital weights w i satisfy the density constraints of the N-body simulation. Due to these density constraints, the w i are not directly proportional to the random ω i , yet they are still very strongly affected by them and thus we call this toy model with randomised orbital weights RANDOM in the following. We adapt r peri,min = r min = 0.05 kpc and r apo,max = r max = 21 kpc. If an orbit conserves the sign of at least one angular-momentum component, then this orbit has a duplicate companion orbit with opposite direction of rotation in the library (this concerns tube orbits and spherical orbits; see also Section C). Let i + and i − denote the orbital indices of such a pair of orbits, where i + refers to the orbit with initial velocity components v r , v θ , v φ > 0. To increase the degree of difficulty for reconstructing the RANDOM model, we modify the weights of orbits in such a pair: By using c w = 10 −3 , Equation D14 leads to increased weights w i + ,ran and almost eliminated weights w i − ,ran , corresponding to an internal rotation of the model (see also Thomas et al. 2007). We then modelled the LOSVDs of the major axis projection of this RANDOM toy model with SMART in order to test our solver algorithm for the orbital weights. Because the χ 2 minimisation is non-unique (cf. Section 5), the solution of eq. 7 will in general not occur at the true weights w i,ran of the RANDOM model. We therefore tried two different configurations: (i) we used the default orbital bias factors ω i = 1 (Shannon entropy) and (ii) we set ω i = exp(1) · w i,ran . In the second case, we push the maximum entropy solution towards the original, randomly generated weights. Only in this case, we can predict the location of the global solution analytically and test our implementation rigorously. For ensuring fair conditions we set the initial weights to w i,ini = 1 in both cases when we start the minimisation. An efficient optimization algorithm should be able to find two equivalent degenerate solutions, one for case (i) and a different one for case (ii). The true orbital weights of the RANDOM model should be recovered (only) in the latter case.
SMART is indeed able to fit the input data for both cases very well. We stopped the minimisation atχ 2 = 2 × 10 −9 for ω i = 1 andχ 2 = 5 × 10 −9 for ω i = exp(1) · w i,ran . Figure D1 shows the comparison between the true orbital weights of the RANDOM model and the recovered weights when using ω i = 1 (case (i) above). As expected, these weights differ from the RANDOM weights, since the χ 2 minimisation and, hence, the recovery of the entire phasespace distribution function, is non-unique. Orbits that contribute to a larger number of LOVSD data points and, in turn, are better constrained by the data have fitted weights that tend to be closer to the true ones. Figure D2 shows the orbital weights of the input model and the modelled fit using ω i = exp(1) · w i,ran . As explained above, in this case we expect the recovered solution to exactly match with the RANDOM model. Figure D2 demonstrates the high quality of our solver algorithm since the fitted weights indeed match exactly with the RANDOM weights (∆w = 4 × 10 −9 ).
The two different fits presented in Figs. D1 and D2 illustrate how the maximum entropy technique with variable ω i can be used to sample different solutions that are equivalent with respect to the χ 2 minimisation.
In Fig. D3 we compare the 'macroscopic' internal velocity structure of these two fits. Even though the detailed population of the orbits (i.e. the exact phase-space distribution function) of both models differ, they both reproduce the internal moments of the RANDOM model very well (rms σ = 5 × 10 −3 for ω i = 1 and rms σ = 3 × 10 −5 for ω i = exp(1) · w i,ran ). As further discussed in Section 5, this suggests, that when fitting the full LOSVDs the remaining degeneracies in the recovery of the distribution function of a triaxial galaxy (with known normalized densities) show only little effect on the 'macroscopic' galaxy parameters of interest.

APPENDIX E: ORBIT FRACTION RECOVERY WHEN TESTING THE UNIQUENESS OF THE ANISOTROPY RECOVERY
We here show the remaining orbit fraction recovery plots for the analysis described in Section 5.2. The recovery of the z-loops was already shown in Fig. 14 of that section. Here, we plot the recoveries of the fraction of orbits classified as x-loops (see Fig. E1), box/chaotic orbits (see Fig. E2) as well as spherical/Kepler orbits (see Fig. E3). Independent of the line of sight (different rows) or chosen entropy method (black lines), the orbit class fractions of the individual input toy-models (different rows) are well recovered. This paper has been typeset from a T E X/L A T E X file prepared by the author. Figure B1: Velocity and surface brightness maps of the simulation (top row) as well as velocity maps of the model and χ 2 -map as deviation from the kinematic input data with the modeled fit (bottom row) for different projections. The individual projection axis can be read from the title (for the major axis projection and a more detailed caption description see Fig. 3). Figure B1: (continued) Figure B2: Velocity and surface brightness maps of the simulation (top row) as well as velocity maps of the model and χ 2 -map as deviation from the kinematic input data with the modeled fit (bottom row) for the noisy and psf convolved intermediate axis projection (for a more detailed description see caption in Fig. 3 and Section 4.5.1). Figure C1: Orbital Phase-space coverage in the x − y plane. We here plot the normalised velocity components v x , v y of all orbits in a single energy shell. This includes the orbital starting points (thick dots) as well as SOS-crossings during the time evolution of the orbits (thin dots). To consider all five dimensions of our starting space we distinguish between several radial and azimuthal annulus sectors of the plane (the chosen radial and azimuthal intervals are labelled above the individual panels). The orbit library produces a representative coverage of the phase-space because the orbital imprints homogeneously fill the area enclosed by the unit circle.  Figure D1: Recovery of the orbital weights of the RANDOM test model described in Section D. Each dot represents one orbit with the color indicating the number of LOSVD data points N constraints that the orbit contributes to. For the recovery of the weights we fitted the kinematics of the RANDOM model using our fiducial constant bias factors ω i = const. (Shannon entropy). Due to the degeneracy in the χ 2 minimisation, the recovered orbital weights differ from the original weights. Larger differences are observed for orbits which are less constrained by data (blue points).  Figure D2: Same as Fig. D1, but here we fitted the RANDOM model using ω i = exp(1) · w i,ran . With these bias factors the unique maximum of eq. 7 should occur at the original weights w i,ran of the RANDOM model. And indeed, even though we start the fit with an initial guess where all the weights are equal, the true solution is reproduced with very high precision. This demonstrates the robustness and high accuracy of our optimization algorithm. The two different fits presented here and in Fig. D1 illustrate how the maximum entropy technique with variable ω i can be used to sample different solutions that minimise the χ 2 . fit : ω i = exp(1) · w i,ran Figure D3: Internal velocity structure of the two fits presented in Fig. D1 (black lines) and D2 (red lines). Despite the different orbital weight distributions, the anisotropy in the second order velocity moments is almost equally well recovered in both cases. This provides evidence that when using the full information contained in the entire LOSVDs, the remaining degeneracies in the recovery of the distribution function of triaxial galaxies (i.e. the differences in the detailed population of the various orbits) have only little impact on the 'macroscopic' internal galaxy properties like the anisotropy β. This, in turn, is crucial for a precise mass reconstruction.  Figure E1: Recovery of the x-loop orbit fractions of different input toy models by models with different entropy methods. We here show the same analysis as in Figure 13 and Figure 14 but now for the reconstruction of the fraction of orbits classified as x-tubes. Independent of the tested projection, the x-loop fractions of the individual input toy models are well recovered by the models using different entropy methods.  Figure E3: Recovery of the spherical/Kepler orbit fractions of different input toy models by models with different entropy methods. For a more detailed caption description see Fig. 13 and E1.