The rotation curve and mass distribution of M31

To gain a better understanding of the Andromeda galaxy M31 and its role in the Local Group, measuring its mass precisely is essential. In this work, we have constructed the rotation curve of M31 out to $\sim$125 kpc using 13,679 M31 objects obtained from various sources, including the LAMOST data release 9 (LAMOST DR9), the DESI survey, and relevant literature. We divide all objects in our sample into bulge, disk and halo components. For the sources in the M31 disk, we have measured their circular velocities by a kinematic model with asymmetric drift corrections. For the bulge and halo objects, we calculate their velocity dispersions and use the spherical and projected Jeans equation to obtain the circular velocities. Our findings indicate a nearly isotropic nature for the M31 bulge, while the halo exhibits tangential anisotropy. The results show that the rotation curve remains constant at $\sim$220 km s$^{-1}$ up to radius $\sim$25 kpc and gradually decreases to $\sim$170 km s$^{-1}$ further out. Based on the newly determined rotation curve, we have constructed a mass distribution model for M31. Our measurement of the M31 virial mass is $M_{\rm vir} = 1.14^{+0.51}_{-0.35} \times 10^{12} M_\odot$ within $r_{\rm vir} = 220 \pm 25$ kpc.


INTRODUCTION
Out of the billions of galaxies in the observable universe, only a few galaxies in the Local Group can be resolved into individual objects and studied in full phase space such as the positions, velocities and element abundances.The Andromeda galaxy M31, being the nearest spiral dynamical system to us and the only large spiral galaxy that can be studied in full and detail, serves as an ideal astrophysical laboratory for exploring the formation and evolution of galaxies.
Mass is a fundamental property of a galaxy.An accurate measurement of the mass of M31 is crucial for understanding its matter distribution, formation and evolution history, and the role played in the Local Group.So far, there are a number of efforts devoted to the determination of the dynamical mass of M31, with a variety of tracers such as emission-line objects (e.g.Rubin & Ford 1970;Kafle et al. 2018), globular clusters (e.g.Galleti et al. 2006;Lee et al. 2008;Veljanoski et al. 2014), tidal stellar streams (e.g.Ibata et al. 2004;Fardal et al. 2013;Dey et al. 2023) and satellite galaxies (e.g.Côté et al. 2000;Watkins et al. 2010;Tollerud et al. 2012); and a variety of methods such as the rotation curve method (e.g.Rubin & Ford 1970;Carignan et al. 2006;Chemin et al. 2009;Sofue 2015), velocity distribution fitting (e.g.Evans & Wilkinson 2000;Watkins et al. 2013), virial theorem (e.g.Hartwick & Sargent 1974;Tollerud et al. 2012), projected mass estimator and tracer mass estimator (e.g.Federici et al. 1990Federici et al. , 1993;;Lee et al. 2008;Veljanoski et al. 2013).
Even with these extensive efforts to determine the mass of M31, an accurate estimation still remains elusive.A summary of previous mass measurements is presented in Fig. 1, indicating significant un-★ E-mail: zhangxw@mail.ynu.edu.cn† E-mail: bchen@ynu.edu.cncertainties in many of the results due to unclear parameter values like the halo anisotropy parameter and limited data sample sizes.Intriguingly, it is still uncertain whether M31 or the Milky Way is the most massive member of the Local Group (e.g.Evans & Wilkinson 2000;Watkins et al. 2010).
Most of the mass measurements of M31 are based on the line-ofsight (LOS) velocity measurements of objects in it.Many projects have been undertaken to observe the spectra of objects in M31.Merrett et al. (2006, hereafter M06) published a catalogue of 3,300 emission-line objects with their positions, magnitudes, and velocities, which were discovered using the Planetary Nebula Spectrograph.Halliday et al. (2006) presented a catalogue of 723 planetary nebulae (PNe), utilizing the WYFFOS fiber spectrograph.Galleti et al. (2009) presented LOS velocities of 118 M31 globular clusters from a series of spectroscopic observations (Galleti et al. 2006(Galleti et al. , 2007(Galleti et al. , 2009)).Caldwell et al. (2008) used the Multi-Mirror Telescope (MMT) Hectospec to observe spectra of 670 cluster candidates.There are also other catalogues of emission-line objects (e.g.Sanders et al. 2012;Martin et al. 2018) and clusters (e.g.Veljanoski et al. 2013;Fan et al. 2011Fan et al. , 2012) ) in M31 that have been established.The Spectroscopic and Photometric Landscape of Andromeda's Stellar Halo (SPLASH; Guhathakurta et al. 2006) survey utilized the DEIMOS instrument on the Keck II Telescope to study the kinematic and dynamical properties of the disk and halo of M31.More than 5,000 red giants in 50 fields have been observed so far (Gilbert et al. 2018).Most recently, the Dark Energy Spectroscopic Instrument (DESI) survey has provided LOS velocities of 7,617 M31 sources (Dey et al. 2023).
We have attempted to compile the most comprehensive catalogue of sources in M31 with LOS velocity measurements by utilizing data from these projects.This promising sample will enable us to accurately determine the mass of M31.This paper is the first in a series of articles in which we will establish the rotation curve of M31.The rotation curve is a crucial tool for constraining the mass distribution of the galaxy, including its dark matter content (e.g.Rubin & Ford 1970;Carignan et al. 2006;Chemin et al. 2009).

