Triaxial Orbit-based Dynamical Modeling of Galaxies with Supermassive Black Holes and an Application to Massive Elliptical Galaxy NGC 1453

Most stellar-dynamical determinations of the masses of nearby supermassive black holes (SMBHs) have been obtained with the orbit superposition technique under the assumption of axisymmetry. However, few galaxies—in particular massive early-type galaxies—obey exact axisymmetry. Here we present a revised orbit superposition code and a new approach for dynamically determining the intrinsic shapes and mass parameters of triaxial galaxies based on spatially resolved stellar kinematic data. The triaxial TriOS code described here corrects an error in the original van den Bosch et al. code that gives rise to incorrect projections for most orbits in triaxial models and can significantly impact parameter search results. The revised code also contains significant improvements in orbit sampling, mass constraints, and run time. Furthermore, we introduce two new parameter-searching strategies—a new set of triaxial shape parameters and a novel grid-free sampling technique—that together lead to a remarkable gain in efficiency in locating the best-fit model. We apply the updated code and search method to NGC 1453, a fast-rotating massive elliptical galaxy. A full 6D parameter search finds p=b/a=0.933−0.015+0.014 and q = c/a = 0.779 ± 0.012 for the intrinsic axis ratios and T = 0.33 ± 0.06 for the triaxiality parameter. Despite the deviations from axisymmetry, the best-fit SMBH mass, stellar mass-to-light ratio, and dark matter enclosed mass for NGC 1453 are consistent with the axisymmetric results. More comparisons between axisymmetric and triaxial modeling are needed before drawing general conclusions.


INTRODUCTION
Elliptical galaxies exhibit a wide range of isophotal shapes and surface brightness profiles.There is an intrinsic uncertainty in inferring the 3D stellar luminosity density from the observed 2D isophotes on the sky.When stellar kinematics from spectroscopic observations are combined with photometric information, stronger constraints can be placed on the intrinsic 3D shapes of elliptical galaxies (e.g., Binney 1985;Franx et al. 1991).An idealized galaxy obeying exact axisymmetry would, by construction, have a regular surface brightness distri-bution without any isophotal twists and have perfectly aligned photometric and kinematic axes.Triaxial systems, on the other hand, can have isophotal twists, misaligned photometric and kinematic axes, and other spatially varying kinematic features absent in an axisymmetric system.This consideration led Binney (1985) to argue that triaxiality is common among elliptical galaxies.
Since then, a more detailed picture has emerged.Elliptical galaxies with lower stellar mass (M * 10 11.5 M ) tend to exhibit properties typical of axisymmetry (e.g., Emsellem et al. 2007;Weijmans et al. 2014;Cappellari 2016;Foster et al. 2017).Comparatively, elliptical galaxies with higher mass (M * 10 11.5 M ) typically exhibit photometric twists, slow or no rota-tion, and misalignments between the photometric and kinematic axes, suggesting triaxial intrinsic shapes (e.g., Veale et al. 2017bVeale et al. ,a, 2018;;Krajnović et al. 2018;Ene et al. 2018;Goullaud et al. 2018;Ene et al. 2020).Thus, it is vital to understand the role of triaxiality in dynamical galaxy modelling, particularly in studying massive elliptical galaxies and their central black holes in the local universe.
The most massive SMBHs observed in the nearby universe lie in centers of some of the most massive nearby elliptical galaxies (Ma et al. 2014).However, few triaxial SMBH mass (M BH ) measurements have been published thus far, perhaps because of the complexity in orbital structures, high-dimensional parameter space, and the associated computational cost required to model stellar orbits in triaxial potentials.To date, all published M BH measurements based on triaxial orbit modeling have been performed using the code initially presented in van den Bosch et al. (2008).This code was first applied to determine the intrinsic shapes and M BH of two fast-rotating elliptical galaxies M32 and NGC 3379 (van den Bosch & de Zeeuw 2010).In this work, M32 was found to be near oblate axisymmetry with M BH = (2.4 ± 1.0) × 10 6 M , fully consistent with M BH from earlier axisymmetric models (van der Marel et al. 1998;Joseph et al. 2001;Verolme et al. 2002).NGC 3379, on the other hand, was found to be moderately triaxial, and the inferred M BH = (4 ± 1) × 10 8 M was double the value derived from axisymmetric models (Gebhardt et al. 2000;Shapiro et al. 2006).In a subsequent application to the S0 galaxy NGC 3998 (Walsh et al. 2012), the best-fit model was found to be moderately triaxial although oblate axisymmetry was not ruled out.Feldmeier-Krause et al. (2017) applied the van den Bosch et al. (2008) code to the nuclear star cluster and SMBH at the Galactic center.The cluster shape was strongly triaxial, and the inferred M BH was consistent within 1 σ of the values inferred from the orbit of the S2 star (Gravity Collaboration et al. 2019;Do et al. 2019).
More recently, den Brok et al. (2021) used the van den Bosch et al. (2008) code to model PGC 046832.This galaxy exhibits dramatic twists, and the resulting models preferred strong variations in triaxiality.However, while axisymmetric models suggested a central black hole mass of 6 × 10 9 M , the triaxial models prefer models with no central black hole.Instead they report an upper bound on the central black hole mass of 2 × 10 9 M .This differs significantly from the value determined from axisymmetric models.
In addition to these published triaxial M BH values, the van den Bosch et al. (2008) code has been used to determine several M BH in the nearly axisymmetric limit (Seth et al. 2014;Walsh et al. 2015Walsh et al. , 2016Walsh et al. , 2017;;Ahn et al. 2018).It has also been used to estimate the intrinsic triaxiality of galaxies under the assumption of a fixed M BH (e.g., van den Bosch et al. 2008;Leung et al. 2018;Zhu et al. 2018a,b;Poci et al. 2019;Yang et al. 2020;Jin et al. 2020).
We have been revamping the van den Bosch et al. (2008) code for a systematic study of the SMBHs and other mass components in the ∼ 100 most massive local early-type galaxies in the MASSIVE survey (Ma et al. 2014).As a first step, we introduced a version of the code capable of achieving the exact axisymmetric limit (Liepold et al. 2020;Quenneville et al. 2021).The original van den Bosch et al. (2008) code was (intentionally) not built to respect axisymmetry, but it had been used to perform (nearly) axisymmetric orbit modeling, leading to unexplained inconsistencies when the resulting M BH values were compared to those from axisymmetric orbit codes (e.g., Ahn et al. 2018).Our axisymmetrized version of the code has bridged this gap and now enables dynamical modeling of galaxies using stellar orbits that properly obey axisymmetry.We applied our axisymmetrized code to NGC 1453, a fast-rotating elliptical galaxy in the MASSIVE survey, and obtained a significant detection of its SMBH with M BH = (2.9 ± 0.4) × 10 9 M (Liepold et al. 2020).Models without black holes were excluded at the 8.7σ level.
For clarity, we refer to the original code (which was unnamed) by the citation van den Bosch et al. (2008), and refer to our versions as the TriOS (Triaxial Orbit Superposition) code.
In this paper, we move beyond the axisymmetric limit of Quenneville et al. (2021), and present a triaxial version of the TriOS code and a first application of this code.This triaxial TriOS code differs in a number of major ways from the original van den Bosch et al. (2008) code.We have implemented these changes to correct a number of bugs and issues that we uncovered during extensive tests of the original code for triaxial potentials.As a start, we correct a major error in the orbit construction part of the code that incorrectly flips some velocity components for the tube orbits.Our tests indicate that for most viewing angles, correcting this mistake has a significant impact on the resulting orbital kinematics and galaxy model parameter recovery within the code.Other major changes include (i) modifying the acceleration table used for orbit integration to gain a significant speedup in runtime, (ii) resolving issues with insufficient orbit sampling that can result in spurious shape preferences, and (iii) using a more uniform mass binning scheme to eliminate frequent problems in sat-isfying mass constraints.Details of these changes are described in Section 4.
In addition to these code changes, we introduce a new set of shape parameters in this paper (Section 3) that are chosen to improve the efficiency of parameter searches in triaxial galaxy shapes and orientations.These parameters strike a balance between sampling in galaxy intrinsic shape and galaxy orientation, and result in fewer unrealistically flat galaxy shapes.To place these new parameters in context, we provide a summary (Section 2) of the parameters used in previous work to describe a triaxial galaxy's intrinsic and observed axis ratios, the relations of viewing angles and sky projections, and how an observed surface brightness is deprojected to obtain a 3D intrinsic shape within the TriOS code.
We apply our triaxial TriOS code to NGC 1453 in the final part of the paper (Section 5).Since triaxial modeling typically involves at least five parameters (three for shapes and at least two for mass parameters), we introduce an efficient new search strategy for sampling this multi-dimensional parameter space.This new strategy does not rely on direct grid searches used in previous orbit modeling studies.Instead, we apply nested Latin hypercube sampling to a 6D parameter space and are able to converge to a best-fit model for NGC 1453 with an order-of-magnitude fewer sample points.The resulting best-fit triaxial model is compared to the best-fit axisymmetric model from Liepold et al. (2020).

