© The Authors 2025

1. Introduction

The Milky Way (MW), Andromeda (M31), and Triangulum (M33) are the most massive members and dominant gravitational components of the Local Group (LG), whose total mass is estimated to be (3 − 5)×10¹² M_⊙ (e.g., van der Marel et al. 2012a; González et al. 2014; Carlesi et al. 2022; Benisty et al. 2022). The LG provides an ideal laboratory for investigating stellar formation and evolution, including studies of star formation histories (e.g., Tolstoy et al. 2009; Weisz et al. 2014) and constraints on evolutionary tracks of massive stars (e.g., Massey & Olsen 2003; Ekström et al. 2012). It also serves as a crucial benchmark for models of hierarchical galaxy formation and cosmology (e.g., McConnachie 2012; Klypin et al. 2002; Bullock & Boylan-Kolchin 2017). Nevertheless, a key aspect remains uncertain: the relative motions of its three most massive galaxies.

Accurate proper motion (PM) measurements of M31, M33, and other LG satellites can significantly improve estimates of the LG’s total mass and gravitational field, as well as the position and motion of its barycenter (e.g., Kahn & Woltjer 1959; Lynden-Bell 1981; van der Marel et al. 2012a; Benisty et al. 2022). This would allow for a more accurate positioning of the LG within the cosmic web and deepen our understanding of cosmological expansion, as the gravitational force competes with the cosmological expansion at LG’s boundary. On a smaller, galactic scale, precise measurements of the relative motion between M31 and the MW provide insights into their orbits and dynamical evolution, helping to distinguish between various formation histories, such as a first encounter or past accretion (Hammer et al. 2013), and predict their eventual merger (van der Marel et al. 2012b, 2019; Schiavi et al. 2020; Sawala et al. 2024). Similarly, accurate PM measurements of M31 and its potential satellite M33 are vital for understanding their orbital evolution, such as M33’s first infall scenario (Patel et al. 2017). These insights may also explain whether features such as the warp in the outer disk of M33 (Corbelli et al. 2014; Kam et al. 2017) and its tidal streams (Bekki 2008; McConnachie et al. 2009) are caused by interactions with M31.

Although of vital importance, M31’s PM has only been measured in recent decades due to its large distance. The first direct measurement was made by Sohn et al. (2012) using the Hubble Space Telescope (HST), based on multi-epoch observations of three deep fields (3.37′×3.37′) over a baseline of 5−7 years. After accounting for internal motions, van der Marel et al. (2012a) derived a systemic PM of ${({\mu }_{\alpha }^{\ast },{\mu }_{\delta })}_{\mathrm{M}31}=(43.8\pm 12.7,-31.7\pm 12.2)\mathrm{\mu }$as yr⁻¹, suggesting a nearly radial orbit between M31 and the MW. However, this estimate is strongly model-dependent (van der Marel et al. 2012a). Alternatively, M31’s PM can be inferred indirectly by modeling the line-of-sight (LOS) velocities of its satellite galaxies, under the assumption that the satellites are pressure-supported and share the bulk motion of M31 (e.g., van der Marel & Guhathakurta 2008; Salomon et al. 2016). These methods, however, tend to favor a non-negligible transverse motion for M31, in contrast to the purely radial scenario.

Most recently, high-precision astrometric data from the Gaia mission have enabled direct measurements of M31’s PM using disk stars as tracers. van der Marel et al. (2019) demonstrated the feasibility of this approach by selecting nearly 1000 disk stars from the color–magnitude diagram (CMD) constructed with Gaia Data Release 2 (Gaia DR2; Gaia Collaboration 2018) broadband photometry. They derived a PM of ${({\mu }_{\alpha }^{\ast },{\mu }_{\delta })}_{\mathrm{M}31}=(65\pm 32,-57\pm 31)\mathrm{\mu }$as yr⁻¹, suggesting a significant transverse component comparable to the radial motion, albeit with large uncertainties due to the limited precision of Gaia DR2. Taking advantage of the substantially improved astrometric accuracy in the Early Data Release 3 (Gaia EDR3; Gaia Collaboration 2021), Salomon et al. (2021) refined this measurement. By incorporating H I rotation models from Chemin et al. (2009) to correct for internal disk motions, they derived PMs using the blue disk sample and reported ${({\mu }_{\alpha }^{\ast },{\mu }_{\delta })}_{\mathrm{M}31}=(48.9\pm 10.5,-36.9\pm 8.1)\mathrm{\mu }$as yr⁻¹, in good agreement with the earlier HST result. However, measurements based on the red sample yielded significantly different values, implying a transversely dominated orbit for M31 relative to the MW. Further analysis by Rusterucci et al. (2024), who established a quasar-based reference frame and re-examined the same samples in four spatial regions, confirmed that the discrepancy between blue and red samples persists. This strongly suggests that the observed differences are not caused by spatial selection effects.

In the case of M33, the PM of its center of mass (COM) was first determined by Brunthaler et al. (2005) through Very Long Baseline Array (VLBA) observations of water masers. Thanks to its high angular resolution, the VLBA provided a robust measurement: ${({\mu }_{\alpha }^{\ast },{\mu }_{\delta })}_{\mathrm{M}33}=(23\pm 7,8\pm 9)\mathrm{\mu }$as yr⁻¹. Later, van der Marel et al. (2019) utilized data from Gaia DR2 and derived an estimate of ${({\mu }_{\alpha }^{\ast },{\mu }_{\delta })}_{\mathrm{M}33}=(31\pm 35,-29\pm 32)\mathrm{\mu }$as yr⁻¹. More recently, following a methodology similar to that used for M31, Rusterucci et al. (2024) updated the PM of M33 using Gaia DR3. Interestingly, the systematic discrepancy between blue and red stellar samples seen in M31 was not observed for M33.

Despite substantial progress over the past decade, the PMs of M31 and M33 remain uncertain. For M31, inconsistencies between different datasets, especially between the HST/Gaia blue samples and the Gaia red samples, have left the relative motion between the two most massive spirals in the LG ambiguous. To resolve this issue, further refined measurements based on purer stellar samples are necessary, as the bright sources selected from the Gaia color–magnitude diagram may be contaminated by foreground stars (Salomon et al. 2021).

In this paper, we resolve the discrepancy between Gaia blue and red samples and provide the most reliable Gaia-based PM measurements for M31 and M33. Our analysis utilizes pure supergiant samples, primarily selected from color-color or color–magnitude diagrams using UBVRI or near-infrared (NIR) photometry, and eliminates contaminants through kinematic and astrometric criteria. The sample selection process is described in detail in Sect. 2. By using supergiants with LOS velocity measurements, we establish the disk’s rotation curve and derive the systemic PMs of M31 and M33, as outlined in Sect. 3. The results, including an explanation for the origin of the discrepancy between Gaia blue and red samples, are presented in Sect. 4. Finally, we summarize and discuss our findings in Sect. 5.

2. Data

2.1. Supergiant samples

To provide accurate PM measurements, we selected reliable supergiant candidates as tracers. In recent years, the number of supergiant candidates in M31 and M33 has significantly increased. This is partly due to the effectiveness of using the Johnson Q-index to identify blue supergiants and NIR CMDs to identify red supergiants, which successfully distinguish supergiants from foreground dwarfs. More importantly, LOS velocities obtained from wide-field multi-object fiber spectrographs and high-precision astrometric measurements from Gaia provide two independent methods for robustly excluding foreground contaminants, especially for yellow supergiants that heavily overlap with foreground stars in the CMD (Drout et al. 2009, 2012; Massey & Evans 2016; Massey et al. 2016, 2021; Wu et al. 2025).

In this work, our sample of supergiants was compiled from four sources:

1. Massey et al. (2016) catalog: This compilation includes previous identifications by Massey’s group and incorporates follow-up spectroscopic observations with the Hectospec spectrograph on the MMT. We selected objects labeled as members, possible members, or unknown, excluding those identified as H II regions or clusters.

2. Red supergiants from Massey et al. (2021): Red supergiant candidates identified via NIR CMDs, with foreground contaminants excluded using Gaia DR2 astrometry.

3. Wu et al. (2025) catalog: This catalog provides systematic identifications of supergiants in M31 and M33. They selected candidates based on photometric criteria, then removed contaminants using LOS velocities from the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST; Cui et al. 2012) and astrometric data from Gaia DR3. We retained objects classified as Rank1 and Rank2.

