Free Access
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/0004-6361/201630335
Published online 30 June 2017

© 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 non-existence of Vulcan finally led to the development of new revolutionary ideas in physics – the geometrical interpretation of gravity – explaining precisely the minute anomalous precession of Mercury (Le Verrier 1859; Einstein 1916).

Orbital characteristics of observed extreme trans-Neptunian objects (ETNOs), objects having perihelion distances, q, greater than the orbit of Neptune, and semimajor axes, a, greater than 150 au and less than roughly 1500 au (definition established and motivated in Trujillo & Sheppard 2014), suggest that these objects have been systematically perturbed by some external force. On September 1, 2016, the Minor Planet Center (MPC) orbit database registered 19 ETNOs1. The outlying positions of 90377 Sedna and 2012 VP113 (hereafter, collectively referred to as Sednoids) in the q versus eccentricity, e, diagram (e.g. Fig. 2 in Sheppard & Trujillo 2016) cannot be addressed by the interaction with the known solar system planets and the current Galactic environment (e.g. Dones et al. 2015). Nevertheless, orbits of ETNOs are also characterised by argument of perihelion, ω, longitude of the ascending node, Ω, and perihelion longitude, L,2 which are distributed in a highly non-uniform and clustered fashion (Trujillo & Sheppard 2014; de la Fuente Marcos & de la Fuente Marcos 2014; Brown & Firth 2016; Batygin & Brown 2016). This is an unexpected observation, as these orbital elements are expected to circulate rapidly under the action of the giant planets.

Tightly clustered are also the eccentricity, e (e = 0.82 ± 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.80.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 inclinations3 (de la Fuente Marcos & de la Fuente Marcos 2014).

Batygin & Brown (2016) recognised that orbits of long-periodic ETNOs with a> 250 au (six such objects were known at that time4: Sedna, 2004 VN112, 2007 TG422, 2010 GB174, 2012 VP113, and 2013 RF98) 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 SR349 and 2013 FT28). These two ETNOs have ω and i clustered in line with the original six, but one of them, 2013 FT28, has Ω and L with 90180 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 present-day 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 Tully-Fisher relation (McGaugh et al. 2000), the relation between stellar and dynamical central surface densities (Lelli et al. 2016), and the mass discrepancy-acceleration correlation (Sanders 1990; McGaugh 2004; Desmond 2017). All these can be encapsulated into one law called the radial acceleration relation (McGaugh et al. 2016; Lelli et al. 2017) – a very tight correlation between the radial acceleration traced by rotation curves and the acceleration predicted by the observed distribution of baryons. There is no reference to dark matter in these laws, even in systems that are supposedly dark matter dominated, i.e. systems showing large mass-discrepancy in Newtonian dynamics. Additionally, “for any feature in the luminosity profile there is a corresponding feature in the rotation curve and vice versa” (Sancisi 2004). In the context of dark matter halos, the existence of these laws means that there is a very tight correlation between the baryonic and dark matter distributions. Such a very tight correlation, independent of properties and histories of individual galaxies, is unexpected in the current cosmological paradigm: the ΛCDM model, where Λ is the positive cosmological constant representing dark energy, whose nature is unknown, and CDM stands for cold (non-baryonic) dark matter. All these laws were predicted a priori and follow naturally from the basic principles of Milgromian dynamics (Milgrom 2014, 2016b).

Milgromian dynamics (MD or MOND = modified Newtonian dynamics; Milgrom 1983b), also known as scale-invariant dynamics, uniquely predicts how the distribution of stars and gas in a galaxy relates to its dynamics. For reviews on MD, see for example the Famaey & McGaugh (2012), which is thorough and technical, or Milgrom (2016a), which contains a discussion on the potential connection of MD to fundamental physics. Milgromian dynamics departs from Newtonian dynamics (and general relativity) at low accelerations with respect to some predefined frame (this could be e.g. the cosmic microwave background rest frame), and the definition of low acceleration is determined by a new constant of nature, a0, having units of acceleration. In the low acceleration regime, unlike Newtonian dynamics and general relativity, MD becomes space-time 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 a0 in the data.

In some versions of MD (typically in modified gravity versions of MD), the dynamics of a small subsystem, freely falling in the field of a large external system, is affected by the external free-fall acceleration through the so-called external field effect (EFE; Milgrom 1983b; Famaey & McGaugh 2012). This is the case of, for instance, a dwarf galaxy in the field of a host galaxy, or the solar system falling freely in the field of the Galaxy. The EFE arises from generic non-linearity of equations of MD and from the fact that it is based on acceleration. McGaugh & Milgrom (2013a) used MD to predict velocity dispersions of dwarf satellite galaxies of the Andromeda galaxy. Later, in McGaugh & Milgrom (2013b), they compared their a priori predictions with newly available data. They reported that, “The data are in good agreement with our specific predictions for each dwarf made a priori... MOND distinguishes between regimes where the internal field of the dwarf, or the external field of the host, dominates. The data appear to recognise this distinction, which is a unique feature of MOND not explicable in Λ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 a0), the quasi-Newtonian behaviour is restored (Famaey & McGaugh 2012).