MODELING A TRIAXIAL GALAXY
In this section we summarize the information relevant for modeling a triaxial galaxy, e.g., coordinate systems, intrinsic and apparent shape parameters, viewing angles, and sky projections.

Intrinsic Shapes and Axis Ratios
To describe the 3D structure of a galaxy, we use a Cartesian coordinate system centered at the galaxy's nucleus, in which the x, y, and z axes are directed along the intrinsic major, intermediate, and minor axes of the galaxy, respectively.The z-axis is therefore the symmetry axis of an oblate axisymmetric galaxy, and the x-axis is the symmetry axis of a prolate axisymmetric galaxy.
It is convenient to use a different coordinate system to describe properties projected on the sky.We follow the standard practice and take the x and y axes of this coordinate system to be along the major and minor axes of the projected surface brightness distribution of a galaxy.The z axis is along the line-of-sight.
We use a, b, c to denote the lengths of the three principal axes of a triaxial ellipsoidal isodensity surface, assuming c ≤ b ≤ a.We use a and b to denote the lengths of the (observed) major and minor axes of the projected ellipse on the sky.Four useful axis ratios are where p is the intrinsic intermediate-to-major axis ratio, q is the intrinsic minor-to-major axis ratio, u represents a compression factor between the intrinsic major axis and the apparent major axis on the sky due to projection, and q is the flattening of the projected shape.These quantities obey the inequalities The upper and lower limits of u correspond to the intrinsic major axis lying in the plane of the sky (u = 1 or a = a) and the intrinsic intermediate axis lying in the plane of the sky (u = p or a = b), respectively.The commonly used triaxiality parameter is which ranges between 0 for an oblate axisymmetric shape (a = b), and 1 for a prolate axisymmetric shape (b = c), with values between 0 and 1 indicating a triaxial shape.

Viewing Angles and Sky Projections
A line of sight between an observer and a galaxy is specified by two viewing angles (θ, φ), where θ and φ are the usual polar angles in the galaxy's intrinsic (x, y, z) coordinate system.Thus, θ = 0 • is for a line of sight along the intrinsic minor axis (i.e., a face-on view down the z-axis), and θ = 90 • is for lines of sight in the x − y plane (i.e., an edge-on view with the intrinsic minor axis in the sky plane).Similarly, φ = 0 • is for lines of sight in the x − z plane (i.e., the intrinsic intermediate axis is in the sky plane), and φ = 90 • is for lines of sight in the y − z plane (i.e., the intrinsic major axis is in the sky plane).
Given a triaxial 3D density stratified on similar concentric ellipsoids, the viewing angle θ and φ are sufficient to project the 3D shape and determine the 2D projected coordinate system (x , y ).To de-project an observed 2D shape on the sky, however, a third angle, ψ, is needed to completely specify the intrinsic coordinate system.This third angle ψ specifies the remaining degree of freedom once θ and φ are fixed -a rotation of the galaxy around the line of sight.More precisely, ψ is defined as the angle between the y axis, and the line defined by the intersection of the x − y and x − y planes.When ψ = 0 • , the x − y plane and x − y plane intersect along the y axis; when ψ = 90 • , the x − y plane and x − y plane intersect along the x axis.
Together, the three angles (θ, φ, ψ) uniquely specify the orientation of the intrinsic axes with respect to the projected axes.If the 3D density is stratified on similar concentric ellipsoidal surfaces, the axis ratios (p, q, u) of Equation ( 1) can be uniquely determined from the projected surface brightness and (θ, φ, ψ) using the equations from Appendix A of de Zeeuw & Franx (1989).

Deprojecting Observed Surface Brightness
Within the TriOS code, the 3D stellar density distribution is described by a sum of multiple Gaussian components of varying widths and axis ratios using the Multi-Gaussian Expansion (MGE) scheme (Cappellari 2002).To determine these components, one first fits a 2D MGE to the observed surface brightness of the galaxy.Each MGE component is allowed to have its own projected flattening q to account for radially varying ellipticity in the observed isophotes.In addition, each MGE component can have a different position angle (PA) to accommodate any observed isophotal twists.
In general, the deprojection of a 2D surface brightness distribution to give a 3D triaxial luminosity density is not unique.MGE is a parametric method of choosing one particular 3D density for a given 2D surface brightness and set of intrinsic axes.Non-parametric deprojection methods have also recently been developed for triaxial galaxies in de Nicola et al. ( 2020), but the TriOS code is not yet capable of using these deprojections.
For a set of (θ, φ, ψ) that specifies the alignment of the galaxy's intrinsic principle axes (x, y, z), one can determine the deprojection of each MGE component that shares these principle axes (if a valid deprojection exists).This deprojection is unique due to the assumption that each 2D gaussian corresponds with a 3D gaussian density with similar concentric ellipsoidal surfaces of constant density.The axis ratios p and q of each deprojected MGE component can have their own values.The triaxiality parameter T , on the other hand, has the convenient property that it is identical for all MGE components when the components share the same PA (i.e., no isophotal twists).1