4. Blue supergiants from LGGS photometry: We supplemented the blue end using photometric data from the Local Group Galaxies Survey (LGGS; Massey et al. 2006, 2016). We selected sources with Johnson Q-index < −0.6, a reddening-free indicator of intrinsic color. Previous studies have shown this selection yields a clean sample with minimal foreground contamination (Massey et al. 2006). To eliminate non-stellar objects, we crossmatched with catalogs including Galleti et al. (2004a, 2007), Merrett et al. (2006), Azimlu et al. (2011), Sanders et al. (2012), and Zhang et al. (2020) for M31, and Sarajedini & Mancone (2007), Ciardullo et al. (2004), Hodge et al. (1999), and Zhang et al. (2020) for M33. Finally, we retained sources with low reddening, defined as G_BP − G_RP < 0.5.

After crossmatching with Gaia DR3 within a radius of 2″, we obtained a total of 10 173 supergiant candidates in M31 and 10 883 in M33. Among them, 8450 candidates in M31 and 8691 in M33 have available PM measurements (i.e., a five-parameter or six-parameter astrometric solution).

2.2. Sample selection

To improve the purity of the disk tracer sample, we applied a cleaning procedure to the supergiant candidates. This involved a series of selection criteria designed to remove potential contaminants and sources with unreliable astrometric measurements. The specific steps of this procedure, along with the number of remaining sources after each step, are summarized in Tables 1 and 2.

For M31, we adopted the following basic parameters: COM position (RA, Dec) = (10.68°, 41.27°) (Evans et al. 2010), inclination i = 77.5° (Walterbos & Kennicutt 1988; de Vaucouleurs et al. 1991), and position angle PA = 37.5° (Chemin et al. 2009; Corbelli et al. 2010). We then applied a series of selection criteria to exclude contaminants and sources with poor measurements, as is detailed below:

Step 1. Taking the angular size of M31’s disk to be 1.8°, its projected shape on the sky is approximated by an ellipse. We excluded all sources lying outside this elliptical boundary, resulting in a sample of 8397 sources for further cleaning.

Step 2. We corrected the parallax zero-point following Lindegren et al. (2021a), and removed sources whose distances exceed that of M31 (784 kpc; Stanek & Garnavich 1998) at the 2σ confidence level. Specifically, we retained sources satisfying ϖ − 2σ_ϖ < 1/784, leaving 7615 sources.

Step 3. To exclude potential non-single stars, we applied three astrometric quality criteria from Fabricius et al. (2021):

1. RUWE < 1.4;

2. ipd_frac_multi_peak < 2;

3. ipd_gof_harmonic_amplitude < 0.1.

This step reduced the sample to 6614 sources.

Step 4. To mitigate the effects of flux excess in crowded or faint regions, we adopted the criterion from Riello et al. (2021), and retained sources whose corrected excess factor (C_corr) satisfies

$$\begin{array}{c}\hfill |C_{\mathrm{corr}}|<3\times (5.9898\times {10}^{-3}\times 8.817481\times {10}^{-12}\times G^{7.618399}).\end{array}$$(1)

After this cut, 4829 sources remained.

Step 5. We excluded sources with G-band magnitudes fainter than 20 mag to ensure the reliability of astrometric measurements, yielding 2880 sources.

Step 6. Finally, we applied PM cuts: $|{\mu }_{\alpha }^{\ast }|<0.19\mathrm{mas}{\mathrm{yr}}^{-1}+2{\sigma }_{{\mu }_{\alpha }^{\ast }}$ and |μ_δ|< 0.19 mas yr⁻¹ + 2σ_μδ, where ${\sigma }_{{\mu }_{\alpha }^{\ast }}$ and σ_μδ are the PM uncertainties in right ascension and declination, respectively. The threshold of 0.19 mas yr⁻¹ follows the criterion adopted in Salomon et al. (2021) and Rusterucci et al. (2024). These limits correspond to rejecting sources with transverse velocities exceeding around 750 km s⁻¹ at the distance of M31 in either direction. This final selection yielded 2462 sources, which are considered a clean sample of disk members for subsequentanalysis.

For M33, we applied a screening process similar to that used for M31, with appropriate adjustments. The adopted basic parameters were: COM position (RA, Dec) = (23.46°, 30.66°) (van den Bergh 2000), inclination i = 54° (de Vaucouleurs et al. 1991), and position angle PA = 23° (de Vaucouleurs et al. 1991; McConnachie et al. 2006). Contaminants and sources with poor-quality measurements were excluded through the following steps:

Step 1. Assuming an angular size of 1.0° for M33’s disk, its projected shape on the sky is approximated by an ellipse, within which all 8691 sources were located.

Step 2. Adopting a distance to M33 of 840 kpc (Galleti et al. 2004b; Breuval et al. 2023), we retained sources that satisfy ϖ − 2σ_ϖ < 1/840, yielding 8076 sources.

Steps 3–5. The same astrometric quality and photometric excess criteria as applied for M31 were used here. The number of sources remaining after each step is listed in Table 2.

Step 6. The PM limits were adjusted to $|{\mu }_{\alpha }^{\ast }|<0.14\mathrm{mas}{\mathrm{yr}}^{-1}+2{\sigma }_{{\mu }_{\alpha }^{\ast }}$ and |μ_δ|< 0.14 mas yr⁻¹ + 2σ_μδ, corresponding to a transverse velocity threshold of about 550 km s⁻¹ at the distance of M33. This threshold is similar to that adopted in Salomon et al. (2021) and Rusterucci et al. (2024), although it is slightly more stringent. We consider it more robust given M33’s lower disk rotation velocity (see Sect. 3.1) and shallower gravitational potential compared to M31. The final cut left 2282 sources, which are considered a clean sample of disk members for subsequent analysis.

3. Analysis

3.1. Rotation curve from disk objects

To accurately assess the contribution of disk rotation to the observed PM of supergiant stars, we constructed rotation curve (RC) models for the disks of M31 and M33. The modeling samples were selected from those identified through the six steps described above (see Tables 1 and 2), with the additional requirement that only stars with available heliocentric LOS velocity measurements were included. This selection ensures that the samples used to construct the RCs are both clean and robust.

3.1.1. M31 rotation curve

For M31, among the 2462 sources listed in Table 1 after applying the six selection steps, 321 have LOS velocity measurements from Drout et al. (2009), Massey et al. (2009), Cordiner et al. (2011), Massey & Evans (2016), and Wu et al. (2025). As young tracers, these supergiant candidates experience negligible asymmetric drift, which can be ignored in the RC calculations. Since the LOS velocities from different catalogs were obtained using different instruments, we applied the zero-point offsets derived by Wu et al. (2025) to adjust the MMT LOS velocities to the LAMOST scale. The spatial distribution of these 321 sources is shown in panel a of Fig. 1.

Fig. 1. Panel (a): Spatial distribution of 321 supergiant candidates for modeling the RC of M31 disk region. The red star denotes the COM of M31, and the dashed black outer line represents the 1.8° ellipse of the projected disk. Inner and outer green ellipses denote the reliable range of our derived RC. Panel (b): RC of the M31 disk region. The orange line and the shaded area represent the best fit and the corresponding standard error from this work, while the blue and green error-bars represent results from previous studies, as labeled in the bottomright corner. Panel (c): Spatial distribution of 235 supergiant candidates for modeling the RC of M33 disk region. The red star denotes the COM of M33, and the dashed black outer line represents the 1° ellipse of the projected disk. Inner and outer green ellipses denote the reliable range of our derived RC. Panel (d): RC of the disk of M33. The orange line and the shaded area illustrate the best fit and its corresponding standard error from this work. The blue and green error-bars represent results obtained from previous studies by H I observations, as labeled in the bottomright corner.

Equipped with observed LOS velocities, we followed a procedure similar to Rubin & Ford (1970) to derive the RC. The LOS velocity of a disk target can be expressed as

$$\begin{array}{c}\hfill V_{\mathrm{LOS}}=V_{\mathrm{LOS}}^{\mathrm{M}31}+V_c(R)sinicos\theta .\end{array}$$(2)

Here, cos θ = X/R, where X represents the distance along the semimajor axis and R is the radial distance of the object within the M31 disk plane, and i is the inclination angle setting to be 77.5° (Walterbos & Kennicutt 1988; de Vaucouleurs et al. 1991). $V_{\mathrm{LOS}}^{\mathrm{M}31}$ is the systemic LOS velocity of M31 and V_c is the circular speed at R.

Instead of assuming a constant V_c, we modeled the RC using a fifth-order polynomial to account for potential variations, following the approach adopted in previous studies (e.g., Rubin & Ford 1970; Rastorguev et al. 2017). A Markov chain Monte Carlo (MCMC) method was employed to simultaneously fit both $V_{\mathrm{LOS}}^{\mathrm{M}31}$ and the coefficients of the polynomial. To facilitate the fitting, the sources were divided into eight radial bins, each containing at least 10 stars. We derived a $V_{\mathrm{LOS}}^{\mathrm{M}31}$ of −300.8 ± 1.2 km s⁻¹, consistent with previous measurements (van der Marel & Guhathakurta 2008; Watkins et al. 2013; Salomon et al. 2016).

