# Dark matter fraction derived from the M31 rotation curve

F Hammer, Y.B Yang, P Amram, L Chemin, G.A Mamon, J.L Wang, I Akib, Y.J Jiao, H.F Wang

► To cite this version:

F Hammer, Y.B Yang, P Amram, L Chemin, G.A Mamon, et al. Dark matter fraction derived from the M31 rotation curve. Astronomy & Astrophysics - A&A, 2025, 694, pp.A16. 10.1051/0004-6361/202452753. hal-04847844

HAL Id: hal-04847844

https://hal.science/hal-04847844v1

Submitted on 16 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Dark matter fraction derived from the M31 rotation curve

F. Hammer<fnref n="1" />,*, Y. B. Yang<fnref n="1" />, P. Amram<fnref n="2" />, L. Chemin<fnref n="3" />, G. A. Mamon<fnref n="4" />, J. L. Wang<fnref n="1" />, I. Akib<fnref n="1" />, Y. J. Jiao<fnref n="1" />, and H. F. Wang<fnref n="5" />

Received 25 October 2024 / Accepted 30 November 2024

## ABSTRACT

Mass estimates of a spiral galaxy derived from its rotation curve must account for the galaxy's past accretion history. There are several lines of evidence indicating that M31 experienced a major merger 2 to 3 Gyr ago. In this work, we generated a dynamical model of M31 as a merger remnant that reproduces most of its properties, including from the central bar to the outskirts. The model accounts for M31's past major merger and reproduces the details of its rotation curve, including its 14 kpc bump and the observed increase of velocity beyond 25 kpc. We find non-equilibrium and oscillatory motions in the gas of the merger-remnant outskirts caused by material in a tidal tail returning to the merger remnant. A total dynamical M31 mass of $4.5 \times 10^{11} M_{\odot}$ within 137 kpc was obtained after scaling it to the observed HI rotation curve. Within this radial distance, we find that 68% of the total dynamical mass is dark.

**Key words.** galaxies: halos – galaxies: interactions – galaxies: kinematics and dynamics – Local Group – dark matter

## 1. Introduction

The Andromeda galaxy (M31) is our nearest giant spiral neighbor. It has been considered as a twin of our Galaxy for decades, and as such, it offers a precious external view to our home. M31 has played a considerable role in establishing the modern cosmology, including in determining the difference between galaxies and Galactic nebulae as well as the presence of dark matter (DM) in the galaxy halo. The latter discovery took advantage of the very large sky angle displayed by the M31 disk, which is seen almost edge-on, allowing for easy measurement of its rotation from radial velocities. Lundmark (1925) was one of the first to measure the rotation curve (RC) of M31, followed by Babcock (1939), and then Mayall (1951). They had already observed that it remained flat up to almost 20 kpc. Later, Rubin et al. (1978) and Bosma (1978) showed that several spiral galaxies had flat RCs, including up to 35 kpc for M31.

Considerable theoretical and observational progress has been accomplished since the late 1970s. Formerly, it was believed that disk galaxies were formed by an initial collapse (Eggen et al. 1962). Because thin disks cannot support significant merger events without being destroyed (Toth & Ostriker 1992), the common consensus was that giant galactic disks have a stable rotation since very early epochs. This conjecture faced the problem of the angular momentum necessary to support disk galaxies (Navarro & Steinmetz 2000), that is, giant disks generated from cosmological fluctuations were too small. However, accurate comparisons between distant and nearby galaxies (Hammer et al. 2009; Delgado-Serrano et al. 2010) indicated that about 50% of spiral galaxy progenitors were in a merger phase 6 Gyr ago. This can be understood because gas-rich

mergers are very efficient at reforming a disk in the remnant (Hopkins et al. 2009), which also solves the disk angular momentum acquisition problem. This especially concerns M31, for which the accumulation of observations have revealed the giant Stream in its halo (Ibata et al. 2001), the past 2–3 Gyr star formation events (Williams et al. 2017), the stationary 10 kpc ring (Dalcanton et al. 2012; Lewis et al. 2015), and the unexpectedly large velocity dispersion of stars older than 2 Gyr (Dorman et al. 2015). All of these exceptional features have been reproduced by a gas-rich 4:1 merger that has occurred 2 to 3 Gyr ago (Hammer et al. 2018, hereafter H18), which confirms the conjecture that the Milky Way and M31 are not twin galaxies because they have very different histories and dynamical properties (Hammer et al. 2007). The robustness of the M31 modeling has been confirmed by its ability to reproduce different observational features discovered later (Bhattacharya et al. 2021, 2023; Escala et al. 2022; Tsakonas et al. 2024).

Having a recent merger in M31 implies a disk in which stars and gas have only recently reached their present orbits, which may warrant disk equilibrium in the inner regions but possibly not in the outskirts. The goal of this paper is to identify the DM content of M31 that is expected from the M31 RC in the context of a recent major merger. Section 2 compares the different M31 RCs extracted from the HI gas motions and shows that they can be reconciled together. Section 3 describes the RCs derived from M31 major merger models, how they can account for asymmetric disk effects, and which requirements on the DM content are necessary to best reproduce the M31 RC. Section 4 describes at which conditions disk stars and gas can be considered at equilibrium, which is based on orbit eccentricities, on the relative amplitude of radial and vertical velocities, and on observed differences between approaching and receding sides of

<footnote n="1">¹ LIRA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, Cergy Paris Université, 92190 Meudon, France</footnote>

<footnote n="2">² Aix-Marseille Univ., CNRS, CNES, LAM, 38 rue Frédéric Joliot Curie, 13338 Marseille, France</footnote>

<footnote n="3">³ Instituto de Astrofísica, Facultad de Ciencias Exactas, Universidad Andres Bello, Fernandez Concha 700, Las Condes, Santiago RM, Chile</footnote>

