© The Authors 2025

1. Introduction

A century ago, astronomers were aware of only two prominent satellite galaxies orbiting the Milky Way, the Large Magellanic Cloud (LMC) and the Small Magellanic Cloud (SMC), a view that had persisted for centuries (von Humboldt & Otte 1997). This census doubled when Shapley (1938) used photographic plates to discover two additional dwarfs, Sculptor and Fornax. With the discovery of Sagittarius in 1994, the number of known Milky Way satellite galaxies reached 11 (Ibata et al. 1994), and has since increased substantially due to advances in observational technology, the deployment of larger telescopes, and systematic sky surveys. Currently, approximately 60 dwarf galaxies have been discovered around the Milky Way (Simon 2019). Notably, these bright satellite galaxies tend to lie along a great circle known as the “Magellanic Plane” (Lynden-Bell 1976). Observations show that Milky Way satellite galaxies are highly anisotropic and kinematically correlated, with most located in a thin plane roughly perpendicular to the disk of the Milky Way (Kroupa et al. 2005; Metz et al. 2008). The most luminous among these galaxies appear to align within a thin, disk-like structure known as the Vast Polar Structure (VPOS; Pawlowski et al. 2012; Taibi et al. 2024). The VPOS exhibits a root mean square (rms) thickness of 20–30 kpc, a radius of about 250 kpc, and an axis ratio c/a of 0.18–0.30 (Pawlowski 2018, 2021b).

Similarly, a comparable thin and anisotropic distribution of satellite galaxies is also found around Andromeda, referred to as the Great Plane of Andromeda (GPoA) (Ibata et al. 2013; Conn et al. 2013), and around Centaurus A, known as the Centaurus A Satellite Plane (CASP) (Tully et al. 2015; Müller et al. 2018; Kanehisa et al. 2023b). These structures are characterized by an axis ratio of c/a = 0.1 (Pawlowski et al. 2013) and a c/a ratio = 0.2 (Tully et al. 2015), respectively. Both the GPoA and CASP planes exhibit significant correlations in the line-of-sight velocities of their satellite members when observed edge-on (Ibata et al. 2013; Müller et al. 2021). This suggests that satellites corotate, with satellites on one side of the host predominantly receding and those on the other side approaching relative to the host. This rotating pattern around its host galaxy, similar to the VPOS, supports the idea that both the GPoA and CASP constitute a nonrandom, coherent structures. Moreover, the first proper motion measurements for three on-plane satellite galaxies of M31 indicate that these galaxies are consistent with co-orbiting along the spatially identified GPoA (Sohn et al. 2020; Pawlowski & Tony Sohn 2021; Casetti-Dinescu et al. 2024). A similar flattened structure with coherent line-of-sight velocities was discovered around the galaxy NGC 4490 and is found to be rare in cosmological simulations (Karachentsev & Kroupa 2024; Pawlowski et al. 2024).

Today, the Lambda cold dark matter (ΛCDM) cosmological model is widely accepted among researchers due to its ability to explain a range of astrophysical phenomena such as Big Bang nucleosynthesis, the accelerated expansion of the Universe (Riess et al. 1998), the power spectrum of the cosmic microwave background (Spergel et al. 2003), large-scale structure (Peacock et al. 2001; Tegmark et al. 2004), and the matter-energy budget of our Universe (Ade et al. 2016). However, on smaller scales, it presents several challenges. Among these, the plane of satellite galaxies has consistently remained a topic of debate (Pawlowski 2018), and such structures are considered highly improbable within the standard ΛCDM cosmological model (Kroupa et al. 2005; Pawlowski 2021b). Several solutions have been proposed to explain this discrepancy; however, no general consensus has been reached. These proposed solutions include accretion along filaments of the cosmic web (Libeskind et al. 2011; Lovell et al. 2011), infall of satellite galaxies in groups (Li & Helmi 2008; D’Onghia & Lake 2008; Metz et al. 2009; Shao et al. 2018; Wang et al. 2013; Júlio et al. 2024), galaxy mergers (Smith et al. 2016; Kanehisa et al. 2023a), effects induced by the infall of a massive LMC (Garavito-Camargo et al. 2021; Pawlowski et al. 2022; Vasiliev 2023), tidal dwarf galaxies (Pawlowski et al. 2011; Fouquet et al. 2012; Banik et al. 2022), and hydrodynamics and baryonic physics (Pawlowski et al. 2015; Ahmed et al. 2017; Müller et al. 2018; Pawlowski et al. 2019). Thus, while each hypothesis contributes valuable insight, a unified explanation remains elusive, highlighting the complexity of galaxy formation and evolution processes.

Precise investigation of the orbits of satellite galaxies requires a thorough understanding of several critical parameters. These parameters encompass precise measurements of a satellite’s position, distance, line-of-sight velocity, proper motions, the underlying gravitational potential of the host galaxy, the Galactocentric distance of the Sun, and its motion relative to the Galactic center. However, many of these key parameters remain inadequately constrained, introducing a significant uncertainty in orbital analyses (Bland-Hawthorn & Gerhard 2016; Joshi 2007; Schönrich et al. 2010). Measurement uncertainties, especially in proper motions, can significantly affect analyses (Pawlowski 2021a). These uncertainties can be broadly classified as statistical or systematic (Li et al. 2021). Statistical uncertainties arising from random measurement variations can be mitigated through statistical methods and the use of larger datasets. Conversely, systematic uncertainties, which result from inherent biases in the measurement process or instruments, pose more fundamental challenges.

Maji et al. (2017) demonstrate that the Milky Way’s disk of satellites (DoS) is dynamically transient. Their orbit modeling shows that the 11 classical satellites move away from the current DoS, causing it to thicken from c/a ∼ 0.18 (height ∼19.6 kpc) to c/a ∼ 0.36 (45 kpc) in 0.5 Gyr, and c/a ∼ 0.42 (64 kpc) in 1 Gyr. This evolution persists across the two tested galactic potentials: one replacing the stellar disk with a single-component Miyamoto–Nagai potential and another using only an Navarro–Frenk–White (NFW) dark matter halo potential. However, Pawlowski et al. (2017) notes that Maji et al. (2017) overlook measurement uncertainties, such as those related to proper motions, and argues that ignoring observational biases further undermines the reliability of their conclusions regarding the plane of satellite galaxies. Similarly, Sawala et al. (2023) argues that the plane of satellites is more likely to be transient rather than rotationally supported, because backward and forward integration of satellite galaxy orbits results in a widening of the inferred plane thickness. However, they do not account for systematic errors in Gaia proper motions. These uncertainties are crucial because they can significantly affect the measured values and potentially alter the results. Without accounting for these uncertainties, the conclusions drawn about the stability and alignment of the Milky Way’s satellite plane might be overly optimistic or pessimistic. This could lead to misinterpretations about the dynamical status of these satellite structures, thereby influencing our understanding of the formation and evolution of the Milky Way and its satellite system. Incorporating measurement uncertainties, particularly in proper motion data, is essential to achieve a more realistic and reliable understanding of the satellite galaxies’ orbital behaviors. This is especially true when studying the joint behavior of an ensemble of satellite galaxies, since errors can conceal out intrinsic correlations (and approaches such as Monte Carlo sampling from the uncertainties effectively apply the errors twice; Pawlowski 2021a).

To address this, we explore how different levels and types of measurement uncertainties affect the inferred stability of planes of satellite galaxies. Specifically, we employed computer simulations to model the orbital evolution of test satellite galaxies around a host galaxy that were set up as intrinsically stable satellite planes. We mock-observed these by applying controlled degrees of observational errors on their evolved positions and velocities, and then backward-integrated them to determine how well the intrinsically stable dynamical evolution can be inferred.

