Towards an explanation of orbits in the extreme trans-Neptunian region: The effect of Milgromian dynamics

TL;DR

This paper investigates whether Milgromian dynamics can explain the unusual orbits of extreme trans-Neptunian objects by modeling the external field effect and comparing predictions with observations.

Contribution

It introduces a dynamical model incorporating Milgromian external field effects to explain ETNO orbital distributions, constrained by Cassini data.

Findings

EFE can produce non-uniform ETNO orbital distributions

Predicted distributions conflict with observations unless EFE parameters are unlikely

Milgromian dynamics alone may not fully explain ETNO orbital characteristics

Abstract

Milgromian dynamics (MD or MOND) uniquely predicts motion in a galaxy from the distribution of its stars and gas in a remarkable agreement with observations so far. In the solar system, MD predicts the existence of some possibly non-negligible dynamical effects, which can be used to constrain the freedom in MD theories. Known extreme trans-Neptunian objects (ETNOs) have their argument of perihelion, longitude of ascending node, and inclination distributed in highly non-uniform fashion; ETNOs are bodies with perihelion distances greater than the orbit of Neptune and with semimajor axes greater than 150 au and less than $\sim1500$ au. It is as if these bodies have been systematically perturbed by some external force. We investigated a hypothesis that the puzzling orbital characteristics of ETNOs are a consequence of MD. We set up a dynamical model of the solar system incorporating the external field effect (EFE), which is anticipated to be the dominant effect of MD in the ETNOs region. We used constraints available on the strength of EFE coming from radio tracking of the Cassini spacecraft. We performed several numerical experiments, concentrating on the long-term orbital evolution of primordial (randomised) ETNOs in MD. The EFE could produce distinct non-uniform distributions of the orbital elements of ETNOs that are related to the orientation of an orbit in space. If we demand that EFE is solely responsible for the detachment of Sedna and 2012 VP$_{113}$, then these distributions are at odds with the currently observed statistics on ETNOs unless the EFE quadrupole strength parameter $Q_{2}$ has values that are unlikely (with probability < 1$\%$) in light of the Cassini data.

Tables

Table 1. Table 1: Results of the KS test comparing distributions of ω 𝜔 \omega , Ω Ω \Omega , L 𝐿 L , and i 𝑖 i of the observational data and those inferred from the benchmark model simulation at t = 4 𝑡 4 t=4 Gyr. The KS test was also carried out for the P9 simulation. Each simulation starts with 400 particles; p ​ ( ξ | ζ ) 𝑝 conditional 𝜉 𝜁 p(\xi\leavevmode\nobreak\ |\leavevmode\nobreak\ \zeta) is the probability p-value, where ξ 𝜉 \xi are the observational data in a quoted a 𝑎 a -range, and ζ 𝜁 \zeta are the simulated data in a quoted a 𝑎 a -range, e.g. p ​ ( 150 − 250 | 150 − 250 ) 𝑝 150 conditional 250 150 250 p(150-250\leavevmode\nobreak\ |\leavevmode\nobreak\ 150-250) is the p-value comparing the observational and simulated data in range 150 < a < 250 150 𝑎 250 150<a<250 au. The p 𝑝 p -values larger than 0.5 are highlighted in boldface. There was only one particle with a > 250 𝑎 250 a>250 au left at the end of the P9 simulation, therefore p ( > 250 | > 250 ) p(>250\leavevmode\nobreak\ |>250) is not reported in this case; [ Q 2 ] delimited-[] subscript 𝑄 2 [Q_{2}] = = 10 − 27 superscript 10 27 10^{-27} s -2 .