<footnote n="4">⁴ Institut d'Astrophysique de Paris (UMR7095: CNRS & Sorbonne Université), 98 bis Bd Arago, 75014 Paris, France</footnote>

<footnote n="5">⁵ Dipartimento di Fisica e Astronomia “Galileo Galilei”, Universita di Padova, Vicolo Osservatorio 3, I-35122 Padova, Italy</footnote>

<footnote>* Corresponding author; francois.hammer@obspm.fr</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Comparison between Chemin et al. (2009) (green lines) and Corbelli et al. (2010) (blue dots) RCs. The former has been corrected by a factor of 0.94 to match the latter as best as possible. The dotted blue line shows the expected RC for a pure axisymmetric disk with a scale length of 6.3 kpc and a mass of $2.31 \times 10^{11} M_{\odot}$.</figcaption></figure>

the RC. Section 5 compares our results with the literature studies on M31's mass, and then with the recently determined Milky Way RC, it also shows the limits of this study.

## 2. Rotation curve of M31

The best HI RCs of M31 have been derived by Chemin et al. (2009) and Corbelli et al. (2010). The former study is based on deep 21 cm observations made at the Synthesis Telescope and the 26 m antenna at the Dominion Radio Astrophysical Observatory, and the latter is based on a deep wide-field H I Western and Green Bank Telescope survey of M31. Both surveys of the atomic HI disk of M31 have been made with a uniformly high sensitivity and resolution. Not only were the data taken at different radio telescopes, but their analysis methods also differ. The former study used terminal velocities derived from a multiple Gaussian peak analysis to get the main velocity field, while the second considered intensity-weighted mean velocities. Both studies applied a tilted-ring model to their respective HI velocity field to infer the RC and warp parameters of the disk. Figure 1 compares the two RCs that show many similarities from 8 to 35 kpc. By construction, the curve of Chemin et al. (2009) exceeds that of Corbelli et al. (2010), but both results remain very consistent within the quoted uncertainties, indicating that the terminal and average velocities trace the same gas component through most of the disk. Differences between the curves also come from small differences in the best fit parameters of the disk warping between the two studies. We find that the best match between the two curves occurs when a correction factor of 0.94 is applied to the Chemin et al. (2009) RC, confirming that terminal velocities exceed intensity-weighted mean velocities.

Chemin et al. (2009) and Corbelli et al. (2010) RCs show a very similar behavior, with a sharp increase from ~8 to 14 kpc, then a peak and a decrease beyond 16 kpc, which resembles a single bump. The major difference between the RCs occurs from ~23 to 27 kpc, while beyond ~28 kpc, both RCs show a further increase in velocity. Figure 1 reveals that the 14–15 kpc

bump can be reproduced by the effect of a single axisymmetric disk, which corresponds to a scale length equal to the bump peak divided by 2.2, that is, from 6.2 to 6.7 kpc. The latter value is in excellent agreement with M31 estimates from photometry (Hammer et al. 2007; Blaña Díaz et al. 2017, and references therein), and it allowed us to suspect that the origin of the bump indicates a quite prominent disk. We note, however, that the value adopted in Figure 1 is not realistic since having all the mass concentrated into a massive disk would generate an extremely unstable disk (Ostriker & Peebles 1973).

Because the two M31 RCs from Chemin et al. (2009) and Corbelli et al. (2010) are so similar after re-scaling, we chose to generate a M31 RC that is the average of the two after correcting the first one by a factor of 0.94. We first interpolated the Chemin et al. (2009) RC to the same intervals as that of Corbelli et al. (2010), reevaluating velocities and errors directly from the data. To estimate the final uncertainties, we quadratically added errors from Corbelli et al. (2010) and the systematic uncertainties coming from the difference between the two studies after applying the correction factor of 0.94 to the Chemin et al. (2009) RC. In the following, we compare this averaged M31 RC (see cyan line and error bars in Figures 2 and 3) to the RC derived from simulated galaxies made by the H18’s hydrodynamical simulations.

## 3. Rotation curves from merger modeling

H18 presented a set of five simulated models of M31 that provided very similar RCs (see their Appendix A), though they did not compare them to observational data. Since the stellar content of M31 is relatively well known ($M_{stars} = 1.1 \times 10^{11} M_{\odot}$), we varied the DM content to best match the observed RC. The top-left panel of Figure 2 indicates that for H18’s model 288¹, the rotational velocity is greater than the one observed (compare the full black line to the cyan points). We compared the rotational velocity of the model 288 with the observed velocity at 20 kpc. The choice of 20 kpc was guided by the fact that (1) there, the RC is not affected by the axisymmetric disk component (2) nor by non-equilibrium effects occurring at $R > 25$ kpc (see Section 4) and (3) that at 20 kpc the DM component fully dominates the baryonic component (compare the dashed black and red lines). This choice leads to a mass within 20 kpc that is close to the quadratic sum of the square velocities of the DM and baryonic components. This provides, quite accurately, a factor of 1.6 to be applied on the DM content of model 288, thus warranting that at 20 kpc the model fits the observed RC. The latter factor led us to generate a new model, namely model 371². The top-right panel of Figure 2 shows that the observed rotational velocity is quite well reproduced by the new model 371 (compare the full black line to the cyan points).

The bottom panels of Figure 2 show the position-velocity diagram provided by the HI gas particles ($T \le 20\,000$ K) for model 288 of H18 and for the new model 371, for which we decreased the DM mass by a factor of 1.6. The red and blue curves in the top panels of Figure 2 give the extremal velocities derived from the phase diagram (see bottom panels of Figure 2) after correcting them for a constant inclination ($i = 77$ degrees).

