The MASSIVE Survey. XVII. A Triaxial Orbit-based Determination of the Black Hole Mass and Intrinsic Shape of Elliptical Galaxy NGC 2693

We present a stellar dynamical mass measurement of a newly detected supermassive black hole ( SMBH ) at the center of the fast-rotating, massive elliptical galaxy NGC 2693 as part of the MASSIVE survey. We combine high signal-to-noise ratio integral ﬁ eld spectroscopy ( IFS ) from the Gemini Multi-Object Spectrograph with wide-ﬁ eld data from the Mitchell Spectrograph at McDonald Observatory to extract and model stellar kinematics of NGC 2693 from the central ∼ 150 pc out to ∼ 2.5 effective radii. Observations from Hubble Space Telescope WFC3 are used to determine the stellar light distribution. We perform fully triaxial Schwarzschild orbit modeling using the latest TriOS code and a Bayesian search in 6D galaxy model parameter space to determine NGC 2693 ʼ s SMBH mass ( M BH ) , stellar mass-to-light ratio, dark matter content, and intrinsic shape. We ﬁ nd


INTRODUCTION
The most massive SMBHs in the local universe have been found at the centers of some of the most massive nearby elliptical galaxies.At stellar masses M * 10 11.5 M targeted by the volume-limited MASSIVE galaxy survey (Ma et al. 2014), a majority of these massive elliptical galaxies exhibit slow or no detectable rotation (Veale et al. 2017a,b;Ene et al. 2018).When selected by environments, massive galaxies in the cores of galaxy clusters are also predominantly slow-or nonrotators (e.g., Brough et al. 2017;Loubser et al. 2018;Krajnović et al. 2018a;Graham et al. 2018).In comparison, the ATLAS 3D project surveyed early-type galaxies at lower masses (10 10 M M 10 11.5 M ) and found most to be fast rotators (Emsellem et al. 2011).
When the kinematic axis of a massive elliptical galaxy can be identified, it is often misaligned from the photometric major axis (e.g., Ene et al. 2018).Detailed IFS kinematic maps also show intricate local twists, and the central and main-body kinematic axes within a galaxy are not always aligned (Ene et al. 2020;Krajnović et al. 2020).All these features are strong indications that local massive elliptical galaxies are triaxial in intrin-sic shape, and not axisymmetric as is often assumed in prior dynamical modeling studies of early-type galaxies, which would only allow rotation about the minor axis (Binney 1985).
The MASSIVE galaxy survey (Ma et al. 2014) is designed to study all the major dynamical components -SMBH, stars, and dark matter halo -in the most massive ∼ 100 M * 10 11.5 M early type galaxies (ETGs) in the local volume (out to ∼ 100 Mpc).For a subset of 20 MASSIVE galaxies, we have completed the stellar kinematic measurements from our IFS observations that cover both the galaxies' central regions with high spatial resolution, and wide fields out to at least one effective radius.In this paper, we focus on NGC 2693, a galaxy with one of the largest ratios of V /σ, where the stellar rotation is V ∼ 160 km s −1 and the velocity dispersion is σ ∼ 300 km s −1 (Veale et al. 2017a;Ene et al. 2020).The only other MASSIVE galaxy that exhibits similarly fast and regular rotation is NGC 1453.NGC 1453 is recently studied both in the axisymmetric limit (Liepold et al. 2020) and with fully triaxial orbit modeling (Quenneville et al. 2022).The triaxial models are found to better recover the input kinematics while also fitting the non-axisymmetric features present in NGC 1453.Given the similarities between NGC 2693 and NGC 1453, we turn to NGC 2693 in this paper.
We use the orbit superposition method to obtain dynamical mass measurements of the components of NGC 2693 with observations taken as part of the MAS-SIVE Survey.We use our latest version of the TriOS code (Quenneville et al. 2021;Quenneville et al. 2022) based on the code by van den Bosch et al. (2008) 1 to perform a full triaxial modeling of the stellar orbits in NGC 2693 and to simultaneously constrain the galaxy's intrinsic shape, M BH , and other mass parameters.
In Section 2, we describe the photometric observations used to model the deprojected stellar mass distribution of NGC 2693, as well as the spectroscopic observations from GMOS (Hook et al. 2004) covering the central kpc and wide-field observations from the McDonald Mitchell IFS (Hill et al. 2008).In Section 3, we describe the triaxial modeling and phase space sampling in our Schwarzschild orbit models.In Section 4, we discuss our search for the best-fit triaxial galaxy model, marginalization scheme for extracting best-fit parameters, and resulting best-fit dynamical model.In Section 5, we compare the triaxial model to axisymmetric Schwarzschild orbit models and Jeans modeling of NGC 2693.
We adopt a distance to NGC 2693 of 71.0 Mpc from the MASSIVE-WFC3 project (Goullaud et al. 2018) using the surface-brightness fluctuation technique (Jensen et al. 2021;Blakeslee et al. 2021).At this distance, 1 is 354 pc, assuming a flat ΛCDM cosmology with a matter density of Ω m = 0.315 and a Hubble parameter of H 0 = 70 km s −1 Mpc −1 .

