The M31’s apparent size in the sky is about 3°10′ by 1°, therefore the factors such as kinematics, position of field, and more need to be considered for high precision astrometry. We use the analytical velocity field model presented by van der Marel et al. (2002) to describe PM field, because these equations are derived through mathematical deductions without other assumptions. For the planar system, the velocity of a tracer can be decomposed into the following three components 1: V = V CM + V pn + V int (1) where V _CM is contributed by the space motion of center of M31, V _pn is contributed by precession and nutation, and it is small and is neglected, V _int is contributed by the internalrotation motion.

There are three coordinate systems as shown in Figure 3. The first Cartesian coordinate system (x, y, z): origin at the center of mass of M31, with the x-axis antiparallel to the right ascension axis, the y-axis parallel to the declination axis, and the z-axis toward the observer. The second Cartesian coordinate system (x′, y′, z′): origin at the center of mass of M31, with the (x′, y′)-plane is inclined with respect to the (x, y)-plane by the inclination angle i, the x′-axis along the intersection of the M31 disk plane and (x, y)-plane with position angle θ between the x-axis. The third Left-handed Cartesian coordinate system (x _l, y _l, z _l): origin at the center of the FOV, with z _l-axis along the line of sight, x _l along the tangent to the great circle at the origin, and the great circle pass the center of the M31 to the origin.

The Eq. 1 can be parameterized as follows in (2) and (3): V CM = v z l v x l v y l CM = v t ⁡ sin ⁡ ρ ⁡ cos ϕ − θ t + v sys ⁡ cos ⁡ ρ v t ⁡ cos ⁡ ρ ⁡ cos ϕ − θ t − v sys ⁡ sin ⁡ ρ − v t ⁡ sin ϕ − θ t V int = v z l v x l v y l int = s V R ′ cos 2 ⁡ i ⁡ cos 2 ϕ − θ + sin 2 ϕ − θ 1 / 2 × − sin ⁡ i ⁡ cos ϕ − θ cos ⁡ i ⁡ cos ⁡ ρ − sin ⁡ i ⁡ sin ⁡ ρ ⁡ sin ϕ − θ sin ⁡ i ⁡ cos ϕ − θ cos ⁡ i ⁡ sin ⁡ ρ + sin ⁡ i ⁡ cos ⁡ ρ ⁡ sin ϕ − θ − cos 2 ⁡ i ⁡ cos 2 ϕ − θ + sin 2 ϕ − θ (2) μ W μ N = cos ⁡ i ⁡ cos ⁡ ρ − sin ⁡ i ⁡ sin ⁡ ρ ⁡ sin ϕ − θ D 0 ⁡ cos ⁡ i × − sin Γ − cos Γ cos Γ − sin Γ v x l v y l (3)

These parameters in these equations are divided into four groups. The first group is the position of the field with respect to the center of M31: ρ is angular distance, ϕ is position angle. The second group describes the direction from which the plane of M31 is viewed: i is the inclination angle, θ is the position angle of the line of nodes. The third group describes the motion information of the center of M31: v _sys is the velocity along the line of sight, v _t and θ _t constitute the transverse velocity vector, D ₀ is the distance to observer. The fourth group describes internal rotation: V(R′) is rotation velocity and s = ±1 determines the direction of rotation. Other quatities: μ W = − μ α * , μ _N = μ _δ , Γ is position angle of v x l relative to the direction of positive declination.

For the simulated data, the above parameters can be categorized into two types, the first type is fixed value as follows: (i, θ) is (73.7°,128.3°), D ₀ is 785 kpc, v _sys is 301 km/s (Karachentsev et al., 2004; McConnachie et al., 2005; Chemin et al., 2009). The second type will be fitted, but with initial values given to generate the simulated data, as follows: v _t is 160 km/s, θ _t is −2 rad and V(R′) is −230 km/s.

The stellar distribution of the disk of M31 is assumed to fit a classical surface density profile of the galactic disk, i.e., (Yin et al., 2009; Courteau et al., 2011), it can be represented by an exponential Eq. 4: M r = M 0 e − r / r d (4) where M ₀ is the stellar density at the galactic center; r is the distance from the galactic center; r _d is the disk scalelength; M _r is the stellar density at the position with r distance from the galactic center. Given the simulation purpose, a typical scalelength, 5.5 kpc (Yin et al., 2009), is chosen to approximate the distribution of M31 member stars.

We assume that the spiral galaxy’s disk has an axisymmetric potential, and rotation curve is a function of disk stars’ circular velocity and distance from the disk center. The M31 rotation curve is obtained from Chemin et al. (2009) to approximate the circular motion of the M31 disk stars.

