A new multiple-arc model of the resonant Kuiper belt objects – Plutinos

ABSTRACT

To incorporate the gravitational influence of Kuiper belt objects (KBOs) in planetary ephemerides, uniform-ring models are commonly employed. In this paper, for representing the KBO population residing in Neptune’s 2:3 mean motion resonance (MMR), known as the Plutinos, we introduce a three-arc model by considering their resonant characteristics. Each ‘arc’ refers to a segment of the uniform ring and comprises an appropriate number of point masses. Then the total perturbation of Plutinos is numerically measured by the change in the Sun–Neptune distance (⁠|$\Delta d_{\mathrm {SN}}$|⁠). We conduct a comprehensive investigation to take into account various azimuthal and radial distributions associated with the resonant amplitudes (A) and eccentricities (e) of Plutinos, respectively. The results show that over a 100-yr period: (1) at the smallest |$e=0.05$|⁠, the Sun–Neptune distance change |$\Delta d_{\mathrm {SN}}$| caused by Plutinos decreases significantly as A reduces. It can deviate from the value of |$\Delta d_{\mathrm {SN}}$| obtained in the ring model by approximately 100 km; (2) as e increases in the medium range of 0.1–0.2, the difference in |$\Delta d_{\mathrm {SN}}$| between the arc and ring models becomes increasingly significant; (3) at the largest |$e\gtrsim 0.25$|⁠, |$\Delta d_{\mathrm {SN}}$| can approach zero regardless of A, and the arc and ring models exhibit a substantial difference in |$\Delta d_{\mathrm {SN}}$|⁠, reaching up to 170 km. Then the applicability of our three-arc model is further verified by comparing it to the perturbations induced by observed Plutinos on the positions of both Neptune and Saturn. Moreover, the concept of the multiple-arc model, designed for Plutinos, can be easily extended to other MMRs densely populated by small bodies.

1 INTRODUCTION

The accuracy of planetary ephemerides is crucial for space exploration, and therefore, continuous development efforts are underway (e.g. Newhall, Standish & Williams 1983; Pitjeva 2001; Fienga et al. 2008, 2020; Folkner et al. 2014; Pitjeva & Pitjev 2014). Traditionally, the ephemeris computation includes mutual perturbations of the Sun, the eight planets, the Pluto system, and the satellites, along with considerations for the solar oblateness J2, the solar pressure, figure and tide effects, relativistic corrections, lunar librations, and so on (Pitjeva & Pitjev 2014). In addition to the influencing factors mentioned above, modern planetary ephemerides have been updated to incorporate the gravitational perturbations of a vast number of asteroids observed in the Solar system (Fienga et al. 2008; Folkner et al. 2014; Pitjeva & Pitjev 2014; Tian 2023). The contributions mainly come from three asteroid populations: main belt asteroids (MBAs), Jupiter Trojans (JTs), and Kuiper belt objects (KBOs).

The MBAs are distributed between the orbits of Mars and Jupiter, with a current estimated population of around 700 000 (Li, Xia & Zhou 2019). To simulate the total perturbation caused by the MBAs, the ‘Bigs + ring’ model is commonly adopted. The ‘Bigs’ refer to the MBAs that have the most significant influence on the Earth–Mars distance (e.g. 1 Ceres, 2 Pallas), and they are individually included in the numerical ephemerides calculations (Kuchynka et al. 2010). Generally, ephemerides contain about 300 ‘Bigs’, for which only the self-gravity of the three largest asteroids (i.e. 1 Ceres, 2 Pallas, and 4 Vesta) is considered (Pitjeva & Pitjev 2018a; Fienga et al. 2020) to speed up the calculation. As for the ‘ring’, which is homogeneous, it represents the remaining numerous MBAs and is simply implemented as a gravitational potential (Krasinsky et al. 2002; Fienga et al. 2008; Pitjeva & Pitjev 2014). In the ephemeris INPOP06, Fienga et al. (2008) estimated that the perturbation induced by this asteroid ring on the Earth–Mars distance can reach up to 150 m during the time interval of 1969–2010. To reduce the time cost in ephemeris calculations, Liu, Fienga & Yan (2022) proposed an alternative approach known as the multiple-ring model. They took into account all MBAs, regardless of size, and classified them into 122 families based on their proper elements (Nesvorný, Brož & Carruba 2015). Eventually, they assigned the MBAs to six rings with varying heliocentric radii and inclinations. With regard to the total perturbation of the MBAs, the six-ring model yields a Mars–Earth distance that deviates from the value obtained in the ‘Bigs + ring’ model by an error below 0.5 m over a 10-yr period. This error is acceptable, given that tracking data from Mars spacecraft are available with metric precision (Kuchynka et al. 2010). The main advantage of this new approach is that it improves computational efficiency, as there is no need to incorporate a large quantity of individual MBAs (i.e. the ‘Bigs’) in the numerical ephemerides.

JTs are asteroids that share the orbit of Jupiter, but leading and trailing Jupiter by about |$60^{\circ }$| in longitude, i.e. around the L4 and L5 triangular Lagrangian points. These asteroids are said to be settled in the 1:1 mean motion resonance (MMR) with Jupiter. As of February 2023, more than 12 000 JTs are registered in the Minor Planet Center (MPC).¹ Since the L4 and L5 swarms are distributed in two separate regions, respectively, this unique configuration of JTs cannot be simulated by a simple ring model like those used for the MBAs. Furthermore, there are many more JTs in the L4 swarm than the L5 swarm, with a number ratio expected to be between 1.3 and 2 (Jewitt, Sheppard & Porco 2004; Szabó et al. 2007; Grav et al. 2011, 2012). According to Li et al. (2023b), this number asymmetry holds for JTs with absolute magnitudes |$H\lt 15$| (equivalent to diameters |$D\gt 7$| km, assuming a typical albedo of 0.04), while the fainter objects are too few to be statistically significant. The number asymmetry of JTs also contradicts the assumption of a homogeneous-ring model. In Li & Sun (2018), besides the 226 largest JTs with |$H\lt 11$|⁠, the authors modelled the remaining objects with |$H\gt 11$| using two arcs located around Jupiter’s L4 and L5 points, respectively. Their findings demonstrated that the total effect of JTs can lead to a change of |$\sim 70$| m in the Earth–Mars distance during the 2014–2114 time interval.

The KBOs are icy celestial bodies located beyond the orbit of Neptune. The orbits of KBOs can be grouped into three classes (Gladman, Marsden & Vanlaerhoven 2008): (1) Resonant KBOs: These objects occupy the MMRs with Neptune, and they have small to large eccentricities ranging from |$e=0.05$| to 0.35. A significant portion of resonant KBOs is found in Neptune’s 2:3 MMR at |$\sim 39.4$| au, sharing this resonance with Pluto. This particular population is called ‘Plutinos’, and Alexandersen et al. (2016) estimated that there are about 9000 Plutinos with diameters |$D\gt 100$| km. Numerous resonant KBOs are also located in Neptune’s 1:2 MMR at |$\sim 47.8$| au, comprising an estimated population of about 4400 objects with |$D\gt 100$| km (Chen et al. 2019). It is worth noting that Neptune’s 2:3 MMR and 1:2 MMR positions are usually considered to be the inner and outer boundaries for the main Kuiper belt (Petit et al. 2023). In addition, the Outer Solar System Origins Survey (OSSOS) has enriched other MMRs. For example, the 2:5 resonant KBOs are more numerous than initially anticipated, and their total population may be as large as that of the Plutinos (Volk et al. 2016). (2) Classical KBOs: These objects are not in Neptune’s MMRs and typically exhibit small to moderate eccentricities and inclinations. Approximately 30 000 non-resonant KBOs with |$D\gt 100$| km are estimated to exist in the main Kuiper belt (Petit et al. 2023). (3) Scattered KBOs: These objects have highly eccentric orbits, with perihelion distances |$q\gt 30$| au and semimajor axes |$a\gt 50$| au (Jewitt, Luu & Trujillo 1998; Fernández, Gallardo & Brunini 2004). Lawler et al. (2018) proposed that there are about 90 000 scattered KBOs with |$D\gt 100$| km.

