On the origin of extreme trans-Neptunian objects within Modified Newtonian Dynamics Free

ABSTRACT

In this work, we investigate the dynamical origin of extreme trans-Neptunian objects (ETNOs) under the action of the External Field Effect (EFE), which is a consequence of Modified Newtonian Dynamics (MOND) applied to gravity around the Sun embedded in the gravitational field of the Galaxy. We perform N-body integrations of known ETNOs treated as massless particles and perturbed by four giant planets and EFE. Backward integrations show that these objects originated in the giant planet region, from where they were scattered and then evolved to their current orbits. A striking example of such evolution is Sedna, which may have been temporarily in a horseshoe orbit with Jupiter and Saturn only 30 Myr ago. Another interesting example is the newly discovered retrograde ETNOs, whose dynamical connection with prograde ETNOs and Centaurs is shown. The EFE is considered as an alternative to Planet Nine in explaining the anomalous distribution of ETNO orbits, namely the orbital plane clustering and apsidal confinement. We also analyse the effect of MOND on the obliquity of the solar spin with respect to the invariant plane of the solar system.

1 INTRODUCTION

The existence of an unobserved distant planet in the solar system, referred to as Planet Nine, has been proposed in response to the unexplained structure of the trans-Neptunian solar system. The anomalous features are the perihelion detachment of Sedna (Brown, Trujillo & Rabinowitz 2004) and other objects with perihelion distances |$q$| beyond the gravitational influence of Neptune, highly inclined or even retrograde (the inclination |$i \gt 90~$|°) orbits of Centaurs (Batygin & Brown 2016b), the apsidal confinement as well as the clustering of orbital planes of extreme trans-Neptunian objects (ETNOs) with |$a \gtrsim 150 \ \mbox{au}$| and |$q \gt 30 \ \mbox{au}$| (Trujillo & Sheppard 2014) or, based on better statistics, for |$a \gtrsim 250 \ \mbox{au}$| and |$q \gtrsim 42 \ \mbox{au}$| (Batygin & Brown 2016a; Brown & Batygin 2021). Following the discovery of new objects, the apsidal confinement proved to be bimodal (Sheppard & Trujillo 2016), with the two modes differing by |$180~$|°. For a thorough review of the Planet Nine hypothesis, see Batygin et al. (2019).

These features, however, are not yet fully confirmed. While Brown (2017) and Brown & Batygin (2019) argue that the probability that the clustering is due to random chance combined with observational bias is very low, Bernardinelli et al. (2020) and Napier et al. (2021) point out that there is no evidence for the clustering of the observed orbits and Clement & Kaib (2020) conclude that further new discoveries of ETNOs are needed to confirm the necessity of an additional planet. We do not attempt to settle the dispute in this article, and we treat the orbital plane clustering and the apsidal confinement as presumed anomalous features of the ETNO orbits. However, even if they are not real, two facts remain problematic for solar system formation scenarios, namely the existence of ETNOs with large |$q$| and the high inclinations of both the Centaurs and the recently discovered retrograde ETNOs with |$q \gt 50 \ \mbox{au}$| and |$i \gt 160~$|°.

However, an additional planet on the periphery of the solar system is not the only possible explanation for the anomalous features of trans-Neptunian solar system. Modified Newtonian Dynamics (MOND; Milgrom 1983a) was originally proposed as an alternative to dark matter to explain the flat rotation curves of galaxies. However, it turned out that MOND is able to explain or even predict other properties of galaxies. The proposed gravity model explains the Tully–Fisher law (Milgrom 1983b; Tully & Fisher 1977) with a power-law coefficient consistent with observations (Milgrom 1983b; Torres-Flores et al. 2011). It was also predicted that there would be a significant discrepancy between the dynamical mass and the luminous mass for galaxies with low surface brightness (Milgrom 1983b; de Blok & McGaugh 1997). It was also shown that the discrepancy between the centripetal acceleration and the gravitational acceleration caused by the luminous matter measured for many galaxies has regularities as expected from MOND (Milgrom 1983b; Lelli et al. 2017).

The original formulation of MOND (Milgrom 1983a) was further developed into a potential theory (Bekenstein & Milgrom 1984; Milgrom 2010) that allows dynamical studies of disc galaxies. It was shown that MOND can solve the problem of the stability of disc galaxies (Milgrom 1989; Brada & Milgrom 1999). For a long time, MOND had a problem with the complete description of the motions of galaxies within galaxy clusters (Sanders 1999), which brought it some criticism. However, recent studies (López-Corredoira et al. 2022) have shown that the discrepancy between luminous mass and virial mass in galaxy clusters can be explained with MOND if appropriate assumptions about hydrostatic equilibrium are used. For a thorough review of dark matter and MOND, see Sanders (2010) and Famaey & McGaugh (2012).

However, the effects of MOND are not limited to galactic and extra-galactic scales, they are important close to the Sun as well. The External Field Effect (EFE; Milgrom 1983a) appears when a massive body, e.g. Sun, is embedded in an external gravitational field, e.g. the field of the Galaxy. It was shown in Milgrom (2009) that the effect can play a role in the Solar system even in the giant planets region. Even more significant is the effect in the trans-Neptunian region. Paučo & Klačka (2016) considered formation of ETNOs as well as highly inclined Centaurs within EFE and found that such objects can be naturally formed. They also demonstrated that both the groups belong to the same trans-Neptunian objects, TNO,¹ population, but in different stages of their evolution.

The other two anomalous features of the ETNOs proved more difficult to explain in the context of EFE. Paučo (2017) showed that the reconstruction of the apsidal confinement as well as the clustering of the orbital planes would require an improbable value of the EFE quadrupole strength parameter, |$Q$|⁠, with respect to the Cassini data (Blanchet & Novak 2011; Hees et al. 2014, 2016).

In this paper we reconsider the possible explanation of the ETNO characteristics within EFE. Our approach differs in some respects from that presented in Paučo (2017). We compute the EFE perturbing acceleration by directly solving the modified Poisson equation (Milgrom 2010) instead of using an analytical model as in Paučo (2017). We include all giant planets in the |$N$|-body model of motion, not just Neptune as in Paučo (2017). We have mainly focused on known objects rather than synthetic populations.

The reason for the latter is that a comparison between the observed and the synthetic distribution of the orbital elements is problematic. The first difficulty is that the final synthetic distribution depends on the initial distribution, which is not known. The second difficulty is the epoch of comparison. Since the EFE-induced acceleration depends on the position of the Galaxy centre in the ecliptic reference frame, one should compare the synthetic distribution with the observations only in the vicinity of the current phase of the solar motion in the Galaxy. For example, it is possible that the orbits are not clustered around a fixed plane, but around the plane whose orientation changes as the Sun orbits around the centre of the Galaxy.

Through forward and backward integrations of known Centaurs and ETNOs, we have shown that the transition between the two groups of objects occurs on typical time-scales of tens to hundreds of Myr, but can even be as short as a few Myr for the widest orbits. Their evolution is chaotic due to close encounters with all giant planets. The simulations show that most of these objects in the past had semimajor axes |$a$| below |$\sim 30 \ \mbox{au}$| as well as moderate eccentricities |$e$| and inclinations |$i$|⁠. This suggests that the ETNOs and Centaurs originate from the unstable region of the solar system, with their initial orbits lying between the orbits of the giant planets. The objects may also return to this region from time to time.

In addition, we investigate the possible origin of the three newly discovered retrograde ETNOs (2022 FN12, 2022 FM12, 2019 EE6) whose |$a \sim 150 \ \mbox{au}$|⁠, |$q \gt 50 \ \mbox{au}$|⁠, and |$i \gt 160\,$|° within EFE, showing that they had lower inclinations as well as |$q \lt 30 \ \mbox{au}$| in the past. This shows that they belong to the same population as Centaurs and prograde ETNOs.

The paper is structured as follows. In Section 2 we compute the EFE-induced acceleration for the Sun embedded in the gravitational field of the Galaxy. In Section 3 we analyse the evolution of the test particles in a simple model without the planetary perturbations and with a fixed direction towards the centre of the Galaxy. In Section 4 we discuss the distribution of orbital elements of TNOs in the context of the Planet Nine hypothesis. In Section 5 we introduce the |$N$|-body model of motion and examine the evolution of the known TNOs. In Section 6 we analyse the clustering of orbits and apsidal confinement of ETNOs within EFE. Section 7 is devoted to the discussion of the possible existence of Planet Nine within MOND, the influence of EFE on the obliquity of the solar spin with respect to the invariant plane of the solar system, and possible implications of the results presented in this paper on the dark matter problem. The conclusions can be found in Section 8.

2 EFE FOR THE SUN–GALAXY SYSTEM

We use a quasi-linear formulation of MOND (Milgrom 2010) in which the modified gravitational potential |$\phi$| is given by the solution of the equation,

where |$\phi _N$| is a Newtonian potential given by the standard Poisson equation |$\nabla ^2 \phi _N = 4 \pi G \rho$|⁠, where |$\rho$| denotes the baryonic mass density. The above equation is linear in |$\phi$| and can be solved in the same way as the standard equation.

The transition between the regime of high acceleration (Newtonian) and the regime of low acceleration (Milgromian) is described by a function |$\nu$| whose particular form is not given by the theory. There are several different classes of functions (see e.g. Famaey & McGaugh 2012). In our analysis, the following function is used:

where |$\tilde{y} \equiv |\nabla \phi _N|/a_0$| and |$a_0$| is a fundamental acceleration of the order of |$\sim 10^{-10}\,$| m s⁻². Both the value of |$a_0$| and the coefficient |$\alpha$| must be determined from observations. We use |$a_0 = 1.2 \times 10^{-10}\,$| m s⁻² (McGaugh, Lelli & Schombert 2016) as well as different values of |$\alpha = 1.5, 2$| or 2.5. The coefficient |$\alpha$| determines how smooth the transition between the Newtonian and Milgromian regimes is. For smaller |$\alpha$|⁠, the transition is smoother.