Section 2 (see also Figure 1) indicates the presence of a large bump in the RC that is likely associated with the impact of the axisymmetric disk component of M31. To account for this, the

<footnote n="1">¹ We first chose model 288 because it provides a HI disk with a size that is comparable to the one observed (see H18’s Figure 1).</footnote>

<footnote n="2">² The properties of the new model are presented at the end of this section.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Derivation of the terminal velocity for two M31 models based on the maximal velocity method of the position-velocity diagram (Chemin et al. 2009) and comparison with the observed M31 RC (see Sect. 2). Top panels: RC derived from the maximal velocity method for models 288 (left) and 371 (right). The blue, red, and black lines represent RCs of the approaching and receding sides and the average. The observed M31 RC is shown by cyan dots with error bars and a cyan line. For illustration, the contribution of stars, DM, and gas to the modeled RCs are given by red, black, and green dashed lines, assuming a spherical symmetry, and the smooth black line corresponds to their summation. Bottom panels: Position-velocity diagram plotted as a kernel density plot for gas particles with $T \le 20\,000$ K (green dots) along the receding (X < 0) and approaching (X > 0) sides, with red and blue lines representing the maximal velocities that are also shown in the top panels. The figures also account for the gas particle densities in the phase diagram, which could be reasonably well compared to Figure 5 of Chemin et al. (2009). In all panels, there are a few spikes that are due to single particles, which have been removed for further analyses (see Figure 3).</figcaption></figure>

top panels of Figure 3 show the radial evolution of the stellar surface density profile (red solid line), which can be decomposed into a bulge (red dashed line) and a disk (blue long-dashed line), together with that of the HI gas (green dashed line). Although the stellar profile of the outskirts is underestimated, particularly for model 371 beyond 22 kpc<fnref n="3" />, we may use this decomposition to predict the M31 model RCs in order to explicitly account for the asymmetric disk component. The bottom-left panel of Figure 3 confirms that model 288 overestimates the observed velocity amplitude of the HI M31 RC (reduced $\chi^2 = 9.3$), while model 371 provides a much better match between the HI observations and gas modeling (reduced $\chi^2 = 2.2$, see bottom-right panel of Figure 3). We further note that the observed increase of the gas rotational velocity beyond 28 kpc (see Section 2) is well reproduced by model 371 but not by model 288.

Another important characteristic of the HI M31 RC is the discrepancies in velocity amplitude between approaching and receding RCs at different radii. This indicates that there are some perturbations of the disk dynamics. Such discrepancies were observed by Chemin et al. (2009) between 8 and 12 kpc (inner discrepancy) and then between 17 and 21 kpc (outer discrepancy). We first noticed that while the Corbelli et al. (2010) RC

shows a similar behavior, the discrepancy between both sides has a smaller amplitude (15 versus 30 km s⁻¹) and slightly different radii: from 11 to 13 kpc and from 22 to 31 kpc for the inner and outer discrepancies, respectively. We suggest that the observational differences are related to the different methodologies used by both studies to recover the RCs (see Section 2). Blaña Díaz et al. (2017) proposed that the 10 kpc ring could be related to the bar’s outer Lindblad resonance, for which the bar pattern speed is in agreement with dynamical results from Blaña Díaz et al. (2018). This suggests that the inner discrepancy that occurs in the same radial range could also be associated with the bar. However, we still have to interpret the observed outer discrepancy between the receding and approaching velocities, which is predicted by both models 288 and 371 to occur between 17 to 30 kpc (see red and blue curves in the top panels of Figure 2). We note that the modeled discrepancy in velocity reaches approximately a maximal amplitude of 40 km s⁻¹, which is quite similar to what has been found by Chemin et al. (2009) but significantly larger than results from Corbelli et al. (2010). This could be due to the fact that we used the same methodology to extract modeled RCs as Chemin et al. (2009).

The above suggests that to best reproduce the observed M31 RC, the required missing mass has to be reduced from that considered by H18 (see also Section 5). However, before such a conclusion could be reached, we had to verify whether the main successes of H18 modeling were kept after reducing the DM content by a single factor of 1.6 (model 371). Figure 4 shows the

<footnote n="3">³ We verified that adding the corresponding missing baryonic mass does not change the RC beyond 22 kpc by more than 2 km s⁻¹ at 35 kpc and can be neglected. We also found negligible changes when passing from a thin to a thick disk for calculating the circular velocity.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Derivation of the M 31 modeled RC after decomposition of the baryonic component into a disk and bulge component for models 288 (left) and 371 (right). Top panels: Stellar mass surface density in logarithmic scale (red solid line) that is decomposed into a bulge component (red dashed line) and a disk (blue long-dashed line) together with that of the HI gas (green dashed line). The bulge and disk mass distributions result from a Plummer and exponential disk profile, respectively, and the characteristic numbers are given in the top right of each panel. Bottom panels: Same as the top panels of Figure 2 except that the black line represent the RC obtained after accounting for the bulge (red dashed line), axisymmetric disk (blue dashed line), and DM (black dashed line) components, respectively. The RC described by the maximal velocity method from Figure 2 is displayed with a magenta line, while the M 31 RC is shown by a cyan line. The blue and black dotted lines illustrate the resulting RCs when the disk is (mis)-represented by a spherical component.</figcaption></figure>