For the RC construction, we consider results outside the radial range of 3.35 to 22.1 kpc from the M31 center to be unreliable due to the limited number of sources. These boundaries are indicated as green ellipses in panel a of Fig. 1. The RC within this range, together with previous measurements from Chemin et al. (2009) and Zhang et al. (2024), is shown in panel b of Fig. 1. Overall, our RC is consistent with previous results within the 3σ confidence level. Specifically, it closely follows the H I observations of Chemin et al. (2009) between around 4 and 15 kpc and then trends toward an intermediate position between their results and the latest measurements from Zhang et al. (2024), which are based on multiple tracers. The final expression for our RC is

$$\begin{array}{cc}\hfill V_c(R)=& -0.0004109\times R^5+0.02668\times R^4\hfill \\ \hfill & -0.5688\times R^3+3.470\times R^2+17.23\times R+92.44.\hfill \end{array}$$

With the RC in hand, we can precisely determine the contribution of disk rotation to the motion of each disk object. Accordingly, we revised the selection criteria from Step 6 in Sect. 2.2 as

$$\begin{array}{cc}& |{\mu }_{\alpha }^{\ast }-{\mu }_{\alpha ,\mathrm{R}C}^{\ast }|<0.14\mathrm{mas}{\mathrm{yr}}^{-1}+{\sigma }_{{\mu }_{\alpha }^{\ast }},\hfill \\ \hfill & |{\mu }_{\delta }-{\mu }_{\delta ,\mathrm{R}C}|<0.14\mathrm{mas}{\mathrm{yr}}^{-1}+{\sigma }_{{\mu }_{\delta }},\hfill \end{array}$$

where ${\mu }_{\alpha ,\mathrm{RC}}^{\ast }$ and ${\mu }_{\delta ,\mathrm{RC}}$ are the predicted PM components arising solely from disk rotation, as derived from our RC models. Given the significant internal kinematic variations across M31 and their contributions to the observed PMs, these revised criteria are more physically motivated than those in Step 6 originally proposed by van der Marel et al. (2019) and Salomon et al. (2021), and effectively remove nonphysical outliers and hyper-velocity stars with deviations exceeding about 500 km s⁻¹ from the expected RC. We note that a 1σ threshold was adopted for the cuts, instead of the 2σ used in the original Step 6, to more effectively exclude contaminants or sources with largeuncertainties.

Using above cuts and selecting only supergiant stars within the radius range of 3.35 to 22.1 kpc (the valid range of the RC model), a total of 1587 sources remain. Among them, 783 are blue sample stars (blue and yellow supergiants) and 804 are red sample stars (red supergiants). The spatial distribution and CMD positions of these sources are shown in the left-hand panels of Fig. 2.

Fig. 2. Spatial and CMD distributions of final M31 and M33 supergiant samples (see Sect. 3.1). The left-hand panels display the spatial and CMD distributions of the M31 sample stars, while the right-hand panels present those of the M33 sample stars. Filled blue and red symbols denote blue (including some yellow) and red supergiants with a 5p astrometric solution, respectively; open cyan and pink circles indicate blue and red supergiants with a 6p solution. The dashed magenta line in the CMD panels marks GBP − GRP = 2.6. Clearly, M31 contains a significantly larger population of extremely red supergiant stars compared to M33.

3.1.2. M33 rotation curve

For M33, 235 out of the 2282 sources remaining after the six selection steps (Table 2) have heliocentric LOS velocity measurements (Drout et al. 2012; Wu et al. 2025). Again, we applied a zero-point correction to the MMT LOS velocities to align their scale with that of LAMOST. The spatial distribution of these sources is shown in panel c of Fig. 1.

Following the procedure similar to that used for M31, we derived the RC of M33’s disk region using the LOS velocities of disk objects, which can be predicted from RC by

$$\begin{array}{c}\hfill V_{\mathrm{LOS}}=V_{\mathrm{LOS}}^{\mathrm{M}33}+V_c(R)sinicos\theta .\end{array}$$(3)

Here the inclination angle i was adopted as 54.0° from de Vaucouleurs et al. (1991). The definitions of $V_{\mathrm{LOS}}^{\mathrm{M}33}$ and θ are similar to those in Eq. (2) for M31. Again, the RC was modeled using a fifth-order polynomial. The best fitting results were then derived using an MCMC method. The derived $V_{\mathrm{LOS}}^{\mathrm{M}33}$ is −179.5 ± 0.6 km s⁻¹, in line with previous studies (van der Marel & Guhathakurta 2008; Kam et al. 2017). The RC is considered reliable within 0.38−8 kpc from M33’s center, as indicated by the green ellipses in panel c of Fig. 1. Within this range, our results agree with previous H I-based RC measurements by Corbelli et al. (2014) and Kam et al. (2017), particularly showing strong consistency between 1 and 8 kpc (panel d of Fig. 1). The best-fit RC is expressed as

$$\begin{array}{cc}\hfill V_c(R)=& -0.01626\times R^5+0.2555\times R^4\hfill \\ \hfill & -0.6677\times R^3-7.518\times R^2+50.59\times R+11.7.\hfill \end{array}$$

For the same reason as for M31, we revised the criteria in Step 6 of Sect. 2.2 for M33 as

$$\begin{array}{cc}& |{\mu }_{\alpha }^{\ast }-{\mu }_{\alpha ,\mathrm{R}C}^{\ast }|<0.11\mathrm{mas}{\mathrm{yr}}^{-1}+{\sigma }_{{\mu }_{\alpha }^{\ast }},\hfill \\ \hfill & |{\mu }_{\delta }-{\mu }_{\delta ,\mathrm{R}C}|<0.11\mathrm{mas}{\mathrm{yr}}^{-1}+{\sigma }_{{\mu }_{\delta }}.\hfill \end{array}$$

After applying these revised cuts, a total of 1429 sources remain within the radial range of 0.38 to 8 kpc (the valid range for the RC model). This includes 772 blue sample stars (blue and yellow supergiants) and 657 red sample stars (red supergiants). Their spatial distribution and CMD positions are shown as blue and red dots in the right-hand panels of Fig. 2.

3.2. Method for proper motion determination

The method adopted here to determine the systemic PMs of M31 and M33 is generally similar to those used in previous studies, such as Salomon et al. (2021) and Rusterucci et al. (2024). The dispersion and asymmetric drift are neglected in our modeling given the young, kinematically cold nature of the massive supergiant stars used in this study. In other words, the motions of these supergiants are assumed to follow perfectly circular orbits within the disks of M31 and M33. Under this assumption, the expected PM of a supergiant consists of two components: one arising from its circular motion within the disk, and the other representing the systemic PM of the galaxy. This can be written as

$$\begin{array}{c}\hfill {\mathit{\mu }}^{\mathrm{model}}={\mathit{\mu }}_{\mathrm{disk}}+{\mathit{\mu }}_{\mathrm{sys}},\end{array}$$(4)

where μ_disk is the contribution from the circular motion, which can be precisely predicted at any position within M31 and M33 based on the RCs derived in Sects. 3.1.1 and 3.1.2. μ_sys denotes the systemic PMs of M31 and M33, respectively, which are the quantities we aim to determine. To infer the systemic PMs of M31 and M33, we used an MCMC method to sample the posterior distribution. The probability density function for one tracer at position (α_k, δ_k) can be written as

$$\begin{array}{cc}\hfill f({\mu }_{\alpha ,k}^{\ast },{\mu }_{\delta ,k})=& \frac{1}{2\pi {\sigma }_{{\mu }_{\alpha ,k}^{\ast }}{\sigma }_{{\mu }_{\delta ,k}}\sqrt{1-{\rho }_k^2}}\hfill \\ \hfill & \times exp(-\frac{1}{2(1-{\rho }_k^2)}[\frac{{\mathrm{\Delta }}_{{\mu }_{\alpha ,k}^{\ast }}^2}{{\sigma }_{{\mu }_{\alpha ,k}^{\ast }}^2}+\frac{{\mathrm{\Delta }}_{{\mu }_{\delta ,k}}^2}{{\sigma }_{{\mu }_{\delta ,k}}^2}-\frac{2{\rho }_k{\mathrm{\Delta }}_{{\mu }_{\alpha ,k}^{\ast }}{\mathrm{\Delta }}_{{\mu }_{\delta ,k}}}{{\sigma }_{{\mu }_{\alpha ,k}^{\ast }}{\sigma }_{{\mu }_{\delta ,k}}}]),\hfill \end{array}$$(5)