Since the late 2000s, the dynamical model of KBOs in planetary ephemerides has undergone gradual refinement. When modelling the Kuiper belt, most of the previous ephemerides mainly focused on the main Kuiper belt (Pitjeva & Pitjev 2018a, b; Fienga et al. 2022). In the ephemeris EPM2008, besides the 21 biggest KBOs, the perturbations caused by the remaining smaller objects were modelled by a 1D ring with a heliocentric radius of 43 au (Pitjeva 2010). By adjusting the mass of the ring to fit the observation data obtained from spacecraft, Pitjeva estimated that the total mass of KBOs should not exceed |$0.027M_{\oplus }$|⁠, where |$M_{\oplus }$| denotes the Earth mass. In the updated ephemeris EPM2013, Pitjeva & Pitjev (2014) improved the KBO model by individually considering the 31 most massive objects, while the smaller objects were still collectively represented by a single ring. Considering that the inner and outer edges of the main Kuiper belt correspond to Neptune’s 2:3 and 1:2 MMRs, respectively, Pitjeva & Pitjev (2018b) introduced an 8 au-wide annulus spanning from 39.4 to 47.8 au (i.e. a 2D ring) to represent the total perturbation of numerous small KBOs. However, such a 2D ring model poses a severe drawback as objects with different heliocentric distance at the annulus rotate with the same angular velocity. This is obviously different from the reality. To overcome this drawback, these authors developed a new model comprising three separate rings: two rings positioned at 39.4 and 47.8 au represent the 2:3 and 1:2 resonant KBOs, respectively, and the third ring is placed at 44 au, symbolizing the ‘core’ of the Kuiper belt predominantly inhabited by the classical KBOs. In a recent study, Di Ruscio et al. (2020) also employed the three-ring model to estimate the mass of the Kuiper belt. By fitting the high-precision measurements of Saturn obtained from the Cassini mission, they provided a total mass for the Kuiper belt of |$0.061 \pm 0.001M_{\oplus }$|⁠.

When dealing with a large number of uniformly distributed asteroids, a single- or multiple-ring model would be appropriate for representing their total perturbation. However, for the resonant KBOs, their motions are restricted to specific regions within the resonance’s phase space, rather than encompassing it entirely (Li et al. 2022), i.e. they may not distributed uniformly. It is evident that Neptune’s MMRs lead to special distributions of these KBOs’ phases with respect to Neptune, and this factor can affect the perturbation induced by the resonant KBOs. Accordingly, we aim to construct a new model to represent the total perturbation of these resonators, and estimate the level of the inaccuracy introduced when using a ring model.

In this paper, we focus on the Plutinos, which constitute the largest resonant population observed in the Kuiper belt. The resonant angle |$\sigma$| associated with Neptune’s 2:3 MMR is defined by

where |$\lambda$| and |$\varpi$| are the mean longitude and the longitude of perihelion of the Plutino, respectively, and |$\lambda _N$| is the mean longitude of Neptune. For a stable Plutino, its resonant angle |$\sigma$| librates around |$180^{\circ }$| with a resonant amplitude of |$A\lt 180^{\circ }$|⁠. Here, we designate the resonant amplitude as a half-width of the variation of |$\sigma$|⁠. Consequently, |$\sigma$| cannot traverse from 0 to |$360^{\circ }$|⁠, indicating that the mean longitudes of Plutinos may not evenly spread across the range of 0–|$360^{\circ }$|⁠. The primary objective of this study is to develop a dynamical model that account for this inhomogeneity in the azimuthal distribution of Plutinos.

In our previous study on the perturbation of JTs (Li & Sun 2018), a similar issue regarding the azimuthal distribution was discussed. In that study, JTs are allowed to be on circular orbits with eccentricities |$e=0$| because the resonant term of the 1:1 Jovian MMR in the expansion of the disturbing function does not contain e (see chapter 6 of Murray & Dermott 2000). However, Neptune’s 2:3 resonance is the eccentricity-type, and its strength is proportional to e. Consequently, Plutinos must possess e values greater than 0. When involving the contribution of e, their orbital configuration becomes more complex. Thus it is necessary to additionally consider the resulting asymmetry in the radial distribution of Plutinos.

The rest of this paper is organized as follows. In Section 2, we design the dynamical model of the Solar system and choose the change in the Sun–Neptune distance as a measurement for the perturbation of the KBOs, including the Plutinos. In Section 3, we construct both the ring model and arc model to simulate the total perturbation of the Plutinos, followed by a comparison of the results obtained from these two models. We conduct a detailed investigation of the plausible azimuthal and radial distributions of the arcs, which are determined by Plutinos’ resonant amplitudes A and eccentricities e, respectively. We then provide a concise analytical expression to describe the contributions of A, e, and the total mass of Plutinos to the change in the Sun–Neptune distance. In Section 4, we verify the applicability of our three-arc model by comparing it to the perturbations induced by observed Plutinos with debiased orbits on the positions of both Neptune and Saturn. The conclusions and discussion are given in Section 5.

2 DYNAMICAL MODEL OF THE SOLAR SYSTEM

Before studying the Plutinos, this section initiates the construction of our dynamical model, focusing on the entire KBO population. The unperturbed model of the Solar system comprises the Sun and eight planets from Mercury to Neptune. The planets’ masses, initial heliocentric positions, and velocities are taken from DE405 (Standish 1998). First, their orbits are adjusted from the mean equatorial system to the J2000.0 heliocentric ecliptic system at epoch 2021 July 5. From this epoch forward, we construct the perturbed model of the Solar system by incorporating the gravitational perturbations of massive KBOs. In the subsequent analysis, our goal is to assess the influence of KBOs on the modern planetary ephemerides. To achieve this, we compare the motion of Neptune, the closest planet to the KBOs, in both the unperturbed and perturbed models.

To quantify the gravitational effect of KBOs on the orbit of Neptune, we examine the change in the distance between the Sun and Neptune, defined as

where |$d_{\mathrm {SN}1}$| and |$d_{\mathrm {SN}0}$| indicate the Sun–Neptune distances calculated with and without perturbations from the KBOs, respectively. We note that although the planetary ephemerides should be parametrized in the barycentric coordinate system, it is common to assess the validity of a perturbation model using the relative distance between two celestial bodies (e.g. Sun–planet, planet–planet), rather than the barycentric distance (Kuchynka et al. 2010; Pitjeva 2013; Li & Sun 2018; Pitjeva & Pitjev 2018b; Liu et al. 2022). This approach is preferred due to the variation of the solar system barycentre (SSB) in different perturbation models. For instance, the inclusion of all the asteroids in the main belt and Kuiper belt can cause a displacement of the SSB of the order of 100 km (Li et al. 2019). Consequently, the barycentric positions of the Sun and all other celestial bodies would change significantly, but the relative coordinates remain the same.

In our numerical simulations of the Solar system’s evolution, we employ the 19th-order Cowell prediction-correction algorithm (PECE) with a time-step of 0.5 d, which is chosen based on the orbital period of the innermost celestial body (i.e. Mercury) in our models (Huang & Zhou 1993; Li & Sun 2018); and the data output step is set to 5 d. In this N-body code, we take into account gravitational interactions among the Sun, planets, and KBOs but ignore the effects of KBO self-gravity.

First of all, we begin by estimating the magnitude of the perturbation caused by the KBO population on the Sun–Neptune distance. For simplicity, we here consider the influence of the 11 most massive KBOs that have their individual satellites, as listed in Table 1. From the motion of their satellites, their masses are well determined by applying Kepler’s third law. The orbital elements of these objects are gathered from the MPC at the specific epoch of 2021 July 5, as mentioned earlier. This way, we can measure the perturbations from these 11 prominent KBOs through the change in the Sun–Neptune distance, denoted as |$\Delta d_{\mathrm {SN}}$|⁠. Fig. 1(a) shows that, between the years 2021 and 2121, the resultant |$\Delta d_{\mathrm {SN}}$| could attain an approximate value of 16.5 km. This substantial value indicates a significant perturbation caused by the KBOs on the Sun–Neptune distance. As a result, modern planetary ephemerides may have to incorporate the gravitational effects of the KBOs.

It is important to note that the 11 selected objects represent only a fraction of the entire KBO population. These objects collectively possess a total mass of about |$0.007 M_{\oplus }$|⁠. While the total mass of KBOs could be an order of magnitude larger, reaching up to |$0.2M_{\oplus }$| as estimated theoretically by Pitjeva & Pitjev (2018b). Therefore, it is reasonable to suppose that the combined influence of all KBOs on Neptune’s position could be even stronger. Moreover, the spatial distribution of KBOs will also play a significant role, as we will discuss later when examining the case of Plutinos.

In this paper, we will use the Sun–Neptune distance to measure the perturbations of KBOs. However, given that the current observations are typically conducted from Earth or its vicinity, it is more customary to determine a planet’s position relative to Earth (Standish 1998; Folkner et al. 2014; Fienga et al. 2020). For example, the Earth–Mars distance has been commonly used to measure perturbations induced by the MBAs (Fienga et al. 2008) or the JTs (Li & Sun 2018). To validate our choice of using the change in the Sun–Neptune distance, with or without considering the 11 massive KBOs listed in Table 1, we also calculate the change in the Earth–Neptune distance:

where the variables have the similar meanings to equation (2), but the subscript ‘E’ indicates Earth. In Fig. 1(b), the red curve illustrates the temporal variation of |$\Delta d_{EN}$| resulting from the perturbation of these 11 KBOs. Over a period of 100 yr, the change in the Earth–Neptune distance can reach a value as large as |$\Delta d_{EN}\sim 17.7$| km, comparable to the largest |$\Delta d_{\mathrm {SN}}$| displayed in Fig. 1(a).