Blanchet & Novak (2011) constrained the MOND transition function by analysing the solar system data and for the class of functions defined in equation (2) they found that |$\alpha \ge 8$| in order to fit the precession rates of the planets. Similarly, Hees et al. (2016) constrained the parameter to be |$\gt 6$|⁠. For the ‘standard’ value of |$\alpha = 2$|⁠, the expected EFE-induced precession rates make only |$0.1-5\,$|ppm (it is larger for the outer planets) of the precession rates caused by the Newtonian planet-planet perturbations. The parameters of the solar system are, however, known with significantly worse precision. The relative standard deviation of the planets’ masses |$\sim 46\,$|ppm (slightly more for Mercury, |$\sim 64\,$|ppm).² The EFE-induced precession can be then balanced by changing the planets’ masses within the uncertainties. Due to the difficulties in determining the |$\alpha$| parameter, we use the ‘standard’ value of 2 and two other close values, 1.5 and 2.5.

Although MOND was formulated to model the rotation curves of galaxies, its consequences extend to smaller scales. When a massive body is embedded in an external gravitational field, e.g. the Sun in the Galaxy’s field, the so-called EFE appears (Milgrom 2009). It stems from the fact that between the Sun and the Galaxy centre there is a point where the resulting Newtonian gravitational acceleration disappears and in a certain area around it the Newtonian acceleration |$g_N \equiv |\nabla \phi _N| \lesssim a_0$|⁠. Solving equation (1) for such a system, it turns out that the gravitational field is modified not only in the region where |$g_N \lesssim a_0$|⁠, but also near the Sun where |$g_N \gg a_0$|⁠. The perturbing acceleration has the form of a quadrupole field in the vicinity of the Sun (Milgrom 2009).

In this work we solve equation (1) without the assumption of a constant external field and beyond the limit |$g_N \gg a_0$|⁠, which would allow us to follow the evolution of trans-Neptunian objects thousands of astronomical units away from the Sun. We treat the field of the Galaxy as a point-source field. There are numerous estimates for the distance and rotational velocity of the Sun with respect to the centre of the Galaxy. We use one of the most recent estimates from (Hunt, Bovy & Carlberg 2016), in which the velocity |$\mathrm{ V}_{\odot } = 239 \pm 9\,$|km s|$^{-1}$| and the distance |$\mathrm{ R}_{\odot } = 7.9 \pm 0.3\,$| kpc.

The rotational velocity and distance can be translated into the mass of the Galaxy inside |$\mathrm{ R}_{\odot }$|⁠, denoted |$M_{{\mbox{G}}}$|⁠. Since the centripetal acceleration of the Sun |$a_{{\mbox{cp}}} \equiv \mathrm{ V}_{\odot }^2/\mathrm{ R}_{\odot } \approx 2.3 \times 10^{-10}\,$| m s⁻², which is only |$\approx 2 a_0$|⁠, we need MOND prescription to find |$M_{{\mbox{G}}}$|⁠, thus

where |$\mu (\tilde{x})$| is a transition function of MOND³ where |$\tilde{x} \equiv g/a_0$| (Milgrom 1983a) of the form

The value of |$\mu$| depends on |$\alpha$| and for |$\tilde{x} \approx 2$| and |$\alpha = 2$| this gives |$\mu \approx 0.89$|⁠, resulting in |$M_{{\mbox{G}}}\approx 9.34 \times 10^{10}\, \mbox{m}_{\odot }$|⁠, while for other |$\alpha$| the mass may be slightly larger or smaller. Naturally, different |$\mathrm{ R}_{\odot }$| and |$\mathrm{ V}_{\odot }$| would give different |$M_{{\mbox{G}}}$|⁠. We have tested various values of these quantities and found that they are not crucial for the overall dynamics of trans-Neptunian objects, although, the evolution of a particular object would change in detail for different |$\mathrm{ R}_{\odot }$| and |$\mathrm{ V}_{\odot }$|⁠.

Another aspect of the Galaxy’s gravitational field must be considered. The Sun does not move in an exactly circular and planar orbit. Its velocity has non-zero components both in the radial direction and in the direction perpendicular to the Galaxy disc. An analysis of the periodicity of the comet flux of the Oort cloud (Matese et al. 1995) shows that the period of the radial motion |$T_R \sim 170\,$| Myr and an amplitude of the radial oscillations |$\sim 0.5\,$| kpc. The variation of |$g_N$| with an amplitude of |$\sim 14\,$| per cent is obtained. The vertical oscillation occurs with a period of |$70\,$| Myr. Since the current perpendicular velocity |$V_z \sim 7\,$|km s|$^{-1}$| (Karim & Mamajek 2017), the maximum deviation of the Sun above or below the midplane of the Galaxy disc, |$z_{{\mbox{max}}} \sim 100\,$| pc. Assuming harmonic oscillations, the maximum acceleration perpendicular to the disc is |$g_{z,{\mbox{max}}} \sim 2.5 \times 10^{-11}\,$| m s⁻², which is an order of magnitude weaker than the centripetal acceleration. The maximum angle by which the total acceleration deviates from the radial direction is |$\sim 7\,$|°, while in the current epoch it is |$\sim 1\,$|°.

A more realistic model of the Galaxy should also take into account the dependence of |$M_{{\mbox{G}}}$| on |$\mathrm{ R}_{\odot }$|⁠, which changes with time. Because of the uncertainties in determining |$\mathrm{ R}_{\odot }$| and |$\mathrm{ V}_{\odot }$|⁠, we omit all these corrections from our analysis and treat the Galaxy’s potential as a point-mass potential and the Sun as moving in a circular orbit around the centre of the Galaxy.

If the masses |$\mbox{m}_{\odot }$| and |$M_{{\mbox{G}}}$| are known, the total potential at the position |${\pmb {r}}$| is given by

where |${\pmb {r}}_{\odot }$| and |${\pmb {r}}_{{\mbox{G}}}$| are the position vectors of the Sun and the Galactic centre. It can be used to calculate the right-hand side of equation (1). Since |$\phi _N$| has axial symmetry, the potential |$\phi$| we are looking for has the same symmetry. We can therefore use a cylindrical coordinate system whose origin is in the Sun. One of the axes, |$z$|⁠, points towards the centre of the Galaxy; the second, |$x$|⁠, is perpendicular to |$z$|⁠. We solve equation (1) numerically using the standard five-point method (e.g. Hoffman 2001). The grid is non-uniform, the size of the grid cells ranging from |$\sim 2.6 \ \mbox{au}$| near the Sun and the critical point (defined by |$g_N = 0$|⁠) to |$\sim 670 \ \mbox{au}$| in the outer parts of the domain. The size of the area is |$x \in [0, 90]\,$| kau and |$z \in [-90, +90]\,$| kau. The boundary condition at |$x = 0$| results from the cylindrical symmetry, i.e.

while at the other three limits the potential corresponds to the Milgromian potential of the point-mass Galaxy alone.

Fig. 1 shows the solution of equation (1) for the transition function coefficient |$\alpha =2$|⁠. The so-called phantom mass density is shown in the upper left-hand panel. With respect to the dark matter hypothesis, the modification of the standard Poisson equation, equation (1), can be interpreted as the addition of dark matter to the right-hand side of the equation, i.e.

where |$\rho _p$| would be the density of dark matter. In terms of MOND, the additional mass is called phantom dark matter, whose distribution is determined by the baryonic mass distribution according to equation (4) (Milgrom 2010).

The phantom mass density can be positive or negative, while it disappears at the position of the Sun (which agrees with the analytical results in Milgrom 2009), as well as at the axially symmetric surface marked with green dashed curves. The remaining three panels of Fig. 1 show the perturbative acceleration |${\pmb {a}}_{{\mbox{pert}}}$| given by the partial derivatives of |$\phi$| (computed numerically), completed with the centrifugal acceleration. The latter results from the fact that the reference frame is not inertial as it relates to the Sun orbiting the centre of the Galaxy. The addition of the centrifugal acceleration is equivalent to the substruction of the Milgromian gravitational acceleration of the Sun by the Galaxy. The total acceleration of a test particle with respect to the Sun is then given by

where the components of |${\pmb {a}}_{{\mbox{pert}}}= (a_x, a_y, a_z)$| read⁴

The acceleration |${\pmb {a}}_{{\mbox{pert}}}$| disappears at the Sun and at two other points on the |$z$|-axis. It is strongest in the lobes extending around the critical point, where its magnitude reaches |$3 \times 10^{-11}\,$|m s⁻². The lower left-hand panelshows the magnitude of the acceleration relative to the Newtonian acceleration. In the region of |$\sim 5000 \ \mbox{au}$| the perturbation can be up to |$10~$| per cent of |$g_N$|⁠, which means that the elongated orbits of the semimajor axes |$\sim 2500 \ \mbox{au}$| are strongly perturbed. For even wider orbits, |$a \gtrsim 10\,$| kau, the relative perturbation can be as high as |$20~$| per cent (this region is not shown). The dynamics of the Oort cloud is thus supposed to be significantly affected by EFE (Iorio 2010; Paučo & Klačka 2017).

The bottom right-hand panelshows the perturbative acceleration in the smaller region around the Sun. The vector field shows that the perturbing force is acting along the |$x$|-axis towards the Sun and along the |$z$|-axis away from it. It has the form of a quadrupole field, as it should be for EFE near the Sun (Milgrom 2009). However, in this plot both EFE and the enhanced gravity effect (EGE; the effect of MOND for an isolated Sun) are shown together. The latter effect is radially symmetric and acts in the direction of the Sun.

For different values of |$\alpha$|⁠, the perturbative acceleration differs both near the Sun and farther from it. Fig. 2 shows the results for |$\alpha = 1.5$| (top row) and |$\alpha = 2.5$| (bottom row) at two different scales. For |$\alpha = 1.5$| the perturbation is generally stronger, but near the Sun it acts as an additional gravitational pull regardless of direction. This means that the EGE dominates the EFE. For |$\alpha = 2.5$| the perturbation is weaker and the picture is more complex. Similar to the case of |$\alpha = 2$|⁠, EFE is more important than EGE.