PHOTOMETRIC AND SPECTROSCOPIC OBSERVATIONS
NGC 2693 is a relatively isolated galaxy, being the only identified member of its galaxy group in the 2MASS "high-density contrast" group catalog (Crook et al. 2007).We obtain photometric observations of NGC 2693 in the F110W filter of HST and spectroscopic observations using GMOS in IFS mode on the 8.1 m Gemini North telescope and the Mitchell IFS on the 2.7m Harlan J. Smith Telescope at McDonald Observatory.In this section, we describe these observations, the data reduction process, modeling the surface brightness profile, and the extraction of the stellar kinematics of NGC 2693.
Table 1.Best-fit MGE parameters to the NGC 2693 HST WFC3 IR photometry.Each Gaussian component (k) is parametrized by a central surface density I k = L k /2πσ 2 k q k , dispersion σ k , and axis ratio q k , where primed variables are projected quantities.The central surface densities have been corrected for a galactic extinction of 0.017 mag and assume an absolute (Vega) magnitude of 3.89 for the Sun.The components all have a PA of 167.9 • east of north.We model the spatial distribution of stars in NGC 2693 with observations from the IR channel of HST WFC3 in the F110W filter (Figure 1).Observations (GO-14219, P.I.J. Blakeslee) were taken over a single orbit and have a total exposure time of 2695 seconds.This orbit was divided into five dithered exposures with a sub-pixel dither pattern to improve measurements of the point-spread function (PSF).We reduce the images using STScI's standard reduction pipeline, a specialized Python program2 to correct for variable background levels, and the Astrodrizzle package (Gonzaga et al. 2012).We additionally perform background subtraction, removing a neighboring galaxy located 55 to the south of NGC 2693, and construct a mask to exclude foreground stars, other galaxies, and detector artifacts.The final F110W image has a pixel scale of 0. 1.For details on the photometric data reduction, see Goullaud et al. (2018).
The HST observations of NGC 2693 show slightly boxy isophotes near the center of the galaxy, which become disky at radii larger than ∼ 5 .There is a small compact dust disk extending 1. 5 in radius from the center.The galaxy's luminosity-weighted ellipticity is nearly constant with radius (beyond the region of the central dust disk) with a mean value L = 0.27 ± 0.002 (Goullaud et al. 2018).Below we parametrize the surface brightness of NGC 2693 as a sum of 2D Gaussians with a common center and PA.
We run the Cappellari (2002) Multi-Gaussian Expansion (MGE) code with regularization to avoid flattened components that artificially restrict the range of inclination angles that can be used during the dynamical modeling.We then tweak the MGE solution using Galfit (Peng et al. 2002).We set a lower boundary on each component's projected axis ratio, q k , of 0.65, which was determined from the previous regularized MGE run.The WFC3 PSF is accounted for by an empirical PSF constructed from extracting, summing, and re-normalizing 11 bright stars within the field of view.We apply the mask of other objects in the field plus a mask for the central dust disk.Since no high-resolution multi-band imaging for NGC 2693 was available, we construct a dust mask by eye initially, conservatively flagging only the most obviously affected pixels, and then we adopt an iterative approach.After the Galfit run converged, we examine the residual image and extend the dust mask to neighboring pixels with residuals above a selected threshold.We continually repeat the process, each time modifying the threshold, until achieving residuals at the ∼ 5% level at the nucleus.
Our best-fit MGE is composed of 10 Gaussians with the central surface brightness I k , projected dispersion σ k , and q k given in Table 1.The model is a good description of the data, as seen in Figure 1, with residuals below ∼ 6% out to a radius of ∼ 70 .
All of the MGE components have the same PA of 167.9 • east of north, but we run Galfit again allowing for the PA to vary between components.The initial guesses for the parameters are set to the best-fit MGE from Table 1, and we use the same empirical PSF and mask.When allowing for PA twists, our best-fit MGE consists of three circular Gaussians at small radii.The next five Gaussians have PA twists of < 3 • relative to a PA of 167.9 • and the outer three components have smaller PAs of ∼ 155 • .While we subtracted the companion galaxy and masked the remaining residuals, the outermost region of NGC 2693 likely remains contaminated by the companion galaxy causing the smaller PAs for the largest three MGE components.Nevertheless, these two MGEs differed by less than 5% at all radii for which we have kinematic data.Beyond 50 , the relative difference approaches ∼ 10% primarily due to the companion.The negligible impact of allowing for a PA twist does not seem to improve our fit but rather fits the remaining contamination of the companion galaxy, and thus we adopt the MGE model with a spatially constant PA.
We further test the companion galaxy's impact on the measured photometric position angle by fitting subsections of the full HST WFC3 image.We perform two additional fits, one to the central 20 × 20 region and another to the central 40 ×40 region as opposed to the full 141 × 125 image.We use our fiducial MGE presented in Table 1 as initial guesses for the ten Gaussians, and we require the PA and center for all ten components to be the same.In both cases, the preferred photometric position angle changes negligibly compared to the fiducal case of 167.9 • E of N. When fitting only to the central 20 × 20 region of the image, the preferred photometric PA is 168.6 • E of N, and when fitting to the 40 × 40 region, the PA of the Gaussians is 168.2 • E of N.

Central GMOS kinematics
We observed the central region of NGC 2693 in the 2016B semester with the two-slit mode of GMOS, providing a 5 × 7 (∼1.7 kpc × 2.4 kpc) field-of-view composed of 1000 hexagonal lenslets, for a projected diameter of 0.2 per lenslet.In total, six science exposures were taken, each of 1200 seconds exposure time.The median seeing was 0.6 FWHM.We used the R400-G5305 grating with the CaT filter for clean coverage from 7800 − 9330 Å.A 5 × 3.5 region of the sky, offset 1 from NGC 2693's central region was simultaneously observed.The spectral resolution was determined from arc lamp lines for each lenslet, with a median value of 2.3 Å FWHM.
The GMOS lenslets are binned to achieve a target S/N of 125 using the Voronoi-binning procedure (Cappellari & Copin 2003).This procedure results in a total of 60 bins.The spectra are co-added from individual lenslets within a single Voronoi bin as described in Ene et al. (2019).Example spectra for three bins at increasing radii are shown in Figure 2 (black curves).
We use the CaII triplet absorption features over a rest wavelength of 8420-8770 Å to extract the stellar line-ofsight velocity distribution (LOSVD) for each bin with the penalized pixel-fitting (pPXF) method of Cappellari (2017).We choose to decompose each LOSVD into a Gauss-Hermite series up to order n = 8 as is done in Liepold et al. (2020).An additive polynomial of degree zero (a constant) and a multiplicative polynomial of degree three are used to model the stellar continuum for the spectra.We compare two sets of stellar templates chosen from the MILES Calcium Triplet (CaT) Library of 706 stars (Cenarro et al. 2001): the first set is limited to the 15 Each panel shows a different moment of the Gauss-Hermite expansion of the line-of-site velocity distribution: the top-left two panels are the velocity V and velocity dispersion σ, with the higher-order hi moments characterizing deviations from Gaussianity.The +x axis of the galaxy is located 167 degrees East of North (North is up and East is to the left).The velocity map shows a prominent rotation pattern with maximal velocities of |V | ∼ 160 km s −1 ; the σ maps shows a central peak.The black lines are surface brightness contours from the best-fitting MGE model in Table 1 and Figure 1.stars from Table 2 of Barth et al. (2002); the second set includes all 706 stars in the library.The resulting kinematics are consistent within measurement errors, and we choose to use the former set of kinematics for dynamical modeling discussed later in this paper.These template spectra cover the wavelength range of 8347-9020 A with a 1.5 A spectral resolution FWHM.
We follow a nearly identical kinematic extraction procedure to that of Liepold et al. (2020).Three example stellar templates broadened by the best-fit LOSVDs are shown in Figure 2 (blue curves).The errors on the kinematic moments are determined with bootstrap methods as described in Section 4 of Ene et al. (2019).Figure 3 shows the kinematic maps of the eight Gauss-Hermite velocity moments for all 60 GMOS bins; Figure 4 shows the corresponding radial profiles of the moments (blue filled circles).The velocity map shows regular rotation with the maximum velocity reaching |V | ∼ 160 km s −1 ; the σ map shows a central peak with an amplitude of ∼ 320 km s −1 .The mean errors on V and σ are 4.3 km s −1 and 5.0 km s −1 , respectively.The errors on higher order moments (h 3 through h 8 ) are similar in amplitude, ranging from 0.010 to 0.013.