Upon examining the red |$\Delta d_{EN}$|-curve in Fig. 1(b), one can observe short-period oscillations on a time-scale of 1 yr, corresponding to Earth’s revolution around the Sun. To extract the low-frequency change in the Earth–Neptune distance, we applied a fast Fourier transform (FFT) low-pass filter to |$\Delta d_{EN}$|⁠. This process yields the smooth black curve, as plotted in Fig. 1(b). Comparing the black curves in Figs 1(a) and (b), we find a remarkable agreement between |$\Delta d_{\mathrm {SN}}$| and the smoothed |$\Delta d_{EN}$|⁠, with a difference between these two quantities being less than 0.4 km. Based on the results obtained here and above, we propose that the Sun–Neptune distance can serve as a representative measurement for evaluating the gravitational perturbations caused by the KBOs.

Adopting the Sun–Neptune distance as a metric offers several advantages. First, it simplifies the description of the secular variation of Neptune’s position over time by effectively avoiding the short-period changes observed in the Earth–Neptune distance due to the Earth’s annual movement. More importantly, another notable advantage is the substantial reduction in computational cost. By using the Sun–Neptune distance, there is no need to provide Earth’s position in locating Neptune. This allows us to remove the four terrestrial planets from the dynamical model of the Solar system, with Jupiter becoming the innermost planet instead of Mercury. Given the Jupiter-to-Mercury period ratio of |$\sim 20$|⁠, we can increase the integration time-step from 0.5 to 10 d (Li, Zhou & Sun 2007). As a result, the numerical integration process would be substantially faster.

In the simplified model of the Solar system, the gravitational effects of the terrestrial planets are implemented by adding their combined masses to the Sun. However, before proceeding with further investigations, we need to validate that this approximation does not substantially influence the change in the Sun–Neptune distance caused by the perturbations of the KBOs. To assess this, we once again calculate the value of |$\Delta d_{\mathrm {SN}}$| resulting from the perturbations of the 11 most massive KBOs within the simplified model. The results reveal that the profile of the time evolution of |$\Delta d_{\mathrm {SN}}$| remains nearly identical to that obtained in the complete model, which includes all eight planets. Indeed, the relative error in the largest value of |$\Delta d_{\mathrm {SN}}$| is as small as 0.019 per cent. Therefore, without compromising our main results, the simplified Solar system model can reduce the calculation expense to merely 1/20 of the original amount.

In a brief summary, our dynamical model for subsequent calculations includes the Sun and the four giant planets, optionally incorporating the KBOs. Furthermore, our attention will be directed towards a specific subset of KBOs, namely the Plutinos in the 2:3 MMR with Neptune. To evaluate the gravitational perturbation induced by the Plutinos, we will utilize the change in the Sun–Neptune distance, i.e. |$\Delta d_{\mathrm {SN}}$|⁠.

3 PERTURBATION OF THE PLUTINOS

When considering a large number of KBOs uniformly distributed in azimuth, their total perturbation should manifest as a continuous ring. By the term ‘continuous’, we mean that the acceleration of a body is computed using an asteroid ring potential (see equation 8 in Fienga et al. 2008). However, for the sake of computational simplicity, it is common to approximate the perturbation of KBOs with a discrete ring. In contrast to the continuous ring, the discrete one consists of a number of moving point masses. In Pitjeva & Pitjev (2018b) and Di Ruscio et al. (2020), 40 point masses were adopted to represent the discrete ring. Although this number may seem very small, it remains reasonable since the masses and positions of these point masses can be adjusted to fit observational data. But our current study aims to theoretically investigate the difference between the ring and arc models, so we have to minimize the impact of such discretization.

From a geometrical perspective, as the number n of point masses increases, the discrete ring approximation becomes increasingly accurate in representing the continuous ring. Similarly, as the arc is a portion of the ring, the discrete arc will also approach the continuous one as the value of n becomes large enough. Determining an appropriate number of point masses within a ring or arc allows us to achieve a balance between computational efficiency and the robustness of our simulations. It is important to note that, unlike a homogeneous ring that remains independent of time, the arcs will rotate with Neptune, meaning their azimuthal positions constantly change with time. Therefore, the use of a discrete arc model offers an advantage: the point masses can naturally revolve alongside Neptune. Another advantage is that, by assigning eccentricities to the orbits of the point masses in the discrete model, we can easily account for the influence of eccentric Plutinos.

In the following investigation of both the ring and arc models, we exclusively concentrate on the zero-inclination case. Adopting the planar model aligns with previous works that modelled KBOs using the 1D ring (Pitjeva 2010) or the 2D annulus (Pitjeva & Pitjev 2018b). Additionally, it is worth noting that Neptune’s 2:3 resonance is, in fact, an eccentric-type resonance (Li et al. 2020). The dynamical features of this resonance are primarily determined by the Plutino’s eccentricity, while its inclination has limited significance.

3.1 The ring model

Considering the ring model, we assume that the mass is distributed evenly across its azimuth. For a continuous ring, its gravitational potential exerting at the location |$(x^{\prime },y^{\prime }, z^{\prime })$| can be described as

where G is the gravitational constant, |$\rho _a$| is the linear density of the ring, and |$(r,\phi)$| are the polar coordinates of the ring relative to the Sun. Notably, the two parameters |$\rho _a$| and r are assumed to be fixed, resulting in equation (4) with only one variable, i.e. |$\phi$|⁠. Consequently, the integration can be performed using the Romberg algorithm.

To determine the linear density |$\rho _a$| of the ring, it is necessary to obtain the total mass of the Plutinos. The total mass of the entire population of KBOs was estimated to be within a wide range of |$M_{kb}=0.01-0.2M_{\oplus }$| (Pitjeva & Pitjev 2018b). Later, through an analysis of high-precision measurements of Saturn from the Cassini mission, Di Ruscio et al. (2020) proposed a more specific mass of |$M_{kb}=0.061M_{\oplus }$|⁠, with Plutinos accounting for approximately |$1/6$| of this mass. Based on these findings, we will assume a total mass of |$M_{\mathrm {plu}}=0.01M_{\oplus }$| for the Plutinos. It should be noted that, to a first-order accuracy, the perturbation caused by the Plutinos is linearly proportional to |$M_{\mathrm {plu}}$| (Li & Sun 2018). This implies that, after selecting a specific value for |$M_{\mathrm {plu}}$| in our calculations, we can easily generalize the results for the perturbation of Plutinos when their total mass is refined in the future observations. Further details on this will be discussed in Section 3.3.

Regarding the radial distance r of the ring, we should set it to be the nominal 2:3 resonance location of 39.4 au. However, in order to determine the minimum number density required for our discrete-arc model, we need to first compare the difference between a uniform continuous ring and a uniform discrete ring. In this case, choosing 39.4 au may pose a potential issue. If the point masses are distributed around 39.4 au, the uniformity of the discrete ring may be disrupted by Neptune’s 2:3 resonance. To avoid this issue when comparing the continuous and discrete ring models, we opt for a slightly larger radial distance of |$r=43.4$| au, which is 1.44 times the Neptune’s semimajor axis of |$a_N=30.1$| au. At the heliocentric distance of 43.4 au, objects with eccentricities smaller than 0.1 are not in the 8th or lower order resonances with Neptune (Li, Lawler & Lei 2023a).

We then calculate the changes in the Sun–Neptune distance induced by both the continuous and discrete rings, denoted as |$\Delta d_{\mathrm {SN}}^{ring}$| and |$\Delta d_{\mathrm {SN}}^{ring}(n)$|⁠, respectively. The difference between the perturbations of these two rings can be expressed as

where n refers to the number of point masses used to construct the discrete ring. As discussed earlier, increasing the value of n helps reduce the effects of discretization, quantified by |$\delta d(n)$| in equation (5). We compute the values of |$\delta d$| by varying n from 500 to 6000, and the results are shown in Fig. 2. It is evident that the difference |$\delta d(n)$| exhibits temporal variation due to the evolution of the system. Therefore, we parametrize this difference by its maximum value, denoted as |$\max \delta d(n)$|⁠, within the considered epoch from 2021 to 2121. For |$n=500$|⁠, the value of |$\max \delta d(n)$| can be as large as about 30 km. However, as n increases to 5000, |$\max \delta d(n)$| decreases by an order of magnitude, reaching only 3.8 km. In addition, with a further increase in n to 6000, we observe that |$\max \delta d(n)$| decreases by as little as |$\sim 0.1$| km. Even if we were to continue increasing n, it would have a limited impact on decreasing the value of |$\max \delta d(n)$|⁠. As a matter of fact, |$\max \delta d(n)\lt 4$| km measures the absolute error between the continuous and discrete ring models. Considering a value as large as 100 km for |$\Delta d_{\mathrm {SN}}^{ring}$|⁠, the relative error is just about 4 per cent.