This paper is structured as follows. In Sect. 1, we describe the methodology, including setting up an artificial satellite plane, the plane-fitting routine, the orbit integrations performed, the application of mock measurement errors, and the subsequent backward integration of the test satellites. Section 3 presents our results for various combinations of measurement errors, assumed Milky Way potentials, and initial orbital eccentricity distributions. We end with a discussion and conclusions in Sect. 4.

2. Methodology

2.1. Plane setup

To study the evolution of artificial planes of satellite galaxies, the test satellites must be appropriately arranged around the host galaxy to initially form a planar structure. To achieve this, we used the Milky Way and its satellite distribution as a model to establish an artificial plane of satellite galaxies. The Milky Way satellite plane, VPOS, spans approximately 250 kpc in radius, with a root mean square thickness ranging from 20 to 30 kpc. Therefore, the radial distribution, R, of N_sat test satellites ranges from 20 to 250 kpc, while the vertical distribution extends from 0 to 20 kpc (Pawlowski 2018, 2021b). To generate a random radial distribution, we used the Von Neumann rejection algorithm, also referred to as rejection sampling (Neumann 1951). This statistical technique facilitates the generation of random samples from a desired probability distribution. Figure 1 shows the face-on and edge-on views of the randomly distributed N_sat = 25 test satellites around a host galaxy, initially set up as a planar structure. Here, the satellite plane is set perpendicular to the host galaxy’s stellar disk. To produce the radial distance distribution, we used a power-law distribution function:

Fig. 1. Face-on (left) and edge-on (right) views of Nsat = 25 randomly generated test satellites, each represented by a different color. The star symbol denotes the host galaxy.

$$\begin{array}{c}\hfill f(r)=r^{\alpha }.\end{array}$$(1)

Here, for α, we used −3. For the velocities, we first defined the unit vectors for each test satellite in the radial, tangential, and perpendicular directions, with respect to the satellite plane. In vector representation, these are written as

$$\begin{array}{cc}\hfill {\mathbf{V}}_{\text{rad; i}}& =|{\mathbf{V}}_{\text{rad; i}}|·{\widehat{r}}_i\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill {\mathbf{V}}_{\text{tan; i}}& =|{\mathbf{V}}_{\text{tan; i}}|·{\widehat{t}}_i\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill {\mathbf{V}}_{\text{perp; i}}& =|{\mathbf{V}}_{\text{perp; i}}|·{\widehat{p}}_i,\hfill \end{array}$$(4)

where V_{rad; i}, V_{tan; i}, and V_{perp; i} represent the radial, tangential, and perpendicular velocity vectors of the i-th test satellite, respectively. The quantities |V_{rad; i}|, |V_{tan; i}|, and |V_{perp; i}| denote the magnitudes of the radial, tangential, and perpendicular velocity, respectively. The unit vectors ${\widehat{r}}_i$, ${\widehat{t}}_i$, and ${\widehat{p}}_i$ correspond to the radial, tangential, and perpendicular directions for the i-th test satellite.

Hammer et al. (2021) and Cautun & Frenk (2017) investigated the relations between radial and tangential velocities of Milky Way satellite galaxies and report a tangential bias. Specifically, Cautun & Frenk (2017) show that nine out of ten satellites with measured proper motions exhibit highly tangentially biased motions, with 80% or more of their orbital kinetic energy attributed to tangential motion. Given the significance of tangential velocity, we initially set both |V_{perp; i}| and |V_{rad; i}| to zero, indicating no motion in the radial and perpendicular directions. Thus, only the tangential component is relevant. Satellites were initialized with a tangential speed equal to the circular velocity of the MW_fiducial gravitational potential (see Table 1), guaranteeing circular orbits. To make the orbits eccentric, we rotated the tangential unit vectors of all satellites by an angle θ around the perpendicular unit vector. For each test satellite, the values for θ were drawn from a uniform distribution with a maximum range of −θ_max and +θ_max, as listed in Table 2.

Once we obtained the radial V_{rad; i}, perpendicular V_{perp; i}, and tangential V_{tan; i} velocity components, we converted them into a Cartesian coordinate system. This allowed us to fully describe the motion of the test satellites in the simulated galaxy system.

2.2. Plane fitting

Once the test satellites were distributed around their host galaxy, the next step was to perform plane fitting to analyze their spatial arrangement. This technique allowed us to determine whether the satellites were arranged in a planar structure or scattered more isotropically around the host galaxy. For this purpose, we used an unweighted plane-fitting technique Metz et al. (2007). This method begins with the calculation of the moment-of-inertia tensor of the satellite distribution, which is followed by its diagonalization. The first step of the process is to determine the centroid r₀, of the data points:

$$\begin{array}{c}\hfill r_0=\frac{1}{N}{\displaystyle \sum _{i=1}^N}{r}_i.\end{array}$$(5)

Subsequently, an eigenvalue analysis of the moment of inertia tensor, T₀, is performed for the position vectors, ${\widehat{r}}_i=r_i-r_0,i=1\dots n$ (Metz et al. 2007; Pawlowski et al. 2013):

$$\begin{array}{c}\hfill T_0={\displaystyle \sum _{i=1}^N}[{({\mathbf{r}}_i-{\mathbf{r}}_0)}^2·\mathbf{1}-({\mathbf{r}}_i-{\mathbf{r}}_0)·{({\mathbf{r}}_i-{\mathbf{r}}_0)}^T].\end{array}$$(6)

Here, 1 denotes the unit matrix and r^T is the transposed version of the vector r. The square root of the eigenvalues of the moment of inertia indicates the extent along three axes (a, b, c) of the fitted ellipsoid to the satellite distribution. These values are proportional to the rms deviation relative to the eigenvectors of T. The eigenvector corresponding to the largest eigenvalue defines the normal of the plane, encompassing the centroid. Meanwhile, the eigenvectors corresponding to theintermediate and smallest eigenvalues indicate the directions of the intermediate and major axes of the distribution, respectively. Consequently, we determined the axial ratios c/a and b/a. A small value of c/a suggests two possibilities: if b/a is large, it indicates an oblate distribution resembling a thin plane; if b/a is similarly small (c/a ≈ b/a), it suggests a narrow, prolate distribution resembling a filament-like shape (Metz et al. 2007).

2.3. Adopted potential

Understanding the gravitational potential and mass distribution of the Milky Way is crucial to numerically integrating the orbits of test satellites. Despite extensive research over the past century, the exact distribution of the Milky Way’s mass remains a topic of debate. Current estimates of the total mass of the Milky Way range between 0.5 − 3 × 10¹² M_⊙ (Wang et al. 2020), varying based on the methodology and the underlying assumptions (Bland-Hawthorn & Gerhard 2016). Therefore, significant uncertainty surrounds the mass of the Milky Way, necessitating consideration of these ambiguities in mass distribution.

For this study, we used one fiducial potential for forward integration and three distinct gravitational potentials, as detailed in Table 1, to investigate their effects on the evolution and inferred stability of the plane of satellite galaxies in backward integration. This allowed us to evaluate the effect of adopting an incorrect Milky Way mass on the inferred satellite plane’s evolution. To numerically integrate orbits, we relied on Galpy (Bovy 2015) to investigate the orbital properties of the test satellites. Within the galpy.potential module, the MWPotential2014 model represents the Milky Way’s gravitational potential and is obtained by fitting to a wide range of observational data. This simplified static potential model comprises the following components:

A spherically symmetric power-law density potential with an exponential cut-off for the bulge:

$$\begin{array}{c}\hfill \rho (r)={\rho }_0{r}^{-\alpha }exp(-{\left(\frac{r}{r_c}\right)}^2),\end{array}$$(7)