overall stellar distribution predicted by model 371, and it indeed reproduces the M 31 disk outskirts, including the giant stream, shells, and clumps, in a manner similar to model 288 (see Figure 8 of H18). We also verified that the age-velocity dispersion relation established by Dorman et al. (2015) is also reproduced. The left panel of Figure 5 shows the distribution of stellar particles in the face-on disk, and it illustrates that the new model 371 also reproduces the 10 kpc ring of star formation (see blue points that identify young stars) and the presence of a bar, which is slightly more prominent than in model 288. Furthermore, as done in H18, we verified that for a snapshot of model 371 taken 1 Gyr later (snapshot 780 instead of 680), both the outskirt morphological features and age-velocity dispersion relation are significantly diluted and merely reproduce the observations. The success of the new model 371 in reproducing M 31's overall properties is not unexpected because these exceptional properties are caused by the interaction (including resonances) of the two disks of the merger interlopers<fnref n="4" />, that is, without much impact from the DM fraction.

**4. The radius to which the mass can be estimated in M31**

Studies of RCs of remnant disks or disks that survive after a merger are not common in the literature. Here, we try to establish some simple figures to understand how such galactic disks can be stabilized over time. As for any dynamical system, the central

parts can rapidly return to equilibrium. However, they can be affected by non-axisymmetric perturbations such as those caused by a bar, and this appears to be the case for M 31 (see Section 3). The situation is different for the gas disk outskirts, which are particularly interesting here because they are the most sensitive to the presence of the DM that is expected to be dominant in the halo.

Non-axisymmetric perturbations of the outer M 31 regions can be attributed to the impact of ancient tidal tails formed during the merger, in which part of the gas and stars may have reached their apocenter before returning to the remnant galaxy with relatively large velocities. Such a mechanism explains well the formation of the giant stream (see H18), which is made by particles forming loops around the M 31 remnant after returning from a tidal tail associated with the secondary galaxy involved in the merger. The above result is a significant success of the M 31 modeling because the latter structure is found to be composed by many substructures associated with the superposition of the loops that are seen edge-on (Hammer et al. 2018; Escala et al. 2021; Tsakonas et al. 2024). Each of the five models selected by H18 to reproduce M 31 and its outskirts form two or even three tidal tails that are associated with the orbital path of the secondary galaxy. Each model also shows a tidal tail that is associated with the main galaxy and possesses a slightly tilted angular momentum compared to that of the remnant disk. The two above properties are shared by model 371 (see Figure 4 and panels d and e of Figure 6).

Figure 6 provides a description of the mechanisms at work for the merger model 371. Panel c shows how the gas orbital motions differ from circularity, indicating an offset angle (ΔAngle) that increases with radii and reaches a median value

<footnote n="4">⁴ We emphasize here that the bar occurs due to the resonance between the initial baryonic disks, while the giant stream is mostly made by particles coming from the secondary interloper.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Sky projection of the stellar density of model 371 (top panel) that can be compared to model 288 (bottom panel) as well as other similar modeling made by H18 (see their Figure 8). The arrows indicate major features, from left to right and top to bottom: NE clump, W shell, G1 clump, and the giant stream.</figcaption></figure>

close to 10 degrees at 30 kpc. Similar radial increases were found for the ratios between radial (panel a) and azimuthal (panel b) to tangential velocities, respectively. Panels d and e show a large scale map of the remnant gas that is dominated by a gigantic tidal tail (see right sides of each panel), which is directly linked to the disk. Figure 6 suggests that the radial increase of the non-circular motions are caused by particles returning from the tail. Moreover, one (the other) side of the disk has negative (positive) radial velocities, which supports well the fact that the disk is permanently shocked by gas particles infalling from the tidal tail. We note that it also explains the oscillatory behavior of the gas particles in panel c as well as the differences between the approaching and receding RCs occurring mostly beyond 17 kpc (see the differences between the green lines in Figure 1 and between the red and blue lines in the top panels of Figure 3). To verify the above, it is useful to investigate the properties of model 371 but at the much later epochs for which it is expected to reach an equilibrium even at the disk outskirts. Figure 7 is the same as Figure 6 but for the remnant taken at 7.8 instead of 2.6 Gyr after the merger. It indicates that most of the system appears relaxed. In particular, there is a flat distribution of the angle offset from circularity, the radial and azimuthal velocities, for which contribution becomes very small at all radii when compared to tangential velocities.

Figure 8 describes the orbital evolution from the merger (4.2 Gyr) to the present (6.8 Gyr) epoch. The leftmost and rightmost panels indicate that gas motions are increasingly eccentric,

from the top to the bottom panels, for which gas particles are passing from pericenter to apocenter. At 10 kpc (top panel, blue) gas motions are mostly circular (same apocenter as pericenter), while for gas particles selected at 25 and 30 kpc (pink and magenta), the apocenter and pericenter may differ by more than 10 kpc. One may assume that the main perturbation occurred at the merger epoch and that the gaseous disk system gradually regains equilibrium after a certain number of orbital times (Gnedin & Ostriker 1999). For example, Gnedin & Ostriker (1999) estimated that after a tidal shock, stars in a globular cluster retrieve the virial relation after three dynamical times, which we can translate as the number of orbital times for the disk. However, according to Gnedin & Ostriker (1999), a shocked system is still oscillating, and 10 to 15 dynamical times are necessary to reach full virial equilibrium. By similarity, this would mean that at 10 kpc, gas dynamics can be considered to fully be at virial equilibrium while still oscillating further away.

However, the gas relaxation process behaves differently than expected based on stellar system interactions. After a major merger, the gas suffers from violent relaxation and a dissipative process, which leads to a faster relaxation process than that of a stellar system. It provides a short timescale for the gas-disk formation, and even in some cases, the gas disk is implemented just before the end of the merger process (Barnes 2002). Figure 8 illustrates how the gas returning from a tidal tail gradually forms a rotating disk. Barnes (2002) has shown that tidal tail gas particles initially side swipe the central disk and then follow eccentric orbits before settling down to the disk. The accreted gas intersects and shocks, which converts kinetic energy into heat that is radiated away through cooling (Barnes & Hernquist 1996). Gas located in outer regions has more energy and less density, conditions in which the dissipation is less efficient than in the inner region. In other words, the gas reaches circular orbits preferentially in the inner region and gradually in the outer regions, as shown in Figures 6 and 8.

