A Slowly Precessing Disk in the Nucleus of M31 as the Feeding Mechanism for a Central Starburst

NIR observations of the nucleus of M31 were taken with the LGS AO system on the W. M. Keck II 10 m telescope (van Dam et al. 2006; Wizinowich et al. 2006). The laser was positioned on the nucleus to correct high-order atmospheric aberrations. Low-order aberrations were corrected using two different tip-tilt stars during different observing runs. During good seeing, a close and faint globular cluster was used for tip-tilt correction; it is located 35'' to the southwest from the nucleus with R = 16.2 at 00 42 42.203 + 41 15 46.01 (J2000). A second globular cluster, NB95, with R = 14.9 and located 53'' southeast from the nucleus at 00 42 47.973 + 41 15 37.07, was used as the tip-tilt reference when seeing was poor or variable or if clouds were present (Battistini et al. 1993; Galleti et al. 2007). These tip-tilt objects are located in a region of M31 with high surface brightness and a strong gradient in the unresolved galaxy light. Thus, in order to properly background-subtract the tip-tilt wavefront sensor, the sky background was measured at a manually selected sky position such that the light from M31's bulge had comparable intensity to that of the tip-tilt object positions.

2.1.1. Keck/OSIRIS

Spectroscopic measurements of the M31 nucleus were made using OSIRIS (Larkin et al. 2006), an IFU spectrograph, on 2008 October 21 (PI: M. Rich), 2010 August 15 (PI: J. Lu), and 2010 August 28–29 (PI: A. Ghez; Table 1). The individual OSIRIS data cubes cover 08 × 32 and were oriented at a P.A. of 56° along the major axis of the eccentric stellar disk. The field was sampled with a pixel scale of 005 pixel⁻¹, and each spatial pixel (spaxel) provides an independent spectrum across the K-broadband filter (Kbb: λ = 2.18 μm, Δλ = 0.4 μm) with an average spectral resolution of R ∼ 3600. A total of 76 individual exposures were taken with t_exp = 900 s in a 3 × 1 mosaic pattern. Observations taken the night of 2010 August 28 were subject to poor weather; the spatial resolution of these frames was generally poor (see Section 3.2), and thus these frames were removed from the analysis, leaving 54 total exposures. Observations on two of the remaining three nights were centered on the field, while the northwest and southeast pointings of the mosaic were observed on the other night. On each night, observations of telluric standard stars and blank sky were also obtained for use in calibration.

Table 1.  M31 Keck/OSIRIS Data

Night	No. of frames	Position
2008 Oct 21	13	center
2010 Aug 15	18	center
2010 Aug 28a	22	NW, SE
2010 Aug 29	23	NW, SE

Note.

^aNight removed from analysis owing to poor spatial resolution.

Download table as:  ASCIITypeset image

Total integration time on the center of the field consists of ∼35 exposures (9 hr), while the edges of the FOV were observed in 5–20 exposures (1.25–5 hr; see Figure 2, bottom panel).

Zoom In Zoom Out Reset image size

2.1.2. Keck/NIRC2

Supporting archival imaging of the nucleus from HST was also used. Optical images of the M31 nucleus were obtained with HST and downloaded from the Hubble Legacy Archive. The nucleus was imaged on 2006 June 15–16 using the Advanced Camera for Surveys (ACS) High Resolution Channel (HRC) with the F330W and F435W filters (PI: T. Lauer; Lauer et al. 2012). Total exposure time with the F330W filter was 8120 s over 12 exposures, while the total exposure time with the F435W filter was 2384 s over eight exposures.

The nucleus was imaged on 1994 October 02 using the Wide Field and Planetary Camera 2 (WFPC2) with the F555W, F814W, and F1042M filters (PI: M. Rich; Rich et al. 1996). Total exposure time for the F555W filter was 1680 s over six exposures, for the F814W filter was 1280 s over five exposures, and for the F1042M filter was 5000 s over 10 exposures.

The OSIRIS spectroscopic data were reduced using the OSIRIS reduction pipeline⁷ (Krabbe et al. 2004) to subtract a dark frame, correct bad pixels and cosmic rays, perform the data cube assembly and wavelength calibration, and remove the background sky and telluric OH emission. The pipeline version utilized for the reduction was v3.2. The pipeline includes internal logic to account for the hardware changes that have occurred since the data were taken and uses the modules appropriate for the date of observations.