with a power-law index of α = 1.8 and a cut-off radius of r_c = 1.9 kpc. The disk component is modeled using the axisymmetric Miyamoto–Nagai disk (Miyamoto & Nagai 1975):

$$\begin{array}{c}\hfill \mathrm{\Phi }(R,z)=-\frac{{\mathrm{\Phi }}_0}{\sqrt{R^2+{(a+\sqrt{z^2+b^2})}^2}},\end{array}$$(8)

where a = 3 kpc and b = 0.28 kpc.

The dark matter halo is characterized by the NFW halo (Navarro et al. 1996):

$$\begin{array}{c}\hfill \rho (r)=\frac{{\rho }_0}{(r/h){(1+r/h)}^2},\end{array}$$(9)

with a characteristic radius h = 16 kpc.

The MWPotential2014 model is scaled such that the circular velocity at a distance of 8 kpc from the galactic center within the disc (z = 0 kpc) is 220 km/s, with a virial mass of approximately 8 × 10¹¹ M_⊙. In this study, we refer to this potential as the fiducial potential, MW_fiducial. For the two remaining Milky Way potentials, we adjusted MW_fiducial by rescaling the halo mass to achieve the desired mass distribution. Specifically, the MW_low model was configured with a halo mass 40% lower than that of MW_fiducial, while MW_high had a halo mass 40% higher than that of MW_fiducial, thereby providing deeper insights into the influence of mismatching the Milky Way’s mass on the inferred behavior and stability of satellite galaxy planes.

2.4. Forward Integration

We conducted forward integration of the satellite system over a 5 Gyr time interval for N = 20 random realizations of the initial setup. We selected three distinct maximum tangential velocity angles, θ_max = 40°, 60°, and 80°, for forward integration, as shown in Table 2. Notably, larger angles allow for more eccentric orbits.

The resulting cumulative distribution functions (CDFs) of eccentricity for the three different values of θ are shown in Fig. 2. We also overlaid, as a solid black curve, the eccentricity data for observed Milky Way satellites from (Li et al. 2021) onto the plot, which contains data for 46 satellite galaxies from Fritz et al. (2018), McConnachie & Venn (2020) and Simon (2019). However, for our comparison, we selected 35 satellite galaxies, as 11 had eccentricities e ≫ 1. This addition allowed us to compare the simulated eccentricities with those of observed Milky Way satellite galaxies. From the comparison, it is apparent that θ_max = 80° aligns closely with the observed eccentricities of satellite galaxies. Each forward integration was conducted using the MW_fiducial potential model. Figure 2 also shows that the CDF range for Li et al. (2021) spans approximately 0.2 to 1.0. For the simulated eccentricities of our initial test runs, the range was between 0.0 and 1.0, because circular orbits were also included. To align the simulated eccentricities with the observed data, we excluded values ranging between θ = −20° and 20° to ensure that perfectly circular orbits were not generated.

Fig. 2. Comparison of eccentricity (e) distributions across different θ. The figure presents the eccentricity CDF for three different θmax: 40°, 60°, and 80°, from left to right, respectively, shown in gray. The solid black line in each panel represents the observed data for Milky Way satellites from (Li et al. 2021), while the overlaid gray CDFs correspond to simulated results for Nsat test satellites.

2.5. Uncertainties

To precisely determine the motion of a satellite galaxy, its position, distance, and 3D motion are crucial. Spectroscopy is used to track the motion of a galaxy via the Doppler effect. However, it is essential to note that radial velocity measures only one dimension. To determine how a galaxy moves tangentially across the sky, proper motions are measured. The Hubble Space Telescope (HST) has measured the proper motion of most of the 11 brightest classical satellite galaxies of the Milky Way (Redd 2018). Since then, largely due to the advent of Gaia, the availability of proper motion measurements has improved significantly (Gaia Collaboration 2016a,b). However, highly accurate proper motion measurements remain limited and are available only for the more massive, luminous, or nearby satellitegalaxies.

According to Gaia Data Release 2 (DR2), the random proper motion errors for Milky Way satellites range from approximately 0.006 to 0.641 mas yr⁻¹ (Fritz et al. 2018), and the corresponding systematic errors are estimated to be approximately 0.035 to 0.066 mas yr⁻¹ (Gaia Collaboration 2018a,b). In Gaia Early Data Release 3 (EDR3), the uncertainties in proper motions, derived from both statistical and systematic errors, have approximately halved (Lindegren et al. 2021; McConnachie & Venn 2020; Li et al. 2021; Battaglia et al. 2022).

In this study, we investigated four different ranges of proper motion uncertainties: ϵ_μ = 0.00 mas yr⁻¹, 0.04 mas yr⁻¹, 0.08 mas yr⁻¹, and 0.12 mas yr⁻¹. We specifically opted for 0.00 mas yr⁻¹ proper motion uncertainties to represent cases where proper motion was measured with the highest accuracy. This condition allowed us to examine the evolution of the plane of satellite galaxies under various parameters, independent of proper motion errors. We randomly sampled these proper motion uncertainties from Gaussian distributions withwidth ϵ_μ.

For distance uncertainties, we chose two specific values: ϵ_dist = 0% and 5%. Distances are crucial factors in understanding the spatial distribution and dynamics of satellite galaxies around the Milky Way, as well as their effects on the inferred stability of the plane of satellite galaxies. Again, distance errors were sampled from a Gaussian with width ϵ_dist.

2.6. Backward integration

After completing the forward integration, we recorded the orbital parameters, such as right ascension, declination, proper motion in right ascension and declination, distance, and line-of-sight velocity, for each test satellite system at the final snapshot of all simulation realizations. We then mocked-observed these systems by adding uncertainties in proper motions and distance such that

$$\begin{array}{cc}\hfill {\mu }_{\alpha {\ast }_{\text{new}}}& ={\mu }_{\alpha \ast }+{ϵ}_{\mu }\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill {\mu }_{{\delta }_{\text{new}}}& ={\mu }_{\delta }+{ϵ}_{\mu }\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill {\text{dist}}_{\text{new}}& =\text{dist}+{ϵ}_{\text{dist}}.\hfill \end{array}$$(12)

Equations (10)–(12) provide the updated proper motions and distance parameters for each test satellite, incorporating their respective uncertainties. The values of ϵ_μ were drawn from a Gaussian distribution. In contrast, for ϵ_dist, a percentage error was drawn from a Gaussian distribution, so that more distant satellites had larger absolute distance uncertainties. Table A.1 details the ranges of uncertainties in these parameters. After performing mock observations, we integrated the system backward for 5 Gyr. We repeated this process 30 times for each of the 20 forward integrations, using different randomly generated errors each time, resulting in a total of 600 backward integrations to return the test satellite system to its initial state. We conducted this procedure using all three Milky Way potential models. Each row in Table A.1 represents a combination of various parameters that could influence the inferred stability of the satellite plane. Overall, with four proper motion uncertainties, two distance uncertainties, and three potential models, there were 24 model combinations for each θ_max. Each model underwent 600 backward integrations, totaling 43 200 integrations in total.

3. Results

In this section we present the main results of the simulation¹. For brevity, and because it most closely resembles the eccentricities of the observed Milky Way satellite galaxies, we show only results for θ_max = 80 °. The results for the other θ_max are qualitatively similar and are provided in Appendix A.

3.1. Effect of proper motion uncertainties

Figure 3 presents results that focus on introducing proper artificial motion and distance uncertainties to the test satellite system before backward integration. These results are divided into four panels, each corresponding to a distinct level of proper motion uncertainty.

Fig. 3. Effect of proper motion uncertainties on the plane of satellite galaxies. Each panel shows Nrealization forward integrations as black curves with the mean indicated by a dot-dashed line, and Mrealization backward integrations as green curves with the mean indicated by a dashed line.