where ${\mathrm{\Delta }}_{{\mu }_{\alpha ,k}^{\ast }}$ and Δ_μδ, k are the differences between the observed PM of the k-th supergiant star and the modeled value given by Eq. (4). Here, ${\sigma }_{{\mu }_{\alpha ,k}^{\ast }}$ and σ_μδ, k are the PM uncertainties in right ascension and declination, respectively, taken from Gaia DR3. In addition, the measured correlation between the two components, ρ_k ≡ pmra_pmdec_corr, is also taken into account. To evaluate the impact of potential residual contamination on the derived PMs, we first adopted a likelihood function that includes a contamination component, following the functional form used in Salomon et al. (2021). In this model, the observed PM distribution is treated as a mixture of two populations, true M31 members and contaminants, with the likelihood given by

$$\begin{array}{c}\hfill ln\mathcal{L}={\displaystyle \sum _{k=1}^n}ln[(1-{\eta }_c)f({\mu }_{\alpha ,k}^{\ast ,t},{\mu }_{\delta ,k}^t)+{\eta }_{c}g({\mu }_{\alpha ,k}^{\ast ,c},{\mu }_{\delta ,k}^c)],\end{array}$$(6)

where η_c is the contamination fraction, $({\mu }_{\alpha ,k}^{\ast ,t},{\mu }_{\delta ,k}^t)$ represent the PMs of true members, and $({\mu }_{\alpha ,k}^{\ast ,c},{\mu }_{\delta ,k}^c)$ represent those of contaminants. The contaminant probability density function, g, has the same form as Eq. (5), but with distinct mean and dispersion parameters from those of the true members. Following Salomon et al. (2021), we imposed uniform priors on the contaminant dispersions: $0\le {\sigma }_{{\mu }_{\alpha }}^{\ast ,c},{\sigma }_{{\mu }_{\delta }}^c\le 2$ mas yr⁻¹. We performed MCMC sampling using this likelihood and found that the inferred contamination fraction is below 0.1%, indicating that our supergiant sample is highly pure and any residual contaminants have negligible impact on the PM results. Therefore, for the subsequent analysis, we adopted a simplified likelihood function without the contamination component to infer the systemic PM via the MCMC technique:

$$\begin{array}{c}\hfill ln\mathcal{L}={\displaystyle \sum _{k=1}^n}lnf({\mu }_{\alpha ,k}^{\ast },{\mu }_{\delta ,k}).\end{array}$$(7)

4. Result and discussion

4.1. Origin of the discrepancy between M31 blue and red samples

Using the carefully selected supergiant star samples and the established methodology described above, we remeasured the PM of M31. As is shown clearly in Fig. 2, the supergiant star sample comprises two distinct populations: blue and red supergiant stars¹. These populations are broadly consistent with those used in previous studies (Salomon et al. 2021; Rusterucci et al. 2024). We first derived the PM of M31 using the two subsamples separately. From the blue supergiants, we obtained ${({\mu }_{\alpha }^{\ast },{\mu }_{\delta })}_{\mathrm{M}31,\mathrm{Blue}}=(37.7\pm 5.2,-17.0\pm 4.4)\mathrm{\mu }$as yr⁻¹, while the red supergiants yielded ${({\mu }_{\alpha }^{\ast },{\mu }_{\delta })}_{\mathrm{M}31,\mathrm{Red}}=(15.3\pm 7.7,-50.3\pm 6.3)\mathrm{\mu }$as yr⁻¹. The PM derived from the red sample is significantly large, indicating a non-negligible transverse motion of M31, whereas the result from the blue sample suggests a nearly negligible tangential component. This apparent discrepancy between the two results echoes similar tensions reported by Salomon et al. (2021) and Rusterucci et al. (2024).

As for the cause of this discrepancy, Rusterucci et al. (2024) investigated the influence of spatially dependent zero-point offsets using quasars (QSOs) and found that the mean motion remains close to 0 μas yr⁻¹. Further analysis, dividing the sky into four quadrants, shows that the discrepancy between the blue and red supergiant samples persists across all regions, suggesting that spatial systematics are unlikely to be the dominant cause.

Instead, we recall two key aspects highlighted by the Gaia team during astrometric validation. The first concerns the nature of the astrometric solutions in Gaia (E)DR3, which include two-parameter (position only), five-parameter (5p, position, parallax, and PM), and six-parameter (6p, 5p plus pseudo-color) solutions (Lindegren et al. 2021b). The 6p solution is typically adopted when the effective wavenumber (ν_eff), which provides essential color information for correcting chromatic effects in the point-spread function, cannot be reliably determined, particularly in crowded fields or areas affected by bright stars. In such cases, a pseudo-color (${\widehat{\nu }}_{\mathrm{eff}}$) is introduced and solved simultaneously with the five standard astrometric parameters. As a result, the 6p solution generally yields less precise astrometry and exhibit both larger random and systematic uncertainties than the 5p solution, due to more complex observational conditions and additional correlations introduced by the pseudo-color term (Lindegren et al. 2021b; Fabricius et al. 2021). The second factor is the color dependence of astrometric quality. Based on open cluster members as calibrators, Fabricius et al. (2021) demonstrated that redder faint sources (G > 18) show significant systematic offsets in PM relative to the cluster median, particularly toward the red end of the color distribution (G_BP − G_RP > 2.0; see their Figs. 17, 18, 24, and 25).

Given the two potential issues discussed above, we began by examining the composition of our samples in terms of 5p and 6p astrometric solutions. In M31, the blue supergiant sample contains 232 stars with a 5p solution and 551 with a 6p one, whereas the red supergiant sample is heavily dominated by 6p sources, with 775 stars having a 6p solution and only 29 having a 5p one. A similar pattern is observed in M33, where the red sample again shows a significantly higher proportion of 6p sources compared to the blue. The spatial and CMD distributions of these stars, categorized by both color and astrometric solution type, are presented in Fig. 2.

To assess whether the observed discrepancies are linked to systematics associated with astrometric solution type or stellar color, we examined the PM distribution as a function of G_BP − G_RP for four subsamples in M31: blue 5p, blue 6p, red 5p, and red 6p, as is shown in the left-hand panels of Fig. 3. In this analysis, the PM values have been corrected for disk rotation. For the blue and red 6p subsamples, each color bin contains at least 130 and 150 sources, respectively, spanning the G_BP − G_RP range in steps of more than 0.2 mag. The blue 5p subsample is similarly binned, with each comprising at least 115 sources. The filled red circle represents the 29 red stars with a 5p solution, while the last open pink circle denotes 155 red 6p sources with G_BP − G_RP > 2.6. Two notable trends can be found:

Fig. 3. Left-hand panels: PM distribution of the blue and red M31 samples as a function of GBP − GRP. The data points in the left (right) panel show the weighted average differences between the observed PMs, ${\mu }_{\alpha ,\mathrm{obs}}^{\ast }$ (μδ, obs), from Gaia DR3 and the predicted rotational components ${\mu }_{\alpha ,\mathrm{RC}}^{\ast }$ (${\mu }_{\delta ,\mathrm{RC}}$), derived from the RC in Sect. 3.1.1. Weights are proportional to 1/σα*2 or 1/σδ2, and the errors, also weighted by uncertainty, are indicated accordingly. Right-hand panels: Similar to the left-hand panels, but for M33.

1. Among the 6p sources, PM values exhibit a marked discontinuity at the red extreme (G_BP − G_RP > 2.6), with significant deviations from the rest of the sample, particularly in the δ direction. The spatial distribution of these sources is relatively uniform, suggesting that the offset is unlikely to originate from spatial effects. Instead, this likely reflects systematic errors in the Gaia 6p astrometric solution. This finding is consistent with Fabricius et al. (2021), who reported substantial PM offsets in faint, red sources, populations typically dominated by the 6p solution given their faint nature (most with G > 18; Lindegren et al. 2021b). Due to the limited number of red stars with a 5p solution, it remains unclear whether they are similarly affected. However, the 29 red 5p stars, despite having a mean color near G_BP − G_RP = 2.6, exhibit PM values comparable to those of bluer stars with either a 5p or 6p solution. This suggests that the 5p solution, regardless of color, is generally reliable in our analysis.

2. A systematic offset is also observed between 5p and 6p solutions: stars with a 5p solution exhibit systematically higher PM values in the α direction compared to their 6p counterparts, independent of color.