To better visualize the small discrepancy between the perturbations induced by the continuous and discrete rings, we focus on the case of |$n=6000$|⁠. Fig. 3 depicts the temporal variations of |$\Delta d_{\mathrm {SN}}^{ring}$| and |$\Delta d_{\mathrm {SN}}^{ring}(n)$|⁠, represented by the black and red curves, respectively. We observe that the relative error between these two curves remains below 4 per cent throughout the 100-yr period. Therefore, we conclude that the discrete ring consisting of 6000 point masses is adequate to achieve a reliable approximation of the continuous ring.

By setting |$n=6000$| for the discrete ring, its number density |$\Sigma _0$| in azimuth can be calculated via

where |$2\pi$| is the radian measure of a complete ring, r is the heliocentric radius of the ring, and |$r_0=43.4$| au is the reference radius. In the subsequent investigation of the arc model, it is essential to ensure that the number density of the discrete arc satisfies the condition given in equation (6). Let l be the length of a discrete arc in azimuth, the minimum number of the point masses required to construct this arc can be calculated as

For the arc-like longitude distribution of Plutinos driven by Neptune’s 2:3 resonance, we will first determine the associated length l in the next section.

3.2 The three-arc model

As the primary focus of this paper, we again note that Plutinos, a distinct group of KBOs trapped in Neptune’s 2:3 resonance, are not uniformly distributed. Therefore, modelling them as a homogeneous ring may introduce a certain level of inaccuracy due to their unique spatial distribution. This distribution will be represented by the arc model, in which the Plutinos are arranged into some arcs. By ‘arc’, we are referring to a segment of the ring. We are now about to demonstrate the azimuthal distribution of Plutinos in physical space. For the 2:3 resonance, its critical resonant angle |$\sigma$| is defined by equation (1). According to the expression for |$\sigma$|⁠, the azimuthal distribution of Plutinos can be represented by the differences in the mean longitude between the Plutinos and Neptune, expressed as

To establish the initial conditions for the Plutinos, we uniformly select the resonant angles |$\sigma$| from |$180^{\circ }-A_{\mathrm {max}}$| to |$180^{\circ }+A_{\mathrm {max}}$|⁠, where |$A_{\mathrm {max}}$| represents their maximum resonant amplitude. And, the longitudes of perihelia |$\varpi$| are randomly chosen between 0 and |$360^{\circ }$|⁠. Then, given Neptune’s mean longitude |$\lambda _N$| at the beginning of the calculation, we can determine the initial mean longitudes |$\lambda$| of the Plutinos using equation (8). It is essential to note that when unfolding the phase space of the 2:3 resonance, three resonant islands emerge, corresponding to |$\sigma \in$| [0, |$360^{\circ }$|], [|$-360^{\circ }$|⁠, 0], and [|$360^{\circ }$|⁠, |$720^{\circ }$|], respectively (Li et al. 2023a). Therefore, besides the commonly considered range of [0, |$360^{\circ }$|], |$\sigma$| is allowed to vary within the other two ranges by associating the angle |$\sigma =\sigma \pm 360^{\circ }$|⁠. Although such displacements in |$\sigma$| preserve the resonant behaviours of Plutinos, the measurements of their longitude positions relative to Neptune, i.e. |$\Delta \lambda$|⁠, would significantly change. For example, with |$A_{\mathrm {max}}=120^{\circ }$| and an initial |$\lambda _N=0$|⁠, the initial values of |$\lambda$| are distributed within three intervals of |$[20^{\circ }, 220^{\circ }]$|⁠, |$[140^{\circ }, 340^{\circ }]$|⁠, and |$[260^{\circ }, 360^{\circ }]\bigcup [0, 100^{\circ }]$|⁠. Fig. 4 sketches the resulting azimuthal distribution of Plutinos, where the three arcs (A1–B1, A2–B2, and A3–B3) corresponds to the three respective |$\lambda$| intervals. These arcs overlap at the regions indicated by the red thicker curves. So obviously, the azimuthal distribution of Plutinos is better described using the three-arc model rather than the ring model.

For Plutinos with zero-inclination orbits, besides |$\lambda$| and |$\varpi$|⁠, two additional orbital elements, the semimajor axis a and the eccentricity e, need to be set. Initially, all Plutinos are placed at |$a=39.4$| au, corresponding to the nominal location of Neptune’s 2:3 resonance, with varying initial values of e. Bearing in mind that the initial |$\lambda$| is derived from the resonant angle |$\sigma$| via equation (8), and given that |$\sigma \in [180^{\circ }-A_{\mathrm {max}}, 180^{\circ }+A_{\mathrm {max}}]$|⁠, these Plutinos form three distinct arcs, as illustrated in Fig. 4. The azimuthal lengths of these arcs are determined by |$A_{\mathrm {max}}$|⁠, while their curvatures (equivalent to heliocentric distances of Plutinos) are influenced by the value of e. Thus, the maximum resonant amplitude |$A_{\mathrm {max}}$| and the eccentricity e of Plutinos serve as the two adjustable parameters of our three-arc model. We examine the orbits of observed Plutinos from the MPC and then choose representative values for these two parameters within the ranges of |$A_{\mathrm {max}}\le 120^{\circ }$| and |$e=0.05-0.35$|⁠. In addition, our earlier work (Li, Zhou & Sun 2014) shows that an upper limit of |$A=120^{\circ }$| is crucial for the stability of Plutinos over the age of the Solar system. Regarding the lower limit of |$e=0.05$|⁠, since the strength of the 2:3 resonance is proportional to e, Plutinos with |$e\lt 0.05$| are generally too weak to sustain their resonant behaviours. In fact, |$A_{\mathrm {max}}$| and e can, respectively, introduce asymmetries in the azimuthal and radial distributions of Plutinos, resulting in different magnitudes of changes in the Sun–Neptune distance. Interestingly, the combined effects of these two orbital parameters can possibly weaken the total perturbation of Plutinos, as we will show below.

To ensure an adequate number of Plutinos in each discrete arc to effectively replace the continuous arc, we consider the case of |$A_{\mathrm {max}}=120^{\circ }$|⁠. In this scenario, the azimuthal length of each arc reaches its maximum of |$200^{\circ }$|⁠, as depicted by the individual arcs Ai-Bi (⁠|$i=1,2,3$|⁠) in Fig. 4. Let |$l=200^{\circ }=10\pi /9$| and |$r=39.4$| au in equations (6) and (7), we can obtain that a minimum number of |$n_{\mathrm {arc}}\approx 3000$| Plutinos is needed for representing a single arc. Accordingly, in subsequent 100-yr calculations, we employ 9000 equal point masses for our three-arc model, regardless of |$A_{\mathrm {max}}$|⁠. The total mass of these point masses is still adopted to be |$M_{\mathrm {plu}}=0.01M_{\oplus }$|⁠, as we did for the ring model.

Before concluding this section, we pause to discuss the model’s centre of mass. For the continuous ring model analytically described by equation (4), the centre of mass can be fixed at the Sun (e.g. Kuchynka et al. 2010). However, for the discrete arc model consisting of a number of point masses, especially for the eccentric one, it is not feasible to simply maintain its centre of mass in the same manner. Instead, the arcs’ orbits are bound to the Sun by gravitation. A similar approach was previously used for the point-mass models in Pitjeva & Pitjev (2018b) and Di Ruscio et al. (2020). As stated at the beginning of this section, the advantage of this approach is that the arcs can naturally revolve alongside Neptune, better mimicking the motion of observed Plutinos. We acknowledge that the centre of mass of our arc model could drift, here by a value of less than 1 km. Similarly, the center of mass of observed Plutinos would also drift. In Section 4, we will see that the arc model indeed effectively represents the perturbations of observed Plutinos.

3.3 Perturbation difference between the ring and arc models

To compare the total perturbations of Plutinos modelled by the ring and arc models, we calculate the corresponding changes |$\Delta d_{\mathrm {SN}}$| in the Sun–Neptune distance. In the ring model, |$\Delta d_{\mathrm {SN}}$| is derived from the continuous ring described by equation (4), where the radius r is fixed at 39.4 au. On the other hand, in the arc model, |$\Delta d_{\mathrm {SN}}$| is obtained from three discrete arcs that we constructed in Section 3.2. These arcs are also located at 39.4 au but have varying orbital parameters |$A_{\mathrm {max}}$| and e to account for the specific distribution of Plutinos.

The temporal evolutions of |$\Delta d_{\mathrm {SN}}$| are plotted in Fig. 5, covering a time-span of 100 yr from 2021 to 2121. In each panel (a)–(d), we present the same |$\Delta d_{\mathrm {SN}}$| variation for the ring model (depicted by the black curve), which will be considered as the reference case. It is observed that the change in the Sun–Neptune distance induced by the Plutino-ring perturbation experiences a substantial increase shortly after the beginning in 2021, reaching a maximum value of |$\Delta d_{\mathrm {SN}}\sim 154$| km around the year 2100. Subsequently, |$\Delta d_{\mathrm {SN}}$| slightly decreases to about 143 km by the end of our calculation at the epoch 2121.

