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

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^{\prime}$- and $y^{\prime}$-axes of this coordinate system to be along the major and minor axes of the projected surface-brightness distribution of a galaxy. The $z^{\prime}$-axis is along the line of sight.

We use a, b, and c to denote the lengths of the three principal axes of a triaxial ellipsoidal isodensity surface, assuming c ≤ b ≤ a. We use $a^{\prime}$ and $b^{\prime}$ to denote the lengths of the (observed) major and minor axes of the projected ellipse on the sky. Four useful axis ratios are

$$\begin{eqnarray}&&p=\displaystyle \frac{b}{a},\quad q=\displaystyle \frac{c}{a},\quad u=\displaystyle \frac{a^{\prime} }{a},\quad q^{\prime} =\displaystyle \frac{b^{\prime} }{a^{\prime} },\end{eqnarray} \tag{ 1 }$$

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^{\prime}$ is the flattening of the projected shape. These quantities obey the inequalities

$$\begin{eqnarray}\begin{array}{rcl} & & 0\leqslant c\leqslant b^{\prime} \leqslant b\leqslant a^{\prime} \leqslant a,\\ \mathrm{or}\quad & & 0\leqslant q\leqslant {uq}^{\prime} \leqslant p\leqslant u\leqslant 1.\end{array}\end{eqnarray} \tag{ 2 }$$

The upper and lower limits of u correspond to the intrinsic major axis lying in the plane of the sky (u = 1 or $a^{\prime} =a$) and the intrinsic intermediate axis lying in the plane of the sky (u = p or $a^{\prime} =b$), respectively.

The commonly used triaxiality parameter is

$$\begin{eqnarray}&&T=\displaystyle \frac{{a}^{2}-{b}^{2}}{{a}^{2}-{c}^{2}}=\displaystyle \frac{1-{p}^{2}}{1-{q}^{2}},\end{eqnarray} \tag{ 3 }$$

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.

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^{\prime} ,y^{\prime} )$. To deproject 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^{\prime}$-axis and the line defined by the intersection of the $x^{\prime} -y^{\prime}$ and x−y planes. When ψ = 0°, the x−y plane and $x^{\prime} -y^{\prime}$ plane intersect along the $y^{\prime}$-axis; when ψ = 90°, the x−y plane and $x^{\prime} -y^{\prime}$ plane intersect along the $x^{\prime}$-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).

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^{\prime}$ 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. Nonparametric 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 principal 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). ³

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, 2018b; Poci et al. 2019). Because 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.

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:

$$\begin{eqnarray}\begin{array}{rcl}T & = & \displaystyle \frac{{a}^{2}-{b}^{2}}{{a}^{2}-{c}^{2}}=\displaystyle \frac{1-{p}^{2}}{1-{q}^{2}},\\ {T}_{{\rm{maj}}} & \equiv & \displaystyle \frac{{a}^{2}-a{{\rm{{\prime} }}}^{2}}{{a}^{2}-{b}^{2}}=\displaystyle \frac{1-{u}^{2}}{1-{p}^{2}},\\ {T}_{min} & \equiv & \displaystyle \frac{b{{\rm{{\prime} }}}^{2}-{c}^{2}}{{b}^{2}-{c}^{2}}=\displaystyle \frac{{\left({{uq}}^{{\rm{{\prime} }}}\right)}^{2}-{q}^{2}}{{p}^{2}-{q}^{2}},\end{array}\end{eqnarray} \tag{ 4 }$$

where T_maj parameterizes the length of the projected major axis, $a^{\prime}$, relative to its allowed limits a and b, and ${T}_{\min }$ parameterizes the length of the projected minor axis, $b^{\prime}$, relative to its allowed limits b and c. It then follows from the inequalities in Equation (2) that $(T,{T}_{\mathrm{maj}},{T}_{\min })$ form a unit cube, i.e.,

$$\begin{eqnarray}&&\begin{array}{l}0\leqslant T\leqslant 1,\\ 0\leqslant {T}_{\mathrm{maj}}\leqslant 1,\\ 0\leqslant {T}_{\min }\leqslant 1.\end{array}\end{eqnarray} \tag{ 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^{\prime} =a$ or u = 1) and the intrinsic intermediate axis lying in the sky plane ($a^{\prime} =b$ or u = p), respectively; (iii) ${T}_{\min }=0$ and 1 correspond to the intrinsic minor axis lying in the sky plane ($b^{\prime} =c$ or ${uq}^{\prime} =q$) and the intrinsic intermediate axis lying in the sky plane ($b^{\prime} =b$ or ${uq}^{\prime} =p$), respectively. While both the T_maj = 1 and ${T}_{\min }=1$ planes correspond to the intrinsic intermediate axis lying in the sky plane, they represent two complementary ranges of viewing angles such that b is equal to the projected major axis $a^{\prime}$ for T_maj = 1, whereas b is equal to the projected minor axis $b^{\prime}$ for ${T}_{\min }=1$.

