2 Model construction

2.1 Deprojection and disk potential

Figure 1: Derivation of the nuclear disk, and its (smooth) gravitational potential. a) Observed sky brightness distribution, from Lauer et al. (1998); the dotted curve has magnitude 14.3, and successive isocontours differ by 0.25 mag. b) Brightness profiles along the P1-P2 line: observed (dotted curve), the Sérsic bulge (dashed-dotted curve), and the bulge-subtracted nuclear disk (solid curve). c) Brightness distribution of the disk, viewed face-on; the isocontours follow the same convention as in a). Isocontours of the disk-potential are displayed in d); the units are such that the deep minimum near P1 has depth equal to unity, and successive isocontours mark increments of 0.05. In c) and d), "X'' marks the location of the center of mass.

Figure 1a shows the nucleus of M 31, plotted from the V-band, HST observations of Lauer et al. (1998). The UV cluster and the MDO are at the origin. P2 is near the MDO, with sky coordinates ( $0 \farcs 023,~ 0\arcsec$), and P1 is located at ( $-0 \farcs 48,~0\arcsec$). The bulge was assumed to be spherical, with a Sérsic brightness profile (Sérsic 1968; Kormendy & Bender 1999) - see Fig. 1b. The center of mass (COM) of the bulge, disk, and MDO was required to coincide with the bulge center; this common location was determined, by an iterative method, at ( $-0 \farcs 0684, 0\arcsec$), in agreement with Kormendy & Bender (1999). With one notable exception, (Bacon et al. 2001), all investigations have assumed that the nuclear disk is coplanar with the much larger galactic disk of M 31. We obtained very poor results with this assumption. Therefore, we resolved to determine the inclination and orientation of the nuclear disk, based on the photometry, similar to Bacon et al. (2001). The disk light covered an approximately elliptical region, with a ragged edge. The plane in which the best-fit ellipse (to the edge) deprojected to a circle was defined to be the disk plane; its inclination (i), and PA of the line of nodes, were $\simeq 51.54{\degr}$, and $\simeq$ $62.66{\degr}$, respectively. The face-on view of the disk, shown in Fig. 1c, has mass $\simeq$ $2.15\times 10^7~M_{\odot}$.

To minimize edge-effects, the self-gravitational potential was evaluated in the disk plane, at 10⁴ grid points within a central square, of side equal to $2 \farcs 28$. However, the Newtonian $\vert{\vec r} - {\vec r'}\vert^{-1}$ contributions from the entire disk of Fig. 1c, which has diameter $\simeq$ $3 \farcs 6$, was used. The grid values were fit to a 20-th order polynomial function of the Cartesian coordinates, $\Phi_{\rm d}(r)$, a contour plot of which is displayed in Fig. 1d. The polynomial form smoothed the potential, facilitated coding of the integrator, and checking of the integrated orbits in the nearly Keplerian limit (Sridhar & Touma 1999). Figures 1c and 1d can be imagined as either snapshots of a rotating configuration, or as steady images in a frame rotating with some angular speed $\Omega$, about an axis normal to the disk plane, and passing through the COM. The forces on a test star include the gravitational attractions of the MDO and disk, as well as centrifugal and Coriolis forces. The contribution of the bulge was ignored, because it is so much smaller than the other forces.

2.2 The orbit library

Figure 2: Orbits in the rotating frame and photometric fits. The axes in all panels are sky positions. a) and b) show prograde and retrograde loop orbits, respectively, as seen on the sky, for $\Omega = 16 \mbox{~km~ s${}^{-1}$ ~ pc${}^{-1}$ }$; the parent (resonant) orbits are overdrawn as the solid curves. The photometry in  c) is from (Lauer et al. 1998), smoothed with a Gaussian beam of $FWHM=0 \farcs 17$. The (bulge-subtracted) light in the region enclosed by the dashed box was employed in our Schwarzschild-type iterative method. d) is our model disk, including the bulge. The dotted lines in both c) and d) have magnitude equal to 14.5, and successive isocountours differ by 0.25 magnitudes. The brightness is displayed in the "negative'' mode, to better emphasize the distribution.

Orbits were computed in the rotating frame by numerically integrating the equations of motion,