DATA
We hereby introduce the data set used in the current work.The sample comprises objects primarily selected from two sources, M31 objects identified from the LAMOST data and existing ones from the literature.

Sources selected from the LAMOST data
In this work, we select objects in M31 from the LAMOST DR9 low resolution catalogue (Luo et al. 2015) released in March 2022.Based on 10 years of spectroscopic surveys, the LAMOST DR9 lowresolution catalogue contains 11.22 million spectra with a resolution of  ∼ 1800 and a wavelength range of 3800 -9000 Å.We briefly describe here how we selected the sources in M31 from LAMOST.For more details we refer the reader to Chen et al. (in prep.).
To identify objects in M31, we initially search for all sources within a sky area of ±15 • centered on M31 from the LAMOST DR9.This yields a sample of 648,485 spectra.To distinguish genuine M31 members, we first use a machine learning algorithm to classify these objects.We cross match the LAMOST sources with the SIMBAD database and catalogues from previous studies.We adopt objects with confirmed types as our training sample.Using the fluxes of the individual wavelengths as the input parameters, we train a random forest classifier to categorize the LAMOST objects into three classes: objects with emission lines, stars without emission lines, and star clusters.Spectra with significant emission lines and no obvious continuum are selected as candidates for emission-line nebulae.The nebulae candidates are then classified into various types such as PNe, H ii regions, and others based on their emission line flux ratios.This classification is carried out using the Baldwin-Phillips-Terlevich (BPT) diagram (Baldwin et al. 1981).Spectra with significant emission lines and continuum are classified as stars with emission lines.Using a straightforward kinematic analysis similar to that of Drout et al. (2009), we select M31 stars based on their positions and LOS velocities, both those with emission lines and those without.For star cluster candidates, we examine their morphology by analyzing their Pan-Andromeda Archaeological Survey (PAndAS; McConnachie et al. 2009) -and -band images.
As a result, the LAMOST M31 source catalogue includes a total of 1,019 sources, consisting of 258 emission-line nebulae, 382 clusters and candidates, and 379 supergiants and candidates.

Sources from the literature
The quantity of tracers has a greater impact on the accuracy and precision of mass measurement than the quality (Hughes et al. 2021).In recent years, the number of spectroscopic observed sources in M31 has significantly increased, thanks in large part to the development of wide-field multi-object spectrographs.To expand the sample size, we have also incorporated previous M31 source catalogues that have LOS velocity measurements in our work.The sources in these catalogues are primarily stars (Drout et al. 2009;Cordiner et al. 2010;Dey et al. 2023, referred to as D09, C10 and D23), emission-line objects (Merrett et al. 2006;Halliday et al. 2006;Sanders et al. 2012;Martin et al. 2018, referred to as M06, H06, S12, and M18), and clusters (Chen et al. 2016;Sakari et al. 2016;Veljanoski et al. 2014, referred to as C16, S16 and V14).Revised Bologna Catalogue Version 5 (RBCV5; Galleti et al. 2004) 1 and Studies of Resolved Objects in M31 2 by Caldwell & Romanowsky (2016) are also included (referred to as RBCV5 and Caldwell).We have compiled all published catalogues from the literature that we are aware of.The information pertaining to all collected catalogues from previous studies is presented in Table 1, including the facilities used in observation, fiber sizes, the wavelength coverage and resolution of the spectra, the types and size of samples.
To combine all these catalogues, both position error and LOS velocity error should be considered.We conduct cross-matching between catalogues, selecting a radius based on the tracer types and fiber sizes.We adopt a matching radius of 1.5 arcsec for stars and emission-line objects in M06, H06, S12, M18, D09, C10, D23, and Caldwell; whereas for clusters in C16, S16, V14, RBCV5, and Caldwell, we use a radius of 3 arcsec.If a source appears in multiple catalogues, we prioritize its LOS velocity measurement and type with smaller velocity uncertainty or better observing condition.It is possible that there are systematic differences between velocities from different catalogues.We perform cross-matching between literature catalogues and the LAMOST M31 source catalogue.We obtain zero-point values from the velocities of the common sources, which are then utilized to correct the velocities of all sources in the literature catalogues and align them with those from LAMOST.In The resulting combined catalogue comprises 13,679 M31 sources, all of which have corresponding LOS velocity measurements.Among these sources, 8,341 are stars, with the majority originating from D23.There are 3,997 emission-line objects, which consist of 3,274 PNe and candidates, 661 H ii regions and candidates, and 62 sources with types listed as 'other type' or 'unknown'.Additionally, there are 1,103 clusters, the majority of which (902) are globular clusters located in the halo of M31.The spatial distribution of all the catalogued sources is shown in Fig. 2.
In our sample, 10,353 sources are accompanied by LOS velocity errors, with 72% of them having errors ≤ 20 km s −1 , and 66% having errors ≤ 10 km s −1 .