Wide-field Mitchell kinematics
NGC 2693 is one of ∼ 100 MASSIVE galaxies observed using the Mitchell IFS.Three dither positions were used, and during each dither, we interleaved two 20 minute science frames with one 10 minute sky frame, for a total exposure time of 2 hours on-source.The Mitchell IFS consists of 246 fibers covering a 107 × 107 field of view.The observations covered a spectral range of 3650- ) and shape parameters (T, Tmaj, Tmin) = (0.39, 0.09, 0.17).The spatial bins have been unfolded such that the bins whose centers lie between −90 • and +90 • of the photometric PA are plotted with positive radius and others with negative radius.Note that the axes are on a linear scale between −1 and 1 and logarithmically scaled otherwise.
5850 A that include the Ca HK region, the G−band region, Hβ, Mgb, and several Fe features.
Each fiber in the central region of NGC 2693 yields a spectrum of S/N 50, while the outer fibers are binned to achieve a S/N ≥ 20.We obtain a total of 33 bins but drop the outermost 4 bins at ∼ 46 due to low S/N .We model the stellar LOSVD in a similar way to the GMOS data, fitting up to order n = 6 due to the lower S/N .We use the MILES library of 985 stellar spectra (Sánchez-Blázquez et al. 2006;Falcón-Barroso et al. 2011) as the templates.Details of the Mitchell data reduction and kinematic measurements are described in Veale et al. (2017a) and Veale et al. (2017b).
The radial profiles of the Mitchell kinematic moments are shown in Figure 4 (pink filled circles).The innermost few Mitchell fibers (each with a 4.1 diameter) overlap the FOV of GMOS.Reassuringly, the Mitchell kinematics in this region are in good agreement with the GMOS values averaged over the GMOS FOV.The mean errors on the Mitchell kinematics are about two times larger than the GMOS kinematics.The mean errors on V and σ are 8.95 km s −1 and 12.2 km s −1 , respectively, with errors on higher order moments (h 3 through h 6 ) ranging from 0.030 to 0.039.

The TriOS Code and Input Kinematics
We use the TriOS code (Quenneville et al. 2021;Quenneville et al. 2022) to perform orbit modeling of NGC 2693.The galaxy is modeled as a stationary, triaxial gravitational potential composed of a SMBH, a stellar component specified by the deprojected MGE, and a dark matter halo.The six parameter galaxy model and the method used to search this parameter space is described in detail in Section 4.1.
For a collection of orbits spanning the phase space, we integrate each orbit for 2000 (loop orbits) or 200 (box orbits) dynamical times.At a large number of steps along the trajectory, we project the orbit onto the sky, keeping track of the projected position and line-of-sight kinematics.For each dataset (GMOS and Mitchell), the projected position is convolved with the respective instrumental PSF, and steps in the trajectory are binned according to the apertures described in Sections 2.2 and 2.3.This produces (i) a measure of the fraction of that orbit's time spent in each aperture and (ii) the LOSVD associated with the orbits within each kinematic bin.As described in Section 3.2, a superposition of these orbital contributions is found that best reproduces the observed kinematic LOSVD.We repeat this process for many gravitational potentials to find the galaxy parametrization that most closely reproduces the observed kinematics.
We fit for eight moments of the GMOS kinematics and six moments of the Mitchell kinematics presented in Section 2. We constrain additional moments up to and including h 12 to be 0.0 ± δ, where δ represents typical errors in the highest odd and even moments of each kinematic dataset.Our previous tests have showed that constraining only the lowest four moments and leaving the higher-order moments free in the orbit model could result in spurious behavior in the higher-order moments and the resulting LOSVDs (Figs. 10 and 11 of Liepold et al. 2020).We test the effects of constraining up to h 16 , but moments h 13 to h 16 are sufficiently close to 0 when constraining the first twelve moments that the added constraint on the last four moments does not change our fits.Thus, we opt to leave moments higher than h 13 unconstrained.We note that leaving moments h 9 to h 12 unconstrained can shift the inferred best fit parameters by ∼ 10%.
Throughout the analysis, we model the GMOS and Mitchell PSFs as single, circularly symmetric Gaussians with FWHM of 0.56 and 1.2 respectively.To keep the potential non-singular at the origin, we use a Plummerstyle potential for the central black hole with a force softening length of 3 × 10 −4 arcsec, which is roughly three orders of magnitude smaller than the central-most GMOS bins.

Orbital Phase-space Sampling and Orbital Weight Optimization
Insufficient orbit sampling can bias the inferred mass and shape parameters, so particular care is needed in choosing initial conditions for the orbits.In the TriOS code, the orbits are initialized in two different spaces, referred to as "start spaces" (Schwarzschild 1993;van den Bosch et al. 2008).A Cartesian coordinate system centered on the galaxy's nucleus is used.The x, y, and z axes are parallel to the intrinsic major axis a, intrinsic intermediate axis b, and intrinsic minor axis c, respectively.
The first start space, called the x-z start space, launches loop orbits from the (x, z) plane and a velocity in the +y direction at N E = 40 different energies (implicitly sampled over radii).The (x, z) positions are determined by dividing the start space into N I2 rays spanning polar angles from 0 to π/2 in the x-z plane; along each ray, we space N I3 orbits.The code allows for additional dithering of orbits, where we group together N dither = 3 adjacent initial conditions for denser and smoother phase-space sampling.Dithering orbits increases the sampling density by a factor of N 3 dither since dithering is performed in all three dimensions.Orbits launched from the same initial position but with velocity in the −y direction are valid orbits.To include these orbits in our modelling, we invert the LOSVD from each orbit in the (x, z) start space and store the resulting orbits in a second 'retrograde' library.
Similar to NGC 1453 (Section 4.3 of Quenneville et al. 2022), we find spurious oscillations in the goodness-of-fit χ 2 landscape while varying T for NGC 2693 when using N I2 = 9 and N I3 = 9 in the (x, z) start space.The spacing between these oscillations matches the spacing between dithered initial conditions, resulting in periodic local minima and thus biased results for the intrinsic galaxy shape.We eliminate these unwanted oscillations by increasing the angular sampling N I2 of orbits in the (x, z) start space from N I2 = 9 to N I2 = 15.The total number of orbits used in our loop library is therefore , where the factor of 2 accounts for the time reversed copy of each orbit.
A second start space, called the stationary start space, is used to generate a library of box orbits in a triaxial system.In this start space, orbits with a given energy E are launched at rest with v x = v y = v z = 0 from starting positions, where the gravitational potential is Φ(R, θ, φ) = E, and θ and φ are the polar and azimuthal angles in the usual spherical coordinate system.Both θ and φ are sampled uniformly between 0 and π/2 at N θ = 9 and N φ = 9 locations.We find no oscillatory behavior in the resulting χ 2 for this start space, so we choose (N θ , N φ , N dither ) = (9, 9, 3).The total number of orbits in our box library is thus 3 3 × 40 × 9 × 9 = 87, 480.
With three integrated orbit libraries consisting of 291, 600 + 87, 480 = 379, 080 orbits for a given galaxy model, we solve for the linear combination of orbital weights that best fits the observed kinematics and surface brightness.As noted, we use a dithering factor of 3, meaning that 3 3 = 27 neighboring orbits are bundled while solving for the orbital weights.Accordingly, there are 379, 080/3 3 = 14, 040 independent weights in each model.We compute these weights to minimize the χ 2 associated with the kinematics using non-negative least squares (Lawson & Hanson 1995), under the constraint that both the projected mass within each aperture and the 3D mass distribution in coarse bins are consistent within 1% of the MGE.We do this for all kinematic bins simultaneously to build a model LOSVD in each bin for each galaxy model.
We note that the choice of 1% mass constraint above is based on the range of values commonly adopted in earlier orbit modeling work (e.g., 2% in van den Bosch & de Zeeuw 2010; 1% in Walsh et al. 2017).Our results do not depend on the exact value used: relaxing the 3D mass constraint from 1% to 10% changes our best-fit M BH reported below by only a few %, which is well within the 1σ uncertainty.Additionally, the typical χ 2 associated with the mass constraint is only a small fraction (∼ 0.5%) of the size of the χ 2 from the stellar kinematics.We obtain identical M BH regardless of whether we use χ 2 from kinematics alone or include χ 2 associated with the mass constraint.