For each data cube, the telluric absorption spectrum was modeled and removed using a combination of empirical and theoretical telluric absorption spectra. The empirical telluric spectrum was created using two standard stars with different spectral types as described in Hanson et al. (1996). A hot, early-type A0V star was observed (HD 209932) each night. This star has a featureless blackbody spectrum in the K band, except for a wide Brγ absorption feature at 2.166 μm. A solar-type star was also observed (HD 218633) each night and divided by a solar spectrum constructed from NSO data by ESO and convolved to match the OSIRIS spectral resolution (Livingston & Wallace 1991; Maiolino et al. 1996) to obtain a telluric spectrum. The solar-type telluric standard was used to replace the region around Brγ, from 2.151 to 2.181 μm, in the A0V spectrum. The corrected A0V spectrum was then divided by a 9500 K blackbody to obtain the final empirical telluric absorption spectrum.

However, changing atmospheric conditions throughout each night and a variable spectral resolution across the OSIRIS FOV led to large residuals when correcting for telluric absorption using the empirical telluric spectrum, which was observed at only a single detector location. This was specifically an issue where blueshifted stellar CO bandheads in the M31 spectra overlap with the large telluric residuals. To better model the changing telluric line depths, we created model telluric spectra with the molecfit package (Kausch et al. 2015; Smette et al. 2015), which generates synthetic telluric absorption spectra. Given an atmospheric profile for a given night and observatory location, molecfit uses a radiative transfer code to obtain absorption line depths and fits the output directly to the telluric features in a science spectrum. Before using molecfit, the OSIRIS M31 science frames were first mosaicked by night and by pointing, roughly corresponding with first and second half-nights, to better capture temporal telluric variations. A subset of 3 × 3 spaxels was extracted from a region of each mosaic where the stellar absorption lines did not overlap with the telluric lines. The median of these spectra was taken to create a single typical spectrum per half-night, and molecfit was run on each of these representative spectra to obtain the telluric absorption line depths. However, these model spectra could not be used directly in the telluric correction, as OSIRIS introduces a shape to the continuum that is not fit well by the molecfit modeling. The resulting molecfit line depths were combined with the continuum from the empirical telluric spectrum to produce a final telluric absorption spectrum, which was used to telluric-correct the M31 data.

The final mosaic was assembled by calculating shifts between frames using cross-correlations between each OSIRIS frame and the NIRC2 K-band image, rebinned to the OSIRIS spatial scale. Frames within a single dithered set (typically four to nine frames) were first mosaicked together, as the offsets between these frames are given by the input dithers. Each of these mosaics was then cross-correlated with the NIRC2 image to obtain its relative shift. The frames were then combined using a mean clipping method. The final OSIRIS data set consists of one fully combined data cube and three subset cubes, each containing 1/3 of the data, used for determining uncertainties.

The formal errors output by the pipeline include real spatial and spectral variations caused by removed OH sky lines, telluric absorption, interpixel sensitivity, and cosmic rays. However, the overall error is too large by nearly an order of magnitude relative to empirical error estimates, determined via the standard error on the mean of three subset cubes at each spaxel and spectral channel. As the magnitude of the flux errors impacts the errors on the kinematic measurements (Section 3.6), the errors need to be properly scaled. To combine the error variation captured by the formal pipeline errors with the more accurate magnitude of the empirical errors, we scale the pipeline errors by the ratio of the median of the two error arrays at every spaxel. We adopt these scaled flux errors for the remainder of the analysis.

The OSIRIS image quality was estimated by first collapsing the cube along the spectral direction. Then the OSIRIS image was compared to the NIRC2 image by convolving the NIRC2 image with an iteratively determined kernel that reduced the difference between the NIRC2 and OSIRIS images. A Gaussian kernel was used to represent the diffraction-limited core. The amplitude was set to 1 and not allowed to vary, while the width of the core Gaussian was allowed to vary. This gives an estimate of the spatial resolution of the OSIRIS images compared to the NIRC2 images. The resulting spatial resolution ranges between 50 and 350 mas FWHM for the individual exposures, and the median spatial resolution across all frames is 260 mas. The majority of the lower-quality frames are in the southeast pointing of the mosaic, due to worse observing conditions during the nights this pointing was observed. We experimented with removing various combinations of the lowest-quality frames, striving to maintain even coverage across the FOV. Ultimately, all frames from 2010 August 28, the night with the worst weather, were removed, resulting in a median spatial resolution of 205 mas among the remaining frames and a resolution of 124 ± 6 mas for the full combined mosaic (Figure 3). We proceeded with the analysis using this clipped set.

Zoom In Zoom Out Reset image size

The old stellar population of the bulge is the dominant source of light at 55 < r < 300'' (Kormendy & Bender 1999, hereafter KB99) and is a source of foreground contamination in observations of the nuclear disk. In addition to being a source of excess surface brightness, the bulge's slow rotation and internal dispersion can distort the kinematic signature of the nuclear disk. Both the surface brightness and the kinematic contribution of the bulge must be removed before extracting the nuclear disk kinematics.