Each panel shows an array of curves, in gray and green, indicating the N_realization and M_realization numbers for forward and backward integration, respectively. The average value at each time step is shown as a dashed-dotted curve, with arrows pointing to the right and left. In this case, we adopt MW_fiducial for both forward and backward integrations, without introducing any distance uncertainties. In the upper left panel, when ϵ_μ = 0.00 mas yr⁻¹, the mean minor-to-major ratio, c/a, obtained from backward integration, aligns perfectly well with the mean minor-to-major ratio obtained from forward integration, as each forward integration curve is overplotted by backward integration. This demonstrates a high degree of predictability and confirms that the orbital plane remains stable and unchanged in scenarios where the proper motion is measured with perfect accuracy. It also implies that, at least in the theoretical framework of the simulation, the inherent dynamics of the satellite system are sufficiently deterministic to allow precise backtracking of its orbital path, provided that the initial conditions are precisely known. These results serve as a control scenario against which the other panels, which include proper motion errors, can be compared. Model BI-01 from Table A.1 shows similar behavior: even after 3 Gyr of backward integration, the difference between forward and backward, Δ_c/a = 0.0, and the fractional change, f_c/a, is 1.00, meaning that both forward and backward integration yield the same output when proper motion uncertainties are zero.

In the upper right panel, we have applied Gaia-level systematic uncertainties ϵ_μ = 0.04 mas yr⁻¹. Here, a slight uncertainty in proper motion leads to a visible divergence between the mean axial ratios of forward and backward integrations. The mean axial ratio from the backward integration does not align as closely with the forward integration as in the first case, when there were no uncertainties. This suggests a decrease in the predictability and inferred stability of the orbital plane with the introduction of proper motion errors. However, the curve appears to retain its shape to some degree. The effect of uncertainties is also evident from Model BI-02 in Table A.1, which shows that with the introduction of ϵ_μ = 0.04 mas yr⁻¹, Δ_c/a shifts from 0.00 to 0.073, and f_c/a becomes 1.69. Thus, the inferred past plane width is typically 70% higher than the true value.

Increasing the proper motion uncertainties to 0.08 mas yr⁻¹ further increases the divergence between the means of backward and forward integration, as shown in the lower left panel. This indicates a further decline in the predictability and stability of the orbital plane as proper motion uncertainty increases. Moreover, Table A.1 also shows that in Model BI-03, Δ_c/a shifts yet again, from 0.07 to 0.17, and f_c/a rises to 2.64, meaning that the typical inferred plane width is almost three times larger than the true one.

Under extreme proper motion uncertainties, 0.12 mas yr⁻¹, the divergence between the mean forward and mean backward integration, shown by dashed curves, is more pronounced. The mean axial ratio from backward integration shows a significant divergence from the forward integration mean, implying a severe reduction in the orbital plane’s predictability and stability with the highest proper motion uncertainty examined.

Moreover, the original shape of the curve is no longer discernible and deviates significantly more from the forward integration curve compared to the other uncertainties examined. This suggests that the actual past evolution of the system cannot be reliably traced under such high uncertainties. This shows that the inferred dynamics of the system are highly sensitive to proper motion uncertainties, making it challenging to draw reliable conclusions about the satellites’ orbital evolution under these conditions. This case is represented by Model BI-04 in Table A.1, which shows that Δ_c/a = 0.24 and f_c/a = 3.32.

3.2. Effect of intrinsic plane height

We also tested whether the effect of proper motion uncertainties on the plane of satellite galaxies was sensitive to our choice of the initial vertical distribution of the intrinsic satellite planes. Specifically, we ran the system with three different plane heights – h = 10 kpc, h = 20 kpc, and h = 30 kpc – bracketing the range of observationally inferred satellite plane heights for the Milky Way and M31 (Ibata et al. 2013; Pawlowski et al. 2013).

Figure 4 displays two panels showing the results for adopted proper motion uncertainties of 0.04 mas yr⁻¹ and 0.08 mas yr⁻¹. In both plots, the dash-dotted curves represent mean axis ratios of the forward integrations for different heights of h = 10, 20, and 30 kpc, indicated by orange, gray and cyan, respectively. As expected, when the vertical distribution is lower, the mean axis ratio of the system is lower, making the system flatter compared to the other values of h. In backward integration, after adding proper motion uncertainties, we observe that all three curves exhibit a similar trend: errors in proper motion lead to an apparent widening of the satellite plane. This effect is substantial in all three cases, indicating that our findings are insensitive to the exact choice of intrinsic plane height.

Fig. 4. Effect of proper motion uncertainties on satellite galaxies under different initial vertical distributions of h = 10, h = 20 kpc, and h = 30 kpc. The left panel shows proper motion uncertainties of 0.04 mas yr−1, while the right panel displays uncertainties of 0.08 mas yr−1. Dashed curves indicate forward integrations and solid curves represent backward integrations.

Table 3 summarizes the quantitative impact of proper motion uncertainties on the inferred change in the mean axis ratio. For a lower intrinsic plane height of h = 10  kpc (model BI-L10 and BI-H10 in Table 3), the relative effect is stronger, with f_c/a = 2.45 (3.94) compared to f_c/a = 1.64 (2.25) for the model with h = 20 kpc, for adopted proper motion uncertainties of 0.04 (0.08) mas yr⁻¹. Similarly, for an intrinsically wider plane (models BI-L30 and BI-H30 in Table 3), the relative change in plane flattening is slightly lower, with f_c/a = 1.29 and 1.72 for proper motion uncertainties of 0.04 and 0.08 mas yr⁻¹, respectively. However, the absolute change in axis ratio, measured via Δ_c/a, does not differ as significantly.

3.3. Effect of distance uncertainties

To further assess the impact of observational uncertainties, we introduced additional 5% distance uncertainties prior to backward integration. This allowed us to study the behavior of the plane under both proper motion and distance uncertainties. In all four panels of Fig. 5, backward integration includes two curves: one without distance uncertainties (shown in green) and one with 5% distance uncertainties (shown in red). As shown in the upper left panel of Fig. 5, when proper motion uncertainties are ignored, distance uncertainties have a noticeable effect: the system does not return to its original state but instead experiences a slight change. Despite this, the mean axis ratios from the backward and forward integrations remain in good agreement, suggesting that the inferred fundamental orbital evolution is not significantly affected, even with the introduction of 5% uncertainties in the distances of the test satellites. This is also evident in Table A.1, where in models BI-01 and BI-13, both Δ_c/a and f_c/a change only slightly after 3 Gyr of backward integration, from 0.00 to 0.02 and from 1.00 to 1.17, respectively. Thus, a 5% error in the distances translates to a widening of the inferred past plane height over the true value, reaching up to 17% at specific times (e.g., after 3 Gyr), although the deviation remains negligible for the majority of the 5 Gyr integration period.

Fig. 5. Similar to Fig. 3, but including a comparison with 5% distance uncertainties in the system. The green curve in each panel represents 0% distance uncertainties, while the red curve shows 5%. Proper motion uncertainties are: 0.00 mas yr−1 (upper left), 0.04 mas yr−1 (upper right), 0.08 mas yr−1 (lower left), and 0.12 mas yr−1 (lower right). Both forward and backward integrations are performed using the MWfiducial potential model.

As we introduce proper motion uncertainties, we observe that the relative effect of distance uncertainties fades away as the proper motion error becomes dominant. This is evident in the lower right panel of Fig. 5, where the highest proper motion uncertainties result in both curves showing similar behavior. Moreover, although not included in this paper, we also performed a test for 10% distance uncertainties. We find that the effect on the system becomes more pronounced with increasing uncertainties. However, this effect diminishes as the uncertainties in proper motions increase.