Next, we analyse the three-arc model, beginning with the case of |$e=0.05$|⁠, as shown in Fig. 5(a). To better understand the trend in the variation of |$\Delta d_{\mathrm {SN}}$| resulting from changes in the azimuthal distribution of Plutinos, we selectively plot colourful curves for representative resonant amplitudes of |$A \le A_{\mathrm {max}}=10^{\circ }$| (red), |$40^{\circ }$| (blue), |$80^{\circ }$| (magenta), and |$120^{\circ }$| (green). For each given |$A_{\mathrm {max}}$|⁠, the Plutinos have their resonant angles |$\sigma$| confined to the range of |$[180^{\circ }-A_{\mathrm {max}}, 180^{\circ }-A_{\mathrm {max}}]$|⁠. One can see that for Plutinos with |$A \le A_{\mathrm {max}}=10^{\circ }$|⁠, around the year 2100, |$\Delta d_{\mathrm {SN}}$| reaches its peak of 89 km, approximately half of the maximum |$\Delta d_{\mathrm {SN}}$| obtained in the ring model (indicated by the black curve). But as |$A_{\mathrm {max}}$| increases from |$10^{\circ }$| to |$120^{\circ }$|⁠, the corresponding curve (from red to green) approaches the black one associated with the ring model, indicating a diminishing asymmetry in the azimuthal distribution of Plutinos. This phenomenon can be attributed to the fact that the |$\sigma$|-range widens at larger |$A_{\mathrm {max}}$| and becomes closer to the complete interval of |$\sigma =[0, 360^{\circ }]$| as in the ring model. Nevertheless, even when |$A_{\mathrm {max}}=120^{\circ }$|⁠, the discrepancy in |$\Delta d_{\mathrm {SN}}$| between the ring (black curve) and arc (green) models can still exceed 10 per cent.

Fig. 5(b) illustrates the outcomes for another representative case of |$e=0.1$|⁠. We can see that the difference between the |$\Delta d_{\mathrm {SN}}$| values derived from the ring and arc models gradually diminishes as |$A_{\mathrm {max}}$| increases. This trend in |$\Delta d_{\mathrm {SN}}$| with increasing |$A_{\mathrm {max}}$| remains apparent and similar to the case of |$e=0.05$| (see Fig. 5a). It suggests that, at relatively small e, the asymmetry of the azimuthal distribution of Plutinos in their gravitational perturbations could be notable.

However, when the eccentricities of Plutinos are large, the asymmetry in their radial distribution becomes important and can significantly impact the Sun–Neptune distance. Figs 5(c) and (d) show the results for the cases of |$e=0.25$| and 0.3, respectively. A remarkable pattern emerges as the colourful curves representing different |$A_{\mathrm {max}}$| values cluster around |$\Delta d_{\mathrm {SN}}=0$|⁠. These high-eccentricity cases highlight our arc model in two key respects: (1) the perturbation caused by Plutinos becomes exceptionally weak due the combined effects of the azimuthal and radial distribution asymmetries of Plutinos, determined by |$A_{\mathrm {max}}$| and e, respectively. Over the entire 100-yr evolution, the resulting effect on the Sun–Neptune distance differs by only 10–20 km compared to the unperturbed Solar system model, where |$\Delta d_{\mathrm {SN}}=0$|⁠. A plausible explanation for the weak perturbation induced by Plutinos with large eccentricities could stem from their resonant configuration. When Plutinos reach their perihelia, they tend to cluster in two segments in longitude, leading and trailing Neptune by |$\sim 90^{\circ }$|⁠, respectively. As their eccentricities increase, the distances between Neptune and the Plutinos would increase during conjunctions (Chiang & Jordan 2002). Therefore, the gravitational effects of more eccentric Plutinos on Neptune become weaker. (2) Regardless of |$A_{\mathrm {max}}$|⁠, the colourful curves deviate prominently from the black curve, exhibiting differences as large as about 170 km in |$\Delta d_{\mathrm {SN}}$|⁠. This significant deviation clearly demonstrates that the ring model is completely unsuitable for accurately representing the perturbations induced by the Plutinos.

Finally, we provide a concise description for the dependence of the perturbation of Plutinos on their parameters: the maximum resonant amplitude |$A_{\mathrm {max}}$|⁠, the eccentricity e, and the total mass |$M_{\mathrm {plu}}$|⁠. Fig. 6 summarizes the changes in the Sun–Neptune distance induced by the three-arc model at the end of our 100-yr calculation, denoted as |$\Delta d^{(\mathrm {arc})}_{\mathrm {SN}}(T=100)$|⁠, for various values of |$A_{\mathrm {max}}$| and e. Through linear fitting of the data points corresponding to different |$(e, A_{\mathrm {max}})$| combinations in this figure, we establish the following measurement:

where the coefficients |$\alpha$| and |$\beta$| are given by

Furthermore, the total mass |$M_{\mathrm {plu}}$| of Plutinos may be a crucial parameter in determining |$\Delta d^{(\mathrm {arc})}_{\mathrm {SN}}$|⁠. Although the precise value of |$M_{\mathrm {plu}}$| currently remains uncertain, the first-order approximation made in previous works (Kuchynka et al. 2010; Li & Sun 2018) suggests that the perturbations induced by point masses are nearly proportional to their masses. By incorporating the mass factor into equation (9), we can derive a comprehensive expression for the perturbation of Plutinos on the Sun–Neptune distance:

To validate equation (11), we selected typical parameters of |$e=0.25$| and |$A_{\mathrm {max}}=120^{\circ }$|⁠, and varied the total mass |$M_{\mathrm {plu}}$| in the range of 0.005–0.025 |$M_{\oplus }$|⁠. First, we numerically recalculate the change in the Sun-Neptune distance, i.e. |$\Delta d_{\mathrm {SN}}$|⁠, using the three-arc model with different |$M_{\mathrm {plu}}$| values. The resulting |$\Delta d_{\mathrm {SN}}$| at the end of the 100-yr integration is indicated by the squares in Fig. 7. Then, for various total masses |$M_{\mathrm {plu}}$|⁠, we predict |$\Delta d^{(\mathrm {arc})}_{\mathrm {SN}}$| at |$T=100$| via equation (11). The prediction is depicted by the solid line in red. Notably, good agreement between the numerical results and our predictions is evident, as the absolute error is less than |$10^{-4}$| km. This error is several orders of magnitude smaller than the obtained |$\Delta d_{\mathrm {SN}}$| value or the difference with the continuous versus discrete ring modularization proposed in Section 3.1. Accordingly, the nearly linear dependence of |$\Delta d_{\mathrm {SN}}$| on |$M_{\mathrm {plu}}$|⁠, as described in equation (11), can be nicely supported.

At this point, we propose that the parameter |$M_{\mathrm {plu}}$| can be decoupled from the other parameters (i.e. |$e, A_{\mathrm {max}}$|⁠) in our arc model. This implies that when the total mass |$M_{\mathrm {plu}}$| of Plutinos is refined in future surveys, one can easily estimate their influence on the Sun–Neptune distance using a linear relation in |$M_{\mathrm {plu}}$| instead of conducting repetitive numerical computations.

4 COMPARISON WITH OBSERVED POINT-MASS PLUTINOS

As we proposed, the impact of the eccentric Plutino orbits has to be modeled by a non-uniformly distributed mass, i.e. the eccentric arcs, over a circular ring. Naturally, we need to verify that the three-arc model is sufficient to represent the perturbations induced by the individual contributions of observed Plutinos on eccentric orbits. To do so, we will compare the positions of the planets perturbed by the three-arc model and those perturbed by the point-mass eccentric Plutinos from observations.

As of May 2024, over 3200 KBOs have been registered in the MPC database.² For the sake of identifying the Plutinos, we first consider the KBOs with semimajor axes of |$39\le a \le 40$| au. By this criterion, 519 KBOs are selected to be Plutino candidates. Then the orbits of these candidates are numerically integrated in a Solar system model consisting of the Sun and four Jovian planets. Finally, an object is regarded as a Plutino if its resonant angle |$\sigma$| librates, i.e. the resonant amplitude of |$A\lt 180^\circ$|⁠. The total duration of the integration is chosen to be 1 Myr. We note that the typical libration period of the 2:3 resonant angle |$\sigma$| is about 20 000 yr; for example, Pluto has a resonant period of 19 670 yr (Cohen & Hubbard 1965; Murray & Dermott 2000). Therefore, the integration time of 1 Myr is approximately equal to 50 resonant periods, which is long enough to identify the Plutinos that can stably exhibit the libration behaviour of |$\sigma$|⁠.