First, we determined the fraction of bulge light in each spaxel. Previous studies have found that both the optical and the NIR surface brightness profile of the bulge is well fit by a Sérsic profile (KB99; Courteau et al. 2011; Dorman et al. 2013) of index ∼2, half-light radius ∼1 kpc, and half-light surface brightness 17.55–17.85 mag arcsec⁻² (optical V band, λ = 555 Å, and I band, λ = 806 Å) or 15.77 mag arcsec⁻² (L band, λ = 3.6 μm). The I-band surface brightness profile from Courteau et al. (2011, model M) was adopted and scaled to match the NIRC2 K-band wide-field image at a radius of 10'', or a distance at which the bulge dominates the surface brightness (Figure 4). The 1D Sérsic profile was converted to a 2D profile using the model M bulge ellipticity. We used the bulge P.A. of 66 from Dorman et al. (2013), as they fit the bulge profile with the same data set; the differences with their best-fit parameters are negligible. Both the 2D bulge profile and the NIRC2 image were smoothed to match the OSIRIS resolution, and the two smoothed images were divided to obtain the ratio of bulge luminosity to total luminosity in the nucleus (Figure 5).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Bulge-dominated spaxels were selected based on the bulge surface brightness ratio map. As no spaxels had a majority contribution from the bulge (Figure 5), those spaxels in which the bulge surface brightness ratio is at least 0.42 were selected to be representative of the bulge. This ensures that enough spaxels, roughly 160, were obtained to derive a high-quality bulge template spectrum.

Next, the intrinsic spectrum of the integrated bulge light was estimated using the penalized pixel fitting (pPXF) package (Cappellari & Emsellem 2004), which fits a linear combination of stellar templates convolved with a line-of-sight velocity distribution (LOSVD). pPXF was also configured to fit and add a fourth-degree Legendre polynomial to the convolved spectrum, to account for continuum shape differences between the stellar templates and the science spectrum (see Section 3.6 for more details). The bulge-dominated spaxels were fit with pPXF. The output best-fit linearly combined stellar templates for each spaxel, along with the additive Legendre polynomial, were used to create a template stellar spectrum for each bulge-dominated spaxel. Each resulting template spectrum reflected the flux level of the corresponding spaxel in the data cube. To normalize these spectra, the template spectrum for each spaxel was divided by its median. The median of all normalized bulge template spectra at each wavelength was then taken to obtain a median stellar template for the bulge.

The median bulge spectrum was dominated by two stellar template spectra: a late K giant and a late K dwarf. The late K giant is likely more representative of the actual bulge population, which is estimated to be roughly the age of the galaxy (Saglia et al. 2010). However, the late K dwarf star has stronger Na lines than that of the giant. The center of M31 is estimated to have enhanced metallicity compared to the Milky Way (Saglia et al. 2010), and thus the inclusion of the dwarf stellar template in the fit may be compensating for this different abundance pattern. Alternatively, varying Na line strengths may represent variation in the initial mass function of the old stellar population (McConnell et al. 2015, and references therein).

Finally, the bulge spectrum was convolved with a Gaussian LOSVD with a fixed velocity of −340 km s⁻¹, at the systemic velocity of the galaxy (R. Bender 2018, private communication), and a dispersion of 110 km s⁻¹ at all points in the field, or the dispersion in the bulge-dominated spaxels before bulge subtraction. No bulge rotation was included since previous work by KB99 found that the bulge is rotating slowly at 2.65(r/1'') km s⁻¹ and that including this rotation does not appreciably affect the bulge subtraction within the compact nuclear region. The final bulge spectrum was multiplied by the median flux in each spaxel and the appropriate surface brightness ratio to create a bulge cube. This cube was then subtracted from the observed spectral cube. An example spectrum, before and after bulge subtraction, and the scaled bulge template spectrum are shown in Figure 6.

Zoom In Zoom Out Reset image size

There are systematic sources of error in the bulge subtraction that we do not fully explore but that should be kept in mind. The edges of the OSIRIS FOV are subject to much more bulge subtraction than the center (Figure 5). Subtracting so much flux ends up amplifying the noisiness of the bulge-subtracted spectrum. This noise is not purely Poissonian, but also includes noise from imperfectly corrected sky and telluric features, which introduce larger systematic errors into the spectra at the edge of the FOV. In addition, the bulge is typically modeled as a Sérsic profile, which is sharply peaked at the origin. However, Bacon et al. (2001) also fit the bulge using a multi-Gaussian expansion (MGE), which models the bulge surface brightness as the sum of three Gaussians. Their MGE bulge profile is much shallower in the inner few arcseconds, with a difference of a few tenths of a magnitude at the origin from the equivalent Sérsic profile. Alternative bulge subtraction methods may subtract more or less flux than does our method, which may impact the kinematics derived from the bulge-subtracted spectra. In Appendix A, we explore in more detail several alternative methods of bulge subtraction and show that their impact does not change the overall conclusions of this paper.