Equation (4), along with the requirement that q > 0, yields the inequality

$$\begin{eqnarray}&&\displaystyle \frac{(1-T){T}_{\min }}{1-{{TT}}_{\mathrm{maj}}}\lt q{{\prime} }^{2},\end{eqnarray} \tag{ 6 }$$

implying that for an observed axis ratio $q^{\prime}$ on the sky, only the $(T,{T}_{\mathrm{maj}},{T}_{\min })$ region satisfying the inequality has valid deprojections. When the projected shape is flattened ($q^{\prime} \lt 1$), some models within the unit cube will result in negative (and thus invalid) values of the squared minor-axis length, c². This volume surrounds the line $(T,{T}_{\mathrm{maj}},{T}_{\min })=(T,1,1)$, which does not have a valid deprojection for any flattened projected shape.

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}_{\mathrm{maj}},{T}_{\min })$ to (p, q, u). To do so, we use these sequential expressions

$$\begin{eqnarray}\begin{array}{rcl}1-{q}^{2} & = & \displaystyle \frac{1-q{{\prime} }^{2}}{1-(1-T){T}_{\min }-q{{\prime} }^{2}{{TT}}_{\mathrm{maj}}},\\ 1-{p}^{2} & = & T(1-{q}^{2}),\\ 1-{u}^{2} & = & {T}_{\mathrm{maj}}(1-{p}^{2}).\end{array}\end{eqnarray} \tag{ 7 }$$

For a given set of $(T,{T}_{\mathrm{maj}},{T}_{\min })$, these equations define the deprojection from an observed MGE component with flattening, $q^{\prime}$, to its 3D shape parameters, (p, q, u).

Similarly, it is useful to convert $(T,{T}_{\mathrm{maj}},{T}_{\min })$ to the angles (θ, ϕ, ψ):

$$\begin{eqnarray}\begin{array}{rcl}{\cos }^{2}\theta & = & {T}_{\min }(1-T\,{T}_{\mathrm{maj}}),\\ {\tan }^{2}\phi & = & \displaystyle \frac{1-{T}_{\mathrm{maj}}}{{T}_{\mathrm{maj}}}\displaystyle \frac{1-{T}_{\min }}{1-{T}_{\min }(1-T)},\\ {\tan }^{2}\psi & = & \displaystyle \frac{\left[1-{T}_{\min }(1-T)\right](1-T\,{T}_{\mathrm{maj}})(1-{T}_{\min })}{{T}^{2}(1-{T}_{\mathrm{maj}}){T}_{\mathrm{maj}}{T}_{\min }}.\end{array}\end{eqnarray} \tag{ 8 }$$

We choose to use the branch where 0° ≤ θ ≤ 90°, 0° ≤ϕ ≤ 90°, and 90° ≤ ψ ≤ 180°, though other equivalent branches exist as well. ⁴ The inverse expressions relating $(T,{T}_{\mathrm{maj}},{T}_{\min })$ and (θ, ϕ, ψ) are given in the Appendix.