Prior practice
As discussed in Section 2, either (p, q, u) or (θ, φ, ψ) can be used to specify the shape of a triaxial galaxy and its sky projections.One can in principle search in either space when running orbit models to determine a galaxy's intrinsic shape and mass parameters.In practice, however, prior triaxial orbit modeling studies favored (p, q, u) over the angles.In these studies, the orbit models were typically run for a grid of regularly spaced values of (p, q, u) (e.g., van den Bosch & van de Ven 2009; van den Bosch & de Zeeuw 2010; Walsh et al. 2012;Jin et al. 2019).In a few other triaxial studies, u was fixed to some value close to 1 while the parameter search was conducted over p and q in a regular 2D grid (e.g., Zhu et al. 2018a,b;Poci et al. 2019).Since u ∼ 1 corresponds to the intrinsic major axis lying close to the sky plane, these studies did not search over all allowed viewing angles.
The argument used by van den Bosch & van de Ven (2009) for favoring conducting parameter searches in (p, q, u) rather than (θ, φ, ψ) is that a change in the angles can result in either a very small or very large change in axis ratios, depending on the angles being explored.We note, however, that the converse is also true: a change in the axis ratios can result in either a very small or very large change in the principal axes' alignment, depending on the values of these ratios.Two models with similar axis ratios, but viewed along very different lines of sight, can result in very different observables.An optimal sampling should consider both the intrinsic shape and the alignment of the line of sight.

Properties of new parameters
Here we propose a new set of variables to parameterize a galaxy's intrinsic triaxial shape and its sky projections.The advantages of conducting parameter searches in these variables over either (p, q, u) or (θ, φ, ψ) during triaxial orbit modeling will be discussed in Section 3.4.
For the first shape parameter, we choose the triaxiality parameter T (Equation 3).We define the next two parameters with forms analogous to T : where T maj parameterizes the length of the projected major axis, a , relative to its allowed limits a and b, and T min parameterizes the length of the projected minor axis, b , relative to its allowed limits b and c.It then follows from the inequalities in Equation ( 2) that (T, T maj , T min ) form a unit cube, i.e., (5) The limiting cases represented by each face of the unit cube has the following physical significance: (i) T = 0 and 1 correspond to oblate axisymmetric (a = b or p = 1) and prolate axisymmetric (b = c or p = q) shapes, respectively; (ii) T maj = 0 and 1 correspond to the intrinsic major axis lying in the sky plane (a = a or u = 1) and the intrinsic intermediate axis lying in the sky plane (a = b or u = p), respectively; (iii) T min = 0 and 1 correspond to the intrinsic minor axis lying in the sky plane (b = c or uq = q) and the intrinsic intermediate axis lying in sky (b = b or uq = p), respectively.While both the T maj = 1 and T min = 1 planes correspond to the intrinsic intermediate axis in the sky plane, they represent two complementary ranges of viewing angles such that b is equal to the projected major axis a for T maj = 1, whereas b is equal to the projected minor axis b for T min = 1.Equation ( 4), along with the requirement that q > 0, yields the inequality implying that for an observed axis ratio q on the sky, only the (T, T maj , T min ) region satisfying the inequality has valid deprojections.When the projected shape is flattened (q < 1), some models within the unit cube will result in negative (and thus invalid) values of the squared minor axis length, c2 .This volume surrounds the line (T, T maj , T min ) = (T, 1, 1), which does not have a valid deprojection for any flattened projected shape.

Relating (T, T maj , T min ) to old parameters
While Equations (4) relate our new parameters to (p, q, u), it is often useful to do the inverse and convert a given set of (T, T maj , T min ) to (p, q, u).To do so, we use these sequential expressions For a given set of (T, T maj , T min ), these equations define the deprojection from an observed MGE component with flattening, q , to its 3D shape parameters, (p, q, u).
Similarly, it is useful to convert (T, T maj , T min ) to the angles (θ, φ, ψ): , We choose to use the branch where 0 and 90 • ≤ ψ ≤ 180 • , though other equivalent branches exist as well. 2 The inverse expressions relating (T, T maj , T min ) and (θ, φ, ψ) are given in Appendix A. Equations ( 7) and ( 8), as well as Equations ( A1) and (A2), follow directly from the definitions in Equation (4), and the general expressions for the deprojection of a triaxial density that is stratified on similar, concentric ellipsoids (e.g., de Zeeuw & Franx 1989, discussed further in appendix A).Furthermore, since Equation (8) and its inverse Equation (A2) make no reference to the observed flattening, the same set of values (T, T maj , T min ) can be used for a density that is composed of multiple such components with different flattening values.Thus, (T, T maj , T min ) and (θ, φ, ψ) are simply different parameterizations of the same space.Equations ( 7) and ( 8), along with the equations listed in Appendix A, make no reference to the MGE formalism and are applicable to any triaxial system that meets these conditions.The existence and uniqueness of a valid deprojection are not affected by the choice of shape space parameterization outside the principal planes.
To illustrate the properties of T maj and T min , we plot a set of lines of constant T maj and T min in a galaxy's intrinsic coordinate system x, y, and z in Figure 1.The corner points (T maj , T min ) = (1, 0), (0, 0), and (0, 1) correspond to viewing angles along the short, intermediate, and long axes, respectively.The point (T maj , T min ) = (1, 1) represents a line of sight lying along the line θ = η = tan −1 ( T /(1 − T )) in the x−z plane, which only results in a valid model for round projected shapes.For flattened shapes, there are no valid deprojections for lines of sight within a solid angle surrounding this direction.This non-deprojectable region increases in size, as the projected shape becomes flatter.

Advantages of T, T maj and T min
The parameters T , T maj , and T min have a number of desirable properties.First, as Figure 1 illustrates, T maj and T min change relatively uniformly with the lineof-sight direction.This is in contrast to the axis ratio space, (p, q, u), in which tiny changes can result in large differences in the angles.For example, models with p = 0.99 and a fixed q would undergo a 90 • rotation in φ when u is varied from 0.99 to 1.
Similarly, the galaxy shape varies much more uniformly with (T, T maj , T min ) than with (θ, φ, ψ).Again, tiny changes in the latter can result in large differences in galaxy shape.For example, when an observed surface brightness (without isophotal twists) is deprojected into a 3D ellipsoidal shape with principle axes defined by (θ, φ, ψ) = (89 • , 45 • , 90 • ), Equation (A2) shows that the resulting 3D shape has T = 0, i.e., it is oblate axisymmetric.As ψ is increased from 90 • by only ∼ 1 • , however, the deprojected shape varies drastically, with oblate axisymmetry at ψ = 90 • to prolate axisymmetry at ψ ∼ 91 • , with the full range of triaxialities lying in between.From Equation (A2) with φ = 45 • , we find prolate axisymmetry (T = 1) to occur when ψ − 90 • = (90 • /π) arctan 2 cos θ/ sin 2 θ on our chosen branch.As θ approaches 90 • , the value of ψ that gives prolate axisymmetry approaches 90 • .For θ = 89 • (and φ = 45 • ), prolate axisymmetry occurs at ψ = 90.99985• .The behavior in the example above arises from coordinate singularities in the (θ, φ, ψ) space.When the line-of-sight is chosen to lie in a principal plane (i.e., cos (θ) = 0, 1 or sin (2φ) = 0), it is impossible for continuous photometric twists to arise in projection as triaxiality is varied.One consequence of this is that the only valid values of ψ are 0 • or 90 • , meaning it is no longer an independent parameter.Thus, (θ, φ, ψ) are insufficient to fully specify the 3D projection.The parameters (T, T maj , T min ), on the other hand, have no such singularity.In the above example, the proximity of the chosen value of θ to 90 • causes the rapid shift in shape with ψ.
Another desirable property of T maj and T min is that, similar to T (see Section 2.3), they do not vary among MGE components with different axis ratios, so long as there are no isophotal twists.
This invariant property can be explained by identifying T maj and T min as the shifted and rescaled versions of the conical coordinates, µ pro and ν pro within the galaxy's intrinsic coordinate system (Franx 1988), where µ pro = a 2 and ν pro = b 2 .Since the coordinate surfaces of µ pro and ν pro are the same for all MGE components, the shifted and scaled quantities T maj and T min do not vary between components.
The advantages of T , T maj , and T min are especially clear for systems not far from axisymmetry.Towards oblate axisymmetry (T ≈ 0), we have T min ≈ cos 2 θ and T maj ≈ cos 2 φ.Thus, a uniform sampling in √ T min and T maj will result in a nearly uniform sampling in the cosines of the inclination and the azimuthal angle.The same behavior holds towards prolate axisymmetry (T ≈ 1) since the roles of T maj and T min are simply switched if the x and z axis labels are interchanged.Thus, for nearly axisymmetric galaxies, a uniform sampling in (T, T maj , √ T min ) results in fewer unrealistically flattened models.