Eight simulated observed fields, numbered from 1 to 8, are shown in Figure 4. Their celestial coordinates are (00h49m08.6s, 42°45′01.6″), (00h44m54.9s 40°57′13.1″), (00h36m37.4s 39°45′58.7″), (00h40m23.8s 41°36′06.3″), (00h48m15.4s, 41°51′54.0″), (00h41m19.4s 40°11′42.0″), (00h36m44.4s 40°26′06.0″), (00h45m00.0s 42°28′48.0″), respectively. The reason for these fields: (1) the field 1, the HST-observed field, can be compared. (2) they are located at approximately the same physical distance from the center of M31 disk, ensuring consistency in the rotational speed.

The Gaia DR3 reference stars are selected for two reasons: First, they are point sources that offer higher astrometric precision, whereas extended sources such as background galaxies are limited by various shapes and sizes. Second, there are a large number of them, unlike the sparsely distributed quasars. There are 1288 QSOs in the 2° region centered on M31 (Liao et al., 2019), and the number in the simulated observed field is 1 or 2 and is approximately 1/100 of the Gaia reference stars. The advantage in terms of number (N) leads to an improvement in astrometric precision proportional to N . When using Gaia catalog as the reference stars, the positional errors from the catalog are almost the same for M31 member stars at any epoch, so the errors do not dominate the error budget of PM result of M31 member stars. However, the PM errors from the catalog are cumulative for M31 member stars and are retained. Therefore, the PM errors have a significant impact on the precision of M31 PMs. It is necessary to modify the PMs of the catalog using observational images and then calculate the plate model parameters. At an observed epoch t, the pixel coordinates x(t), y(t) of reference stars are generated from the positions, positional errors, PMs, PM errors, parallaxes of the Gaia DR3 catalog, telescope parameters, etc. The celestial coordinates α _C(t), δ _C(t) of reference stars are extrapolated using the positions, PMs and parallaxes from the Gaia DR3 catalog, where the subscript “C” stands for catalog. The standard coordinates ξ _C(t), η _C(t) are obtained by performing a gnomonic projection equation of α _C(t), δ _C(t), as defined in Eq. 5. ξ = cos ⁡ δ ⁡ sin α − α 0 sin δ 0 ⁡ sin ⁡ δ + cos δ 0 ⁡ cos ⁡ δ ⁡ cos α − α 0 η = cos δ 0 ⁡ sin ⁡ δ − sin δ 0 ⁡ cos ⁡ δ ⁡ cos α − α 0 sin δ 0 ⁡ sin ⁡ δ + cos δ 0 ⁡ cos ⁡ δ ⁡ cos α − α 0 (5) where α, δ denotes celestial spherical coordinates, with the subscript ’0’ (α ₀, δ ₀) meaning the celestial position of the plate center or “tangent-point” of the Cartesian coordinates, in which the η-axis is parallel to the declination axis and the ξ-axis is parallel to right ascension axis.

Then the plate model can be calculated by Eq. 6 ξ C t = a x t + b y t + c η C t = d x t + e y t + f (6) where a, b, c, d, e, f is the affine transformation parameters between standard coordinates (ξ, η) and pixel coordinates (x, y), also known as the plate model. This plate model can be applied to the x(t), y(t) to obtain α _O(t), δ _O(t), where subscript “O” stands for observation. And, besides the simulated epochs, J2030.0 to J2039.0, the J2016.0 epoch from the Gaia DR3 catalog can be an additional observed epoch for modifying the PMs of the catalog. The PMs and errors, μ α * O , μ δ O , Δ μ α * O , Δ μ δ O , are obtained by solving Eq. 7 (Perryman et al., 1997) using the least squares method. ξ O t = μ α * O t − t r e f − p r e f b O t π / A 1 − r r e f b O t π / A + Δ α * η O t = μ δ O t − t r e f − q r e f b O t π / A 1 − r r e f b O t π / A + Δ δ r r e f = cos δ r e f ⁡ cos α r e f cos δ r e f ⁡ sin α r e f sin δ r e f p r e f = − sin α r e f cos α r e f 0 q r e f = − sin δ r e f ⁡ cos α r e f − sin δ r e f ⁡ sin α r e f cos δ r e f (7) where subscript “ref” and “O” stands for reference and observation respectively; t _ref is the reference epoch; b O t is the barycentric position of observer; π is parallax; A is astronomical unit; r _ref, p _ref and q _ref are orthogonal unit vectors, with original point (α _ref, δ _ref); Δα* = Δα cos δ _ref and Δδ are the residuals at the reference epoch in α _ref and δ _ref.