We identified the position of the SMBH in our OSIRIS and NIRC2 data using the position of P3, the very compact (r ∼ 006; Lauer et al. 1998) young nuclear cluster that contains the SMBH (Bender et al. 2005). Unlike the old stellar population, which is bright in the NIR, P3 is NIR-dark and thus effectively invisible in our OSIRIS and NIRC2 data. However, P3 is bright in the UV; we adopted as the SMBH position the location of source 11, or the brightest UV source in the Lauer et al. (2012) analysis of HST/ACS F330W and F435W images. The NIR image cannot be directly aligned to the F330W image because there are no compact sources that are bright in both the blue and NIR filters that can be used to determine the transformation. Instead, we aligned images at successively longer wavelengths using an evolving set of compact, bright sources to fit the translation, scale, and rotation (Table 2), aligning adjacent-bandpass frames using only the sources common to both frames (see Appendix B). At least 16 bright compact sources are used for aligning each pair of frames. After alignment, each frame was transformed into the reference frame, F330W. The error on each frame's alignment was taken as the residual error, converted to arcseconds. The total error in alignment from F330W to the NIRC2 K-band image is 0033, which is smaller than the 005 spaxels in the OSIRIS data. The OSIRIS data cube was aligned to the NIRC2 K-band image by cross-correlating the OSIRIS cube, collapsed in the wavelength direction, with the final registered NIRC2 K-band image.

Figure 7 shows the region within 03 of P3 in each of the eight wavelength images after alignment. The pixel scaling was adjusted in each frame to match that of the F330W image. The SMBH coordinates were taken from the background-subtracted, subsampled F435W image (Lauer et al. 2012, private communication) and converted to the aligned K'-band image coordinates. The SMBH position is marked with the black cross in each panel in Figure 7.

Zoom In Zoom Out Reset image size

Fitting the precession rate of the disk requires a very high signal-to-noise ratio (S/N). The S/N of the data varies over the FOV and was calculated empirically. For each spaxel, we fit a straight line in a region of the spectrum dominated by the stellar continuum between 2.270 and 2.280 μm. After dividing by the fitted line, the S/N was estimated as the mean of the normalized flux divided by the standard deviation. The S/N map is shown in the middle panel of Figure 2.

The data were spatially binned using a Voronoi tessellation algorithm (Cappellari & Copin 2003), which optimally sums spaxels until a target S/N is reached while enforcing a roundness criterion on the resulting bins to ensure that the resulting bins are as compact as possible. Spaxels already at or above the target S/N remain unbinned. We used the S/N map thus generated to spatially bin the spaxels to a target S/N of 40. The resulting bins and their S/N are shown in Figure 8. Spaxels near P1, P2, and the SMBH remained unbinned, while the fainter outer regions of the eccentric disk and the innermost bulge were optimally binned. Spectra in the binned spaxels were summed and their errors combined in quadrature. Spaxels with a very low initial S/N (<2) on the very edge were masked and discarded in the following analysis.

Zoom In Zoom Out Reset image size

Kinematics were extracted using the pPXF method described in Cappellari & Emsellem (2004). This method fits the observed spectrum with a linear combination of spectral templates, convolved by an LOSVD, which we parameterized as a fourth-order Gauss–Hermite expansion of a Gaussian profile (Appendix A of van der Marel & Franx 1993). A fourth-degree Legendre polynomial was fit and added to the convolved spectrum to correct for mismatch between the continuum of the template and that of the science spectrum. The moments of the best-fit LOSVD (including v, σ, and the higher-order moments h₃ and h₄) and the weights of the spectral templates were returned. The Gemini GNIRS spectral template library (Winge et al. 2009), a library of evolved and main-sequence late-type stars observed in the K' band at R ∼ 5900, was used as inputs to pPXF. Only those templates observed in both the blue (2.15–2.33 μm) and the red (2.24–2.43 μm) modes were used, for a total of 23 template stars. The final template library covers G4–G5 II, F7–M0 III, K0 IV, and G3–K8 V stars. Independent fits were performed for each spaxel over the range 2.185–2.381 μm to obtain measurements of v, σ, h₃, and h₄. Example fits for spaxels near P1, in P2, and in the disk away from these regions are shown in Figure 9.

The templates preferred by the pPXF fits on the bulge-subtracted stellar population are nearly identical to those found in the fits to the bulge-dominated spaxels. Similar to the bulge population, the old eccentric disk is best fit by a late K giant and a late K dwarf template spectrum. This implies, as found by previous work (Saglia et al. 2010), that the old eccentric disk has a similar stellar population to that of the bulge. However, we defer a full analysis of the eccentric disk stellar population to future work.