Near the Sun, EFE has the same functional form as tidal acceleration at position |${\pmb {r}}$| by a distant perturber, i.e.

where |$A_{{\mbox{tidal}}}$| is the tidal parameter that can be expressed with the mass |$m_{{\mbox{pert}}}$| and the distance |$\Delta _{{\mbox{pert}}}$| of the perturber as |$A_{{\mbox{tidal}}}= m_{{\mbox{pert}}}/\Delta _{{\mbox{pert}}}^3$|⁠. The unit vector pointing from the Sun to the perturber is denoted by |$\hat{{\pmb {r}}}_{{\mbox{pert}}}$|⁠. Within the EFE, |$\hat{{\pmb {r}}}_{{\mbox{pert}}}$| would point towards the centre or anticentre of the Galaxy. Using the numerically solved modified Poisson equation, for |$\alpha = 2$| we can find the tidal parameter |$A_{{\mbox{tidal}}}= 1.77 \times 10^{-13}\, \mbox{m}_{\odot } \, \mbox{au}^{-3}$|⁠, while for |$\alpha = 2.5$|⁠, |$A_{{\mbox{tidal}}}= 1.30 \times 10^{-13}\, \mbox{m}_{\odot } \, \mbox{au}^{-3}$|⁠. As mentioned earlier, for |$\alpha = 1.5$| the EGE dominates the EFE, hence the total perturbation deviates from the tidal form.

Holman & Payne (2016a) analysed the astrometry of Pluto and a few TNOs to constrain the tidal parameter due to the hypothetical Planet Nine and found that the upper limit of |$A_{{\mbox{tidal}}}$| is between |$10^{-12}$| and |$10^{-10}\, \mbox{m}_{\odot } \, \mbox{au}^{-3}$|⁠, depending on the position of the perturber in the sky. Below this limit, the perturbation cannot be detected based on the currently available observations. Both values of the EFE–induced |$A_{{\mbox{tidal}}}$| are well below the upper limit.

The values can be expressed in terms of |$m_{{\mbox{pert}}}$| and |$\Delta _{{\mbox{pert}}}$|⁠. Taking a representative mass of Planet Nine, |$m_{{\mbox{pert}}}= 10\, \mbox{m}_{\oplus }$| (Batygin & Brown 2016a), we find that the values of |$A_{{\mbox{tidal}}}= 1.77 \times 10^{-13}\, \mbox{m}_{\odot } \, \mbox{au}^{-3}$| and |$A_{{\mbox{tidal}}}= 1.30 \times 10^{-13}\, \mbox{m}_{\odot } \, \mbox{au}^{-3}$| can be obtained for |$\Delta _{{\mbox{pert}}}= 554 \ \mbox{au}$| and |$\Delta _{{\mbox{pert}}}= 614 \ \mbox{au}$|⁠, respectively. For smaller |$m_{{\mbox{pert}}}= 6.2\, \mbox{m}_{\oplus }$| (Brown & Batygin 2019) the distances are |$472 \ \mbox{au}$| and |$524 \ \mbox{au}$|⁠, respectively. All distances are consistent with the current constraints for Planet Nine’s orbit. According to Brown & Batygin (2019), the most likely position of P9 in the sky is near the anticentre of the Galaxy and the most likely distance is |$\sim 400-500 \ \mbox{au}$|⁠.

While the long-term dynamics of TNOs within the two hypotheses, i.e. Planet Nine and MOND, differ from each other, when modelling an additional unexplained acceleration in the outer parts of the solar system, Milgrom’s gravity corrections cannot be distinguished from the acceleration due to Planet Nine. Fienga et al. (2016) analysed the Cassini radio ranging data to constrain the true anomaly of Planet Nine, assuming the orbital elements according to Batygin & Brown (2016a). They found that for the true anomaly |$\nu \in [108, 129]\,$|° the Cassini residuals are reduced. This range of |$\nu$| corresponds to the direction being shifted by |$50-60\,$|° with respect to the direction towards the centre of the Galaxy. The latter lies within the uncertainty zone for which the perturbative acceleration is below the detection limit. Iorio (2017) analysed the perturbation of the Saturn’s orbit inferred from the Cassini data and found that the true anomaly of Planet Nine could be constrained to |$\nu \in [130, 240]\,$|°, i.e. close to the aphelion, which itself lies close to the anticentre of the Galaxy. Holman & Payne (2016b) has relaxed constraints on the mass and orbit of Planet Nine in their analysis of the Cassini data. It turns out that the direction towards the centre of the Galaxy is on the boundary between decreasing and increasing residuals. Therefore, the MOND perturbation of gravity cannot be ruled out by the Cassini data.

3 A SECULAR MODEL OF ISOLATED AND FIXED EFE

Once the modified Poisson equation is solved for given parameters of the Galaxy, the perturbative acceleration for any |$(z,x)$| point can be determined by interpolation between the grid nodes. Before turning to the |$N$|-body integrations of an asteroid under the combined effect of Newtonian acceleration due to the Sun, the Galaxy and the giant planets, as well as the MOND perturbation, in this section we examine the secular dynamics of an asteroid without the gravitational attraction of the planets and assuming that the position of the Galaxy centre, which generally varies in time due to the motion of the Sun around the Galaxy centre, is fixed in the ecliptic reference frame.

The Hamiltonian of a test particle under the MOND effect of the Sun and the Galaxy is as follows:

where |$\mathcal {H}_1$| is the Keplerian part, while |$\mathcal {H}_2$| is the perturbative Hamiltonian. The perturbative Hamiltonian is equal to

where |$\phi$| is the solution of equation (1) and |$\phi _{{\mbox{cf}}}$| is the centrifugal potential, which is as follows:

The constant value is added to the Hamiltonian in equation (7) so that |$\mathcal {H}_2 = 0$| at the Sun.

The canonical Delaunay angle-action variables are defined as (e.g. Brouwer & Clemence 1961)

where |$\mathcal {M}, \omega , \Omega$| are the mean anomaly, the argument of the perihelion and the longitude of the ascending node, respectively, while |$a, e, i$| denote the semimajor axis, the eccentricity, and the inclination of the orbit. The Keplerian Hamiltonian depends only on |$I_1$|⁠, i.e.

The perturbing Hamiltonian generally depends on all variables. However, because of the axial symmetry of the perturbation, it is useful to choose the |$z$|-axis of the reference frame to coincide with the axis of symmetry (also denoted |$z$|⁠). With such a choice, |$\mathcal {H}_2$| does not depend on |$\varphi _3$|⁠, so |$I_3$|⁠, the projection of angular momentum onto the |$z$|-axis, is the integral of motion. The Hamiltonian therefore has two degrees of freedom. To further simplify the Hamiltonian, the method of averaging (e.g. Arnold, Kozlov & Neishtadt 2006) is used. The perturbation |$\mathcal {H}_2$| is much smaller than the Keplerian part, so we can average the Hamiltonian over the mean anomaly, which is a fast variable

Averaging over the mean anomaly can be replaced by averaging over the eccentric anomaly |$E$|⁠, using the Kepler equation. We obtain

After the averaging, |$\langle \mathcal {H} \rangle$| does not depend on |$\varphi _1$|⁠, thus |$I_1$| is an integral of motion. The secular Hamiltonian of a test particle thus has one degree of freedom and is integrable. The averaging is done numerically using the 100th order Gauss-Legendre quadrature. For a given position of a test particle in its Keplerian orbit, the value of |$\mathcal {H}_2$| (⁠|$\mathcal {H}_1$| does not need to be averaged as it is constant after averaging) is determined using bicubic interpolation on the grid nodes.

Once the integrals |$I_1$| and |$I_3$| are set, the secular Hamiltonian depends on |$(\varphi _2, I_2)$|⁠, which translates into the dependence on the argument of pericentre and eccentricity. Since the reference frame whose |$z$|-axis points towards the centre of the Galaxy (the G-frame from now on)⁵ differs from the ecliptic reference frame normally used to describe the orbits of planets and asteroids in the solar system, the argument of pericentre and the inclination are referred to as |$\omega _G$| and |$i_G$|⁠, respectively, when expressed in the reference frame related to the Galaxy.

Since the secular Hamiltonian has one degree of freedom, it is sufficient to construct a phase diagram using contours with constant values of |$\langle \mathcal {H} \rangle$|⁠. Fig. 3 shows such contours for different |$I_1$| and |$I_3$|⁠, which can be translated into |$a$| and |$i_0$|⁠, the latter being the inclination of a circular orbit, which gives |$I_3 = I_1 \cos i_0$|⁠. Each panel was obtained for different |$\alpha$|⁠. The plots show representative phase diagrams. The dynamics of the system from the left-hand panel is simple. There are oscillations of |$e$| accompanied by rotations of |$\omega _G$|⁠. There is a stable equilibrium near the origin of the diagram. In the middle panel the dynamics is more complex as there are two stable equilibria and one unstable equilibrium. The diagram is asymmetric, i.e. the equilibrium for |$\omega _G = \pi /2$| occurs for lower |$e$| than that for |$\omega _G = -\pi /2$|⁠. In the right-hand panel, the situation is reversed. Due to the axial symmetry of the problem, the equilibria only exist for |$\omega _G = \pm \pi /2$|⁠.

The equilibria are examined more systematically in Fig. 4. Each panel corresponds to a particular pair |$(a, \alpha)$| and shows families of equilibria parametrized by |$i_0 \in [0, \pi ]$|⁠. Since |$\omega _G = \pm \pi /2$|⁠, the ordinate axis is |$e \sin \omega _G$| and a positive value means that |$\omega _G = \pi /2$|⁠, while a negative value means that |$\omega _G = -\pi /2$|⁠.

For |$\alpha = 1.5$| there are only stable equilibria (apart from a tiny part of the rightmost panel, i.e. |$a = 5000 \ \mbox{au}$| and |$e \sin \omega _G \approx -1, i_G \approx \pi /2$|⁠), the one for |$e \sim 0$| (increasing for larger |$a$|⁠) in the whole range of |$i_G$|⁶ and for |$i_G \sim 0$| or |$\pi$| (or |$\sin i_G \sim 0$|⁠) in the whole range of |$e$|⁠. The first family is characterized by |$\omega _G = \pi /2$|⁠, while the other families have |$\omega _G = -\pi /2$|⁠. The families with |$\sin i_G \sim 0$| are not seen in the diagrams in Fig. 3.