In the central regions, the stellar bar plays an important role in draining the angular momentum to the gaseous disk; this has been extensively studied by Hopkins et al. (2009). The stellar bar has a different rotation compared to the gaseous disk, and this small offset exerts a torque on the gas that transfers its angular momentum into the system toward the central region. Since there is a bar in the M 31, the process should also play a role in regulating the inside gaseous disk. In the outskirts, gas takes more time to adopt circular motions because it is permanently fed by the gas particles returning from the tail. Similar conclusions can be derived from the middle and right panels of Figure 8. The bottom panels indicate that gas particles have not fully reached orbital and circular motions inside the thin disk, for which $J_z/J$ and eccentricity are equal to one and zero, respectively. The middle panels of Figure 8 show that the gas has difficulty reaching precisely $J_z/J = 1$ at the outskirts. This is because the cold gas has a sound speed of ~16 km/s (for a gas temperature of 20 000 K), which is much smaller than the rotation speed, so the gas pressure force is much smaller, and the gas will roughly follow the ballistic trajectory. Because of the above arguments, Figure 8 confirms the conjecture of Gnedin & Ostriker (1999) that virial equilibrium is not reached if the gas has orbited less than three (and five) times around the M 31 center (i.e., beyond 25 kpc). This makes a calculation of mass from the virial theorem beyond that radius inapplicable.

The top-middle panel of Figure 8 also shows a very interesting property at 10 kpc because a few gas particles show almost periodic variations of their angular momentum during the whole 2.6 Gyr elapsed time since the merger. We suggest
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Stellar particle distribution of the face-on disk from model 371 within 40 kpc for snapshot 680 (left) and 1200 (right). The color of the image is based on assuming each newborn stellar particle evolves as a simple stellar population and has been observed with the filters of *u*, *g*, and *z* of the Sloan Digital Sky Survey. For example, the blue colors in the images indicate a stellar age of about 100 Myr.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Oscillation behavior of gas particles in the merger model 371 and snapshot 680 (i.e., 2.6 Gyr after the merger). Panels a and b give the ratio of the radial (and azimuthal) velocities to tangential velocities, respectively. In both panels, points and error bars give the median of the absolute values at five different radii. The colors in panel b indicate whether $\Delta\text{Angle}$ is negative (magenta) or positive (cyan). Panel c shows the radial evolution of $\Delta\text{Angle}$, which is the difference between the disk-projected velocity vector and that of a pure circular motion. Points and error bars give the median of the absolute values at five different radii. Panels d and e provide the gas particle distribution within 300 × 600 kpc² and show the presence of a tidal tail on the right sides whose angular momentum is slightly tilted versus that of the remnant disk. The colors are the same as in panels a and b, respectively.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Same as Figure 6 but for snapshot 1200 (7.8 Gyr after the merger) of model 371. The gas particles forming an almost vertical line with very negative values at R = 27 kpc are caused by a tidal dwarf that is found superimposed to the disk (see right panel of Figure 5).</figcaption></figure>

that this reveals the role of the bar, which has been proposed to be at the origin of gas shocks observed in M 31's central regions (Feng et al. 2024).

Future works investigating the properties of the bar and how to model it would be very useful, as such investigations

are beyond the scope of this paper. Preliminary results have shown that the modeled bar and its associated velocity field show strong variations with time, which could be associated with non-equilibrium conditions below 10 kpc. It is also possible that the central regions have an oscillatory behavior, suggesting a pseudo-equilibrium, that is, a system that is varying but retrieving similar conditions during periodic elapsed epochs.

**5. Discussion**

In this paper, we have shown that the M 31 RC is reproduced by model 371, which is similar to M 31 models developed by H18 but after reducing the DM mass by a factor 1.6. We have
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Temporal evolution from T = 4.2 Gyr (merger epoch) to 6.8 Gyr (present-day of the remnant M31) of randomly selected gas particle orbits at five different areas of the disk. The areas correspond to 2 kpc wide rings selected in the gas disk at T = 6.8 Gyr (snapshot 680) with radii of 10 (blue), 15 (cyan), 20 (green), 25 (pink), and 30 (magenta) kpc, respectively. Left panel: Radial orbital evolution of the gas particles from 10 kpc (top panel) to 30 kpc (bottom panel). The black dotted lines illustrate the orbital paths of characteristic gas particles in each panel. Middle panel: Ratio of the disk angular momentum ($J_z/J$) to the total angular momentum in the same rings as in the left panel. The black solid lines give the median, and the black dashed lines give the 1σ variations. Right panel: Same as the other panels but for the eccentricity evolution.</figcaption></figure>

also shown that the former H18 models systematically overestimate the M31 RC because of their too large DM fractions (see left panels of Figures 2 and 3). One may notice from Figure 3 that the baryonic mass is also larger for model 288 than for model 371. We verified that model 288’s RC still overestimates the observed RC after adopting the smaller baryonic mass associated with model 371.