We also performed additional tests to examine the correlation between proper motion errors and distance errors. We find that, although there is a noticeable change, the overall effect is not extreme. This is discussed in Appendix B.1.

3.4. Effect of adopted potential

In the upper left panel of Fig. 6, without the uncertainties in proper motion and distance, the results provide a clear picture of the fundamental dynamical behavior of the satellite plane when subjected to different potential models. The colored curves depict the mean backward integration under different Milky Way potentials. The backward integration, particularly the green curve representing the MW_fiducial potential, demonstrates perfect consistency with its forward counterpart, implying a stable evolution when the potential model is the same. This case has been discussed in the previous sections.

Fig. 6. Evolution of the plane of satellite galaxies with proper motion uncertainties included in each panel. The black dashed-dotted curves represent the mean axis ratio for forward integration, while the colored dashed curves show the mean axis ratio for backward integration under different potential models characterized by halo mass. Specifically, the blue, green, and red dashed curves correspond to the c/a values under the MWlow, MWfiducial, and MWhigh potential models, respectively. Note that these results do not account for distance uncertainties.

The red curve, representing the MW_high model, exhibits a reasonable degree of stability, although it does not completely resemble the forward integration, suggesting some sensitivity to the host galaxy’s mass assumptions. This difference is shown in BI-09 of Table A.1, which shows Δ_c/a = 0.01 after 3 Gyr, and f_c/a = 1.13. In contrast, the blue curve, representing the MW_low mass model, shows a markedly different behavior. Its axis ratio initially increases before flattening, which could imply a dynamical history distinct from the correct one for the plane of satellites within such a lower mass potential model. Here, for the lower mass Milky Way potential, Δ_c/a = 0.03 after 3 Gyr, and f_c/a = 1.27.

The other three panels of Fig. 6 show the effect of both potential models and proper motion uncertainties on the stability of the plane of satellites. In the upper right panel, we have proper motion uncertainties of ±0.04 mas yr⁻¹. The effect under MW_fiducial was already discussed in the previous section. For the MW_high model, the trend is similar to MW_fiducial, but with slightly larger axis ratio values. The MW_low model again shows a distinct behavior, where it initially rises and then flattens out over time.

When adding more uncertainties, up to Gaia levels and beyond, the effect strengthens in all the remaining panels, and the plane of the system becomes even wider. Under 0.12 mas yr⁻¹ proper motion uncertainty, all three systems show a similar trend, suggesting that higher proper motion uncertainties overshadow the effect of using an incorrect Milky Way model. This observation indicates that, at higher levels of proper motion uncertainty, the distinctions between the various Milky Way potential models become less pronounced, as the increased uncertainty tends to dominate the inferred dynamics of the satellite plane.

Building upon this, we introduced an additional layer of complexity by incorporating distance uncertainties. A comparison between Fig. 6 and Fig. 7, reveals a noticeable effect of distance uncertainties across all three potential models, but the evolution of the plane does not change drastically. However, in the presence of higher proper motion uncertainties, the effect of distance uncertainties is again overshadowed under all three potential models. Models BI-13 to BI-24 show this behavior when 5% uncertainties are induced.

Fig. 7. Similar to Fig. 6, but here we include a 5% uncertainty in distances.

We also verified that our results remain qualitatively unchanged when using different Milky Way gravitational potentials. Specifically, we performed the analysis with MilkyWayPotential2022 from Gala (Price-Whelan 2017), including variations in halo mass and proper motion uncertainties, finding no substantial differences in the overall trends (see Appendix B.2).

3.5. Analysis of individual axes

The behavior of the minor-to-major axial ratio in the MW_low potential model differs from that in other models, warranting a closer examination of how the extent of the test satellite system changes along its axes, from major to minor. Figure 8 illustrates the influence of proper motion uncertainties on the mean values of the three axes: a (major), b (intermediate), and c (minor). The influence of proper motion uncertainties is organized by row, while the axes are displayed by column in the figure. Each panel includes curves in three colors, depicting the impact of different potential models. The solid curves represent the results when ϵ_dist = 0%, while the dashed-dotted curves show the results when ϵ_dist = 5%. The first row shows ±0.00 mas yr⁻¹ proper motion uncertainties and 0% distance uncertainties. Here, both MW_fiducial and MW_high show similar behavior across the axes, both with and without distance uncertainties. Over time, the axes tend to decrease, particularly under the MW_high potential model. In contrast, MW_low shows a completely different and opposite behavior. During backward integration, all the axes under MW_low increase over time.

Fig. 8. Comparative analysis of galactic orbital axes under proper motion uncertainties and different potential models at θtan = 80°. Panels shows the mean values of 600 major, minor, and intermediate axes calculated through backward integration. Each panel corresponds to a unique combination of proper motion uncertainty and potential models. Dash-dotted curves show the effect of distance uncertainties, while solid curves show results without distance uncertainties. Rows correspond to different levels of proper motion uncertainty, while columns represent the three axes: major (left), intermediate (center), and minor (right).

In these models, the variations in axis lengths are a direct consequence of changes in gravitational pull due to mass differences. In the MW_high model, increased halo mass leads to more compact axes as the structure becomes more gravitationally bound. Conversely, in the MW_low model, a reduction in mass leads to an elongation of the axes, particularly in a and b. This suggests that in the MW_low model, more satellites become unbound, causing them to drift away from their host galaxies. We emphasize that the MW_fiducial model serves as a baseline for the MW_high and MW_low models, with the halo mass adjusted by 40% (increased and decreased, respectively). However, the response of the axis extents is striking, with the Milky Way potential with a 40% reduction in halo mass having a more pronounced effect on the system than a 40% increase. Moreover, including distance errors causes each axis in every potential model to extend slightly, indicating a noticeable, though moderate, effect on the orbital axes.

We also tested and recalculated the c/a ratio by excluding test satellites beyond 300 kpc, considering only those within this threshold, as discussed in Appendix B.3. We find that this does exclusion not result in a major influence.

3.6. Inferred escaping satellites

Figure 9 shows the average number of test satellites that surpass a distance of 300 kpc from the galactic center under proper motion uncertainties and various potential models. The plot is divided into four panels that show three-color curves representing the impact of potential models. The solid curves represent the results for ϵ_dist = 0%, while the dashed-dotted curves correspond to ϵ_dist = 5%. In the upper left panel, when there are no proper motion uncertainties, the MW_low curve initially shows a relatively stable number of satellites exceeding 300 kpc in backward integration, both with and without distance uncertainties. Over time, the number gradually increases before reaching a plateau. This outcome is expected, as reducing the halo mass by 40% lowers the gravitational binding force, allowing more satellites to become unbound and drift beyond 300 kpc.

Fig. 9. Average number of test satellite exceeding a radial distance of 300 kpc under varying proper motion uncertainties. The four panels correspond to different level of proper motion uncertainty. Each panel displays six curves, representing the number of test satellites beyond 300 kpc of distance under three different potential models. Dash-dotted curves show results with distance uncertainties included, while solid curves reflect only proper motion errors.

Unlike the MW_low model, the other two potential models exhibit opposite behavior. In MW_fiducial, the number of test satellites that go beyond 300 kpc decreases over time, and at 5 Gyr, all N_sat test satellites remain under the 300 kpc limit. This trend is even more pronounced in the MW_high potential model, where the number of distant test satellites progressively decreases over time due to the higher halo mass. Furthermore, the other three panels demonstrate the effect of proper motion uncertainties on satellite dynamics, showing that as uncertainties increase, more test satellites remain at distances greater than 300 kpc across all models. Proper motion errors preferentially increase the inferred velocities, and thus satellite energies, effectively ejecting more test satellites.