Measurement errors on the LOSVD moments were assessed via a Monte Carlo (MC) analysis using 100 simulations. In each simulation, the input spectrum was randomly perturbed within the 1σ flux errors and the pPXF fits were recalculated. The uncertainty of each LOSVD moment (i.e., mean v, σ, h3, h4) was taken as the standard error of the mean of the fitted moments from the MC trials. Error maps are shown in Figure 10. The median errors are 11.3 km s⁻¹ on v and 16.7 km s⁻¹ on σ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition to measurement errors, mismatch between the spectral templates and the intrinsic stellar population can provide a source of systematic error. The magnitude of this error was assessed via a jackknife analysis. One spectral template at a time was dropped from the library used in the fitting routine. pPXF fits were run for each of these 23 template library subsets, and the systematic errors from template mismatches were calculated as the standard error of the mean of the resulting outputs. Error maps are shown in Figure 11, and the median errors are 1.3 km s⁻¹ on v and 1.6 km s⁻¹ on σ. These errors are nearly an order of magnitude lower than the measurement errors and thus negligible compared to the observational error. We disregard this source of error in the remainder of the analysis.

Zoom In Zoom Out Reset image size

The results of the pPXF kinematic fitting are shown in Figures 12 and 13. The velocity map shows a clear rotation signature, but with a slight asymmetry in the morphology of the velocity peaks. In addition, the peak of the velocity is offset from the P1 peak flux (Figure 14). There is also an asymmetry in the magnitude of the velocities on either side of the disk, with maxima of 223 ± 7 km s⁻¹ on the northeast side and −312 ± 16 km s⁻¹ on the southwest side. In addition, the zero of the velocity gradient does not coincide with the SMBH position, but is offset by 015 toward P1, which is consistent with previous long-slit measurements (Bender et al. 2005).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The dispersion is peaked close to but not at the position of the SMBH as is shown in the zoomed kinematic maps in Figure 13 and in the contour plot in Figure 14. The dispersion reaches a maximum of 381 ± 55 km s⁻¹ at an offset in (x, y) of (+009, −01) from the SMBH, or a distance of 013 on the P2 side. The wider dispersion peak (σ > 200 km s⁻¹) is concentrated on the P2 side and extends toward P1.

Overall, the maps for the higher-order moments h3 and h4 are quite noisy, and trends are difficult to discern. There is a possible enhancement of h4 on the P1 side of the SMBH, but the results are not robust. We find that the significance of the h3 and h4 maps is particularly sensitive to the bulge subtraction.

Closer examination of Figure 7 shows that there are clear differences in the structure of the stellar population at each wavelength; notably, P3 is bright and compact only in the F330W and F435W frames and is dark in the NIR. Figure 15 shows the inner 02 in the F330W, background-subtracted F435W (Lauer et al. 2012, private communication) and K' images, with the SMBH position marked. In contrast with the bluer images, there is little to no flux at the position of P3 and the SMBH in the NIR, indicating a separation of the young and old stellar populations in the nucleus.

Zoom In Zoom Out Reset image size

A clearer comparison of the difference between the flux maps at two different wavelengths can be seen in Figure 16. Lauer et al. (1998) observed the nuclear eccentric disk with HST/WFPC2 in three filters: F300W, F555W, and F814W. We show their F555W flux map, convolved to the OSIRIS spatial resolution, in comparison with the OSIRIS flux map. There are clear structural differences visible between the two flux maps, particularly around the P1 peak. In addition, the morphology differences seen around the SMBH in Figure 7 can been seen here: the secondary P2 peak is more northerly in the OSIRIS data than in the F555W image.

Zoom In Zoom Out Reset image size

This structural difference can also be seen in 1D cuts across each image (Figure 17). A 1D cut has been taken across each frame in Figure 7 at the P.A. of the HST/STIS long-slit observations (see Section 4.3) and passing through the SMBH. Each 1D cut has been flux-calibrated according to the published instrumental zero-points. The nucleus is much brighter in the NIR than at UV or optical wavelengths, so offsets in magnitudes (shown in the figure legend) have been added to each cut to equalize them on the P1 side. The UV peak at the origin in the F330W and F435W frames is apparent. The optical bandpasses show an extended, shallower flux peak roughly coincident with the SMBH position and extending slightly to the southwest. However, the H and K' flux cuts do not show this shallower peak and instead are flat from the SMBH position to the cutoff at the southwest edge.

Zoom In Zoom Out Reset image size

The nucleus of M31 has been well studied with both ground- and space-based spectroscopic observations. Several recent studies are summarized in Table 3. In this section, we compare the OSIRIS kinematics with two high-resolution data sets: HST/STIS long-slit observations reported in Bender et al. (2005, hereafter B05), and ground-based IFS observations reported in Bacon et al. (2001, hereafter B01).

