{'fraction_non_alphanumeric': 0.045665423451287185, 'fraction_numerical': 0.02740380241972164, 'mean_word_length': 4.006033697632058, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Three-dimensional Keplerian orbit-superposition models of the nucleus of M31

We present three-dimensional eccentric disc models of the nucleus of M31, modelling the disc as a linear combination of thick rings of massless stars orbiting in the potential of a central black hole. Our models are nonparametric generalisations of the parametric models of Peiris&Tremaine. The models reproduce well the observed WFPC2 photometry, the detailed line-of-sight velocity distributions from STIS observations along P1 and P2, together with the qualitative features of the OASIS kinematic maps. We confirm Peiris&Tremaine's finding that nuclear discs aligned with the larger disc of M31 are strongly ruled out. Our optimal model is inclined at 57 degrees with respect to the line of sight of M31 and has a position angle of 55 degrees. It has a central black hole of mass 10^8 solar masses, and, when viewed in three dimensions, shows a clear enhancement in the density of stars around the black hole. The distribution of orbit eccentricities in our models is similar to Peiris&Tremaine's model, but we find significantly different inclination distributions, which might provide valuable clues to the origin of the disc.


INTRODUCTION
As the nearest spiral galaxy to our own, M31 allows the study of a galactic centre in unmatched detail. A particularly striking feature of high-resolution V -and I-band photometry of the central few arcseconds of M31 is its double nucleus (Light et al. 1974;Lauer et al. 1993;Bacon et al. 1994;King et al. 1995): there are two peaks in surface brightness, the brighter of which (known as P1) is extended and offset ∼ 0. ′′ 5 from the fainter peak (known as P2), which is elongated and close to the photometric centre of the bulge. Early spectroscopic observations resolved steep gradients in the stellar rotation velocities and a prominent peak in the stellar velocity dispersion, which hints at the presence of a massive black hole (Dressler & Richstone 1988;Kormendy 1988;Bacon et al. 1994;van der Marel et al. 1994). Tremaine (1995) put forward an elegant explanation for these observations: the nucleus is a massive disc of red stars on eccentric, nearly Keplerian, approximately aligned orbits around a central black hole located at P2; P1 is generated by orbital crowding of stars lingering at apocentre. Subsequent observations have been consistent with his model. Lauer et al. (1998) observed the nucleus with the Wide Field Planetary Camera 2 (WFPC2) on the corrected Hubble Space Telescope (HST), confirming the bimodal structure. Kormendy & Bender (1999) obtained spectroscopy with the Subarcsecond Imaging Spectrograph (SIS) finding an asymmetric rotation curve and a constant colour across the nucleus, showing P1 and * Email: ckb@thphys.ox.ac.uk, magog@thphys.ox.ac.uk P2's colours to be consistent with each other but not the bulge or a globular cluster. Further kinematics were recorded with the Faint Object Camera , the integral field spectrograph OASIS (Bacon et al. 2001) and the Space Telescope Spectroscopy Imaging Spectrograph (STIS) on HST (Bender et al. 2005). All of these observations show that the kinematic centre of the nucleus is very close to P2, as predicted by T95's eccentric disc model.
Further refinement of this picture has come from studying the nucleus at ultraviolet wavelengths. Almost all of the UV emission from the nucleus comes from a tiny (< 0. ′′ 1) source located at P2 (King et al. 1995;Lauer et al. 1998;Brown et al. 1998), whose optical-UV colours and spectra are consistent with a population of A stars (Lauer et al. 1998;Bender et al. 2005), distinct from the red K-type spectrum of the rest of the nucleus (Bender et al. 2005). Bender et al. (2005) find that this compact, young, blue population has a maximum velocity dispersion of σ = 1183 ± 201 km s −1 , significantly higher than the red stars' ≈ 250 km s −1 . They label the UV peak P3 and have found that its photometry and kinematics are well modelled by a separate, almost-circular disc of blue stars around a central black hole of mass (1.1-2.3) × 10 8 M ⊙ . This provides very strong evidence in support of the presence of a supermassive black hole, which is the most fundamental requirement of T95's model. Most recently, Lauer et al. (2012) has found that the surface brightness profile of the young population is described by an exponential profile of scale length 0. ′′ 075 ± 0. ′′ 01.
Tremaine's original eccentric disc model consisted of three Keplerian orbits, coplanar with the main disc of M31, projected onto the sky and convolved with Gaussian point spread functions c 2013 RAS (PSFs). The original model included neither disc self-gravity nor a realistic treatment of the disc's internal velocity dispersions but was still able to broadly reproduce kinematic features. To be stable over many dynamical times such a disc would require apsidal alignment to be maintained; T95 proposed that this could be achieved by the disc self-gravity and argued that two-body relaxation would lead to a disc thickness of ≈ 0.3 times the disc radius. Estimates of the disc mass from mass to light ratios place it at around 10 7 M ⊙ . This is substantial enough to affect the dynamics of the disc although the nucleus falls within the sphere of influence of the black hole, which dominates the orbits and imposes regularity. A self-gravitating eccentric disc mode maintaining orbital alignment would also precess under a uniform pattern speed. Since the T95 model, several two-dimensional self-gravitating models with these properties have been constructed (Statler 1999;Bacon et al. 2001;Salow & Statler 2001;Sambhus & Sridhar 2002;Salow & Statler 2004). These massive models have found a variety of disc masses, pattern speeds and orbital distributions. Peiris & Tremaine (2003) took a different approach. They revised the 1995 model and constructed fully three-dimensional models with a parametric distribution function that ignored the self gravity of the disc; the gravitational potential in their models is due solely to the central black hole, which greatly simplifies the modelling procedure. Their models are the most successful to date at fitting the observed kinematics. They also found that thickened disc models that were misaligned with respect to the large-scale M31 disc produced significantly better fits than coplanar models, echoing a result seen in the 2d models.
The current picture of the nucleus (see figure 1) has the black hole (hereafter BH) embedded in P3, which is explained by a flat, circular exponential disc of blue stars 0. ′′ 033 from the photometric bulge centre. This young, blue disc is surrounded by the larger, red, eccentric disc. P1 is made of stars crowding at apoapsis, while P2 is now identified as stars at pericenter in the elongated region on the anti-P1 side of P3. The system remains of great interest: the origin of the nuclei and the relationship between the red and blue populations are unexplained and the mass and pattern speed of precession of the disc have not been pinned down consistently. Understanding the dynamics of the disc will help determine the BH mass more accurately which is of use in better determining the relationship between a BH and the host galaxy. It is also of interest to confirm whether the eccentric disc is aligned with the main disc of M31.
In this paper we present a natural development of the modelling approach started by Peiris & Tremaine (2003). Like them, we ignore the self-gravity of the disc, leaving the construction of fully self-gravitating models for a subsequent paper. But instead of considering parameterised functional forms for the phase-space distribution function (hereafter DF) of the nucleus, we model the DF non-parametrically as a mixture of Gaussian rings whose amplitudes are allowed to vary: one of our goals is to "let the data speak for themselves" and then examine closely the structure of the DF. The paper is organised as follows. In section 2 we summarise the data we use and additional post-processing applied. Section 3 describes our modelling procedure and section 4 our results. Section 5 sums up.
Throughout this paper we adopt a distance of 770kpc to M31.  (Lauer et al. 1998) annotated to show the relative positions of P1, P2 and P3 in arcseconds. The pair of white lines represents the positioning and width of the STIS slit (Bender et al. 2005). The brightest point in the image is at P1 and is marked with cross. Centred in the image and marked with a cross is the location of the black hole and P3. P2 is on the anti-P1 side of the black hole.

OBSERVATIONAL DATA
Our models use observations from three instruments: the HST photometry recorded with WFPC2, detailed in Lauer et al 1998 (hereafter L98), the high resolution kinematics from STIS (Bender et al 2005, B05) and the kinematic maps from the integral field spectrograph OASIS (Bacon et al 2001, hereafter B01).

Photometry
We use the deconvolved V -band (F555W) WFPC2 photometry of L98. This image is 1024×1024 pixels in size with a pixel scale of 0. ′′ 011375. P3 is located at (513, 518). We note that the orientation of the image has North 82.3 • clockwise from the y-axisthe data has been rotated away from the original alignment of the CCD, so that the corners of the image are blank -but we make no attempt to reorient the image. This image has already been reduced and deconvolved but we perform some processing on the image in order to reduce this data to a tractable number of observables for our fitting procedure. Bright foreground stars are first masked out and we crop the image to a 512×512 pixel region centred on P3. This is a region of side 5. ′′ 824 and incorporates the nucleus and the inner parts of the bulge and removes the blank corners of the image. This cropped image is then rebinned to "super-pixels" of side 2 l pixels containing the mean of the 2 2l − n pixels, with n the number of masked pixels falling in the super-pixel. We consider two schemes for choosing super-pixels: in the first we simply take l = 2 everywhere so that super-pixels contain a 4 × 4 pixel region. In the other l is chosen based on the surface brightness of the image scaled by r −1/2 , where r is the radius from P3, to provide a good balance of detail around P2 and P1. This gives l = 0 (the resolution of the deconvolved image) at P1 and P2 and l = 4 in the outermost parts of the image. A total of 4096 super-pixels are generated in this scheme. Ajhar et al (1997) and Lauer et al (1998) showed that that at the resolution of WFPC2 the underlying number of stars in M31 shows strong surface brightness fluctuations. The limited number of stars is the dominant residual between the structure of the nucleus and any model we fit. We treat the noise from the fluctuations as the sole source of our errors in our photometry and ignore photon noise so that the fractional error in a spatial bin is given by σ =N −1/2 , whereN is the effective number of stars within the spatial bin. L98 gives the effective magnitude of each SBF "star" as m I = 23.4, or 224 L ⊙ . The error associated with a spatial bin containing n L ⊙ is therefore √ 224n L ⊙ .

STIS kinematics
The STIS kinematics are taken from Table 5 of Bender et al. (2005). These consist of the bulge-subtracted Gauss-Hermite coefficients derived with the Fourier Correlation Quotient Method at the 22 positions listed for each of V , σ, h 3 and h 4 and their accompanying errors. We adopt the same slit widths of 0.1 ′′ and position angle of 39 • . The quoted positions are given with respect to the location of P3. We use equation (21) of Peiris & Tremaine (2003) to model the STIS PSF. When this is convolved with a a 0. ′′ 1-square top hat to represent the slit, the result is well approximated by a double Gaussian with amplitudes I 1 = 0.24, I 2 = 0.76 and dispersions σ 1 = 0. ′′ 042 and σ 2 = 0. ′′ 087 respectively. We assume that equation (1) as the effective PSF of the STIS observations. A more sophisticated treatment might take account of asymmetries in the PSF and the variations in spatial binning along the slit.

OASIS kinematics
We also make use of the kinematics of the integral field spectrograph OASIS, which were kindly provided by Eric Emsellem. We opt to use the higher S/N data set, M8. This consists of V , σ, h 3 and h 4 values derived from spectra taken at 1123 positions, spaced by 0. ′′ 09. We registered the image with the WFPC2 data using a similar process to B01. We found a close match in angle (0.7 • ) and a small offset of (-0. ′′ 02, -0. ′′ 03) between the two images. Like B01 we assumed that the OASIS measurements have a PSF that can be described by a sum of three Gaussians and allowed the parameters of these Gaussians to float freely in the registration process, taking the form The resulting PSF differs from that found by B01 with σ 1 = 0. ′′ 230, σ 2 = 0. ′′ 587, σ 3 = 0. ′′ 440 and I 2 /I 1 = 0.836 and I 3 /I 1 = 0.057.

MODELLING PROCEDURE
Our models are straightforward generalisations of those of PT. In particular, we ignore the self gravity of the disc and the gravitational influence of the bulge and assume that the potential is purely Keplerian. The mass of the central black hole M • is the single free parameter in our model potential. We assume that the BH is located at P3.

Coordinate systems
Following PT03 we use three different coordinate systems: "orbit plane", "disc plane" and "sky plane". All three coordinate systems have origin O coincident with the BH. Our models have biaxial symmetry. This provides a natural definition of the "disc plane" (x, y, z) coordinate system: the model is symmetric under reflections in the (x, y) and (x, z) planes. The orbit-plane and sky-plane coordinate systems are defined as follows. In the potential of the BH all orbits are Keplerian ellipses. An orbit with semi-major axis a and eccentricity e defines a coordinate system (x ′ , y ′ , z ′ ) in which the Oz ′ axis is parallel to the orbit's angular momentum vector and the Ox ′ axis points towards pericentre. That is, where the eccentric anomaly E is related to the mean anomaly M through Kepler's equation M = E − e sin E. In this coordinate system the apocentre of the orbit is located at (x ′ , y ′ , z ′ ) = (−a(1 + e), 0, 0) and the pericentre at (a(1 − e), 0, 0). Each star defines its where the angles ω, I and Ω are the star's argument of periapsis, the inclination and longitude of the ascending node, respectively. Any orbit in the Keplerian potential of the BH can be labelled by the five integrals of motion (a, e, ω, I, Ω), but it proves convenient to replace e and ω by the eccentricity vector e ≡ (e cos ϖ, e sin ϖ, 0) = (e x , e y , 0), where the longitude of periapsis ϖ = ω + Ω. The vector e points from the BH towards the projection of the pericentre onto the z = 0 disc plane. Its magnitude is the scalar eccentricity e of the orbit. Projected, sky-plane coordinates (X,Y, Z) are related to disc plane coordinates via in which (θ a , θ i , θ l ) are the three Euler angles specifying the orientation of the disc with respect to the observer's reference frame.
The (X,Y ) plane is the sky plane, with the positive X axis pointing west and positive Y axis north. The Z axis then points towards the observer; the line-of-sight velocity is therefore V los = −Ż, following the usual convention that that receding objects have V los > 0.

Distribution function
We model only the old stellar population of the disc; the young stars around P3 and the bulge are treated as contaminants (see §3.4 below). Assuming that the old stars are collisionless and homogenous, their dynamics can be completely described by a distribution . We further assume that the system is in a steady state. By Jeans' theorem the DF must be expressible as a function of the integrals of motion (a, e, I, Ω). The number of stars within a small phase-space volume element d 3 xd 3 v is then the 1 2 (GM • ) 3/2 a 1/2 sin I factor coming from the Jacobian relating (x, v) to (a, e, I, Ω, M).
We assume that the distribution function can be decomposed into a weighted sum of rings, in which each ring has a uniform distribution of Ω ∈ [0, 2π) and Gaussian distributions in a, e, I, and it is understood that f j (a, e, I) = 0 if a < 0 or |e| > 1. The normalisation factor N j is included to give each ring unit total luminosity. We use a Monte Carlo method to construct the rings, which avoids explicit calculation of N j .
Our rings span a (20 × 9 × 4) grid in (a j , e x, j , σ I, j ), with mean semimajor axis a j running logarithmically from 0. ′′ 03 to 10 ′′ , mean x component of eccentricity vector e x, j from −0.8 to +0.8 in steps ∆e = 0.2, and the dispersion in inclination σ I, j drawn from {12 • , 24 • , 36 • , 48 • }. For the spreads in (a, e) we take σ a, j = 0.8a j ∆ log a and σ e, j = 0.8∆e = 0.16. All rings have mean e y, j = 0. Together with the zero mean inclination of each ring this makes our models symmetric under reflection in the y = 0 and z = 0 planes. Rings with mean e x > 0 are "aligned" in the sense that their apocentre (and therefore peak density) lies somewhere along the negative x axis. Rings with mean e x < 0 have the opposite orientation.

Observables of each ring
We compute the projected properties of each ring f j (equ. 8) by sampling it with 10 6 points, drawing N o = 5000 values of (a, e, I, Ω, ω) from the distribution f j (a, e, I)a 1/2 sin I that appears on the right-hand side of (6), then sampling 200 points equispaced in mean anomaly M along each of these orbits. To avoid sampling artefacts we select the initial value of M for each orbit at random.

WFPC
We use this sample of 10 6 points to estimate the projected surface density (or, equivalently, the zeroth velocity moment) on a fine grid of 0. ′′ 011375 × 0. ′′ 011375 pixels on the (X,Y ) sky plane by simply counting the number of points that fall within each pixel. The contribution of the j th ring to each WFPC2 "superpixel" is found by taking the mean value of µ 0,fine j (X,Y ) over the relevant range of (X,Y ). We do not carry out any PSF convolution for the WFPC2 photometry because the image of Lauer et al. (1998) has already been deconvolved.

OASIS
The OASIS kinematics of Bacon et al. (2001) consist of measurements of mean velocity V (X,Y ) and velocity dispersion σ(X,Y ), together with Gauss-Hermite coefficients h 3 (X,Y ) and h 4 (X,Y ). We ignore their h 3 and h 4 measurements and assume that their measured V (X,Y ) and (V 2 + σ 2 )(X,Y ) distributions probe directly the (PSF-convolved) first-and second-order moments of the lineof-sight velocity distribution (LOSVD) at the point (X,Y ) on the sky. This is a good approximation provided the underlying PSFconvolved LOSVDs are reasonably close to Gaussian. For modelling purposes it is more natural to consider luminosity-weighted moments where I(X,Y ) is the underlying surface brightness. We obtain I(X,Y ) by convolving the WFPC2 image with the OASIS PSF (2). The observational errors on these first-and second-order luminosity-weighted velocity moments are obtained by adding the uncertainties on I, V and σ in quadrature in the obvious way: Following (9) above, the contribution of the j th ring to the luminosity-weighted first and second moments of the line-of-sight velocity distribution, are estimated by weighting each of the 10 6 sample points by V Z and V 2 Z , respectively. Having these µ i,fine j (X,Y ) distributions we convolve with the OASIS PSF (2) to obtain the contribution the ring makes to the model's predictions for the first and second velocity moments of the OASIS kinematics.

STIS
Our first attempt at fitting the STIS kinematics was based on the same assumption that we could use the STIS V and σ measured by B+05 as direct estimates of the (STIS PSF-convolved) first-and second-order moments of the LOSVD. That did not work well: at STIS resolution the LOSVDs are far from Gaussian, as we show on Figure 2, and the strong high-velocity wings caused by the BH mean that it is not possible to obtain reliable estimates of the first and second moments from the observed spectra. This last point was one of the motivations for the introduction of Gauss-Hermite series (van der Marel & Franx 1993) to parametrize LOSVDs. Therefore we fit our models to B+05's Gauss-Hermite parametrizations of the LOSVDs along the STIS slit. Our method for fitting the Gauss-Hermite coefficients follows the same lines used in other orbit-superposition models (e.g., Cretton et al. (1999)), but taking extra care to treat the normalisation of the LOSVDs correctly.
We recall some details of Gauss-Hermite expansions. Suppose that we are given an LOSVD L 0 (v), normalised such that L 0 (v) dv = 1. The Gauss-Hermite expansion of this L 0 (v) is where w ≡ (v − V )/σ, α(w) = e −w 2 /2 / √ 2π is the standard Gaussian and the H j (w) are Hermite polynomials. We adopt vdMF93's normalisation for the latter. Using the orthogonality properties of the H j , it is easy to show that the Gauss-Hermite coefficients h l are given by That is, there is a different Gauss-Hermite series (13) for each choice of (γ,V, σ). The "line-strength" parameter γ simply scales all the h j , but it proves important as we shall now see. A particularly natural choice of the parameters (γ,V, σ) are those that minimise in which case it can be shown that the first few Gauss-Hermite coefficients from (14) become h 0 = 1 and h 1 = h 2 = 0. In other words, if we choose (γ,V, σ) to be the parameters of the best-fit Gaussian to the LOSVD L(v), then the LOSVD can be written as with h 3 , h 4 , ... given by the integral (14). Conversely, if we adopt the parametrization (16) and fit (γ,V, σ, h 3 , h 4 , ...) simultaneously to L 0 (v), then we get back the same parameters we would obtain by first fitting (γ,V, σ) by minimising (15) and then using (14) to find the h i . For a strongly non-Gaussian L 0 (v) the parameters (γ,V, σ) obtained by minimising (15) need not be close to zeroth-, first-and second-order moments of L 0 (v), as shown on Figure 2. In particular, γ need not be close to 1. Gauss-Hermite fits to the LOSVDs of real galaxies, including B+05's measurements of M31, adopt the parametrization (16) and fit (γ,V, σ, h 3 , h 4 , ...). Unfortunately, the parameter γ is rarely reported, presumably because it is strongly affected by systematic effects in the fitting procedure, such as template mismatch, and because it does not affect the shape or width of the LOSVDs. Nevertheless, we note that it is an essential part of any Gauss-Hermite expansion. For now we assume that γ is known. Our iterative scheme for reconstructing it from the models is described in section 4.3 below.
Notice that equation (14) shows that the Gauss-Hermite coefficients h l can be thought of as modified moments: each h l is the integral over velocity space of some linear combination of the classical moments 1, w, w 2 , ..., w l , but weighted by the Gaussian factor α(w). Therefore, given a Gauss-Hermite fit (γ,V, σ, h 3 , h 4 ) to the line-of-sight velocity distribution at projected position (X,Y ), we treat (γ,V, σ) as being known perfectly and take the (luminosityweighted) modified moments (μ 0 obs , ...,μ 4 obs ) ≡ I × (1, 0, 0, h 3 , h 4 ) as our observables, where the surface brightness I(X,Y ) is obtained by convolving the WFPC2 image by the STIS PSF (1). We use equations (10) of vdMF93 to propagate B+05's quoted uncertainties (∆V i , ∆σ i , ∆h 3,i , ∆h 4,i ) to our observational errors (∆μ 0 i,obs , ..., ∆μ 4 i,obs ). Just as for the classical velocity moments, for each ring j = 1, 2, ... we use our Monte Carlo sample of 10 6 positions and velocities to evaluate the expressions for the modified moments, in which the rescaled velocity Then we convolve each of these distributions with the STIS PSF (1) and read off the values at (X,Y ) = (X i ,Y i ), giving the contribution of the j th ring to the k th modified moment of the i th STIS data point, µ k i,obs .

Modelling the effects of the bulge and P3
Our ring system is designed to model only the old red stars of the eccentric disc, but some of the light observed in the central few arcsec of M31 comes from other sources. The two main contaminants are M31's bulge and the compact young stellar cluster at P3. We follow Kormendy & Bender (1999) in modelling the surface brightness of the bulge as a Sersic profile I(R) = I 0 exp(−(R/R n ) 1/n ) with index n = 2.19, scale radius R n = 14. ′′ 0 and central V -band surface brightness I 0 = 15.40 mag. Our kinematic model for the bulge is very simple: it is non rotating and has a constant velocity dispersion σ bulge = 120 km s −1 . We add the contribution from this simple bulge model to our models' predictions for the WFPC2 photometry and the OASIS V and σ maps. We do not add it to the STIS predictions; we assume that B05 have successfully removed the bulge contribution from their STIS kinematics.
To model the contribution the young stars from P3 make to the V -band light we include a another component having surface brightness Σ(X,Y ) = Σ 3 exp(−R/R 0 ) in which R 2 = X 2 + Y 2 and the scale length R 0 = 0. ′′ 075 (L12). We follow B05 in assuming that the kinematics extracted from the red spectra are unaffected by the young stars; these stars affect only the WFPC2 photometry, not the OASIS or STIS kinematics.

Fitting the weights
For given BH mass M • and disc orientation (θ a , θ i , θ l ), the model's prediction for any observable O i can be be written as ∑ j P i j w j , in which w j is the weight given to the j th ring (equ. 7) and the matrix P i j gives the contribution that the j th ring makes to the i th observable, calculated using the method described in §3.3 above. An observable O i can be the light within a "superpixel" (given by WFPC2 photometry), a classical first-or second-order velocity moment (OASIS kinematics) or a 0 th ...4 th -order modified moment (STIS kinematics). We do not include the zeroth-order moments of the OASIS kinematics as these contain no additional information over the WFPC2 photometry.
Having a vector of observables (O 1 , O 2 , ....) and associated uncertainties (∆ 1 , ∆ 2 , . . .), we use a non-negative linear least squares algorithm (Lawson & Hanson 1974) to find the vector of non-negative weights w that minimises We then take this minimum value of χ 2 as a measure of the goodness of fit of the model with parameters (M • , θ l , θ i , θ a ).

RESULTS
Ideally, we would like to include all data sets (WFPC2, STIS, OA-SIS) in our set of O i and fit simultaneously to them. However, calculating the contribution the j th ring makes to a particular modified moment is computationally intensive and a full study of the parameter space spanned by M • , θ l , θ i and θ a that includes the STIS and OASIS data sets is not viable. We therefore conduct our modelling in two stages. We first fit models to the WFPC2 photometry and the OASIS velocity distribution in order to find the best set of orientation angles (θ a , θ i , θ l ). Then, having these angles, we fit to the WFPC2 photometry and STIS kinematics to obtain our estimates of the BH mass and the structure of the phase-space DF of the nucleus. As a final test of this model, we "observe" it at OASIS resolution and compare it (by eye) to the real M31. Before embarking on any of this model fitting, however, we first confirm that our models are indeed able to reproduce some of PT03's results.

A test: reproducing PT03's model
An immediate test of our modelling machinery is to reproduce some of PT03's results by using a sum of rings (8) that approximates one of their DFs. We focus here on their favoured nonaligned model, but we found comparable results with the poorerfitting model that is forced to be aligned with the main M31 disc. PT03's non-aligned models have DF (20) where σ e = 0.307. The parametric form for the backbone eccentricities is given by with α = 0.197 pc −1 , a e = 4.45 pc, a g = 1.71 pc and w = 1.52 pc. For the dispersion in inclination the form is where σ 0 I = 24. • 6 and a I = 31.5 pc. The function g(a) that sets the semi-major axis distribution is with a 0 = 1.37 pc, M • = 10.2 × 10 8 M ⊙ , c 1 = 4 pc −1 and c 2 = 4.24 pc. We treat the overall normalisation Σ 0 as a free parameter. We approximate this DF by a sum of 720 rings of the form (8) whose semi-major axes a j are spaced logarithmically in radius from 0.1 pc to 19 pc, with σ a, j = 0.8a j ∆ log a. We set σ e, j = 0.307 for all rings and set the mean eccentricity e j and dispersion in inclination σ I, j of the j th ring by evaluating (21) and (22) above at a = a j . The weights w j are set proportional to a j g(a j ), the factor of a j coming from the fact that our rings are equispaced in log a. Figure 3 shows the resulting model predictions for the WFPC2 photometry when viewed at the same angles as PT03's best-fit nonaligned model (θ a = −34. • 5, θ i = 54. • 1, θ l = −42.8. • 2) and including the contribution of the bulge model from sec. 3.4. Our reconstruction of the projected surface brightness in their model agrees closely with their figure 3: the model produces a nucleus that is broader than the observations and has an overly extended flat profile at P1.

Determining the orientation of the disc
We conduct an exhaustive scan over the space of orientation angles (θ a , θ i , θ l ) and black hole mass to obtain our best guess for the orientation of the disc. We first do this for a model that is forced to be aligned (θ i = 77. • 5) with the larger-scale M31 disc, before letting θ i vary freely. We include in these fits the 1123 OASIS V data points as well as a broad field spanning 5. ′′ 824 × 5. ′′ 824 from the WFPC2 data: this utilises our varied super-pixel scheme and contains 4096 data points. We do not include the OASIS σ maps in the fits, as it is dangerous to assume that the measured V 2 + σ 2 is a good estimate of the true second moment of the LOSVD ( figure 2).
For the aligned model we fix θ i = 77. • 5 and take θ a and θ l in the intervals [−90 • , 0 • ] and [−55 • , −25 • ] with spacing ∆θ = 2 • . A contour plot depicting the shape of the χ 2 goodness of fit for the inner region of this range is shown in figure 4. As expected the aligned model constrains θ l tightly but there is a weak dependence on θ a , with the best fitting model falling at θ a = −60 • and θ l = −35 • .
The non-aligned model looks at θ a from −58 • to −20 • , θ i from 50 • to 69 • and θ l from −39 • to −29 • , with ∆θ a = 2 • and ∆θ l = ∆θ i = 1 • . The best fit appears at θ a = −34 • , θ i = 57 • , θ l = −35 • for a black hole mass M • ≃ 1.25 × 10 8 M ⊙ . The shape of the χ 2 goodness of fit for the non-aligned model is shown in figure 5. The presented scan is for a single random distribution of stars projected with different angles: scans using distributions drawn from a different random seed found broadly similar results with slight variations ( figure 6). This presents issues with determining a single best fitting set of angles. Using a much larger (10 8 )  number of stars per disc in the region around the best fit angles informed our final selection of specific angles, however even then the random seed affected results. While in principle the angles could be determined more accurately, in practice shot noise from the model (and also the finite number of stars in the real nucleus!) limits the accuracy in each angle to the order of ≈ 1 • . Our final selection of angles was determined from the projection of the likelihood exp(− 1 2 χ 2 ) onto each of the θ a , θ i and θ l axes. We have experimented with allowing the bulge surface brightness I 0 and the central surface brightness Σ 3 of the young stellar disc described in sec. 3.4 to float our fitting procedure by including them as additional "weights" w j in the model and adding two ad-ditional columns to the projection matrix P i j , but we find that this makes little difference to our results.

Fitting WFPC2 photometry and STIS kinematics
Having the orientation angles we now drop the OASIS V maps and focus on using WFPC2 photometry together with the STIS LOSVDs to further constrain the model. The region of the WFPC2 photometry we use is restricted to an ellipse of semi-major axis 1. ′′ 6 and axis ratio 0.6 centred on P2. This ellipse is just large enough to encompass all STIS positions. We rebin Lauer et al. (1998)'s dithered image 4 by 4 into "superpixels" of side 0. ′′ 0455. Our vector of observables O i consists of the WFPC2 fluxes in all 2335 such superpixels that lie within the ellipse, together with the 5 modified moments (h i,0 , ...,h i,4 ) obtained from (γ i ,V i , σ i , h 3,i , h 4 i ) for each of the i = 1...22 LOSVDs measured by Bender et al. (2005) using the procedure described earlier in section 3.3.3.

Reconstruction of the γ(R) profile
The only complication in this is that we do not know the line strengths γ for any of our LOSVDs. We do, however, expect γ to be reasonably close to one and so we first fit a model in which all γ i = 1. Then, having the weights w j , we construct a realization of this model and "observe" it convolved with the STIS PSF (1). We fit Gauss-Hermite coefficients (γ,V, σ, h 3 , h 4 ) to the model LOSVDs at each point along the STIS slit, giving us a more informed estimate of how γ varies along the slit. We then repeat the whole fitting procedure, replacing our original γ i = 1 guesses with values read off from this reconstructed γ(R) distribution. We find that the resulting model converges after only a couple of iterations of this scheme. Figure 7 plots representative γ profiles obtained by this procedure for models with black hole masses M • /10 8 M ⊙ = 1.0, 0.8 and 1.2.   This shows that γ is significantly depressed close to the black hole where the LOSVDs are least well described by simple Gaussians.

Comparison of best-fit model against observations
The  Figure 8 shows the Gauss-Hermite coefficients of the reconstructed models along the STIS slit. The agreement with B05's observed kinematics is good, but there some features the best M • = 1.0 × 10 8 M ⊙ model cannot reproduce: the model does not fit the detailed shape of the V (R) and h 3 (R) profiles between R = −0. ′′ 9 and 0. ′′ 3 well and it predicts a central σ(R) that is slightly too high. Figure 9 shows how well the M • = 1.0 × 10 8 M ⊙ model fits to the WFPC2 photometry. Results for the other two black hole masses are similar. The agreement is good in the central regions, but beyond about 1 arcsec from P2 the model's surface brightness profile falls off too slowly compared to the observations. One possible explanation for this is that our model is simply too coarse; the thickness of each of the n a = 20 rings (equ. 8) in these models is σ a = 0.23a, which sets the models' characteristic radial spatial resolution. Another is that our model for the contribution of the bulge light (sec. 3.4) may be wrong within the innermost couple of arcsec.
Based on these comparisons, we interpret the relatively low values of χ 2 of our models not as a indication of the outstanding quality of our model fits, but instead as a sign that the treatment of the observational uncertainties -particularly of the WFPC2 photometry -could be improved. Nevertheless, we believe that the reader will agree that simple "chi-by-eye" tests indicate that our models produce the best fits to date of the M31 eccentric disc system. Of course, this is to be expected given that we have > 700 free parameters to play with, which is at least an order of magnitude more than most previous models of M31's nucleus.

DF of best-fitting model
What can we learn from all these free parameters, specifically the orbit weights w j ? Figure 10 shows two views of the DF f (a, e x , e y , I) of our best-fit model. This model has a strong negative eccentricity gradient between 0. ′′ 5 and 1. ′′ 2 (corresponding to P1), which is very similar to PT03's best-fit model; even the dispersion in eccentricity (figure 11) is similar to their value of 0.307. There are three important differences between our model and theirs, however.
(i) Beyond 1. ′′ 2 the mean eccentricity in our models becomes negative, meaning that the rings become mildly antialigned. We find that the strength of this feature depends on the details of our bulge model and so it is hard to judge its significance, but we note that just such a feature was predicted by Statler (1999) in his analysis of thin, self-gravitating discs.
(ii) Whereas PT03's parametrised model had an exponentially declining σ I profile, our model fits a σ I profile that increases with radius for a > 0. ′′ 15. This last point is qualitatively consistent with the predictions of collisional models of disc evolution (Stewart & Ida 2000;Peiris & Tremaine 2003). The detailed agreement is not so good though: whereas the collisional models predict σ I /σ e ≃ 0.5, our models fit σ I /σ e ≃ 2. It is not immediately clear, however, how far one can apply these calculations that assume almost circular e = 0 discs to the strongly eccentric disc in M31.
(iii) f (a) in our model falls steeply towards the centre from a =  This distinctive change in the distribution of the red stars is almost cospatial with the young, A-star population that make up P3. Figure 12 shows the face-on and edge-on projected density distributions of our best-fit model. Our machinery fits a density distribution which, when viewed face on, is broadly similar to the density distribution adopted by Salow & Statler (2004) in their self-gravitating razor-thin models. Apart from our neglect of self gravity, the most significant difference between our models and theirs is that ours have a secondary density peak around P2 and also have significant thickness.

Other projections of the best-fit model
Although the models we present in this section are not fit to the OASIS maps, we can nevertheless compare our reconstructed models' predictions against the real OASIS maps. Figure 13 shows the results for the M • = 1.0×10 8 M ⊙ model. The model reproduces the shape and orientation of velocity and dispersion maps very well. This provides an independent test of the orientation angles and the broad-brush features of the DF inferred by our models.

A brief experiment: retrograde orbits
Finally, we note that one proposed origin of the eccentric disc in M31 is from an m = 1 instability in an initially circular stellar disc that contains a small counter-rotating population of stars (Touma 2002;Jog & Combes 2009;Kazandjian & Touma 2012). N-body simulations of this instability (Kazandjian & Touma 2012) show that the dominant, prograde population becomes the eccentric disc, while the minor, retrograde population is puffed up into a strongly triaxial distribution. As a quick experiment to test whether this is easily detectable, we have tried doubling up our ring distribution by including in our models the retrograde ring corresponding to each of the 720 prograde rings considered above, giving a total of 1440 rings. This produces a noticeably better fit to the observations, with χ 2 = 784 instead of 812, with the fitting procedure picking out three main counter-rotating components (figure 14) at (a, e) = (0. ′′ 2, 0.8), (1. ′′ 2, 0.6) and (3, −0.1). None of these is easily identifiable with the more diffuse component predicted by Kazandjian & Touma (2012), although we suspect that that might be hard to detect with naive models such as ours. We interpret the innermost retrograde ring as suggesting that our ring-based decomposition fails at radii r < 0. ′′ 2, perhaps by not being fine enough. Similarly, we take the presence of the other two rings as a hint that our bulge model could perhaps be improved.

CONCLUSIONS
We have constructed an eccentric disc models of the nucleus of M31 by modelling it as a linear combination of fattened rings of stars orbiting in a purely Keplerian potential. Our models are an obvious generalisation of those of PT03 and -as expected -are noticeably (albeit perhaps not significantly) more successful at reproducing the features observed in the inner arcsec or so of M31.  One of the most fundamental parameters of any model of this system is its orientation. Like PT03, we assume that the disc has biaxial symmetry, but we find some differences in the Euler angles (θ a , θ i , θ l ) that specify the orientation of these symmetry axes. Our values of θ a = −34 • and θ i = 57 • agree reasonably well with their θ a = −34. • 5 and θ i = 54. • 51, but our value of θ l = −35 • (which directly controls the position angle on the sky) differs significantly from their θ l = −42. • 8, more than can be accounted for by the effects of shot noise in either the models or in the distribution of stars in the real disc.
Our orientation is also consistent with the N-body model of B01 (θ i = 55 • ± 5 • and θ l = −36 • ) and the models of Salow & Statler (2004) (θ i = 63. • 51 ± 10. • 80 and assumed θ l = −33. • 6) and Sambhus & Sridhar (2002) (θ i = 51. • 54 and θ l = 27. • 34), though it should be noted all these models are 2d which imposes constraints on their geometry; Sambhus & Sridhar (2002) obtained their orientation by de-projecting the photometry of the disc such that the outer isophotes became circular. We note that these values of θ i found for the old, red distribution of stars are close to the inclination of the young population in P3 (θ i = 55 • ± 2 • ) measured by B+05.
The most interesting result of our models is the DF they infer from the data, bearing in mind that they contain no prior "wisdom" about which DFs are dynamically plausible. The models suggest the presence of a distinct, compact disc of red stars within 0. ′′ 15 of the black hole. This disc is anti-aligned (with eccentricity e < 0) with respect to the larger-scale eccentric disc (which has e > 0). Outside this compact region the eccentricity distribution of our models is very similar to PT03's. The main difference between our models and theirs is in the inclination distribution: our models fit an RMS inclination profile σ I (a) that closely tracks the dispersion in eccentricity, with σ I (a)/σ e (a) ≃ 2. This is interesting: models of the collisional evolution of (circular) discs predict that this ratio should be ∼ 0.5.
Another difference between our models and PT03 is that we find evidence for an anti-aligned feature (e < 0) at a > 1. ′′ 5. How- Figure 14. The distribution of orbits in a model that includes both prograde and retrograde orbits. The top two panels show the DF f (a,e x ) (left) and the rms inclination (right) for the prograde stars. The bottom two panels show the corresponding distributions for the retrograde population. ever, the details of this feature depend on our assumed bulge model, which merits further investigation.
Our models prefer black hole masses M • of the order of 1.0 × 10 8 M ⊙ ; masses higher than about 1.2 × 10 8 M ⊙ are weakly ruled out. Although it would be possible to use our machinery to carry out a full scan of black hole masses and orientation angles, followed by a systematic investigation of the degeneracies in the DF, we believe that a more pressing task is to include the self gravity of the disc. The stars contribute a significant fraction (∼ 20%) of the mass of the BH+eccentric disc system, which means that it is dangerous to read too much into our present, purely Keplerian models. Past 2d (Bacon et al. 2001;Sambhus & Sridhar 2002;Salow & Statler 2004, e.g.,) models have shown it is possible to get plausibly good fits to the kinematics for large disc masses. The space of 3d disc distributions that project to yield the observed surface brightness profile is degenerate, but our orbital ring system serves as a good starting point for this investigation. Work on self-gravitating selfconsistent 3d disc models is now underway, and we hope that such models will provide further insight into the origin of the eccentric disc in M31 and possibly elsewhere (Lauer et al. 2005).