For |$\alpha = 2$| and |$\alpha = 2.5$| there is a family of unstable equilibria. They exist for sufficiently high inclinations, e.g. for |$a = 500 \ \mbox{au}$| the unstable equilibria appear for |$i_G \gtrsim 40-50~$|° depending on |$\alpha$|⁠. They appear when the stable equilibrium |$e \sim 0$| bifurcates. As |$i_0$| (and |$i_G$|⁠) increases, the stable equilibria move towards higher eccentricities, reaching 1 for |$i_G = \pi /2$|⁠. For larger |$a$|⁠, the families of equilibria change. The unstable equilibria exist for a narrower range of |$i_G$| or may even disappear.

Since the centre of the Galaxy is near the ecliptic, the orbits with low inclination |$i$| (in the ecliptic reference frame) have a high inclination |$i_G$| in the reference frame defined by the centre of the Galaxy. Furthermore, the orbits whose perihelions point in the direction of the Galaxy centre or anticentre have |$\omega _G \sim \pi /2$| and |$-\pi /2$|⁠, respectively.

The presence of an unstable equilibrium for |$e \sim 0$| leads to high-amplitude oscillations of |$e$| between |$\sim 0$| and |$\sim 1$|⁠, which means that an asteroid with a large |$a$| can have its perihelion during evolution both close to the Sun and very far from it.

In order to follow the evolution of the mean system, one can use one of two approaches. In the first case, the secular Hamiltonian, equation (8), is used to construct the canonical equations of motion. In this case, one needs to calculate the partial derivatives of the Hamiltonian numerically. In the second approach, the Gauss’s planetary equations (e.g. Murray & Dermott 1999) are averaged over the Keplerian motion. We have chosen the second method, although both are good.

An example evolution is shown in Fig. 5. The perihelion distance (the upper left-hand panel) varies between |$\sim 1000 \ \mbox{au}$| and |$\lesssim 5 \ \mbox{au}$|⁠. The asteroid remains detectable over a short part of the periodic oscillations of |$q$| (⁠|$q \lesssim 80 \ \mbox{au}$| for the current statistics of known TNOs), in particular the asteroid can be classified as a Centaur (the grey area denotes |$q \in [5.2, 30] \ \mbox{au}$|⁠) or as an ETNO (the yellow area denotes |$q \in [30, 80] \ \mbox{au}$|⁠).

The upper middle panel of Fig. 5 shows the evolution of the inclination in two different reference frames, i.e. the red curve corresponds to the G-frame, while the black curve corresponds to the frame in which the coordinates of the Galaxy centre are |$(0, -1, 0)$|⁠. The latter are close to the current position of the Galaxy centre in the ecliptic frame and can be referred to as the E-frame for short. The inclination in the G frame, |$i_G$|⁠, varies between |$\sim 90$|° and |$\sim 60~$|°, while the inclination in the E frame, |$i_E$|⁠, varies throughout the whole range. The variation is faster when |$q$| is small. The orbit is alternately prograde and retrograde. The rapid reorientation of the orbital plane results from the precession of the orbit in the G-frame (the red curve in the upper right-hand panel). For a given torque acting on the orbit, the precession is faster with a lower angular momentum of the orbit, i.e. higher |$e$|/lower |$q$| when |$a$| is constant. In the E-frame, the longitude of the ascending node oscillates.

The bottom left-hand panel of Fig. 5 shows the evolution in the |$(i_E, q)$| diagram. The initial position is marked with the blue symbol. The system evolves along a U-shaped trajectory, visiting the regions with |$i_E$| either |$\sim 0$|° or |$\sim 180~$|° for moderate and large |$q$|⁠, and with |$i_E$| between these values for small |$q$|⁠. The evolution in this diagram is not periodic since there are two fundamental frequencies of motion. One is related to the variation of |$\omega$|⁠, the second to |$\Omega$|⁠. While |$i_G$| and |$q$| are related by the conservation of projection of angular momentum on the |$z$|-axis, this relationship does not exist in the E-frame. As a result, the phase trajectory is not closed and with a longer integration (more cycles of the |$q$| oscillations), the trajectory would fill the diagram more evenly and the U-shaped structure might be less clear.

The two remaining panels of Fig. 5 show the evolution of the argument of pericentre |$\omega _G$| and |$e$| (bottom middle) and |$q$| (bottom right). The direction of the motion is marked with arrows. The phase trajectory is closed since |$\Omega _G$| is a cyclic variable in the G-frame. The eccentricity varies between low and high values along an 8-shaped trajectory. The variation of the perihelion distance is shown in the logarithmic scale. As was mentioned earlier, the maximum |$e$|/minimum |$q$| is reached for |$\omega _G = \pm \pi /2$|⁠. The example system has a large initial |$q$|⁠, but we could also start the simulation with a small |$q$|⁠. If the initial parameters are |$q \sim 5 \ \mbox{au}$| (which would correspond to an asteroid scattered by Jupiter) and |$\omega _G \sim 90~$|°, the value of |$\omega _G$| reached for |$q \in [30, 80] \ \mbox{au}$| (the upper limit corresponds to the current detection limit) would be close to |$\sim 45~$|°. Other values of |$\sim -45, -135, 135~$|° are reached over the whole cycle. This would mean that the apsidal line could be shifted by about |$45~$|° with respect to the Galaxy centre-anticentre line.

The presence of the unstable equilibrium at low eccentricity leads to a large-amplitude variation of |$q$|⁠. For |$\alpha = 1.5$| there is no such equilibrium and the variation occurs with smaller amplitude. However, it is still possible to drift perihelia of scattered asteroids outside the orbit of Neptune. Fig. 6 shows the energy levels for |$\alpha = 1.5$| and |$\alpha = 2$|⁠, which is for a moderate semimajor axis |$a = 500 \ \mbox{au}$| and |$i_0 = 90\,$|° (implying that |$I_3 = 0)$|⁠. For both values of |$\alpha$|⁠, the perihelia, which are initially close to Jupiter’s orbit (⁠|$q \sim 5.2 \ \mbox{au}$|⁠), are driven to greater distances. Therefore, it can be difficult to establish a specific value of |$\alpha$| based on the observed properties of ETNOs. Most of the simulations presented further in this paper were performed for |$\alpha = 2$|⁠.

In a more realistic model, the EFE-induced cycle competes with planetary perturbations. While the former speeds up with larger orbits, the latter depends inversely on |$a$|⁠. Fig. 7 illustrates these dependencies. The time-scales for |$\omega$| and |$\Omega$| are shown as functions of |$a$| with the red and green curves, respectively. The solid curves denote the EFE, while the dashed curves correspond to the planetary perturbation. The time-scales depend not only on |$a$|⁠, but also on |$e$| and |$i$|⁠, as well as |$\omega _G$| for the EFE perturbation. Therefore, for illustration, we have chosen |$e = 0.8$| and |$i_G = 70~$|° (which corresponds to |$i_E = 20~$|°). These numbers correspond to typical values for known ETNOs. The initial values for |$\omega _G$| and |$\Omega _G$| are set to 0.

The orbits of the planetary perturbers are assumed to be circular and coplanar with each other. The perturbations are expanded to third order in a Taylor series with respect to a small quantity |$a/a_P$|⁠, where |$a_P$| is the size of planets’ orbits, and then averaged over the Keplerian motions. As is known, planetary perturbations lead to the precession of an asteroid orbit with respect to the axis perpendicular to the orbital plane of the planets, as well as to the rotation of the apsidal line.

It is clear that for |$a \sim 250 \ \mbox{au}$| the time-scales for the two effects are approximately equal. For smaller |$a$| the planetary perturbations dominate, while for larger |$a$| the MOND effect is more significant. For the semimajor axis beyond a few hundred astronomical units, the MOND-induced evolution dominates over the planetary perturbations by several orders of magnitude. Moreover, for |$a \gtrsim 350 \ \mbox{au}$|⁠, the Milgromian evolution occurs over a period shorter than the Sun’s orbital period around the centre of the Galaxy, which is |$\sim 200\,$| Myr according to Hunt et al. (2016). For |$a \gtrsim 3000 \ \mbox{au}$|⁠, the EFE-induced cycle is only an order of magnitude slower than Keplerian motion, so the secular model is not sufficient.

In this paper we extend the model in two ways. We abandon averaging because we want to model possible close approaches between the planets and an asteroid. Thus, the interaction between planets and asteroids is modelled within the restricted |$N$|-body model. We also take into account the motion of the Sun in the Galaxy (a circular orbit is assumed), which causes the axis of symmetry of EFE to change in time. Before that, in the next section, we describe observational evidence for the hypothesis of an unseen planet in the periphery of the solar system.

4 EXTREME TRANS-NEPTUNIAN OBJECTS AND THE PLANET NINE HYPOTHESIS

The hypothesis of Planet Nine is justified by the anomalous structure of trans-Neptunian solar system. The first indication of the planet’s existence was the detachment of the perihelia of Sedna and several other ETNOs from the gravitational influence of the known planets (Brown et al. 2004; Gomes, Matese & Lissauer 2006). Nevertheless, the existence of Planet Nine would be problematic if the detachment of the perihelia were the only evidence for it, since the proposed orbit of Planet Nine is itself detached from the giant planets (Batygin & Brown 2016a) and could not be explained by planet–planet scattering alone. In order to solve this problem, three mechanisms are proposed: (i) in situ formation in an extended protoplanetary disc (Kenyon & Bromley 2016); (ii) scattering from the region of giant planets followed by orbital circularization (Bromley & Kenyon 2016; Eriksson, Mustill & Johansen 2018); (iii) capture of a free-floating planet or a planet from another system during close stellar encounters (Li & Adams 2016).