We then examined the PM distribution for M33 using the same approach. As is shown in the right-hand panels of Fig. 3, the blue sample includes at least 90 (125) stars per bin for the 5p (6p) subsamples, with G_BP − G_RP binned over intervals of 0.2 (0.125) mag. For the red 6p sample with G_BP − G_RP ≤ 2.6, each bin contains at least 125 stars in 0.2 mag intervals. Additionally, the final open pink circle represents the 23 red 6p stars with G_BP − G_RP > 2.6. The overall patterns observed in M33 are consistent with those found in M31.

From these tests, we conclude that the discrepancy between the blue and red samples primarily arises from systematic issues associated with the 6p astrometric solution. These issues are twofold: (1) unreliable astrometry for sources at the extreme red end, particularly those with G_BP − G_RP > 2.6, and (2) systematic offsets between 5p and 6p solutions, regardless of color. These effects are further amplified by the differing proportions of 5p and 6p solutions in the two samples: the red sample is dominated by 6p sources, while the blue sample includes a substantial fraction of 5p sources (approximately 30% in M31 and 35% in M33). As a result, the red sample is more vulnerable to the systematics of the 6p solution.

To address these two issues, we first mitigated the impact of unreliable astrometry by restricting the 6p sample to sources with G_BP − G_RP ≤ 2.6. Second, to account for systematic offsets, we examined potential differences in the zero-point corrections between the 5p and 6p solutions using QSOs, as described below.

4.2. Proper motion zero-point offset

Previous studies have shown that zero-point offsets still exist in Gaia PM measurements (Fabricius et al. 2021; Salomon et al. 2021; Rusterucci et al. 2024). Building on the discussion above in Sect. 4.1, we used QSOs with both 5p and 6p astrometric solutions, whose PMs are expected to be zero due to their extragalactic nature, to determine the corresponding zero-point offsets in the directions of M31 and M33 for each solution type, and to assess whether systematic differences exist between the 5p and 6p solutions. The QSO sample was constructed by crossmatching sources from the Gaia DR3 in_qso_candidates table with the QSO catalog compiled by Liao et al. (2019). This crossmatch serves to ensure the purity of the selected QSO sample. To further ensure the statistical reliability of the 6p solution analysis, we supplemented the QSO sample with additional Gaia QSO candidates flagged with astrometric_selection_flag == 1 for the 6p solution, which have an estimated purity of 98% (Gaia Collaboration 2023).

For M31, we selected background QSOs within an angular radius of 10° from the galaxy center. This selection was motivated by the 5 − 10 μas yr⁻¹ zero-point variations reported over angular scales of 4° −20° by Salomon et al. (2021). Limiting the sample to within 10° reduces spatial systematics while preserving a sufficient number of QSOs for a robust zero-point estimate. This yields approximately 4200 5p and 2900 6p sources. To further ensure the sample’s quality and to robustly estimate the zero-point offsets, we apply the following five filtering criteria–similar to those used for the supergiant sample:

1. ω − 2σ_ω < 0;

2. ${\mu }_{\alpha }^{\ast }-{\sigma }_{{\mu }_{\alpha }^{\ast }}<0<{\mu }_{\alpha }^{\ast }+{\sigma }_{{\mu }_{\alpha }^{\ast }}$;

3. μ_δ − σ_μδ < 0 < μ_δ + σ_μδ;

4. RUWE < 1.4, ipd_frac_multi_peak < 2, ipd_gof_harmonic_amplitude < 0.1;

5. astrometric_matched_observations ≥ 10.

After applying these criteria, 1670 QSOs with a 5p solution and 931 with a 6p solution remain. Considering the spatial dependence of zero-point offsets noted earlier (e.g., Salomon et al. 2021; Rusterucci et al. 2024), typically on the order of 5 − 10 μas yr⁻¹, these variations cannot explain the much larger PM discrepancy of about 50 μas yr⁻¹ observed between the blue and red samples (Rusterucci et al. 2024). We therefore examined the color dependence of the zero-point offsets as a function of G_BP − G_RP. To ensure statistical robustness, the full sample of 2601 sources was divided into 13 color bins, each containing at least 200 sources and spanning a minimum color width of 0.04 mag. We inferred the zero-point offsets in each bin using the method described in Sect. 3.2, assuming the QSOs have no intrinsic motion, i.e., ${\mathit{\mu }}^{\mathrm{model}}={\mathit{\mu }}_{\mathrm{sys}}$. Under this assumption, we applied Eq. (7) to derive the posterior distributions in each color bin. The results, shown as gray symbols in the left-hand panels of Fig. 4, indicate no significant color-dependent trend within the relatively narrow color range of the QSO sample, which spans from the red end of the blue supergiant sample to the blue end of the red supergiant sample. For the entire unbinned QSO sample, the overall zero-point offsets are ${({\mu }_{\alpha }^{\ast },{\mu }_{\delta })}_{\mathrm{zpt}}^{\mathrm{M}31}=(-2.6\pm 5.1,-2.7\pm 4.1)\mathrm{\mu }$as yr⁻¹, estimated using the same posterior inference method, and are largely consistent with those reported previously (Salomon et al. 2021; Rusterucci et al. 2024).

Fig. 4. Left-hand panels: Inferred PM values of M31 background QSOs in right ascension (${\mu }_{\alpha }^{\ast }$; left) and declination (μδ; right) as a function of GBP − GRP color. The upper sub-panels show the inferred PMs in each color bin for the full QSO sample (gray points), with the dashed gray line and shaded region indicating the overall zero-point offset and its 1σ uncertainty from the MCMC analysis. The lower sub-panels show results for QSOs with different astrometric solutions: blue and orange points correspond to the inferred PMs in each color bin for 5p and 6p sources, respectively. Dashed lines and shaded areas denote the overall PMs and their 1σ uncertainties for the two subsamples. Right-hand panels: Same as the left-hand panels, but for M33 background QSOs.

Then, using the same sample, we separately examined the potential color dependence of the zero-point offsets for the 5p and 6p QSO subsamples, and assessed whether any systematic differences exist between the two solution types. Following the same approach as above, the sources were binned by G_BP − G_RP color, with each bin containing at least 150 sources for the 5p sample and 100 for the 6p sample, and a minimum color width of 0.04 mag. The zero-point offsets in each bin were inferred using the same method described previously and are shown as blue (5p) and orange (6p) points in the left-hand panels of Fig. 4. No clear color-dependent trend is observed within either solution type; however, the 5p and 6p sources exhibit notably different zero-point offsets. The offsets for the 6p solution are substantially larger in magnitude and have larger uncertainties compared to those of the 5p solution. The overall zero-point offsets for each solution type were then derived using the posterior inference method, yielding ${({\mu }_{\alpha }^{\ast },{\mu }_{\delta })}_{\mathrm{zpt},5\mathrm{p}}^{\mathrm{M}31}=(-1.5\pm 5.4,-2.1\pm 4.3)\mathrm{\mu }$as yr⁻¹ and ${({\mu }_{\alpha }^{\ast },{\mu }_{\delta })}_{\mathrm{zpt},6\mathrm{p}}^{\mathrm{M}31}=(-11.5\pm 16.8,-6.8\pm 13.0)\mathrm{\mu }$as yr⁻¹. As was expected, the uncertainties in the 6p solution are approximately three to four times larger than those of the 5p solution.

For M33, we applied a similar procedure to determine the zero-point offsets, using approximately 7000 5p and 3800 6p background QSO candidates located within 10° from M33 center. Following the same cleaning process as for M31, 1961 QSOs with a 5p solution and 840 with a 6p solution were finally selected. Similar to M31, we first examined whether the zero-point offsets for M33 are dependent on G_BP − G_RP color. Accordingly, the full sample of 2801 sources was divided into 14 bins, each containing more than 200 sources ranging in magnitude by over 0.04 mag. The inferred PM values in each bin are shown as the gray symbols in the right-hand panels of Fig. 4. As with M31, no significant color-dependent trend is observed across the QSO sample’s color range for M33. The overall zero-point offsets were derived as ${({\mu }_{\alpha }^{\ast },{\mu }_{\delta })}_{\mathrm{zpt}}^{\mathrm{M}33}=(2.8\pm 6.4,2.0\pm 5.3)\mathrm{\mu }$as yr⁻¹. We then derived the zero-point offsets separately for the 5p and 6p subsamples, with the results again presented in the right-hand panels of Fig. 4. The pattern closely resembles that of M31: both the 5p and 6p sources exhibit nonzero, color-independent PM zero-points, but the offsets differ significantly, with the 6p solution generally showing larger offsets and greater uncertainties. For M33, the inferred overall zero-point offsets for the 5p and 6p solutions are: ${({\mu }_{\alpha }^{\ast },{\mu }_{\delta })}_{\mathrm{zpt},5\mathrm{p}}^{\mathrm{M}33}=(4.4\pm 6.7,2.6\pm 5.5)\mathrm{\mu }$as yr⁻¹ and ${({\mu }_{\alpha }^{\ast },{\mu }_{\delta })}_{\mathrm{zpt},6\mathrm{p}}^{\mathrm{M}33}=(-17.8\pm 23.5,-7.9\pm 20.2)\mathrm{\mu }$as yr⁻¹, respectively. Again, the uncertainties in the 6p solution are approximately three to four times larger than those of the 5p solution.