Table 3.  Data Quality Comparison

Instrument	Filter	FWHM	IFU?	Reference
HST FOC	∼B	∼005	No	Statler et al. (1999)
CFHT SIS	∼I	063	No	KB99
CFHT OASIS	∼I	041–05	Yes	B01
HST STIS	∼I	012	No	B05
Keck OSIRIS	K	012	Yes	This work

Download table as:  ASCIITypeset image

4.3.1. Comparison with STIS

HST/STIS long-slit spectroscopy was obtained for the nucleus of M31 in 1999 July; the full analysis was published in B05. Spectra were taken in both blue (2900–5700 Å) and red (8272–8845 Å, the Ca triplet) modes; we confine our discussion solely to the red mode. The slit was oriented at a P.A. of 39°, along the major axis of the optical disk. The slit width was 01, and the estimated resolution is FWHM = 012.

We compare the kinematics reported in B05 with our OSIRIS kinematics in Figure 18. In addition to the original STIS kinematics, we show the STIS kinematics smoothed along the slit by a Gaussian kernel derived using the OSIRIS point-spread function (PSF). The comparison to the OSIRIS kinematics was derived by taking a cut across the kinematic maps in Figure 12 at the STIS P.A. Each point along the STIS slit is represented by the mean of three OSIRIS spaxels orthogonal to the slit to match the STIS slit width, and the errors are taken as the Monte Carlo errors of the same spaxels, added in quadrature.

Zoom In Zoom Out Reset image size

Though overall the kinematics reported by STIS and by OSIRIS are similar, there are several inconsistencies, particularly in the velocity and dispersion. Overall, while the velocity gradient from the OSIRIS data is well matched to that from STIS, the OSIRIS data show lower peak velocities on both the red- and blueshifted sides of the nuclear disk than the STIS kinematics. In addition, the wings of the velocity profile from OSIRIS are lower than those from STIS, which may point to residuals arising from differences in the bulge subtraction methods. The peak dispersion values are the same, within the errors, but the OSIRIS dispersion peak is broader, with a possibility of a secondary dispersion peak at 05, on the P2 side. This extended enhanced dispersion in the NIR-bright old stellar population could be probing the apoapse of the eccentric disk.

4.3.2. Comparison with OASIS

AO-corrected integral field spectroscopy was obtained by the OASIS instrument on CFHT and first presented in B01. Two different mosaics of different spatial resolution were presented, with resolution ranging from 041 to 05. Similar to the STIS data, the chosen wavelength range was 8351–9147 Å to optimize observations of the Ca triplet. OASIS provides a 4'' × 3'' FOV.

We show the bulge-subtracted OASIS kinematics in Figure 19. These kinematics were extracted along the OASIS kinematic major axis (P.A. = 564) and were initially presented in B01 (their Figure 10, solid line). We also show a cut across the OSIRIS kinematics at the same P.A. The red points are the original OSIRIS kinematics, and the red curve is the OSIRIS kinematics smoothed to the OASIS resolution. The B01 velocity and dispersion profiles have been shifted by −01, to better align the velocity gradient with that of the OSIRIS data and to account for differing SMBH positions. Overall, the OASIS velocity profiles agrees well with the smoothed OSIRIS data and are similar to the correspondence seen between the OASIS and STIS kinematics, as reported in B01 (their Figure 13). However, the OASIS and OSIRIS kinematics are not quite aligned along the radial direction. This offset is potentially due to a difference in SMBH position between the OASIS and OSIRIS data; shifting the comparison cut in the OSIRIS data by up to 01 from our SMBH position brings the data into closer alignment along the radial direction, though the correspondence is still not perfect. Also, smoothing the OSIRIS kinematics to match the OASIS resolution brings the two velocity profiles into agreement, but the reported OASIS dispersion peak is higher than that seen in the smoothed OSIRIS data. The peak dispersion in the smoothed 2D OSIRIS dispersion map is less than 250 km s⁻¹, or significantly lower than that seen by OASIS. It is unclear how the OASIS data are able to probe the dispersion at such high resolution given that beam smearing should lower the peak resolved dispersion. Differences in bulge subtraction may be the culprit, or else the OASIS resolution may be better than reported.

Zoom In Zoom Out Reset image size