Eventually, 460 Plutinos have been identified, and their eccentricity distribution is shown in Fig. 8. Considering that the arc model shows a greater difference from the ring model at high eccentricity, we will focus on the eccentricity range of |$e\gt 0.2$| below. There are 278 Plutinos that belong to this e range, accounting for more than 60 per cent of the entire Plutino population. It is evident that the number of currently observed Plutinos is much smaller than the theoretically estimated population, which is over 9000 for diameters |$D\gt 100$| km, as mentioned in Section 1. This number deficit is because surveys for KBOs are usually limited in relatively small areas (⁠|$\sim 100$| square degrees) due to their faint magnitudes (⁠|$m_R\sim 23$| with |$D\sim 100$| km) (Bannister et al. 2018; Bernardinelli et al. 2020; Trujillo et al. 2022). Therefore, the observational incompleteness of the available observed Plutinos can lead to a biased spatial distribution. The spatial distribution of Plutinos is crucial as it significantly affects the estimation of their total perturbation. Furthermore, for the Plutinos with |$e\gt 0.2$|⁠, they are more likely to be observed when they reach their perihelia and are at their brightest. However, given the orbital period of the Plutino population, which is about 250 yr long, those objects that are near their aphelia, i.e. at relatively larger heliocentric distance, would be more difficult to observe currently and in the near future.

In order to offset the observational bias, we have constructed two methods to simulate the realistic spatial distribution of the overall Plutinos:

(1) Dispersion in space (DISpa): Since the number of Plutinos currently observed is much less than theoretically predicted, we enlarge the sample size following the idea of the Monte Carlo method (Li et al. 2019; Li & Xia 2020). It is obvious that the spatial dispersion of the Plutino samples strongly depends on the three angles of their orbital elements: longitude of ascending node (⁠|$\Omega$|⁠), argument of perihelion (⁠|$\omega$|⁠), and mean anomaly (M). Thus, for each of the 278 Plutinos with |$e\gt 0.2$|⁠, we generate 10 synthetic samples with the same a, e, and i, as well as with the same resonant angle |$\sigma$|⁠, but their |$\Omega$| and |$\omega$| are randomly selected between |$0^{\circ }$| and |$360^{\circ }$|⁠. The mean anomaly M of a synthetic sample is determined according to equation (1). In this method, we produced a more ‘complete’ Plutino population, consisting of 2780 synthetic Plutino samples, referred to as DISpa samples.

(2) Dispersion in time (DITim): Regarding the current 278 Plutinos with |$e\gt 0.2$|⁠, we numerically integrated their trajectories for |$10^7$| yr, under the gravitational perturbations of four Jovian planets. Since the |$10^7$|-yr evolution time is much longer than the secular procession period of the orbits of Plutinos (Li et al. 2014), this method can naturally disperse the observed Plutinos within the 6D orbital distribution, even with a sample size as small as 278. These samples will be denoted by the DITim samples.

It is worth mentioning that the inclinations of the 278 considered Plutinos are taken to be the observed values, meaning that both the DISpa and DITim samples generated from this population follow an inclination distribution similar to the observations. This may further allow us to test the validity of the theoretical arc model, which is co-planar in principle.

Subsequently, we calculated the perturbation on the Sun–Neptune distance induced by the point-mass Plutinos from the two methods described above.The results are shown by the red curves in Fig. 9, with the left and right panels corresponding to the DISpa and DITim samples, respectively. This figure also displays the Sun–Neptune distance difference |$\Delta d_{\mathrm {SN}}$| induced by the ring model (black curve) and the arc models with eccentricities |$e=0.2$|⁠, 0.25, and 0.3 (indicated by different shades of blue). It can be seen that, in general, the time evolution of |$\Delta d_{\mathrm {SN}}$| for the point-mass Plutinos deviates significantly from that of the ring model, but it is much closer to that of the arc models. Specifically, the arc model with |$e=0.25$| nicely mimics the point-mass Plutinos, as the corresponding regular blue and black curves have a difference in |$\Delta d_{\mathrm {SN}}$| of less than 2 km for the first 50 yr. This could be understood by acknowledging that, as shown in Fig. 8, the eccentricity |$e=0.25$| corresponds approximately to the average eccentricity of the considered Plutinos with |$e\gt 0.2$|⁠.

However, after 50 yr of evolution, the difference in |$\Delta d_{\mathrm {SN}}$| between the cases of the point-mass Plutinos and the arc model with |$e=0.25$| increases, reaching about 20 and 40 km for the DISpa and DITim samples, respectively. This noticeable discrepancy may be due to the fact that observational bias has not been completely eliminated. Nevertheless, we find that the difference in |$\Delta d_{\mathrm {SN}}$| between the cases of the point-mass Plutinos and the ring model is much larger, exceeding 140 km. We note that observations of the farthest planet, Neptune, primarily rely on astronomical observations. Due to limitations imposed by Earth’s atmosphere and the accuracy of star catalogues, the orbital accuracy of Neptune is currently at the level of a few thousand kilometers (Stone 2001; Folkner et al. 2014; Pitjeva & Pitjev 2018b). While the 140 km difference in Neptune’s position that we derived may seem small compared to current observations, it is noteworthy because a similar magnitude of orbital variation, |$\Delta d_{\mathrm {SN}}$|⁠, will be introduced once the ring model developed in previous works incorporates the planetary ephemerides. Furthermore, our arc model will become increasingly valuable as future in situ observations lead to significant improvements in the positional accuracy of Neptune.

Having already taken into account the comparison of the 3-arc model to the observed Plutinos, represented by the DISpa and DITim samples. As a supplement, a new experiment has been conducted, letting the observed Plutinos be represented by the population of 2:3 resonant KBOs from the L7 synthetic model (Kavelaars et al. 2009; Petit et al. 2011; Gladman et al. 2012). The L7 synthetic model, based on the Canada France Ecliptic Plane Survey (CFEPS), gives a distribution of KBOs extrapolated down to small sizes, with |$H_g$| magnitude of 8.5. Consistent with the DISpa and DITim samples, here we consider the L7 model samples with |$e\gt 0.2$|⁠, resulting in a total number of about 1500 objects. These objects were then used to calculate the perturbation on the Sun–Neptune distance, and the results are shown in Fig. 9 (see the light magenta curves). We note that, in relation to the right panel, we consider the L7 model samples after a |$10^7$|-yr evolution. In general, the |$\Delta d_{\mathrm {SN}}$|-curve resulting from the L7 model is similar to that from the DISpa (or DITim) samples, tending to align more closely with the 3-arc model while diverging significantly from the ring model. One may notice a slight discrepancy: at an earlier epoch of about 30 yr, the time evolution of |$\Delta d_{\mathrm {SN}}$| for the L7 model starts to appear different from that for the 3-arc model. This discrepancy could be accounted for by the L7 synthetic model being derived solely from the CFEPS samples, which contain only 24 observed Plutinos. In fact, the CFEPS project website³ states that the L7 model may not represent the intrinsic orbit distribution of the KBOs, and direct comparisons with other KBO models should be avoided. Nevertheless, it is worth remarking that even within this 30-yr window, when referring to the L7 model, the ring model induces a maximum difference in |$\Delta d_{\mathrm {SN}}$| of over 40 km, which is ten times larger than the difference of about 4 km observed with the 3-arc model. This suggests that the 3-arc model is still a better representation of the observed Plutinos compared to the ring model.

Throughout the above analysis, three points need to be noted. First, to maintain consistency with the ring and arc models considered in Section 3, we assume the total mass of the eccentric Plutino samples to be equivalent to that of the entire Plutino population, i.e. |$0.01M_{\oplus }$|⁠. In reality, however, the former should be less massive than the latter. Since |$\Delta d_{\mathrm {SN}}$| nearly linearly depends on the total mass of Plutinos, a smaller total mass could further reduce the absolute discrepancy between the point-mass Plutinos and our arc model. Secondly, in the construction of the DISpa samples, the initial conditions of the system are adjusted to be at the epoch of the year 2024, which is the same as the epoch of the KBOs registered in the MPC data base. Therefore, the initial epoch in Fig. 9 is 3 yr later than that in Fig. 5. Thirdly, the x-coordinates are different in the left and right panels in Fig. 9. This is due to the distinct initial epochs of the modelled planetary systems. As mentioned above, in the left panel, the initial epoch of the DISpa samples – in the year 2024 – represents the current time of the observed Plutinos. Meanwhile, in the right panel, the DISpa samples refer to the observed Plutinos after a |$10^7$|-yr evolution, meaning the initial epoch corresponds to |$10^7$|-yr into the future. Consequently, since the arc models in these two panels refer to planetary systems at different epochs, the resulting profiles of the |$\Delta d_{\mathrm {SN}}$| curves are not the same. We believe this epoch difference can further support the comparison between the point-mass Plutinos, the arc model, and the ring model. A similar rationale also applies to Fig. 10.