Galaxy Model
We use a six-parameter model to describe the triaxial potential of NGC 2693.Three parameters are for the mass components: black hole mass M BH , stellar massto-light ratio M * /L F 110W (hereafter M * /L), and enclosed dark matter mass at 15 kpc, M 15 .A logarithmic dark matter halo is assumed, where the enclosed dark matter mass at radius r is where R c is the scale radius and V c is the circular velocity of the halo.We fit for the enclosed dark matter mass at 15 kpc, which is approximately the outer edge of our outermost mass bin.
The other three parameters of our galaxy model encode the intrinsic shape of NGC 2693.Previous dynamical studies used angles (θ, φ, ψ), or axis ratios (u, p, q) = (a /a, b/a, c/a), to relate the projected and intrinsic shapes of a triaxial galaxy.Here, u is a compression factor relating the intrinsic (unprimed) and projected (primed) length scales; p is a ratio between the intrinsic intermediate axis b and the intrinsic major axis a (with projected major axis a ); and q is the ratio between the intrinsic minor axis c and the intrinsic major axis.Additionally, q = b /a describes the projected flattening of the MGE component.In general these axis ratios are different for each MGE component as those components have different projected flattenings.
We instead adopt the new parameters (T, T maj , T min ) introduced in Quenneville et al. ( 2022): where T is the commonly used triaxiality parameter, T maj parametrizes the length of the projected major axis relative to the minimum value b and maximum value a, and T min parametrizes the length of the projected minor axis relative to the minimum c and maximum b.The three parameters T, T maj and T min form a convenient unit cube, each with an allowed range of 0 to 1.They have a number of additional desirable properties when compared to the axis ratios (u, p, q) or angles (θ, φ, ψ); see Section 3.4 of Quenneville et al. (2022).Note that while these shape parameters are expressed in terms of the axis ratios in Equation (2), they are constant for different MGE components when the PAs of the MGE components are identical.As shown in Equations ( 8) and (A2) in Quenneville et al. (2022), the angles (θ, φ, ψ) can be directly computed from a (T, T maj , T min ) triplet and vice versa.

Latin Hypercube Sampling and Bayesian Search
As in Quenneville et al. (2022), we use the grid-free Latin hypercube sampling method (McKay et al. 1979) to search the 6D model parameter space: M BH , M * /L, M 15 , T , T maj , and T min .Latin hypercube sampling is becoming increasingly common in computer-designed experiments due to the simplicity of the sampling algorithm and the desirable space-filling properties in high dimensional parameter spaces.We adopt the Latin hypercube scheme described in Jin et al. (2005) and implemented in the Python package SMT (Bouhlel et al. 2019).
In this method, we first divide each dimension of our search-space into N cells, where N is the number of galaxy models in a Latin hypercube batch.We then uniformly sample points in each dimension until each cell contains a point, uniformly filling the space.We note that we use the "center" dispersal criterion in SMT, which places new points in the center of each hypercube cell.The result is a set of model points spanning six dimensions more uniformly than a regular grid and allowing for a more representative sampling of the likelihood function as a function of the six model parameters.T = 0.39 +0.04 −0.04 0 .9 8 0 .9 9 1 .0 u = 0.991 +0.003 −0.004 u 0 .8 8 0 .9 0 0 .9 2 p 0 .6 9 0 .7 1 0 .7 3 0 .7 5 0 .9 8 0 .9 9 We marginalize over a smoothed 6D landscape generated with Gaussian process regression.The 1σ, 2σ, and 3σ contours are represented by the curves in red, green, and blue respectively, and as different shade of grey in the 1D panels.Above each 1D posterior distribution are the extracted best-fit values and 1σ confidence intervals.(Upper right) 1D and 2D marginalized posteriors in the axis-ratio space of (u, p, q).Below each 1D posterior are the best-fit values and 1σ confidence intervals.
When drawing points from intrinsic-shape space, we opt to sample uniformly in T maj and √ T min rather than in T maj and T min .For nearly axisymmetric galaxies, this sampling space results in fewer unrealistically flat models, increasing the efficiency of our parameter search.
After running preliminary searches over broad ranges of parameters, we choose the uniform prior ranges of M BH ∈ [0, 4] × 10 9 M , M * /L ∈ [2.0, 2.85] M /L , M 15 ∈ [1, 13] × 10 11 M , T ∈ [0.15, 0.65], T maj ∈ [0.0, 0.6], and T min ∈ [0.0, 0.4] for our parameters.After each hypercube realization of roughly 1000 galaxy mod-els, we model the χ 2 likelihood landscape as a function of the six parameters using Gaussian process regression.We construct posterior distributions of the space and estimate of the best-fit values for each parameter using dynamic nested sampling (Speagle 2020).To test for convergence in our model sampling, we perform jackknife resampling where the regression and parameter inference is repeatedly performed using subsets of the full suite of models.In total, we generate 8530 galaxy model points which, when jackknife-resampled, converge on the same best-fit galaxy parameters.
To map the high-likelihood region in fine detail and to serve as an additional test for convergence from our hypercube iterations, we perform one more independent check of our best-fit parameters.We again perform a 6-D gaussian process regression function to the χ 2 landscape produced by the 8530 models described above.We then sample 1000 additional model points with another hypercube, rejecting those which fall outside of the 3σ confidence volume, as estimated by that regression fitted to the previous 8530 model points.The rejection sampling procedure very efficiently populates the χ 2 minimum.Of the 1000 points in that sample, ∼ 9% lie within the 3σ confidence volume for 6 parameters (∆χ 2 ≈ 20.06), compared to just ∼ 4% from the uninformed 8530-point sample.
Both the uniform hypercube of 8530 models and the 1000 rejection-sampled hypercube models converge on the same best fit parameters.We include both sets here.The resulting 6D posterior distribution of our 9530 production run models is shown in Figure 5.We determine the best-fit value and uncertainties by fitting the χ 2 landscape with Gaussian process regression with a squared-exponential covariance kernel and sampling that landscape with the dynamic nested sampler described in Speagle (2020).A uniform prior is assumed for all parameters.The 1σ, 2σ, and 3σ confidence regions in 1D and 2D are computed from the posterior distribution marginalized over all other dimensions.These confidence levels correspond to ∆χ 2 = 1, 4, and 9 in 1-D and ∆χ 2 ≈ 2.3, 6.2, 11.8 in 2D.