CODE CORRECTIONS AND IMPROVEMENTS
In this section, we describe the key corrections, improvements, and speedups made to the van den Bosch et al. ( 2008) code.See Section 4 of Quenneville et al. (2021) for other general changes that we had implemented (regardless of axisymmetry).

Correct orbital mirroring mistakes
The TriOS code is written for a static triaxial potential that is symmetric under reflection along each of the three principal axes of a triaxial system.Under this assumption, any orbital property only needs to be calculated in one octant of the orbit space; it can then be "mirrored" into the other seven octants by symmetry.2008) code.We plot the fractional error between the incorrect and corrected schemes (see Table 1) in the kinematic map of the line-of-sight velocity dispersion, σ, for a single orbit.The orbit is chosen from the x − z start space of a triaxial model with T = 6 × 10 −6 for NGC 1453, but it is representative of typical short-axis tubes in a triaxial potential.Each panel represents a different viewing inclination angle θ.The fractional error is largest near θ = 45 • , reaching beyond 50% for some parts of the orbit.
Taking advantage of this symmetry, the code initializes orbits in only one octant (x, y, z > 0) and integrates only these orbits.Seven additional copies of each orbit are then created by simply mirroring along the three axes.The recipe for how to flip the signs of the velocity components is given in Table 2 of van den Bosch et al. (2008).The exact procedure depends on whether the orbit is a short-axis tube, long-axis tube, or box.These orbits are classified as follow: throughout its trajectory, an orbit is labelled a box orbit if all three components of its angular momentum (L x , L y , L z ) change sign, and a tube orbit if exactly one component of angular momentum maintains its sign.The tube orbits are further classified according to the angular momentum component that maintains its sign, i.e., a long-axis (i.e.x-axis) tube maintains the sign of its L x , an intermediate-axis (y-axis) tube maintains the sign of L y , and a short-axis (z-axis) tube maintains the sign of L z .Orbits that don't fall into the tube or box orbit classifications are flipped in the same way as box orbits.
We discovered that the tube orbits are incorrectly flipped for four of the eight octants in Table 2 of van den Bosch et al. (2008).We indicate the incorrect components in boldface and give the corrected expressions in Table 1.The mistakes are such that the mirrored positions and velocities are inconsistent with one another, and the two do not combine to give a valid trajectory.A consequence of these mistakes is that the magnitude of each component of L is not always preserved by the mirroring, as it should be, and the resulting | L| is also not preserved.For instance, for the short-axis tube flip, the original recipe would change the amplitudes of L x and L y for 4 of the 8 copies, and the resulting total L would not be preserved.Similarly, L y and L z are incorrect for 4 copies of the long-axis tubes, and L x and L z are incorrect for 4 copies of the intermediate-axis tubes.
To illustrate the impact of the incorrect orbital flips, we plot the error in the line-of-sight velocity dispersion, σ, for a single short-axis tube orbit for three different viewing angles in Figure 2. We first integrate the trajectory of this orbit within the potential and then compute the 7 mirrored copies using the original and corrected flips in Table 1.The fractional difference in the projected σ between the two schemes is then plotted for three different values of viewing inclination angle θ.The errors vary across the plane of the sky and exceed 50% for θ = 45 • .For this orbit, the incorrect flip scheme tends to under-predict σ along the galaxy's projected major axis and over-predict σ near the edges.The orbit shown in Figure 2 is typical of short-axis tubes in triaxial potentials.Long-axis tube orbits exhibit similar error patterns when the appropriate axis labels are switched.While the pattern of velocity dispersion error is different for each orbit, systematic errors with magnitudes of 10% − 100% are typical, with peak errors of over 1600% in some cases for orbital inclinations near 45 • .
To assess further the impact of the incorrect flips, we perform full orbit modeling for a grid of triaxial models for NGC 1453 using the original and then the corrected scheme.Overall, when the correct flips are used, we find that χ 2 is lowered by a wide range of values depending on the triaxiality and viewing angles.For instance, the value of χ 2 can decrease by more than 100 for strongly triaxial models, while it can change by less than 5 or even increase slightly for other models.The overall χ 2 landscape is therefore significantly altered by our corrections.
Due to the symmetry of the tube orbits, the errors in the orbital flips can cancel out when the galaxy is viewed along a principal axis.Nearly axisymmetric models that are viewed edge-on or face-on will be similarly unaffected.Outside of these special cases, the orbital kinematics have significant errors.The incorrect flips were not used in our axisymmetric modelling of NGC 1453 (Liepold et al. 2020;Quenneville et al. 2021) since we used an axisymmetrization procedure in place of the flips in the TriOS code.
The discussion above is relevant only for tube orbits.For box orbits, we find the flips given in Table 2 of van den Bosch et al. (2008) to be correct.However, in addition to this set of 8 mirrored orbits, we choose to include 8 more orbits for each point in the stationary start space (defined in section 4.3) that correspond to enforcing time reversal symmetry for the box orbits.This addition ensures that box orbits have the expected even parity in their line-of-sight velocity distributions (LOSVDs).In the cases that we have examined, these orbits already have small enough odd LOSVD components that this change makes very little difference.