Furthermore, when comparing the curves with and without distance uncertainties, we observe a clear effect: these uncertainties increase the number of test satellites exceeding 300 kpc. In Table A.1, the column N_sat > 300 kpc lists models BI-05 to BI-08 without induced distance errors in MW_low, and models BI-17 to BI-20, with distance errors. These results indicate that the inclusion of distance errors at 3 Gyr leads to a greater inferred loss of satellites in the system.

4. Conclusion and outlook

We investigated the impact of different types of measurement uncertainties on the inferred stability of artificial, intrinsically stable planes of satellite galaxies. Our key results are as follows.

• When the test satellite galaxies’ proper motions are measured with perfect accuracy (or ϵ_μ = 0.00 mas yr⁻¹) and no other errors are introduced, we infer that the system remains highly stable.

• At the level of minimum Gaia uncertainties (with ϵ_μ = 0.04 mas yr⁻¹, similar to Gaia systematics), the inferred plane widths begin to diverge from the true value toward larger values. Although the thin initial plane remains relatively well constrained, its inferred width relative to the true plane almost doubles within 3 Gyr. As we increase the applied proper motion uncertainties (ϵ_μ = 0.08 and 0.12 mas yr⁻¹), the inferred plane widths increase substantially, and one might infer that the intrinsically stable satellite planes are highly unstable. The c/a shows a linear increase with the addition of more uncertainties.

• When performing backward integration with three different intrinsic plane heights (h = 10, 20, and 30 kpc), we find that applying realistic uncertainties in the proper motion introduces substantial broadening in the plane of the satellites in all cases. The apparent broadening effect we identified is therefore not sensitive to the exact intrinsic plane height assumed for the plane.

• When applying 5% distance uncertainties, we observe a noticeable effect that increases the plane width during backward integration relative to the intrinsic plane width. However, this effect does not dominate the inferred evolution of the satellite plane. Additionally, the impact of distance uncertainties becomes less important with higher levels of proper motion uncertainties.

• The adoption of incorrect Milky Way potential models leads to an apparent widening of the inferred past satellite planes. This impact is especially significant for MW_low, which has a 40% smaller halo mass compared to the fiducial model considered for our intrinsic satellite planes. Due to its lower mass, MW_low struggles to maintain cohesion among the test satellites, leading to increased radial distances and an extension of its major, minor, and intermediate axes. However, under larger proper motion uncertainties, the influence of adopting incorrect Milky Way potential models diminishes, as the system deviates from its true configuration across all potential models.

Thus, this study emphasizes the necessity of precise determination of proper motion, distances of satellite galaxies, and the Milky Way’s gravitational potential. Variations in these quantities can significantly shift our inference of the evolution of intrinsically stable satellite planes, especially with respect to proper motion and halo mass of the host galaxy. We argue that underestimating the halo mass also leads to the inference that more satellites escape from the host (and might thus appear unbound) than in the true fiducial model.

Accurate measurement of the Milky Way potential is therefore vital for a comprehensive study of its satellite system and for broader insights into galactic evolution and interactions. We also find that the results are not very sensitive to the number of integrations; in our study, 600 backward integrations per model are sufficient.

However, it is crucial to acknowledge certain limitations in this research. A significant caveat is the omission of the gravitational effects of massive satellite galaxies such as the LMC. These celestial bodies, owing to their substantial mass, can affect the mass estimates of the Milky Way galaxy (Erkal et al. 2020) and exert significant gravitational forces that profoundly influence the Milky Way and its surroundings (Garavito-Camargo et al. 2021; Pawlowski et al. 2022; Vasiliev 2023), and thus the stability of the satellite plane. This omission represents a notable limitation in our analysis. Additionally, this study assumes a spherical dark matter halo; however, it is important to acknowledge that the halo could be triaxial (Law et al. 2009), which might have implications for our findings.

Future studies could enhance our understanding by incorporating the gravitational influence of massive satellite galaxies, such as the LMC and SMC, as well as by considering different shapes of the dark matter halo. Such an approach would allow for a more comprehensive analysis of the satellite plane dynamics, taking into account the complexities of real-world scenarios. Additionally, in this study, we assume a static dark matter halo over the 5 Gyr integration period. However, this simplification overlooks important evolutionary effects. Over time, the halo mass is expected to grow as it accretes dark matter, resulting in a smaller mass in the past and a larger mass in the future. This affects satellite orbits, particularly at larger radii. Pace et al. (2022) focuses on observations and finds that the orbital computations of 46 dwarf spheroidal galaxies (dSphs) using a static Milky Way potential reveal that, while pericenter and apocenter distances can be derived, the inclusion of the LMC’s potential alters these parameters by more than 25% for 40% of the sample. Similar studies have also been conducted using cosmological simulations D’Souza & Bell (2022), Santistevan et al. (2024). Additionally, the influence of the baryonic disk leads to adiabatic contraction, increasing the dark matter density at smaller radii, which affects the inner orbits of satellites (Cautun et al. 2020). However, in general, such additional dynamical effects are expected to increase the uncertainty in past satellite trajectories, as we have only limited constraints on their exact implementation. We therefore expect them to exacerbate the trend of inferring less stable satellite plane dynamics when backward integrating. Moreover, our models does not include an important dynamical process: dynamical friction. This effect, especially for massive satellites orbiting within a host potential, can lead to the gradual decay of orbits as they lose energy and angular momentum to the dark matter background (Esquivel & Fuchs 2007; Ogiya & Burkert 2016; Taylor & Babul 2001). Over time, this process can significantly alter satellite trajectories and contribute to the disruption of coherent planar structures.

Our method can be extended to examine the observed planes of satellites of M31, Cen A, and others. However, the application to these systems requires careful consideration of observational uncertainties, such as proper motions, distances, and gravitational potentials, which are significantly greater compared to our current analysis (in the case of M31) or thus far inaccessible (for more distant systems). In addition to these observational limitations, the evolution of such planes is also sensitive to the initial distribution of the positions and velocities of the test satellites. Variations in the initial conditions could lead to notable differences in the structure and persistence of the planes over time. By accounting for these factors, our approach can provide valuable insights into the dynamics of satellite galaxies in diverse environments, further broadening the scope of our findings.

Our numerical experiments demonstrate how intrinsically stable, co-orbiting planes of satellite galaxies display a significant increase in their plane height when backward integrated after applying measurement errors. This effect is already evident for proper motion errors comparable to the current Gaia systematics of ∼0.033, mas yr⁻¹ (Lindegren et al. 2021; McConnachie & Venn 2020; Li et al. 2021). Sawala et al. (2023) argued that such an increase in plane height indicates that the Milky Way’s Vast Polar Structure is more likely a transient rather than a rotationally supported feature, and thus consistent with ΛCDM expectations. Similar arguments were previously put forward by Maji et al. (2017), based on less accurate proper motion measurements. In their studies, the authors integrated the orbits of observed Milky Way satellite galaxies in a static potential, akin to our numerical experiments, and measured the time evolution of the inferred thickness of plane fits. However, our work shows that such an increasing thickness does not necessarily imply an absence of a kinematically coherent, dynamically stable, or long-lived satellite plane. In the case of an underlying, intrinsically stable satellite plane, realistic measurement errors alone, as well as assuming a mismatching Galactic potential, can already cause apparent thickening. Furthermore, care must be taken in how errors are applied. For example, the measured most-likely proper motions of observed satellite galaxies are already subject to measurement error. Sampling the errors in a Monte Carlo fashion (e.g., as demonstrated by Sawala et al. 2023) then effectively applies these errors twice. In our study, this means that our applied proper motion error of ϵ_μ = 0.04 mas yr⁻¹ corresponds to a Monte Carlo sampling of orbits based on observed data with a measurement error of ${ϵ}_{\mu }/\sqrt{2}=0.028{\mathrm{mas}\mathrm{yr}}^{-1}$. Monte Carlo sampling therefore enhances the effect we identified in our study of artificially increasing the inferred satellite plane height when backward integrating, because it typically adds more dispersion to an ensemble of satellites. Merely finding an increase in the inferred plane height of an observed satellite structure under orbit integration–without accounting for the impact of measurement errors (especially in proper motions) or mismatches in the assumed Galactic gravitational potential–is therefore insufficient to claim that the structure is not dynamically stable or transient, or that it is therefore consistent with ΛCDM.