Recently, McGaugh (2016) predicted velocity dispersion, σvel, of Crater II in MD. Although internal accelerations are very low (≪ a0), 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 (Boylan-Kolchin 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 produce5 enough torque to deliver primordial solar system bodies to Sedna-like orbits even in the present-day Galactic environment. Here we reassessed this problem and checked whether the perihelia of both Sednoids can be acquired by the action of EFE. We used recent constraints available on the strength of EFE in the ETNOs region coming from the radio tracking of the Cassini spacecraft (Blanchet & Novak 2011; Hees et al. 2014, 2016). Assuming that the statistics on the orbital elements of ETNOs are representative of the real population, we also investigated whether the spatial orientation of the orbits of ETNOS can be explained by this effect of MD as well.

We further discuss the observational data on ETNOs in Sect. 2. Milgromian dynamics and the dynamical effect it has on the solar system are introduced in Sect. 3. We present a dynamical model of the solar system obeying laws of MD that are applicable in the ETNOs region in Sect. 4. The comparison between the P9 scenario in Newtonian dynamics and the MD scenario without an additional planet is made in 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 q-filter removes five objects (2013 UH15, 2010 VZ98, 2001 FP185, 2013 FS28, and 2015 SO20) from the ETNOs sample. One of these objects, 2013 FS28, is an outlier in terms of ω, 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 well-determined orbits of ETNOs (dashed histogram); ETNOs were defined in Trujillo & Sheppard (2014) as objects with 150 <a< 1500 au, q> 30 au. We have to bear in mind that this choice of the boundary q may be just our means of forcing the data to behave as we expect.

The two q-populations are distinct in terms of ω 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 SS422 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.

thumbnail 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 well-determined 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 well-determined 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.

thumbnail 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 well-determined orbits of TNOS and Centaurs fulfilling the criteria 100 <a< 1500 au, q< 35 au (colour histograms). The distributions inferred from all well-determined orbits of ETNOs (150 <a< 1500 au, q> 30 au) are indicated as dashed bars.

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

On September 1, 2016, the MPC orbit database listed 40 TNOs with a> 100 au, 19 ETNOs, and 8 long-periodic ETNOs with a> 250 au. When we restrict ourselves to q> 35 au TNOs only, these numbers reduce to 32 TNOs with a> 100 au, 14 ETNOs, and the same number of 8 long-periodic ETNOs. Clearly, we are dealing with small-number statistics. However, as pointed out elsewhere (Brown & Firth 2016; Batygin & Brown 2016), it is unlikely that the observed clustering is purely due to chance. Assuming the ω, Ω, 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, 106 times, and calculated the circular variance of each set, σi2=1R\hbox{$\sigma^{2}_{i}=1-R$}, R=(S2+C2)1/2\hbox{$R=(\overline{S}^{2}+\overline{C}^{2})^{1/2}$}, S=(1/8)n=18sinθn\hbox{$\overline{S}=\left(1/8\right)\sum_{n=1}^{8}\sin\theta_{n}$}, C=(1/8)n=18cosθn\hbox{$\overline{C}=\left(1/8\right)\sum_{n=1}^{8}\cos\theta_{n}$}. The probability is given as a frequency of sets for which σi2\hbox{$\sigma^{2}_{i}$} 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 off-ecliptic 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.

thumbnail 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 best-fit 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 survey6, 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 oddly7 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 KBOs8. 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 19.1+3.9-3.6\hbox{$19.1\substack{+3.9 \\ -3.6}$} deg, and with a width of 6.9+4.1-2.7\hbox{$6.9\substack{+4.1 \\ -2.7}$} 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 = gN but fromg=ν(gN/a0)gN,\begin{eqnarray} \label{MOND_basic} {\vec{g}} = \nu( g^{\rm N}/a_{0}){\vec{g}}^{\rm N}, \end{eqnarray}(1)where gN is expected Newtonian gravitational acceleration, gN is its norm, ν(y), ygN/a0 is a modification factor, the so-called interpolating function, fulfilling behaviour ν(y) → 1 for y ≫ 1 (Newtonian limit) and ν(y) → y− 1/2 for y ≪ 1 (deep-MD limit), and a0 ~ 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) galaxies9, 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) ·[μ(|∇Φ|a0)∇Φ]=4πGρ=2ΦN,\begin{eqnarray} \label{AQUAL} \nabla\cdot\left[\mu\left(\frac{\vert\nabla\Phi\vert}{a_{0}}\right)\nabla\Phi\right]~=~4\pi G\rho~=~\nabla^{2}\Phi_{\rm N}, \end{eqnarray}(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 quasi-linear MOND (QUMOND), states that the governing gravitational potential Φ is given by (Milgrom 2010) 2Φ=·[ν(|ΦN|a0)ΦN].\begin{eqnarray} \label{QUMOND} \nabla^{2}\Phi~=~\nabla\cdot\left[\nu\left(\frac{\vert\nabla\Phi_{\rm N}\vert}{a_{0}}\right)\nabla\Phi_{\rm N}\right]. \end{eqnarray}(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 full-fledged 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 rrM, where rM ≡ (GM/a0)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 = ggN, can be written as δg=[ν(gN/a0)1]gN.\begin{eqnarray} \label{mond_effect1} \delta {\vec{g}} = \left[\nu(g^{\rm N}/a_{0}) - 1\right]{\vec{g}}^{\rm N}. \end{eqnarray}(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 a0, 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 MD10 and it comes from the non-linearity of the modified Poisson equation. An embedding external field that is constant and uniform, with magnitude ge, gea0, 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 gev02/R02402\hbox{$g_{\rm e}\approx v^{2}_{0}/R_{0}\approx 240^{2}$} km2 s-2/8.0 kpc ≐ 2.0 × 10-10 m s-2 (for values of the Galactic parameters v0 and R0, we refer to McMillan & Binney 2010; Schönrich 2012). Close to the Sun, at distances rrM, rM ≡ (GM/a0)1/2 ~ 7000 au, where rM 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 Φ=ΦNQ22rirj(eiej13δij)=ΦN+δΦ,\begin{eqnarray} \label{mond_effect2} \Phi &=& \Phi_{\rm N} - \frac{Q_{2}}{2} r^{i}r^{j}\left(e_{i}e_{j}-\frac{1}{3}\delta_{ij}\right)\nonumber \\ &=& \Phi_{\rm N} + \delta\Phi, \end{eqnarray}(5)where ege/ge is a unitary vector pointing towards the Galactic centre11 (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 Q2, which depends on the specific modified Poisson formulation, form of the interpolating function, and strength of the external field in natural units a0 (Milgrom 2009a; Blanchet & Novak 2011). If not stated otherwise, Q2 is assumed to be positive. All routinely used forms of the interpolating function yield12Q2> 0 (Milgrom 2009a; Blanchet & Novak 2011). If Q2 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 Q2, 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 Q2 close to the theoretical maximum allowed by the Earth-Cassini spacecraft range data (hereafter Cassini data), Q2 ~ 10-27 s-2 (Hees et al. 2014), then EFE is ~ 103 times stronger than the Galactic tides (Hees et al. 2014). Instantaneous action of EFE is the same as that of P9 lying in the direction of the GC-anticentre, characterised by the tidal parameter GM′/r′3 = Q2/ 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 present-day Galactic anticentre (the preferred direction in Eq. (5)). Comparison between orbital characteristics of ETNOs induced on long timescales by EFE and those induced by P9 is a topic of this study.

  • Asphericity effect

This effect appears only in modified gravity formulations of MD. It affects every system with non-spherically symmetric matter distribution. Suppose we have an isolated high-acceleration system S, which assumes μ = ν = 1 to the desired accuracy everywhere within S, of total mass M and extension R, RrM, where rM=GM/a0\hbox{$r_{\rm M}=\sqrt{GM/a_{0}}$}. Suppose also that mass distribution within S, ρ(r), varies on timescales much larger than rM/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>rM. 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 Φ=ΦNαGrM5rirjQij=ΦN+δΦ,Qij=12ρ(r)(r2δij3rirj)d3r,\begin{eqnarray} \label{mond_effect3} &&\Phi = \Phi_{\rm N} - \frac{\alpha G}{r^{5}_{\rm M}} r^{i}r^{j}Q_{ij}\nonumber \\ && \ \ \ = \Phi_{\rm N} + \delta\Phi,\nonumber \\ &&Q_{ij} = \frac{1}{2}\int\rho({\vec{r'}})\left(r'^{2}\delta_{ij}-3r'_{i}r'_{j}\right){\rm d}^{3}r', \end{eqnarray}(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 system13, MM and Qij 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 gN ~ a0, but beyond gN ~ κa0, α 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, a0 and a0κa0\hbox{$a^{*}_{0}\equiv\kappa a_{0}$}, 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 constraints14 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 regime15. Referring to Table 2 in Hees et al. (2016), we chose νn(y)=[1exp(yn)]1/(2n)+[11/(2n)]exp(yn)\hbox{$\overline\nu_{n}(y)=\left[1-\exp(-y^{n})\right]^{-1/(2n)}+\left[1-1/(2n)\right]\exp(-y^{n})$} interpolating function family as most representative. The Cassini data restrict νn\hbox{$\overline\nu_{n}$} to n ≥ 2. The enhanced gravity effect is completely negligible in the whole solar system16 for νn\hbox{$\overline\nu_{n}$}, n ≳ 3. Thus, in the first approximation, we neglect the enhanced gravity effect in the benchmark model.

The parameter Q2 of Eq. (5) is position dependent in general (Milgrom 2009a). However, as shown by Blanchet & Novak (2011) in the framework of AQUAL, Q2 can be considered a constant in a large sphere centring on the Sun. The variation of Q2 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 Q2, Q2 = Q2(0). All higher multipoles are omitted. The value of Q2 was constrained by the Cassini data to be Q2 = 3 ± 3 × 10-27 s-2 (nominal ± 1σ value) (Hees et al. 2014). From now on, if we refer to Q2 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 exg=+cos(wt+τ)eyg=sin(wt+τ)ezg=0ex=0.054876exg+0.494109eygey=0.993821exg0.110991eygez=0.096477exg+0.862286eyg,\begin{eqnarray} \label{unit_vector} &&e_{x-g} = +\cos(w t + \tau) \nonumber \\ &&e_{y-g} = -\sin(w t + \tau)\nonumber \\ &&e_{z-g} = 0 \nonumber \\ &&e_{x} = -~0.054876~e_{x-g} + 0.494109~e_{y-g}\nonumber \\ &&e_{y} = -~0.993821~e_{x-g} - 0.110991~e_{y-g}\nonumber \\ &&e_{z} = -~0.096477~e_{x-g} + 0.862286~e_{y-g}, \end{eqnarray}(7)where the first three lines are components of e in the Galactic coordinates17 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 GC18, lying in the Galactic mid-plane, 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 long-living orbital confinement in the solar system ruled by the dynamical laws of MD. For this purpose we made use of the benchmark model. Additionally, we model the secular effect of the giant planets by considering the Sun of radius R = aU, where aU is the semimajor axis of Uranus and J2 is the moment of magnitude J2=imiai2/(2MR2),\hbox{$J_{2}=\sum_{i}m_{i}a^{2}_{i}/(2M R^{2}),$} and where we sum over the giant planets with masses mi, semimajor axes ai, 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 Le plane, considering Q2 taken from a set Q2 = { 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 Q2.

thumbnail Fig. 4

Evolution of test particles in Le 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 long-living particles in the last Gyr of the simulation. The more transparent tracks refer to evolution in the first 3 Gyr of the simulation and also comprise unstable orbits. Three starting semimajor axes were used as follows: a = 150 au (first column), a = 350 au (second column), and a = 550 au (third column). N-body 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 J2 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 present-day GC-anticentre. The ecliptic latitude of the GC is only minus few degrees. Here our intention is to demonstrate the existence of the long-living libration regions; in the next section we examine development of the orbital clustering, starting with completely randomised orbits.

Every orbital integration in this paper was carried out using mercury6 N-body integration software package (Chambers 1999). The hybrid symplectic-Bulirsch-Stoer algorithm with a timestep equal to a tenth of the orbital period of Neptune was employed. Phase τ in Eq. (7) was chosen in a way that at t = 4 Gyr the vector e points towards the present-day 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 simulation20. 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 = aU, 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 LL equal to 0 (20 particles for each a) and 180 deg (20 particles for each a), and nil inclinations21. 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.

thumbnail Fig. 5

Evolution of test particles in ωe (first row), Ω−e (second row), and Le (third row) plane as inferred from the benchmark model simulation with Q2 = 1.

thumbnail Fig. 6

Evolution of test particles in ωe (first row), Ω−e (second row), and Le (third row) plane as inferred from the benchmark model simulation with negative Q2, Q2 = −1.

We depict evolution of test particles in Le plane as inferred from the first P9 simulation in Fig. 4. All particles meet the condition ii′ = 0 by design. The discussed orbital elements always refer to the heliocentric ecliptic J2000 reference frame. Two transparency levels are used as a proxy for two different time windows of the simulation. The less transparent evolutionary tracks represent long-living particles in the last Gyr of the simulation. The more transparent tracks refer to evolution in the first 3 Gyr of the simulation and also comprise unstable orbits. Orbits with a = 150 au typically experience trivial apsidal circulation. Only one small libration island exists at L ≈ 240 deg (aligned with P9) for near-circular orbits. For a = 350 au, two long-term stable libration regimes exist (Batygin & Brown 2016): first, the resonant dynamics regime concerning eccentric orbits (with eccentricities similar to those of the observed ETNOs) and centred at L ≈ 60 deg and anti-aligned 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 Le 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 long-term stable orbits.

The picture in MD is rich as well. Orbits with a = 150 au are characterised by trivial circulation of the angles (L, Ω,ω). No stable libration develops at this a-level for realistically high value of Q2. 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 Q2 = 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 Q2. In the case of lower eccentricities, we get a circulation for models with Q2 ≲ 4, and neither circulation nor strong confinement using models with Q2 ≳ 4. The specialty of values L = 90 and 270 deg is still noticeable.

The case Q2 = 1, a = 550 au is especially interesting. The motion in Le 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 Le plane as inferred from the benchmark model simulations using Q2 = 1 in Fig. 5.

Despite the rich structure in ωe, Ω−e, and Le planes in the MD, the simulated evolutionary paths are in sharp contrast with present-day observations. The simulations predict that the vast majority of ETNOs should have Ω ∈ (180,360) deg. Hence, the vast majority of ETNOs should have ascending node lying in the same half-space, with the GC lying in the same half-space. This is not an observed trend. Of 19 ETNOs, 14 have Ω ∈ (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 Q2). When Q2 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 that22 although Q2 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 Q2 = 3 ± 3 (nominal ± 1σ value).

We substituted Q2 → −Q2 and repeated the benchmark model simulations. This means we considered Q2s that are 13σ off of the nominal value of Hees et al. (2014), with corresponding probabilities for Q2 ranging from 16 (Q2< 0) down to 0.1% (Q2< −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 Le planes without qualitative change. This can be seen in the case Q2 = −1 in Fig. 6. Looking at EER at a = 550 au, Ω lies sometimes preferentially in range (0,180) deg (Q2 = −1), sometimes in range (180,360) deg (Q2 = −4 and − 6), and sometimes there is preference for neither of the two regions (Q2 = −2).

6. Add Neptune and stir

Neptune may have influenced the distribution of the orbital elements of ETNOs through orbital filtering (e.g. by preferentially ejecting ETNOs in specific orbits) and resonant effects. As a next step, we added Neptune orbiting the Sun to the dynamical model. We consider the gravitation of the Sun and planet Neptune in the usual Newtonian way as a direct N-body interaction. The secular effect of the three remaining giant planets is modelled by considering the Sun of radius R = aU, where aU is the semimajor axis of Uranus, and J2 moment of magnitude J2=imiai2/(2MR2)\hbox{$J_{2}=\sum_{i}m_{i}a^{2}_{i}/(2M R^{2})$}, where we sum over the remaining giant planets with masses mi, semimajor axes ai, 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 Q2 taken from a set Q2 = { 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 Q2 = 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 same23a, 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 VP113 (Trujillo & Sheppard 2014) was surprising as their orbits did not fit into the standard picture of the solar system, settled in the current Galactic environment, containing a set of known planetary bodies. The literature discussed scenarios in which these orbits were shaped in environments that, in the past, were different from the present-day environment, for instance the possibility that the Sun was born in a star cluster (Kaib & Quinn 2008; Brasser et al. 2012) and/or was traveling much closer to the GC in the past (Kaib et al. 2011).

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

Evolution in ae plane is depicted during the last Gyr of the benchmark model simulation in Fig. 7 for various values of Q2. Figures showing ae diagrams are always based on the simulations starting with 800 particles. With Q2 = 0 and no P9, the particles would be stuck between the two rightmost/bottom-most dashed lines (q in range between 30 and 50 au). Moreover, assuming Q2 = 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 Q2 to Q2> 0.5. The orbit of 2012 VP113, with similar perihelion distance but much lower semimajor axis than Sedna, requires Q2> 4, if we demand that EFE is solely responsible for the elevation of perihelia. Even with Q2 = 4−6, 2012 VP113 still lies at the low-end border of the achievable eccentricities at a ≈ 250 au. Objects with higher eccentricities are easier to detect. The orbit of 2012 VP113 is thus only marginally compatible with the EFE origin scenario.

6.1.1. Beyond the benchmark model

What effect can be expected from the next lowest-order correction (octupole) neglected in Eq. (5) and in the benchmark model? This correction (adds to RHS of Eq. (5)) can be written as24 (Blanchet & Novak 2011) δΦoct=Q36rirjrk[eiejek15(eiδjk+ejδik+ekδij)].\begin{eqnarray} \label{octupole} \delta\Phi_{\rm oct} &=& \frac{Q_{3}}{6} r^{i}r^{j}r^{k}\left[e_{i}e_{j}e_{k}-\frac{1}{5}\left( e_{i}\delta_{jk}+e_{j}\delta_{ik}+e_{k}\delta_{ij}\right)\right]. \end{eqnarray}(8)The octupole strength parameter Q3 (typically Q3< 0) depends on the specific modified Poisson formulation, form of the interpolating function, and strength of the external field in natural units a0. The parameter Q3 is position dependent in general. The variation of Q3 within our volume of interest is up to several percent (Blanchet & Novak 2011). We treat Q3 as position independent, Q3 = Q3(0).

Blanchet & Novak (2011) tested some of the commonly used interpolating functions μ, for instance the μn(x) = x/ (1 + xn)1 /n family, in the framework of AQUAL, and found that | Q3 | < 1.2 × 10-40 m-1 s-2 close to the Sun. From now on, when we refer to Q3 without explicitly stating its units, the unit is 10-40 m-1 s-2. The highest | Q3 | 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 Q2 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 | Q3 | ~ 10-2 (Blanchet & Novak 2011). We found that with such a low value of Q3, the incorporation of the octupole term, Eq. (8), has a negligible effect. We note that the most stringent constraints on Q2 comes from Hees et al. (2016), where QUMOND framework was applied. The values of Q3 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).

thumbnail Fig. 7

ae diagram from the last Gyr of the benchmark model simulation incorporating gravitational interaction with Neptune. Q2 = 1 (top left panel), Q2 = 2 (top right panel), Q2 = 4 (bottom left panel), and Q2 = 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 long-living 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 left-most/upper-most lines (q in range between 100 and 550 au) would hardly be observed if they really exist.

thumbnail Fig. 8

ae 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 ν2\hbox{$\overline{\nu}_{2}$}. ν2\hbox{$\overline{\nu}_{2}$} yields Q2 = 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).

thumbnail Fig. 9

ae 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 νn\hbox{$\overline{\nu}_{n}$} family and QUMOND, ν2\hbox{$\overline{\nu}_{2}$} yields Q2 = 5.9 when one assumes ηge/a0 = 2.3, a0 = 0.815 × 10-10 m s-2 (Hees et al. 2016). The quoted value of a0 comes from rotation curves fits using ν2\hbox{$\overline{\nu}_{2}$}. With ν2\hbox{$\overline{\nu}_{2}$} one gets, besides EFE, as well as a non-negligible enhanced gravity effect. We incorporated the enhanced gravity effect, Eq. (4), into the model and repeated the simulation. We used ν2(y)\hbox{$\overline{\nu}_{2}(y)$}, y=|gN+geN|/a0\hbox{$y=\vert {\vec{g}}^{\rm N}_{\odot}+{\vec{g}}^{\rm N}_{\rm e}\vert/a_{0}$}, where gN\hbox{${\vec{g}}^{\rm N}_{\odot}$} is Newtonian gravitation from the Sun, and geN\hbox{${\vec{g}}^{\rm N}_{\rm e}$} was found from ge=ν2(|geN|/a0)geN\hbox{${\vec{g}}_{\rm e}=\overline{\nu}_{2}(\vert{\vec{g}}^{\rm N}_{\rm e}\vert/a_{0}){\vec{g}}^{\rm N}_{\rm e}$} knowing ge. Direction of fields ge and geN\hbox{${\vec{g}}^{\rm N}_{\rm e}$} rotate as prescribed by Eq. (7). An outcome is in Fig. 8. 2012 VP113 still lies at the low-end border of the achievable eccentricities at a ≈ 250 au. We remind that νn\hbox{$\overline{\nu}_{n}$}, n< 2, is disfavoured by the Cassini data (Hees et al. 2016), and for νn\hbox{$\overline{\nu}_{n}$}, 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 (Q2 = 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.

thumbnail Fig. 10

ae diagram from the last Gyr of the benchmark model simulation incorporating gravitational interaction with Neptune. Negative values of Q2, Q2 = −1 (left panel) and Q2 = −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 long-living 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 Q2

We show ae diagrams that are analogous to the previous diagrams, but now with negative values of Q2, Q2 = −1, and Q2 = −4, in Fig. 10. The orbits of both Sednoids are produced in the MD model as long as | Q2 | > 4. With Q2 = −4, 2012 VP113 lies at the low-end 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 Kolmogorov-Smirnov (KS) test25 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 Q2. Additionally, we also did the KS test based on results from Q2 = 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.

Table 1

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), p-values, 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 (105), interacting gravitationally with five planets (known giant planets + P9) in a full N-body manner for 4 Gyr, leading to no angular clustering at the end of the simulation.

In the semimajor axis range 150 <a< 250 au, the MD model test yields p-values, p(150−250 | 150−250), similar to the control test, except for the distribution of Ω when even lower p-values 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 Q2 ≳ 4. However, the compatibility between the observational data and the MD model remains low in the case of ω and L, though p-values, p( > 250 | > 250), are slightly higher than in the control test if, e.g. Q2 = 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.

Table 2

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

Table 3

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 La plane in the last Gyr of the Newtonian simulation with P9. We initialised P9 in the orbit (#1). Initialising P9 in the orbit (#2) would have yielded a similar picture. Evolutionary tracks of the particles are shown every time their orbits met the criteria, i.e. i< 50 deg, q< 80 au. The evolutionary tracks continue down to a = 100 au, but we do not show them below a = 200 au as they are just a forest of full-width horizontal lines. Clustering in ω 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.

thumbnail 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 full-width 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 long-living particles at t = 4 Gyr. The known ETNOs with q> 35 au (circles) and q< 35 au (triangles) are plotted. TNO 2003 SS422 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.

thumbnail Fig. 12

Evolution of test particle orbits, taken from the P9 simulation, in Ω−a (left panel) and La 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 Q2 = 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 Q2 = 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 Q2 ≳ 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 Q2 = 2 simulation. Strong confinement in Ω develops for stable orbits in a wide range of Q2 values, Q2 ∈ (0.5,6), and also for potentially unstable orbits if Q2 ≳ 2.

thumbnail Fig. 13

Evolution of test particle orbits, taken from the benchmark model simulation, in ωa plane. The strength of EFE was assumed to be Q2 = 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 full-width horizontal lines.

The distribution of L shows a coherent picture for Q2 in the range Q2 ∈ (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 La plane, as taken from Q2 = 2 simulation. As Q2 increases, clustering develops at lower and lower semimajor axes. For example, for Q2 = 1, the clustering develops only beyond a = 500 au, while the border is at a = 300 au for Q2 = 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 (Q2 = 1) and 16 (Q2 = 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 Q2 values, Q2 ∈ (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 Q2, 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 Q2 ≈ 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 Q2, one gets rid of low-inclination ETNOs gradually. Objects with higher inclinations are more resilient. With Q2 = 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.

thumbnail Fig. 14

Evolution of test particle orbits, taken from the benchmark model simulation, in Ω−a (left panel) and La (right panel) plane. The strength of EFE was assumed to be Q2 = 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 full-width horizontal lines. Increasing Q2, the clustering in L sets in at lower semimajor axes, e.g. with Q2 = 6; similar clustering as in the figure develops for a> 300 au.

thumbnail 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 Q2 = 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.

thumbnail 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 Q2 = 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 Q2 = 0.2 could lift the perihelion of P9 from the edge of the planetary region26 to few hundred au, while with P9 in the orbit (#1) and EFE of strength Q2 = 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 ae 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 Q3 ~ 10-2; the inclusion of the enhanced gravity effect in the case27Q2 = 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 Q2

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

Evolution in Ω−a and La plane, inferred from the last Gyr of Q2 = −2 simulation, is visualised in Fig. 18. Varying Q2 in the range Q2 ∈ (−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 Q2 = −1 beyond a = 450 au, while for Q2 = −6 already beyond a = 300 au. We found no apparent clustering in ω in our simulations with Q2< 0.

Histogram distributions of ω, Ω, L, and i inferred at t = 4 Gyr from the Q2 = −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 | Q2 | = 4−6, Q2< 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 Q2 values is highly improbable (probability < 1%) in light of the Cassini data.

thumbnail Fig. 17

Perihelion distance as a function of time for 4 selected orbits from the benchmark model simulation characterised by high-q amplitude. Assumed values of Q2 and initial (final) values of a and q are indicated.

Table 4

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 trans-Neptunian objects (ETNOs), have their orbits oriented in space in highly non-isotropic fashion. A general trend in numerical simulations is that any primordial clustering of ETNOs orbits is quickly washed away by the perturbations of giant planets. If the distributions of ETNOs angular orbital elements ω, Ω, 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 Galaxy28 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 Q2, and constrained by the Earth-Cassini range measurements to be Q2 = 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 Q2.

thumbnail Fig. 18

Evolution of test particle orbits, taken from the benchmark model simulation in Ω−a (left panel) and La (right panel) plane. A negative Q2, Q2 = −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 full-width horizontal lines. Decreasing Q2, Q2 < 0, the clustering in L sets in at lower semimajor axes, e.g. with Q2 = −6 similar clustering as in the figure develops for a > 250 au.

thumbnail 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 Q2, Q2 = −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| Q2 | ≳ 1 × 10-27 s-2.

  • The orbit of the other Sednoid, 2012 VP113, seems only marginally compatible with the EFE scenario of its origin; 2012 VP113 requires Q2> 6 × 10-27 s-2, which is more than 1σ off of the nominal value of Hees et al. (2014), or Q2< 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 a-range is bimodal (peaking around ω = 0 and 180 deg); moreover, the confinement fades away when Q2 is not close to Q2 = 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 Q2; if Q2< 0, the clustering centres close to 90 and 270 deg, while when Q2> 0, the clustering centres close to 0 (360) and 180 deg.

  • In the stable population, the clustering in L looms for | Q2 | ≳ 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 Q2 values, and traced by both stable and potentially unstable orbits with a> 150 au; depending on the sign of Q2, either Ω ∈ (0,180) deg (Q2< 0) or Ω ∈ (180,360) deg (Q2> 0) is preferred by the simulated ETNOs, especially when we look at a> 250 au orbits.

  • If | Q2 | ≳ 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 present-day 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 Q2. There are not any interpolating function families known, giving arise to Q2< 0 in AQUAL or QUMOND at the present time, though maybe these are possible to construct. Assuming Q2 to be negative solves the problem. For instance Q2 = −1 × 10-27 s-2 (1.3σ off of the nominal Q2 of Hees et al. 2014) yields clustering in Ω in line with the observations. The observed clustering in Ω (if it is of dynamical origin) constrains Q2 to Q2< 0. Perihelia of Sednoids and ETNOs distributions of L and i prefer | Q2 | > 4 × 10-27 s-2. The Q2s lower than Q2 = −4 × 10-27 s-2 (2.3σ off of the nominal Q2 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 Q2. Also, provided that P9 were found, Q2 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 ephemerides29. 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 Q2, 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.


1

All ETNOs have q> 33 au, 14 of 19 ETNOs have q> 35 au, and 8 of 19 have q> 40 au.

2

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.

3

We add to this that many surveys are performed preferentially close to the ecliptic.

4

Asteroid 2000 CR105 with a = 228 au is very close to the six in terms of ω, Ω, and L, however, asteroid 2001 FP185, with a = 226 au, is about 90 deg ahead in L. Hence, the boundary was identified to be at a = 250 au.

5

This depends strongly on the shape of MD interpolating function and the value of a0. 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).

6

These are spring and autumn preferentially. By inspecting dates of ETNOs discoveries, one finds that they were discovered primarily in spring and autumn.

7

This is only odd if there is no substantial external perturbation present.

8

ETNOs are a subset of the scattered KBOs.

9

Systems in deep-MD regime when isolated.

10

This effect does not appear in the modified inertia formulations of MD.

11

In AQUAL, e is pointing in the direction of the real (Milgromian) gravitational field of the Galaxy. In QUMOND, e is pointing in the direction of the Newtonian gravitational field. For the approximately axially symmetric Galaxy, these directions are nearly the same.

12

In Milgrom (2009a), a dimensionless parameter q, q=2Q2(GM)1/2/(3a03/2)\hbox{$q=-2Q_{2}(GM_{\odot})^{1/2}/(3a^{3/2}_{0})$}, was employed instead of Q2.

13

Solar system is not an isolated system, so we have a superposition of EFE and the asphericity effect.

14

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).

15

Also, these are the functions for which the enhanced gravity effect is negligible in the planetary region, compared to EFE.

16

Except for small regions, where internal and external gravitational fields add up into a vector of small magnitude.

17

Principal axis points from the Sun to the GC.

18

The Sun is moving in the clockwise direction, looking from the north Galactic pole perspective.

19

In the case a = 150 au, the particles in e = 0.90 and e = 0.95 orbits were doomed as q<aU holds for them at the time of their initialisation.

20

In fact, none of the particles experienced such large change in a. The discarded particles typically came down to r = aU.

21

In this case, it is meaningful to examine evolution in the Le plane only, as all inclinations were set to zero.

22

Milgrom used q, q2Q2(GM)1/2/(3a03/2)\hbox{$q\equiv-2Q_{2}(GM_{\odot})^{1/2}/(3a^{3/2}_{0})$}, instead of Q2, to quantify the strength of EFE.

23

Starting values of a, e, and i were altered only negligibly during the time of the simulation.

24

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.

25

To be clear in statistical terms and notation, we give a following example: the probability p-value being equal to 0.05 means that the null hypothesis is not rejected at 5% level; in our case the null hypothesis is always the hypothesis in which the two compared datasets have the same distribution.

26

We did not examine how migration through the giant planets’ region, whence P9 would probably originate from, proceeds. The drawback can be the slowness of the EFE induced migration, as P9 in an eccentric orbit crossing orbits of the giant planets, is endangered by their energetic kicks.

27

We compared this with the benchmark model Q2 = 6 simulation.

28

Do not confuse this with the very different tidal effect, arising from differential gravity. Existence of EFE means that the strong equivalence principle is broken.

29

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

  1. Batygin, K., & Brown, M. E. 2016, AJ, 151, 22 [NASA ADS] [CrossRef] [Google Scholar]
  2. Bekenstein, J. D. 2004, Phys. Rev. D, 70, 083509 [NASA ADS] [CrossRef] [Google Scholar]
  3. Bekenstein, J., & Milgrom, M. 1984, ApJ, 286, 7 [NASA ADS] [CrossRef] [Google Scholar]
  4. Blanchet, L., & Novak, J. 2011, MNRAS, 412, 2530 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bosma, A. 1981, AJ, 86, 1825 [NASA ADS] [CrossRef] [Google Scholar]
  6. Boylan-Kolchin, M., Bullock, J. S., & Kaplinghat, M. 2012, MNRAS, 422, 1203 [NASA ADS] [CrossRef] [Google Scholar]
  7. Brasser, R., Duncan, M. J., Levison, H. F., Schwamb, M. E., & Brown, M. E. 2012, Icarus, 217, 1 [NASA ADS] [CrossRef] [Google Scholar]
  8. Brown, M. E. 2001, AJ, 121, 2804 [Google Scholar]
  9. Brown, M. E., & Batygin, K. 2016, ApJ, 824, L23 [NASA ADS] [CrossRef] [Google Scholar]
  10. Brown, R. B., & Firth, J. A. 2016, MNRAS, 456, 1587 [NASA ADS] [CrossRef] [Google Scholar]
  11. Brown, M. E., Trujillo, C., & Rabinowitz, D. 2004, ApJ, 617, 645 [NASA ADS] [CrossRef] [Google Scholar]
  12. Caldwell, N., Walker, M. G., Mateo, M., et al. 2017, ApJ, 839, 20 [NASA ADS] [CrossRef] [Google Scholar]
  13. Chambers, J. E. 1999, MNRAS, 304, 793 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  14. Cowan, N. B., Holder, G., & Kaib, N. A. 2016, ApJ, 822, L2 [NASA ADS] [CrossRef] [Google Scholar]
  15. de la Fuente Marcos, C., & de la Fuente Marcos, R. 2014, MNRAS, 443, L59 [NASA ADS] [CrossRef] [Google Scholar]
  16. de la Fuente Marcos, C., & de la Fuente Marcos, R. 2016, MNRAS, 459, L66 [NASA ADS] [CrossRef] [Google Scholar]
  17. de la Fuente Marcos, C., & de la Fuente Marcos, R., & Aarseth, S. J. 2016, MNRAS, 460, L123 [NASA ADS] [CrossRef] [Google Scholar]
  18. Desmond, H. 2017, MNRAS, 464, 4160 [CrossRef] [Google Scholar]
  19. Dones, L., Brasser, R., Kaib, N., & Rickman, H. 2015, Space Sci. Rev., 197, 191 [NASA ADS] [CrossRef] [Google Scholar]
  20. Einstein, A. 1916, Annalen der Physik, 354, 769 [NASA ADS] [CrossRef] [Google Scholar]
  21. Famaey, B., & McGaugh, S. S. 2012, Liv. Rev. Rel., 15, 10 [Google Scholar]
  22. Fienga, A., Laskar, J., Manche, H., & Gastineau, M. 2016, A&A, 587, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Fortney, J. J., Marley, M. S., Laughlin, G., et al. 2016, ApJ, 824, L25 [NASA ADS] [CrossRef] [Google Scholar]
  24. Gomes, R. S., Matese, J. J., & Lissauer, J. J. 2006, Icarus, 184, 589 [NASA ADS] [CrossRef] [Google Scholar]
  25. Gomes, R. S., Soares, J. S., & Brasser, R. 2015, Icarus, 258, 37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Gulbis, A. A. S., Elliot, J. L., Adams, E. R., et al. 2010, AJ, 140, 350 [NASA ADS] [CrossRef] [Google Scholar]
  27. Haghi, H., Bazkiaei, A. E., Zonoozi, A. H., & Kroupa, P. 2016, MNRAS, 458, 4172 [NASA ADS] [CrossRef] [Google Scholar]
  28. Hees, A., Folkner, W. M., Jacobson, R. A., & Park, R. S. 2014, Phys. Rev. D, 89, 102002 [NASA ADS] [CrossRef] [Google Scholar]
  29. Hees, A., Famaey, B., Angus, G. W., & Gentile, G. 2016, MNRAS, 455, 449 [NASA ADS] [CrossRef] [Google Scholar]
  30. Holman, M. J., & Payne, M. J. 2016a, AJ, 152, 80 [NASA ADS] [CrossRef] [Google Scholar]
  31. Holman, M. J., & Payne, M. J. 2016b, AJ, 152, 94 [NASA ADS] [CrossRef] [Google Scholar]
  32. Iorio, L. 2010, The Open Astron. J., 3, 1 [NASA ADS] [CrossRef] [Google Scholar]
  33. Iorio, L. 2013, Classical and Quantum Gravity, 30, 165018 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  34. Jílková, L., Portegies Zwart, S., Pijloo, T., & Hammer, M. 2015, MNRAS, 453, 3157 [Google Scholar]
  35. Kaib, N. A., & Quinn, T. 2008, Icarus, 197, 221 [NASA ADS] [CrossRef] [Google Scholar]
  36. Kaib, N. A., Roškar, R., & Quinn, T. 2011, Icarus, 215, 491 [NASA ADS] [CrossRef] [Google Scholar]
  37. Lawler, S. M., Shankman, C., Kaib, N., et al. 2017, AJ, 153, 33 [NASA ADS] [CrossRef] [Google Scholar]
  38. Lelli, F., McGaugh, S. S., Schombert, J. M., & Pawlowski, M. S. 2016, ApJ, 827, L19 [NASA ADS] [CrossRef] [Google Scholar]
  39. Lelli, F., McGaugh, S. S., Schombert, J. M., & Pawlowski, M. S. 2017, ApJ, 836, 152 [Google Scholar]
  40. Le Verrier, U. J. 1846, Astron. Nachr., 25, 65 [NASA ADS] [Google Scholar]
  41. Le Verrier, U. J. 1859, Annales de l’Observatoire de Paris, 5 [Google Scholar]
  42. McGaugh, S. S. 2004, ApJ, 609, 652 [NASA ADS] [CrossRef] [Google Scholar]
  43. McGaugh, S. S. 2016, ApJ, 832, L8 [NASA ADS] [CrossRef] [Google Scholar]
  44. McGaugh, S., & Milgrom, M. 2013a, ApJ, 766, 22 [NASA ADS] [CrossRef] [Google Scholar]
  45. McGaugh, S., & Milgrom, M. 2013b, ApJ, 775, 139 [NASA ADS] [CrossRef] [Google Scholar]
  46. McGaugh, S. S., Schombert, J. M., Bothun, G. D., & de Blok, W. J. G. 2000, ApJ, 533, L99 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  47. McGaugh, S. S., Lelli, F., & Schombert, J. M. 2016, Phys. Rev. Lett., 117, 201101 [NASA ADS] [CrossRef] [Google Scholar]
  48. McMillan, P. J., & Binney, J. J. 2010, MNRAS, 402, 934 [NASA ADS] [CrossRef] [Google Scholar]
  49. Milgrom, M. 1983a, ApJ, 270, 371 [Google Scholar]
  50. Milgrom, M. 1983b, ApJ, 270, 365 [NASA ADS] [CrossRef] [Google Scholar]
  51. Milgrom, M. 1994, Annals of Physics, 229, 384 [NASA ADS] [CrossRef] [Google Scholar]
  52. Milgrom, M. 2009a, MNRAS, 399, 474 [NASA ADS] [CrossRef] [Google Scholar]
  53. Milgrom, M. 2009b, ApJ, 698, 1630 [NASA ADS] [CrossRef] [Google Scholar]
  54. Milgrom, M. 2010, MNRAS, 403, 886 [NASA ADS] [CrossRef] [Google Scholar]
  55. Milgrom, M. 2012, MNRAS, 426, 673 [NASA ADS] [CrossRef] [Google Scholar]
  56. Milgrom, M. 2014, MNRAS, 437, 2531 [NASA ADS] [CrossRef] [Google Scholar]
  57. Milgrom, M. 2016a, ArXiv e-prints [arXiv:1605.07458] [Google Scholar]
  58. Milgrom, M. 2016b, Phys. Rev. Lett., 117, 141101 [NASA ADS] [CrossRef] [Google Scholar]
  59. Morbidelli, A., & Levison, H. F. 2004, AJ, 128, 2564 [NASA ADS] [CrossRef] [Google Scholar]
  60. Paučo, R., & Klačka, J. 2016, A&A, 589, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Rubin, V. C., Ford, Jr., W. K., Thonnard, N., & Burstein, D. 1982, ApJ, 261, 439 [NASA ADS] [CrossRef] [Google Scholar]
  62. Sancisi, R. 2004, in Dark Matter in Galaxies, eds. S. Ryder, D. Pisano, M. Walker, & K. Freeman, IAU Symp., 220, 233 [Google Scholar]
  63. Sanders, R. H. 1990, A&ARv, 2, 1 [NASA ADS] [CrossRef] [Google Scholar]
  64. Schönrich, R. 2012, MNRAS, 427, 274 [NASA ADS] [CrossRef] [Google Scholar]
  65. Shankman, C., Kavelaars, J. J., Lawler, S. M., Gladman, B. J., & Bannister, M. T. 2017, AJ, 153, 63 [NASA ADS] [CrossRef] [Google Scholar]
  66. Sheppard, S. S., & Trujillo, C. 2016, AJ, 152, 221 [NASA ADS] [CrossRef] [Google Scholar]
  67. Trujillo, C. A., & Sheppard, S. S. 2014, Nature, 507, 471 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1

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.

Table 2

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

Table 3

Results of the KS test after doubling the initial number of test particles.

Table 4

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

thumbnail 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 well-determined 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 well-determined 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
thumbnail 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 well-determined orbits of TNOS and Centaurs fulfilling the criteria 100 <a< 1500 au, q< 35 au (colour histograms). The distributions inferred from all well-determined orbits of ETNOs (150 <a< 1500 au, q> 30 au) are indicated as dashed bars.

In the text
thumbnail 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 best-fit 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
thumbnail Fig. 4

Evolution of test particles in Le 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 long-living particles in the last Gyr of the simulation. The more transparent tracks refer to evolution in the first 3 Gyr of the simulation and also comprise unstable orbits. Three starting semimajor axes were used as follows: a = 150 au (first column), a = 350 au (second column), and a = 550 au (third column). N-body 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 J2 moment of the Sun. There were no stable orbits at a = 550 au level.

In the text
thumbnail Fig. 5

Evolution of test particles in ωe (first row), Ω−e (second row), and Le (third row) plane as inferred from the benchmark model simulation with Q2 = 1.

In the text
thumbnail Fig. 6

Evolution of test particles in ωe (first row), Ω−e (second row), and Le (third row) plane as inferred from the benchmark model simulation with negative Q2, Q2 = −1.

In the text
thumbnail Fig. 7

ae diagram from the last Gyr of the benchmark model simulation incorporating gravitational interaction with Neptune. Q2 = 1 (top left panel), Q2 = 2 (top right panel), Q2 = 4 (bottom left panel), and Q2 = 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 long-living 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 left-most/upper-most lines (q in range between 100 and 550 au) would hardly be observed if they really exist.

In the text
thumbnail Fig. 8

ae 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 ν2\hbox{$\overline{\nu}_{2}$}. ν2\hbox{$\overline{\nu}_{2}$} yields Q2 = 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
thumbnail Fig. 9

ae 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
thumbnail Fig. 10

ae diagram from the last Gyr of the benchmark model simulation incorporating gravitational interaction with Neptune. Negative values of Q2, Q2 = −1 (left panel) and Q2 = −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 long-living 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
thumbnail 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 full-width 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 long-living particles at t = 4 Gyr. The known ETNOs with q> 35 au (circles) and q< 35 au (triangles) are plotted. TNO 2003 SS422 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
thumbnail Fig. 12

Evolution of test particle orbits, taken from the P9 simulation, in Ω−a (left panel) and La plane (right panel).

In the text
thumbnail Fig. 13

Evolution of test particle orbits, taken from the benchmark model simulation, in ωa plane. The strength of EFE was assumed to be Q2 = 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 full-width horizontal lines.

In the text
thumbnail Fig. 14

Evolution of test particle orbits, taken from the benchmark model simulation, in Ω−a (left panel) and La (right panel) plane. The strength of EFE was assumed to be Q2 = 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 full-width horizontal lines. Increasing Q2, the clustering in L sets in at lower semimajor axes, e.g. with Q2 = 6; similar clustering as in the figure develops for a> 300 au.

In the text
thumbnail 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 Q2 = 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
thumbnail 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 Q2 = 4.

In the text
thumbnail Fig. 17

Perihelion distance as a function of time for 4 selected orbits from the benchmark model simulation characterised by high-q amplitude. Assumed values of Q2 and initial (final) values of a and q are indicated.

In the text
thumbnail Fig. 18

Evolution of test particle orbits, taken from the benchmark model simulation in Ω−a (left panel) and La (right panel) plane. A negative Q2, Q2 = −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 full-width horizontal lines. Decreasing Q2, Q2 < 0, the clustering in L sets in at lower semimajor axes, e.g. with Q2 = −6 similar clustering as in the figure develops for a > 250 au.

In the text
thumbnail 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 Q2, Q2 = −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 (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.