Each of these scenarios has its difficulties (see Batygin et al. 2019 for a brief overview on this topic). In the second scenario, for example, it was proposed that the planet was first scattered and then its orbit was circularized by planetesimals forming an extended disc up to |$\sim 1000 \ \mbox{au}$| (Eriksson et al. 2018). The chain of proposed solutions to the ETNOs detachment problem is then as follows: to explain detached ETNOs – a distant detached Planet Nine is proposed; to explain detached Planet Nine – an extended detached disc of planetesimals is proposed. The addition of more and more objects only seems to shift the problem, and without more evidence, the Planet Nine hypothesis would be in trouble.

However, there are other features of ETNOs that support this hypothesis. Planet Nine has been shown to explain the high inclinations of Centaurs (Batygin & Brown 2016b), which cannot otherwise be reproduced due to scattering by known planets. However, the most telling feature of ETNO orbits is the clustering of orbital planes and the apsidal confinement reported for |$a \gtrsim 250 \ \mbox{au}$| (Batygin & Brown 2016a). As mentioned in the introduction, the robustness of these features has been questioned by several authors (Bernardinelli et al. 2020; Clement & Kaib 2020; Napier et al. 2021). In order to verify the features, we have presented the statistics of the TNO orbital elements in Fig. 8. In the upper panel, the longitudes of the ascending node are plotted against the semimajor axes. The colours of the symbols as well as their size encode the perihelion distance (see caption). The grouping of |$\Omega$| is not clear, although for |$a \gtrsim 200 \ \mbox{au}$| it can be noted that the range |$(0, 180)\,$|° is preferred over the range |$(-180, 0)\,$|°, especially for the most detached orbits (magenta symbols with |$q \gt 45 \ \mbox{au}$|⁠). The distribution of |$\varpi$| (middle panel) shows a similar feature. The grouping of |$\omega$|⁠, on the other hand, is more likely to be found in the range |$(-90, 90)\,$|°.

By choosing certain criteria for |$a$| and |$q$|⁠, one can achieve a more or less clear grouping of the elements. Brown & Batygin (2021) suggested to consider only objects with |$q \gt 42 \ \mbox{au}$| and |$a \in (150, 1000) \ \mbox{au}$| which are on the one hand sufficiently detached from Neptune and on the other hand not too far away to be influenced by Planet Nine. Fig. 9 illustrates the clustering of orbital planes in the |$(i \cos \Omega , i \sin \Omega)$|-diagram for different selection criteria. Each row represents a different minimum |$a$|⁠, while each column represents a different minimum |$q$|⁠. The upper half-plane of each panel is favoured over the lower half-plane for all criteria, which means that the clustering of orbital planes is independent of a particular choice of |$a$| and |$q$| limits.

The apsidal confinement for prograde ETNOs is shown in Fig. 10 for different selection criteria. The perihelia group around the ecliptic, which can be attributed to observational bias. For |$q \gt 40 \ \mbox{au}$|⁠, the ecliptic longitude in the range |$(0, 180)\,$|° is preferred over the range |$(180, 360)\,$|°. This corresponds to a grouping around the direction of the anticentre of the Galaxy, although there are also perihelia closer to the centre of the Galaxy.

The statistics of the perihelia and inclinations are shown in Fig. 11. Apart from the previously selected objects, this set is extended by Centaurs, |$q \in (5.2, 30) \ \mbox{au}$| and objects with smaller |$a$| down to |$100 \ \mbox{au}$|⁠. Centaur statistics are limited to |$a \lt 2000 \ \mbox{au}$|⁠. The inclinations of the Centaurs are distributed over the whole range, but for |$i \gtrsim 45\,$|° only objects with |$q \lesssim 10 \ \mbox{au}$| are observed. Moreover, for the objects with |$a \gt 500 \ \mbox{au}$|⁠, the perihelion distances are close to |$5.2 \ \mbox{au}$|⁠, indicating that these Centaurs were scattered by Jupiter. The objects with |$q \gt 30 \ \mbox{au}$| have inclinations of |$\lesssim 60\,$|°, but there are three recently discovered objects with |$q \gt 50 \ \mbox{au}$| and retrograde orbits (they are labelled in the graph). Also, we can see that the objects with |$q \gt 30 \ \mbox{au}$| and |$a \gt 500 \ \mbox{au}$| (magenta dots) generally have lower inclinations than the objects with |$a \lt 500 \ \mbox{au}$|⁠.

As shown in the previous section and already demonstrated by Paučo & Klačka (2016), the detachment of the perihelia as well as highly inclined Centaurs can be naturally explained within EFE. The connection between Centaurs and ETNOs is also natural in EFE and was also shown by Paučo & Klačka (2016). The lower left-hand panel in Fig. 5 shows that they belong to the same population of objects, but at two different stages of the evolution. We can also see that the retrograde ETNOs can be explained within EFE as well.

The apsidal confinement as well as the clustering of the orbital planes are more difficult to explain. Paučo (2017) concluded that the EFE quadrupole strength parameter required to reconstruct the observed properties of the orbits is not consistent with observational constraints.⁷ As we mentioned in the introduction, a comparison between the synthetic and the observed orbits is problematic. The difficulty lies in the choice of the initial set of orbits and in the way the final set is constructed for comparison with the observations. To understand the first problem, in the next section we trace the evolution of the observed ETNOs backwards in time and show that they originate from the region of the giant planets of the solar system.

5 EVOLUTION OF KNOWN ETNOS AND CENTAURS

The restricted |$N$|-body model with MOND perturbation is used to study the evolution of the known ETNOs and Centaurs. The model of motion is the following. The system consists of the Sun, the four giant planets and a given number of massless asteroids. The model is thus restricted in the sense that the asteroids do not affect the motions of the planets. Apart from the Newtonian gravitational interaction between the bodies, each of them is subject to the perturbation |${\pmb {a}}_{{\mbox{pert}}}$|⁠, which is determined for each position in space using bicubic interpolation, as described earlier in this paper.

The perturbation is computed on the assumption that there are no massive bodies other than the Sun and the point-like Galaxy, which is not necessarily fulfilled due to the existence of massive planets. However, a more general Milgromian model of motion is much more complex and would be impractical for the problem studied here. The simplified model can be validated in the following way. In the regime of giant planets, the planetary contribution to |$\phi _N$| is not negligible, but the MOND perturbation, |${\pmb {a}}_{{\mbox{pert}}}$|⁠, is very weak. On the other hand, MOND effects are important for large distances from the Sun, where |$\phi _N$| used to compute |${\pmb {a}}_{{\mbox{pert}}}$| can be well approximated by the Sun alone, the planetary contribution to |$\phi _N$| is of the order of |$10^{-3}$|⁠. In both ranges, an additional contribution to |${\pmb {a}}_{{\mbox{pert}}}$| resulting from the presence of giant planets is a higher order effect and is omitted in this analysis.

5.1 Prograde ETNOs

Following the original selection criteria in Batygin & Brown (2016a), we selected 25 ETNOs with |$a \gt 250 \ \mbox{au}$| and |$q \gt 30 \ \mbox{au}$| and randomly chose 100 clones for each of them according to their orbital uncertainties. Such a set of 2525 massless particles were numerically integrated backwards in time. All the asteroids evolve into orbits with perihelions well inside the orbit of Neptune. For ETNOs with initial |$a \gtrsim 1000 \ \mbox{au}$| a few Myr were enough to reach the low |$q$| range, while for the ones with |$a \sim 250 \ \mbox{au}$| a few hundreds Myr was necessary.

Once a given object reaches the giant planets region, its evolution becomes strongly chaotic due to close encounters with the planets. Close encounters result in the energy and angular momentum gain or loss of the object, thus not only |$e$| but also |$a$| is modified. If the energy is lost, the asteroid orbit becomes smaller and may reach the region |$a \lt 30 \ \mbox{au}$|⁠. If also the eccentricity decreases, the entire orbit of the object can reside in the giant planets region.

The evolution of a particular object cannot be traced because of the strong chaos and one should only investigate the past orbits of ETNOs statistically. Nevertheless, to illustrate the possible origin of a particular ETNO, we present the evolution of Sedna (one of its clones) in Fig. 12. Once its perihelion reaches the orbits of Neptune, Uranus, Saturn, and Jupiter (upper left-hand panel), the semimajor axis changes chaotically towards smaller values down to the sizes of the orbits of Saturn and Jupiter (upper right-hand panel). This occurs in the epoch around |$-30\,$| Myr and is accompanied by a decrease in eccentricity (the middle left-hand panel), while the inclination remains moderate (the middle right-hand panel) and the |$\Omega$| and |$\omega$| angles vary throughout (the lower panels).

Further past evolution results in the ejection of the object. However, this does not mean that this object entered the solar system from outside, if the evolution of the clone was the real evolution of Sedna. It is the result of a strong chaos and sooner or later every object is ejected from the solar system, regardless of the direction of evolution.

Fig. 13 shows the temporary evolution of the same clone of Sedna about the epoch |$-30\,$| Myr. The clone remains close to the orbit of Jupiter or Saturn and is temporarily in a 1:1 mean motion resonance (MMR) with one or other of the giant planets. It may also be in a transient 3:2 MMR with Jupiter. Since Jupiter and Saturn are themselves near 5:2 MMR, the clone is in the region of overlapping resonances with the two giant planets. This is illustrated in the lower panels where the evolution of the resonant angles can be seen. The proximity to the orbits of Jupiter or Saturn means that the respective resonant angles oscillate around |$180\,$|°, but with a half amplitude of |$\sim 180\,$|°. The resonant angles for 1:1 MMR are simply differences between the mean longitudes of the asteroid and one of the planets. Such oscillations of the angles are characteristic of horseshoe orbits in the restricted three-body problem. For the proximity to 3:2 MMR with Jupiter, the resonant angle oscillates around |$\sim 30\,$|° with a moderate amplitude of |$\sim 120\,$|°. In all cases, the oscillations occur for only a few cycles.