Best-Fit Triaxial Model
The best-fit galaxy model is an excellent fit to the observed stellar kinematics, shown in Figure 4.The bestfit parameters from the 6D posterior distribution in Figure 5 are listed in Table 2.The total χ 2 of the best-fit triaxial model is 265.9, with 19.7 from the higher-order moments which are constrained to be 0.0 ± δ.As a test, we have repeated the regression using only four Gauss-Hermite moments measured from the spectra while setting h 5 and beyond to 0.0 ± δ described above.The best-fit M BH is within 0.5σ confidence interval of the fiducial model in Table 2, but the uncertainties on M BH in this case increase by ∼ 10%.The same trend is reported in Table 2 of Liepold et al. (2020).
As we will discuss further in Section 6.1, our inferred M BH = (1.7 ± 0.4) × 10 9 M for the SMBH in NGC 2693 is within the intrinsic scatter of the SMBH-galaxy scaling relations.For the intrinsic shape of NGC 2693, we can compare our inferred axis ratios (p, q) = (0.902 ± 0.009, 0.721 +0.011  −0.010 ) with the mean galaxy shapes for 49 slow-rotating galaxies in the MAS- Table 2. Summary of best-fit galaxy models for NGC 2693.
For each parameter, we marginalize over the other dimensions and report the 1σ uncertainties.The axisymmetric orbit models and JAM models have fixed inclination of 70 • .In orbit models, θ is the inclination angle in the oblate axisymmetric limit (ψ = 90 • , or equivalenly p = 1), with θ = 90 • being edge-on and θ = 0 • being face-on.† We measure βz in the orbit model as a function of radius, shown in the bottom panel of Figure 6.The best-fit JAM value of βz = 0.07±0.01 is consistent with the range of βz values measured from this best-fit model, with values ranging from βz = −0.27 at small radii to βz = 0.28 at large radii in both the triaxial and axisymmetric Schwarzschild models.
SIVE survey from Ene et al. (2018).In that work, based on the observed ellipticity and misalignment between the kinematic and photometric axes, the mean flattening of the galaxy sample was estimated statistically to be (p, q) = (0.88, 0.65), with 56% of galaxies having p > 0.9.The triaxial shape of NGC 2693 is therefore quite close to the mean values.NGC 1453, the other MASSIVE galaxy for which we have performed triaxial orbit modeling (Quenneville et al. 2022), has best-fit shape parameters of (p, q) = (0.93, 0.78), indicating a slightly less flattened shape than NGC 2693 and the mean MASSIVE galaxy.Additionally, the intrinsic shapes of NGC 2693 and NGC 1453 are consistent with the distribution of shapes of fast rotators found in the IllustrisTNG50 and IllustrisTNG100 simulations (Pulsoni et al. 2020), where the mean axis ratios are (p, q) ∼ (0.9, 0.52) and the dispersion is σ ∼ 0.15 for the most massive fast rotating elliptical galaxies.We note that by construction, the axis ratios (p, q, u) of each MGE component obey the relation 0 ≤ q ≤ uq ≤ p ≤ u ≤ 1 (see, e.g., Sec.2.1 of Quenneville et al. 2022).For the MGE of NGC 2693, the most flattened component has q = 0.684.For physically useful deprojections, we may expect q 0.2 (Binney & de Vaucouleurs 1981).The allowed ranges of (p, q, u) are therefore quite narrow, in particular for u, which is constrained to be between ∼ 0.9 and 1.The errors on these parameters in Table 2, while appearing small on absolute terms, are on the order of ∼ 5 − 10% of the allowed ranges.
We also highlight that the recovered viewing angle θ = 66 •+4 • −3 • , which corresponds to the galaxy's inclination in the oblate axisymmetric limit, is consistent with the inclination of the galaxy estimated from the nuclear dust disk at NGC 2693's center, which we measure to be i ≈ 70 We have also run two additional galaxy models using the best-fit parameters shown in Table 2, but with M BH = 0 and M BH twice the best-fit value to assess what features in the kinematics provide the black hole mass constraint.Both models are a worse fit to the kinematic data, giving a ∆χ 2 = 38.1 when there is no black hole present, and a ∆χ 2 of 17.9 when the black hole is twice as massive.As expected, the inner kinematic data provide significant constraints on M BH , with ∼ 30% and 50% of the additional ∆χ 2 coming from the inner ∼ 1 data for the two test cases, respectively.
We use the computed orbit libraries to calculate the orbital composition, as well as the radial velocity dispersion σ r and tangential velocity dispersion σ t ≡ σ 2 θ + σ 2 φ /2 of NGC 2693.We present the orbital fractions and two anisotropy parameters β and β z as a function of radius in Figure 6.Long-axis tubes and box orbits, both of which are only present in triaxial potentials, make up ∼ 40% of the orbits at small radii and 35% of the orbits at outer parts of the galaxy.Near the center of the galaxy, the orbits are mostly tangential with β < 0 but become radially anisotropic beyond ∼ 1 kpc.

AXISYMMETRIC DYNAMICAL MODELING OF NGC 2693
For a comparison study, we have performed axisymmetric modeling of NGC 2693 using both the orbit superposition method and Jeans modeling.We describe the results from each method below.