In a brief summary, we believe that our three-arc model has a significant advantage in representing the total perturbation induced by observed Plutinos, especially for those on the eccentric orbits, which comprise a majority of the Plutino population.

Having constructed the Plutino models, we apply them to evaluate the impact of the Plutinos on Neptune’s longitude and latitude. This could be useful because observational constraints for Neptune, the farthest planet, are often given in these two parameters. By making simple coordinate transformations on our previous results from the ring model, the arc model, and the DISpa Plutino samples, in Table 2, we provide the maximum changes in Neptune’s longitude (⁠|$\Delta \lambda _N$|⁠) and latitude (⁠|$\Delta \phi _N$|⁠). For reference, the corresponding maximum changes in the Sun–Neptune distance (⁠|$\Delta d_{\mathrm {SN}}$|⁠) are also included. The numbers in Table 2 are measured within the first 50 yr of evolution, rather than a century. This shorter time-scale is chosen because, as shown in Fig. 9, the arc model better mimics the observed Plutinos during the first 50 yr; afterwards, the difference in |$\Delta d_{\mathrm {SN}}$| may become more evident. One may notice that, in comparison with the observational accuracy of geocentric longitudes and latitudes, the numbers given in Table 2 appear pretty small for all the modelling. However, as we argued above, this does not undermine the applicability of our arc model relative to the ring model, as the latter could yield much larger differences in Neptune’s position when compared to the observed Plutinos. Additionally, the consideration of Saturn’s positions will further validate the arc model, as discussed in the following section.

4.1 Perturbation on Saturn positions

Compared to observations of Neptune, which rely mainly on astronomical observations, the accuracy of Saturn positions has been significantly enhanced by in situ observations. Based on a series of astrometric very long baseline interferometry (VLBI) observations of the Cassini spacecraft, Jones et al. (2011) pointed out that the positional accuracy for Saturn is about 2 km. Later, based on Cassini tracking data over the period between 2004 and 2017, the accuracy of the geocentric positions of Saturn was considerably improved to about 0.025 km (Fienga et al. 2021; Park et al. 2021). Therefore, present high-accuracy observations of Saturn positions may further verify the applicability of our three-arc model.

Results are shown in Fig. 10, comparing the point-mass Plutinos, the arc model, and the ring model as in Fig. 9, but this time plotted for the induced changes in the Sun–Saturn distance, denoted by |$\Delta d_{SS}$|⁠. Similar to Fig. 9, this figure also shows that the time evolutions of |$\Delta d_{SS}$| for the point-mass Plutinos (red curve) and the arc models (blue curves) generally match well over the 100-yr evolution. However, there are considerable differences when compared to the ring model (black curve), in both the exact values and the curves’ profiles.

Specifically, as seen in the left panel of Fig. 10, when we consider the perturbation of the DISpa Plutino samples as the standard, the maximum difference in |$\Delta d_{SS}$| for the arc model with |$e=0.25$| is only about 0.96 km, while the difference for the ring model is more than twice as large, at approximately 2.15 km. Comparable results are achieved when taking the perturbation of the DITim Plutino samples as the standard, as shown in the right panel of Fig. 10. The maximum differences in |$\Delta d_{SS}$| are 0.87 and 1.80 km for the arc model and ring model, respectively, with the latter being over twice as large. We again note that the different profiles of the corresponding |$\Delta d_{SS}$|-curves (e.g. the red curves) between the two panels in Fig. 10 are also due to the different initial epochs of the DISpa and DITim samples, as illustrated for Fig. 9. Nevertheless, the comparisons among various models would not be affected.

Furthermore, to compare with the VLBI and Cassini observations, Fig. 11 presents a zoom-in view of Fig. 10 on the evolution of |$\Delta d_{SS}$| over time. The 15-yr period shown corresponds approximately to the interval of Cassini observations, as we mentioned above. Notably, during this period, the Sun–Saturn distance change, |$\Delta d_{SS}$|⁠, has already reached a local maximum once. When comparing with the DISpa samples, as shown in the left panel of Fig. 11, the maximum difference in |$\Delta d_{SS}$| induced by the arc model with |$e=0.25$| is about 0.82 km, which is less than half of the 1.69 km difference produced by the ring model. And, the arc model exhibits even greater accuracy when compared with the DITim samples (right panel of Fig. 11), yielding a maximum |$\Delta d_{SS}$|-difference of only about 0.20 km. We additionally note that, although the L7 synthetic model was previously advised against direct comparison with other KBO models, its perturbation on Saturn’s positions – indicated by the light magenta curve in Fig. 11 – also seems to match the arc model (blue curves) quite well, with maximum differences in |$\Delta d_{SS}$| ranging from 0.2 to 0.3 km. In contrast, the |$\Delta d_{SS}$| difference between the L7 synthetic model and ring model consistently exceeds 1 km. Thus, the L7 synthetic model may further emphasize the significance of incorporating the arc model into planetary ephemerides instead of the ring model.

As our previous consideration of the changes in Neptune’s longitude and latitude, Table 3 presents similar results for Saturn, covering both the 50 yr and 15 yr evolutionary time-scales. Over 50 yr of evolution, when compared to the point-mass Plutinos represented by either the DISpa samples or the L7 synthetic model, the arc model provides a noticeably better fit than the ring model for Saturn’s longitude (⁠|$\Delta \lambda _S$|⁠), latitude (⁠|$\Delta \phi _S$|⁠), and the Sun–Saturn distance (⁠|$\Delta d_{SS}$|⁠). While for the 15-yr evolution, due to the shorter time span, the differences in |$\Delta \lambda _S$| and |$\Delta \phi _S$| are not significant, remaining at the level of |$\lesssim 0.3$| mas. However, the maximum change in |$\Delta d_{SS}$| induced by the ring model still differs from other models by more than 0.8 km, which is much larger than the current observational accuracy of about 0.025 km for the geocentric Saturn distance.

In summary, when predicting Saturn’s heliocentric distance, compared to the observed Plutinos, the arc model provides a better representation, as it results in a discrepancy of about 1 km smaller than that of the ring model. Given that the observational accuracy of the geocentric positions of Saturn is as small as |$\sim 0.025$| km, we suggest that future ephemeris computations, when modeling Plutinos, take into account their unique distribution driven by Neptune’s 2:3 MMR, rather than assuming a homogeneous ring model.

5 CONCLUSIONS AND DISCUSSION

This paper is devoted to developing an appropriate dynamical model to represent the Plutinos, which are KBOs trapped in the 2:3 MMR with Neptune, and to estimate their total perturbation by including the peculiar features of this MMR. Previously, the total perturbation of Plutinos was simulated using a 1D, homogeneous ring model. However, since the 2:3 MMR involves three stable equilibrium points, Plutinos are expected to exhibit a unique spatial distribution. Consequently, the ring model may introduce a certain level of inaccuracy.

We analysed the potential azimuthal distribution of Plutinos resulting from the libration of the resonant angle |$\sigma$| of Neptune’s 2:3 MMR. Since |$\sigma$| does not complete a full circulation from 0 to |$360^{\circ }$| (i.e. the resonant amplitude |$A\lt 180^{\circ }$|⁠), the mean longitudes |$\lambda$| of Plutinos relative to Neptune will be confined to three overlapping arcs. The length of each arc is equal and determined by the maximum resonant amplitude |$A_{\mathrm {max}}$|⁠, with |$120^{\circ }$| serving as an upper limit for stable Plutinos. To capture this azimuthal distribution, we introduced a new three-arc model to represent the perturbations induced by Plutinos. Each arc is represented by a number of point masses, and a minimum number density of |$6000/2\pi$| per radian in the azimuth has been deduced for effectively representing their gravitational perturbations. Accordingly, our three-arc model comprises a total of 9000 point masses, each arc possessing an equal mass.

We then assessed the total perturbation of Plutinos by calculating the change in the Sun–Neptune distance, |$\Delta d_{\mathrm {SN}}$|⁠, between the perturbed (including Plutinos) and unperturbed (excluding Plutinos) Solar systems over the 100-yr period. For comparison, we employed both the ring and arc models to simulate the perturbations of Plutinos. Given Plutinos having eccentricities of |$e=0.05$|⁠, our three-arc model with |$A_{\mathrm {max}}=10^{\circ }$| yielded a peak |$\Delta d_{\mathrm {SN}}$| of 89 km. In contrast, the reference case of a homogeneous ring resulted in a significantly larger |$\Delta d_{\mathrm {SN}}$|⁠, reaching approximately 154 km. This discrepancy is anticipated since a small |$A_{\mathrm {max}}$| can lead to a severe asymmetry in the azimuthal distribution of Plutinos. Nevertheless, this asymmetry noticeably decreases when considering the largest |$A_{\mathrm {max}}=120^{\circ }$|⁠, with the relative difference in |$\Delta d_{\mathrm {SN}}$| between the arc and ring models being about 10 per cent over the 100-yr period.