Acknowledgments

MSP acknowledges funding via a Leibniz-Junior Research Group (project number J94/2020). Moreover, this work made use of the following software packages: astropy (Astropy Collaboration 2013, 2018, 2022), matplotlib (Hunter 2007), numpy (Harris et al. 2020), python (Van Rossum & Drake 1995), scipy (Virtanen et al. 2020; Gommers et al. 2025), and galpy (Bovy 2015). Software citation information aggregated using https://www.tomwagg.com/software-citation-station/ (Wagg & Broekgaarden 2024; Wagg et al. 2024).

References

Appendix A: Analysis at θ_max = 40° and 60°

This appendix presents the simulations involving less eccentric orbits. Results and analysis for θ_tan = 80% have been discussed in detail in Section 3. Following a similar pattern, we here present the results which are obtained for θ_tan = 40°and 60°.

Figure A1 and A2 show results analogous to those in Section 3. These scenarios exhibit qualitatively comparable patterns. The figures are segmented into four panels, each representing different levels of proper motion uncertainties. The mean c/a ratios for forward integration are depicted with black curves, whereas those for backward integration are shown in colored curves. Similar to the previously discussed case, it becomes increasingly difficult to accurately determine the stability of the satellite plane as the level of proper motion errors escalates. Furthermore, in these results too, the behavior of the system under MW_low is different from the other two potential models.

Figure 2 shows that as θ_max increases, the orbits become more eccentric, and this effect is clearly observed when we consider the results under θ_max. By comparing the results at θ_max = 40%, 60%, 80%, as shown in Tables A2, A3, and A.1, we see that at a given level of proper motion uncertainty, the satellite plane widens more for larger θ_max. For instance, in Model BI-02, when ϵ_μ = 0.04 mas yr⁻¹, Δ_c/a is 0.04, 0.06, and 0.07 for θ_tan = 40%, 60%, 80%, respectively, and increases gradually. Furthermore, the influence of eccentric orbits is also apparent when comparing the curves in Fig. A1, A2, and 6. In cases with less eccentric orbits (θ_tan = 40%), the plane appears relatively smoother and flatter. In contrast, as the orbits become more eccentric, the evolution of the plane flattening becomes increasingly irregular and pronounced, indicating a significant impact of orbital eccentricity on the satellite’s dynamics. This suggests that as θ_max increases, not only does the eccentricity of the orbits intensify, but the overall structure and stability of the orbital plane are also affected, leading to a more complex and varied orbital configuration.

Similarly, Fig. A3 and A4 display the effects of distance uncertainties on the test satellite system for θ_tan= 40°and 60°, respectively. The analysis, as discussed in Section 3, indicates that the inclusion of distance errors has a noticeable effect on the c/a. However, this does not significantly impact the overall evolution of c/a.

Finally, similar to Fig. 9, here Fig. A5 and A6 display the number of test satellites that exceeded the 300 kpc radius threshold under θ_max = 40°and 60°, respectively. The effect of proper motion errors is depicted in each of the four panels, with each line within a panel representing the influence of different Milky Way potential models. Dash-dotted curves show the effect of 5% distance uncertainties, while solid curves show 0% distance uncertainties. The influence of eccentricities is clearly apparent in these figures. For a given model, less eccentric orbits result

Fig. A.1. Analogous to Fig 6, the results for θtan= 40°: Evolution of the plane of satellite galaxies with proper motion uncertainties included in each panel. The black dashed-dotted curves represent the mean axis ratio for forward integration, while the colored dashed curves indicate the mean axis ratio for backward integration under different potential models characterized by halo mass. Specifically, the blue, green, and red dashed curves correspond to the c/a values under the MWlow, MWfiducial, and MWhigh potential models, respectively. Note that these results do not account for distance uncertainties.

Fig. A.2. Analogous to Fig 6, results for θtan= 60°: Evolution of the plane of satellite galaxies with proper motion uncertainties included in each panel. The black dashed-dotted curves represent the mean axis ratio for forward integration, while the colored dashed curves indicate the mean axis ratio for backward integration under different potential models characterized by halo mass. Specifically, the blue, green, and red dashed curves correspond to the c/a values under the MWlow, MWfiducial, and MWhigh potential models, respectively. Note that these results do not account for distance uncertainties.

Fig. A.3. Analogous to Fig 7, results for θtan= 40°: Evolution of the plane of satellite galaxies with proper motion uncertainties included in each panel. The black dashed-dotted curves represent the mean axis ratio for forward integration, while the colored dashed curves show the mean axis ratio for backward integration under different potential models characterized by halo mass. Specifically, the blue, green, and red dashed curves correspond to the c/a values under the MWlow, MWfiducial, and MWhigh potential models, respectively. Note that in these results we account 5% of distance uncertainties.

Fig. A.4. Analogous to Fig 7, results for θtan= 60°: Evolution of the plane of satellite galaxies with proper motion uncertainties included in each panel. The black dashed-dotted curves represent the mean axis ratio for forward integration, while the colored dashed curves indicate the mean axis ratio for backward integration under different potential models characterized by halo mass. Specifically, the blue, green, and red dashed curves correspond to the c/a values under the MWlow, MWfiducial, and MWhigh potential models, respectively. Note that in these results we account 5% of distance uncertainties.

Fig. A.5. Average number of test satellite exceeding 300 kpc of radial distance under proper motion uncertainties. Similar to Fig. 9, but under θtan= 40°.

Fig. A.6. Average number of test satellite exceeding 300 kpc of radial distance under proper motion uncertainties. Similar to Fig. 9, but under θtan= 60°.

in fewer test satellites escaping the system compared to those with higher eccentricities. This trend is further corroborated by the data presented in Tables A.1, A2, and A3. In these tables, under all models, the column representing the number of satellites with distances greater than 300 kpc (N_sat > 300 kpc) consistently increases as the eccentricity of the orbits increases.

Appendix B: Additional checks

Apart from the main results, we also checked some additional checks to see the behavior of the satellite plane.

B.1. Correlated PM and Dist errors

We also test whether a correlation between proper motion and distance errors affects the inferred evolution of the satellite plane, as fainter test satellites could have both larger proper motion and distance uncertainties. To investigate this, we performed an additional analysis, comparing two cases (see Fig. B1): one with a proper motion error of 0.04 mas yr⁻¹ and a 10% distance error, and another with a PM error of 0.1 mas yr⁻¹ and a 10% distance error. Each plot includes two backward integrations–one with a PM-distance error correlation and one without. In both cases, the same errors are used, but in the ’correlated’ case, the errors are ranked and assigned to the test satellites such that satellites farther away have progressively larger distance and proper motion errors. We find that the effect is not extreme but leads to an increase in the inferred widening of the satellite plane. This suggests that the approach of using uncorrelated errors is conservative and at most underestimates the impact of uncertainties. Moreover, as expected, with higher PM error, the effect becomes more pronounced, aligning with what we observed in our main results.

B.2. Backward integration with MilkyWayPotential2022