Based on the above analysis, we find no significant color dependence in the zero-point offsets. However, a substantial systematic difference exists between the zero-points of the 5p and 6p solutions, which leads to the observed PM discrepancies between the blue and red samples, given their significantly different 5p and 6p fractions. Therefore, it is necessary to apply separate zero-point corrections for sources of each solution type. This correction will help further reduce the PM discrepancies observed between the blue and red samples in M31, as indicated below.

4.3. Proper motion of M31

Based on the above discussion, we now present our approach to resolving the discrepancy in M31 PM measurements between the blue and red samples, which arises from two issues associated with the 6p solution. We first revisited the PM results derived from the initial samples constructed in Sect. 3.1.1, as is shown in Fig. 5. Following Salomon et al. (2021), we applied the overall zero-point correction ${({\mu }_{\alpha }^{\ast },{\mu }_{\delta })}_{\mathrm{zpt}}^{\mathrm{M}31}$, as determined in Sect. 4.2, to both the red and blue samples. As is clearly shown in panel a of Fig. 5, a significant discrepancy persists between the two samples, in agreement with previous results of Salomon et al. (2021) and Rusterucci et al. (2024). To address this issue, and guided by the tests presented in Sect. 4.1, we first removed extremely red stars with G_BP − G_RP > 2.6 from the red 6p sample, as their PM measurements are likely unreliable. After this cut, the updated PM measurements from both the blue and red samples, shown in panel b of Fig. 5, are consistent within 2σ uncertainties, significantly reducing the discrepancy observed in the initial samples. Building on this color cut, we further derived the M31 PM values by applying zero-point corrections separately to the 5p and 6p solutions, using background QSOs matched to each solution type, as described in Sect. 4.2. The resulting measurements are shown in panel c of Fig. 5, where the PM values from the blue and red samples agree within 1σ uncertainties. Following the two-step refinement, the tension in M31 PM measurements between the blue and red samples is effectively resolved.

Fig. 5. Procedure for addressing the M31 PM discrepancy measured between the blue and red samples in the heliocentric frame assuming a distance of 784 kpc taken from Stanek & Garnavich (1998). Each panel is labeled in the topleft corner to indicate the corresponding sample. In panels (a)–(c), our PM results from the blue and red samples are connected by solid black lines. For comparison, we also show two other measurements based on Gaia (E)DR3 from Salomon et al. (2021) (gray squares) and Rusterucci et al. (2024) (gray triangles), connected by dashed gray lines. In panel (d), we compare our final results with previous measurements. Here, the black inverted triangle indicates a purely radial orbit of M31 relative to the MW. The filled magenta and open blue and red points show our new results from 5p, 6p blue, and 6p red samples (GBP − GRP ≤ 2.6), respectively. In addition to Gaia (E)DR3, other measurements include results from HST observations (Sohn et al. 2012; van der Marel et al. 2012a) (olive diamond) and Gaia DR2 (van der Marel et al. 2019) (black triangle).

Finally, we examined the results using both color and solution type: 5p (261 stars), blue 6p (551 stars), and red 6p stars with G_BP − G_RP ≤ 2.6 (620 stars). The 5p sample is not further divided by color, as red 5p sources are too scarce for meaningful analysis. After applying the zero-point corrections separately for the 5p and 6p solutions, we determined the corrected systemic PM and the corresponding transverse velocity in the heliocentric frame. These results are presented in the first three rows of Table 3 and shown in panel d of Fig. 5. Overall, the results from the three subsamples are consistent within 1σ, further demonstrating the effectiveness of our procedure in resolving the PM discrepancy between the blue and red samples. Nevertheless, the larger uncertainties associated with the 6p results compared to the 5p measurement suggest that the 5p sample, despite having fewer sources, yields more accurate and robust result due to its significantly smaller uncertainties. Compared with previous measurements, as also shown in panel (d), our 5p result is in excellent agreement with that obtained from HST (Sohn et al. 2012; van der Marel et al. 2012a), but offers improved precision.

After deriving the results, we provide a more intuitive illustration of the PM kinematic of M31 sample stars and its COM in the left panel of Fig. 6. The full sample consists of 1432 sources (5p stars and 6p stars with G_BP − G_RP ≤ 2.6), for which the observed PMs were corrected for zero-point offsets separately according to their solution type. The tracers were then divided into six bins in position angle, each containing at least 210 sources. The differences between the average observed PMs and the systemic PM derived from the 5p sample (first row in Table 3, denoted by the magenta segment at the left panel of Fig. 6) for these six bins exhibit a clear circular motion pattern across the disk, consistent with the predictions from the RC. This agreement further validates the robustness of our results.

Fig. 6. PM kinematics of M31 (left) and M33 (right). All 5p and 6p sample stars with GBP − GRP ≤ 2.6 are marked as red dots. The blue segments represent their predicted PMs from disk rotation, as determined in Sects. 3.1.1 and 3.1.2, without considering the systemic PM component. The derived M31 and M33 COM PMs (after zero-point correction) using the 5p sample in the heliocentric frame are shown as the magenta segments at the bottomleft corner of each panel. Green segments illustrate the average observed PMs in six (four) position angle bins for M31 (M33), after zero-point correction and subtraction of the COM PM. Each bin contains more than 210 (300) sources for M31 (M33), with bin boundaries indicated by dashed black lines. In each case, the PM direction starts at the dot and points along the line segment.

To transform the derived PMs and transverse velocity from the heliocentric frame to the galactocentric rest frame, and to determine the relative motion between M31 and the MW, we must correct for the solar reflex motion. In this work, we used the solar motion of (U_⊙,  V_⊙,  W_⊙) = (7.01,  10.13,  4.95) km s⁻¹ (Huang et al. 2015) and the circular speed at the solar position V_c(R₀) = 234.04 km s⁻¹ (Zhou et al. 2023). This leads to a reflex motion of ${({\mu }_{\alpha }^{\ast },{\mu }_{\delta })}_{\odot }=(36.4,-20.2)\mathrm{\mu }$as yr⁻¹ for M31 at 784 kpc, represented by the black inverted triangle in panel d of Fig. 5. If the derived PMs in the heliocentric frame match these values, it indicates zero transverse and a purely radial motion between M31 and the MW. After correcting for this solar reflex motion, we derived the PMs and corresponding transverse velocities in the galactocentric rest frame, which are provided in the last three rows of Table 3. The result from the 5p sample corresponds to a transverse velocity of V_trans = 46.7 ± 24.0 km s⁻¹, in agreement with recent cosmological simulations (Sawala et al. 2023). This supports the scenario of a nearly radial orbit for the relative motion between M31 and the MW, suggesting that the LG is gravitationally bound.

By combining the PM measurement from the 5p sample with the systemic LOS velocity of M31 derived in Sect. 3.1.1, both corrected for solar reflex motions, we derived the galactocentric 3D velocity of the M31 COM, ${\mathit{v}}_{\mathrm{DR}3}^{\mathrm{M}31}=(V_X,V_Y,V_Z)=(24.1\pm 25.9,-107.9\pm 17.4,38.8\pm 22.8)$ km s⁻¹. The uncertainties were propagated from the individual components using a Monte Carlo method.

Notably, the CMD distributions of our blue and red samples differ slightly from those used by Salomon et al. (2021). To assess the robustness of our results to sample selection, we reanalyzed an alternative CMD-selected sample following Salomon et al. (2021), applying the same quality cuts described in Sects. 2.2 and 3.1.1, together with our correction procedure. The PM results at each step show excellent agreement with those from our supergiant samples, confirming the robustness of our measurements and the reliability of our methodology in addressing the discrepancy between blue and red samples.

4.4. Proper motion of M33

