Issue 
A&A
Volume 603, July 2017



Article Number  A11  
Number of page(s)  19  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201630335  
Published online  30 June 2017 
Towards an explanation of orbits in the extreme transNeptunian region: The effect of Milgromian dynamics
Faculty of Mathematics, Physics and Informatics, Comenius University in Bratislava, Mlynská dolina, 842 48 Bratislava, Slovakia
email: pauco@fmph.uniba.sk
Received: 22 December 2016
Accepted: 16 March 2017
Context. 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 nonnegligible dynamical effects, which can be used to constrain the freedom in MD theories. Known extreme transNeptunian objects (ETNOs) have their argument of perihelion, longitude of ascending node, and inclination distributed in highly nonuniform fashion; ETNOs are bodies with perihelion distances greater than the orbit of Neptune and with semimajor axes greater than 150 au and less than ~1500 au. It is as if these bodies have been systematically perturbed by some external force.
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 longterm orbital evolution of primordial (randomised) ETNOs in MD.
Results. The EFE could produce distinct nonuniform 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.
Key words: minor planets, asteroids: general / Kuiper belt: general / gravitation
© ESO, 2017
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 nonexistence 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 transNeptunian 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 ETNOs^{1}. The outlying positions of 90377 Sedna and 2012 VP_{113} (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, ω, longitude of the ascending node, Ω, and perihelion longitude, L,^{2} which are distributed in a highly nonuniform 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 ± 0.06), and inclination, i (i = 20 ± 7 deg) of the ETNOs, where in brackets we present mean ± 1σ 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 inclinations^{3} (de la Fuente Marcos & de la Fuente Marcos 2014).
Batygin & Brown (2016) recognised that orbits of longperiodic ETNOs with a> 250 au (six such objects were known at that time^{4}: Sedna, 2004 VN_{112}, 2007 TG_{422}, 2010 GB_{174}, 2012 VP_{113}, and 2013 RF_{98}) are very tightly clustered in terms of ω, Ω, 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 ~ 10M_{⊕}, perihelion longitude 180 deg ahead of the six, orbiting in an eccentric, a ~ 700 au, e ~ 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 SR_{349} and 2013 FT_{28}). These two ETNOs have ω and i clustered in line with the original six, but one of them, 2013 FT_{28}, has Ω 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 ω, Ω, 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 presentday planetary architecture plus an additional 10 M_{⊕} 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 TullyFisher relation (McGaugh et al. 2000), the relation between stellar and dynamical central surface densities (Lelli et al. 2016), and the mass discrepancyacceleration 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 massdiscrepancy 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 ΛCDM model, where Λ is the positive cosmological constant representing dark energy, whose nature is unknown, and CDM stands for cold (nonbaryonic) 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 scaleinvariant 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 spacetime scale 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 freefall acceleration through the socalled 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 nonlinearity 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 ΛCDM”. The EFE can also explain the fact that rotation curves in the outer parts of galaxies are sometimes declining (Haghi et al. 2016). In the outer parts, when the internal accelerations are smaller than the external fields (which are ≲ a_{0}), the quasiNewtonian behaviour is restored (Famaey & McGaugh 2012).
Recently, McGaugh (2016) predicted velocity dispersion, σ_{vel}, of Crater II in MD. Although internal accelerations are very low (≪ a_{0}), EFE reduces σ_{vel} down to σ_{vel} ≈ 2 km s^{1}. For isolated Crater II (EFE is not present e.g. in modified inertia theories of MD) McGaugh obtained σ_{vel} ≈ 4 km s^{1}. Significantly higher σ_{vel} would mean falsification of MD. In ΛCDM, higher σ_{vel} is expected (BoylanKolchin et al. 2012; McGaugh 2016). The resolution was found recently, σ_{vel} = 2.7 ± 0.3 km s^{1} (Caldwell et al. 2017). 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 produce^{5} enough torque to deliver primordial solar system bodies to Sednalike orbits even in the presentday 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 Sects. 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 qfilter removes five objects (2013 UH_{15}, 2010 VZ_{98}, 2001 FP_{185}, 2013 FS_{28}, and 2015 SO_{20}) from the ETNOs sample. One of these objects, 2013 FS_{28}, is an outlier in terms of ω, with ω = 101.5 deg (Sheppard & Trujillo 2016). Histogram distributions of the orbital elements ω, Ω, 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 ≈ 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 welldetermined 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 qpopulations are distinct in terms of ω and Ω when a> 150 au. In the case of ω, there are no ETNOs with q> 35 au, but many objects with a> 150 au, q< 35 au, occupying the region around ω = 180 deg. Object 2003 SS_{422} is listed in the NASA Jet Propulsion Laboratory (JPL) small body database with orbital elements (and their 1σ uncertainties) at Epoch 2457600.5 (July 31, 2016): a = 193.5 ± 47.8 au, q = 39.4 ± 1.0 au, i = 16.8 ± 0.1 deg, ω = 209.8 ± 17.2 deg, and Ω = 151.1 ± 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 ω 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 Ω 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.
Fig. 1 Histogram distributions of argument of perihelion (top left panel), longitude of ascending node (top right panel), perihelion longitude (bottom left panel), and inclination (bottom right panel) inferred from the welldetermined orbits of TNOs fulfilling the criteria 100 <a< 1500 au, q> 35 au (colour histograms). The orbits were discerned according to their a and divided into 3 groups: 100 <a< 150 au (yellow histogram), 150 <a< 250 au (blue histogram), and 250 <a< 1500 au (green histogram). The distributions inferred from all welldetermined orbits of ETNOs (150 <a< 1500 au, q> 30 au) are indicated as dashed bars. The blue and green bars in a given bin can, but do not have to, sum up to the height of a dashed bar in that bin, as the dashed bar also comprises objects with q in the range from 30 to 35 au. The histogram bars are aligned to the left, for example, in the bottom right panel, 0 indicates inclination range (0,10) deg, 10 indicates inclination range (10,20) deg, and so on. 
Fig. 2 Histogram distributions of argument of perihelion (top left panel), longitude of ascending node (top right panel), perihelion longitude (bottom left panel), and inclination (bottom right panel) inferred from the welldetermined orbits of TNOS and Centaurs fulfilling the criteria 100 <a< 1500 au, q< 35 au (colour histograms). The distributions inferred from all welldetermined orbits of ETNOs (150 <a< 1500 au, q> 30 au) are indicated as dashed bars. 
Looking at TNOs in Fig. 1, the nonuniformity of the distributions is apparent, with a nonuniformity degree being adependent, 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 longperiodic 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 longperiodic ETNOs. Clearly, we are dealing with smallnumber 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 ω, Ω, 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 ω, Ω, 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, θ_{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, , , , . The probability is given as a frequency of sets for which is smaller than the circular variance of the observational data.
Brown & Firth (2016) plotted the distribution of ω 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 ω 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 ω 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 ω is twofold: (1) if one does observations close to the ecliptic (as many surveys do), ωs close to ω = 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 ω ∈ (180,360) deg, and by contrast, ω ∈ (0,180) is preferred when observations are made from the northern hemisphere. By employing offecliptic observations, made primarily from Chile, and taking into account the fact that there are no q> 35 au ETNOs with ω close to 180 deg, Sheppard & Trujillo (2016) concluded that the clustering in ω for q> 35 au ETNOs is real with high probability.
Fig. 3 Unbiased cumulative distribution function of the scattered KBOs (dashed line) vs. observed cumulative distribution of ETNOs (solid line), both as a function of inclination. The thick dashed line represents the bestfit model of Gulbis et al. (2010), while the thinner dashed lines are its 1σ variations. The dots correspond to the observed inclinations of ETNOs and are equidistantly spaced between 0 and 1 on the vertical axis. 
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 survey^{6}, 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 ω clustered around ω ≈ 340 deg and L biased, Ω is biased as well. The Galactic plane avoidance then suppresses ETNOs with Ω around 110, where oddly^{7} the majority of ETNOs reside, and around 290 deg. The asymmetry of the Ω distribution cannot be explained by the intrinsic declination bias, as this is symmetric with respect to the preferred value Ω = 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 KBOs^{8}. 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 sini (Brown 2001), centred at deg, and with a width of 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 ω, Ω, 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 g = g^{N} but from(1)where g^{N} is expected Newtonian gravitational acceleration, g^{N} is its norm, ν(y), y ≡ g^{N}/a_{0} is a modification factor, the socalled interpolating function, fulfilling behaviour ν(y) → 1 for y ≫ 1 (Newtonian limit) and ν(y) → y^{− 1/2} for y ≪ 1 (deepMD limit), and a_{0} ~ 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) galaxies^{9}, 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) (2)where the latter equality is the classical Poisson equation with Newtonian potential Φ_{N}, (baryonic) matter density ρ, and Newton’s gravitational constant G, μ is inverse of ν from Eq. (1), and the equation of motion is g = −∇Φ. The second modified Poisson approach, known as quasilinear MOND (QUMOND), states that the governing gravitational potential Φ is given by (Milgrom 2010) (3)Having μ as inverse of ν, 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.
Equation (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 fullfledged modified inertia theory developed yet (but see discussion in Sect. 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 ≪ r_{M}, where r_{M} ≡ (GM_{⊙}/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, δg = g−g^{N}, can be written as (4)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 MD^{10} and it comes from the nonlinearity of the modified Poisson equation. An embedding external field that is constant and uniform, with magnitude g_{e}, g_{e} ≲ 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 (μ and ν arbitrarily close to 1). The solar system is falling freely in the external gravitational field of the Galaxy. The external field has magnitude km^{2} s^{2}/8.0 kpc ≐ 2.0 × 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 ≪ r_{M}, r_{M} ≡ (GM_{⊙}/a_{0})^{1/2} ~ 7000 au, where r_{M} is characteristic radius in MD solar system, EFE expresses itself dominantly as an anomalous quadrupole, δΦ, so a governing potential entering the equation of motion is written (5)where e ≡ g_{e}/g_{e} is a unitary vector pointing towards the Galactic centre^{11} (GC) and Φ_{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 yield^{12}Q_{2}> 0 (Milgrom 2009a; Blanchet & Novak 2011). If Q_{2} is taken to be a constant in some region of interest, the value δΦ 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 theoretical maximum allowed by the EarthCassini spacecraft range data (hereafter Cassini data), Q_{2} ~ 10^{27} s^{2} (Hees et al. 2014), then EFE is ~ 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 GCanticentre, characterised by the tidal parameter GM′/r^{′3} = Q_{2}/ 3, where M′ is mass of P9 and r′, r′ ≫ 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 presentday 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 nonspherically symmetric matter distribution. Suppose we have an isolated highacceleration system S, which assumes μ = ν = 1 to the desired accuracy everywhere within S, of total mass M and extension R, R ≪ r_{M}, where . Suppose also that mass distribution within S, ρ(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 (6)where Φ_{N} is the expected Newtonian potential, δΦ is a dominant MD correction from the asphericity effect, a quadrupole field, and α, typically α ~ 1, is a numerical factor that depends on the form of the interpolating function. In the case of the solar system^{13}, M ≈ M_{⊙} 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} ~ a_{0}, but beyond g^{N} ~ κa_{0}, α in Eq. (6) becomes enhanced by κ^{2} approximately (Milgrom 2012). An example of such a theory is, for instance TeVeS (Bekenstein 2004), where κ 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 , seems rather unnatural. The anomalous perihelion precession of Saturn constrains κ to be κ ≤ 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 ~ 900 au.
By combining an analysis of rotation curves of galaxies and solar system constraints^{14} 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^{15}. Referring to Table 2 in Hees et al. (2016), we chose interpolating function family as most representative. The Cassini data restrict to n ≥ 2. The enhanced gravity effect is completely negligible in the whole solar system^{16} for , n ≳ 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 ± 3 × 10^{27} s^{2} (nominal ± 1σ 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 e of Eq. (5), pointing towards the GC, is assumed to circulate uniformly with time according to (7)where the first three lines are components of e in the Galactic coordinates^{17} 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 GC^{18}, lying in the Galactic midplane, with angular speed w; τ 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 longliving 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 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 ω−e, Ω−e, and L−e plane, considering Q_{2} taken from a set Q_{2} = { 1.0, 2.0, 4.0, 6.0 }. In MD, the external perturbation, substituting a distant planet, is the quadrupole anomaly arising from EFE; see Eqs. (5) and (7). The direction of the perturbation is known and we alter only the parameter Q_{2}.
Fig. 4 Evolution of test particles in L−e plane as inferred from the P9 simulation. The particles and P9 orbit the Sun in the same plane by design. The eccentricity range of the observed ETNOs is indicated by the dashed horizontal lines. Two transparency levels are used as a proxy for two different time windows of the simulation. The less transparent evolutionary tracks represent longliving 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. Three starting semimajor axes were used as follows: a = 150 au (first column), a = 350 au (second column), and a = 550 au (third column). Nbody encounters are completely omitted in this case, hence semimajor axes do not change during the simulation. The secular effect of the giant planets is modelled with a strong J_{2} moment of the Sun. There were no stable orbits at a = 550 au level. 
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 ^{19}, i drawn randomly from a uniform distribution in range (0, 5) deg, ω drawn randomly from a uniform distribution in range (0, 360) deg, and Ω 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 presentday GCanticentre. The ecliptic latitude of the GC is only minus few degrees. Here our intention is to demonstrate the existence of the longliving 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 Nbody integration software package (Chambers 1999). The hybrid symplecticBulirschStoer algorithm with a timestep equal to a tenth of the orbital period of Neptune was employed. Phase τ in Eq. (7) was chosen in a way that at t = 4 Gyr the vector e points towards the presentday GC. We tested that the value of τ 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 simulation^{20}. 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′ = 10 M_{⊕} was assumed to orbit in a′ = 700 au, e′ = 0.6, i′ = 0 deg, and L′ = 240 deg orbit (Batygin & Brown 2016). The primed quantities always refer to P9. The test particle orbits were initialised with L−L′ equal to 0 (20 particles for each a) and 180 deg (20 particles for each a), and nil inclinations^{21}. 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′ = 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.
Fig. 5 Evolution of test particles in ω−e (first row), Ω−e (second row), and L−e (third row) plane as inferred from the benchmark model simulation with Q_{2} = 1. 
Fig. 6 Evolution of test particles in ω−e (first row), Ω−e (second row), and L−e (third row) plane as inferred from the benchmark model simulation with negative Q_{2}, Q_{2} = −1. 
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′ = 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 longliving 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 ≈ 240 deg (aligned with P9) for nearcircular orbits. For a = 350 au, two longterm 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 ≈ 60 deg and antialigned with P9 and, second, the secular dynamics regime concerning less eccentric orbits, centring at L ≈ 240 deg and aligned with P9. These absolute positions come from our choice of L′, which was tuned to match the observed clustering position; see Fig. 1. Batygin & Brown (2016) always refer to the relative longitude of perihelion ϖ−ϖ′ 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′ 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′ = 30 deg and inclined test particle orbits, there were no longterm stable orbits.
The picture in MD is rich as well. Orbits with a = 150 au are characterised by trivial circulation of the angles (L, Ω,ω). No stable libration develops at this alevel 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 Ω in all 4 MD simulations. Both eccentric (ETNOs) and less eccentric particles, with both a = 350 and 550 au, strongly prefer Ω ∈ (180,360) deg. In ETNOs eccentricity range (EER), e ∈ (0.69,0.93), stable, and often tight, libration around Ω ≈ 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} ≲ 4, and neither circulation nor strong confinement using models with Q_{2} ≳ 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 ≈ 90 and L ≈ 270 deg in a broad range of eccentricities, spanning from low eccentricities to EER. For eccentricities e ≳ 0.9, the circulation of the angles sets in. In ω−e plane, there are orbits clustered near ω ≈ 0 and ω ≈ 180 deg, below and at the low eccentricity end of EER. We depict the evolution of test particles in ω−e, Ω−e, and L−e plane as inferred from the benchmark model simulations using Q_{2} = 1 in Fig. 5.
Despite the rich structure in ω−e, Ω−e, and L−e planes in the MD, the simulated evolutionary paths are in sharp contrast with presentday observations. The simulations predict that the vast majority of ETNOs should have Ω ∈ (180,360) deg. Hence, the vast majority of ETNOs should have ascending node lying in the same halfspace, with the GC lying in the same halfspace. This is not an observed trend. Of 19 ETNOs, 14 have Ω ∈ (0,180) deg, while 3 of them have Ω very close to 180 deg (Sheppard & Trujillo 2016). Moreover, the tightness of the observed clustering in ω 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.
The fact that the MD models predict the clustering region that lies roughly 180 deg ahead of the observed clustering in Ω 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 that^{22} although Q_{2} was positive for all interpolating functions μ he studied, he was “not able to prove that this is always the case from basic properties of μ” (Milgrom 2009a). Hees et al. (2014) showed that the Cassini data yield Q_{2} = 3 ± 3 (nominal ± 1σ value).
We substituted Q_{2} → −Q_{2} and repeated the benchmark model simulations. This means we considered Q_{2}s that are 1−3σ 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 Ω, while leaving evolution in ω−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, Ω 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 Nbody 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 , 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 (ω,Ω). 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, 0.5, 1.0, 2.0, 4.0, 6.0 }. For the sake of comparison, simulations with P9 in Newtonian dynamics (no EFE ⇔Q_{2} = 0) were carried out. Two initial P9 orbits were considered: (#1) a′ = 700 au, e′ = 0.6, i′ = 30 deg, ω′ = 150 deg, Ω′ = 113 deg (Batygin & Brown 2016; Brown & Batygin 2016; Fienga et al. 2016); (#2) the same^{23}a′, e′, i′ as for orbit (#1), but ω′ = 100 deg, Ω′ = 140 deg. The orbit (#2) orbital elements account for variation in ω′ and Ω′, caused by the giant planets, and at t = 4 Gyr the final ω′ and Ω′ were very close to the orbit (#1) starting values.
6.1. Sednoids
The discovery of Sedna (Brown et al. 2004) and 2012 VP_{113} (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 presentday 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 planetmass 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/bottommost 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 VP_{113}, 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 VP_{113} still lies at the lowend border of the achievable eccentricities at a ≈ 250 au. Objects with higher eccentricities are easier to detect. The orbit of 2012 VP_{113} 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 lowestorder correction (octupole) neglected in Eq. (5) and in the benchmark model? This correction (adds to RHS of Eq. (5)) can be written as^{24} (Blanchet & Novak 2011) (8)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 μ, for instance the μ_{n}(x) = x/ (1 + x^{n})^{1 /n} family, in the framework of AQUAL, and found that  Q_{3}  < 1.2 × 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 μ_{1}(x) = x/ (1 + x). It is known that μ_{1}, and its inverse ν_{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 μ_{n} family, allowed by the solar system data (Blanchet & Novak 2011; Hees et al. 2016), yield typically  Q_{3}  ~ 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 μ and ν, and a given external field strength. Approximation of this kind is common; see e.g. Milgrom (2009a).
Fig. 7 a−e diagram from the last Gyr of the benchmark model simulation incorporating gravitational interaction with Neptune. Q_{2} = 1 (top left panel), Q_{2} = 2 (top right panel), Q_{2} = 4 (bottom left panel), and Q_{2} = 6 (bottom right panel) were assumed. Two colours (transparency levels) are used as a proxy for two inclination ranges. The blue (less transparent) evolutionary tracks represent orbits with inclinations spanning i ∈ (10,30) deg, the interval the inclinations of Sednoids lie in. The grey (more transparent) tracks refer to all the other longliving orbits. Sedna and 2012 VP 113 are indicated as large dots. The dashed lines indicate q = 30, 50, 100, and 550 au levels (from right/bottom to left/upwards). Objects between the two leftmost/uppermost lines (q in range between 100 and 550 au) would hardly be observed if they really exist. 
Fig. 8 a−e diagram from the last Gyr of the benchmark model simulation incorporating gravitational interaction with Neptune. The enhanced gravity effect is accounted for in the model using the interpolating function . yields Q_{2} = 5.9. Only the evolution of test particles in orbits with inclinations spanning i ∈ (10,30) deg, which is the interval that the inclinations of Sednoids lie in, is shown. Sedna and 2012 VP 113 are indicated as large dots. The dashed lines indicate q = 30, 50, 100, and 550 au levels (from right/bottom to left/upwards). 
Fig. 9 a−e diagram from the last Gyr of the P9 simulation incorporating gravitational interaction with Neptune. Only evolution of test particles in orbits with inclinations spanning i ∈ (10,30) deg, which is the interval that the inclinations of Sednoids lie in, is shown. Sedna and 2012 VP 113 are indicated as large dots. The dashed lines indicate q = 30, 50, 100, and 550 au levels (from right/bottom to left/upwards). 
Switching back to family and QUMOND, yields Q_{2} = 5.9 when one assumes η ≡ g_{e}/a_{0} = 2.3, a_{0} = 0.815 × 10^{10} m s^{2} (Hees et al. 2016). The quoted value of a_{0} comes from rotation curves fits using . With one gets, besides EFE, as well as a nonnegligible enhanced gravity effect. We incorporated the enhanced gravity effect, Eq. (4), into the model and repeated the simulation. We used , , where is Newtonian gravitation from the Sun, and was found from knowing g_{e}. Direction of fields g_{e} and rotate as prescribed by Eq. (7). An outcome is in Fig. 8. 2012 VP_{113} still lies at the lowend border of the achievable eccentricities at a ≈ 250 au. We remind that , n< 2, is disfavoured by the Cassini data (Hees et al. 2016), and for , n ≳ 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.
Fig. 10 a−e diagram from the last Gyr of the benchmark model simulation incorporating gravitational interaction with Neptune. Negative values of Q_{2}, Q_{2} = −1 (left panel) and Q_{2} = −4 (right panel) were used. Two colours (trasparency levels) are used as a proxy for two inclination ranges. The blue (less transparent) evolutionary tracks represent orbits with inclinations spanning i ∈ (10,30) deg, which is the interval that the inclinations of Sednoids lie in. The grey (more transparent) tracks refer to all the other longliving orbits. Sedna and 2012 VP 113 are indicated as large dots. The dashed lines indicate q = 30, 50, 100, and 550 au levels (from right/bottom to left/upwards). 
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 τ and angular speed of the Sun w, present in Eq. (7). 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 VP_{113} lies at the lowend border of the achievable eccentricities at a ≈ 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 KolmogorovSmirnov (KS) test^{25} comparing distributions of ω, Ω, 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.
Results of the KS test comparing distributions of ω, Ω, L, and i of the observational data and those inferred from the benchmark model simulation at t = 4 Gyr.
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), pvalues, p( > 250  > 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 (~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 Nbody 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 pvalues, p(150−250  150−250), similar to the control test, except for the distribution of Ω when even lower pvalues are reported. This means some clustering in Ω 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} ≳ 4. However, the compatibility between the observational data and the MD model remains low in the case of ω and L, though pvalues, p( > 250  > 250), are slightly higher than in the control test if, e.g. Q_{2} = 4. The distribution of Ω 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.
Results of the KS test taking into account only the orbits, of ETNOs and test particles, with q> 35 au.
Results of the KS test after doubling the initial number of test particles.
We now investigate the simulated data visually. Figures 11, 12 show evolution of test particle orbits in ω−a, Ω−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 fullwidth horizontal lines. Clustering in ω emerges only for a ≳ 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 ω 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 ω ≈ 0 to the number of particles with ω ≈ 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.
Fig. 11 Evolution of test particle orbits, taken from the P9 simulation, in ω−a plane. Only particles surviving at least 3 Gyr are taken into account. The evolutionary tracks of the particles are shown in the last Gyr of the simulation every time their orbits met the criteria: i< 50 deg, q< 80 au. The tracks below a = 200 au are not shown as these are just a forest of fullwidth horizontal lines. The X symbols indicate the positions of potentially unstable particles, which are discarded at some time between t = 3 and 4 Gyr, at the time of their destruction. The plus signs indicate positions of longliving particles at t = 4 Gyr. The known ETNOs with q> 35 au (circles) and q< 35 au (triangles) are plotted. TNO 2003 SS_{422} is denoted as a square and its large uncertainty in a is indicated by the error bar. Shaded regions indicate the observed clustering regions for a> 250 au ETNOs. 
Fig. 12 Evolution of test particle orbits, taken from the P9 simulation, in Ω−a (left panel) and L−a plane (right panel). 
We turn our attention to MD. Figure 13 shows the evolution of test particle orbits in ω−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 ω 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 ω near 0 (360) as well as near 180 deg is preferred by the test particles. With Q_{2} ≳ 2, the clustering in ω fades away.
The test particle orbits strongly cluster around Ω ≈ 270 deg. This is illustrated in the left panel of Fig. 14, where we show results from Q_{2} = 2 simulation. Strong confinement in Ω develops for stable orbits in a wide range of Q_{2} values, Q_{2} ∈ (0.5,6), and also for potentially unstable orbits if Q_{2} ≳ 2.
Fig. 13 Evolution of test particle orbits, taken from the benchmark model simulation, in ω−a plane. The strength of EFE was assumed to be Q_{2} = 1. Only particles surviving at least 3 Gyr are taken into account. The evolutionary tracks of the particles are shown in the last Gyr of the simulation every time their orbits met the criteria, i.e. i< 50 deg, q< 80 au. The tracks below a = 500 au are not shown as these are just fullwidth horizontal lines. 
The distribution of L shows a coherent picture for Q_{2} in the range Q_{2} ∈ (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 ω, Ω, 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 ω. Strong clustering emerges in terms of Ω across wide range of Q_{2} values, Q_{2} ∈ (0.5,6), as majority of the orbits prefer Ω ∈ (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} ≈ 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 lowinclination 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 ∈ (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.
Fig. 14 Evolution of test particle orbits, taken from the benchmark model simulation, in Ω−a (left panel) and L−a (right panel) plane. The strength of EFE was assumed to be Q_{2} = 2. Only particles surviving at least 3 Gyr are taken into account. Evolutionary tracks of the particles are shown in the last Gyr of the simulation every time their orbits met the criteria, i.e. i< 50 deg, q< 80 au. The tracks below a = 400 au are not shown as these are just fullwidth horizontal lines. Increasing Q_{2}, the clustering in L sets in at lower semimajor axes, e.g. with Q_{2} = 6; similar clustering as in the figure develops for a> 300 au. 
Fig. 15 Histogram distributions of ω (top left panel), Ω (top right panel), L (bottom left panel), and i (bottom right panel) of test particle orbits at the end of the benchmark model simulation, i.e. at t = 4 Gyr. The strength of EFE was assumed to be Q_{2} = 1. Orbits are discerned according to their semimajor axis and divided into three groups: 150 <a< 250 au (blue histograms), 250 <a< 550 au (green histograms), and 550 <a< 1500 au (red histograms). The observed distributions for ETNOs with 250 <a< 1500 au, multiplied by factor 4, are shown as well (dashed histograms). Histograms are aligned to the left. 
Fig. 16 Histogram distributions of ω (top left panel), Ω (top right panel), L (bottom left panel), and i (bottom right panel) of test particle orbits at the end of the benchmark model simulation, i.e. at t = 4 Gyr. The strength of EFE was assumed to be Q_{2} = 4. 
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 region^{26} 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 Sect. 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} ~ 10^{2}; the inclusion of the enhanced gravity effect in the case^{27}Q_{2} = 5.9; the variation of w and τ in Eq. (7); 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 ω, Ω, L, and i, of simulated and observed ETNOs. The pvalues 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σ range given by the Casssini data). In the case of ω, L, and i distributions, the reported pvalues are similar to their Q_{2}> 0 counterparts as long as models with the same  Q_{2}  are compared. In the case of Ω, the pvalues are much higher.
Evolution in Ω−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} ∈ (−6,−1) yields a qualitatively similar picture as in Fig. 18. Potentially unstable orbits cluster around L ≈ 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 ω in our simulations with Q_{2}< 0.
Histogram distributions of ω, Ω, L, and i inferred at t = 4 Gyr from the Q_{2} = −1 simulation are shown in Fig. 19. The clustering in Ω 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 ∈ (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.
Fig. 17 Perihelion distance as a function of time for 4 selected orbits from the benchmark model simulation characterised by highq amplitude. Assumed values of Q_{2} and initial (final) values of a and q are indicated. 
Results of the KS test comparing distributions of ω, Ω, L, and i of the observational data and those inferred from the benchmark model simulation at t = 4 Gyr.
7. Summary and discussion
Solar system bodies in orbits with perihelia beyond Neptune and semimajor axes 150 <a< 1500 au, dubbed extreme transNeptunian objects (ETNOs), have their orbits oriented in space in highly nonisotropic 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 ω, Ω, 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 ~ 10 M_{⊕}); it has to orbit the Sun in a distant, eccentric, and moderately inclined orbit (a ~ 400−700 au, q ~ 250 au, i ~ 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 Galaxy^{28} 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 EarthCassini range measurements to be Q_{2} = 3 ± 3 × 10^{27} s^{2} (nominal ± 1σ 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}.
Fig. 18 Evolution of test particle orbits, taken from the benchmark model simulation in Ω−a (left panel) and L−a (right panel) plane. A negative Q_{2}, Q_{2} = −2 was assumed. Only particles surviving at least 3 Gyr are taken into account. Evolutionary tracks of the particles are shown in the last Gyr of the simulation every time their orbits met the criteria, i.e. i < 50 deg, q< 80 au. The tracks below a = 400 au are not shown as these are just fullwidth horizontal lines. Decreasing Q_{2}, Q_{2} < 0, the clustering in L sets in at lower semimajor axes, e.g. with Q_{2} = −6 similar clustering as in the figure develops for a > 250 au. 
Fig. 19 Histogram distributions of ω (top left panel), Ω (top right panel), L (bottom left panel), and i (bottom right panel) of test particle orbits at the end of the benchmark model simulation, i.e. at t = 4 Gyr. A negative Q_{2}, Q_{2} = −1 was assumed. Orbits are discerned according to their semimajor axis and divided into three groups: 150 < a < 250 au (blue histograms), 250 < a < 550 au (green histograms), and 550 < a < 1500 au (red histograms). The observed distributions for ETNOs with 250 < a < 1500 au, multiplied by factor 4, are shown as well (dashed histograms). Histograms are aligned to the left. 
Our results can be summarised as follows:

The orbit of Sedna can be explained by EFE torquing if Q_{2}  ≳ 1 × 10^{27} s^{2}.

The orbit of the other Sednoid, 2012 VP_{113}, seems only marginally compatible with the EFE scenario of its origin; 2012 VP_{113} requires Q_{2}> 6 × 10^{27} s^{2}, which is more than 1σ off of the nominal value of Hees et al. (2014), or Q_{2}< 4 × 10^{27} s^{2}, which is more than 2.3σ 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 ω−a plane until a ≳ 600 au, the clustering in ω is traced only by the potentially unstable orbits, and the inferred distribution of ω at this arange is bimodal (peaking around ω = 0 and 180 deg); moreover, the confinement fades away when Q_{2} is not close to Q_{2} = 1 × 10^{27} s^{2}.

Concerning perihelion longitude L, visually appealing clustering traced by the potentially unstable orbits develops for a ≳ 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}  ≳ 4 × 10^{27} s^{2}, with distribution of L peaking close to 90 and 270 deg.

The clustering in longitude of ascending node Ω 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 Ω ∈ (0,180) deg (Q_{2}< 0) or Ω ∈ (180,360) deg (Q_{2}> 0) is preferred by the simulated ETNOs, especially when we look at a> 250 au orbits.

If  Q_{2}  ≳ 4 × 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 presentday observations. First of all, the simulations did not reproduce the tightness of the clustering in ω down to a = 150 au. Second, the clustering region in Ω 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 × 10^{27} s^{2} (1.3σ off of the nominal Q_{2} of Hees et al. 2014) yields clustering in Ω in line with the observations. The observed clustering in Ω (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 × 10^{27} s^{2}. The Q_{2}s lower than Q_{2} = −4 × 10^{27} s^{2} (2.3σ 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 ephemerides^{29}. 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 ω seems problematic in both frameworks; with the tuned orbit of Planet 9, ω is still not confined until a ≳ 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.
Longitude of an object at perihelion in the heliocentric ecliptic coordinates. It is not equal to the longitude of perihelion, ϖ ≡ ω + Ω, 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 ω, Ω, and i give L and B unambiguously: sin(L−Ω) = sinωcosi/ cosB, cos(L−Ω) = cosω/ cosB.
This 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).
In Milgrom (2009a), a dimensionless parameter q, , was employed instead of Q_{2}.
The 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).
Blanchet & Novak (2011) state an opposite sign in the RHS of Eq. (8), as they have been using ∇^{2}Φ_{N} = −4πGρ and g = ∇Φ. We used ∇^{2}Φ_{N} = 4πGρ and g = −∇Φ throughout the paper, hence the difference.
To be clear in statistical terms and notation, we give a following example: the probability pvalue 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.
Fienga et al. (2016) found that the true anomaly, φ, most probably lies in the interval (108,129) deg, using the same orbit of P9 as we considered in this paper.
Acknowledgments
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.
References
 Batygin, K., & Brown, M. E. 2016, AJ, 151, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Bekenstein, J. D. 2004, Phys. Rev. D, 70, 083509 [NASA ADS] [CrossRef] [Google Scholar]
 Bekenstein, J., & Milgrom, M. 1984, ApJ, 286, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Blanchet, L., & Novak, J. 2011, MNRAS, 412, 2530 [NASA ADS] [CrossRef] [Google Scholar]
 Bosma, A. 1981, AJ, 86, 1825 [NASA ADS] [CrossRef] [Google Scholar]
 BoylanKolchin, M., Bullock, J. S., & Kaplinghat, M. 2012, MNRAS, 422, 1203 [NASA ADS] [CrossRef] [Google Scholar]
 Brasser, R., Duncan, M. J., Levison, H. F., Schwamb, M. E., & Brown, M. E. 2012, Icarus, 217, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, M. E. 2001, AJ, 121, 2804 [Google Scholar]
 Brown, M. E., & Batygin, K. 2016, ApJ, 824, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, R. B., & Firth, J. A. 2016, MNRAS, 456, 1587 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, M. E., Trujillo, C., & Rabinowitz, D. 2004, ApJ, 617, 645 [NASA ADS] [CrossRef] [Google Scholar]
 Caldwell, N., Walker, M. G., Mateo, M., et al. 2017, ApJ, 839, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Chambers, J. E. 1999, MNRAS, 304, 793 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Cowan, N. B., Holder, G., & Kaib, N. A. 2016, ApJ, 822, L2 [NASA ADS] [CrossRef] [Google Scholar]
 de la Fuente Marcos, C., & de la Fuente Marcos, R. 2014, MNRAS, 443, L59 [NASA ADS] [CrossRef] [Google Scholar]
 de la Fuente Marcos, C., & de la Fuente Marcos, R. 2016, MNRAS, 459, L66 [NASA ADS] [CrossRef] [Google Scholar]
 de la Fuente Marcos, C., & de la Fuente Marcos, R., & Aarseth, S. J. 2016, MNRAS, 460, L123 [NASA ADS] [CrossRef] [Google Scholar]
 Desmond, H. 2017, MNRAS, 464, 4160 [CrossRef] [Google Scholar]
 Dones, L., Brasser, R., Kaib, N., & Rickman, H. 2015, Space Sci. Rev., 197, 191 [NASA ADS] [CrossRef] [Google Scholar]
 Einstein, A. 1916, Annalen der Physik, 354, 769 [NASA ADS] [CrossRef] [Google Scholar]
 Famaey, B., & McGaugh, S. S. 2012, Liv. Rev. Rel., 15, 10 [Google Scholar]
 Fienga, A., Laskar, J., Manche, H., & Gastineau, M. 2016, A&A, 587, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fortney, J. J., Marley, M. S., Laughlin, G., et al. 2016, ApJ, 824, L25 [NASA ADS] [CrossRef] [Google Scholar]
 Gomes, R. S., Matese, J. J., & Lissauer, J. J. 2006, Icarus, 184, 589 [NASA ADS] [CrossRef] [Google Scholar]
 Gomes, R. S., Soares, J. S., & Brasser, R. 2015, Icarus, 258, 37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gulbis, A. A. S., Elliot, J. L., Adams, E. R., et al. 2010, AJ, 140, 350 [NASA ADS] [CrossRef] [Google Scholar]
 Haghi, H., Bazkiaei, A. E., Zonoozi, A. H., & Kroupa, P. 2016, MNRAS, 458, 4172 [NASA ADS] [CrossRef] [Google Scholar]
 Hees, A., Folkner, W. M., Jacobson, R. A., & Park, R. S. 2014, Phys. Rev. D, 89, 102002 [NASA ADS] [CrossRef] [Google Scholar]
 Hees, A., Famaey, B., Angus, G. W., & Gentile, G. 2016, MNRAS, 455, 449 [NASA ADS] [CrossRef] [Google Scholar]
 Holman, M. J., & Payne, M. J. 2016a, AJ, 152, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Holman, M. J., & Payne, M. J. 2016b, AJ, 152, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Iorio, L. 2010, The Open Astron. J., 3, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Iorio, L. 2013, Classical and Quantum Gravity, 30, 165018 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Jílková, L., Portegies Zwart, S., Pijloo, T., & Hammer, M. 2015, MNRAS, 453, 3157 [Google Scholar]
 Kaib, N. A., & Quinn, T. 2008, Icarus, 197, 221 [NASA ADS] [CrossRef] [Google Scholar]
 Kaib, N. A., Roškar, R., & Quinn, T. 2011, Icarus, 215, 491 [NASA ADS] [CrossRef] [Google Scholar]
 Lawler, S. M., Shankman, C., Kaib, N., et al. 2017, AJ, 153, 33 [NASA ADS] [CrossRef] [Google Scholar]
 Lelli, F., McGaugh, S. S., Schombert, J. M., & Pawlowski, M. S. 2016, ApJ, 827, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Lelli, F., McGaugh, S. S., Schombert, J. M., & Pawlowski, M. S. 2017, ApJ, 836, 152 [Google Scholar]
 Le Verrier, U. J. 1846, Astron. Nachr., 25, 65 [NASA ADS] [Google Scholar]
 Le Verrier, U. J. 1859, Annales de l’Observatoire de Paris, 5 [Google Scholar]
 McGaugh, S. S. 2004, ApJ, 609, 652 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S. S. 2016, ApJ, 832, L8 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S., & Milgrom, M. 2013a, ApJ, 766, 22 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S., & Milgrom, M. 2013b, ApJ, 775, 139 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S. S., Schombert, J. M., Bothun, G. D., & de Blok, W. J. G. 2000, ApJ, 533, L99 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 McGaugh, S. S., Lelli, F., & Schombert, J. M. 2016, Phys. Rev. Lett., 117, 201101 [NASA ADS] [CrossRef] [Google Scholar]
 McMillan, P. J., & Binney, J. J. 2010, MNRAS, 402, 934 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 1983a, ApJ, 270, 371 [Google Scholar]
 Milgrom, M. 1983b, ApJ, 270, 365 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 1994, Annals of Physics, 229, 384 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 2009a, MNRAS, 399, 474 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 2009b, ApJ, 698, 1630 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 2010, MNRAS, 403, 886 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 2012, MNRAS, 426, 673 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 2014, MNRAS, 437, 2531 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 2016a, ArXiv eprints [arXiv:1605.07458] [Google Scholar]
 Milgrom, M. 2016b, Phys. Rev. Lett., 117, 141101 [NASA ADS] [CrossRef] [Google Scholar]
 Morbidelli, A., & Levison, H. F. 2004, AJ, 128, 2564 [NASA ADS] [CrossRef] [Google Scholar]
 Paučo, R., & Klačka, J. 2016, A&A, 589, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rubin, V. C., Ford, Jr., W. K., Thonnard, N., & Burstein, D. 1982, ApJ, 261, 439 [NASA ADS] [CrossRef] [Google Scholar]
 Sancisi, R. 2004, in Dark Matter in Galaxies, eds. S. Ryder, D. Pisano, M. Walker, & K. Freeman, IAU Symp., 220, 233 [Google Scholar]
 Sanders, R. H. 1990, A&ARv, 2, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Schönrich, R. 2012, MNRAS, 427, 274 [NASA ADS] [CrossRef] [Google Scholar]
 Shankman, C., Kavelaars, J. J., Lawler, S. M., Gladman, B. J., & Bannister, M. T. 2017, AJ, 153, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Sheppard, S. S., & Trujillo, C. 2016, AJ, 152, 221 [NASA ADS] [CrossRef] [Google Scholar]
 Trujillo, C. A., & Sheppard, S. S. 2014, Nature, 507, 471 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Results of the KS test comparing distributions of ω, Ω, L, and i of the observational data and those inferred from the benchmark model simulation at t = 4 Gyr.
Results of the KS test taking into account only the orbits, of ETNOs and test particles, with q> 35 au.
Results of the KS test comparing distributions of ω, Ω, L, and i of the observational data and those inferred from the benchmark model simulation at t = 4 Gyr.
All Figures
Fig. 1 Histogram distributions of argument of perihelion (top left panel), longitude of ascending node (top right panel), perihelion longitude (bottom left panel), and inclination (bottom right panel) inferred from the welldetermined orbits of TNOs fulfilling the criteria 100 <a< 1500 au, q> 35 au (colour histograms). The orbits were discerned according to their a and divided into 3 groups: 100 <a< 150 au (yellow histogram), 150 <a< 250 au (blue histogram), and 250 <a< 1500 au (green histogram). The distributions inferred from all welldetermined orbits of ETNOs (150 <a< 1500 au, q> 30 au) are indicated as dashed bars. The blue and green bars in a given bin can, but do not have to, sum up to the height of a dashed bar in that bin, as the dashed bar also comprises objects with q in the range from 30 to 35 au. The histogram bars are aligned to the left, for example, in the bottom right panel, 0 indicates inclination range (0,10) deg, 10 indicates inclination range (10,20) deg, and so on. 

In the text 
Fig. 2 Histogram distributions of argument of perihelion (top left panel), longitude of ascending node (top right panel), perihelion longitude (bottom left panel), and inclination (bottom right panel) inferred from the welldetermined orbits of TNOS and Centaurs fulfilling the criteria 100 <a< 1500 au, q< 35 au (colour histograms). The distributions inferred from all welldetermined orbits of ETNOs (150 <a< 1500 au, q> 30 au) are indicated as dashed bars. 

In the text 
Fig. 3 Unbiased cumulative distribution function of the scattered KBOs (dashed line) vs. observed cumulative distribution of ETNOs (solid line), both as a function of inclination. The thick dashed line represents the bestfit model of Gulbis et al. (2010), while the thinner dashed lines are its 1σ variations. The dots correspond to the observed inclinations of ETNOs and are equidistantly spaced between 0 and 1 on the vertical axis. 

In the text 
Fig. 4 Evolution of test particles in L−e plane as inferred from the P9 simulation. The particles and P9 orbit the Sun in the same plane by design. The eccentricity range of the observed ETNOs is indicated by the dashed horizontal lines. Two transparency levels are used as a proxy for two different time windows of the simulation. The less transparent evolutionary tracks represent longliving 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. Three starting semimajor axes were used as follows: a = 150 au (first column), a = 350 au (second column), and a = 550 au (third column). Nbody encounters are completely omitted in this case, hence semimajor axes do not change during the simulation. The secular effect of the giant planets is modelled with a strong J_{2} moment of the Sun. There were no stable orbits at a = 550 au level. 

In the text 
Fig. 5 Evolution of test particles in ω−e (first row), Ω−e (second row), and L−e (third row) plane as inferred from the benchmark model simulation with Q_{2} = 1. 

In the text 
Fig. 6 Evolution of test particles in ω−e (first row), Ω−e (second row), and L−e (third row) plane as inferred from the benchmark model simulation with negative Q_{2}, Q_{2} = −1. 

In the text 
Fig. 7 a−e diagram from the last Gyr of the benchmark model simulation incorporating gravitational interaction with Neptune. Q_{2} = 1 (top left panel), Q_{2} = 2 (top right panel), Q_{2} = 4 (bottom left panel), and Q_{2} = 6 (bottom right panel) were assumed. Two colours (transparency levels) are used as a proxy for two inclination ranges. The blue (less transparent) evolutionary tracks represent orbits with inclinations spanning i ∈ (10,30) deg, the interval the inclinations of Sednoids lie in. The grey (more transparent) tracks refer to all the other longliving orbits. Sedna and 2012 VP 113 are indicated as large dots. The dashed lines indicate q = 30, 50, 100, and 550 au levels (from right/bottom to left/upwards). Objects between the two leftmost/uppermost lines (q in range between 100 and 550 au) would hardly be observed if they really exist. 

In the text 
Fig. 8 a−e diagram from the last Gyr of the benchmark model simulation incorporating gravitational interaction with Neptune. The enhanced gravity effect is accounted for in the model using the interpolating function . yields Q_{2} = 5.9. Only the evolution of test particles in orbits with inclinations spanning i ∈ (10,30) deg, which is the interval that the inclinations of Sednoids lie in, is shown. Sedna and 2012 VP 113 are indicated as large dots. The dashed lines indicate q = 30, 50, 100, and 550 au levels (from right/bottom to left/upwards). 

In the text 
Fig. 9 a−e diagram from the last Gyr of the P9 simulation incorporating gravitational interaction with Neptune. Only evolution of test particles in orbits with inclinations spanning i ∈ (10,30) deg, which is the interval that the inclinations of Sednoids lie in, is shown. Sedna and 2012 VP 113 are indicated as large dots. The dashed lines indicate q = 30, 50, 100, and 550 au levels (from right/bottom to left/upwards). 

In the text 
Fig. 10 a−e diagram from the last Gyr of the benchmark model simulation incorporating gravitational interaction with Neptune. Negative values of Q_{2}, Q_{2} = −1 (left panel) and Q_{2} = −4 (right panel) were used. Two colours (trasparency levels) are used as a proxy for two inclination ranges. The blue (less transparent) evolutionary tracks represent orbits with inclinations spanning i ∈ (10,30) deg, which is the interval that the inclinations of Sednoids lie in. The grey (more transparent) tracks refer to all the other longliving orbits. Sedna and 2012 VP 113 are indicated as large dots. The dashed lines indicate q = 30, 50, 100, and 550 au levels (from right/bottom to left/upwards). 

In the text 
Fig. 11 Evolution of test particle orbits, taken from the P9 simulation, in ω−a plane. Only particles surviving at least 3 Gyr are taken into account. The evolutionary tracks of the particles are shown in the last Gyr of the simulation every time their orbits met the criteria: i< 50 deg, q< 80 au. The tracks below a = 200 au are not shown as these are just a forest of fullwidth horizontal lines. The X symbols indicate the positions of potentially unstable particles, which are discarded at some time between t = 3 and 4 Gyr, at the time of their destruction. The plus signs indicate positions of longliving particles at t = 4 Gyr. The known ETNOs with q> 35 au (circles) and q< 35 au (triangles) are plotted. TNO 2003 SS_{422} is denoted as a square and its large uncertainty in a is indicated by the error bar. Shaded regions indicate the observed clustering regions for a> 250 au ETNOs. 

In the text 
Fig. 12 Evolution of test particle orbits, taken from the P9 simulation, in Ω−a (left panel) and L−a plane (right panel). 

In the text 
Fig. 13 Evolution of test particle orbits, taken from the benchmark model simulation, in ω−a plane. The strength of EFE was assumed to be Q_{2} = 1. Only particles surviving at least 3 Gyr are taken into account. The evolutionary tracks of the particles are shown in the last Gyr of the simulation every time their orbits met the criteria, i.e. i< 50 deg, q< 80 au. The tracks below a = 500 au are not shown as these are just fullwidth horizontal lines. 

In the text 
Fig. 14 Evolution of test particle orbits, taken from the benchmark model simulation, in Ω−a (left panel) and L−a (right panel) plane. The strength of EFE was assumed to be Q_{2} = 2. Only particles surviving at least 3 Gyr are taken into account. Evolutionary tracks of the particles are shown in the last Gyr of the simulation every time their orbits met the criteria, i.e. i< 50 deg, q< 80 au. The tracks below a = 400 au are not shown as these are just fullwidth horizontal lines. Increasing Q_{2}, the clustering in L sets in at lower semimajor axes, e.g. with Q_{2} = 6; similar clustering as in the figure develops for a> 300 au. 

In the text 
Fig. 15 Histogram distributions of ω (top left panel), Ω (top right panel), L (bottom left panel), and i (bottom right panel) of test particle orbits at the end of the benchmark model simulation, i.e. at t = 4 Gyr. The strength of EFE was assumed to be Q_{2} = 1. Orbits are discerned according to their semimajor axis and divided into three groups: 150 <a< 250 au (blue histograms), 250 <a< 550 au (green histograms), and 550 <a< 1500 au (red histograms). The observed distributions for ETNOs with 250 <a< 1500 au, multiplied by factor 4, are shown as well (dashed histograms). Histograms are aligned to the left. 

In the text 
Fig. 16 Histogram distributions of ω (top left panel), Ω (top right panel), L (bottom left panel), and i (bottom right panel) of test particle orbits at the end of the benchmark model simulation, i.e. at t = 4 Gyr. The strength of EFE was assumed to be Q_{2} = 4. 

In the text 
Fig. 17 Perihelion distance as a function of time for 4 selected orbits from the benchmark model simulation characterised by highq amplitude. Assumed values of Q_{2} and initial (final) values of a and q are indicated. 

In the text 
Fig. 18 Evolution of test particle orbits, taken from the benchmark model simulation in Ω−a (left panel) and L−a (right panel) plane. A negative Q_{2}, Q_{2} = −2 was assumed. Only particles surviving at least 3 Gyr are taken into account. Evolutionary tracks of the particles are shown in the last Gyr of the simulation every time their orbits met the criteria, i.e. i < 50 deg, q< 80 au. The tracks below a = 400 au are not shown as these are just fullwidth horizontal lines. Decreasing Q_{2}, Q_{2} < 0, the clustering in L sets in at lower semimajor axes, e.g. with Q_{2} = −6 similar clustering as in the figure develops for a > 250 au. 

In the text 
Fig. 19 Histogram distributions of ω (top left panel), Ω (top right panel), L (bottom left panel), and i (bottom right panel) of test particle orbits at the end of the benchmark model simulation, i.e. at t = 4 Gyr. A negative Q_{2}, Q_{2} = −1 was assumed. Orbits are discerned according to their semimajor axis and divided into three groups: 150 < a < 250 au (blue histograms), 250 < a < 550 au (green histograms), and 550 < a < 1500 au (red histograms). The observed distributions for ETNOs with 250 < a < 1500 au, multiplied by factor 4, are shown as well (dashed histograms). Histograms are aligned to the left. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.