Equations (7) and (8) 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 the Appendix). Furthermore, because Equation (8) and its inverse given in the Appendix make no reference to the observed flattening, the same set of values $(T,{T}_{\mathrm{maj}},{T}_{\min })$ can be used for a density that is composed of multiple such components with different flattening values. Thus, $(T,{T}_{\mathrm{maj}},{T}_{\min })$ and (θ, ϕ, ψ) are simply different parameterizations of the same space. Equations (7) and (8), along with the equations listed in the Appendix, 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}_{\mathrm{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 $\theta =\eta ={\tan }^{-1}(\sqrt{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 nondeprojectable region increases in size, as the projected shape becomes flatter.

Zoom In Zoom Out Reset image size

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 line-of-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}_{\mathrm{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 $\psi -90^\circ =(90^\circ /\pi )\arctan \left(2\cos \theta /{\sin }^{2}\theta \right)$ 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 (\theta )=0,1$ or $\sin (2\phi )=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}_{\mathrm{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 ${\mu }_{\mathrm{pro}}=a{{\prime} }^{2}$ and ${\nu }_{\mathrm{pro}}=b{{\prime} }^{2}$. Because 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. Toward oblate axisymmetry (T ≈ 0), we have ${T}_{\min }\approx {\cos }^{2}\theta$ and ${T}_{\mathrm{maj}}\approx {\cos }^{2}\phi$. Thus, a uniform sampling in $\sqrt{{T}_{\min }}$ and $\sqrt{{T}_{\mathrm{maj}}}$ will result in a nearly uniform sampling in the cosines of the inclination and the azimuthal angle. The same behavior holds toward prolate axisymmetry (T ≈ 1) because 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,\sqrt{{T}_{\mathrm{maj}}},\sqrt{{T}_{\min }})$ results in fewer unrealistically flattened models.

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.

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 follows: throughout its trajectory, an orbit is labeled 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 do not 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 $\vec{L}$ is not always preserved by the mirroring, as it should be, and the resulting $| \vec{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 four of the eight copies, and the resulting total L would not be preserved. Similarly, L_y and L_z are incorrect for four copies of the long-axis tubes, and L_x and L_z are incorrect for four 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 seven 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 overpredict σ 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°.

Zoom In Zoom Out Reset image size

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 χ² is lowered by a wide range of values depending on the triaxiality and viewing angles. For instance, the value of χ² 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 χ² 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 modeling of NGC 1453 (Liepold et al. 2020; Quenneville et al. 2021) because 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 eight mirrored orbits, we choose to include eight 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.

In order to speed up orbital integration, the orbit code precomputes 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 slowdown. 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 precomputed over a grid spanning the radial range

$$\begin{eqnarray}\begin{array}{rcl}{r}_{\mathrm{interp},\min } & = & \min \left[0.1\times \min (\sigma {{\prime} }_{i}),0.01\,{r}_{\min }\right],\\ {r}_{\mathrm{interp},\max } & = & \max \left[6\times \max ({\sigma }_{i}),1.05\,{r}_{\max }\right],\end{array}\end{eqnarray} \tag{ 9 }$$

where $\sigma {{\prime} }_{i}$ is the length of the semimajor 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 ${\rm{\Phi }}(x={r}_{\min },y=0,z=0)$ and ${\rm{\Phi }}(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}_{\mathrm{interp},\min }=0.01\,{r}_{\min }$ and ${r}_{\mathrm{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}_{\mathrm{interp},\max }$. The inner boundary of ${r}_{\mathrm{interp},\min }=0.01\,{r}_{\min }$, however, can be problematic because centrophilic box orbits can pass well within $0.01\,{r}_{\min }$. The DOP853 Runge–Kutta integrator in the TriOS code uses adaptive time steps, tuning them to minimize errors in the position and velocity between time steps. In this scheme many acceleration evaluations are required in regions of the trajectory where the time step 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\buildrel{\prime\prime}\over{.} 1$ and ${r}_{\mathrm{interp},\min }\,=0.01\,{r}_{\min }=0\buildrel{\prime\prime}\over{.} 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 more efficient use of the acceleration table, we choose to decouple ${r}_{\mathrm{interp},\min }$ and ${r}_{\mathrm{interp},\max }$ from r_min and r_max, which are used to determine the range of orbital energy sampling. When ${r}_{\mathrm{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 speedup 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}_{\mathrm{interp},\min }$. For the example shown in Figure 3, choosing ${r}_{\mathrm{interp},\min }\sim 0\buildrel{\prime\prime}\over{.} 0001$ would reduce the orbit integration time by a factor of ≳2 compared with the original setting of ${r}_{\mathrm{interp},\min }=0.01\,{r}_{\min }=0\buildrel{\prime\prime}\over{.} 001$. Because energy conservation prevents orbits from passing outside r_max, setting ${r}_{\mathrm{interp},\max }$ to be slightly larger than r_max minimizes the integration time while maximizing the interpolation accuracy.

Zoom In Zoom Out Reset image size

Because we do not 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.

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}_{{{\rm{I}}}_{2}}=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}_{{{\rm{I}}}_{3}}=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}_{{{\rm{I}}}_{2}},{N}_{{{\rm{I}}}_{3}},{N}_{\mathrm{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.

Zoom In Zoom Out Reset image size

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 a line at an angle

$$\begin{eqnarray}&&\eta ={\tan }^{-1}\sqrt{\displaystyle \frac{T}{1-T}}.\end{eqnarray} \tag{ 10 }$$

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 long-axis tubes (for a prolate axisymmetric potential).

When orbits are well sampled, model properties such as the goodness of fit (χ²) 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, χ² varies periodically with T with a frequency matching the spacing between dithered orbits, $(\pi /2)/{N}_{{{\rm{I}}}_{2}}$, 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 artificial oscillations in χ² with that same period. Because 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 χ² versus T by increasing ${N}_{{I}_{2}}$, 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}_{{I}_{2}}$ from the default value of 9 to 15 and beyond removes the oscillations and also yields convergent results. We choose ${N}_{{I}_{2}}=18$ for the x − z start space.

We do not find similar issues for the other start space. Nonetheless, we increase ${N}_{{I}_{2}}$ 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}_{{I}_{2}},{N}_{{I}_{3}},{N}_{\mathrm{dither}})=(40,18,9,3)$ for both start spaces. This results in a total 40 × 18 ×9 × 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).

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 prespecified 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 of mass bins would have difficulty satisfying the 1% precision and/or contribute disproportionately high values to the total χ² 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 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 χ² landscape, as illustrated in the right panel of Figure 5. The more oblate (T ≲ 0.1) and round ($q\gtrsim 0.9q^{\prime}$) systems are more prone to this issue.

Zoom In Zoom Out Reset image size

We find that this mass misfitting problem can be easily resolved by using mass bins linearly spaced in $\cos (\theta )$ 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 bin-to-bin variations representing the galaxy's intrinsic deviation from spherical symmetry. Correspondingly, the prespecified 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 χ² returned by the code includes contributions from fits to the masses as well as kinematics. In this case, models with significant mass misfits due to uneven binning schemes would have disproportionately larger χ² 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, 2020) and wide-field McDonald Mitchell IFS data (Veale et al. 2017b, 2017a, 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 5) 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 surface-brightness 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₀ = 70 kms⁻¹ Mpc⁻¹.

We conduct the search for the best-fit galaxy shape in the new triaxial parameters $(T,{T}_{\mathrm{maj}},{T}_{\min })$ introduced in Section 3. The dark matter halo is modeled as a logarithmic potential. We parameterize it through its mass within 15 kpc, M₁₅, 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₁₅, we sample the 6D parameter space of galaxy models.

We determine the best-fit parameters by minimizing a χ² 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 is used to calculate the χ² 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. However, 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 six 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₁₅ ∈ [3.5, 10.5] × 10¹¹ M_⊙, and M_BH ∈ [1, 5] × 10⁹ M_⊙, and the full range between 0 and 1 for $(T,\sqrt{{T}_{\mathrm{maj}}},\sqrt{{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 six-dimensional χ² surface is computed from the previously run models. We use this interpolation to estimate the χ² for O(10⁴) points chosen using the LHS scheme described above in the original volume and select points where the estimated χ² is within Δχ² = 20.06 (3σ for six parameters) of the estimated global minimum. To avoid premature optimization we perform this routine 10 times where random subsets of half of all previously run models are used to build the interpolation function. With this scheme we select roughly 1000 model points that are expected to lie near the global χ² 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 χ² 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,\sqrt{{T}_{\mathrm{maj}}},\sqrt{{T}_{\min }})$ to $(T,{T}_{\mathrm{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}_{\mathrm{maj}},{T}_{\min })$ do not have the strong degeneracy apparent in (p, q, u).

Zoom In Zoom Out Reset image size

The standard values of Δχ² = 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 two degrees of freedom, giving Δχ² ≈ 2.3, 6.2, 11.8. This is different from most previous work using the van den Bosch et al. (2008) code, where typically ${\rm{\Delta }}{\chi }^{2}=\sqrt{2{N}_{\mathrm{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 χ² values for each model and is much larger than our values. However, while this is true when the input data are varied according to their noise level as discussed in Vasiliev & Valluri (2020), the noise level in the χ² values between models are significantly smaller when the input data are fixed.

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 modeling 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

Table 2. Best-fit Triaxial Model Parameters for NGC 1453 from the 6D Likelihood Landscape in Figure 6

Parameter	Value
MBH(109 M⊙)	2.9 ± 0.4
M*/LF110W(M⊙/L⊙)	2.02 ± 0.07
M15(1011 M⊙)	${7.0}_{-0.5}^{+0.6}$
T	0.33 ± 0.06
Tmaj	0.89 ± 0.03
${T}_{\min }$	0.12 ± 0.03
u	${0.941}_{-0.013}^{+0.012}$
p	${0.933}_{-0.015}^{+0.014}$
q	0.779 ± 0.012
θ (°)	73 ± 3
ϕ (°)	19 ± 3
ψ (°)	${92.7}_{-0.8}^{+0.7}$

Note. For each parameter, all other dimensions have been marginalized over.

Download table as:  ASCIITypeset image

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

Zoom In Zoom Out Reset image size

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

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 χ² 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 nonaxisymmetric models are favored, we examine the 2D maps of V and the lowest three 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 antisymmetric for points mirrored across the projected minor axis. Any observed systemic deviation from bisymmetry would then indicate triaxiality.

Zoom In Zoom Out Reset image size

For this reason, we decompose each GMOS moment map into a bisymmetrized component (second column) and a nonbisymmetrized 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 285. The residual pattern persists and cannot 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 nonbisymmetric 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 to those in the data (third column). In comparison, the best-fit triaxial model is able to fit these nonbisymmetric features to a large extent, producing essentially random residuals (fifth column).

Figure 8 indicates that the preference for triaxiality is driven by the nonaxisymmetric 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 being 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⁹ 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.