Following a similar procedure as for M31, we next derived the systemic PM of M33 in the heliocentric frame. We began by remeasuring the PM using the initial sample constructed in Sect. 3.1.2, applying the overall zero-point correction ${({\mu }_{\alpha }^{\ast },{\mu }_{\delta })}_{\mathrm{zpt}}^{\mathrm{M}33}$ derived in Sect. 4.2. The results are shown in panel a of Fig. 7. The discrepancy between the blue and red samples is notably smaller than that observed in M31, particularly in the declination component, consistent with the findings of Rusterucci et al. (2024). This is primarily due to the limited number of extremely red supergiants with G_BP − G_RP > 2.6 in M33; that is, the PM measurements in M33 are less affected by the first issue of the 6p solution found in Sect. 4.1. For the right ascension component, the observed difference between the blue and red samples mainly arises from the systematic offsets between the 5p and 6p astrometric solutions, further amplified by the differing proportions of 5p and 6p stars in the two samples. Similar to M31, we first excluded these extremely red 6p stars, and the corresponding results are shown in panel b of Fig. 7. On top of this color cut, we further derived the systemic PM of M33 by applying separate zero-point corrections to the 5p and 6p stars, using background QSOs corresponding to each solution type, as described in Sect. 4.2. These results are shown in panel c of Fig. 7. Clearly, this step effectively mitigates the discrepancy between the blue and red samples in the right ascension component observed in panel (a). Finally, we determined the M33 COM PM using three subsamples: 5p (334 stars), blue 6p (502 stars), and red 6p stars with G_BP − G_RP ≤ 2.6 (570 stars). After applying the respective zero-point corrections for the 5p and 6p solutions, the final results are listed in the first three rows of Table 4 and illustrated in panel d of Fig. 7. Overall, the three measurements are consistent within 1σ, with the 5p result being the most robust and accurate due to its smaller uncertainties. Compared to previous studies, our 5p measurement agrees with the VLBA result within 1.5σ, with comparable uncertainties.

Fig. 7. Procedure for determining M33’s systemic PM in the heliocentric frame, assuming a distance of 840 kpc (Galleti et al. 2004b; Breuval et al. 2023). Each panel is labeled in the topleft corner to indicate the corresponding sample. The results at each step are indicated using the same symbols as in Fig. 5. For comparison, results from Rusterucci et al. (2024) based on Gaia DR3 are shown as gray triangles connected by dashed gray lines. In panel (d), we compare our final results with previous measurements. The blue inverted triangle indicates a purely radial orbit of M33 relative to the MW. The filled magenta and open blue and red points show our new results from 5p, 6p blue, and 6p red samples (GBP − GRP ≤ 2.6), respectively. The green cross represents the measurement by VLBA (Brunthaler et al. 2005) and the black triangle shows the result from Gaia DR2 (van der Marel et al. 2019).

The right panel of Fig. 6 displays a more intuitive PM kinematics of M33 sample stars and its COM. After applying zero-point corrections separately for 5p and 6p solutions, all 1406 sources (5p stars and 6p stars with G_BP − G_RP ≤ 2.6) were grouped into four position angle bins of over 300 sources each. The differences between the average observed PM values and the systemic PM derived from the 5p sample (first row in Table 4, denoted by the magenta segment at the right panel of Fig. 6) in these bins, exhibit a clear circular motion pattern. This aligns well with the predictions from the RC, thereby further confirming the reliability of our measurement.

To transform the derived PMs and transverse velocity from the heliocentric frame to the galactocentric rest frame, we then corrected for the solar reflex motion. Using the same solar motion parameters adopted for M31, we applied a reflex motion of ${({\mu }_{\alpha }^{\ast },{\mu }_{\delta })}_{\odot }=(38.2,-31.4)\mathrm{\mu }$as yr⁻¹ at the position of the M33 COM by adopting its distance of 840 kpc. After correcting for this solar reflx motion, the results are given in the last three rows of Table 4. Combined with the M33 systemic LOS velocity derived in Sect. 3.1.2, we obtained the galactocentric 3D velocity ${\mathit{v}}_{\mathrm{DR}3}^{\mathrm{M}33}=(-50.3\pm 30.6,76.3\pm 27.6,215.4\pm 24.5)$ km s⁻¹, with uncertainties estimated following the same procedure as in Sect. 4.3.

Fig. 8. Numerical orbital integration backward for M33 in the galactocentric frame centered on M31. The 3D galactocentric velocity of M31 is taken from this work. Orbits are integrated backward for 6 Gyr from the present-day position and velocity. Panels (a) and (b) show the orbit of M33 using its 3D velocity measured solely in this work, while panels (c) and (d) adopt the velocity from this work combined with the VLBA measurement. Similar to Patel et al. (2017) and van der Marel et al. (2019), we consider four mass combinations for M31 and M33, with two corresponding virial radii of M31: Rvir = 299 kpc (solid gray lines) for an M31 mass of 1.5 × 1012 M⊙, and Rvir = 329 kpc (dashed gray lines) for an M31 mass of 2 × 1012 M⊙. Panels (a) and (c) display the orbit of M33 projected onto the Y − Z plane, while panels (b) and (d) show the evolution of M33’s distance from M31 as a function of lookback time.

Based on the galactocentric 3D velocities of the M31 and M33 COM derived in this work, the relative motion between M33 and M31 can be described by a radial component of V_rad, DR3 = −169.7 km s⁻¹ and a transverse component of V_trans, DR3 = 204.6 km s⁻¹. For M33, if we combine the result in this work (5p sample) with that from VLBA (Brunthaler et al. 2005), a weighted average PM result in the heliocentric frame can be derived as ${({\mu }_{\alpha }^{\ast },{\mu }_{\delta })}_{\mathrm{DR}3+\mathrm{VLBA}}^{\mathrm{M}33}=(30.6\pm 5.7,19.0\pm 5.7)\mathrm{\mu }$as yr⁻¹. Correcting for the solar reflex motion, this corresponds to a galactocentric 3D velocity ${\mathit{v}}_{\mathrm{DR}3+\mathrm{VLBA}}^{\mathrm{M}33}=(1.5\pm 18.5,97.4\pm 17.9,181.8\pm 19.6)$ km s⁻¹. The relative motion between M33 and M31 then comprises a radial component of V_{rad, DR3 + VLBA} = −185.0 km s⁻¹ and a transverse component of V_{trans, DR3 + VLBA} = 170.0 km s⁻¹.

To investigate whether a first infall scenario is supported for M33 given the newly derived 3D velocities of M31 and M33, we followed a similar procedure to that of Patel et al. (2017) and van der Marel et al. (2019), numerically integrating M33’s orbit backward over 6 Gyr while neglecting the influence of the MW. The adopted potential models and parameters were identical to those used in Patel et al. (2017) and van der Marel et al. (2019). The distances to M31 and M33 were set to 784 and 840 kpc, respectively, and their velocities were updated using the values derived in this work. For M33, we considered two cases: using only our Gaia DR3-based result (${\mathit{v}}_{\mathrm{DR}3}^{\mathrm{M}33}$, shown in the left-hand panels of Fig. 8) and combining it with the VLBA result (${\mathit{v}}_{\mathrm{DR}3+\mathrm{VLBA}}^{\mathrm{M}33}$, shown in the right-hand panels of Fig. 8). In both cases, regardless of the assumed mass of M31 and M33, the results consistently support a first infall scenario; namely, that M33 is currently falling into M31’s virial radius for the first time.

5. Summary

In this work, we revisited the systemic PMs of M31 and M33 using massive supergiant stars as tracers. These tracers have been selected from color-color or color–magnitude diagrams, excluding foreground contaminants based on kinematic and astrometric properties, making them purer than those selected solely from CMD of Gaia broadband photometry. To further ensure sample purity, we applied several criteria to retain reliable disk objects. To assess the contribution of the supergiants’ disk rotation to the observed PMs, we established RC models for the disk of M31 and M33 using sources with LOS velocity measurements from the reliable disk objects selected earlier.

Using the astrometric measurements from Gaia DR3, we derived the systemic PMs for M31 and M33. For M31, the results from initial samples show significant discrepancy between blue and red samples, in line with previous studies. To address this, we identified its origin as arising from systematic astrometric offsets between sources with a 6p astrometric solution and those with a 5p one. The discrepancy was further amplified by the differing fractions of 5p and 6p sources in the blue and red samples. The red sample is dominated by 6p sources while the blue sample includes a substantial fraction of 5p sources. Specifically, we identified two main issues associated with 6p sources: (1) unreliable astrometric measurements for extremely red stars, specifically those with G_BP − G_RP > 2.6, which exhibit significant deviations from the bulk of the 6p population. Combined with the findings of Fabricius et al. (2021), we conclude that the astrometry of these red 6p sources is likely unreliable; (2) a color-independent systematic offset between 5p and 6p solutions, which is caused by the different zero-point offsets between these two solution types. To correct for these effects, we applied the following steps: (1) excluding 6p sources with extreme red colors by restricting the sample to those with G_BP − G_RP ≤ 2.6; (2) correcting the PM zero-point offsets separately for 5p and 6p sources using background QSOs, thereby accounting for systematic differences between the two solution types. After applying these corrections, the PM measurements derived from M31 blue and red samples become consistent within 1σ, and both indicate a radially dominated orbit of M31 relative to the MW. This resolves the previously found discrepancy between PM measurements using blue and red samples and provides robust evidence of a radially dominated M31–MW relative motion.