$$% \ddot{\vec r} = -\frac{GM(\vec {r} - \vec {r}_{\rm MDO})}{\... ...ga^2{\vec r} -2\Omega~\left(\hat{z}\times \dot{\vec r}\right),$$ (1)

using a 4th-order, adaptive step size, Runge-Kutta scheme. The global structure of orbits was explored by studying Poincaré surfaces of section. The principal families of orbits were lenses and loops. Lens orbits change the sign of their orbital angular momentum (Sridhar & Touma 1997, 1999). Stars on such orbits will collide with the MDO, in the time it takes an orbit to precess. These time scales do not exceed a million years, even for quite large orbits; dwarf, as well as giant stars on lens orbits will be lost to the MDO (if not tidally disrupted before). Hence lens orbits were excluded from our modeling. Other orbits that were also omitted included chaotic orbits, and those parented by higher order resonances. The loops orbits were of two kinds: prograde and retrograde, some of which are shown in Figs. 2a and 2b; these were the only orbits included in our orbit library. The kinematic model of Tremaine (1995), the Kepler-averaged dynamics of Sridhar & Touma (1999), and studies of slow, linear modes by Tremaine (2001), all suggest the use of the prograde loops as the back bones of the orbit library. The necessity of including retrograde loops is less obvious, and was stimulated by the investigations of Touma (2001). It turned out that the retrograde loops significantly improved fits near P2.

For each value of energy, loops of two or three different "thicknesses'' (i.e. deviations from the parent loop) were computed. Each orbit was sampled, and populated with "stars'', spaced apart uniformly in time. All stars in an orbit are accorded the same (unknown) mass; this is not a restriction, because in a collisionless system, the relevant physical quantities are the mass per orbit. The numbers of "stars'' in an orbit was chosen proportional to the inverse square of the energy (approximately, square of the "semi-major axis'') of the parent loop; thus 25stars sufficed for an orbit with $a=0 \farcs 02$, whereas an orbit with $a=0 \farcs 6$ was sampled by more than 10 000 stars. Altogether, the positions and velocities of $\sim$237 000 stars, populating 50 prograde orbits and 20 retrograde orbits, were recorded.

2.3 Richardson-Lucy deconvolution

Orbit masses were determined by iteratively imposing on the model, consistency with the bulge-subtracted sky brightness of a region covered by the orbits; the dashed box of Fig. 2c encloses this region. The box was divided into 112 equal square cells, each of side $0 \farcs 09$; each cell was small enough to give good resolution, and large enough (16 pixels) to keep pixel noise levels low. The "observed'' mass per cell, $\mu_j$ (for $j=1\ldots 112$), was obtained from the observed light, by multiplication with $\Upsilon_V$; these numbers composed our basic data. We defined m_i (for $i=1\ldots 70$) as the mass of orbit i, that also lies within the box; the total mass in the orbit exceeds m_i. We normalized $\sum_{i =1\ldots 70} m_i = \sum_{j=1\ldots 112}\mu_j = 1$, to unit mass in the box. A linear relationship, $\mu_j ~=~ \sum_{i=1\ldots 70} K(j\vert i)~m_i$, exists between the "observed masses'' $\mu_j$, and the unknown masses m_i. The positive kernel, K(j|i), is known from the orbit library. It has the property, $\sum_{j=1\ldots 112}K(j\vert i)=1$, for all $i=1\ldots 70$. An initial guess, $\{m_i^{\rm g}\}$, was iterated by the RL algorithm (Richardson 1972; Lucy 1974). The problem being overdetermined, about 5000iterations ensured good convergence to some $\{m_i^{\rm f}\}$. Velocities were then transformed to the inertial frame. Rescaling of $\{m_i^{\rm f}\}$ to physical values, and including the portions of orbits outside the box, provided a numerical distribution function. The entire process, beginning from the selection of an orbit-library, was repeated for several values of $\Omega$, between 5 and $25\mbox{~km~ s${}^{-1}$ ~ pc${}^{-1}$ }$. For any chosen value of $\Omega$, the final set of orbit masses, $\{m_i^{\rm f}\}$, corresponds to a prediction for the cell masses, $\mu_j^{\rm f} ~=~\sum_{i=1\ldots 70} K(j\vert i)~m_i^{\rm f}$, which should be compared with the data, $\{\mu_j\}$. For models with $\Omega ~=~ \{15, 16, 17\}$, the root-mean-squared deviation in mass per cell are, $\{0.26, 0.20, 0.28\}$; other values of $\Omega$ resulted in very poor models.