Schwarzschild Orbit Modeling in the Axiysmmetric Limit
We use the axi-symmeterized version of the TriOS code first described in Quenneville et al. (2021), with further improvements in mass binning and acceleration table discussed in Quenneville et al. (2022).Liepold et al. (2020) first applied this code to NGC 1453; here we use similar settings to achieve axisymmetry within the triaxial TriOS code.We ensure the low L z space is well sampled by tube orbits and do not include the box orbit library (which has L z = 0) explicitly.We set the viewing angle ψ sufficiently close to 90 • in the input parameter file, i.e., |ψ − 90 • | = 10 −9 , to ensure no long-axis tube orbits are present.For the remaining short-axis tube orbits, we enforce axisymmetry by making 40 copies of each orbit, each copy rotated successively by 2π/40 about the intrinsic minor axis of the galaxy (Section 3 of Quenneville et al. 2021).These three precautions are necessary to run the triaxial code in the axisymmetric limit and obtain robust parameter constraints.We choose (N I2 , N I3 , N dither ) = (9, 9, 3) for the phase space sampling, and include two copies of the integrated orbit library in our minimization.This gives a total number of 174, 960 orbits for our axisymmetric galaxy models.
We search for the best-fit galaxy model using the Latin hypercube scheme outlined in Section 3 over three dimensions: M BH , M * /L, and M 15 , with a fixed inclination angle i = 70 • , estimated from the nuclear dust disk.Our hypercube consists of 2000 galaxy models drawn from a range M BH = [0, 4] × 10 9 M , M * /L = [1.8,2.8] M /L , and M 15 = [1, 13] × 10 11 M .The best-fit axisymmetric model parameters are listed in Table 2.
Our best fit axisymmetric model of NGC 2693 prefers a ∼ 40% larger M BH compared to the triaxial case, though the recovered best-fit M * /L and dark matter halo are consistent with triaxial modeling at the 1σ level.By construction, axisymmetric models produce bisymmetric kinematic maps, meaning that the LOSVDs are symmetric about the photometric major axis and antisymmetric for points mirrored about the photometric minor axis.By contrast, LOSVDs in triaxial models are only point-symmetric about the origin.The apparent minor-axis rotation in "Non-bisymmetric Data Component" panel of Figure 7 therefore can be fit by triaxial models but not axisymmetric models.Our best-fit axisymmetric model fails to account for this component (lower middle panel), while our best-fit triaxial model captures well the full velocity features and produces featureless and nearly zero residuals (lower right panel).Our best-fit axisymmetric model fails to account for this component (middle, left panel), while our best-fit triaxial model captures well the full velocity features and produces featureless and nearly zero residuals (middle, right panel).
We have run an additional test to verify that the nonbisymmetric component of the kinematics show in Figure 7 is due to a physical non-alignment between the photometric and kinematic axes and can not be simply "rotated away".In this test, instead of using the best-fit photomertic PA 167.9 • given by the MGE (see Sec. 2.1), we inflate the photometric PA to be 175 • , a value that would minimize the magnitude of the non-bisymmetric component.Using this inflated PA, we then refit the MGE and re-compute the non-bisymmetric component of our input kinematics.Figure 8 compares the result for this test (bottom panel) with that of our fiducial PA (top panel).While the non-bisymmetric feature is indeed much reduced for PA = 175 • , the MGE isophotes for this PA provide a noticeably worse fit to the observed surface brightness profile at both at large and small radii.Lipka & Thomas (2021) recently argued that edgeon axisymmetric models have a larger model flexibility than face-on projections and thus can fit observational data better, biasing the recovered inclination to- wards i ∼ 90 • .The rationale is that in edge-on models, the prograde and retrograde orbits have opposite velocities along the line of sight and contribute uniquely to the model's LOSVDs, whereas in face-on models, the two sets of orbits have negligible line-of-sight velocities, making them virtually interchangeable and effectively reducing the number of unique orbits used in superposition and minimization routines.They reported a ∆χ 2 ∼ 30 bias, favoring edge-on inclinations.We have performed a parameter search including inclination as a fourth model parameter, sampling 1000 values from i = [68 • , 89 • ] in the hypercube.Our regression find a best fit value i = 87.6 •+0.9 • −1.8 • , with a ∆χ 2 between the lowest and highest inclinations of ∼ 25, slightly smaller than that reported in Lipka & Thomas (2021).Despite this preference for edge-on inclinations in the axisymmetric models, our best-fit M BH and M * /L barely change when we include inclination as a free parameter: M BH = (2.4 ± 0.5) × 10 9 M and M * /L = (2.23 ± 0.1)M /L .These are both consistent within the confidence intervals of the i = 70 • results, so our results are robust to choice of inclination angle.

Jeans Anisotropic Models
We further model the stellar kinematics of NGC 2693 as an axisymmetric system using Jeans anisotropic modeling (JAM ;Cappellari 2008aCappellari , 2020)).JAM solves the Beyond ∼ 1 , it is difficult to distinguish between the models.The best-fit model is a good match to the observations and the model without a SMBH underestimates Vrms at the nucleus.
Jeans equations assuming a velocity ellipsoid that is aligned with a cylindrical coordinate system (R, z, φ) or a spherical coordinate system (R, θ, φ).We adopt a cylindrically aligned velocity ellipsoid, which is flattened along the z-axis and is characterized by the anisotropy parameter β z = 1 − (σ z /σ R ) 2 , where σ z and σ R are the velocity dispersions parallel to the rotation axis and in the radial direction.JAM has the advantage of being computationally inexpensive and previous studies generally have found similar results between (axisymmetric) Schwarzschild models and JAM (e.g., Seth et al. 2014;Krajnović et al. 2018b;Thater et al. 2019).
In our model, the gravitational potential comes from the BH, stars, and dark matter.The galaxy's surface brightness from Table 1 is deprojected into a 3D stellar mass density given an inclination angle, i = 70 • , and M * /L, while the dark matter halo is parametrized by the logarithmic profile in Equation (1).Thus, the free parameters in our model are: M BH , M * /L, M 15 (the dark matter mass enclosed within 15 kpc), and β z .Given these parameters and the GMOS PSF, JAM predicts the second moment, which we compare to the observed V rms , with V rms = (V 2 + σ 2 ).We use the point-symmetrized V and σ from the GMOS and Mitchell observations, excluding the innermost Mitchell kinematics that spatially overlap with the GMOS kine-matics.We additionally exclude the outermost four Mitchell bins as was done in the Schwarzschild models in earlier sections.
The model parameters are optimized using Bayesian inference and the nested sampling code Dynesty (Speagle 2020), which estimates posteriors and evidences.We adopted a likelihood L ∝ exp(−χ 2 /2) where χ 2 = Σ i (D i − M i ) 2 /σ 2 i and D i and M i are the observed and model V rms , respectively, and σ i is the V rms uncertainty for each spatial bin.When running with Dynesty, we use 500 live points and stop the initial sampling stage once reaching a threshold of 0.05, which is the logratio between the current estimated Bayesian evidence and the remaining evidence.The batch sampling stage is stopped when the fractional error on the posterior reaches 0.02.We assume uniform priors, with all free parameters sampled linearly.The best-fit values and 1σ uncertainties are taken to be the median and 68% confidence intervals of the posterior distributions, respectively.
The results are shown in Table 2 and the comparison between the best-fit model and the observed V rms is given in Figure 9.The model reproduces the data well, with a reduced χ 2 of 1.07.Figure 9 also displays three models with M BH set to 0 M , 2.1 × 10 9 M (the 3σ lower bound), and 3.7 × 10 9 M (the 3σ upper bound), with M * /L, M 15 , and β z fixed to the values from the best-fit JAM model in Table 2.The M BH = 0 M case fails to match the kinematics in the inner region and further demonstrates the need for a black hole in the galaxy potential.The M BH and M * /L from the best-fit JAM model are consistent within 1σ of the axisymmetric Schwarzschild model results in Section 5.1, and β z falls within the range of values extracted from the best-fit axisymmetric Schwarzschild model.The M 15 value from JAM is lower than the Schwarzschild model result, but it remains consistent at the 2σ level.Liepold et al. (2020) found a similar result in the analysis of NGC 1453, with JAM favoring a value of M 15 half that inferred from the axisymmetric Schwarzschild models.
We complete additional JAM runs to test assumptions made during the modeling.In our fiducial model, we fix i = 70 • , which is the inclination angle inferred from the dust disk, but we also test allowing i to be a free parameter.We find a preference for i = 81 • , however all the angles for which the MGE could be deprojected fall within the 3σ uncertainties.In addition, we test using a spatially varying anisotropy, with a parameter (β z,in ) assigned to the MGE components with σ k < 3 and another parameter (β z,out ) attributed to the remaining MGE components.We then examine a case where β z,in corresponded to the MGE components with σ k < 6 .These choices were motivated by the previously run Schwarzschild models, which suggested a change in the anisotropy between a radius of ∼ 3 − 6 .Next, we fit to only the GMOS data, which extend to a radius of 3. 4, and we assume a spatially constant anisotropy.Finally, we test including the Mitchell kinematics from the outer spatial bins during the fit, adopting an modified MGE constructed using a dust mask with fewer central pixels flagged, and increasing the number of live points and applying different sampling thresholds in Dynesty.
Even when changing the model in these various ways, we nearly always find consistent results at the 1σ level with the fiducial model.The exception is when we fit to only the GMOS data with a spatially constant β z ; we find that M * /L is consistent within the 2σ uncertainties while the remaining parameters are in agreement at the 1σ level with the fiducial model.
The above work assumed a cylindrically aligned velocity ellipsoid, but we also examine using spherically aligned JAM.In this case, we find a large anisotropy, with β = 1−(σ θ /σ r ) 2 = 0.39 and an order-of-magnitude smaller M BH with 3σ uncertainties that extend to 0 M .We also fit spherically aligned JAM to only the GMOS kinematics and recover the same results.When repeating the run and fixing β = 0.0, we find that M BH is constrained with a best-fit value of 3.4 × 10 9 M .In this case, the M BH and remaining parameters are consistent with the fiducial (cylindrically aligned) JAM model given the 1σ uncertainties.
Despite the assumptions of cylindrically aligned JAM, the inferred M BH and M * /L match (at the 1σ level) the results from the more complex axisymmetric orbit model in Section 5.1 (Table 2).The enclosed dark matter mass from JAM is ∼ 40% lower than that from the axisymmetric orbit model, but it is within 2σ uncertainties of the orbit model.As Table 2 shows, the uncertainties in the best-fit parameters from JAM tend to be much smaller than those from the axisymmetric orbit model.We continue to see a shift in the M BH compared to the best-fit value from the triaxial Schwarzschild model, with the JAM value being 75% more massive than the value predicted from the triaxial modeling; see Section 6.2 for further discussion.2013) and ∼ 5% above the relation in Saglia et al. (2016); it is within the intrinsic scatter of both relations, with values of 0.38 dex and 0.417 dex, respectively.
For the M BH − M bulge relation, we use the total stellar mass of NGC 2693 from our best-fit triaxial model, M * = 7.2×10 11 M , as the bulge mass. 3The NGC 2693 M BH is ∼ 25% smaller than the value predicted by the mean M BH − M bulge relation of McConnell & Ma (2013) and ∼ 18% smaller than the value predicted by the Saglia et al. (2016) relation.Again, this SMBH is within the intrinsic scatter of both relations, with values of 0.34 and 0.535 dex, respectively.