Coordinate transformation
In this study, we adopt the coordinates of the M31 center as ( M31 ,  M31 ) = (10.6847083• , 41.26875 • ) (de Vaucouleurs 1958), the position angle and inclination angle of the M31 disk as 38 • (Kent 1989) and 77 • (Walterbos & Kennicutt 1987), respectively, and the distance of M31 as  M31 = 780 kpc (Holland 1998).To simplify analysis, we establish a Cartesian coordinate system projected in M31 disk, with the X-axis coinciding with the M31 major axis in the northeast direction and the Y-axis coinciding with the M31 minor axis in the northwest direction.
For each source, we are able to obtain its major axis projection distance () and minor axis projection distance ( ).However, determining its radii, which is the distance between the object and the center of M31, is challenging as the actual distance from the Sun is unknown.Therefore, we use different methods to calculate radii values for objects in the disk and in the bulge/halo of M31.For disk objects which are assumed to be concentrated on the disk plane, we use their projected distance in the disk.This gives us elliptical contours of radii ( = √︁  2 + ( /cos 77 • ) 2 ).For sources in the halo/bulge which are generally spherically distributed, we use their two-dimensional (2D) projected distance in the sky.This leads to round contours of radii (  = √  2 +  2 ).In the current work, we use a combination of these radii, denoted by R, as shown in Fig. 1.Furthermore, the symbol  is utilized to represent the three-dimensional (3D) radii.