Modify acceleration table for significant speedup
In order to speed up orbital integration, the orbit code pre-computes a lookup table of acceleration values over a spatial grid and performs a trilinear interpolation to closely approximate the true acceleration.If an orbit passes outside the radial range of this grid, the acceleration is then computed from scratch, which is multiple orders of magnitude slower than interpolating values from the lookup table.It is therefore prudent to choose the extent of the grid wisely because even a small number of orbits passing outside the table's coverage can dominate the total runtime and unnecessarily increase the computation time of the entire orbit library.
We have noticed that some orbits can indeed pass outside the radial range used in the original code and result in a significant slow down.To eliminate this situation, we have made a simple modification to the radial range used for the acceleration table.In van den Bosch et al. (2008), the acceleration is pre-computed over a grid spanning the radial range where σ i is length of the semi-major axis of the ith projected MGE component, σ i is the length of the semimajor axis of the corresponding intrinsic MGE component, and r min and r max are the innermost and outermost orbital equipotential radii in the model.Thus, the lowest and highest energy orbits included in the model have energies Φ(x = r min , y = 0, z = 0) and Φ(x = r max , y = 0, z = 0), where Φ is the gravitational potential of the model.
In practice, we find that the second conditions in Equation ( 9) typically determine the range of the acceleration table, i.e., r interp,min = 0.01 r min and r interp,max = 1.05 r max .The outer boundary is never exceeded because energy conservation prevents orbits from passing outside r max and therefore r interp,max .The inner boundary of r interp,min = 0.01 r min , however, can be problematic because centrophilic box orbits can pass well within 0.01 r min .The DOP853 Runga-Kutta integrator in the TriOS code uses adaptive timesteps, tuning them to minimize errors in the position and veloc- ity between timesteps.In this scheme many acceleration evaluations are required in regions of the trajectory where the timestep is smaller, namely, when the trajectory passes closest to the central black hole where the orbits are most likely to reach below 0.01 r min .The fraction of acceleration evaluations within this boundary is somewhat model-dependent and may be higher when box orbits are launched from well within the SMBH's sphere of influence because the potential felt by those orbits is largely spherical and supportive of highly centrophilic box orbits.For a typical case of r min = 0.1 and r interp,min = 0.01 r min = 0.001 , we find as many as a sixth of the acceleration evaluations during the box orbit integrations to lie outside the lookup table.This minority of acceleration evaluations take up more than 50% of the total time when constructing the orbit library.
To enable a more efficient use of the acceleration table, we choose to decouple r interp,min and r interp,max from r min and r max which are used to determine the range of orbital energy sampling.When r interp,min is allowed to be smaller than 0.01 r min , we find the total time to integrate orbits can be reduced by a factor of a few, with a negligible change in accuracy.This speed-up is illustrated in Figure 3.As the acceleration table is extended to smaller radii, fewer orbits fall outside the radial coverage of the table, and the average integration time for box orbits drops significantly with decreasing r interp,min .For the example shown in Figure 3, choosing r interp,min ∼ 0.0001 would reduce the orbit integration time by a factor of 2 compared with the original setting of r interp,min = 0.01 r min = 0.001 .Since energy conservation prevents orbits from passing outside r max , setting r interp,max to be slightly larger than r max minimizes the integration time while maximizing the interpolation accuracy.
Since we don't typically vary the interpolation boundaries by more than 1 dex, the density of points in the interpolation grid does not change dramatically, and we find that the accuracy of the interpolated potential is sufficient.However, if the boundaries are changed more drastically, the number of radial interpolation points should be adjusted to maintain the desired accuracy.

Resolve issues with insufficient orbit sampling
The TriOS code samples orbit initial conditions from two separate spaces, referred to as start spaces (Schwarzschild 1993;van den Bosch et al. 2008).In the first start space ("stationary start space"), all orbits start from rest on the equipotential surface for a given energy.This start space contains only box, box-like, and chaotic orbits.
The second start space ("x − z start space") contains mainly tube orbits and samples orbits in the x−z plane, with velocity vectors pointing along the y-axis.As illustrated in Figure 4, orbits of a given energy in this space are sampled over the region bounded by the equipotential and thin-orbit curves.Typically, N I2 = 9 rays of orbits are sampled uniformly in polar angles from 0 to π/2 in the positive x and z quadrant; along each ray, N I3 = 9 orbits are uniformly spaced between thin-orbit curves and equipotential curve.Additionally, the code allows for dithering, where orbits with N dither adjacent initial conditions in each dimension are integrated and then bundled together to form each of the 9 × 9 orbits in order to improve phase space sampling.Figure 4 illustrates the case of (N I2 , N I3 , N dither ) = (9, 9, 3), where 27×27 tube orbits are launched in the positive quadrant of the x − z start space for a given energy.
For a triaxial model, the short-axis tubes (red points in Figure 4) and long-axis tubes (blue points) occupy two regions of the x − z start space separated by the focal curve.As derived in Appendix A of Quenneville et al. (2021), the focal curve is roughly approximated by Thus, as T increases from 0 to 1, the focal curve moves smoothly from the z-axis to the x-axis, and the composition of the tube orbits changes from being all short-axis tubes (for an oblate axisymmetric potential) to all longaxis tubes (for a prolate axisymmetric potential).
When orbits are well sampled, model properties such as the goodness-of-fit (χ 2 ) should vary smoothly as η (and hence T ) is varied.In our test runs for NGC 1453, however, we find that on top of a smooth variation, χ 2 varies periodically with T with a frequency matching the spacing between dithered orbits, (π/2)/N I2 , resulting in multiple spurious local minima at different values of T .Further testing reveals that these local minima arise from insufficient orbit sampling: as T increases, the focal curve approximated with η crosses rays of orbits in a periodic manner, resulting in the artificial oscillations in χ 2 with that same period.Since the periodic behavior is coherent as other model parameters are changed, it can have a significant impact on the recovered value of T and its uncertainty.Other parameter values are mainly impacted through their correlations with T .
We are able to eliminate the spurious oscillations in χ 2 vs. T by increasing N I2 , which increases the number of radial rays in the x − z start space and therefore improves the sampling in the polar angle.For the models presented in Section 5, we find that increasing N I2 from the default value of 9 to 15 and beyond removes the oscillations and also yields convergent results.We choose N I2 = 18 for the x − z start space.
We do not find similar issues for the other start space.Nonetheless, we increase N I2 to 18 for the stationary start space as well so as to maintain equal sizes for the tube and box orbit libraries.In summary, we use (N E , N I2 , N I3 , N dither ) = (40, 18, 9, 3) for both start spaces.This results in a total 40 × 18 × 9 × 3 3 × 3 = 524, 880 integrated orbits in each galaxy model, where the last factor of 3 accounts for the 3 orbit libraries (the x − z start space, its time-reversed copy, and the stationary start space).