Since model 371 reproduces most of M31’s properties as well as H18’s modeling, we investigated how the M31 RC can predict the dynamical mass and to which limiting radius the virial relation can be retrieved. Both observations (see Figure 1) and simulations (see left panels of Figures 2 and 3) show a significant velocity increase beyond 25 kpc. However, Figure 6 indicates that at these radii, rotational velocities are affected by instabilities caused by gas particles returning to the merger remnant from a tidal tail. These particles also affect the orbital eccentricity as well as the radial and azimuthal velocities of the gas particles in the disk outskirts. The temporal evolution of gas particle orbits since the merger epoch (4.2 Gyr) confirms the above (see left panel of Figure 8). During the 2.6 Gyr elapsed time shown in Figure 8, gas particles at disk outskirts have experienced only a small number of orbits, two to four at 30 and 25 kpc, respectively. These numbers are too small to warrant virial equilibrium with the M31 remnant gravitational potential. We consider that model 371 represents well the matter distribution in M31. Figure 9 gives the radial distribution of the DM and of the baryonic components as well as the $M_{DM}/M_{bar}$ ratio. It indicates a baryonic component that represents 32% of the M31 dynamical mass when considering it in a volume limited by $R_{200}$ = 137 kpc, a DM, and a total mass

of M31 of $M_{200}$ = $2.95 \times 10^{11} M_{\odot}$ and $M_{tot}$ = $4.5 \times 10^{11} M_{\odot}$, respectively<fnref n="5" />.

However, there is an important concern related to the baryonic mass, which may be different from one study to another. The bulge mass chosen for fitting the stellar mass distribution of model 371 is very close to that of the accurate modeling of M31’s central region (Blaña Díaz et al. 2018): 4.5 versus $4.25 \times 10^{10} M_{\odot}$. The total stellar mass of M31 in model 371 is $12 \times 10^{10} M_{\odot}$, which is an intermediate value compared to Chemin et al. (2009, $9.5 \times 10^{10} M_{\odot}$), Kafle et al. (2018, $10.3 \times 10^{10} M_{\odot}$), and Corbelli et al. (2010, $12.6 \times 10^{10} M_{\odot}$). We verified that changing the baryonic content to values adopted in other studies would not change our conclusions about the DM total mass simply because the latter dominates the RC beyond 20 kpc (see bottom-left panel of Figure 9). Given the above, we may consider that our mass determination accounts for all substructures of M31 and of its RC as well as for equilibrium conditions and that it supersedes the mass estimates made by Chemin et al. (2009) and Corbelli et al. (2010).

Bhattacharya (2023) reviewed the different mass estimates for M31. There is quite a large variety of estimates because they are based on many different probes. However, very few of them (Escala et al. 2022; Dey et al. 2023) have accounted for the fact that M31 has experienced a recent merger, which affects equilibrium conditions. For example, it might be difficult to consider orbits of probes lying at distances significantly larger than

<footnote n="5">⁵ In this paper, $R_{200c}$ is the virial radius for which the enclosed DM mass, dubbed $M_{200c}$, corresponds to an overdensity of 200 times the critical density. The term $M_{tot}$ accounts for DM and baryons, including $0.35 \times 10^{11} M_{\odot}$ of the neutral gas (see Figure 9)</footnote>
<figure><img src="image_10.png" /><figcaption>Fig. 9. Radial distribution of DM and baryonic matter. The full lines show the radial DM distribution of model 371 (snapshot 680), within 35 kpc (left) and 300 kpc (right). Long-dashed lines give the same for model 288 and snapshot 762. Top panels: Ratio of $M_{DM}/M_{bar}$ versus radius. Bottom panels: Evolution of DM (black lines), baryon (red lines), stellar (blue lines), and gas (green lines) masses.</figcaption></figure>

30–40 kpc, as they would not have had time to experience more than one or two orbits since the merger time (see bottom-left panel of Figure 8). This especially concerns mass determinations based on dwarf galaxy satellites or distant globular clusters that predict total masses from $11.4 \times 10^{11} M_{\odot}$ (Zhang et al. 2024, star clusters at distances larger than 50 and even 100 kpc), $13.5 \times 10^{11} M_{\odot}$ (Veljanoski et al. 2013, globular clusters up to 130 kpc), and $14 \times 10^{11} M_{\odot}$ (Watkins et al. 2010, dwarf galaxies within 300 kpc) to $19.9 \times 10^{11} M_{\odot}$ (Sofue 2015, all probes). In the following, we restrict comparisons to previous works that accounted for kinematic probes at distances smaller than ~35–40 kpc.

Dey et al. (2023) proposed an innovative method to estimate the mass of M31 based on the extent of the giant stream and associated shells (see also the pioneering studies of Fardal et al. 2013 and Escala et al. 2022). They derived a mass for the DM halo of $6.3 \times 10^{11} M_{\odot}$ within 125 kpc, which is twice our value ($2.95 \times 10^{11} M_{\odot}$ within 137 kpc). The study of Dey et al. (2023) is based on a progenitor of the giant stream with a much smaller mass ($10^8 M_{\odot}$) than in the H18 modeling, and it is unclear whether this can affect their determined total M31 mass. Furthermore, Dey et al. (2023) acknowledged that their determined mass may overestimate the mass by a factor of two to three (see a detailed description in Sanderson & Helmi 2013), which would render their mass estimate consistent with that from model 371.

Kafle et al. (2018) studied M31’s planetary nebulas (PNes) in an attempt to derive the escape velocity profile of M31 and then its total mass. They found an escape velocity of $470 \text{ km s}^{-1}$ after assuming a distribution of stars with absolute values of radial velocities larger than $300 \text{ km s}^{-1}$ (see their Figure 2). Using the stellar radial velocity distribution from model 371 (snapshot 680), we found a very similar figure, especially for large velocity values (see top panel of Figure 10), but an escape velocity of only $321 \text{ km s}^{-1}$. The DM mass derived by Kafle et al. (2018) is $M_{200} = 7 \times 10^{11} M_{\odot}$, which is slightly more than twice the value derived from model 371. We suspect that part of the discrepancy is related to an extrapolation of the Navarro-Frenk-White (NFW) profile (Navarro et al. 1997), which requires a very shallow profile at the outskirts (see a discussion about the impact of this extrapolation in Jiao et al. 2021). To verify this, we fit model 371’s DM-mass radial distribution with a Dehnen profile (Dehnen 1993):