Comparison of Triaxial and Axisymmetric Models
There are few studies that compare M BH determination from fully triaxial stellar dynamical models to axisymmetric models of the same galaxy.The best-fit M BH for both M32 (van den Bosch & de Zeeuw 2010) and NGC 1453 (Liepold et al. 2020;Quenneville et al. 2022) were unchanged when relaxing the assumption of axisymmetry, whereas M BH in NGC 3379 increased by a factor of ∼ 2 in the triaxial case (van den Bosch & de Zeeuw 2010).We note that the mass modeling performed for M32 and NGC 3379 did not simultaneously model the dark matter halo of the two galaxies.In comparison, we make no assumptions on the dark matter halo of NGC 2693, and instead constrain the dark matter mass at 15 kpc directly as was done for NGC 1453 (Quenneville et al. 2022).Furthermore, the triaxial code of van den Bosch et al. (2008) had an incorrect scheme for mirroring orbits, which we fixed in the TriOS code used for NGC 1453 and NGC 2693 here.
In the case of NGC 3998, Walsh et al. (2012) applied the triaxial code of van den Bosch et al. (2008) and considered different dark matter halos.The grid-based parameter search did not allow for simultaneously varying all parameters in their model.While NGC 3998 was not modelled in the axisymmetric regime, the gasdynamical measurement of M BH disagreed with the stellar-dynamical value by a factor of ∼ 4 (De Francesco, G. et al. 2006).
Recently, den Brok et al. ( 2021) applied the van den Bosch et al. (2008) code to the brightest cluster galaxy PGC 046832 to determine its intrinsic shape, central black hole mass, and orbital composition.The galaxy has a unique velocity map, exhibiting both a kinemat-ically decoupled core and dramatic twists in the velocity field, suggesting a non-axisymmetric intrinsic shape.Their triaxial models prefer prolate galaxy shapes in the inner 10 arcseconds of the galaxy, and oblate shapes beyond 10 arcseconds, though these models only provide an upper bound on the black hole mass of M BH 2 × 10 9 M .While this disagrees considerably with the results from their best-fit axisymmetric models, which prefer M BH ∼ 6 × 10 9 M , it remains to be seen if their triaxial result would change after the incorrect orbit mirroring in the van den Bosch et al. (2008) code and other issues discussed in Quenneville et al. (2022) are fixed.
In the case of NGC 2693, the best-fit orbit model in the axisymmetric limit and the best-fit JAM model favor M BH that is 40%-70% higher than the triaxial orbit model, but the difference is within ∼ 2σ confidence level (see Table 2).Similar comparison studies are needed from more galaxies to assess whether any systematic difference exists in M BH values determined from different methods.

SUMMARY
We have reported detection of a SMBH with M BH = (1.7 ± 0.4) × 10 9 M at the center of the massive, fastrotating galaxy NGC 2693 targeted by the MASSIVE survey.Using HST stellar light profiles and extensive IFS kinematic data covering a FOV from ∼ 150 pc to 15 kpc as constraints (Section 2), we have performed triaxial orbit modeling with the TriOS code to determine the galaxy's internal stellar orbit structure, M BH , M * /L, dark matter content, and intrinsic 3D shape (Section 3).We modeled the gravitational potential of NGC 2693 with 6 parameters and performed a 6D Bayesian search using Latin hypercube sampling of ∼ 10, 000 galaxy models to find the model that best matches our input data (Section 4).
Despite NGC 2693 exhibiting properties typically indicating an intrinsic axisymmetric shape, we find the best-fit model to be triaxial with T = 0.39 ± 0.04 and intrinsic axis ratios p = b/a = 0.902 ± 0.009 and q = c/a = 0.721 +0.011 −0.010 .We find that triaxial models are needed to account for non-axisymmetric features seen in the residuals of our accompanying axisymmetric models (Figure 7).When limiting ourselves to axisymmetry, we find 40% larger best-fit black hole mass of M BH = (2.4 ± 0.6) × 10 9 M from axisymmetric orbit modeling, and 75% larger best-fit black hole mass of M BH = (2.9 ± 0.3) × 10 9 M from JAM modeling (Section 5); both values are within ∼ 2σ confidence level of M BH determined from triaxial modeling (Table 2).
We have examined orbit flexibility in our galaxy models to assess possible effect of "generalized degrees of freedom" (Ye 1998;Spiegelhalter et al. 2002) on parameter determinations.Using a similar measure as Lipka & Thomas (2021) to estimate the effective number of parameters, we find that our models in the axisymmetric limit have a similar behavior as Lipka & Thomas (2021), in which edge-on orientations tend to have higher model flexibility (Section 5.1).Such varying model flexibility can be attributed to varying degeneracy between prograde and retrograde short-axis loop orbits as the line-of-sight approaches the symmetry axis.For triaxial models, however, we find the model flexibility to vary much less in the region around the best-fit models, and our best-fit triaxial shape parameters change by less than 1σ in a number of preliminary tests.It is possible that the additional presence of box and long-axis tube orbits in triaxial potentials has led to a weaker dependence of model flexibility on viewing angles.We will report the full results in a subsequent paper.
This paper adds to only a handful of other stellar dynamical modeling studies not limited to axisymmetric galaxy shapes (Section 6.2).Most of the remaining galaxies in the MASSIVE survey exhibit more prominent kinematic and photometric twists and less rotation compared to NGC 2693, further providing evidence that massive early-type galaxies have triaxial intrinsic shapes.More stellar dynamical measurements beyond the axisymmetric limit will inform whether the systematic differences in M BH seen for NGC 2693 in this paper is a common occurrence.