Improve intrinsic mass binning scheme
In addition to kinematic constraints, the TriOS code enforces self-consistency of the mass model by requiring that the orbital weights be chosen to reproduce an input mass distribution (e.g., deprojected surface brightness profile of a galaxy).This is done by binning the mass in spherical coordinates (r, θ, φ), and requiring that the mass in each bin be reproduced to within a pre-specified precision (typically 1%).van den Bosch et al. (2008) uses linearly spaced bins between 0 • and 90 • for θ and φ, and logarithmically spaced bins between r min and r max /2 for r, where r min and r max are the innermost and outermost equipotential radii discussed in Section 4.2.
In the axisymmetrized TriOS code (Quenneville et al. 2021), we changed the radial binning scheme above to ensure sufficient orbits are used to represent the innermost and outermost mass bins.During our subsequent tests for triaxial systems, however, we noticed occasional problems with mass misfits in which a handful mass bins would have difficulty satisfying the 1% precision and/or contribute disproportionately high values to the total χ 2 of the galaxy model under examination.We are able to trace the problem to uneven bin sizes in θ used in the original code: the bins near the poles contained much less mass, as shown in the left panel of Figure 5.Because of this, the mass near the z axis was subject to The color scale here indicates the fraction of mass that falls within a given angular bin, summed over radius.The bottom row shows an example of the resulting χ 2 in the mass fits for a triaxial galaxy for the two binning schemes.The color scale here indicates χ 2 from attempting to fit a particular mass model, summed over radius.Only the 3D mass distribution is fit, with an error of 1% assumed on each bin.The most significant contributions to the mass χ 2 are from bins near the z-axis that contain very little mass.The triaxial mass model shown here has MBH = 2.9 × 10 9 M , M * /LF110W = 2.0, T = 0.10, q = 0.96q , Tmaj = 0.95, and Tmin = 0.12.much more stringent constraints than elsewhere, leading to frequent difficulties in satisfying the 1% fitting criterion.Even in the absence of kinematic constraints, spurious variations would arise in the χ 2 landscape, as illustrated in the right panel of Figure 5.The more oblate (T 0.1) and round (q 0.9q ) systems are more prone to this issue.
We find that this mass misfitting problem can be easily resolved by using mass bins linearly spaced in cos (θ) and φ, rather than in θ and φ.The resulting bins at a given radius then occupy the same volume, and the mass in each bin is much more uniform, with the binto-bin variations representing the galaxy's intrinsic deviation from spherical symmetry.Correspondingly, the pre-specified mass constraint criterion is enforced more uniformly.
For clarity, we have chosen to illustrate the mass misfitting issue in Figure 5 without imposing any kinematic constraints.When kinematic constraints are added in full orbit modeling (see Section 5), the total χ 2 returned by the code includes contributions from fits to the masses as well as kinematics.In this case, mod-els with significant mass misfits due to uneven binning schemes would have disproportionately larger χ 2 values, leading to potential biases in the recovered galaxy parameters.We apply the updated TriOS code described in the previous section to NGC 1453, a massive elliptical galaxy targeted by the MASSIVE survey (Ma et al. 2014).In Liepold et al. (2020), we performed orbit modeling of NGC 1453 using the axisymmetrized TriOS code.We refer the reader to that paper for a detailed description of the input kinematic and photometric data.In brief, the stellar kinematics are measured over 135 spatial bins from our high-spatial resolution Gemini GMOS IFS data (Ene et al. 2019(Ene et al. , 2020) ) and wide-field McDonald Mitchell IFS data (Veale et al. 2017b(Veale et al. ,a, 2018)).The first eight Gauss-Hermite moments are measured from the IFS spectra and used to constrain the stellar LOSVD in each kinematic bin; see Figure 4 of Liepold et al. (2020).
The MGE components representing the galaxy's mass distribution (see Section 2.3) are obtained from deprojections of our HST WFC3 photometry (Goullaud et al. 2018).Here we use the same input data but relax the assumption of axisymmetry in the orbit models.In order to ensure that all trajectories within the model are representative of their equilibrium distributions, we integrate each orbit in the x − z start space for 2000 times the orbital period for a thin tube orbit of the same energy.For orbits in the stationary start space, we integrate for 200 times the orbital period, as is typical of previous studies using the van den Bosch et al. (2008) code.
Due to the regular isophotes of NGC 1453 (Figure 5 of Liepold et al. 2020), we use the same PA for all MGE components and do not model isophotal twists.This is a common simplifying assumption (e.g., van den Bosch & de Zeeuw 2010; Walsh et al. 2012;Feldmeier-Krause et al. 2017) and it enables us to explore the galaxy's shape using the new scheme outlined in Section 3.
For the distance to NGC 1453, we adopt our new determination of 51.0 Mpc from the MASSIVE-WFC3 project (Goullaud et al. 2018) using the surfacebrightness fluctuation technique (Jensen et al. 2021).At this distance, 1 is 245 pc for a flat ΛCDM model with a matter density of Ω m = 0.315 and a Hubble parameter of H 0 = 70 km s −1 Mpc −1 .
We conduct the search for the best-fit galaxy shape in the new triaxial parameters (T, T maj , T min ) introduced in Section 3. The dark matter halo is modelled as a logarithmic potential.We parameterize it through its mass within 15 kpc, M 15 , which is roughly the central radius of the outermost kinematic bins, following Liepold et al. (2020).As in Liepold et al. (2020), we fix the scale radius of the dark matter halo to 15 kpc.Combining the three shape parameters with the three mass parameters M BH , M * /L F110W , and M 15 , we sample the 6D parameter space of galaxy models.
We determine the best-fit parameters by minimizing a χ 2 that includes terms for each LOSVD moment within each aperture, the projected light within each aperture, as well as the binned 3D mass density in order to enforce self-consistency for the stellar density.For each model, the best-fit set of weights are used to calculate the χ 2 differences between models.Lipka & Thomas (2021) recently suggested that recovery of the inclination of axisymmetric models can be biased unless the intrinsic flexibility of the models is accounted for.How-ever, a triaxial exploration of model flexibility is beyond the scope of the present study.
Instead of conducting model searches on a regular grid as was done in previous studies, we use the more efficient method of Latin hypercube sampling (McKay et al. 1979).There are many techniques for ensuring spatial uniformity in multidimensional spaces.We adopt the scheme described in Deutsch & Deutsch (2012), as implemented in the LHSMDU python package (Moza 2020).This procedure results in models that span a more continuous range of values than a regular grid, and are more uniformly spaced than random sampling.This approach allows a more representative sampling of the 6 dimensions with many fewer points than a regular grid.
We initially use a hypercube consisting of 1000 models spanning the range of M * /L F110W ∈ [1.7, 2.3], M 15 ∈ [3.5, 10.5]×10 11 M , and M BH ∈ [1, 5]×10 9 M , and the full range between 0 and 1 for (T, T maj , √ T min ).Of these models, 927 resulted in valid deprojections.We then use a rejection-based scheme to choose subsequent sets of model points.A Gaussian process interpolation of the 6-dimensional χ 2 surface is computed from the previously-run models.We use this interpolation to estimate the χ 2 for O(10 4 ) points chosen using the LHS scheme described above in the original volume and select points where the estimated χ 2 is within ∆χ 2 = 20.06(3σ for 6 parameters) of the estimated global minimum.To avoid premature optimization we perform this routine 10 times where random subsets of half of all previouslyrun models are used to build the interpolation function.With this scheme we select roughly 1000 model points which are expected to lie near the global χ 2 minimum to evaluate with the TriOS code.We perform two iterations of this rejection scheme, yielding roughly 3000 total model evaluations.
The resulting 6D likelihood landscape is shown in Figure 6.To determine the best-fit value and uncertainties, we fit the χ 2 landscape using Gaussian process regression with a squared-exponential covariance function (Pedregosa et al. 2011).To make the 2D contours shown in Figure 6, we transform this smoothed surface from (T, T maj , √ T min ) to (T, T maj , T min ), or (p, q, u).The marginalized 1D likelihood is also shown for each parameter.The shapes of the 2D contours in Figure 6 clearly demonstrate that (T, T maj , T min ) do not have the strong degeneracy apparent in (p, q, u).
The standard values of ∆χ 2 = 1, 4, 9 are used to define the 1σ, 2σ, and 3σ confidence intervals for 1 degree of freedom when considering the marginalized landscape for each variable individually.For the 2D contours, we use the values for 2 degrees of freedom, giving ∆χ 2 ≈ 2. 3, 6.2, 11.8.This is different from most previ-ous work using the van den Bosch et al. (2008) code, where typically ∆χ 2 = √ 2N obs is used to define the 1σ confidence interval, where N obs is the number of apertures on the sky, multiplied by the number of moments fitted within each aperture.This value is chosen to represent the intrinsic noise in the χ 2 values for each model, and is much larger than our values.However, while this is true when the input data are varied according to its noise level as discussed in Vasiliev & Valluri (2020), the noise level in the χ 2 values between models are significantly smaller when the input data are fixed.