The original disk models consisted of ∼10⁷ particles in disk-plane coordinates, which were rotated to two different orientations as described in PT03: the aligned model, matched to the orientation of the large-scale galactic disk, and the nonaligned model, with orientation parameters left free (Table 4, second and third columns). The aligned models were a poor fit to both our data and those of PT03; thus, we drop the aligned models from the discussion and focus on the nonaligned models. Throughout the model fitting, we adopt the best-fit parameters of the PT03 nonaligned model, including a black hole mass of M_• = 1.02 × 10⁸ M_⊙, and keep the parameters fixed unless otherwise specified. In order to fit the model to the OSIRIS data, the particles from the nonaligned model were first rotated to an orientation of our choosing as specified by three angles: the inclination of the disk with respect to the sky (θ_i), the angle of the ascending node in the sky plane (θ_l), and the angle from the ascending node to the periapse vector (θ_a). A solid-body rotation was also introduced by adding a fixed rotation speed, Ω_P (in units of km s⁻¹ pc⁻¹), to the disk-plane model velocities. A right-handed coordinate system was assumed, so positive precessions are counterclockwise. For a given set of model parameters, the particles were randomly perturbed in their sky-plane positions using a Gaussian kernel that matched the OSIRIS resolution in order to smooth out features at higher resolutions. The particles were then spatially binned using the tessellated pattern used for the data (Section 3.5), and the line-of-sight velocities were binned in increments of 5 km s⁻¹. The resulting model LOSVD was fit using the Gauss–Hermite expansion in the same manner as the observations.

Table 4.  Model Fitting Results

Parameter	Aligneda	Nonaligneda	Best Fitb
θl	−523	−428	−33° ± 4°
θi	775	541	44° ± 2°
θa	−110	−345	−15° ± 15°
ΩPc	...	...	0.0 ± 3.9

Notes.

^a PT03. ^bThis work. ^cNot fit in PT03; units are km s⁻¹ pc⁻¹.

Download table as:  ASCIITypeset image

In order to find the best-fit orientation and precession with our new data, a grid search was performed over Ω_P, θ_i, θ_l, and θ_a. The grid of angles was centered on the original best-fit nonaligned values and spread over a range of ±15° with a step of 5° for each angle. Precession values in the coarse grid ranged from −30 to +30 km s⁻¹ pc⁻¹ in steps of 5 km s⁻¹ pc⁻¹.

The goodness of fit was tested by calculating the χ², or the sum of the squared residuals between the OSIRIS flux or kinematics and that of the model, weighted by the squared flux or MC errors. Only the flux, velocity, and dispersion moments were used in the fitting process, as the h3 and h4 measurements have lower S/N. For the flux, the OSIRIS data and model were first normalized by dividing by the sum of the flux in each image. Only the inner 13 was used for the calculation. The resulting ${\chi }_{m}^{2}$ for each moment, m, was divided by the number of spaxels (1746) within the inner region minus the number of free parameters (4), to obtain a reduced ${\tilde{\chi }}_{m}^{2}$. The best-fitting model across all moments was selected based on a weighted sum ${\tilde{\chi }}^{2}={\sum }_{m}{\tilde{\chi }}_{m}^{2}{w}_{m}$ over the flux, velocity, and dispersion moments. The weight used was the inverse of the minimum ${\tilde{\chi }}_{m}^{2}$ value for each moment, which is equivalent to error rescaling, in order to prevent any one moment from dominating the fit.

The errors on the best-fit parameters were estimated by performing two Monte Carlo analyses. First, errors were calculated using 100 samples of the data. In each, the flux and kinematic maps were randomly perturbed using a normal distribution matched to the MC errors for the data, and the ${\tilde{\chi }}^{2}$ values were recalculated for every model in the grid. Notably, all Monte Carlo samples yielded the same best-fit model parameters as the original, nonperturbed data; however, the best-fit ${\tilde{\chi }}^{2}$ for each MC sample varied. To define a 68% (1σ) confidence interval from the data error analysis, we adopted the standard deviation of the best-fit ${\tilde{\chi }}^{2}$ values, or ${\rm{\Delta }}{\tilde{\chi }}_{\mathrm{data}}^{2}=0.01$. MC errors were also calculated for the models using 100 simulations. In each realization, the model x and y sky-plane position of each particle was perturbed randomly within a Gaussian kernel that matched the OSIRIS resolution. For each MC simulation, the model particles in the best-fit orientation were randomly perturbed within the same Gaussian kernel. To define a 68% (1σ) confidence interval from the model stochasticity, we adopted the standard deviation of the ${\tilde{\chi }}^{2}$ values from the best-fit models, or ${\rm{\Delta }}{\tilde{\chi }}_{\mathrm{model}}^{2}=0.08$. The errors from the data and the model MC simulations were added in quadrature, for a total error of ${\rm{\Delta }}{\tilde{\chi }}^{2}=0.08$. Errors on each of the fitted model parameters were derived using this error.