$$ \rho(r) = \frac{M_0 (3-n)}{4 \pi h^3} \frac{h}{r^n (r+h)^{4-n}} \quad (1) $$

<figure><img src="image_11.png" /><figcaption>Fig. 10. Comparison with Kafle et al. (2018). Top panel: Distributions of PNe radial velocities (blue histogram) from Kafle et al. (2018) compared to predictions from model 371 (black histogram). The two red-dashed vertical lines show the limits chosen by Kafle et al. (2018) to calculate M31’s escape velocity from PNe with absolute velocity values in excess of $300 \text{ km s}^{-1}$. Bottom panel: Cumulative mass $M_{200}$ for Model 371 (red points) fit by the function of Dehnen (1993) (red line, see Eq. (1)) compared to that from Kafle et al. (2018, green-dashed line). Vertical dotted lines indicate the corresponding $R_{200}$ radii: 137 and 240 kpc, respectively.</figcaption></figure>

We found $M_0 = 5.93 \times 10^{11} M_{\odot}$, $h = 21.64$ kpc, and $n = 1.41$ (see red line in the bottom panel of Figure 10). At a large distance, this density profile follows $\rho \propto r^{-4}$ instead of $r^{-3}$ for the NFW profile, which could partly explain the mass estimate discrepancy.

The bottom panel of Figure 10 further investigates the origin of the discrepancy between the two mass estimates, and it confirms the fact that even below 15 kpc, Kafle et al. (2018) overestimated the mass (see also their Fig. 6). Roche et al. (2024) summarized the major difficulties in determining mass from escape velocities that are related to (1) the small amount of high-velocity PNe stars (a few tens of stars) used for modeling the tail of the stellar speed distribution and (2) systematics due to strong prior distributions on the fitting parameters, including strong degeneracy between parameters.

The above discussion is necessary for deriving a proper estimate of the DM total mass and its fraction to the dynamical mass. For all methods based on a limited volume such as the HI RC, the substructures surrounding the disk, or the escape velocity from PNe, the DM mass that is calculated has to be within a relatively
limited volume, and this follows the definition of what is the missing mass. The latter mass has to be evaluated in a given volume, as does the DM or missing mass fraction. Any extrapolation has to be model dependent, and it should be well described, especially to calculate the DM fraction. It means that within 35 kpc, model 371 has 1.3 and $1.7 \times 10^{11} M_{\odot}$ for DM and baryonic masses, respectively. Not surprisingly, limiting the volume in which mass is estimated automatically reduces the DM fraction.

A similar high baryonic fraction has been found by Jiao et al. (2023) when estimating the Milky Way mass from the Gaia DR3 RC. They found a total mass of $2.06 \times 10^{11} M_{\odot}$ within a radius of 122 kpc, which is made of 1.44 and $0.62 \times 10^{11} M_{\odot}$ of DM and baryons, respectively. It leads to a baryonic fraction of 30%, which does not change even when accounting for a smaller radius. Conversely, for M 31 our prediction of a baryonic fraction of 32% within a volume limited by $R_{200} = 137$ kpc<fnref n="6" />, increases to 43% if it is estimated within 35 kpc. The main reason for this difference is related to the Keplerian decline claimed to be detected by Jiao et al. (2023), which implies an accurate fit by an Einasto (Einasto 1965) profile, while the detection of a Keplerian decline in such a perturbed galaxy such as M 31 is excluded.

## 6. Conclusions

In this paper, we have generated a dynamical model of M 31 that accounts for a large set of its observed properties, from the bar to the galaxy outskirts. This modeling also succeeds in reproducing its RC detailed properties, namely, the 14 kpc bump and the velocity increase beyond 25 kpc. The former is associated with the axisymmetric disk (see Figures 1 and 3), and the latter is caused by returning material from a tidal tail (see Figure 6) that is out of equilibrium. The above is supported by the fact that for a recent (i.e., ~2.6 Gyr old) major merger, gas and stars beyond 30 kpc have no time to experience more than two orbits in the disk or in the halo.

The DM mass content of our modeling has been calibrated on the M 31 RC without affecting its accuracy in reproducing M 31's galaxy properties. Our mass estimate supersedes those assuming equilibrium for globular clusters or satellite galaxies at distances well beyond 30–40 kpc, as these probes cannot perform even a single orbit since the merger epoch. In other words, kinematical probes at large radii in such a M 31-like galaxy, well beyond 25 kpc, do not trace the underlying mass distribution because they are still not at equilibrium with M 31's remnant gravitational potential and consequently cannot be relied on to determine its total mass. Given that our modeling accurately reproduces M 31's overall properties, including its detailed RC, it also improves estimates based on former analyses of the RC. We find that the M 31 RC is consistent with a total baryonic fraction of 32%, which is far larger than that derived from former studies of external galaxies, including M 31. It is also larger than expectations from the cosmic background (Planck Collaboration VI 2020, baryonic fraction of 15.6%), which casts significant doubt on the existence of a missing baryon problem.

Finally, we have compared our results to those based on dynamics of substructures in the M 31 disk outskirts and to the escape velocity of PNes, and we find that the results are consistent after accounting for their large associated systematics. Further works for consolidating the M 31 modeling are very much wel-

comed. For example, improving the modeling of the central region kinematics could lead to interesting constraints on the DM halo content, which is somewhat linked to the bar properties. A detailed study of kinematical features discovered in the giant stream and in shells (Dey et al. 2023) is underway (see Tsakonas et al. 2024 and a paper in prep.). Other tests are likely to be done by different teams, including about the new *Hubble* Space Telescope observations of the M 31 southern disk (Chen et al. 2024).