Stationary reference sources, such as QSOs and background galaxies, can be used to align images from different observed epochs. To transform moving Gaia DR3 reference stars into “stationary”, we adopt the approach as shown in Figure 5: Select an epoch as reference epoch and the corresponding image as the reference image. Here, the intermediate epoch (J2035.0) is reference epoch, and represented by “mid”. The pixel coordinates of Gaia stars at other epochs can be adjusted to the intermediate epoch by applying the modified PMs Eq. 8. x t m i d = x t + t m i d − t v x y t m i d = y t + t m i d − t v y (8) where t _mid is intermediate epoch; v _x , v _y are stellar velocities in the pixel coordinates; x(t)_mid, y(t)_mid are adjusted pixel coordinates. This adjustment ensures that each image at different epochs can be aligned with the reference image throught “stationary” Gaia stars, and the transformation as shown in Eq. 9. x t m i d = a x t m i d + b y t m i d + c y t m i d = d x t m i d + e y t m i d + f (9) where a, b, c, d, e, f are transformation parameters between images of one epoch and intermediate epoch. For M31 member stars in each image, applying affine transformation and the plate model at the intermediate epoch will convert pixel coordinates to celestial coordinates. Subsequently, the PMs of the M31 member stars can be calculated.

The adopted theoretical PMs of member stars within the first simulated observed field are calculated from J2030.0 to J2039.0, as shown in Figure 6. We assume that all member stars in the field have the same PM is unreasonable due to the significant degree of dispersion (∼3.0 μas/yr) relative to the expected observed PMs’ precision (∼10.0 μas/yr), and the field is divided into 2 × 2 blocks for reducing the degree. The distribution of one of the blocks is shown in Figure 7, where the dispersion of member stars is reduced to 1.5 μas/yr.

The hypothetical strategy is to observe the all 8 fields quasi-simultaneously from J2030.0 to J2039.0, and the number of observed epochs per year is 1, 5, 10, 20 and 50. Figure 8 shows the relationship between the uncertainties of the average PMs within the first field and the number of observed epochs per year. It can be seen that uncertainty of 10 μas/yr is achievable with 10 epochs per year, and of 5 μas/yr with 50 epochs.

The uncertainty decreases with the increase in the number of epochs per year. But the improvement effect also keeps decreasing, 5 to 10 epochs achieves a good level of precision. There is simple Eq. 10 to estimate the uncertainty, which predicts only the Poisson noise’s behaviour and ignores the unknown systematic errors. δ ≈ σ r e f 2 / N r e f + σ member 2 / N member / N epoch * 2 / T span (10) where σ _ref, σ _member are positional uncertainties for reference stars and M31 member stars; N _ref, N _member are corresponding star numbers; N _epoch is number of epochs per year; T _span is time span of PM calculation and 2 means displacement. Corresponding to 10 observed epochs in Figure 8, the values as follows: σ _ref, σ _member are 2.5 mas, 7.5 mas; N _ref, N _member are 165, ≈ 15000; N _epoch is 10; T _span is 9. And the result is 10.1 mas.

Using the simulated data provided above, the simulated PMs under 10 observed epochs per year are compared with the adopted theoretical values, as shown in Table 2. The differences between adopted theoretical and simulated observed values are introduced by PM errors in Gaia catalog.

From Eqs 2, 3 in the velocity field model, the simulated observed PM results (μ _W , μ _N ) or ( μ α * , μ δ ) can be used to fit the information about M31’s kinematics. The parameters (i, θ, D ₀, v _sys) are known, and parameters (ρ, ϕ) can be calculated based on the celestial coordinates of the field and M31’s center. The motion of M31’s center (v _t , θ _t ), and internal rotational velocity (V(R′)) are obtained through the least-squares fit of PM results, as shown in Table 3.

Analyzing the differences between adopted theoretical values v _t , θ _t , V(R′) (160 km/s, −2 rad, −230 km/s) and all fitting values, it could be found that these systematic errors are introduced by the PM errors in the Gaia catalog. Analyzing the differences in precision among fitted values from different fields, it could be found that with an increasing number of fields (1 to 8 fields), the precision improves. Under the condition of the same number of fields, if the fields are distributed on the same side (1, 5, 2, 6 fields or 3, 7, 4, 8 fields), the results are inferior to the symmetrically distributed fields. Under the condition where the observed fields are all symmetrically distributed, the fields away from the two sides of M31 (1, 3 fields) show better results in the internal rotation V(R′), and this may be the internal rotational motion contributes to the line-of-sight velocity at both ends. Furthermore, results of 5, 6, 7, 8 fields and 2, 6, 4, 8 fields show that having fields more dispersed might yield better results.