If this particular clone of Sedna was the actual object, Sedna may even have come from the group of stable Jupiter Trojans after being perturbed by a collision with another member of the group. Naturally, due to the strong chaos, such predictions cannot be made and one can only get the statistical information from the backward integration experiment. Not all of the 25 ETNOs reach orbits close to Jupiter or Saturn. Some of them reach the range of |$a \sim 20-40 \ \mbox{au}$| and |$e \gtrsim 0.3$| with moderate-to-high inclinations.

Fig. 14 shows the statistical results of the experiment and the comparison with the observed distribution of asteroids. On the left side, the known solar system objects are shown as two-dimensional scans mapping the number density of the objects. They are shown in two diagrams |$(a,e)$| in the top row and |$(i,e)$| in the bottom row. A darker shade of red indicates a higher number density of objects. Most asteroids are located in two regions: Trojans and Plutinos. There are also traces of scattered objects from both groups. Asteroids with perihelia between |$5.2 \ \mbox{au}$| and |$30 \ \mbox{au}$| (Centaurs) are much less numerous. There is an almost empty gap for |$a \gt 40 \ \mbox{au}$| and |$q \in (5.2, 30) \ \mbox{au}$|⁠. The inclinations are mostly moderate |$\lesssim 40\,$|° with less numerous objects with higher values of |$i$| accompanied by higher values of |$e$|⁠.

The right-hand panel of Fig. 14 shows the results of backward integrations of 25 ETNOs with 100 clones for each. Clearly, most of the objects in the |$(a,e)$|-diagram are in the gap mentioned above. Therefore, ETNOs could have come from the disc between Jupiter and Neptune after being scattered. Depending on the semimajor axis reached after the scattering, the detachment from Neptune was slower or faster. If |$a$| is moderate (⁠|$\sim 200-300 \ \mbox{au}$|⁠), the next visit to the region of giant planets occurs after one or more revolutions of the Sun around the centre of the Galaxy. If |$a$| is large (⁠|$\sim 1000 \ \mbox{au}$|⁠), the next visit takes place even within a few Myr.

The |$(i,e)$|-diagram shows that there are very few objects with simultaneously low |$i$| and |$e$|⁠. This area is mainly occupied by Plutinos and Trojans (bottom left), which are less likely to get scattered and join the population of trans-Neptunians. Most of the past ETNOs have |$i$| and |$e$| at the boundary of the grouping of Plutinos and Trojans shown in the lower left-hand panel of Fig. 14. The majority of objects are in the |$i \in (0, 60)\,$|° and |$e \in (0.4, 0.9)$| region therefore ETNOs more likely had prograde orbits in the past.

Every time perihelion is in the planetary region, i.e. an ETNO becomes a Centaur, the evolution is chaotic, but the object is not necessarily ejected from the solar system. However, its orbit is significantly perturbed. An example of such an evolution is shown in Fig. 15 for 2014 FE72 (one of the clones). For |$t = -100\,$| Myr |$a \sim 500 \ \mbox{au}$|⁠, while |$q \sim 200 \ \mbox{au}$| and decreases. After reaching |$q \sim 80 \ \mbox{au}$|⁠, it could be detected (according to current detection capabilities) and classified as ETNO. For |$t \sim -75\,$| Myr, the object has |$a \sim 30-50 \ \mbox{au}$| and |$q \sim 10 \ \mbox{au}$|⁠. After a series of close planetary encounters, the semimajor axis for |$t \sim -25\,$| Myr reached |$1000-2000 \ \mbox{au}$|⁠, as observed today. After several more EFE cycles and close encounters with the planets, the semimajor axis grows to |$\sim 10000 \ \mbox{au}$|⁠.

According to MOND, ETNOs are ephemeral objects rather than stable asteroids held together by a distant planet. Their semimajor axes, perihelion distances as well as inclinations vary significantly over time. A given object can evolve between a low-|$a$| Centaur, a high-|$a$| Centaur, an ETNO, in a prograde or retrograde orbit. Recently, the first three retrograde ETNOs were discovered and their orbits are explained in the next subsection within MOND.

5.2 Retrograde ETNOs

The orbital parameters of the three retrograde ETNOs are not well determined, since for two of them (2022 FN12 and 2022 FM12) there are only six observations covering 2 d, while for the third (2019 EE6) there are seven observations covering 63 d. We use the orbital parameters listed in the JPL Small-Body Database (accessed 2022 October 13; Table 1). All objects are detached from Neptune, |$q \gt 50 \ \mbox{au}$|⁠, and retrograde, |$i \gt 160\,$|°. The formal uncertainties of the orbital elements are not determined.

The question naturally arises whether or not these objects are dynamically connected with prograde ETNOs and Centaurs. To be related to the former, the inclinations should decrease, while to be related to the latter means that the perihelion distances should decrease. The lower left-hand panel of Fig. 5 suggests that this is possible. The only potential problem lies in the time required to evolve into the low-|$q$| or low-|$i$| region of parameter space. The semimajor axes of all three asteroids are relatively low, i.e. |$136.3 \ \mbox{au}$| (2022 FN12), |$158.8 \ \mbox{au}$| (2022 FM12), and |$165.5 \ \mbox{au}$| (2019 EE6). Their EFE-induced evolution thus takes place on time-scales longer than the Sun’s orbital period around the centre of the Galaxy (see Fig. 7). For sufficiently low |$a$|⁠, it could be that planetary perturbations strongly dominate the evolution and significant |$q$| and |$i$| variations are not possible.

In order to verify that, we performed backward integrations with 10 clones for each object chosen according to the Gaussian distribution with assumed uncertainties of the orbital elements of |$\sigma _a = 20 \ \mbox{au}$|⁠, |$\sigma _q = 1 \ \mbox{au}$|⁠, and |$\sigma _{{\mbox{angle}}} = 1\,$|° for the angles. Fig. 16 illustrates the results of integrating selected clones of each retrograde ETNO.

The left column shows the results for 2019 EE6. The perihelion distance decreases and for |$t \sim -150\,$| Myr the asteroid becomes a Centaur. A further decrease of |$q$| down to |$10 \ \mbox{au}$| results in strong perturbations by Saturn. At time |$t \sim -300\,$| My the inclination decreases to |$\sim 0$| and the orbit becomes prograde. Shortly thereafter, the orbit becomes retrograde again. The change in |$i$| is accompanied by a significant change in |$a$| and |$q$|⁠. The left-hand panel of Fig. 17 shows the evolution of the same object in the |$(q, i)$|-diagram, where the U-shaped trajectory can be seen, which is similar to that shown in the lower left-hand panel of Fig. 5. Clearly, 2019 EE6 could have been a prograde ETNO or a Centaur in the past.

In the middle column of Fig. 16 the results for 2022 FM12 are shown. At epoch |$\sim -200\,$| Myr the object becomes a Centaur. The orbit is perturbed by the planets. In particular, the perturbation by Saturn leads to a significant increase in |$a$|⁠, which in turn leads to a faster evolution of |$q$|⁠. The perihelion distance first increases up to |$\sim 200 \ \mbox{au}$| and then decreases again down to the planetary region. At the same time, the inclination decreases to moderate values. The perihelion leaves the planetary region and returns, with the inclination decreasing further. In the epoch |$\sim -600\,$| Myr the inclination is as low as a few degrees. The evolution of this object can be followed in the middle panel of Fig. 17. The trajectory is more complicated than that of the previous object, but a characteristic U-shape can still be seen. Similar to 2019 EE6, 2022 FM12 may also have originated in the Centaurs region and passed through the prograde ETNO phase, with the semimajor axis of a few hundred astronomical units.

The third object 2022 FN12 has the smallest semimajor axis among the retrograde ETNOs, i.e. |$a = 135.9 \ \mbox{au}$|⁠, and therefore it is more difficult to experience a significant variation of |$q$|⁠. More clones (we tested 50 clones) were needed to find a configuration that evolves to a Centaur (the right column of Figs 16 and 17). The initial semimajor axis of this clone is |$162 \ \mbox{au}$|⁠, which is significantly larger than the nominal value. However, the uncertainties in the orbital parameters for the retrograde ETNOs are not known and it is possible that the orbit of 2022 FN12 is wider than currently reported in the database. About 800 Myr are needed to reach |$q \lt 30 \ \mbox{au}$| and about |$3.5\,$|Gyr to reach |$i \lt 90~$|°. The clone reaches the region inside Earth’s orbit and the perihelion distance increases again. The object is eventually scattered by Uranus.

Despite the difficulties of reaching a low-|$q$|/low-|$i$| region by 2022 FN12, the tests described above suggest that EFE could in principle explain the origin of retrograde ETNOs. Future improvements of their orbits and possibly new objects of this type would help to test the MOND hypothesis even better.

5.3 Centaurs

In the previous section it was shown that retrograde ETNOs belong to the same group of trans-Neptunians as prograde ETNOs and Centaurs, only at a different stage of evolution. In this section we show the dynamic relationship between the three subgroups through statistical analysis rather than by looking at the evolution of individual objects. We start with the orbits of 32 Centaurs with |$a \in (250, 2000) \ \mbox{au}$|⁠. For each object, apart from the nominal configuration, we have chosen 10 clones according to the uncertainties of the orbital elements. Such a set of 352 objects is integrated forward and backward for |$20\,$|Myr. The choice of this relatively short period was dictated by the comparison between the synthetic and the observed ETNOs, which must be performed for a given epoch of the Sun’s motion in the Galaxy. We performed the simulations for three different values of |$\alpha = 1.5, 2, 2.5$|⁠.

The left-hand panel of Fig. 18 shows the results for |$\alpha = 1.5$|⁠. The shades of blue in the |$q \gt 30 \ \mbox{au}$| part of the diagram depicts the density of the trajectories in phase space, where darker means higher density. We take into account the observational bias for the ecliptic latitudes |$|b| \lt 30\,$|° and only orbits that meet this criterion are used to calculate the density. The densest region is clearly for |$i \lt 40\,$|°, with a maximum around |$20\,$|°, while the second peak is for |$i \gt 140\,$|°. The smaller histogram above the 2-dimensional scan shows density as a function of |$i$|⁠, assuming |$q \in (30, 90) \ \mbox{au}$|⁠.