In Fig. B2, we re-ran some of our simulations and show a comparison between the backward integration done by MilkyWayPotential2022 from Gala. We found that our results remain qualitatively unchanged. In Fig. B2, the first column corresponds to MWPotential2014, and the second to MilkyWayPotential2022. In each plot, we include one example with no proper motion error and one with a 0.04 mas/yr proper motion error. While some minor differences are noticeable, the overall behavior and results of our study remain unaffected by the choice of potential. We also tested a +40% halo mass variation for both MWPotential2014 and MilkyWayPotential2022, and the results remained consistent.

B.3. Re-computing c/a but with test satellites within 300 kpc

Figure B3 shows the evolution of c/a for the MW_fiducial potential model, comparing cases with and without the exclusion of test satellites beyond 300 kpc. The dashed green line represents the backward evolution when these satellites are excluded, while the solid green line shows the backward evolution including all test particles. The black lines indicate the forward integration. The exclusion of satellites beyond 300 kpc does not significantly alter the overall evolution of c/a. While removing some test satellites has a minor effect on the backward integration of c/a, this effect is not substantial if satellites beyond 300 kpc are excluded.

Fig. B.1. Comparison of correlated and uncorrelated PM-distance errors on satellite plane evolution. Two cases are shown with PM errors of 0.04 and 0.1 mas yr−1, both with a 10% distance error. Correlated errors lead to a slightly increased plane widening.

Fig. B.2. Comparison between MilkyWayPotential2022 and MWPotential2014

Fig. B.3. Evolution of c/a for the MWfiducial model with and without excluding test satellites beyond 300 kpc. Dashed and solid green lines show backward evolution, black lines indicate forward integration. Exclusion has minimal impact on the overall trend.

All Tables

All Figures

	Fig. 1. Face-on (left) and edge-on (right) views of Nsat = 25 randomly generated test satellites, each represented by a different color. The star symbol denotes the host galaxy.
In the text

Fig. 2. Comparison of eccentricity (e) distributions across different θ. The figure presents the eccentricity CDF for three different θmax: 40°, 60°, and 80°, from left to right, respectively, shown in gray. The solid black line in each panel represents the observed data for Milky Way satellites from (Li et al. 2021), while the overlaid gray CDFs correspond to simulated results for Nsat test satellites.

In the text

	Fig. 3. Effect of proper motion uncertainties on the plane of satellite galaxies. Each panel shows Nrealization forward integrations as black curves with the mean indicated by a dot-dashed line, and Mrealization backward integrations as green curves with the mean indicated by a dashed line.
In the text

	Fig. 4. Effect of proper motion uncertainties on satellite galaxies under different initial vertical distributions of h = 10, h = 20 kpc, and h = 30 kpc. The left panel shows proper motion uncertainties of 0.04 mas yr−1, while the right panel displays uncertainties of 0.08 mas yr−1. Dashed curves indicate forward integrations and solid curves represent backward integrations.
In the text

Fig. 5. Similar to Fig. 3, but including a comparison with 5% distance uncertainties in the system. The green curve in each panel represents 0% distance uncertainties, while the red curve shows 5%. Proper motion uncertainties are: 0.00 mas yr−1 (upper left), 0.04 mas yr−1 (upper right), 0.08 mas yr−1 (lower left), and 0.12 mas yr−1 (lower right). Both forward and backward integrations are performed using the MWfiducial potential model.

In the text

Fig. 6. Evolution of the plane of satellite galaxies with proper motion uncertainties included in each panel. The black dashed-dotted curves represent the mean axis ratio for forward integration, while the colored dashed curves show the mean axis ratio for backward integration under different potential models characterized by halo mass. Specifically, the blue, green, and red dashed curves correspond to the c/a values under the MWlow, MWfiducial, and MWhigh potential models, respectively. Note that these results do not account for distance uncertainties.

In the text

	Fig. 7. Similar to Fig. 6, but here we include a 5% uncertainty in distances.
In the text

Fig. 8. Comparative analysis of galactic orbital axes under proper motion uncertainties and different potential models at θtan = 80°. Panels shows the mean values of 600 major, minor, and intermediate axes calculated through backward integration. Each panel corresponds to a unique combination of proper motion uncertainty and potential models. Dash-dotted curves show the effect of distance uncertainties, while solid curves show results without distance uncertainties. Rows correspond to different levels of proper motion uncertainty, while columns represent the three axes: major (left), intermediate (center), and minor (right).

In the text

Fig. 9. Average number of test satellite exceeding a radial distance of 300 kpc under varying proper motion uncertainties. The four panels correspond to different level of proper motion uncertainty. Each panel displays six curves, representing the number of test satellites beyond 300 kpc of distance under three different potential models. Dash-dotted curves show results with distance uncertainties included, while solid curves reflect only proper motion errors.

In the text

Fig. A.1. Analogous to Fig 6, the results for θtan= 40°: Evolution of the plane of satellite galaxies with proper motion uncertainties included in each panel. The black dashed-dotted curves represent the mean axis ratio for forward integration, while the colored dashed curves indicate the mean axis ratio for backward integration under different potential models characterized by halo mass. Specifically, the blue, green, and red dashed curves correspond to the c/a values under the MWlow, MWfiducial, and MWhigh potential models, respectively. Note that these results do not account for distance uncertainties.

In the text

Fig. A.2. Analogous to Fig 6, results for θtan= 60°: Evolution of the plane of satellite galaxies with proper motion uncertainties included in each panel. The black dashed-dotted curves represent the mean axis ratio for forward integration, while the colored dashed curves indicate the mean axis ratio for backward integration under different potential models characterized by halo mass. Specifically, the blue, green, and red dashed curves correspond to the c/a values under the MWlow, MWfiducial, and MWhigh potential models, respectively. Note that these results do not account for distance uncertainties.

In the text

Fig. A.3. Analogous to Fig 7, results for θtan= 40°: Evolution of the plane of satellite galaxies with proper motion uncertainties included in each panel. The black dashed-dotted curves represent the mean axis ratio for forward integration, while the colored dashed curves show the mean axis ratio for backward integration under different potential models characterized by halo mass. Specifically, the blue, green, and red dashed curves correspond to the c/a values under the MWlow, MWfiducial, and MWhigh potential models, respectively. Note that in these results we account 5% of distance uncertainties.

In the text

Fig. A.4. Analogous to Fig 7, results for θtan= 60°: Evolution of the plane of satellite galaxies with proper motion uncertainties included in each panel. The black dashed-dotted curves represent the mean axis ratio for forward integration, while the colored dashed curves indicate the mean axis ratio for backward integration under different potential models characterized by halo mass. Specifically, the blue, green, and red dashed curves correspond to the c/a values under the MWlow, MWfiducial, and MWhigh potential models, respectively. Note that in these results we account 5% of distance uncertainties.

In the text

	Fig. A.5. Average number of test satellite exceeding 300 kpc of radial distance under proper motion uncertainties. Similar to Fig. 9, but under θtan= 40°.
In the text

	Fig. A.6. Average number of test satellite exceeding 300 kpc of radial distance under proper motion uncertainties. Similar to Fig. 9, but under θtan= 60°.
In the text

	Fig. B.1. Comparison of correlated and uncorrelated PM-distance errors on satellite plane evolution. Two cases are shown with PM errors of 0.04 and 0.1 mas yr−1, both with a 10% distance error. Correlated errors lead to a slightly increased plane widening.
In the text

	Fig. B.2. Comparison between MilkyWayPotential2022 and MWPotential2014
In the text

	Fig. B.3. Evolution of c/a for the MWfiducial model with and without excluding test satellites beyond 300 kpc. Dashed and solid green lines show backward evolution, black lines indicate forward integration. Exclusion has minimal impact on the overall trend.
In the text