Using 5p sample stars, which provide the most accurate and robust result, we derived a galactocentric transverse velocity of V_trans = 46.7 ± 24.0 km s⁻¹ for M31. This is in excellent agreement with previous HST-based measurement, but with improved precision. For M33, our measurement is consistent with previous VLBA measurements within 1.5σ. Combining both M31 and M33 measurements, M33’s first infall scenario is supported after numerical integration, regardless of the mass model adopted.

Overall, we present the most robust Gaia-based PM measurements for M31 and M33 to date, conclusively demonstrating that the relative motion between M31 and the MW is radially dominated. Our work also highlights the systematic offsets between Gaia 5p and 6p astrometric solutions, as well as the biases observed in extreme red 6p sources. Given that approximately twothirds of the sources in Gaia DR3 lack 5p solution, especially significant for faint sources (the fraction of 5p sources decreases from 88% at G < 18 to only 1.4% at 20 < G < 21; Lindegren et al. 2021b), we recommend additional caution when using the 6p solution in cases requiring high-precision astrometry. This work also provides valuable insights that may help mitigate these issues in the upcoming Gaia DR4.

Acknowledgments

We are grateful to Prof. Haibo Yuan for his invaluable comments and insightful discussions. This work acknowledges the supports from National Key R&D Programme of China (Grant No. 2024YFA1611903) and the National Science Foundation of China (NSFC Grant No. 12422303, 12090040 and 12090044). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This work has made use of data products from the Guo Shou Jing Telescope (the Large Sky Area Multi-Object Fibre Spectroscopic Telescope, LAMOST). LAMOST is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences.

References

All Tables

All Figures

Fig. 1. Panel (a): Spatial distribution of 321 supergiant candidates for modeling the RC of M31 disk region. The red star denotes the COM of M31, and the dashed black outer line represents the 1.8° ellipse of the projected disk. Inner and outer green ellipses denote the reliable range of our derived RC. Panel (b): RC of the M31 disk region. The orange line and the shaded area represent the best fit and the corresponding standard error from this work, while the blue and green error-bars represent results from previous studies, as labeled in the bottomright corner. Panel (c): Spatial distribution of 235 supergiant candidates for modeling the RC of M33 disk region. The red star denotes the COM of M33, and the dashed black outer line represents the 1° ellipse of the projected disk. Inner and outer green ellipses denote the reliable range of our derived RC. Panel (d): RC of the disk of M33. The orange line and the shaded area illustrate the best fit and its corresponding standard error from this work. The blue and green error-bars represent results obtained from previous studies by H I observations, as labeled in the bottomright corner.

In the text

Fig. 2. Spatial and CMD distributions of final M31 and M33 supergiant samples (see Sect. 3.1). The left-hand panels display the spatial and CMD distributions of the M31 sample stars, while the right-hand panels present those of the M33 sample stars. Filled blue and red symbols denote blue (including some yellow) and red supergiants with a 5p astrometric solution, respectively; open cyan and pink circles indicate blue and red supergiants with a 6p solution. The dashed magenta line in the CMD panels marks GBP − GRP = 2.6. Clearly, M31 contains a significantly larger population of extremely red supergiant stars compared to M33.

In the text

Fig. 3. Left-hand panels: PM distribution of the blue and red M31 samples as a function of GBP − GRP. The data points in the left (right) panel show the weighted average differences between the observed PMs, ${\mu }_{\alpha ,\mathrm{obs}}^{\ast }$ (μδ, obs), from Gaia DR3 and the predicted rotational components ${\mu }_{\alpha ,\mathrm{RC}}^{\ast }$ (${\mu }_{\delta ,\mathrm{RC}}$), derived from the RC in Sect. 3.1.1. Weights are proportional to 1/σα*2 or 1/σδ2, and the errors, also weighted by uncertainty, are indicated accordingly. Right-hand panels: Similar to the left-hand panels, but for M33.

In the text

Fig. 4. Left-hand panels: Inferred PM values of M31 background QSOs in right ascension (${\mu }_{\alpha }^{\ast }$; left) and declination (μδ; right) as a function of GBP − GRP color. The upper sub-panels show the inferred PMs in each color bin for the full QSO sample (gray points), with the dashed gray line and shaded region indicating the overall zero-point offset and its 1σ uncertainty from the MCMC analysis. The lower sub-panels show results for QSOs with different astrometric solutions: blue and orange points correspond to the inferred PMs in each color bin for 5p and 6p sources, respectively. Dashed lines and shaded areas denote the overall PMs and their 1σ uncertainties for the two subsamples. Right-hand panels: Same as the left-hand panels, but for M33 background QSOs.

In the text

Fig. 5. Procedure for addressing the M31 PM discrepancy measured between the blue and red samples in the heliocentric frame assuming a distance of 784 kpc taken from Stanek & Garnavich (1998). Each panel is labeled in the topleft corner to indicate the corresponding sample. In panels (a)–(c), our PM results from the blue and red samples are connected by solid black lines. For comparison, we also show two other measurements based on Gaia (E)DR3 from Salomon et al. (2021) (gray squares) and Rusterucci et al. (2024) (gray triangles), connected by dashed gray lines. In panel (d), we compare our final results with previous measurements. Here, the black inverted triangle indicates a purely radial orbit of M31 relative to the MW. The filled magenta and open blue and red points show our new results from 5p, 6p blue, and 6p red samples (GBP − GRP ≤ 2.6), respectively. In addition to Gaia (E)DR3, other measurements include results from HST observations (Sohn et al. 2012; van der Marel et al. 2012a) (olive diamond) and Gaia DR2 (van der Marel et al. 2019) (black triangle).

In the text

Fig. 6. PM kinematics of M31 (left) and M33 (right). All 5p and 6p sample stars with GBP − GRP ≤ 2.6 are marked as red dots. The blue segments represent their predicted PMs from disk rotation, as determined in Sects. 3.1.1 and 3.1.2, without considering the systemic PM component. The derived M31 and M33 COM PMs (after zero-point correction) using the 5p sample in the heliocentric frame are shown as the magenta segments at the bottomleft corner of each panel. Green segments illustrate the average observed PMs in six (four) position angle bins for M31 (M33), after zero-point correction and subtraction of the COM PM. Each bin contains more than 210 (300) sources for M31 (M33), with bin boundaries indicated by dashed black lines. In each case, the PM direction starts at the dot and points along the line segment.

In the text

Fig. 7. Procedure for determining M33’s systemic PM in the heliocentric frame, assuming a distance of 840 kpc (Galleti et al. 2004b; Breuval et al. 2023). Each panel is labeled in the topleft corner to indicate the corresponding sample. The results at each step are indicated using the same symbols as in Fig. 5. For comparison, results from Rusterucci et al. (2024) based on Gaia DR3 are shown as gray triangles connected by dashed gray lines. In panel (d), we compare our final results with previous measurements. The blue inverted triangle indicates a purely radial orbit of M33 relative to the MW. The filled magenta and open blue and red points show our new results from 5p, 6p blue, and 6p red samples (GBP − GRP ≤ 2.6), respectively. The green cross represents the measurement by VLBA (Brunthaler et al. 2005) and the black triangle shows the result from Gaia DR2 (van der Marel et al. 2019).

In the text

Fig. 8. Numerical orbital integration backward for M33 in the galactocentric frame centered on M31. The 3D galactocentric velocity of M31 is taken from this work. Orbits are integrated backward for 6 Gyr from the present-day position and velocity. Panels (a) and (b) show the orbit of M33 using its 3D velocity measured solely in this work, while panels (c) and (d) adopt the velocity from this work combined with the VLBA measurement. Similar to Patel et al. (2017) and van der Marel et al. (2019), we consider four mass combinations for M31 and M33, with two corresponding virial radii of M31: Rvir = 299 kpc (solid gray lines) for an M31 mass of 1.5 × 1012 M⊙, and Rvir = 329 kpc (dashed gray lines) for an M31 mass of 2 × 1012 M⊙. Panels (a) and (c) display the orbit of M33 projected onto the Y − Z plane, while panels (b) and (d) show the evolution of M33’s distance from M31 as a function of lookback time.

In the text