**Acknowledgements.** We warmly thank the referee for very useful comments. We are grateful for the support of the International Research Program Tianguan, which is an agreement between the CNRS in France, NAOC, IHEP, and the Yunnan Univ. in China. Y.-J.J. acknowledges financial support from the China Scholarship Council (CSC, No.202108070090). LC is grateful to the financial support from Agencia Nacional de Investigación y Desarrollo of Chile, through the program Fondecyt Regular 1210992. HFW is supported in this work by the Department of Physics and Astronomy of Padova University through the 2022 ARPE grant: Rediscovering our Galaxy with machines.

## References

Babcock, H. W. 1939, Lick Obs. Bull., 498, 41

Barnes, J. E. 2002, MNRAS, 333, 481

Barnes, J. E., & Hernquist, L. 1996, ApJ, 471, 115

Bhattacharya, S. 2023, ArXiv e-prints [arXiv:2305.03293]

Bhattacharya, S., Arnaboldi, M., Gerhard, O., et al. 2021, A&A, 647, A130

Bhattacharya, S., Arnaboldi, M., Hammer, F., et al. 2023, MNRAS, 522, 6010

Blaňa Díaz, M., Wegg, C., Gerhard, O., et al. 2017, MNRAS, 466, 4279

Blaňa Díaz, M., Gerhard, O., Wegg, C., et al. 2018, MNRAS, 481, 3210

Bosma, A. 1978, Ph.D. Thesis, University of Groningen, The Netherlands

Chemin, L., Carignan, C., & Foster, T. 2009, ApJ, 705, 1395

Chen, Z., Williams, B., Lang, D., et al. 2024, Am. Astron. Soc. Meeting Abstr., 243, 428.06

Corbelli, E., Lorenzoni, S., Walterbos, R., Braun, R., & Thilker, D. 2010, A&A, 511, A89

Dalcanton, J. J., Williams, B. F., Lang, D., et al. 2012, ApJS, 200, 18

Dehnen, W. 1993, MNRAS, 265, 250

Delgado-Serrano, R., Hammer, F., Yang, Y. B., et al. 2010, A&A, 509, A78

Dey, A., Najita, J. R., Koposov, S. E., et al. 2023, ApJ, 944, 1

Dorman, C. E., Guhathakurta, P., Seth, A. C., et al. 2015, ApJ, 803, 24

Eggen, O. J., Lynden-Bell, D., & Sandage, A. R. 1962, ApJ, 136, 748

Einasto, J. 1965, Trudy Astrofiz. Inst. Alma-Ata, 5, 87

Escala, I., Gilbert, K. M., Wojno, J., Kirby, E. N., & Guhathakurta, P. 2021, AJ, 162, 45

Escala, I., Gilbert, K. M., Fardal, M., et al. 2022, AJ, 164, 20

Fardal, M. A., Weinberg, M. D., Babul, A., et al. 2013, MNRAS, 434, 2779

Feng, Z.-X., Li, Z., Shen, J., et al. 2024, ApJ, 963, 22

Gnedin, O. Y., & Ostriker, J. P. 1999, ApJ, 513, 626

Hammer, F., Puech, M., Chemin, L., Flores, H., & Lehnert, M. D. 2007, ApJ, 662, 322

Hammer, F., Flores, H., Puech, M., et al. 2009, A&A, 507, 1313

Hammer, F., Yang, Y., Wang, J., et al. 2018, MNRAS, 475, 2754

Hopkins, P. F., Cox, T. J., Younger, J. D., & Hernquist, L. 2009, ApJ, 691, 1168

Ibata, R., Irwin, M., Lewis, G., Ferguson, A. M. N., & Tanvir, N. 2001, Nature, 412, 49

Jiao, Y., Hammer, F., Wang, J., & Yang, Y. 2021, A&A, 654, A25

Jiao, Y., Hammer, F., Wang, H., et al. 2023, A&A, 678, A208

Kafle, P. R., Sharma, S., Lewis, G. F., Robotham, A. S. G., & Driver, S. P. 2018, MNRAS, 475, 4043

Lewis, A. R., Dolphin, A. E., Dalcanton, J. J., et al. 2015, ApJ, 805, 183

Lundmark, K. 1925, MNRAS, 85, 865

Mayall, N. U. 1951, Publ. Michigan Obs., 10, 19

Navarro, J. F., & Steinmetz, M. 2000, ApJ, 538, 477

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Ostriker, J. P., & Peebles, P. J. E. 1973, ApJ, 186, 467

Planck Collaboration VI. 2020, A&A, 641, A6

Roche, C., Necib, L., Lin, T., Ou, X., & Nguyen, T. 2024, ApJ, 972, 70

Rubin, V. C., Ford, W. K. J., & Thonnard, N. 1978, ApJ, 225, L107

Sanderson, R. E., & Helmi, A. 2013, MNRAS, 435, 378

Sofue, Y. 2015, PASJ, 67, 75

Toth, G., & Ostriker, J. P. 1992, ApJ, 389, 5

Tsakonas, C., Bhattacharya, S., Arnaboldi, M., et al. 2024, in EAS2024, 184

Veljanoski, J., Ferguson, A. M. N., Mackey, A. D., et al. 2013, ApJ, 768, L33

Watkins, L. L., Evans, N. W., & An, J. H. 2010, MNRAS, 406, 264

Williams, B. F., Dolphin, A. E., Dalcanton, J. J., et al. 2017, ApJ, 846, 145

Zhang, X., Chen, B., Chen, P., Sun, J., & Tian, Z. 2024, MNRAS, 528, 2653

<footnote n="6">⁶ However, we held to this 32% baryonic fraction because it directly comes from our modeling that predicts a DM distribution extending up to 137 kpc.</footnote>