This is consistent with the observations of both the prograde and retrograde ETNOs. However, it is noticeable that the two ETNOs with |$i \sim 45\,$|° are in the relatively low density region. The density is slightly higher in this region for |$\alpha = 2$| and 2.5 (see middle and right-hand panels). However, the overall results are qualitatively the same for the three values of |$\alpha$|⁠. This experiment shows that EFE can explain the transition between different subgroups of trans-Neptunians, but the statistics of the orbital elements cannot be used to distinguish between different values of the interpolation parameter |$\alpha$|⁠.

5.4 EFE–induced flows of trans-Neptunians

The uneven distribution of the orbital elements |$\Omega , \omega , \varpi$|⁠, which seem to cluster around certain values, have been interpreted as indicating the existence of an additional planet in the solar system (Batygin & Brown 2016a). Before checking whether the Milgromian dynamics can explain these facts without Planet Nine, we show that there are additional features in the distribution of orbital elements that point to EFE as the main factor for the statistics of ETNOs.

Fig. 11 shows the distribution of |$i$| and |$q$| for TNOs, while the colours and sizes of the symbols encode the values of |$a$|⁠. After dividing the objects into four groups according to |$a$|⁠, i.e. |$a \in (500, 2000) \ \mbox{au}$|⁠, |$a \in (250, 500) \ \mbox{au}$|⁠, |$a \in (150, 250) \ \mbox{au}$| and |$a \in (100, 150) \ \mbox{au}$|⁠, we find that the distributions for each group differ from each other (see Fig. 19). In addition, we can add arrows indicating the direction of evolution of each object when both the EFE and planetary perturbations are considered.

The upper left-hand panel of Fig. 19 shows the |$(i, q)$|– diagram of the asteroids with the widest orbits, i.e. |$a \in (500, 2000) \ \mbox{au}$|⁠. We use the upper limit of |$a$| to exclude extremely wide orbits of Centaurs and keep the maximum |$a$| for them in the same order of magnitude as the maximum |$a$| for ETNOs. The points form a narrow strip that is nearly vertical for |$q \gt 30 \ \mbox{au}$| and nearly horizontal for |$q \lt 10 \ \mbox{au}$|⁠, with a transition between these limiting cases. The arrows show that the flow of objects is consistent with the observed distribution. In the vertical part, the arrows point mainly in the up/down direction, while in the horizontal part, the arrows point mainly in the left/right direction. For wide orbits, the perturbations from the giant planets are much weaker than the EFE perturbations (see Fig. 7) and the EFE should dominate. Such a correspondence between flow and distribution clearly indicates that EFE could be responsible for the distribution.

In the upper middle panel of Fig. 19 the |$(i, q)$|-diagram for |$a \in (250, 500) \ \mbox{au}$| is shown. The velocity distribution is less ordered, especially in the Centaurs region. In addition to the main flow, which is similar to the previous one, a considerable amount of random motion can be seen. The distribution of orbital elements is also less tight. In the ETNOs region (⁠|$q \gt 30 \ \mbox{au}$|⁠), both the distribution of orbital elements and the velocities are relatively well constrained, especially for |$i \lt 30~$|°. The picture shown in this panel is to be expected for the intermediate range of |$a$|⁠. The EFE still dominates over the planetary perturbations for |$a \in (250, 500) \ \mbox{au}$|⁠, but for |$a \sim 250 \ \mbox{au}$| the two perturbations are already of equal magnitude (see Fig. 7).

For even smaller |$a$|⁠, both the distribution and the flow become less confined. The top right-hand and bottom left-hand panels of Fig. 19 show the diagrams for |$a \in (150, 250) \ \mbox{au}$| and |$(100, 150) \ \mbox{au}$|⁠, respectively. The bottom middle and bottom right-hand panels show the close-ups of the ETNOs region. For |$a \in (150, 250) \ \mbox{au}$| the points form a clump with random motions. This is to be expected in the planet–dominated range of |$a$|⁠. For |$a \in (100, 150) \ \mbox{au}$| a positive correlation between |$i$| and |$q$| appears, consistent with the flow. In this range of |$a$|⁠, the planetary perturbations strongly dominate over EFE and the correlation and the flow can be attributed to the conservation of the projection of angular momentum onto the |$z$|-axis. This is expected for the perturbation by the giant planets, which is almost axisymmetric. The absence of low-|$i$| and high-|$q$| objects among TNOs with |$a \lt 150 \ \mbox{au}$| is consistent with the formation scenario described in Nesvorný, Vokrouhlický & Roig (2016), Kaib & Sheppard (2016), and Anderson & Kaib (2021), combining the crossing of MMRs with Neptune during its early migration and the Lidov–Kozai mechanism (Lidov 1962; Kozai 1962).

The experiments described above show that the observed distributions of the orbits agree with the Milgromian model of gravity. In the next section we try to verify whether MOND can explain the non-uniform distribution of the longitudes of the ascending nodes as well as the longitudes of the perihelia of ETNOs.

6 ORBITAL PLANE CLUSTERING AND APSIDAL CONFINEMENT

In the previous section we showed that the evolution of trans-Neptunians is strongly chaotic. Moreover, depending on its evolutionary phase, a given object can be classified as a Centaur or as a prograde/retrograde ETNO with smaller or greater semimajor axis. Some of the known ETNOs may have visited the giant planet region only a few Myr ago, resulting in a significant and unpredictable change in their orbits due to the chaos and uncertainties of the orbits. Other objects were probably scattered from the giant planet region at least a few hundred Myr ago.

This means that it is very difficult, if not impossible, to choose a specific initial epoch and a synthetic initial set of orbits for the numerical simulations, which are then compared with the observational data. Instead, we chose the currently observed objects as initial orbits. The orbits were then integrated forward and backward for |$\pm 400\,$| Myr. The difference from the previous experiment is that we selected all objects with |$a \gt 50 \ \mbox{au}$| and |$q \gt 5 \ \mbox{au}$|⁠. There were 978 objects in our sample. The reason for this selection is that while we use the criteria |$q \in [42, 85] \ \mbox{au}$| and |$a \in [150, 2000] \ \mbox{au}$| for comparison with observations, the initial sample must include all objects whose orbits could possibly evolve into these ranges (at least for a limited time). Furthermore, we are not looking for a particular fixed cluster in parameter space, since the position of the cluster in |$(\Omega , \varpi)$|-space may vary over time and ETNO orbits at a particular epoch may nevertheless be clustered.

Another feature of the distribution of TNOs studied in this experiment was the clustering of orbital planes of objects with |$q\gt 30 \ \mbox{au}$| and relatively small semimajor axes, |$a \in [50, 100] \ \mbox{au}$| (Volk & Malhotra 2017). Unfortunately, our simulations show that MOND does not prevent the randomization of orbits with |$a \in [50, 100] \ \mbox{au}$| after a few Myr. Interestingly, (Brown & Dahlke 2018) show that this is also true for Planet Nine. We discuss a possible solution to this problem in the next section. Here we focus on the clustering of orbits with |$a \gt 150 \ \mbox{au}$|⁠.

Fig. 20 shows the distribution of orbital elements for a few selected epochs in the evolution of the trans-Neptunians. The top row shows the current epoch. Most of the objects have perihelia near the anticentre of the Galaxy (the cluster is shifted by |$\sim 45\,$|° from this direction), only two of them have perihelia near the Galaxy centre. There are gaps between the two groups. A similar picture is repeated in different past and future epochs. The three selected epochs |$-302, -198$|⁠, and |$391\,$| Myr correspond to the situation where the centre/anticentre of the Galaxy is close to the ecliptic. In each of these cases, the perihelia confinement is clearly visible. For higher ecliptic latitudes of the Galaxy centre/anticentre (the bottom two rows), the confinement may be less clear (especially for |$t = 60\,$| Myr).

The right column of Fig. 20 shows the distributions in the respective epochs in the |$(i \cos \Omega , i \sin \Omega)$|-diagram. Clearly, the inclinations in past/future epochs are generally larger than today. Nevertheless, the majority of orbits have |$i \lesssim 60\,$|°. The experiments with Planet Nine show a similar excitation of inclinations (Shankman et al. 2017; Batygin & Brown 2021). One can also observe that |$\Omega$| is unevenly distributed, which is consistent with the observations.

The results presented in Fig. 20 show only selected epochs where both the clustering of orbits and the confinement of perihelia are relatively clear. Fig. 21 shows the full results in the form of density maps for the angles |$\Omega , \varpi , \omega , i$| as a function of epoch. The upper panel shows the map for |$\Omega$|⁠. It is centred at |$90\,$|°, which is close to the current cluster. Within |$\pm 100\,$| Myr most ETNOs are clustered, but the position of the cluster centre varies with time. At more distant epochs, most objects have |$\Omega \in (-\pi /2, +\pi /2)$|⁠, while at past epochs |$\Omega \in (\pi /2, 3 \pi /2)$|⁠. For more distant epochs, the clustering is less clear. However, the sample of objects considered in this study cannot be complete, as some of the objects that are not observed today (e.g. due to |$q \gtrsim 85 \ \mbox{au}$|⁠) may become detectable in a distant epoch and objects observed today are only temporarily detectable.

The second panel from the top shows the density map for |$\varpi$|⁠. It is less clear than that for |$\Omega$|⁠, although configurations with |$\varpi \sim \pm \pi /2$| are preferred. The confinement in |$\varpi$| corresponds to the confinement in the ecliptic longitude of perihelia. Fig. 22 shows a density map in the ecliptic coordinates. The confinement can also be analysed in a one-dimensional representation (upper panel). The configurations with |$\lambda$| near the centre/anticentre of the Galaxy are favoured, but there are also a significant number of objects with |$\lambda \sim 0$|⁠, which is consistent with the observations (the blue symbols). The maximum near the Galaxy centre is at least as strong as that near the anticentre. The observed asymmetry is probably due to the observational bias.