Best-fit Triaxial Model
The best-fit values and the uncertainties for each NGC 1453 parameter are listed in Table 2.For each parameter, all other dimensions have been marginalized over.The best-fit M BH is consistent with the value determined from axisymmetric modelling in Liepold et al. (2020).The value of M * /L F110W has shifted down slightly, but is still consistent within 2σ of the axisymmetric value.
The best-fit shape, on the other hand, is inconsistent with axisymmetry.It is useful to compare our best-fit values of p = 0.93 and q = 0.78 with those inferred statistically from the observed distributions of ellipticity and misalignment angle between the kinematic and photometric axes for 49 slowly-rotating massive elliptical galaxies with measurable kinematic axes in the MAS-SIVE survey (Ene et al. 2018).In that sample, 56% of the galaxies have p > 0.9 with a mean value of 0.88, and the mean value of q is 0.65.Our best-fit shape for NGC 1453 indicates this fast-rotating galaxy is relatively oblate like the MASSIVE slow rotators and is slightly less flattened than the mean of that population.
The orbital composition of the best-fit triaxial model is shown in Figure 7 (top panel).Long-axis tubes and box orbits -two orbit types that are present only in triaxial potentials -together account for ∼ 30% of the orbital weights in the inner part and ∼ 45% in the outer part of NGC 1453.Quasi-planar orbits account for a small fraction of the total mass at small and large radii and are excluded from the plot.While long-axis tubes contribute a significant fraction of the mass, the projected model has fairly little minor axis rotation, due in part to the LOS being close to the intrinsic major axis.
The orbital velocity anisotropy of the best-fit model (bottom panel of Figure 7) is mildly tangential (β < 0) in the inner part and becomes increasingly radial outward.The radial profile has a similar shape to the axisymmetric model presented in Liepold et al. (2020).

Triaxial vs. Axisymmetric Best-fit Models
Table 2. Best-fit triaxial model parameters for NGC 1453 from the 6D likelihood landscape in Figure 6.For each parameter, all other dimensions have been marginalized over.The best-fit triaxial model presented above matches the observed kinematics significantly better than the best-fit axisymmetric model in Liepold et al. (2020).
Even though the best-fit χ 2 values in the two cases -493.0 for axisymmetric versus 382.7 for triaxial -differ by ∼ 110, they should not be compared directly because triaxial potentials require a new library of box orbits, and different numbers of orbits are used (6480 independent weights for axisymmetric versus 19440 for triaxial).Nonetheless, within triaxial modeling, our best-fit triaxiality of T = 0.33 is preferred over nearly oblate axisymmetric models with T ≈ 0 at a confidence level of about 5σ.To understand why non-axisymmetric models are favored, we examine the 2D maps of V and the lowest 3 even Gauss-Hermite moments in the GMOS data in Figure 8 (first row).We recall that axisymmetric models by construction produce only bisymmetric kinematics about the photometric major axis on the sky, meaning that the LOSVDs would be symmetric for points mirrored across the projected major axis and anti-symmetric for points mirrored across the projected minor axis.Any observed systemic deviation from bisymmetry would then indicate triaxiality.
For this reason, we decompose each GMOS moment map into a bisymmetrized component (second column) and a non-bisymmetrized component (third column).The latter exhibits clear systemic deviations from bisymmetry.The most obvious feature is the residual minor axis rotation indicative of kinematic misalignment.These maps assume a bisymmetrization along the projected photometric major axis used by our dynamical models, with a PA of 28.5 • .The residual pattern persists and can not be "rotated away" even if the PA is within uncertainties in the PA determination determined from the isophotal profile from (Goullaud et al. 2018).An axisymmetric model (consistent with the photometry) would be incapable of fitting these non-bisymmetric features in the data.To confirm this point, we plot the residual maps (fourth column) between the GMOS data and the best-fit axisymmetric model of Liepold et al. (2020).Indeed, the axisymmetric model exhibits similar residual patterns as in the data (third column).In comparison, the best-fit triaxial model is able to fit these non-bisymmetric features to a large extent, producing essentially random residuals (fifth column).
Figure 8 indicates that the preference for triaxiality is driven by the non-axisymmetric features in the NGC 1453 kinematics.Even though the nonbisymmetric features are somewhat subtle, they lead to detectable triaxiality, which we find to be best fit with p = 0.933, q = 0.779, and T = 0.33.Thus, despite be- ing a fast rotator with regular isophotal and kinematic features, NGC 1453 is best fit by a triaxial model.This is further evidence for widespread triaxiality in massive elliptical galaxies.
Importantly, however, the best-fit black hole mass M BH = 2.9 × 10 9 M is unchanged from that in the axisymmetric model.The stellar mass-to-light ratio and dark matter mass within 15 kpc agree to within a 1σ confidence level.