$Q_2=0$	$\omega$	$\Omega$	$L$	$i$
$p\left(150-250|150-250\right)$	0.17	0.31	0.66	0.00
$p(>250|>250)$	0.00	0.03	0.04	0.00
P9 (#1)	$\omega$	$\Omega$	$L$	$i$
$p\left(150-250|150-250\right)$	0.95	0.25	0.07	0.98
P9 (#2)	$\omega$	$\Omega$	$L$	$i$
$p\left(150-250|150-250\right)$	0.48	0.48	0.83	1.00
$Q_2=1$	$\omega$	$\Omega$	$L$	$i$
$p\left(150-250|150-250\right)$	0.19	0.03	0.71	0.00
$p(>250|>250)$	0.02	0.00	0.01	0.00
$Q_2=2$	$\omega$	$\Omega$	$L$	$i$
$p\left(150-250|150-250\right)$	0.25	0.01	0.95	0.00
$p(>250|>250)$	0.01	0.01	0.01	0.02
$Q_2=4$	$\omega$	$\Omega$	$L$	$i$
$p\left(150-250|150-250\right)$	0.35	0.00	0.66	0.00
$p(>250|>250)$	0.03	0.00	0.04	0.49
$Q_2=6$	$\omega$	$\Omega$	$L$	$i$
$p\left(150-250|150-250\right)$	0.26	0.01	0.80	0.07
$p(>250|>250)$	0.00	0.00	0.01	0.89

Table 2. Table 2: Results of the KS test taking into account only the orbits, of ETNOs and test particles, with q > 35 𝑞 35 q>35 au.

P9 (#1)	$\omega$	$\Omega$	$L$	$i$
$p\left(150-250|150-250\right)$	0.71	0.71	0.28	0.53
$Q_2=1$	$\omega$	$\Omega$	$L$	$i$
$p\left(150-250|150-250\right)$	0.11	0.18	0.66	0.00
$p(>250|>250)$	0.02	0.00	0.01	0.00
$Q_2=6$	$\omega$	$\Omega$	$L$	$i$
$p\left(150-250|150-250\right)$	0.33	0.02	0.99	0.22
$p(>250|>250)$	0.01	0.00	0.01	0.64

Table 3. Table 3: Results of the KS test after doubling the initial number of test particles. We only considered the orbits of ETNOs and test particles with q > 35 𝑞 35 q>35 au. There was only one particle with a > 250 𝑎 250 a>250 au left at the end of the P9 simulation, therefore p ( > 250 | > 250 ) p(>250\leavevmode\nobreak\ |>250) is not reported in the P9 case.

P9 (#1)	$\omega$	$\Omega$	$L$	$i$
$p\left(150-250|150-250\right)$	0.98	0.57	0.98	0.57
$Q_2=1$	$\omega$	$\Omega$	$L$	$i$
$p\left(150-250|150-250\right)$	0.09	0.27	0.85	0.00
$p(>250|>250)$	0.01	0.00	0.01	0.00
$Q_2=6$	$\omega$	$\Omega$	$L$	$i$
$p\left(150-250|150-250\right)$	0.19	0.09	0.86	0.01
$p(>250|>250)$	0.00	0.00	0.03	0.51

Table 4. Table 4: Results of the KS test comparing distributions of ω 𝜔 \omega , Ω Ω \Omega , L 𝐿 L , and i 𝑖 i of the observational data and those inferred from the benchmark model simulation at t = 4 𝑡 4 t=4 Gyr. Q 2 subscript 𝑄 2 Q_{2} is assumed to be negative, fulfilling | Q 2 | < 6 subscript 𝑄 2 6 |Q_{2}|<6 . [ Q 2 ] delimited-[] subscript 𝑄 2 [Q_{2}] = = 10 − 27 superscript 10 27 10^{-27} s -2 .

$Q_2=0$	$\omega$	$\Omega$	$L$	$i$
$p\left(150-250|150-250\right)$	0.17	0.31	0.66	0.00
$p(>250|>250)$	0.00	0.03	0.04	0.00
P9 (#1)	$\omega$	$\Omega$	$L$	$i$
$p\left(150-250|150-250\right)$	0.95	0.25	0.07	0.98
P9 (#2)	$\omega$	$\Omega$	$L$	$i$
$p\left(150-250|150-250\right)$	0.48	0.48	0.83	1.00
$Q_2=-1$	$\omega$	$\Omega$	$L$	$i$
$p\left(150-250|150-250\right)$	0.15	0.15	0.81	0.00
$p(>250|>250)$	0.00	0.85	0.01	0.00
$Q_2=-2$	$\omega$	$\Omega$	$L$	$i$
$p\left(150-250|150-250\right)$	0.19	0.27	0.98	0.00
$p(>250|>250)$	0.02	0.95	0.07	0.00
$Q_2=-4$	$\omega$	$\Omega$	$L$	$i$
$p\left(150-250|150-250\right)$	0.30	0.48	0.65	0.01
$p(>250|>250)$	0.00	0.90	0.04	0.19
$Q_2=-6$	$\omega$	$\Omega$	$L$	$i$
$p\left(150-250|150-250\right)$	0.19	0.53	0.66	0.01
$p(>250|>250)$	0.01	0.68	0.01	0.40

Equations

$$\mathbf{g}=\nu (g^N/a_0){\mathbf{g}}^N\text{leavevmode},$$

$$\nabla \cdot \left[\mu \left(\frac{|\nabla \Phi |}{a_0}\right)\nabla \Phi \right]\text{leavevmode}=\text{leavevmode}4\pi G\rho \text{leavevmode}=\text{leavevmode}{\nabla }^2{\Phi }_N\text{leavevmode},$$

$${\nabla }^2\Phi \text{leavevmode}=\text{leavevmode}\nabla \cdot \left[\nu \left(\frac{|\nabla {\Phi }_N|}{a_0}\right)\nabla {\Phi }_N\right]\text{leavevmode}.$$

$$\delta \mathbf{g}=\left[\nu (g^N/a_0)-1\right]{\mathbf{g}}^N\text{leavevmode}.$$

$$\Phi$$

$$Q_{ij}$$

$$e_{x-g}$$

$$e_{y-g}$$

$$e_{z-g}$$

$$e_x$$

$$e_y$$

$$e_z$$

$$\delta {\Phi }_{oct}$$

Full text

*Aims. *We investigated a hypothesis that the puzzling orbital characteristics of ETNOs are a consequence of MD.

*Methods. *We set up a dynamical model of the solar system incorporating the external field effect (EFE), which is anticipated to be the dominant effect of MD in the ETNOs region. We used constraints available on the strength of EFE coming from radio tracking of the Cassini spacecraft. We performed several numerical experiments, concentrating on the long-term orbital evolution of primordial (randomised) ETNOs in MD.

Results. The EFE could produce distinct non-uniform distributions of the orbital elements of ETNOs that are related to the orientation of an orbit in space. If we demand that EFE is solely responsible for the detachment of Sedna and 2012 VP113, then these distributions are at odds with the currently observed statistics on ETNOs unless the EFE quadrupole strength parameter $Q_{2}$ has values that are unlikely (with probability ¡ 1$\%$) in light of the Cassini data.

1 Introduction

History has taught us that we should pay attention to orbital anomalies in the solar system. Two of the best historical examples are the discovery of Neptune and nonexistence of the planet Vulcan. Neptune was found right where it was predicted to erase an anomaly in the orbit of Uranus (Le Verrier, 1846), while the non-existence of Vulcan finally led to the development of new revolutionary ideas in physics – the geometrical interpretation of gravity – explaining precisely the minute anomalous precession of Mercury (Le Verrier, 1859; Einstein, 1916)

Orbital characteristics of observed extreme trans-Neptunian objects (ETNOs), objects having perihelion distances, $q$, greater than the orbit of Neptune, and semimajor axes, $a$, greater than 150 au and less than roughly 1500 au (definition established and motivated in Trujillo & Sheppard, 2014), suggest that these objects have been systematically perturbed by some external force. On September 1, 2016, the Minor Planet Center (MPC) orbit database registered 19 ETNOs111All ETNOs have $q>33$ au, 14 of 19 ETNOs have $q>35$ au, and 8 of 19 have $q>40$ au.. The outlying positions of 90377 Sedna and 2012 VP113 (hereafter, collectively referred to as Sednoids) in the $q$ versus eccentricity, $e$, diagram (e.g. Fig. 2 in Sheppard & Trujillo, 2016) cannot be addressed by the interaction with the known solar system planets and the current Galactic environment (e.g. Dones et al., 2015). Nevertheless, orbits of ETNOs are also characterised by argument of perihelion, $\omega$, longitude of the ascending node, $\Omega$, and perihelion longitude, $L$,222Longitude of an object at perihelion in the heliocentric ecliptic coordinates. It is not equal to the longitude of perihelion, $\varpi\equiv\omega+\Omega$, except for the case of nil inclination. Similarly, when we refer to perihelion latitude, $B$, we mean latitude of an object at perihelion in the heliocentric ecliptic coordinates. Orbital elements $\omega$, $\Omega$, and $i$ give $L$ and $B$ unambiguously: $\sin(L-\Omega)=\sin\omega\cos i/\cos B$, $\cos(L-\Omega)=\cos\omega/\cos B$. which are distributed in a highly non-uniform and clustered fashion (Trujillo & Sheppard, 2014, de la Fuente Marcos & de la Fuente Marcos, 2014, Brown & Firth, 2016, Batygin & Brown, 2016). This is an unexpected observation, as these orbital elements are expected to circulate rapidly under the action of the giant planets.

Tightly clustered are also the eccentricity, $e$ ($e=0.82\pm 0.06$), and inclination, $i$ ($i=20\pm 7$ deg) of the ETNOs, where in brackets we present mean $\pm 1\sigma$ values. The observed distribution of eccentricities is the result of a simple bias. We preferentially find objects close to the Earth, hence objects with smaller perihelia and hence larger eccentricities are preferred. Such bias means that most objects should be found with eccentricities in the range 0.8 – 0.9 (de la Fuente Marcos & de la Fuente Marcos, 2014). However, the clustering in $i$ should be of dynamical origin as the intrinsic bias in declination owing to our observations made from the Earth, prefers lower inclinations333We add to this that many surveys are performed preferentially close to the ecliptic. (de la Fuente Marcos & de la Fuente Marcos, 2014).

Batygin & Brown (2016) recognised that orbits of long-periodic ETNOs with $a>250$ au (six such objects were known at that time444Asteroid 2000 CR105 with $a=228$ au is very close to the six in terms of $\omega$, $\Omega$, and $L$, however, asteroid 2001 FP185, with $a=226$ au, is about 90 deg ahead in $L$. Hence, the boundary was identified to be at $a=250$ au.: Sedna, 2004 VN112, 2007 TG422, 2010 GB174, 2012 VP113, and 2013 RF98) are very tightly clustered in terms of $\omega$, $\Omega$, and $L$, and hence, aligned in physical space. This observational fact, plus the large perihelia of Sednoids, has led Batygin & Brown (2016) (see also Brown & Batygin, 2016) to hypothesise the existence of the ninth planet of the solar system, dubbed Planet 9 (hereafter P9), with mass of $\sim 10$ $M_{\oplus}$, perihelion longitude 180 deg ahead of the six, orbiting in an eccentric, $a\sim 700$ au, $e\sim 0.6$, and moderately inclined orbit. Discussion in the astronomical community settled in a consensus that the orbital characteristics of ETNOs are not result of an observational bias and/or smallness of the sample, and the quest of finding P9 was prompted (Brown & Batygin, 2016; Cowan et al., 2016; de la Fuente Marcos & de la Fuente Marcos, 2016; de la Fuente Marcos et al., 2016; Fortney et al., 2016; Fienga et al., 2016; Holman & Payne, 2016b, a). Recently, Sheppard & Trujillo (2016) found four new ETNOs, two of which have $a>250$ au (2014 SR349 and 2013 FT28). These two ETNOs have $\omega$ and $i$ clustered in line with the original six, but one of them, 2013 FT28, has $\Omega$ and $L$ with 90 – 180 deg offset from the original clustering positions.

Shankman et al. (2017) demonstrated a drawback of the P9 scenario; they showed that P9 generically drives ETNOs into giant planets crossing orbits, which subsequently decouples them from any shepherding mechanism of P9 and randomises their $\omega$, $\Omega$, and $L$. Hence, these randomised ETNOs should contaminate the observed sample. However, no such population is observed. In Lawler et al. (2017), the authors evolved primordial planetesimal disk for 4 Gyr in the solar system with the present-day planetary architecture plus an additional 10 $M_{\oplus}$ planet. Both circular (Trujillo & Sheppard, 2014) and eccentric (Batygin & Brown, 2016) orbits of the ninth planet were considered. They found no evidence for orbital clustering in the ETNOs region at the end of the simulation neither for a circular nor an eccentric ninth planet. These findings bring some controversy into the P9 scenario.

Knowing the distribution of stars and gas (baryons), the matter we see, the motion in galaxies is determinable. This is implied by the existence of the Galactic laws, such as the baryonic Tully-Fisher relation (McGaugh et al., 2000), the relation between stellar and dynamical central surface densities (Lelli et al., 2016), and the mass discrepancy-acceleration correlation (Sanders, 1990; McGaugh, 2004; Desmond, 2017). All these can be encapsulated into one law called the radial acceleration relation (McGaugh et al., 2016; Lelli et al., 2017) – a very tight correlation between the radial acceleration traced by rotation curves and the acceleration predicted by the observed distribution of baryons. There is no reference to dark matter in these laws, even in systems that are supposedly dark matter dominated, i.e. systems showing large mass-discrepancy in Newtonian dynamics. Additionally, “for any feature in the luminosity profile there is a corresponding feature in the rotation curve and vice versa” (Sancisi, 2004). In the context of dark matter halos, the existence of these laws means that there is a very tight correlation between the baryonic and dark matter distributions. Such a very tight correlation, independent of properties and histories of individual galaxies, is unexpected in the current cosmological paradigm: the $\Lambda$CDM model, where $\Lambda$ is the positive cosmological constant representing dark energy, whose nature is unknown, and CDM stands for cold (non-baryonic) dark matter. All these laws were predicted a priori and follow naturally from the basic principles of Milgromian dynamics (Milgrom, 2014, 2016b).

Milgromian dynamics (MD or MOND = modified Newtonian dynamics; Milgrom, 1983b), also known as scale-invariant dynamics, uniquely predicts how the distribution of stars and gas in a galaxy relates to its dynamics. For reviews on MD, see for example the Famaey & McGaugh (2012), which is thorough and technical, or Milgrom (2016a), which contains a discussion on the potential connection of MD to fundamental physics. Milgromian dynamics departs from Newtonian dynamics (and general relativity) at low accelerations with respect to some predefined frame (this could be e.g. the cosmic microwave background rest frame), and the definition of low acceleration is determined by a new constant of nature, $a_{0}$, having units of acceleration. In the low acceleration regime, unlike Newtonian dynamics and general relativity, MD becomes space-timescale invariant (Milgrom, 2009b). While MD was invented as an alternative to the dark matter hypothesis, which lacks direct observational evidence, over the years MD acquired strong empirical support; see Famaey & McGaugh (2012), which summarises the ubiquitous presence of $a_{0}$ in the data.

In some versions of MD (typically in modified gravity versions of MD), the dynamics of a small subsystem, freely falling in the field of a large external system, is affected by the external free-fall acceleration through the so-called external field effect (EFE; Milgrom, 1983b; Famaey & McGaugh, 2012). This is the case of, for instance, a dwarf galaxy in the field of a host galaxy, or the solar system falling freely in the field of the Galaxy. The EFE arises from generic non-linearity of equations of MD and from the fact that it is based on acceleration. McGaugh & Milgrom (2013a) used MD to predict velocity dispersions of dwarf satellite galaxies of the Andromeda galaxy. Later, in McGaugh & Milgrom (2013b), they compared their a priori predictions with newly available data. They reported that, “The data are in good agreement with our specific predictions for each dwarf made a priori… MOND distinguishes between regimes where the internal field of the dwarf, or the external field of the host, dominates. The data appear to recognise this distinction, which is a unique feature of MOND not explicable in $\Lambda$CDM.” The EFE can also explain the declination of rotation curves in the outer parts of galaxies (Haghi et al., 2016). In the outer parts, when the internal accelerations are smaller than the external fields (which are $\lesssim a_{0}$), the quasi-Newtonian behaviour is restored (Famaey & McGaugh, 2012).

Recently, McGaugh (2016) predicted velocity dispersion, $\sigma_{vel}$, of Crater II in MD. Although internal accelerations are very low $(\ll a_{0})$, EFE reduces $\sigma_{vel}$ down to $\sigma_{vel}\approx 2$ km s*-1*. For isolated Crater II (EFE is not present e.g. in modified inertia theories of MD) McGaugh obtained $\sigma_{vel}\approx 4$ km s*-1*. Significantly higher $\sigma_{vel}$ would mean falsification of MD. In $\Lambda$CDM, higher $\sigma_{vel}$ is expected (Boylan-Kolchin et al., 2012; McGaugh, 2016). The resolution was found recently to be $\sigma_{vel}=2.7\pm 0.3$ km s*-1* (Caldwell et al., 2016). This nicely demonstrates predictive power of MD in galaxies.

It was recognised that EFE may act even in systems with high internal accelerations (Milgrom, 2010). Paučo & Klačka (2016) reconsidered the hypothesis of the Oort cloud of comets in the framework of MD. They also discussed that EFE could produce555This depends strongly on the shape of MD interpolating function and the value of $a_{0}$. These are free parameters in the theory for the time being (they should follow from a more profound theory), but they are constrained by the effect MD produces in the inner solar system (Hees et al., 2016). enough torque to deliver primordial solar system bodies to Sedna-like orbits even in the present-day Galactic environment. Here we reassessed this problem and checked whether the perihelia of both Sednoids can be acquired by the action of EFE. We used recent constraints available on the strength of EFE in the ETNOs region coming from the radio tracking of the Cassini spacecraft (Blanchet & Novak, 2011; Hees et al., 2014, 2016). Assuming that the statistics on the orbital elements of ETNOs are representative of the real population, we also investigated whether the spatial orientation of the orbits of ETNOS can be explained by this effect of MD as well.

We further discuss the observational data on ETNOs in Sect. 2. Milgromian dynamics and the dynamical effect it has on the solar system are introduced in Sect. 3. We present a dynamical model of the solar system obeying laws of MD that are applicable in the ETNOs region in Sect. 4. The comparison between the P9 scenario in Newtonian dynamics and the MD scenario without an additional planet is made in Sections 5 and 6 in the context of the Sednoid perihelia and the spatial orientation of the orbits of ETNOs. We conclude with a summary and discussion of our results in Sect. 7.

2 More on data

Here we discuss clues to the existence of the external perturbation, which give us clues to the existence of the external perturbation. Sheppard & Trujillo (2016) noted that the value $q=35$ au divides visually the observed population of objects with $a>150$ au and perihelia among and beyond giant planets in two distinct populations. This $q$-filter removes five objects (2013 UH15, 2010 VZ98, 2001 FP185, 2013 FS28, and 2015 SO20) from the ETNOs sample. One of these objects, 2013 FS28, is an outlier in terms of $\omega$, with $\omega=101.5$ deg (Sheppard & Trujillo, 2016). Histogram distributions of the orbital elements $\omega$, $\Omega$, and $i$, as well as perihelion longitude $L$ for $100<a<1500$ au, $q>35$ au TNOs listed in the MPC orbit database on September 1, 2016 (colour histograms) are depicted in Fig. 1. The largest $a$ of all ETNOs has Sedna, $a\approx 500$ au, i.e. no TNOs with $500<a<1500$ au, $q>30$ au are known. Figure 2 considers the latter group of $100<a<1500$ au, $q<35$ au orbits. These could be substantially influenced by Neptune and the other giant planets. In order to illustrate the effect of setting the boundary $q$ to be $q=35$ au, we also show, in both figures, the distributions inferred from all well-determined orbits of ETNOs (dashed histogram); ETNOs were defined in Trujillo & Sheppard (2014) as objects with $150<a<1500$ au, $q>30$ au. We have to bear in mind that this choice of the boundary $q$ may be just our means of forcing the data to behave as we expect.

The two $q-$populations are distinct in terms of $\omega$ and $\Omega$ when $a>150$ au. In the case of $\omega$, there are no ETNOs with $q>35$ au, but many objects with $a>150$ au, $q<35$ au, occupying the region around $\omega=180$ deg. Object 2003 SS422 is listed in the NASA Jet Propulsion Laboratory (JPL) small body database with orbital elements (and their $1\sigma$ uncertainties) at Epoch 2457600.5 (July 31, 2016): $a=193.5\pm 47.8$ au, $q=39.4\pm 1.0$ au, $i=16.8\pm 0.1$ deg, $\omega=209.8\pm 17.2$ deg, and $\Omega=151.1\pm 0.1$ deg. This object was not counted because of its large uncertainty in $a$. If this object would be qualified as ETNO with future observations, then it would be the first ETNO with $q>35$ au with $\omega$ near 180 deg, which is about 180 deg away from the presently established clustering position. Additionally, it seems that $q<35$ au objects prefer higher values of $\Omega$ than their $q>35$ au counterparts. There is no obvious distinction between the two populations in terms of $L$. The group of lesser perihelia also comprises objects on highly inclined and retrograde orbits.

Looking at TNOs in Fig. 1, the non-uniformity of the distributions is apparent, with a non-uniformity degree being $a$-dependent, as can be expected when an external perturbation is present. Orbits with low inclinations (¡10 deg) and $q>35$ au are very rare, regardless of $a$, which is in stark contrast with the expected observational bias (e.g. de la Fuente Marcos & de la Fuente Marcos, 2014).

On September 1, 2016, the MPC orbit database listed 40 TNOs with $a>100$ au, 19 ETNOs, and 8 long-periodic ETNOs with $a>250$ au. When we restrict ourselves to $q>35$ au TNOs only, these numbers reduce to 32 TNOs with $a>100$ au, 14 ETNOs, and the same number of 8 long-periodic ETNOs. Clearly, we are dealing with small-number statistics. However, as pointed out elsewhere (Brown & Firth, 2016; Batygin & Brown, 2016), it is unlikely that the observed clustering is purely due to chance. Assuming the $\omega$, $\Omega$, and $L$ of ETNOs are uniformly distributed, the circulation of these angles is ensured by the perturbations from the giant planets, the probability of the observed clustering in $\omega$, $\Omega$, and $L$, beyond $a=250$ au level, is $0.1\%$, $3.0\%$, and $5.0\%$, respectively. These estimates come from a following Monte Carlo simulation a la Batygin & Brown (2016): we drew 8 angles, $\theta_{n}$ , which we call a set, randomly from a uniform distribution in range from 0 to 360 deg, $10^{6}$ times, and calculated the circular variance of each set, $\sigma^{2}_{i}=1-R$, $R=\left(\overline{S}^{2}+\overline{C}^{2}\right)^{1/2}$, $\overline{S}=\left(1/8\right)\sum_{n=1}^{8}\sin\theta_{n}$, $\overline{C}=\left(1/8\right)\sum_{n=1}^{8}\cos\theta_{n}$. The probability is given as a frequency of sets for which $\sigma^{2}_{i}$ is smaller than the circular variance of the observational data.

Brown & Firth (2016) plotted the distribution of $\omega$ for the whole known TNOs population (taking the data from the JPL small body database) and for a subset of TNOs that might not be in mean motion resonances with Neptune. These authors found that the distribution of $\omega$ is well modelled, especially in the latter case, by overlapping normal distributions, peaking at 0 (360) and 180 deg, with standard deviations of 60 deg. This could be the effect of observational bias (de la Fuente Marcos & de la Fuente Marcos, 2014; Sheppard & Trujillo, 2016). Using the mentioned distribution made of the overlapping normal distributions, the Monte Carlo estimate of the probability of the observed clustering in $\omega$ is $0.5\%$ for $a>250$ au orbits.

Observational bias can be an important aspect of the puzzle. A thorough examination of biases and discussion on this topic can be found, for example in Sheppard & Trujillo (2016). The known bias in $\omega$ is twofold: (1) if one does observations close to the ecliptic (as many surveys do), $\omega$s close to $\omega=0$ (360) deg and 180 deg are preferred because then the perihelion of an object, which is the most likely position in which an object is discovered, is close to the ecliptic and (2) observations from the southern hemisphere tend to prefer $\omega\in(180,360)$ deg, and by contrast, $\omega\in(0,180)$ is preferred when observations are made from the northern hemisphere. By employing off-ecliptic observations, made primarily from Chile, and taking into account the fact that there are no $q>35$ au ETNOs with $\omega$ close to 180 deg, Sheppard & Trujillo (2016) concluded that the clustering in $\omega$ for $q>35$ au ETNOs is real with high probability.

The parameter $L$ can be biased as well. New objects are most likely to be discovered near perihelion, and most surveys avoid the star polluted Galactic plane. Hence, there are two main time windows of the year to carry out a survey666These are spring and autumn preferentially. By inspecting dates of ETNOs discoveries, one finds that they were discovered primarily in spring and autumn., which possibly leads to bias in $L$. This bias should suppress objects near the ecliptic with right ascension, and hence $L$, around 90 deg. This is where all ETNOs with $a>250$ au, $q>35$ au, but one are found; there is another one found around 270 deg, see Fig. 3 in Sheppard & Trujillo (2016). If the observed population is biased in this way in $L$, then the external perturbation scenario is further promoted. With $\omega$ clustered around $\omega\approx 340$ deg and $L$ biased, $\Omega$ is biased as well. The Galactic plane avoidance then suppresses ETNOs with $\Omega$ around 110, where oddly777This is only odd if there is no substantial external perturbation present. the majority of ETNOs reside, and around 290 deg. The asymmetry of the $\Omega$ distribution cannot be explained by the intrinsic declination bias, as this is symmetric with respect to the preferred value $\Omega=180$ deg (de la Fuente Marcos & de la Fuente Marcos, 2014).

Observational bias prefers low inclinations (de la Fuente Marcos & de la Fuente Marcos, 2014). Hence, the clustering in $i$ is inconsistent with being an observational effect. Interestingly, the (biased) inclination distribution of the observed ETNOs resembles the unbiased inclination distribution of the population of scattered Kuiper belt objects (KBOs), high $a$ and $e$ KBOs888ETNOs are a subset of the scattered KBOs.. We illustrate this in Fig. 3, where we show an unbiased cumulative distribution function of the scattered KBOs (dashed line) versus observed cumulative distribution of ETNOs (solid line) as a function of inclination. For scattered KBOs, we used the model from Gulbis et al. (2010), modelling the inclination distribution as a Gaussian multiplied by $\sin i$ (Brown, 2001), centred at $19.1\begin{subarray}{c}+3.9\\ -3.6\end{subarray}$ deg, and with a width of $6.9\begin{subarray}{c}+4.1\\ -2.7\end{subarray}$ deg. Cumulative distribution for ETNOs was constructed in a simplified manner. The dots in Fig. 3 correspond to the observed inclinations of ETNOs and are equidistantly spaced between 0 and 1 on the vertical axis. We conclude that ETNOs inclination distribution is probably not significantly biased. But still, we have to bear in mind the preference of low inclinations and that the width of the real ETNOs inclination distribution is somewhat higher than that inferred from Fig. 3.

In what follows, we consider statistics about the orbital elements $\omega$, $\Omega$, and $i$ of the ETNOs as not significantly biased and representative of the real (unbiased) population.

Batygin & Brown (2016) and also Sheppard & Trujillo (2016) suggested taking into account only ETNOs that are dynamically stable for the age of the solar system. This is a good idea in general. However, it is misleading to investigate the stability of ETNOs in a dynamical model with giant planets alone and no external perturbation, as carried out by Sheppard & Trujillo (2016) and Batygin & Brown (2016), and then infer properties of the external perturbation from such a “stable” sample. We do not distinguish between stable and unstable ETNOs in this way, but we do distinguish them on $q>35$ au and $q<35$ au ETNOs.

3 Milgromian dynamics

In 1983 Milgrom proposed a simple modification of the standard dynamics attributed to low accelerations (Milgrom, 1983b). He did so to explain dynamical discrepancies in the Universe manifested at that time mainly through asymptotically flat rotation curves (Bosma, 1981; Rubin et al., 1982). He assumed acceleration to be calculated not from ${\bf g}={\bf g^{N}}$ but from

[TABLE]

where ${\bf g}^{N}$ is expected Newtonian gravitational acceleration, $g^{N}$ is its norm, $\nu(y)$, $y\equiv g^{N}/a_{0}$ is a modification factor, the so-called interpolating function, fulfilling behaviour $\nu(y)\rightarrow 1$ for $y\gg 1$ (Newtonian limit) and $\nu(y)\rightarrow y^{-1/2}$ for $y\ll 1$ (deep-MD limit), and $a_{0}\sim 10^{-10}$ m s*-2* is a new constant of nature. He realised that such modification leads to a whole new phenomenology, so that low acceleration systems, such as galaxies, should display some special, systematic properties (Milgrom, 1983a, b). At that time, when the amount of data on galaxies, and especially on low surface brightness (LSB) galaxies999Systems in deep-MD regime when isolated., was limited, only one observation, which is actually a premise of writing down Eq. (1) – asymptotically flat rotation curves – was recognised as one of these special properties. Today, when various types of galaxies, including LSB galaxies, were discovered and scrutinised, we now know that many Milgrom’s predictions from his 1983 papers were fulfilled (reviewed e.g. in Famaey & McGaugh, 2012).

Milgrom and others have already developed generally covariant modified gravity theories, where Eq. (1) is implemented as a static weak field limit in the case of systems with spherically symmetrically distributed mass. In these theories the weak field regime is governed by a modified Poisson equation. Two types of the modified Poisson approach are common. In the first type, known as aquadratic Lagrangian theory (AQUAL), the modified Poisson equation is written (Bekenstein & Milgrom, 1984)

[TABLE]

where the latter equality is the classical Poisson equation with Newtonian potential $\Phi_{N}$, (baryonic) matter density $\rho$, and Newton’s gravitational constant $G$, $\mu$ is inverse of $\nu$ from Eq. (1), and the equation of motion is ${\bf g}=-\nabla\Phi$. The second modified Poisson approach, known as quasi-linear MOND (QUMOND), states that the governing gravitational potential $\Phi$ is given by (Milgrom, 2010)

[TABLE]

Having $\mu$ as inverse of $\nu$, the two theories coincide in the spherically symmetric case, giving rise to Eq. (1) (Milgrom, 2010). Outside of the spherical symmetry the two theories are not equivalent.

Eq. (1) can equivalently come from modified inertia instead of the modified gravity discussed above, where modified is the kinetic portion of the classical action, hence the equation of motion. Milgrom (1994) showed that such theories always have to be time nonlocal to be Galilean invariant. As of this complication there is no full-fledged modified inertia theory developed yet (but see discussion in section 7.10 of Famaey & McGaugh, 2012). In modified inertia approach, Eq. (1) holds exactly in the case of circular orbits alone.

3.1 Effect of Milgromian dynamics in the solar system

Three effects of MD are recognised to act in the solar system. The three effects are separable as long as $r\ll r_{M}$, where $r_{M}\equiv(GM_{\odot}/a_{0})^{1/2}$ (Milgrom, 2009a, 2012).

• •

Enhanced gravity effect

This classical MD effect results from departure of the MD interpolating function from unity. In the case of a spherical system, the anomalous acceleration due to this effect, $\delta{\bf g}={\bf g}-{\bf g}^{N}$, can be written as

[TABLE]

This effect successfully encapsulates various aspects of Galactic phenomenology (Milgrom, 1983b; Famaey & McGaugh, 2012). This effect strongly depends on the choice of the interpolating function and value of $a_{0}$, and can be made arbitrarily small by a suitable choice of the interpolating function.

• •

External field effect (EFE)

This effect is inherent to modified gravity formulations of MD101010This effect does not appear in the modified inertia formulations of MD. and it comes from the non-linearity of the modified Poisson equation. An embedding external field that is constant and uniform, with magnitude $g_{e}$, $g_{e}\lesssim a_{0}$, influences internal dynamics in an embedded system (for more details see e.g. Paučo & Klačka, 2016; Famaey & McGaugh, 2012). Contrary to the enhanced gravity effect, EFE appears even in the case of high internal accelerations ($\mu$ and $\nu$ arbitrarily close to 1). The solar system is falling freely in the external gravitational field of the Galaxy. The external field has magnitude $g_{e}\approx v^{2}_{0}/R_{0}\approx 240^{2}$ km2 s*-2* $/\leavevmode\nobreak\ 8.0$ kpc $\doteq 2.0$ $\times 10^{-10}$ m s*-2* (for values of the Galactic parameters $v_{0}$ and $R_{0}$, we refer to McMillan & Binney, 2010; Schönrich, 2012). Close to the Sun, at distances $r\ll r_{M}$, $r_{M}\equiv(GM_{\odot}/a_{0})^{1/2}\sim 7000$ au, where $r_{M}$ is characteristic radius in MD solar system, EFE expresses itself dominantly as an anomalous quadrupole, $\delta\Phi$, so a governing potential entering the equation of motion is written

[TABLE]

where ${\bf e}\equiv{\bf g}_{e}/g_{e}$ is a unitary vector pointing towards the Galactic centre111111In AQUAL, ${\bf e}$ is pointing in the direction of the real (Milgromian) gravitational field of the Galaxy. In QUMOND, ${\bf e}$ is pointing in the direction of the Newtonian gravitational field. For the approximately axially symmetric Galaxy, these directions are nearly the same. (GC) and $\Phi_{N}$ is the expected Newtonian potential (Milgrom, 2009a; Blanchet & Novak, 2011). The strength of the anomaly at a given position is represented by the magnitude of the parameter $Q_{2}$, which depends on the specific modified Poisson formulation, form of the interpolating function, and strength of the external field in natural units $a_{0}$ (Milgrom, 2009a; Blanchet & Novak, 2011). If not stated otherwise, $Q_{2}$ is assumed to be positive. All routinely used forms of the interpolating function yield121212In Milgrom (2009a), a dimensionless parameter $q$, $q=-2Q_{2}(GM_{\odot})^{1/2}/(3a^{3/2}_{0})$, was employed instead of $Q_{2}$. $Q_{2}>0$ (Milgrom, 2009a; Blanchet & Novak, 2011). If $Q_{2}$ is taken to be a constant in some region of interest, the value $\delta\Phi$ has the same functional form as the tidal potential coming from the hypothetical P9 or that from the Galaxy. We do not know the exact form of the interpolating function, therefore we do not know the value of $Q_{2}$, but to be dynamically important in the ETNOs region it would have to be much stronger than the quadrupole coming from the Galactic tides. If we consider values of $Q_{2}$ close to the the theoretical maximum allowed by the Earth-Cassini spacecraft range data (hereafter Cassini data), $Q_{2}\sim 10^{-27}$ s*-2* (Hees et al., 2014), then EFE is $\sim 10^{3}$ times stronger than the Galactic tides (Hees et al., 2014). Instantaneous action of EFE is the same as that of P9 lying in the direction of the GC-anticentre, characterised by the tidal parameter $GM^{\prime}/r^{\prime 3}=Q_{2}/3$, where $M^{\prime}$ is mass of P9 and $r^{\prime}$, $r^{\prime}\gg r$ is its heliocentric distance (Iorio, 2010). Thus, EFE and the dynamical effect of P9 are indistinguishable on short timescales in the ETNOs region. Intriguingly, perihelia of ETNOs with $a$ beyond 250 au cluster close to the direction of the present-day Galactic anticentre (the preferred direction in Eq. (5)). Comparison between orbital characteristics of ETNOs induced on long timescales by EFE and those induced by P9 is a topic of this study.

• •

Asphericity effect

This effect appears only in modified gravity formulations of MD. It affects every system with non-spherically symmetric matter distribution. Suppose we have an isolated high-acceleration system S, which assumes $\mu=\nu=1$ to the desired accuracy everywhere within S, of total mass $M$ and extension $R$, $R\ll r_{M}$, where $r_{M}=\sqrt{GM/a_{0}}$. Suppose also that mass distribution within S, $\rho({\bf r})$, varies on timescales much larger than $r_{M}/c$, where $c$ is the speed of light. The classical Poisson equation coincides with MD equations in S. But, the Poisson equation is not valid everywhere up to infinity. It is definitely not valid for $r>r_{M}$. Hence, the solution of the Poisson equation in S is slightly different from the Newtonian solution. In QUMOND, Milgrom (2012) showed that dynamics in S is governed by the potential

[TABLE]

where $\Phi_{N}$ is the expected Newtonian potential, $\delta\Phi$ is a dominant MD correction from the asphericity effect, a quadrupole field, and $\alpha$, typically $\alpha\sim 1$, is a numerical factor that depends on the form of the interpolating function. In the case of the solar system131313Solar system is not an isolated system, so we have a superposition of EFE and the asphericity effect., $M\approx M_{\odot}$ and $Q_{ij}$ is the quadrupole moment of the solar system.

In the standard formulation of MD and in many parent relativistic theories of MD, this effect is too weak to be detected in the solar system (Milgrom, 2012). In special theories where Newtonian regime is not reached beyond $g^{N}\sim a_{0}$, but beyond $g^{N}\sim\kappa a_{0}$, $\alpha$ in Eq. (3.1) becomes enhanced by $\kappa^{2}$ approximately (Milgrom, 2012). An example of such a theory is, for instance TeVeS (Bekenstein, 2004), where $\kappa$ must have very high value to become compatible with general relativity (Milgrom, 2012). However, the existence of two different acceleration constants in the theory, $a_{0}$ and $a^{*}_{0}\equiv\kappa a_{0}$, seems rather unnatural. The anomalous perihelion precession of Saturn constrains $\kappa$ to be $\kappa\leq 3500$ (Iorio, 2013).

4 Benchmark model

We start with construction of a simple yet realistic dynamical model of the solar system in MD. We call this model a benchmark model. The region of our interest is a sphere of radius 2000 au that is centred on the Sun. This region is large enough to investigate the orbit of ETNO in its entirety and also small enough for Eq. (5) to be applicable; Sedna, a body with the largest aphelion distance, cruises up to $\sim 900$ au.

By combining an analysis of rotation curves of galaxies and solar system constraints141414The solar system constraints in Hees et al. (2016) come from EFE quadrupole anomaly, all the other MD effects are negligible in the planetary zone (Hees et al., 2014, 2016). provided by the Cassini spacecraft, Hees et al. (2016) showed that many popular MD interpolating functions are not allowed by the data. These authors also listed examples of the allowed interpolating functions (their Table 2). Generally speaking, these allowed functions are characterised by extremely rapid transition into the Newtonian regime.151515Also, these are the functions for which the enhanced gravity effect is negligible in the planetary region, compared to EFE. Referring to Table 2 in Hees et al. (2016), we chose $\overline{\nu}_{n}(y)=\left[1-\exp(-y^{n})\right]^{-1/(2n)}+\left[1-1/(2n)\right]\exp(-y^{n})$ interpolating function family as most representative. The Cassini data restrict $\overline{\nu}_{n}$ to $n\geq 2$. The enhanced gravity effect is completely negligible in the whole solar system161616Except for small regions, where internal and external gravitational fields add up into a vector of small magnitude. for $\overline{\nu}_{n}$, $n\gtrsim 3$. Thus, in the first approximation, we neglect the enhanced gravity effect in the benchmark model.

The parameter $Q_{2}$ of Eq. (5) is position dependent in general (Milgrom, 2009a). However, as shown by Blanchet & Novak (2011) in the framework of AQUAL, $Q_{2}$ can be considered a constant in a large sphere centring on the Sun. The variation of $Q_{2}$ within our volume of interest is up to few percent (Blanchet & Novak, 2011).

In the benchmark model, the only MD effect we take into account is EFE quadrupole anomaly. We model EFE as in Eq. (5), with position independent $Q_{2}$, $Q_{2}=Q_{2}(0)$. All higher multipoles are omitted. The value of $Q_{2}$ was constrained by the Cassini data to be $Q_{2}=3\pm 3\times 10^{-27}$ s*-2* (nominal $\pm 1\sigma$ value) (Hees et al., 2014). From now on, if we refer to $Q_{2}$ without explicitly stating its units, the unit is $10^{-27}$ s*-2*. The unit vector ${\bf e}$ of Eq. (5), pointing towards the GC, is assumed to circulate uniformly with time according to

[TABLE]

where the first three lines are components of ${\bf e}$ in the Galactic coordinates171717Principal axis points from the Sun to the GC. and the last three lines represent transformation from the Galactic to the heliocentric ecliptic coordinates at J2000. We assumed that the Sun is moving in the circular orbit around the GC181818The Sun is moving in the clockwise direction, looking from the north Galactic pole perspective., lying in the Galactic mid-plane, with angular speed $w$; $\tau$ is the phase of the motion.

5 Orbital clustering: Initial numerical exploration

Our aim is to investigate whether orbits of primordial ETNOs are driven to long-living orbital confinement in the solar system ruled by the dynamical laws of MD. For this purpose we made use of the benchmark model. Additionally, we model the secular effect of the giant planets by considering the Sun of radius $R=a_{U}$, where $a_{U}$ is the semimajor axis of Uranus and $J_{2}$ is the moment of magnitude $J_{2}=\sum_{i}m_{i}a^{2}_{i}/(2MR^{2}),$ and where we sum over the giant planets with masses $m_{i}$, semimajor axes $a_{i}$, and $M$ is mass of the Sun (Batygin & Brown, 2016). Close encounters with the giant planets were completely omitted in this initial numerical exploration.

We traced the evolution of primordial ETNOs in $\omega-e$, $\Omega-e$, and $L-e$ plane, considering $Q_{2}$ taken from a set $Q_{2}=\{1.0,\leavevmode\nobreak\ 2.0,\leavevmode\nobreak\ 4.0,\leavevmode\nobreak\ 6.0\}$. In MD, the external perturbation, substituting a distant planet, is the quadrupole anomaly arising from EFE; see Eqs. (5) and (4). The direction of the perturbation is known and we alter only the parameter $Q_{2}$.

The ETNOs were modelled as test particles and followed for 4 Gyr. We consider three values of initial semimajor axis: $a=150$, 350, and 550 au. For each value of $a$, 20 test particle orbits were initialised with an $e$ spanning range (0, 0.95) with a constant increment of 0.05 191919In the case $a=150$ au, the particles in $e=0.90$ and $e=0.95$ orbits were doomed as $q<a_{U}$ holds for them at the time of their initialisation., $i$ drawn randomly from a uniform distribution in range (0, 5) deg, $\omega$ drawn randomly from a uniform distribution in range (0, 360) deg, and $\Omega$ calculated so that $L=90$ deg. Similarly, another 20 particles with $L=270$ deg were set up for each value of $a$. The particles were started in perihelion. The values $L=90$ and 270 deg were identified as stable libration regions with trial and error method. We note that $L=270$ – 90 deg is also roughly the direction of the present-day GC–anticentre. The ecliptic latitude of the GC is only minus few degrees. Here our intention is to demonstrate the existence of the long-living libration regions; in the next section we examine development of the orbital clustering, starting with completely randomised orbits.

Every orbital integration in this paper was carried out using mercury6 N-body integration software package (Chambers, 1999). The hybrid symplectic-Bulirsch-Stoer algorithm with a timestep equal to a tenth of the orbital period of Neptune was employed. Phase $\tau$ in Eq. (4) was chosen in a way that at $t=4$ Gyr the vector ${\bf e}$ points towards the present-day GC. We tested that the value of $\tau$ actually had no influence on the results. If a particle had reached the distance 2000 au from the Sun, it was immediately discarded from the simulation202020In fact, none of the particles experienced such large change in $a$. The discarded particles typically came down to $r=a_{U}$.. There is no reason to discard particles beyond 2000 au in the Newtonian simulation, but we wanted to apply the same measure in both frameworks. If a particle had come closer than $r=a_{U}$, it was discarded from the simulation as well.

For the sake of comparison, we ran also two P9 simulations in Newtonian dynamics. In the first simulation, P9 of mass $M^{\prime}=10\leavevmode\nobreak\ M_{\oplus}$ was assumed to orbit in $a^{\prime}=700$ au, $e^{\prime}=0.6$, $i^{\prime}=0$ deg, and $L^{\prime}=240$ deg orbit (Batygin & Brown, 2016). The primed quantities always refer to P9. The test particle orbits were initialised with $L-L^{\prime}$ equal to 0 (20 particles for each $a$) and 180 deg (20 particles for each $a$), and nil inclinations212121In this case, it is meaningful to examine evolution in the $L-e$ plane only, as all inclinations were set to zero.. All the other orbital elements were set up as in the MD simulation. P9 was initialised in aphelion. We also did a simulation assuming $i^{\prime}=30$ deg and inclined test particle orbits with $i$ drawn randomly from a uniform distribution in range (0, 5) deg. In addition to the previous rules for discarding particles, we also discarded all particles that come within one Hill radius from P9.

We depict evolution of test particles in $L-e$ plane as inferred from the first P9 simulation in Fig. 4. All particles meet the condition $i-i^{\prime}=0$ by design. The discussed orbital elements always refer to the heliocentric ecliptic J2000 reference frame. Two transparency levels are used as a proxy for two different time windows of the simulation. The less transparent evolutionary tracks represent long-living particles in the last Gyr of the simulation. The more transparent tracks refer to evolution in the first 3 Gyr of the simulation and also comprise unstable orbits. Orbits with $a=150$ au typically experience trivial apsidal circulation. Only one small libration island exists at $L\approx 240$ deg (aligned with P9) for near-circular orbits. For $a=350$ au, two long-term stable libration regimes exist (Batygin & Brown, 2016): first, the resonant dynamics regime concerning eccentric orbits (with eccentricities similar to those of the observed ETNOs) and centred at $L\approx 60$ deg and anti-aligned with P9 and, second, the secular dynamics regime concerning less eccentric orbits, centring at $L\approx 240$ deg and aligned with P9. These absolute positions come from our choice of $L^{\prime}$, which was tuned to match the observed clustering position; see Fig. 1. Batygin & Brown (2016) always refer to the relative longitude of perihelion $\varpi-\varpi^{\prime}$ in their analysis, as they did not know the orbital elements of the hypothesised P9 a priori. Here, we build on the results of their pioneering numerical experiments. The shift of the libration centres in $L-e$ plane with time follows from the shift of $L^{\prime}$ caused by the giant planets. At 550 au level, no stable libration develops. In the second P9 test case, where we consider an inclined perturber with $i^{\prime}=30$ deg and inclined test particle orbits, there were no long-term stable orbits.

The picture in MD is rich as well. Orbits with $a=150$ au are characterised by trivial circulation of the angles $(L,\leavevmode\nobreak\ \Omega,\leavevmode\nobreak\ \omega)$. No stable libration develops at this $a$-level for realistically high value of $Q_{2}$. The values $a=350$ and 550 au are more interesting. The orbits of particles are strongly confined in longitude of ascending node $\Omega$ in all 4 MD simulations. Both eccentric (ETNOs) and less eccentric particles, with both $a=350$ and 550 au, strongly prefer $\Omega\in(180,360)$ deg. In ETNOs eccentricity range (EER), $e\in(0.69,0.93)$, stable, and often tight, libration around $\Omega\approx 270$ deg occur, for instance where $Q_{2}=1$ or 2, $a=350$ au.

Concerning perihelion longitude $L$, for $a=350$ au, we get a circulation of the angles in EER irrespective of $Q_{2}$. In the case of lower eccentricities, we get a circulation for models with $Q_{2}\lesssim 4$, and neither circulation nor strong confinement using models with $Q_{2}\gtrsim 4$. The specialty of values $L=90$ and 270 deg is still noticeable.

The case $Q_{2}=1$, $a=550$ au is especially interesting. The motion in $L-e$ plane is confined to narrow regions near $L\approx 90$ and $L\approx 270$ deg in a broad range of eccentricities, spanning from low eccentricities to EER. For eccentricities $e\gtrsim 0.9$, the circulation of the angles sets in. In $\omega-e$ plane, there are orbits clustered near $\omega\approx 0$ and $\omega\approx 180$ deg, below and at the low eccentricity end of EER. We depict the evolution of test particles in $\omega-e$, $\Omega-e$, and $L-e$ plane as inferred from the benchmark model simulations using $Q_{2}=1$ in Fig. 5.

Despite the rich structure in $\omega-e$, $\Omega-e$, and $L-e$ planes in the MD, the simulated evolutionary paths are in sharp contrast with present-day observations. The simulations predict that the vast majority of ETNOs should have $\Omega\in(180,360)$ deg. Hence, the vast majority of ETNOs should have ascending node lying in the same half-space, with the GC lying in the same half-space. This is not an observed trend. Of 19 ETNOs, 14 have $\Omega\in(0,180)$ deg, while 3 of them have $\Omega$ very close to 180 deg (Sheppard & Trujillo, 2016). Moreover, the tightness of the observed clustering in $\omega$ is not reproduced in our simple MD models. The discrimination of the orbits of particles according to their $q$ is of no meaning here as the giant planets were omitted in these initial simulations. In order to compare MD models with Newtonian models including P9, we should switch to a more realistic approach, also taking the gravitational interaction between ETNOs and Neptune into account.

5.1 The case of negative $Q_{2}$

The fact that the MD models predict the clustering region that lies roughly 180 deg ahead of the observed clustering in $\Omega$ motivates us to question the direction of the perturbation coming from EFE (effectively the sign of $Q_{2}$). When $Q_{2}$ is positive, a test particle is repelled from the Sun in the direction of the external field, which means that effectively the gravity of the Sun is reduced in that direction, while the particle is attracted to the Sun in the direction perpendicular to the external field, i.e. effectively the gravity of the Sun is enhanced in that direction. In Milgrom (2009a), the author studied the effect MD has in the planetary region of the solar system. In the last section, he states that222222Milgrom used $q$, $q\equiv-2Q_{2}(GM_{\odot})^{1/2}/(3a^{3/2}_{0})$, instead of $Q_{2}$, to quantify the strength of EFE. although $Q_{2}$ was positive for all interpolating functions $\mu$ he studied, he was “not able to prove that this is always the case from basic properties of $\mu$” (Milgrom, 2009a). Hees et al. (2014) showed that the Cassini data yield $Q_{2}=3\pm 3$ (nominal $\pm 1\sigma$ value).

We substituted $Q_{2}\rightarrow-Q_{2}$ and repeated the benchmark model simulations. This means we considered $Q_{2}$s that are 1 - 3$\sigma$ off of the nominal value of Hees et al. (2014), with corresponding probabilities for $Q_{2}$ ranging from 16 ($Q_{2}<0$) down to 0.1$\%$ ($Q_{2}<-6$).

At $a=350$ au level, the substitution causes a shift of the clustering region by roughly 180 deg in $\Omega$, while leaving evolution in $\omega-e$ and $L-e$ planes without qualitative change. This can be seen in the case $Q_{2}=-1$ in Fig. 6. Looking at EER at $a=550$ au, $\Omega$ lies sometimes preferentially in range $(0,180)$ deg ($Q_{2}=-1$), sometimes in range $(180,360)$ deg ($Q_{2}=-4$ and -6), and sometimes there is preference for neither of the two regions ($Q_{2}=-2$).

6 Add Neptune and stir

Neptune may have influenced the distribution of the orbital elements of ETNOs through orbital filtering (e.g. by preferentially ejecting ETNOs in specific orbits) and resonant effects. As a next step, we added Neptune orbiting the Sun to the dynamical model. We consider the gravitation of the Sun and planet Neptune in the usual Newtonian way as a direct N-body interaction. The secular effect of the three remaining giant planets is modelled by considering the Sun of radius $R=a_{U}$, where $a_{U}$ is the semimajor axis of Uranus, and $J_{2}$ moment of magnitude $J_{2}=\sum_{i}m_{i}a^{2}_{i}/(2MR^{2})$, where we sum over the remaining giant planets with masses $m_{i}$, semimajor axes $a_{i}$, and $M$ is mass of the Sun (Batygin & Brown, 2016).

We performed several numerical simulations, aiming to compare the dynamical effects of P9 and MD on orbits of primordial TNOs. Each simulation started with a thin disk of 400 test particles, where 40 particles had semimajor axis $a=100$ au, 40 particles $a=150$ au, and so on up to 40 particles with $a=550$ au. Perihelion distances of the particles with a given $a$ covered range $(30,50)$ au with a constant increment. The inclination, argument of perihelion, and longitude of ascending node of the test particle orbits were randomly drawn from a uniform distribution in range $(0,5)$ deg $(i)$, and $(0,360)$ deg $(\omega,\Omega)$. The particles were initialised at their perihelia. We also ran all the simulations with twice as many particles with an additional 400 particles initialised by the same procedure as delineated above to test the effect of the test particle number.

In the benchmark model simulations, we considered various values of $Q_{2}$ taken from a set $Q_{2}=\{0.2,\leavevmode\nobreak\ 0.5,\leavevmode\nobreak\ 1.0,\leavevmode\nobreak\ 2.0,\leavevmode\nobreak\ 4.0,\leavevmode\nobreak\ 6.0\}$. For the sake of comparison, simulations with P9 in Newtonian dynamics (no EFE $\Leftrightarrow$ $Q_{2}=0$) were carried out. Two initial P9 orbits were considered: (#1) $a^{\prime}=700$ au, $e^{\prime}=0.6$, $i^{\prime}=30$ deg, $\omega^{\prime}=150$ deg, $\Omega^{\prime}=113$ deg (Batygin & Brown, 2016; Brown & Batygin, 2016; Fienga et al., 2016); (#2) the same232323Starting values of $a^{\prime}$, $e^{\prime}$, and $i^{\prime}$ were altered only negligibly during the time of the simulation. $a^{\prime}$, $e^{\prime}$, $i^{\prime}$ as for orbit (#1), but $\omega^{\prime}=100$ deg, $\Omega^{\prime}=140$ deg. The orbit (#2) orbital elements account for variation in $\omega^{\prime}$ and $\Omega^{\prime}$, caused by the giant planets, and at $t=4$ Gyr the final $\omega^{\prime}$ and $\Omega^{\prime}$ were very close to the orbit (#1) starting values.

6.1 Sednoids

The discovery of Sedna (Brown et al., 2004) and 2012 VP113 (Trujillo & Sheppard, 2014) was surprising as their orbits did not fit into the standard picture of the solar system, settled in the current Galactic environment, containing a set of known planetary bodies. The literature discussed scenarios in which these orbits were shaped in environments that, in the past, were different from the present-day environment, for instance the possibility that the Sun was born in a star cluster (Kaib & Quinn, 2008; Brasser et al., 2012) and/or was traveling much closer to the GC in the past (Kaib et al., 2011).

Another possibility is that the torque necessary to elevate perihelia of these bodies comes from an unknown planet-mass body, or maybe multiple bodies, beyond Neptune (Gomes et al., 2006; Batygin & Brown, 2016). In MD, similar orbits could be induced by EFE with no need for an additional external body (Paučo & Klačka, 2016). Both options have in common that they operate continuously. Hence, they produce a “live” population in which perihelia are repeatedly, and slowly, elevated and depressed. This is in contrast to a “fossil” population induced by stellar encounters (Morbidelli & Levison, 2004), where $q$ changes are almost instantaneous, or capture of these bodies from a passing star (Jílková et al., 2015). An excess of classical and large semimajor axis Centaurs, which are icy bodies with $q$ among the giant planets and $a>100$ au, would be natural in the continuously disturbed population (Gomes et al., 2015). The influx of comets into the inner solar system would be altered as well (Gomes et al., 2006; Paučo & Klačka, 2016).

Evolution in $a$–$e$ plane is depicted during the last Gyr of the benchmark model simulation in Fig. 7 for various values of $Q_{2}$. Figures showing $a$–$e$ diagrams are always based on the simulations starting with 800 particles. With $Q_{2}=0$ and no P9, the particles would be stuck between the two rightmost/bottom-most dashed lines ($q$ in range between 30 and 50 au). Moreover, assuming $Q_{2}=0$ and no P9, inclinations of the particles remain low, typically $i<10$ deg.

From Fig. 7 we can conclude that the orbit of Sedna constrains $Q_{2}$ to $Q_{2}>0.5$. The orbit of 2012 VP113, with similar perihelion distance but much lower semimajor axis than Sedna, requires $Q_{2}>4$, if we demand that EFE is solely responsible for the elevation of perihelia. Even with $Q_{2}=4$ – 6, 2012 VP113 still lies at the low-end border of the achievable eccentricities at $a\approx 250$ au. Objects with higher eccentricities are easier to detect. The orbit of 2012 VP113 is thus only marginally compatible with the EFE origin scenario.

6.1.1 Beyond the benchmark model

What effect can be expected from the next lowest-order correction (octupole) neglected in Eq. (5) and in the benchmark model? This correction (adds to RHS of Eq. (5)) can be written as242424Blanchet & Novak (2011) state an opposite sign in the RHS of Eq. (8), as they have been using $\nabla^{2}\Phi_{N}=-4\pi G\rho$ and ${\bf g}=\nabla\Phi$. We used $\nabla^{2}\Phi_{N}=4\pi G\rho$ and ${\bf g}=-\nabla\Phi$ throughout the paper, hence the difference. (Blanchet & Novak, 2011)

[TABLE]

The octupole strength parameter $Q_{3}$ (typically $Q_{3}<0$) depends on the specific modified Poisson formulation, form of the interpolating function, and strength of the external field in natural units $a_{0}$.  The parameter $Q_{3}$ is position dependent in general. The variation of $Q_{3}$ within our volume of interest is up to several percent (Blanchet & Novak, 2011). We treat $Q_{3}$ as position independent, $Q_{3}=Q_{3}(0)$.

Blanchet & Novak (2011) tested some of the commonly used interpolating functions $\mu,$ for instance the $\mu_{n}(x)=x/(1+x^{n})^{1/n}$ family, in the framework of AQUAL, and found that $|Q_{3}|<1.2\times 10^{-40}$ m*-1* s*-2* close to the Sun. From now on, when we refer to $Q_{3}$ without explicitly stating its units, the unit is $10^{-40}$ m*-1* s*-2*. The highest $|Q_{3}|$ was found for the simple interpolating function $\mu_{1}(x)=x/(1+x)$. It is known that $\mu_{1}$, and its inverse $\nu_{1}$, are excluded by the solar system tests (Blanchet & Novak, 2011; Hees et al., 2016) because they give a $Q_{2}$ that is too large. The higher $n$ functions from the $\mu_{n}$ family, allowed by the solar system data (Blanchet & Novak, 2011; Hees et al., 2016), yield typically $|Q_{3}|\sim 10^{-2}$ (Blanchet & Novak, 2011). We found that with such a low value of $Q_{3}$, the incorporation of the octupole term, Eq. (8), has a negligible effect. We note that the most stringent constraints on $Q_{2}$ comes from Hees et al. (2016), where QUMOND framework was applied. The values of $Q_{3}$ that we refer to come from Blanchet & Novak (2011) and were calculated in AQUAL. Thus, when we link the two we assume that their values are similar in both modified gravity frameworks, for a given pair of interpolating functions $\mu$ and $\nu$, and a given external field strength. Approximation of this kind is common; see e.g. Milgrom (2009a).

Switching back to $\overline{\nu}_{n}$ family and QUMOND, $\overline{\nu}_{2}$ yields $Q_{2}=5.9$ when one assumes $\eta\equiv g_{e}/a_{0}=2.3$, $a_{0}=0.815\times 10^{-10}$ m s*-2* (Hees et al., 2016). The quoted value of $a_{0}$ comes from rotation curves fits using $\overline{\nu}_{2}$. With $\overline{\nu}_{2}$ one gets, besides EFE, as well as a non-negligible enhanced gravity effect. We incorporated the enhanced gravity effect, Eq. (4), into the model and repeated the simulation. We used $\overline{\nu}_{2}(y)$, $y=|{\bf g}^{N}_{\odot}+{\bf g}^{N}_{e}|/a_{0}$, where ${\bf g}^{N}_{\odot}$ is Newtonian gravitation from the Sun, and ${\bf g}^{N}_{e}$ was found from ${\bf g}_{e}=\overline{\nu}_{2}(|{\bf g}^{N}_{e}|/a_{0}){\bf g}^{N}_{e}$ knowing ${\bf g}_{e}$. Direction of fields ${\bf g}_{e}$ and ${\bf g}^{N}_{e}$ rotate as prescribed by Eq. (4). An outcome is in Fig. 8. 2012 VP113 still lies at the low-end border of the achievable eccentricities at $a\approx 250$ au. We remind that $\overline{\nu}_{n}$, $n<2$, is disfavoured by the Cassini data (Hees et al., 2016), and for $\overline{\nu}_{n}$, $n\gtrsim 3$, the enhanced gravity effect is completely negligible in the ETNOs region.

We can now directly compare the MD and P9 scenarios. The last Gyr of the Newtonian simulation $(Q_{2}=0)$ with P9 is shown in Fig. 9 for P9 starting orbit (#1). Starting in orbit (#2), P9 yields very similar picture. Apparently, P9 accounts well for the existence of the two. Although orbits close to Sednoids in terms of $a$, $e$, and $i$ were not induced (or not preserved) in the simulation, probably because of the low number of test particles, there were orbits induced with lower semimajor axes and higher perihelion distances than those of Sednoids, which means that the effect of P9 should be sufficient to shape orbits of Sednoids. This is not a new result, P9 was designed to do this.

We also tested how Figs. 7 and 8 change with longer integration time of 4.5 Gyr and found that the differences are minute. Additionally, we varied also phase $\tau$ and angular speed of the Sun $w$, present in Eq. (4). We used the fact that $w$ was found to lie in bounds $29.8-31.5$ km s*-1* kpc*-1* (McMillan & Binney, 2010). It reveals that the effect of these variations is of little importance too.

6.1.2 The case of negative $Q_{2}$

We show $a-e$ diagrams that are analogous to the previous diagrams, but now with negative values of $Q_{2}$, $Q_{2}=-1,$ and $Q_{2}=-4$, in Fig. 10. The orbits of both Sednoids are produced in the MD model as long as $|Q_{2}|>4$. With $Q_{2}=-4$, 2012 VP113 lies at the low-end border of the achievable eccentricities at $a\approx 250$ au.

6.2 Orbital clustering

Can statistics on ETNOs orbits be understood with the aid of MD without additional planets in the solar system? At first, we analysed our simulations without explicit visualisation. We performed the Kolmogorov-Smirnov (KS) test252525To be clear in statistical terms and notation, we give a following example: the probability $p$-value being equal to 0.05 means that the null hypothesis is not rejected at $5\%$ level; in our case the null hypothesis is always the hypothesis in which the two compared datasets have the same distribution. comparing distributions of $\omega$, $\Omega$, $L$, and $i$ of the observational data with those inferred from the benchmark model simulation (MD model test) at the time instant $t=4$ Gyr, considering various values of $Q_{2}$. Additionally, we also did the KS test based on results from $Q_{2}=0$ (no effects of MD in the solar system) simulation with (P9 test) and without (control test) P9. Only test particle orbits fulfilling the criteria – $i<50$ deg, $q<80$ au at $t=4$ Gyr, were considered to account for the observational preference of such orbits. The orbits are discerned according to their semimajor axis and divided into two groups: one with $150<a<250$ and the other with $a>250$ au. The results are presented in Table 1 for simulations starting with 400 particles.

As there was only one test particle with $a>250$ au left in the P9 simulation at $t=4$ Gyr (for both considered initial orbits of P9), $p$-values, $p(>250\leavevmode\nobreak\ |>250)$, are not reported in the P9 case. The P9 test shows that an incorporation of P9 in the Newtonian model increases its likelihood for $a<250$ au, though the compared samples are small ($\sim$10). A possible caveat of the P9 scenario was revealed by Lawler et al. (2017). They performed a similar simulation starting with many more test particles $(10^{5})$, interacting gravitationally with five planets (known giant planets + P9) in a full N-body manner for 4 Gyr, leading to no angular clustering at the end of the simulation.

In the semimajor axis range $150<a<250$ au, the MD model test yields $p$-values, $p(150-250\leavevmode\nobreak\ |\leavevmode\nobreak\ 150-250)$, similar to the control test, except for the distribution of $\Omega$ when even lower $p$-values are reported. This means some clustering in $\Omega$ in the simulated data in MD but at different positions than the observed data.

At the $a>250$ au level, the observed distribution of $i$ might be compatible with EFE if $Q_{2}\gtrsim 4$. However, the compatibility between the observational data and the MD model remains low in the case of $\omega$ and $L$, though $p$-values, $p(>250\leavevmode\nobreak\ |>250)$, are slightly higher than in the control test if, e.g. $Q_{2}=4$. The distribution of $\Omega$ seems to be the most problematic distribution, with the simulated data clustered but having significant offset from the observed clustering position.

We tested how the results in Table 1 change when we restrict the orbits of particles to those with $q>35$ au. An example is given in Table 2. The difference is small. There is prevalence of $q>35$ au orbits at $t=4$ Gyr. Doubling the number of test particles at $t=0$ does not change the results qualitatively. An example can be seen in Table 3. Even after doubling the initial number of particles, there was only one particle with $a>250$ au left at the end of the P9 simulation for both considered starting orbits of P9.

We now investigate the simulated data visually. Figures 11 - 12 show evolution of test particle orbits in $\omega-a$, $\Omega-a$, and $L-a$ plane in the last Gyr of the Newtonian simulation with P9. We initialised P9 in the orbit (#1). Initialising P9 in the orbit (#2) would have yielded a similar picture. Evolutionary tracks of the particles are shown every time their orbits met the criteria, i.e. $i<50$ deg, $q<80$ au. The evolutionary tracks continue down to $a=100$ au, but we do not show them below $a=200$ au as they are just a forest of full-width horizontal lines. Clustering in $\omega$ emerges only for $a\gtrsim 500$ au, as could be seen already in Batygin & Brown (2016). Potentially unstable particles, discarded at some time between $t=3$ and 4 Gyr, have $\omega$ preferentially close to 0 (360) or 180 deg at the time of their destruction, which is another feature already noticed by Brown & Batygin (2016). The relative number of particles with $\omega\approx 0$ to the number of particles with $\omega\approx 180$ can be influenced by altering an initial orbit of Planet 9 (Brown & Batygin, 2016). The best agreement between observations and the P9 model can be seen in terms of $L$. However, the clustering in orbital angles is always traced only by the potentially unstable particles.

We turn our attention to MD. Figure 13 shows the evolution of test particle orbits in $\omega-a$ plane in the last Gyr of the benchmark model simulation assuming $Q_{2}=1$. Again, only particles surviving at least 3 Gyr, orbiting in the observationally preferred orbits, are taken into account. Similarly with the P9 picture, clustering in $\omega$ develops only for large $a$ orbits, in the case of $Q_{2}=1$, $a>600$ au, and it is only traced by the potentially unstable orbits. The value of $\omega$ near 0 (360) as well as near 180 deg is preferred by the test particles. With $Q_{2}\gtrsim 2$, the clustering in $\omega$ fades away.

The test particle orbits strongly cluster around $\Omega\approx 270$ deg. This is illustrated in the left panel of Fig. 14, where we show results from $Q_{2}=2$ simulation. Strong confinement in $\Omega$ develops for stable orbits in a wide range of $Q_{2}$ values, $Q_{2}\in(0.5,6)$, and also for potentially unstable orbits if $Q_{2}\gtrsim 2$.

The distribution of $L$ shows a coherent picture for $Q_{2}$ in the range $Q_{2}\in(1,6)$. Clustering in $L$ is noticeable around $L=0$ (360) and $L=180$ deg and is traced again mainly by potentially unstable orbits. In the right panel of Fig. 14, we show the evolution of test particle orbits in $L-a$ plane, as taken from $Q_{2}=2$ simulation. As $Q_{2}$ increases, clustering develops at lower and lower semimajor axes. For example, for $Q_{2}=1$, the clustering develops only beyond $a=500$ au, while the border is at $a=300$ au for $Q_{2}=6$.

We depict histogram distributions of $\omega$, $\Omega$, $L$, and $i$ of test particle orbits at the end of the simulation, i.e. at $t=4$ Gyr, in Figs. 15 ($Q_{2}=1$) and 16 ($Q_{2}=4$). Only particles with $a>150$ au and $q>35$ au are considered. Stable particles show no apparent clustering in terms of $\omega$. Strong clustering emerges in terms of $\Omega$ across wide range of $Q_{2}$ values, $Q_{2}\in(0.5,6)$, as majority of the orbits prefer $\Omega\in(180,360)$ deg at $t=4$ Gyr. This is the trend already identified in the simple numerical exploration of Sec. 5: with positive $Q_{2}$, we get a simulated clustering position shifted by 180 deg with respect to the observed clustering position.

Stable orbits also slightly prefer $L$ around 90 and 270 deg if $Q_{2}\approx 4$. This distinction between stable and unstable orbits could be linked to the distinction between observed $q>35$ au and $q<35$ au orbits, looming for $a>150$ objects (Sheppard & Trujillo, 2016). By increasing the value of $Q_{2}$, one gets rid of low-inclination ETNOs gradually. Objects with higher inclinations are more resilient. With $Q_{2}=4-6$, majority of the particles in $a>250$ au orbits have $i\in(10,30)$ deg at $t=4$ Gyr, as observed. But the simulated inclination distribution still does not account for the lack of the observed ETNOs with inclinations lower than 10 deg, bearing in mind their observational preference.

Perihelion distance as a function of time for four test particle orbits from the benchmark model simulation, which are characterised by a high amplitude of migration in $q$, is depicted in Fig. 17. This example is relevant also for a more complicated dynamical model, where both EFE and a distant massive planet that is placed on its orbit with the aid of EFE act on ETNOs. For example, EFE with $Q_{2}=0.2$ could lift the perihelion of P9 from the edge of the planetary region262626We did not examine how migration through the giant planets’ region, whence P9 would probably originate from, proceeds. The drawback can be the slowness of the EFE induced migration, as P9 in an eccentric orbit crossing orbits of the giant planets, is endangered by their energetic kicks. to few hundred au, while with P9 in the orbit (#1) and EFE of strength $Q_{2}=0.2$, we get similar picture of the ETNOs region as in Figs. 11 - 12 with pure P9 model. We do not develop these ideas further in this paper.

6.2.1 Beyond the benchmark model

We tested the robustness of the $a-e$ diagrams inferred from the benchmark model simulation in Sec. 6.1.1. We have carried out similar robustness tests again, paying attention to the orbital clustering. These tests include the incorporation of the octupole term, Eq. (8), in the governing potential, assuming $Q_{3}\sim 10^{-2}$; the inclusion of the enhanced gravity effect in the case272727We compared this with the benchmark model $Q_{2}=6$ simulation. $Q_{2}=5.9$; the variation of $w$ and $\tau$ in Eq. (4); and the enhancement of the simulation time to 4.5 Gyr. None of these changes of the model matters; the results remain qualitatively the same.

6.2.2 The case of negative $Q_{2}$

We carried out the KS test that compares distributions of $\omega$, $\Omega$, $L$, and $i$, of simulated and observed ETNOs. The $p$-values of the test are listed in Table 4. It is analogous to Table 1, but now, MD models with $Q_{2}<0$ are investigated. We restricted ourselves to $|Q_{2}|<6$ (within $3\sigma$ range given by the Casssini data). In the case of $\omega$, $L$, and $i$ distributions, the reported $p$-values are similar to their $Q_{2}>0$ counterparts as long as models with the same $|Q_{2}|$ are compared. In the case of $\Omega$, the $p$-values are much higher.

Evolution in $\Omega-a$ and $L-a$ plane, inferred from the last Gyr of $Q_{2}=-2$ simulation, is visualised in Fig. 18. Varying $Q_{2}$ in the range $Q_{2}\in(-6,-1)$ yields a qualitatively similar picture as in Fig. 18. Potentially unstable orbits cluster around $L\approx 90$ and 270 deg. The clustering in $L$ sets in for $Q_{2}=-1$ beyond $a=450$ au, while for $Q_{2}=-6$ already beyond $a=300$ au. We found no apparent clustering in $\omega$ in our simulations with $Q_{2}<0$.

Histogram distributions of $\omega$, $\Omega$, $L$, and $i$ inferred at $t=4$ Gyr from the $Q_{2}=-1$ simulation are shown in Fig. 19. The clustering in $\Omega$ is now aligned with the observed clustering position and traced by both stable and potentially unstable orbits. With $|Q_{2}|=4-6$, $Q_{2}<0$, the majority of the particles in $a>250$ au orbits have $i\in(10,30)$ deg at $t=4$ Gyr, as observed. However, this range of $Q_{2}$ values is highly improbable (probability $<1\%$) in light of the Cassini data.

7 Summary and discussion

Solar system bodies in orbits with perihelia beyond Neptune and semimajor axes $150<a<1500$ au, dubbed extreme trans-Neptunian objects (ETNOs), have their orbits oriented in space in highly non-isotropic fashion. A general trend in numerical simulations is that any primordial clustering of ETNOs orbits is quickly washed away by the perturbations of giant planets. If the distributions of ETNOs angular orbital elements $\omega$, $\Omega$, and $i$ are not mainly products of an observational bias, then these statistics suggest that there should exist some external perturbation that acts continuously on ETNOs and forces the observed orbital clustering. A natural first guess is that the external perturbation comes from an unseen planet. It turns out that such a planet – Planet 9 (P9)– has to be massive $(M\sim 10\leavevmode\nobreak\ M_{\oplus})$; it has to orbit the Sun in a distant, eccentric, and moderately inclined orbit ($a\sim 400-700$ au, $q\sim 250$ au, $i\sim 30$ deg) (Batygin & Brown, 2016; Brown & Batygin, 2016). Moreover, it has to be far from perihelion right now (Fienga et al., 2016). In our contribution, we investigated the hypothesis that the external perturbation does not come from a massive planet in Newtonian dynamics, but from a subtle effect of more general theory of Milgromian dynamics (MD). The dominant effect of MD (formulated as modified gravity) in the ETNOs region is the external field effect (EFE); the external gravitational acceleration from the Galaxy282828Do not confuse this with the very different tidal effect, arising from differential gravity. Existence of EFE means that the strong equivalence principle is broken. does not decouple from the internal dynamics, it enters the equation of motion. We modelled EFE in the solar system as a quadrupole correction to the Newtonian potential along the direction of the external Galactic field; see Eq. (5). We accounted for a variation of the direction of the external field with time. The strength of EFE at a given position was parametrised with the parameter $Q_{2}$, and constrained by the Earth-Cassini range measurements to be $Q_{2}=3\pm 3\times 10^{-27}$ s*-2* (nominal $\pm 1\sigma$ value) (Hees et al., 2014). The action of the giant planets was incorporated into the model. Starting with isotropically oriented orbits of primordial ETNOs, we numerically investigated the evolution of these orbits during the subsequent 4 Gyr, considering various values of the parameter $Q_{2}$.

Our results can be summarised as follows:

• •

The orbit of Sedna can be explained by EFE torquing if $|Q_{2}|\gtrsim 1\times 10^{-27}$ s*-2*.

• •

The orbit of the other Sednoid, 2012 VP113, seems only marginally compatible with the EFE scenario of its origin; 2012 VP113 requires $Q_{2}>6\times 10^{-27}$ s*-2* , which is more than $1\sigma$ off of the nominal value of Hees et al. (2014), or $Q_{2}<4\times 10^{-27}$ s*-2* , which is more than $2.3\sigma$ off of the nominal value of Hees et al. (2014), to securely state that its perihelion was lifted purely by EFE.

• •

Similar to the P9 scenario, EFE does not induce strong confinement in $\omega-a$ plane until $a\gtrsim 600$ au, the clustering in $\omega$ is traced only by the potentially unstable orbits, and the inferred distribution of $\omega$ at this $a$-range is bimodal (peaking around $\omega=0$ and 180 deg); moreover, the confinement fades away when $Q_{2}$ is not close to $Q_{2}=1\times 10^{-27}$ s*-2*.

• •

Concerning perihelion longitude $L$, visually appealing clustering traced by the potentially unstable orbits develops for $a\gtrsim 300$ au and wide range of values of $Q_{2}$; if $Q_{2}<0$, the clustering centres close to 90 and 270 deg, while when $Q_{2}>0$, the clustering centres close to 0 (360) and 180 deg.

• •

In the stable population, the clustering in $L$ looms for $|Q_{2}|\gtrsim 4\times 10^{-27}$ s*-2*, with distribution of $L$ peaking close to 90 and 270 deg.

• •

The clustering in longitude of ascending node $\Omega$ is the most visually appealing and robust, spanning wide range of permitted $Q_{2}$ values, and traced by both stable and potentially unstable orbits with $a>150$ au; depending on the sign of $Q_{2}$, either $\Omega\in(0,180)$ deg $(Q_{2}<0)$ or $\Omega\in(180,360)$ deg $(Q_{2}>0)$ is preferred by the simulated ETNOs, especially when we look at $a>250$ au orbits.

• •

If $|Q_{2}|\gtrsim 4\times 10^{-27}$ s*-2*, $a>250$ au ETNOs with very low inclinations are suppressed, while the population with inclinations between 10 and 30 deg is relatively enhanced.

• •

High inclination and retrograde orbits of large semimajor axis Centaurs can be produced by EFE.

These results are in some tension with the present-day observations. First of all, the simulations did not reproduce the tightness of the clustering in $\omega$ down to $a=150$ au. Second, the clustering region in $\Omega$ is shifted by 180 deg with respect to the observed region, when one assumes, as usual, positive $Q_{2}$. There are not any interpolating function families known, giving arise to $Q_{2}<0$ in AQUAL or QUMOND at the present time, though maybe these are possible to construct. Assuming $Q_{2}$ to be negative solves the problem. For instance $Q_{2}=-1\times 10^{-27}$ s*-2* (1.3$\sigma$ off of the nominal $Q_{2}$ of Hees et al. (2014)) yields clustering in $\Omega$ in line with the observations. The observed clustering in $\Omega$ (if it is of dynamical origin) constrains $Q_{2}$ to $Q_{2}<0$. Perihelia of Sednoids and ETNOs distributions of $L$ and $i$ prefer $|Q_{2}|>4\times 10^{-27}$ s*-2*. The $Q_{2}$s lower than $Q_{2}=-4\times 10^{-27}$ s*-2* (2.3$\sigma$ off of the nominal $Q_{2}$ of Hees et al. (2014)) are very unlikely, with a probability of about $1\%,$ in light of the Cassini data.

Our results predict the effect MD could have on ETNOs and should be compared with future unbiased and robust data, when they are acquired. This means further constraints on the interpolating function as it is directly linked to the value of $Q_{2}$. Also, provided that P9 were found, $Q_{2}$ and the interpolating function would be further constrained in this case. Unfortunately such a test is not of the Popper’s standard of the hypotheses testing. We will not be able to rule out MD completely this way. It is always possible to construct interpolating function giving unmeasurable EFE. Moreover, there are flavours of MD completely without any EFE.

Does P9 do the job better? With P9, one has more freedom to fit the data. Planet 9 in an eccentric orbit can be made massive, hence its effect on ETNOs is strong, while it can be hidden close to its aphelion at the present time in order not to violate the Cassini data or planetary ephemerides292929Fienga et al. (2016) found that the true anomaly, $\phi$, most probably lies in the interval $(108,129)$ deg, using the same orbit of P9 as we considered in this paper.. There is much less freedom in the MD scenario. The external field direction rotates uniformly and very slowly (with period of one Galactic year) around the Sun, while $Q_{2}$, the only free parameter, stays approximately constant and bounded from above with aid of the Cassini data. The clustering in $\omega$ seems problematic in both frameworks; with the tuned orbit of Planet 9, $\omega$ is still not confined until $a\gtrsim 500$ au; see our Fig. 11 or Fig. 8 in Batygin & Brown (2016). Moreover, the results of Shankman et al. (2017) and Lawler et al. (2017) pose a problem for the P9 scenario, as already discussed above. It is unclear whether similar issues would also arise in the scenario where EFE, not P9, shapes the orbits of ETNOs.

In this paper we applied the Occam’s razor approach; we tried to explain the clustered orbits of ETNOs with the effect of MD only. Maybe a more complicated dynamical model, where a distant planetary perturber exists while EFE of MD is also relevant, could perform better in explaining ETNOs orbits, the overall structure of the Kuiper belt, and the origin of the distant planet at the same time.

Acknowledgements.

I am thankful to Jozef Világi and Leonard Kornoš for providing a computational facility. I also thank Vladímir Balek, Jozef Klačka, and Leonard Kornoš for discussions and Aurélien Hees for useful correspondence. This work was supported by Comenius University Grant no. UK/419/2016.