The density map for the argument of perihelion is shown in the third panel of Fig. 21. Configurations with |$\omega$| near 0° or |$180\,$|° are preferred, due to the observational limitation |$|\beta | \lt 30\, \deg$|⁠. In the lower panel of Fig. 21 the map for |$i$| is shown. Most of the objects are in orbits of low-to-moderate inclination. However, there are some ETNOs whose inclination exceeds |$90\,$|° and whose orbits become retrograde. This is consistent with the observations and analysis presented in this paper.

7 DISCUSSION AND FUTURE WORK

The problem of modelling the clustering of orbital planes as well as the apsidal confinement remains not fully solved. On the one hand, the observations are not yet completely conclusive. The features may still be artefacts or are weak. On the other hand, the simulations do not show beyond doubt that the model produces strong non-uniformities in the distribution of the orbital elements, although the results presented above show moderate clustering and confinement. More data are definitely needed to verify the properties of the ETNOs. The model should also be improved.

Such an improvement could be based on the consideration of non-zero asteroid masses. Madigan & McCourt (2016) have shown that the inclusion of interactions between the small masses significantly affects the dynamics. The significance of self-gravitation between the trans-Neptunians has been also demonstrated by Sefilian & Touma (2019), who showed that a massive and moderately eccentric disc of TNOs can balance the giant planets–induced apsidal precession of highly eccentric ETNOs, resulting in their apsidal confinement. Such an extended model of massive TNOs will be investigated within MOND in a future work.

The clustering of orbits is not the only feature that could possibly be explained by the self-gravitating TNOs. The obliquity of the Sun with respect to the invariant plane of the solar system is another puzzle to be solved. The solar spin is tilted by |$I_{\odot } = 7.155 \pm 0.002\,$|° and the orientation of the solar equator is given by |$\Omega _{\odot } = 73.5 \pm 0.1\,$|° (Beck & Giles 2005) with respect to the ecliptic-equinox J2000.0, while the invariant plane has the parameters |$I_{{\mbox{inv}}} = 1.58\,$|°, |$\Omega _{{\mbox{inv}}} = 107.58\,$|° (Souami & Souchay 2012).

There have been several attempts in the literature (Bailey, Batygin & Brown 2016; Gomes, Deienno & Morbidelli 2017) in which authors tried to reconstruct the obliquity as a result of the gravitational interaction of the giant planets with Planet Nine. If the giants are initially in the same plane as the solar equator, the obliquity can be achieved if Planet Nine is in an inclined orbit. The inclinations of the giant planets’ orbits increase due to the exchange of angular momentum with Planet Nine. However, there is a problem with such a scenario. The initial invariant plane, defined as the total angular momentum of five planets, i.e. the known giants and Planet Nine, must initially be inclined with respect to the solar equator. However, it could work if Planet Nine was a free-floating planet captured by the Sun. If Planet Nine was formed together with the known planets in a common disc, the initial inclination between P9 and the four giants would be small and could not increase.

Another way of solving this puzzle is to notice that the planets alter the solar spin (Lai 2016), which precesses around the vector of the total angular momentum of the planets. However, even in this case, the initial obliquity is required, which cannot be generated from initial values |$\sim 0$|⁠. In contrast, Milgrom’s gravity can lead to inclination excitation even if the inclinations are initially zero.

As explained earlier in this paper, the change in inclination results from the precession of the orbit with respect to the axis of symmetry of EFE (the Galaxy centre-anticentre axis). The mechanism is very efficient for trans-Neptunians, but also works for giant planets, just much less efficiently. Fig. 23 shows the backward evolution of |$i$| and |$\Omega$| of the Laplace plane of Jupiter, Saturn, Uranus, and Neptune perturbed by EFE for |$\alpha = 1.5$| (red) and |$\alpha = 2$| (blue). Clearly, the inclination increases by |$\sim 3\,$|° over the lifetime of the solar system. This is still not enough to reconstruct the observed solar obliquity, but more importantly, the evolution of |$\Omega$| is not consistent with the observed solar system parameters. In the right-hand panel of Fig. 23, the evolution is shown in the |$(i \cos \Omega , i \sin \Omega)$|-diagram. The spin of the Sun also evolves due to the interaction with the planets (we used a simple model in Lai 2016). The obliquity does not decrease for past epochs, but even increases slightly.

Trans-Neptunian objects with non-zero mass could hopefully solve this problem, but non-restricted large-|$N$|-body simulations are beyond the scope of this article and we defer them to future studies. The scenario has a few elements: (i) EFE effectively increases the inclinations of TNOs and much weaker the inclinations of giant planets; (ii) TNOs interact with giant planets, especially strongly in epochs of low |$q$|⁠; (iii) giant planets alter the spin of the Sun. All elements are important and should be modelled in a self-consistent way.

Since MOND is presented here as an alternative to Planet Nine, one might ask whether both hypotheses can be correct. If gravity is Milgromian, then Planet Nine would obey chaotic evolution similar to other trans-Neptunian objects. Its perihelion would repeatedly return to the planetary region every few tens or hundreds of Myr, depending on the semimajor axis, showing behaviour similar to that in Fig. 15. The semimajor axis of Planet Nine would change each time perihelion |$q \lesssim 30 \ \mbox{au}$|⁠, but the giant planets would also be significantly perturbed. If both hypotheses are correct, the stability of the solar system would be in question.

An example of the backward and forward evolution of the five-planet solar system within Milgromian dynamics is shown in Fig. 24. Planet Nine visits the inner solar system in epochs |$\pm 200\,$| Myr, which leads to significant changes in its orbit (destabilization in epoch |$-200\,$| Myr) and/or the orbits of the giant planets, i.e. in epoch |$+200\,$| Myr the orbits of Uranus and Neptune are strongly perturbed.

However, the stability problem in the past could be overcome if Planet Nine was a free-floating planet (Li & Adams 2016) or a primordial black hole (Scholtz & Unwin 2020) captured by the Sun relatively recently, i.e. less than |$\sim 200\,$| Myr ago. On the other hand, the future fate of the solar system could still be potentially catastrophic in this scenario.

8 CONCLUSIONS

The hunt for Planet Nine continues by selecting the most likely positions in the sky (Brown & Batygin 2021; Socas-Navarro 2023), calculating theoretical fluxes in various bands (Wright et al. 2010; Cowan, Holder & Kaib 2016; Ginzburg, Sari & Loeb 2016; Linder & Mordasini 2016; Planck Collaboration 2014, 2016) and conducting/planning surveys (Trilling, Bellm & Malhotra 2018; Naess et al. 2021; Rowan-Robinson 2022). So far, Planet Nine has not been discovered. According to the flux models, the planet is detectable, but a large area of the sky must be searched with large telescopes to detect or rule out the putative planet of our solar system.

In this paper we have investigated the effects of MOND on the evolution and dynamical origin of ETNOs. We found that the observed objects could originate from the planetary region of the solar system, with semimajor axes on the order of a few to a few tens of astronomical units. A striking example is Sedna, which could have been in a horseshoe orbit with Jupiter and Saturn about |$30\,$| Myr ago. According to this scenario, ETNOs were scattered by the giant planets into orbits with high-|$a$| and |$q \in (5, 30) \ \mbox{au}$|⁠, which were then detached from the influence of the planets by EFE. According to MOND, both ETNOs and Centaurs are ephemeral objects whose evolution is strongly chaotic.

Milgromian dynamics explains the formation of ETNOs in a natural way. Particularly interesting objects are the newly discovered retrograde ETNOs with relatively small |$a$| and large |$q$|⁠, which have been shown to belong to the same dynamical class as the already known prograde ETNOs and Centaurs. These objects are dynamically unconnected when Newtonian dynamics is used. MOND produces the clustering of orbital planes as well as the apsidal confinement, but these features are not particularly strong in this model. On the other hand, the observations are still inconclusive. More new objects are clearly needed to enrich the sample and confirm or rule out the anomalous features of ETNOs.

We have shown that MOND can produce solar obliquity, but the observed parameters cannot be fully reconstructed if the giant planets are only modified by EFE. We suggest that accounting for the masses of the TNOs could hopefully resolve this discrepancy.

We have considered whether the two alternative explanations for the distribution of ETNOs, namely Planet Nine and MOND, can be correct at the same time. We have shown that the dynamical stability of the solar system with Planet Nine can be problematic if the gravity is Milgromian. On the other hand, Planet Nine could have disrupted the orbits of the known giant planets and broken the chain of mean-motion resonances that are thought to have formed early in the history of the solar system, but not necessarily destabilize the entire system. A scenario in which the initially resonant and more compact solar system has been significantly disrupted by planetesimals has been proposed in the Nice model (Gomes et al. 2005; Morbidelli et al. 2005; Tsiganis et al. 2005). Naturally, if gravity is Milgromian and Planet Nine exists, the history of the solar system would have to be reconsidered. In any case, the Milgromian dynamics of the solar system with Planet Nine is much more complicated than without it and requires further investigation.

MOND has been considered here as an alternative for the Planet Nine hypothesis. However, MOND is also a well-known alternative to the dark matter hypothesis. Since falsification of dark matter is problematic when considering galactic and extra-galactic systems, the dynamics of the solar system can serve as a testing ground for the dark matter hypothesis. In other words: If Planet Nine is not found and the anomalous features of ETNOs are confirmed by future observations and fully explained by MOND, then the dark matter hypothesis could be in trouble.

Acknowledgement

I would like to thank Dan Fabrycky for bringing the Planet Nine hypothesis to my attention and Dr Man Ho Chan for helpful comments that improved the paper.

DATA AVAILABILITY

The orbital parameters of TNOs used in this work were taken from the JPL Small-Body Database https://ssd.jpl.nasa.gov/tools/sbdb_query.html (accessed 2022 October 13). The orbital parameters were computed by JPL' Solar System Dynamics group. The actual values used in the simulations as well as other data underlying this article will be shared on reasonable request to the corresponding author. The uncertainties of the solar system planets masses were taken from the JPL database (Planetary Physical Parameters section), https://ssd.jpl.nasa.gov/planets/phys_par.html (accessed 2023 July 6). The physical parameters were collected by JPL's Solar System Dynamics group.

Footnotes

References