CONCLUSIONS
In this paper we have presented a revised code and a revamped approach for performing dynamical modeling of triaxial galaxies and their central SMBHs using the orbit superposition technique.We discussed a new triaxial version of the TriOS code that is capable of modeling triaxial systems while avoiding several shortcomings of the original van den Bosch et al. (2008) code.As a first application of this code, we performed triaxial orbit modeling of the massive elliptical galaxy NGC 1453 and presented the best-fit galaxy shape and mass parameters.This work complements Liepold et al. (2020) and Quenneville et al. (2021), in which we introduced a properly axisymmetrized version of the TriOS code.
We discovered and corrected a major error in the orbit kinematics in the van den Bosch et al. (2008) code: the tube orbits had wrong signs in certain mirrored velocity components in the orbit library (Table 1), resulting in incorrect projected kinematics.The magnitude of the kinematic errors varies spatially and depends on the viewing angles (Figure 2).This issue impacts all triaxial models that are not viewed along a principal axis, and all nearly axisymmetric models that are not viewed edgeon.How this error affects the best-fit galaxy shapes and mass parameters would have to be assessed on a galaxyby-galaxy basis by re-running the models with the corrected orbital flips in Table 1.In the case of NGC 1453, we find the χ 2 landscapes to be altered drastically, with χ 2 values changing non-uniformly by more than 100 for some models.
Following Quenneville et al. (2021), we continued to find ways to speed up the code.In this updated version of the TriOS code, we achieved another significant speedup (of up to ∼ 50 %; Figure 3) in orbit integration time by a simple extension of the interpolation table used to evaluate orbit accelerations (Section 4.2).The reduction in integration time is particularly pronounced for centrophilic orbits.
We have made two other adjustments in the code that significantly improve the sampling of long-axis tube orbits (Section 4.3) and enforce more uniformly the 3D mass constraints (Section 4.4).After these changes, the behavior of χ 2 vs. T (triaxiality parameter) no longer exhibits spurious oscillations, and the orbit code is able to find reasonable solutions for some mass models that were previously strongly disfavored.
The rest of this paper is devoted to new and improved strategies for searching the multi-dimensional parameter space required to specify triaxial galaxy models.We introduced a new set of shape parameters (Section 3) as well as a novel sampling technique (Section 4.3), which together lead to a remarkable gain in parameter searching efficiency.Searching in the new parameters T , T maj , and T min (Equation 3) avoids significant nonuniformities associated with other parameters used in earlier work.Our Latin hypercube sampling scheme results in an order-of-magnitude reduction in needed sampling points compared with conventional grid searches.
We applied the TriOS code and triaxial sampling scheme to the fast-rotating massive elliptical galaxy NGC 1453 in the MASSIVE survey (Section 5).NGC 1453 has a relatively small twist in the isophotes, and the kinematic and photometric axes are nearly aligned.Despite these properties that are typically in-voked to justify the use of axisymmetric orbit codes, we find the best-fit model to have a triaxiality value of T = 0.33, with intrinsic axis ratios p = 0.933 and q = 0.779.This best-fit triaxial model is able to match the observed kinematic maps significantly better than the best-fit axisymmetric model in Liepold et al. (2020).The improvement is mainly due to the ability of triaxial models to account for non-bisymmetric features in the data (Figure 8).Most other galaxies in the MASSIVE survey exhibit less (or no) rotation and more twists in their photometric and kinematic maps compared to NGC 1453.This is further evidence that massive elliptical galaxies have triaxial intrinsic shapes.M BH in the best-fit triaxial model for NGC 1453 is unchanged from the value measured with the axisymmetrized TriOS code from Liepold et al. (2020).Among the many dozens of stellar dynamical M BH measurements in local galaxies (e.g., McConnell & Ma 2013), NGC 1453 is only one of a handful galaxies whose central SMBH is studied with the full triaxial orbit modeling technique not limited to axisymmetry.In four other galaxies (Section 1), M32 had consistent M BH from axisymmetric and triaxial modeling, the NGC 3379 M BH increased by a factor of ∼ 2 when axisymmetry was relaxed, the PGC 046832 M BH decreased enough to be consistent with 0, while NGC 3998 was only modeled with the triaxial code so no comparison can be made.All four systems were modeled with the original van den Bosch et al. 2008 code, which used the incorrect mirroring scheme.Triaxial orbit modeling of more galaxies is needed for a full assessment of the systematic effects on stellar dynamical M BH measurements when the commonly-made assumption of axisymmetry is relaxed.
We thank Shaunak Modak and Jonelle Walsh for useful discussions.M.E.Q.acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), PGSD3-517040-2018.C.-P.M. acknowledges support from NSF AST-1817100, HST GO-15265, HST AR-14573, the Heising-Simons Foundation, the Miller Institute for Basic Research in Science, and the Aspen Center for Physics, which is supported by NSF grant PHY-1607611.This work used the Extreme Science and Engineering Discovery Environment (XSEDE) at the San Diego Supercomputing Center through allocation AST180041, which is supported by NSF grant ACI-1548562.

Figure 1 .
Figure 1.Isocontours of the new shape parameters, Tmaj and Tmin, in a galaxy's coordinate system, where the x, y, and z axes are chosen to be the intrinsic major, intermediate, and minor axes, respectively.The triaxiality parameter, T , is assumed to be 0.35 here.The parameters Tmaj and Tmin are seen to change relatively uniformly with the line-of-sight direction, resulting in fewer unrealistically flattened models near non-deprojectable regions (see text).

Figure 2 .
Figure 2. Illustration of the impact of the incorrect mirroring scheme in the van den Bosch et al. (2008) code.We plot the fractional error between the incorrect and corrected schemes (see Table1) in the kinematic map of the line-of-sight velocity dispersion, σ, for a single orbit.The orbit is chosen from the x − z start space of a triaxial model with T = 6 × 10 −6 for NGC 1453, but it is representative of typical short-axis tubes in a triaxial potential.Each panel represents a different viewing inclination angle θ.The fractional error is largest near θ = 45 • , reaching beyond 50% for some parts of the orbit.

Figure 3 .
Figure3.Average orbital integration time (per orbit) as a function of the inner interpolation radius, rinterp,min, used to tabulate the accelerations.The stationary start space contains mostly box orbits that pass near the galaxy center.The box orbit integration time increases drastically with rinterp,min, and the value used in the van denBosch et al. (2008) code is typically not small enough to minimize the integration time.

Figure 4 .
Figure 4.An example of the initial orbit locations in the x − z start space for a single energy value in the triaxial TriOS code.Orbits are launched from within the thin-orbit curve (inner grey arc) and equipotential curve (outer grey arc).The orbit initial conditions are sampled with NI 2 = 9 radial rays uniformly spaced in the polar angle from the z-axis to the x-axis, NI 3 = 9 points along each ray, and N dither = 3 to further improve the sampling, resulting in a total of 27×27 orbits.Each of the 27×27 color dots indicates the initial locations of an orbit (color coded by the type of orbits).The black line at angle η (see text) approximates the boundary between long-axis and short-axis tube orbits within this start space.Model χ 2 values are sensitive to the alignment between the angle η and orbit cell boundaries.

Figure 5 .
Figure 5.Comparison of the original (left) and new (right) mass binning scheme in the TriOS code.The top row shows that the bins near the x − y plane contain far more mass than the bins near the z axis due to the significant difference in bin volume in the original scheme (top left).Our new binning scheme evens out the mass considerably (top right).The color scale here indicates the fraction of mass that falls within a given angular bin, summed over radius.The bottom row shows an example of the resulting χ 2 in the mass fits for a triaxial galaxy for the two binning schemes.The color scale here indicates χ 2 from attempting to fit a particular mass model, summed over radius.Only the 3D mass distribution is fit, with an error of 1% assumed on each bin.The most significant contributions to the mass χ 2 are from bins near the z-axis that contain very little mass.The triaxial mass model shown here has MBH = 2.9 × 10 9 M , M * /LF110W = 2.0, T = 0.10, q = 0.96q , Tmaj = 0.95, and Tmin = 0.12.

Figure 7 .
Figure 7. Orbital composition (top) and velocity anisotropy (bottom) of the best-fit triaxial model of NGC 1453 as a function of radius.Short-axis tubes (solid) are dominant throughout the model, with significant contributions from long-axis tubes (dashed) and box orbits (dotted) that are present only in triaxial potentials.The velocity anisotropy parameter, β, has a similar radial profile for the best-fit triaxial (solid) and axisymmetric (dashed) models, being mildly tangentially anisotropic in the inner part and becoming more radially anisotropic in the outer part.

Figure 8 .
Figure 8. Maps of the stellar kinematics from the Gemini GMOS IFS in 135 spatial bins of the central 5 ×7 of NGC 1453.Four velocity moments are shown (from top down): V , σ, h4 and h6.The maps are oriented such that the horizontal and vertical axes are aligned with the galaxy's projected major and minor photometric axes, respectively.The data (first column) are decomposed into a bisymmetric component (second column) and a non-bisymmetric component (third column).To accentuate systematic patterns, we plot the non-bisymmetric component normalized by the moment uncertainty.Since an axisymmetric model can only produce bisymmetric kinematic maps, the residuals from the best-fit axisymmetric model (fourth column) show similar patterns to the bisymmetrized residuals.h6 shows additional residuals that are consistent with bisymmetry, but unable to be fit by an axisymmetric model.A triaxial model (right column) is able to capture most of the systematic behaviour in the input map, resulting in largely random residuals.The residuals have been normalized by the moment uncertainty.