Figure 1 .
Figure 1.(Left) The F110W-band HST image of NGC 2693 used for our photometry.A companion galaxy which is masked from photometric analysis (see text) can be seen ∼ 50 south of NGC 2693.(Left inset) NGC 2693 has a dust disk extending approximately 1.5 (in radius) from the center.This feature is masked from our MGE fitting.(Right) Isophotes of the HST WFC3 IR image of NGC 2693 (black) and the best fit MGE model (red).The inset shows the central region containing the nuclear dust disk.The gray regions are masked when performing the fit as described in the text.

Figure 2 .
Figure 2. CaII triplet region of the GMOS IFS spectrum (black) of NGC 2693 for three example bins located at increasing distance from the nucleus: (top) central bin with S/N = 204, (middle) bin 1.50 from center with S/N = 144, and (bottom) bin 3.44 from center with S/N = 100.The spectrum is fit by broadening a series of stellar templates by the best-fit LOSVD (blue) over a wavelength range of 8420 − 8770 A. The grey shaded regions are excluded from the fit to account for improperly subtracted sky lines.The red points are the (best-fit minus observation) residuals offset by constants for clarity.The typical residual is ∼ 0.5%.

Figure 3 .
Figure 3. Spatial maps of the stellar kinematics extracted from the Gemini GMOS IFS spectra for 60 bins in the central 5 × 7 region of NGC 2693.Each panel shows a different moment of the Gauss-Hermite expansion of the line-of-site velocity distribution: the top-left two panels are the velocity V and velocity dispersion σ, with the higher-order hi moments characterizing deviations from Gaussianity.The +x axis of the galaxy is located 167 degrees East of North (North is up and East is to the left).The velocity map shows a prominent rotation pattern with maximal velocities of |V | ∼ 160 km s −1 ; the σ maps shows a central peak.The black lines are surface brightness contours from the best-fitting MGE model in Table1 and Figure 1.

Figure 4 .
Figure 4. Radial profiles of the velocity moments determined from the GMOS (blue circles) and Mitchell (pink circles) IFS observations.The observed moments are well matched by those predicted from our best-fit triaxial galaxy model (black open squares) with mass parameters (MBH, M * /L, M15) = (1.7 × 10 9 M , 2.35, 7.1 × 10 11 M ) and shape parameters (T, Tmaj, Tmin) = (0.39, 0.09, 0.17).The spatial bins have been unfolded such that the bins whose centers lie between −90 • and +90 • of the photometric PA are plotted with positive radius and others with negative radius.Note that the axes are on a linear scale between −1 and 1 and logarithmically scaled otherwise.

Figure 5 .
Figure 5. (Lower left) 1D and 2D marginalized posteriors for the triaxial orbit models of NGC 2693 described in the text.We marginalize over a smoothed 6D landscape generated with Gaussian process regression.The 1σ, 2σ, and 3σ contours are represented by the curves in red, green, and blue respectively, and as different shade of grey in the 1D panels.Above each 1D posterior distribution are the extracted best-fit values and 1σ confidence intervals.(Upper right) 1D and 2D marginalized posteriors in the axis-ratio space of (u, p, q).Below each 1D posterior are the best-fit values and 1σ confidence intervals.

Figure 6 .
Figure6.(Top) Fraction of orbital weights in each orbital family for the best-fit triaxial galaxy model.The orbital structure is dominated by short-axis tubes at all radii, with a non-zero fraction of the weights occupied by long-axis tubes and box orbits, both of which are present only in triaxial potentials.For axisymmetric models, short-axis tubes are the only allowed orbital family.(Middle) Velocity anisotropy β ≡ 1−σ 2 t /σ 2 r profile of the best-fit triaxial model (pink, solid line) and best-fit axisymmetric model (green, dashed line) of NGC 2693.Inner orbits are tangential out to ∼ 1 kpc and are increasingly radially anisotropic at larger radii in both the axisymmetric and triaxial cases.(Bottom) Anisotropy parameter βz = 1 − (σz/σR) 2 , where σz and σR are the velocity dispersions parallel to the rotation axis and in the radial direction, for the best-fit axisymmetric and triaxial models described in the text.

Figure 7 .
Figure 7. Observed velocity map of the central 5 × 7 of NGC 2693 (upper left), oriented such that the observed photometric major and minor axes are horizontal and vertical, respectively, where the best-fit MGE PA is 167.9 • .We decompose the map into a bisymmetrized component (upper right) and a non-bisymmetric component (lower row), where bisymmetry means symmetry for points mirrored about the photometric major axis and anti-symmetry for points mirrored across the photometric minor axis.The nonbisymmetric component (normalized by measurement uncertainty) shows a prominent apparent minor-axis rotation, a telltale sign of triaxiality.Since axisymmetric models can only produce bisymmetric velocity maps by construction, the residuals from our best-fit axisymmetric model (lower middle) shows a similar pattern as the non-bisymmetrized map.By contrast, our best-fit triaxial model (lower right) is able to reproduce the full observed velocity structure, and the residuals scatter randomly about 0.

Figure 8 .
Figure8.Illustration of the non-alignment between the photometric PA and kinematic features.The upper panel repeats the HST (black) and MGE model isophotes (red) in Figure1and the non-bisymmetric component of the GMOS data in Figure7.The best-fit MGE PA of the photometric major axis is 167.9 • in this fiducial case.In the lower panel, we inflate the MGE PA to 175 • and plot the resulting model fits and non-bisymmetric map assuiming this PA.While the non-bisymmetric velocity pattern is minimized, this inflated PA provides a poor fit to the observed surface brightness profile.

Figure 9 .
Figure 9. Line-of-sight rms velocity (Vrms) determined from the GMOS (blue dots) and Mitchell (pink dots) IFS kinematics.The best-fit JAM model is shown with black open squares.Also shown are three JAM models with MBH fixed to 0M (dotted line), 2.1×10 9 M (the 3σ lower bound; dashed line), and 3.7 × 10 9 M (the 3σ upper bound; dashed line).The three models extend over all radii, although only the model predictions within the central region are plotted.Beyond ∼ 1 , it is difficult to distinguish between the models.The best-fit model is a good match to the observations and the model without a SMBH underestimates Vrms at the nucleus.
Black Hole Scaling Relations To place the NGC 2693 SMBH on the M BH − σ relation, we use the luminosity-weighted velocity dispersion within R e , σ = 296 km s −1 , from Veale et al. (2017b) for NGC 2693.This measurement was obtained from the same Mitchell IFS data used in this paper.The mass of the NGC 2693 SMBH is within 15% of the value pre-dicted by the mean M BH − σ relation in McConnell & Ma (