Velocity Transformation
To transform velocity, we first convert the LOS velocity of each object into the Galactocentric coordinate system, and then into the M31centric coordinate system.This allows us to decompose the peculiar velocity ( pec ) of an object into three components, as: where  ★→⊙ is the motion of the object relative to the Sun, which is the LOS velocity that can be derived from its spectrum,  ⊙→MW is the motion of the Sun relative to the Galactic center, and  MW→M31 is the motion of the Galactic center relative to the center of M31.
For a given object, we first calculate its motion relative to the Galactic centre, i.e.,  ★→MW .In the current work, we adopt the Galactic rotation velocity at the position of the Sun as 239.3 km s −1 (McMillan 2011) and the Solar peculiar motion as ( ⊙ ,  ⊙ ,  ⊙ ) = (11.1,12.24, 7.25) km s −1 (Schönrich et al. 2010).We then have where  and  are longitude and latitude of the object in the Galactic coordinates.
We then compute the motion of the object relative to the center of M31, i.e.,  pec .To accomplish this, we decompose the velocity of M31 relative to the Galactic center in the Galactocentric frame into the Galactocentric radial velocity ( M31r ) and Galactocentric transverse velocity ( M31t ).We adopt ( M31r , where  is the angular separation between the object and the center of M31, and  is the position angle of the object with respect to the M31 center.

KINEMATICS MODELING: THE DISK OF M31
To properly discuss the kinematics of the disk of M31, it is crucial to first consider the range of the disk and the relevant tracers.The commonly accepted range where the disk dominates is between 1.2 and 9 kpc (Courteau et al. 2011).Rotation curves based on H i gas have been observed at much farther distances, such as 38 kpc (Chemin et al. 2009).Moreover, extensive kinematic surveys (e.g.Ibata et al. 2005) and wide-field imaging surveys (e.g.Bernard et al. 2015) have shown that M31 has an extended disk, spanning radii from 15 kpc to 40 kpc and exhibiting extended disk-like substructures at around 60 kpc.In our study, we have explored various disk limits, favoring a more conservative approach to minimize contamination from halo objects and obtain a relatively pure sample of the disk.We plot the radial velocity dispersion values for objects located close to the minor axis of M31, using radii limits of 20, 25, and 30 kpc, respectively, in the upper panel of Fig. 3. Within the radius bins of 20 kpc, we observe a decrease in the radial velocity dispersion values as the radius increases.However, when we consider disk limits of 30 and 25 kpc, we observe an increase in radial velocity dispersion beyond 20 kpc, indicating contaminations from halo objects.Among the objects in our sample, 89% of the H ii regions and candidates, which are believed to be young disk objects, are located within the range of 1.2-20 kpc.H ii regions and candidates beyond this range mainly come from the M06 catalogue and are considered candidates rather than confirmed H ii regions.Therefore, we adopt a maximum radius limit of 20 kpc for the kinematic analysis of the M31 disk.For our analysis of the M31 disk, we have focused on tracers such as stars (including OB stars, supergiants, and other types) and emission-line objects (including PNe, H ii regions, and unclassified objects) located within the selected disk radii range, i.e., from 1.2 to 20 kpc.For the subsequent kinematic analysis, we utilize the cylindrical coordinate system centered on the centre of M31, which is similar to the one in Section 2.3.1 but with slight differences.For disk objects, we disregard their vertical velocity component in the -direction, which is perpendicular to the M31 disk.We measure the angle coordinates  in a counterclockwise direction from the negative semi-axis of the X-axis and the radial coordinate  increasing outward.The positive direction of the tangential velocity   is opposite to the direction of  increase.By adopting this convention, we ensure that the disk sources have positive circular velocities.The positive direction of the radial velocity   points away from the center of M31.We can then decompose  pec into two components:   and   , given by, where  is the inclination angle of the disk of M31.
The tangential velocity   is not equivalent to the circular velocity   due to the non-negligible effect of asymmetric drift (  ), which arises from the heating effect of the disk.We assume the disk is in a steady state and ignore the variations of covariance     with Z.The density of tracers in the disk and their radial velocity dispersion,   , both decrease exponentially with , and the scale lengths for the decrease in number density and radial velocity dispersion are   and   , respectively.The asymmetric drift is then given by (e.g.Binney & Tremaine 2008) where   is the tangential velocity dispersions.In this work we adopt a value of   = 5.3 ± 0.5 kpc based on the work of (Courteau et al. 2011).For the rest of the parameters, we will derive them from our sample as some are not available in the literature at present, and some are accompanied by considerable uncertainty.

Radial and tangential velocity dispersions
The radial velocity dispersion   () of the disk source in M31 is still uncertain.It is very difficult to measure the radial velocity   of an M31 disk source directly.However, we can approximate the LOS velocity and dispersion of the source located on the minor axis of the M31 disk as a radial velocity.We select objects located in a fan-shaped region near the minor axis of the M31 disk, as given by: |/| < 0.1.For these sources, their tangential velocity   is perpendicular to our LOS, so the contribution of the LOS velocity to   can be neglected.The value of the radial velocity   of these objects can be calculated by, The selected tracers are divided into seven radius bins, each containing more than 35 sources.The mean and dispersion uncertainties of the individual bins are estimated from 300 bootstrapping samples.
The top panel of Fig. 3 show the dispersion values for these seven bins.
As in the case of the Galactic disk (e.g.Huang et al. 2016), we assume that the radial velocity dispersions   of the objects in the disk of M31 decrease exponentially with respect to the radii , given by, where  LOS 0 is the radial velocity dispersion at the reference position  0 .In the current work, we adopt  0 = 10 kpc considering the size of the disk of M31.We use Monte Carlo Markov Chain (MCMC) analysis to extract the best-fit values for   0 and   .As a result, we obtain:   0 = 62.7 +3.4 −3.4 km s −1 and   = 19.8+6.5 −4.0 kpc.The best-fit curve is illustrated in top panel of Fig. 3.
In a similar manner, for the tangential velocity   , we select objects that are located close to the major axis of the disk of M31, where the radial velocity   is perpendicular to our LOS.Only objects that meet the criterion | /| < 0.07 are chosen, and their tangential velocity   can be calculated by, We group the selected objects into ten radius bins, with tracer size of each bin larger than 45.The bottom panel of Fig. 3 presents the tangential velocity dispersion values for each bin.We have compared our results with those obtained from previous studies, including M06, H06, and Zou et al. (2011), which are over-plotted in the Fig. 3.The results demonstrate a good agreement with the previous studies.We used an exponential function similar to equation ( 7) to model the tangential velocity dispersion, given by, where   0 is the tangential velocity dispersion at  =  0 and   the scale length of the tangential velocity dispersion.Similarly, we have applied the MCMC analysis to estimate the best-fit parameters, which yields values of   0 = 69.6 +3.3 −3.3 km s −1 and   = 21.0 +5.2 −3.5 kpc.The corresponding best-fit curve is shown in bottom panel of Fig. 3.

Tangential velocity and asymmetric drift correction
In order to determine the circular velocity of objects in M31 disk, we need to make the asymmetric drift corrections to their tangential velocity.We previously computed the tangential velocities of a small subset of the disk objects (508 out of 5261 tracers) that are located near the major axis of the M31 disk.Now, we aim to estimate the tangential velocities of all the selected disk objects.To do so, we assume a radial velocity function for all disk objects.Ideally, the radial velocities of objects in a stable, symmetric disk would be zero.However, as depicted in upper panel of Fig. 4, the mean radial motion   oscillates with , much like what has been observed in the Milky Way (Williams et al. 2013;Huang et al. 2016;Wang et al. 2022).The mean radial velocity varies around zero, indicating either an elliptical mean orbit of the disk objects or the process of radial migration.To smooth out this oscillation relation, we apply a cubic polynomial fit, which is also shown in Fig. 4. Based on the fitted relation, we calculate the radial velocity of all the disk objects, then the tangential velocity of each object can be derived using equation ( 4).The distribution of tangential velocities calculated for all disk sources in M31 is shown in the lower panel of Fig. 4. In the diagram we also show the mean values of the tangential velocities for the objects located close to the major axis of the disk of M31.The tangential velocities vary with radius in a manner very similar to that observed for sources located near the major axis of M31, demonstrating the robustness of our results.
With the tangential velocities of all the selected disk objects and all parameters except for the circular velocity to calculate the asymmetric drift corrections, we can now determine the circular velocities of the individual objects in the disk of M31.This is achieved by combining equation ( 5) with the relation   =   +   .Tracers located near the minor axis (|/| < 0.1) are excluded from the analysis, as their peculiar velocities have minimal contribution to the circular velocities, leading to a larger uncertainty.To account for the errors we use the Monte Carlo (MC) method.

KINEMATICS MODELLING: THE BULGE AND HALO OF M31
We first select objects that belong to the bulge and halo of M31.
According to Section 3, we set the radii range of the bulge as  < 1.2 kpc, and that of the halo as  > 30 kpc.Halo stars from DESI observations, PNe, unclassified emission-line objects, and clusters that fall within these radius ranges are selected.Due to the different methods used to calculate the radii of disk and halo objects, there is overlap between the disk and the halo in the resulting rotation curve (see Fig. 6).This overlap will be utilized to constrain the halo anisotropy parameter.In the subsequent calculations of this section, we will employ the 3D radius () and the 2D projected radius (  ) in the sky (see Section 2.3.1).Unlike the disk, both the bulge and the halo of M31 are assumed to be spherical gravitational potential systems.To model the kinematics of these systems, we use the Jeans equation given by (Jeans 1919;Binney & Tremaine 2008), where  is the number density,   is the radial velocity dispersion,  is the anisotropy parameter, and Φ is the total gravitational potential.The circular velocity is associated with the gravitational potential by, The number density  of the bulge and halo of M31, as utilized in this study, is based on the results of Courteau et al. (2011).The bulge is characterized by a projected profile that can be described by the Sérsic form, with a Sérsic shape index of  = 2.2 ± 0.3, a Sérsic shape parameter   = 1.9992 − 0.3271, and an effective radius   = 1.0 ± 0.2 kpc.On the other hand, the halo exhibits a power law profile in both 2D and 3D.This can be expressed as a 2D surface density profile, Σ = Σ 0   2D  , and a 3D density profile,  =  0   3D , where Σ 0 and  0 represent the 2D and 3D densities at 1 kpc, respectively.The power law indices are chosen as  2D = −2.5 ± 0.2 and  3D = −3.5 ± 0.2 (Courteau et al. 2011).
Following the works of Veljanoski et al. (2014) and Gilbert et al. (2018), we adopt a power law for the radial velocity dispersion,   , such that   ∝   .Now the circular velocities of the bulge and halo are respectively given by, and However, since we are unable to directly measure the actual distance and 3D velocity of individual objects in the M31 halo, we can only obtain projected 2D kinematic parameters from observations.In this study, we provide the deprojected circular velocities of the halo as a replacement for equation ( 13), given by, where  0 Σ 0 = 0.57 kpc −1 is the ratio of number density to surface density at 1 kpc.The derivation of this equation is described in Appendix A. Previous literature does not offer comprehensive 3D information on the circular velocities of the bulge, making it difficult to establish a corresponding equation.Therefore, we continue to use equation ( 12) with 2D parameters to approximate the 3D circular velocities of the bulge.Within the bulge region, there are only two available measurements of circular velocity, and their impact on the virial mass is minimal.Furthermore, the circular velocity values in the bulge are primarily influenced by  bulge .Instead of opting for more complex models that incorporate missing information, substituting the 2D scenario is a preferable choice.
To calculate  LOS for the halo of M31, we divide the halo objects into 20 radial bins.In the case of the inner 15 bins, as the radius increases, the sky area covered by each radial bin widens progressively, resulting in additional dispersion beyond the intrinsic value.To address this, we partition the objects within each bin into four quadrants using the  and  axes of M31.Subsequently, we calculate the LOS velocity dispersion for each quadrant and derive a weighted average as the velocity dispersion value for the bin, considering the number of objects in each quadrant.For the outer 5 bins, we do not differentiate between quadrants due to the limited number of objects in our sample and calculate the LOS velocity dispersion directly.Fig. 5 displays  LOS of the halo objects as a function of   , with previous results from Gilbert et al. (2018) over-plotted.Our results are generally consistent with Gilbert et al. (2018), though slightly higher.We fit  LOS as a power law (Veljanoski et al. 2014;Gilbert et al. 2018), given by, where we keep   0 fixed at 30 kpc and use MCMC analysis to determine the parameters  LOS 0 and .Our analysis yields  LOS 0 = 107.88+1.77 −1.74 km s −1 and  = −0.09± 0.03.The best-fit curve is displayed in Fig. 5.
To obtain the  LOS of the bulge, we separate the bulge objects into two radius bins and compute the radii and LOS velocity dispersions for each bin.However, bulge velocity dispersion is still assumed in power law form as equation ( 15).As there are only two  LOS measurements within 0 <  < 1.2 kpc for the bulge, we have not attempted to fit a velocity dispersion relation. = −0.09± 0.03 obtained from the halo fits is adopted for both the bulge and halo analysis.
Finally, given the lack of knowledge about the anisotropic parameter of the bulge and halo of M31 ( bulge and  halo ), we treat them as free kinematic parameters in conjunction with the dynamical parameters determined in Section 6.However, as observed in the Milky Way, the anisotropy parameters exhibit spatial variations.For the bulge and halo of M31, accounting for the influence of radius on  presents challenges, so we only constrain an spatially averaged value for these parameters.

THE ROTATION CURVE OF M31
Using the most extensive collection of velocity measurements of M31 objects to date, we have generated a rotation curve for M31 that encompasses all of its components, including the bulge, disk, and halo.Hereafter, we adopt the terminology of "rotation velocity" and "rotation curve" rather than "circular velocity" to maintain consistency with prior studies.However, it is important to note that these concepts are only equivalent in the context of a cold disk.The rotation velocity of stars is typically lower than the circular velocity, and in some literature, it is even assumed to be effectively zero in the halo.The resulting rotation curve extends out to a distance of ∼125 kpc, which is presented in Fig. 6.Table 2 displays the values of our calculated rotation curve.Table 3 presents a summary of all the kinematic parameters utilized in the calculation of the rotation curve for the current study.

Rotation curve of the disk of M31
The resulting rotation curve of the disk of M31, along with a comparison to previous H i studies (Carignan et al. 2006;Chemin et al. 2009), is shown in Fig. 7.The disk objects are divided into 20 radius bins, each containing over 220 objects.The typical rotation velocity in the disk of M31 is about 220 km s −1 .Our results generally agree with the previous rotation curves obtained from H i observations.However, for regions at radius  = 7 -18 kpc, our estimates for the rotation curve are slightly lower than those from the H i studies.Since the difference is relatively small and still falls within the 3 range, it is probable that the dissimilarity is caused by systematic differences in the tracers and methods used to measure the rotation curve.
Our calculated rotation velocities in the disk of M31 have a typical uncertainty of about 17 km s −1 .The most significant contributor to this uncertainty is the tangential velocity error.Since we cannot directly measure the tangential velocities of the individual objects in the disk, we have assumed a radial velocity and radius relationship to derive the tangential velocities.Fig. 4 shows that the uncertainties of the fitted radial velocities can range from 9 to 30 km s −1 , leading to an error of around 5 to 38 km s −1 in the obtained rotation velocity.Typically, this results in an error of approximately 16 km s −1 .Aside from this, there are other sources of uncertainty, including errors in the parameters used.For instance, the scale length of number density   can have a maximum value of 6.5 kpc (Kormendy & Bender 1999).Using a different value of   would result in a systematic rotation velocity error of about 5 km s −1 .Finally, the uncertainty of rotation curve caused by the error of the scale length of radial velocity dispersion   is about 2 km s −1 , while that from the error of the scale length of tangential velocity dispersion   is about 1 km s −1 .

Rotation curve of the bulge and halo of M31
The objects in the bulge of M31 has been divided into two radius bins.The rotation velocities are 163.24 and 210.95 km s −1 at distances of 0.22 and 0.65 kpc, respectively.These values are lower than those of the disk.
The uncertainty of the rotation velocities of the bulge is about 33 km s −1 .This can be attributed to our lack of understanding of the bulge parameters, which themselves have large errors.The Sérsic effective radius of the bulge number density distribution   has a maximum reported value of 1.93 by Seigar et al. (2008), which could result in a systematic error of approximately 30 km s −1 for the rotation velocity.measurement of the Sérsic shape index of the bulge number density  has an uncertainty that could lead to a systematic error of approximately 12 km s −1 for the rotation velocity.Additionally, for the velocity dispersion power law index , Gilbert et al. (2018) has measured a value of −0.12, which could result in a systematic error of about 5 km s −1 for the rotation velocity.
Rotation velocities are calculated in 20 radius bins for the halo of M31.Since the way for measuring radii in the disk differs from that used for objects in the halo (as described in Section 2.3.1),there are some overlapped bins in the rotation curve.In these overlapped radii, the rotation curve of the disk agrees well as that of the halo, which strengthens the reliability of our results.The rotation curve of the halo exhibits a radial decline, decreasing from around 230 km s −1 at a distance of 30 kpc to approximately 170 km s −1 at 125 kpc.Equation ( 14) demonstrates high sensitivity to kinematic parameters, making it crucial to obtain precise values.Error propagation analysis reveals that a mere 0.2 error in  2D and  3D can introduce a substantial error of 30 km s −1 into the rotation curve.The surface density power law index of the halo, denoted as  2D and documented in Courteau et al. (2011), exhibits a considerable range from −3.5 to −2.2.This wide variation in  2D significantly impacts the halo's rotation curve.Additionally, the power law index , as reported in Gilbert et al. (2018), is associated with an error of approximately 18 km s −1 in the rotation velocities.Fortunately, due to the degeneracy between these kinematic parameters and  halo , the uncertainty in halo rotation velocity is somewhat mitigated.However, it remains imperative to achieve more precise measurements of these kinematic parameters, particularly .

The gravitational potential model
We have developed a parameterized potential model to fit the rotation curve of M31 and obtain its mass distribution.The model consists of three components: the bulge, disk, and halo, and the total gravitational potential is obtained by summing these components.
In the current work, we adopt the potential models proposed by Hernquist (1990) and Miyamoto & Nagai (1975) for the bulge and disk, respectively.Therefore, the bulge and disk potentials can be respectively given by, and where  is the gravitational constant,  b and  d are the mass of the bulge and disk, respectively, and the parameters , , and  correspond to the scale lengths.If  = 0, the potential of the disk (Φ disk ) simplifies to Plummer's spherical potential (Plummer 1911), while if  = 0, it simplifies to Kuzmin's potential for a razor-thin disk (Binney & Tremaine 2008).The ratio / characterizes the flatness of the M31 disk.Depending on the values of  and , the potential model can range from a nearly infinitesimally thin disk to a spherical system.
For the dark matter halo of M31, we adopt the NFW model (Navarro et al. 1996), given by, where  vir is the virial radius,  is concentration parameter and () = ln(1 + ) −  1+ . vir is the virial mass given by the equation  vir = 4  3  3 vir Δ  , where Δ is the virial overdensity parameter, and   = 3 2 0 8  is the critical density of the universe, where  0 = 70 km s −1 Mpc −1 .We set Δ = 200, which is commonly used in the definition of halo mass in simulations and observations (Seigar  . et al. 2012;Fardal et al. 2013;Veljanoski et al. 2014).By combining equation ( 16), ( 17), and (18), we can create a model for the total potential of M31.This potential model has seven free parameters.In addition to these parameters, we also have the bulge and halo anisotropy parameters (see Section 4), which brings the total number of free parameters to nine.We denote these parameters as Θ = {, log( b ), , , log( d ), log(),  vir ,  bulge ,  halo }.We use a Bayesian approach to model the posterior probability distribution of these free parameters, as given by, (Θ|R,   ,  err ,  LOS ,  LOSerr ) ∝ (Θ) L, where  err is the error of rotation velocity,  LOSerr is the error of LOS velocity dispersion and (Θ) is the prior.We adopt flat priors with the following ranges: (), (), and () are uniformly distributed between 0 and 10 kpc, while [log( b /M ⊙ )] and  [log( d /M ⊙ )] are uniformly distributed between 9 and 12.We also constrain log( b /M ⊙ ) > log( d /M ⊙ ).We consider different previous estimates of the concentration parameters of the Milky Way and M31, as reported in previous studies (e.g.Kafle et al. 2014;Huang et al. 2016;Kafle et al. 2018), and set  [log()] to be uniformly distributed between −2 and 10.Additionally, we set  [log( vir /M ⊙ )] to be uniformly distributed between 10 and 14, so the distribution of  vir /10 12 M ⊙ is  [ =  vir /10 12 M ⊙ ] = 1/(4 ln 10 • ) between 0.01 and 100.L is likelihood, given by, for disk, and for bulge and halo, where N is a Gaussian function, Φ(R) is the potential.Model(R) is given by equation 11, while the bulge and the halo   is given by equation 12 and equation 14.We employ an MCMC analysis to sample prior distribution of these free parameters using the emcee python package (Foreman-Mackey et al. 2013).40 walkers are employed and we run the chains for 10,000 steps to ensure convergence, with the first 3,000 steps are discarded as burn-in.

The best-fitted potential and mass distribution
Fig. 8 displays the marginalized (1D) and joint (2D) probability density functions (PDFs) for the model parameters in both one-and two-dimensional (1D and 2D) forms.The joint probability density reveal correlations among some of the parameters, such as a strong correlation between  b and  as well as anticorrelations between  d and .Other pairwise parameters are generally independent.The best-fit values for the model parameters are estimated using the median values of their marginalized PDFs, while uncertainties are calculated from the 68% probability intervals of the marginalized PDFs.It should be noted that the rotation curve of the disk constrains the rotation of the bulge and halo, effectively determining the values of  bulge and  halo .As a result,  bulge and  halo have independent joint distributions with each other and with  vir .However, this does not mean that we can eliminate the degeneracy between mass and anisotropy, as there are uncertainties in the disk rotation curve.For the same reason,  bulge and  halo in the fitting are still recommended values.
According to our analysis, the bulge of M31 has a anisotropy parameter of  bulge = 0.07 +0.13  −0.29 .Other kinematic parameters would influence value of  bulge , making it varying between −0.12 and 0.30.The narrow down process for the uncertainty of  bulge remains challenging.Nevertheless, our calculation leads us to an approximate value of 0, suggesting that M31 likely possesses an isotropic bulge.Regarding the halo, our conclusion indicates a tangential anisotropy parameter, with  halo = −0.54+0.38  −1.29 .The MCMC sampling probability distribution of  bulge displays a prolonged negative tail and a corresponding large lower error, indicating a diminishing slope of rotation velocity associated with decreasing  bulge .
In our gravitational potential model, the M31 galaxy is divided into the bulge, the disk, and the halo.However, recent studies have suggested that M31 is a barred spiral galaxy (e.g.Athanassoula & Beaton 2006;Beaton et al. 2007;Blaña Díaz et al. 2016;Blaña Díaz et al. 2018).Objects associated with the bar do not follow perfectly circular orbits, which can be misleading when calculating the rotation curve.This can result in a noticeable peak in the rotation curve of the H i gas at around 1-2 kpc, as shown in Fig. 7.It is worth noting that our analysis did not reveal similar structures.The investigations and studies on the bar of M31 primarily focus on the observations in the infrared bands, whereas this particular study employs stars, emission line objects, and star clusters as tracers, utilizing their optical wavelengths spectra.This distinction in observational methods may be a potential factor contributing to the absence of bar structures.Additionally, the rotation velocities of the bulge are determined based on velocity dispersion measurements, rather than being directly derived

SUMMARY
We have obtained a large sample of M31 sources with LOS velocity measurements, consisting of a total of 13,679 objects selected from LAMOST DR9 and the literature.With the largest sample available to date, we have explored the kinematic properties of M31.We have obtained measurements for the radial and tangential velocity dispersion profiles in the disk of M31.Specifically, at a projected radius of  = 10 kpc, our measurements indicate   0 = 62.7 +3.4  −3.4 km s −1 and   0 = 69.6 +3.3 −3.3 km s −1 .We also provide the constraints on the scale lengths of both the radial and tangential velocity dispersion, with   = 19.8+6.5  −4.0 kpc and   = 21.0 +5.2 −3.5 kpc.Our measurements suggest that the velocity dispersion of M31 halo stars exhibits a mild decrease with projected radius, with a power-law index of  = −0.09± 0.03.We also reveal evidence of a isotropic bulge, with a recommended value  bulge = 0.07 +0.13  −0.29 , and tangential anisotropy in the M31 halo, with a recommended value  halo = −0.54+0.38  −1.29 .Based on these derived kinematic properties, together with parameters from existing literature, we have derived the rotation curve of M31 from kinematic models.The rotation curve which encompasses all of the components of M31, including the bulge, disk, and halo, extends up to a distance of about 125 kpc.The resulting rotation curve of M31 is generally flat at 220 km s −1 in the disk, gradually decreasing to 170 km s −1 at ∼125 kpc.We have also provided a detailed discussion of the sources of error in our rotation curve measurements, clearly identifying each source of error.Based on the newly derived rotation curve of M31, we have obtained the best estimate of its virial mass to be  vir = 1.14 +0.51  −0.35 × 10 12  ⊙ within  vir = 220 ± 25 kpc.Our measurement is consistent with the results from earlier rotation curve and escape velocity curve analyses.
However, our understanding of many of the kinematic parameters of M31 remains limited.As large scale multi-object spectroscopic surveys like LAMOST, SDSS, and DESI continue, and new projects like China Space Station Telescope (CSST), WEAVE (William Herschel Telescope Enhanced Area Velocity Explorer Instrument) , and Prime Focus Spectrograph (PFS) in Subaru to be initiated, we anticipate that we will gain a deeper understanding of our neighboring galaxy, much in the same way that we have with the Milky Way.

Figure 1 .
Figure 1.A summary of recent M31 mass measurements from previous works.The cumulative mass profile calculated in this study is represented by the red solid line with grey shade indicating uncertainties.

Figure 2 .
Figure2.The spatial distribution of all sources in our final catalogue.The upper left panel illustrates the distribution of stars in the disk of M31, which includes OB stars (red dots), supergiants (blue dots), and other types (gray dots).The upper right panel depicts the distribution of emission-line objects, including PNe (blue dots), H ii regions (red dots), and unclassified objects (gray dots).The bottom panel shows the distribution of clusters (red dots) and stars in the halo of M31 (gray dots).The black dashed ellipses represent circles with a disk projected radius  = 20 kpc, while the green and purple dashed circles represent sky projected radii of   = 50 kpc and 100 kpc, respectively.

Figure 3 .
Figure3.Dispersions of the radial velocity (upper panel) for the selected objects close to the minor axis and the tangential velocity (bottom panel) for the selected objects close to the major axis.The measured dispersion values for individual radius bins (using 20 kpc as the disk limit) are represented by black dots with error bars.The best-fit curves are shown as black lines, and the 1 regions of the best fits are depicted with light grey shadings.In the upper panel, additional blue and red dots with error bars represent disk limits of 25 and 30,kpc.The black dashed curve, black dots and green dots in the bottom panel show the results from M06, H06, andZou et al. (2011), respectively.

)Figure 4 .
Figure 4. Upper panel:Mean values of the radial velocities for the selected objects close to the minor axis (blue dots).The red line and light grey shading represent the best cubic polynomial fit curve and the 1 region error respectively, used to infer the radial velocity of all the disk objects.Bottom panel: Mean values of the tangential velocities for the selected objects close to the major axis (blue dots) and the density distribution of estimated tangential velocities for all the disk sources of M31 (grey scales).The darker the color, the greater the density.

Figure 5 .
Figure5.The relationship between  LOS and   for the halo of M31.The blue dots correspond to the  LOS values for each of the radius bins.The red curve and the grey shading represent the best-fit curve and the 1 region of the fit, respectively.The black dashed line and black points with errorbars denotes the velocity dispersion fromGilbert et al. (2018).

Figure 6 .Figure 7 .
Figure 6.Final combined 3-part rotation curve of M31, comprising the disk (blue dots), bulge and halo (red dots).The different line styles, as indicated in the top-left corner of the diagram, represent the best-fit rotation curve contributions to the individual components of M31.The black line represents the sum of contributions from all the mass components, while the grey shading corresponds to the 1 region of the fit.

Table 1 .
Information of our collected literature catalogues.

Table 2 .
The rotation curve of M31.

Table 3 .
Kinematic parameters used to calculate the rotation curve.