Figure 20 maps the sum of the weighted, reduced ${\tilde{\chi }}^{2}$ values for the full range of model parameters that were explored. Figure 21 shows ${\tilde{\chi }}^{2}$ values as a function of each parameter in 1D, with the other parameters held constant at their best-fit values. The best-fit ${\tilde{\chi }}_{m}^{2}$ for each moment and the combined ${\tilde{\chi }}^{2}$ are shown in Table 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our flux, kinematics, and residuals from the original nonaligned models of PT03 are shown in the left panels of Figures 22–24. We note some discrepancies, particularly in the flux and velocity comparisons. (Following PT03, we note that comparisons between the data and models are only valid within the inner 13.) The flux residuals may be due to the difference in morphology of the stellar population at different wavelengths (Section 3.4), as the flux models were originally fit to imaging data taken by HST in filter F555W (originally reported in Lauer et al. 1998) and the kinematics were fit to Ca triplet long-slit data from KB99.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The best-fit models and residuals are shown in the right panels of Figures 22–24. The best-fit orientation is [θ_l, θ_i, θ_a] = [−33° ± 4°, 44° ± 2°, −15° ± 15°]. The 1σ errors are derived using the ${\rm{\Delta }}{\tilde{\chi }}^{2}$ above. Overall, the fitting preferred smaller values of θ_i and larger values of θ_l and θ_a than those found in PT03 (see Table 6 for a comparison).

Table 6.  Literature Parameters

Parameter	[1]	[2]	[3]	[4]	[5]	[6]	[7]
θl	...	−48° b	...	...	−2734 b	−48° a,b	−35°
θi	77° a	55° ± 5°	...	77° a	5154	525	57°
θa	...	...	−21° a,b	...	...	...	−34°
ΩPc	≲47.7d	3	16	13.6	16	36.5	...

Notes.

^aAssumed, from the literature. ^bConverted from format given in literature. ^cUnits are km s⁻¹ pc⁻¹. ^dGiven as a function of sin θ_i; adjusted to match our best-fit value.

References. [1] Sambhus & Sridhar 2000; [2] B01; [3] Jacobs & Sellwood 2001; [4] Salow & Statler 2001; [5] Sambhus & Sridhar 2002; [6] Salow & Statler 2004; [7] Brown & Magorrian 2013.

Download table as:  ASCIITypeset image

The reduced ${\tilde{\chi }}^{2}$ values for the flux and velocity residuals preferred a positive precession, while the dispersion residuals preferred a negative precession, but overall the preferred precession is Ω_P = 0.0 ± 3.9 km s⁻¹ pc⁻¹. We show the data, model, and residuals for the best-fitting orientation and precession in the right panels of Figures 22–24.

We also fit the precession using the 1D method formulated in Tremaine & Weinberg (1984, hereafter TW84) and modified by Sambhus & Sridhar (2000). Briefly, the original method requires 1D profiles of both the surface brightness and the LOS velocity along the line of nodes, or the intersection of the sky and disk planes. The modified method allows for the profiles to be taken along the P1−P2 line. In either, the formulation is

$$\begin{eqnarray*}&&{{\rm{\Omega }}}_{P}\,\sin \,{\theta }_{i}{\int }_{-\infty }^{\infty }r\,{\rm{\Sigma }}(r)\,{dr}={\int }_{-\infty }^{\infty }{v}_{\mathrm{LOS}}(r)\,{\rm{\Sigma }}(r)\,{dr},\end{eqnarray*}$$

where r is the projected distance from the black hole, Σ(r) is the surface brightness, and v_LOS(r) is the line-of-sight velocity.

We take the surface brightness and LOS velocity profiles along the P1−P2 line, or parallel to the long axis of our FOV, and calculate the precession along the profile intersecting with the SMBH position. The method allows determination of the precession along strips parallel to this line. The precession along two strips on either side of the profile intersecting the SMBH position is also calculated, for a total of five determinations. The standard deviation is taken as the error. Using the best-fit value for θ_i, we derive Ω_P =−18 ± 5 km s⁻¹ pc⁻¹.

However, the original method by TW84 assumes that the disk is thin; the models from PT03 show that the best-fit disk is quite thick (h/r ∼ 0.4). We check the results of the TW84 method by calculating the precession for the models of the same orientation as the best-fit model into which a known precession value has been injected. We find that the values obtained using the TW84 method for these models are systematically too high. A line is fit to the output TW84 precessions to calibrate the method for the best-fit orientation. We extrapolate the fit to obtain the calibrated TW84 precession for our data: 62 ± 5 km s⁻¹ pc⁻¹. A model is generated using this precession, along with the best-fit orientation angles, and the goodness of fit to the data is calculated (Table 5); this precession gives a significantly worse fit to the data than the best-fit precession does, particularly to the velocity map. A comparison of the data and the models derived using the TW84 precession is shown in Appendix C.