As a matter of fact, the azimuthal distribution of Plutino arcs is governed by |$A_{\mathrm {max}}$|⁠, while their radial distribution depends on e. It is interesting to discover that, with increasing e from 0.05 to 0.1–0.2, the difference in |$\Delta d_{\mathrm {SN}}$| between the arc and ring model becomes more prominent. Notably, for the largest |$e \gtrsim 0.25$|⁠, we emphasize two key findings: (1) the perturbation caused by Plutinos on the Sun–Neptune distance could consistently remain small for any |$A_{\mathrm {max}}$|⁠, with the resulting |$\Delta d_{\mathrm {SN}}$| being only 10–20 km within the considered 100-yr period. (2) The |$\Delta d_{\mathrm {SN}}$| value obtained in the arc model can deviate by as much as 170 km from that in the ring model.

To summarize the above results, we have provided a concise analytical expression in equation (11) to estimate the change in the Sun–Neptune distance caused by the perturbations of Plutinos. This expression takes into account the resonant amplitudes, eccentricities, and total mass of Plutinos.

Subsequently, we further verified the applicability of our three-arc model by comparing it to the perturbations induced by observed Plutinos. After reducing observational bias by simulating the realistic spatial distribution of the Plutinos, the results show that the arc model can mimic Plutinos very well, especially in the first 50 yr of evolution, characterized by a different in |$\Delta d_{\mathrm {SN}}$| as small as |$\lt 2$| km. Considering that the observational accuracy of Saturn positions is much higher and could better constrain the Plutino distribution, we then calculated the change in the Sun–Saturn distance |$\Delta d_{SS}$| induced by different perturbation models. Taking the perturbation of observed Plutinos as the standard, we find that the difference in |$\Delta d_{SS}$| for the arc model is of the order of 1 km, while for the ring model it is about twice as large, at approximately 2 km. Given the |$\sim 0.025$| km observational accuracy of the geocentric positions of Saturn, we propose that our arc model is more appropriate for representing the observed Plutinos. In addition, we re-examined the comparison between the continuous and discrete ring models and found that the difference in |$\Delta d_{SS}$| between these two models is only about 0.03 km. This discrepancy is notably smaller than the differences between the considered models. Specifically, the ring model consistently differs from other models by more than 1 km in |$\Delta d_{SS}$|⁠. Thus, we have confirmed that the discrete model provides an adequate approximation of the continuous one.

As a supplement, we also considered a separate population of point-mass Plutinos from the L7 synthetic model (Kavelaars et al. 2009; Petit et al. 2011; Gladman et al. 2012). When comparing the perturbations induced by the L7 Plutinos on both Neptune’s and Saturn’s positions, the differences observed in the arc model are consistently much smaller than those in the ring model. This may further emphasize the importance of incorporating the arc model into planetary ephemerides rather than the ring model.

Although the obtained results suggest that our three-arc model is more appropriate as it takes into account the A and e distributions of Plutinos, there are still several aspects that can be further refined:

(1) The arc model could be somewhat rough due to the limited observations of observed Plutinos. In our theoretical study, although the arcs are confined to azimuthal ranges narrower than |$0^{\circ }$| to |$360^{\circ }$|⁠, they are indeed evenly distributed within these specific ranges. We identified that currently only about 460 Plutinos have been observed, while there may be as many as 9000 Plutinos with diameters larger than 100 km (Alexandersen et al. 2016). Although we tried to reduce observational bias, a more intrinsic arc model could be constructed for the Plutino population from future observations. For example, the Plutinos may be sparse in certain parts of the arcs while being dense in others. This potential inhomogeneity in the distribution of Plutinos would further diminish the applicability of the ring model. To overcome this issue, we may develop our arc model by considering a series of sub-arcs, each with different |$A_{\mathrm {max}}$| values, embedded in the three main arcs.

(2) We considered the point masses in the arc model with orbital inclinations of |$i=0$|⁠. As motioned in Section 1, this simplicity is reasonable for our initial study because the 2:3 MMR is an eccentricity-type resonance, and its feature is controlled by eccentricity rather than inclination. However, observations show that a significant fraction of observed Plutinos has |$i\sim 5^{\circ }$|–|$35^{\circ }$|⁠. It should be noted that, in Section 4, when evaluating the influence on either the Sun–Neptune or Sun–Saturn distance, the point-mass Plutinos are not distributed in the plane, but have an inclination distribution similar to the observed Plutinos. Given the inclined Plutinos, the planer arc model can still fit their perturbations well. Nevertheless, a more advanced arc model could be developed for higher accuracy, which would involve multiple layers, designed to accommodate various i values. Within each layer, a 1D arc model would be implemented.

In recent years, several projects for surveying the KBOs have been proposed, such as the ecliptic Deep Drilling Field by the Large Synoptic Survey Telescope (LSST; Jones, Jurić & Ivezić 2015; Trilling et al. 2018). The discovery of a larger number of small and faint KBOs is anticipated, enabling us to impose even tighter constraints on the spatial distribution and the total mass of Plutinos. Consequently, a more refined KBO model can further be developed.

Finally, we note that the positions of Neptune are mainly obtained through astronomical observations. Due to limitations posed by the Earth’s atmosphere and the accuracy of star catalogs, Neptune’s positional accuracy is currently at the level of a few thousand kilometers (Stone 2001; Folkner et al. 2014; Pitjeva & Pitjev 2018b). While the difference in computed positions of Neptune between the ring and arc models is a few hundred kilometers, which appears smaller than the observational uncertainty. But, it’s crucial to consider that the perturbation of Plutinos on Neptune’s position is at a similar magnitude, around one hundred kilometers. Therefore, if the improvement of the ephemeris is to be carried out by the individual population of Plutinos, it may be necessary to adopt the arc model to represent their total perturbation. The improvement of accuracy for the positions of other planets (e.g. Saturn) may also leads to such a need. Moreover, as introduced in Section 1, there are several other resonant populations among the KBOs, with the number of each population as large as that of the Plutinos. Our multiple-arc model can be easily generated for these resonant KBOs. And, the combined gravitational effect of all the resonant KBOs should be more noticeable. We anticipate future space missions, such as NASA’s Uranus Orbiter and Probe mission, to greatly improve the positional accuracy of Neptune. After the observational improvement, the contribution of our arc model in developing the ephemeris calculation for resonant KBOs may lead to more accurate predictions of Neptune’s positions.

ACKNOWLEDGEMENTS

This work was supported by the National Natural Science Foundation of China (Nos. 12473061, 11973027, 11933001, 12150009), and National Key R&D Program of China (2019YFA0706601). We would also like to express our sincere thanks to the anonymous referee for the valuable comments.

DATA AVAILABILITY

The data underlying this article are available in the article and in its online supplementary material.

Footnotes

REFERENCES

APPENDIX A: EARTH POSITIONS IN THE BARYCENTRIC COORDINATE SYSTEM

The introduction of the Plutinos in planetary ephemerides will induce a drift in the SSB, which is defined by the centre of mass in each proposed model. Although the relative coordinates between two celestial bodies, such as Sun–planet or planet–planet, remain the same across different models (Pitjeva & Pitjev 2014), it is also very interesting to provide a glimpse of the SSB-Earth positions.

In Fig. A1, we show the effect of three representative models on the SSB-Earth distance (⁠|$\Delta d_{\mathrm {SSB}-\mathrm {Earth}}$|⁠), represented by the solid curves. The amplitudes of the |$\Delta d_{\mathrm {SSB}-\mathrm {Earth}}$| variations consistently range between 2 and 4 km. Nevertheless, one may notice that the |$\Delta d_{\mathrm {SSB}-\mathrm {Earth}}$| values induced by the DISpa samples (in red) are about 30 km larger than those of the arc (blue) and ring (in black) models. This difference is due to the drift of the SSB in these models, as indicated by the horizontal dashed lines. It is worth noting again that the inclusion of all asteroids in the Solar system can cause a displacement of the SSB by as much as 100 km (Li et al. 2019).

With respect to Earth's longitude (⁠|$\Delta _{\mathrm {Longitude}}$|⁠) and latitude (⁠|$\Delta _{\mathrm {Latitude}}$|⁠) in the barycentric coordinate system, Fig. A2 shows that the variations in these angular coordinates are less than 1 mas in both the arc and ring models. While for the point-mass Plutinos generated by the DISpa method, the maximum changes in |$\Delta _{\mathrm {Longitude}}$| and |$\Delta _{\mathrm {Latitude}}$| reach approximately 40 mas and 20 mas, respectively. We believe these relatively large changes can be attributed to the fact that the observational bias has not been completely eliminated, as discussed in the main text. According to the calculations in Li et al. (2019), the contribution of the faintest and least massive asteroids can change the ascending node of the invariable plane of the Solar system in the barycentric frame by about 40 mas. This seems to align with the results concerning the impact of observational bias.