# The imprint of the MOND external field effect

on the Oort cloud comets aphelia distribution

###### Key Words.:

comets: general — gravitation — Oort cloud###### Abstract

Context:Milgromian dynamics (MD or MOND) is a promising physical description excelling especially in galaxies. When formulated as a modified gravity theory, it leads to the so called external field effect (EFE). In the case of the solar system this means that bodies orbiting the Sun are influenced, beyond its tidal effect, by the external gravitational field of the Galaxy with magnitude m s and time-varying direction. Aphelia of intermediate outer Oort cloud (OC) comets (, where is semimajor axis) are distributed non-uniformly on the celestial sphere, showing an apparent concentration around a great circle centered at Galactic longitudes and 135 deg. Such non-uniformity is beyond that attributable to the classical injectors of comets, stellar encounters and the Galactic tides, as well as the expected observational biases.

Aims:We investigated a hypothesis that the great circle concentration of aphelia is a consequence of the long-term action of EFE in the framework of MD.

Methods:We considered exclusively quasi-linear MOND (QUMOND) theory. We built our model of the OC in MD on an analytical approximation of the QUMOND potential for a point mass in the dominant external field of constant magnitude. The model is well applicable at heliocentric distances au. Constraints on the strength of EFE found by the analysis of the Cassini radio-tracking data were taken into account.

Results:We demonstrated characteristic imprint of the EFE on the distribution of aphelia of candidate outer OC comets that migrated down to au. By both analytical and numerical calculations, we showed that the combined effect of EFE and the Galactic tides could qualitatively account for the characteristic features seen in the observed distribution of aphelia of the outer OC comets.

Conclusions:

## 1 Introduction

Spatial orientation of original orbits (i.e. before entering the planetary zone) of new Oort cloud (OC) comets can teach us about the dynamical processes that delivered them from the OC into the inner solar system.
The standard picture is that the Galactic tides and impulses from the stars passing by the OC cooperate in decreasing of perihelion distances of OC bodies, making them eventually observable as new comets (Heisler & Tremaine, 1986; Rickman et al., 2008). The role of passing stars is important especially in the repopulating of the tidal infeed trajectories (Rickman et al., 2008). Aphelia of new outer Oort cloud (OOC) comets (having semimajor axes, , 10 000 au) are distributed non-uniformly on the celestial sphere.
In the Galactic coordinates, the latitudes of aphelia, , are depleted in the equatorial and the polar regions, a feature first noted by (Delsemme, 1987).
This is in accord with what is expected when the Galactic tides make the comets discernible (Torbett, 1986).
When we restrict ourselves to comets with accurately known original orbital energies (class 1A orbits of Marsden & Williams, 2008), distribution^{1}^{1}1For random s, the distribution would be flat. peaks around 0.5 (Matese & Whitmire, 2011). The position of the peak agrees well with the one given by the combined model (stellar encounters + Galactic tides) of Rickman et al. (2008); see their Fig. 5.

Matese et al. (1999) identified an excess of intermediate OOC comets aphelia lying on the celestial sphere close to a great circle centered at Galactic longitudes, , and 135 deg. The great circle excess is still present also in newer data (Matese & Whitmire, 2011). The excess is the most apparent for high-quality orbits, i.e. for the comets about which we are most certain that they are the first-time entrants into the planetary region. The non-uniformity in seems to be beyond that attributable to the known observational biases (Horner & Evans, 2002) and the selection induced by the conventional cometary injection mechanisms (Heisler & Tremaine, 1986; Rickman et al., 2008; Matese et al., 1999). A potential explanation put forward by Matese et al. (1999) is that the excess is induced by a perturbation coming from an unseen massive body with mass of few Jupiter masses, interacting with comets impulsively at a mean distance of au. There is not much room for such perturber in the solar system as implied by the analysis of the Wide-field Infrared Survey Explorer (WISE) data (Luhman, 2014).

Feng & Bailer-Jones (2014) stressed that the aphelia cluster around the preferred plane the normal to which is the solar apex. The authors argue that the non-uniform distribution is an effect of the solar apex motion and that a few strong stellar encounters with proper geometry could be able to produce it. Although, they were not be able to identify such past encounter in their sample of past stellar encounters experienced by the Sun; see also Feng & Bailer-Jones (2015).

Modified Newtonian dynamics (MOND) (Milgrom, 1983; Famaey & McGaugh, 2012 for a review), is a very successful paradigm for galaxies. It predicted existence of the radial acceleration relation (McGaugh et al., 2016), which generalizes well-known systematic properties of galaxies – the mass discrepancy-acceleration relation (Sanders, 1990; McGaugh, 2004) and the baryonic Tully-Fisher relation (McGaugh et al., 2000). In galaxies everything happens as if the gravitation was to be calculated according to Milgrom, not Newton with aid of dark matter (Kroupa et al., 2010; Kroupa, 2012; Famaey & McGaugh, 2012; McGaugh et al., 2016). More than 30 years ago, Milgrom proposed that dynamics in the low acceleration regime – when accelerations are much smaller than some fundamental acceleration , m s – departs from Newtonian dynamics. A particle in the low acceleration regime subject to influence of a spherically symmetric and isolated mass distribution, moves with acceleration of magnitude

(1) |

where is the expected Newtonian gravitational acceleration due to that mass distribution. Eq. (2) is the cornerstone of MOND, meaning that MOND is space-time scale invariant in the low acceleration regime (Milgrom, 2009), contrary to Newtonian dynamics or general relativity. To ensure better clarity, we call fully-fledged classical theories incorporating MOND – Milgromian dynamics (MD).

As first noted by Milgrom (1986), the gravitational potential seen by objects in the OC in MD is not spherically symmetric but dilated along the direction of the external field of the Galaxy. The dilation of the potential is beyond that introduced by the tidal effects. Milgrom speculated that this, plus the enhanced gravitational field of the Sun in MD, could have some impact on our concept of the Oort cloud. Recently, impact of MD on the Oort cloud (Paučo & Klačka, 2016) and the Kuiper belt (Paučo, 2017) was studied. When constraints on the theory based on the Cassini radio-tracking data (Hees et al., 2014, 2016) (the latter hereafter referred to as H+16) are applied, the Oort cloud has a similar size and the planetary barrier operates similarly to its Newtonian counterpart (Paučo & Klačka, 2016). But the delivery of comets could be substantially influenced by the so called external field effect (EFE). EFE could also produce enough torque to detach the orbits of Sedna and 2012 VP (Paučo, 2017). EFE brings in a preferred direction, the time-varying direction of the external field of the Galaxy, which in MD does not decouple from the internal dynamics (Milgrom, 1983; Famaey & McGaugh, 2012). Hence, there are preferred regions in space where EFE causes larger angular momentum changes of OC objects. This in turn could lead to the observed non-uniformity in the distribution of the OOC comets. An imprint of the EFE would be the most easily observed in the regime where tides do not dominate over EFE, i.e. on orbits of intermediate and low- comets. In this paper, we investigate a hypothesis that the great circle concentration (GCC) of aphelia identified by Matese et al. (1999) is caused by the EFE of MD.

## 2 Distribution of outer Oort cloud comets aphelia

The aphelia distribution for the highest-quality (1A) orbits of OOC comets is shown in Fig. 1. This is an equal-area Hammer-Aitoff projection in the Galactic coordinates. The data were taken from Matese & Whitmire (2011), and originally from the 17th Catalogue of Cometary Orbits (Marsden & Williams, 2008). The comets are discerned according to their energies and divided into three groups: i) those with ( au), ii) , iii) ( au). We consider only bounded comets. Comets with aphelion longitudes, , between -180 and -120 deg and latitudes, , between 0 and 30 deg, are probably a residual part of a weak comet shower (Biermann et al., 1983; Matese et al., 1998). This shower is evidenced mainly by the low- (high-) comets.

A strong north-south (NS) bias in the equatorial coordinates of aphelia is prevalent in the cometary catalogue as a whole (Horner & Evans, 2002). The NS bias is a consequence of the simple fact that more observations have been carried out from the northern than from the southern hemisphere in the past. However, when we restrict ourselves to the 1A category orbits, the NS bias is not so apparent. This can be seen in Fig. 1. There are 15 1A comets with aphelion declination higher than 30 deg north, while 19 1A comets have declination higher than 30 deg south. We note that the comets discovered in the north have generally perihelia in the north and vice versa (Horner & Evans, 2002); perihelion in the north means aphelion in the south. Interestingly, there appear to be a preference for the southern aphelia not in the equatorial but in the Galactic coordinates when we look at deg region. If it was not for this NS asymmetry in the Galactic coordinates, the data would look more or less uniformly distributed for deg, with polar holes for deg.

The cometary data are burdened with various selection effects. These include seasonal and diurnal biases, directional and orbital biases, meaning that comets in some directions on the sky or in some specific orbits are more likely to be discovered, as well as some sociological biases (Horner & Evans, 2002). Many of these biases cancel each other, while none of them is able to induce systematics of the observed aphelion distribution (Horner & Evans, 2002).

We show aphelion positions in the Galactic coordinates and the histogram distribution of for 1A-class comets only in Fig. 2. The two comets in the region with longitudes between -180 and -120 deg and latitudes between 0 and 30 deg were discarded as these are probably associated with a comet shower (Biermann et al., 1983; Matese et al., 1998). We discarded also one additional intermediate OOC lying close, in terms of , to the previously discarded couple, and whose aphelion is nearly in the plane of the Galaxy. Therefore, this comet is a good candidate for a comet made discernible due to a stellar encounter. The overpopulated great circle of aphelia identified by Matese et al. (1999) is apparent in the figure, especially in the Galaxy’s equatorial region. These are the puzzling data we aim to explain with the aid of MD.

## 3 Model and analytic calculations

Two classical modified gravity theories, giving arise to MOND phenomenology in galaxies, are in circulation – fully non-linear aquadratic-lagrangian theory (AQUAL; Bekenstein & Milgrom, 1984) and quasi-linear MOND (QUMOND; Milgrom, 2010). To ensure better clarity, we call such fully-fledged classical theories – Milgromian dynamics (MD). AQUAL and QUMOND coincide for spherically symmetric problems (Bekenstein & Milgrom, 1984; Milgrom, 2010), giving the equation of motion in the form of the simple MOND formula (Milgrom, 1983):

(2) |

where is the expected Newtonian gravitational acceleration, is its norm, , , is a modification factor, the so called interpolating function, fulfilling behaviour for (Newtonian limit) and for (deep-MD limit – Eq. (1)), and m s is a new constant of nature. The two theories are distinct outside of the spherical symmetry, although they often give a similar qualitative picture.

We build our model of the OC on the simpler QUMOND theory. In QUMOND, the governing gravitational potential, , is given by (Milgrom, 2010):

(3) |

The problem can be reformulated with aid of phantom matter density (algebraically calculable from the Newtonian potential ) (Milgrom, 2010), hence one does not have to solve Eq. (3) but two linear Poisson equations instead (see e.g. Sect. 6.1.3 of Famaey & McGaugh (2012) for details).

In our exploratory analysis of the problem in MD, we will restrict ourselves to the region of the Oort cloud where the external field from the Galaxy dominates over the internal one from the Sun. This means heliocentric distances larger than , where is gravitational constant, is mass of the Sun, and is magnitude of the external Newtonian gravitational acceleration. Under such condition, one can avoid the direct numerical solution of the Poisson equation. The governing QUMOND potential can be in this case well approximated by (Banik & Zhao, 2015):

(4) |

where , ,

(5) |

and , , , i.e. is the angle between the direction of the external Newtonian field and the direction towards the Sun. is to be calculated from , where can be estimated for instance from the solar motion. We use a simple model of the Sun’s motion – the Sun orbiting in a circular orbit, lying in the Galactic midplane, with angular velocity . Then, is the centripetal acceleration of the Sun in units of .

, as the direction of the external field, coinciding with the direction towards the Galactic center, varies with time. We use

(6) |

The Sun orbits the Galactic center clockwise from the north-Galactic-pole perspective.

For an OC body orbiting in the Galactic midplane , Eqs. (4), (3), and the Gauss planetary equations lead to secular variations:

(7) |

where and are semimajor axis and eccentricity of the body’s orbit, and and are functions involving only eccentricity; details can be found in appendix A. We have averaged over true anomaly. Eqs. (3) form a closed set of equations for , , and . As , the secular change , where is perihelion distance of the orbit, reads

(8) |

As and for , is negative and extremal when or . Moreover, in such case . If we have , where is an integer, decreases the most rapidly on a great circle crossing and 135 deg, same one as the one indicated by the observations of the intermediate OOC comets. However, this does not correspond to the present-day orientation of the external field ( does). In any case, as of the non-stationarity of EFE and combination of the action of EFE and tides, we opt to do some numerical experiments; see Sect. 4.

Eqs. (3) and (8) lead to , i.e. is a decreasing function of . The fact that we observe preference for some values of in the intermediate comets sample, while no preferred values for comets with higher , can be explained by the dominance of the Galactic tide as a comet injector at larger . For angular momentum magnitude changed due to the effect of the Galactic tide, we have (e.g. Torbett, 1986), and we have averaged over true anomaly again. holds approximately for near-parabolic comets, hence we come to .

### 3.1 Tides

Provided that the Galactic mass is distributed axisymmetrically and the Galactic rotation curve is flat at the Sun’s location, the Newtonian tidal acceleration of an OC object orbiting the Sun is (Heisler & Tremaine, 1986; Fouchard, 2004):

(9) |

where , , , and is local mass density in the solar neighborhood. To keep the model in MD simple, we add the tidal acceleration of Eq. (3.1) to the acceleration calculated from Eq. (4). The hint that this might be a good approximation comes from the fact that in QUMOND the problem can be formulated in the Newtonian manner with the aid of phantom mass calculated from the distribution of stars and gas in the Galaxy (Milgrom, 2010; Famaey & McGaugh, 2012). Moreover, at the Sun’s location, an enhancement of the gravity due to MD is low, or equivalently, there is not an obvious need for dark matter at the Sun’s location.

### 3.2 Interpolating function and Galactic parameters

The Cassini radio-tracking data put constraints on the MD interpolating function families as these are related to the strength of EFE which is the dominant effect of MD in the planetary region of the solar system (H+16). Interpolating functions come in pairs with values of , based on the rotation curve fits of galaxies (H+16). As a representative of the allowed interpolating functions we choose

(10) |

family.

In Table 1, we show parameters and of Eq. (4) for a given set – interpolating function, , . Two values of the external field strength are adopted, m s and m s, the same as in H+16. In addition to family, defined in Eq. (3.2), with and 3, we also list and , defined as in H+16, with their boundary s (larger s are allowed by the Cassini data). With m s, the region where the external field dominates over the internal one is au, while with m s, we get au, assuming interpolating function which pairs with m s (H+16). We note that , based on Eq. (2), gives for all allowed interpolating functions.

The values of the angular speed of the Sun used in our models with , which is present in Eqs. (3) and (3.1), are listed in Table 2 for a given external field strength. We determined in two different ways: (i) assuming fixed kpc, we calculated from , and (ii) by making use of the result of McMillan & Binney (2010), i.e. the most probably lies in the interval yr. In (ii), we associate the lowest from this interval with m s and the highest from this interval with m s. The table also lists the corresponding values of the circular speed of the Sun . For the second free parameter in Eq. (3.1), , we assume (Holmberg & Flynn, 2004) throughout the paper.

IF | |||||
---|---|---|---|---|---|

0.815 | 1.9 | 2.3 | 5.04575 | -5.31581 | |

0.815 | 2.4 | 2.9 | 0.22264 | -0.37441 | |

0.743 | 1.9 | 2.6 | 0.00002 | -0.00012 | |

0.743 | 2.4 | 3.2 | 0.00000 | -0.00000 | |

1.450 | 2.4 | 1.7 | 3.36631 | -2.27218 | |

1.460 | 2.4 | 1.6 | 2.80043 | -3.46535 |

m s | ||||
---|---|---|---|---|

1.9 | 2.3 | 8.3 | 221 | 2.723 |

2.4 | 2.9 | 8.3 | 248 | 3.056 |

1.9 | 2.3 | 6.6 | 197 | 3.058 |

2.4 | 2.9 | 7.4 | 234 | 3.232 |

## 4 Simulation

We are interested in the past orbital evolution of OOC comets with intermediate two-body binding energies, ( au). The non-uniformity in the distribution is the most apparent in this -range. This is likely because the tides in this -range are not strong enough to wash away the anisotropy induced by an unknown dynamical mechanism. Also, in the case of MD, EFE is important in this region. The OOC with seem more or less consistent with the trend given by the intermediate comets, though it is only 9 objects; see Fig. 1.

The great circle excess was previously attributed to the impulsive action of a hypothetical massive jovian planet orbiting in the OOC (Matese et al., 1999; Matese & Whitmire, 2011). We hypothesise that the anisotropy, best seen for deg data, comes from the EFE of MD. As a benchmark Milgromian model we use the model with interpolating function, ( m s), yr, and pc, see Tables 1 and 2.

We performed several test-particle simulations in Newtonian and MD, concentrating on the particles’ migration. In Newtonian simulations, gravity of the Sun and tides from the Galaxy (both radial and vertical component) as in Eq. (3.1) are taken into account. In Milgromian simulations, acceleration calculated from the potential in Eq. (4) and tides from the Galaxy as in Eq. (3.1) were simply added together. 180 000 test particles were started with aphelia randomly distributed on the celestial sphere. We started particles with randomly drawn from a uniform distribution between 30 and 60 ( au). Perihelion distances were assigned randomly from a uniform distribution between /60 and au, where is the semimajor axis of a given particle. , , , and (we use the standard notation for the orbital elements, and the angles may be defined with respect to an arbitrary reference frame) were drawn from a flat distribution.

As applicability of Eq. (4) is limited to the region where the gravity of the Sun is smaller than the external field , we had followed the particles until they reached the heliocentric distance au. Then, the particles were discarded from the simulation and their positions, , and velocities, , were recorded. For each recorded particle, the eccentricity vector

(11) |

was calculated. vector was taken as the aphelion direction and the spherical Galactic coordinates of the aphelion were calculated. The orbital integration was carried out in REBOUND (Rein & Liu, 2012). We used WHFast integrator (Rein & Tamayo, 2015) with a timestep of years. The particles were followed for 2 Gyr.

## 5 Simulation results and discussion

In our simulations in MD, we found characteristic imprint of the EFE on the distribution of aphelia of particles recorded at a given time (let us call these injected particles).

In Figs. 3 and 4, the distribution of aphelia in the Galactic coordinates is depicted for the injected particles. When Newtonian, tides only, model is considered (middle row of Fig. 4), the aphelia cluster around deg, without a preference in , as expected (Rickman et al., 2008). As time goes, the tidal infeed trajectories are continually less and less populous and the scatter of aphelia around the values deg continually decreases. This is where the effect of passing stars would step in, repopulating the tidal infeed trajectories (Rickman et al., 2008).

The MD benchmark model simulation shows an apparent non-uniformity not only in the distribution of but as well. In Fig. 4, aphelia of the particles that came down to au during time intervals , , and , where Myr is the period of the Sun’s rotation around the Galactic center in the model, are shown. is set to Myr in order to obtain a meaningfully rich sample. We note that during , the variation of the external field direction is negligible. is the time corresponding to the phase at which the simulated great circle position corresponds to the observed one. The particles’ aphelia cluster around a great circle in . The clustering is the most apparent in the equatorial region. The position of the great circle rotates with the rotating direction of the external field with period . The trend is that after few , the GCC in begins to stand out relatively to the horizontal stripes centered at deg, and the simulated distribution resembles the observed one in Fig. 2; see also Fig. 3. We tested that this picture does not depend on the heliocentric distance where we record the eccentricity vector (let us call it the injection distance); but we have to remain in the region where the gravity of the Sun is smaller than the external field, i.e. au. We repeated the simulations with the injection distance set also to 12 000 and 8000 au, getting qualitatively the same picture.

The simulated GCC matches the position of the observed one at the time when the external field is about 60 deg behind the present-day orientation. The external field needs additional about 38 Myr to align with the present-day orientation. As our model does not allow for tracing of the particles much below the heliocentric distance of 10 000 au, we are not able to determine the orientation of the external field at the time of an actual injection into the inner solar system. In any case, if the delineated position of the GCC is preserved, what is indicated by the fact that the dominant change in perihelion distance is made in the regions where our model applies, then this poses a problem for the EFE scenario. Another complication is that we do not know on what timescale stellar encounters reshuffles the cloud in MD, if we presume that the effect of passing stars is similar as in Newtonian dynamics. Thus, we are not able to make definite quantitative prediction.

The dependence on the interpolating function and the strength of the external field was tested. Restricting to combinations in Table 1, there is no combination but the benchmark triplet , m s, , that would yield a rotating GCC of aphelia of the injected particles. In all the other models, the equatorial region between about deg north and deg south is empty. In the case of and , the aphelia distribution resembles that of the Newtonian, tides-only, model, with some concentrations and dilutions in the horizontal stripes, lying at s whose positions change with time. In all the other remaining combinations, the effect of MD is too weak and is washed away by the Galactic tides. These findings are in concordance with what we find in Paučo (2017) for Sednoids, i.e. we could elucidate both, the orbits of Sednoids and the GCC of intermediate- OOC comets, if we use parameters similar to those in the benchmark model.

We also investigated the dependence on initial semimajor axis. If we initialise particles with lower energies, between 15 and 30 ( au), while leaving the procedure intact otherwise, the injected particles bear only weak imprint of MD. The aphelia distribution resembles that of the Newtonian, tides-only, model, with some barely-noticeable concentrations and dilutions in the horizontal stripes, lying at s whose positions change with time. This is in line with the fact that the GCC, or any other form of -non-uniformity, is not seen in the high- data.

## 6 Summary

Aphelia of intermediate OOC are distributed non-uniformly on the celestial sphere, showing an apparent concentration around the great circle centered at Galactic longitudes and 135 deg. In the framework of MD, we have studied if there exist preferred orientations of orbits of OC bodies with intermediate semimajor axes where their perihelia are more effectively reduced. At heliocentric distances where the gravity of the Sun is much smaller than the external field of the Galaxy, i.e. for au, MD potential incorporating EFE can be well approximated with the analytical formula of Eq. (4). Using Eq. (4), we have shown analytically for orbits in the Galactic equatorial region that such preferred locations, though varying with time, indeed exist. The fact that we do not observe non-uniformity also in the distribution of high- comets can be explained as the consequence of the fact that the influence of EFE weakens as increases, while the influence of the Galactic tides grows stronger rapidly with increasing . We have backed up this claim by numerical experiments as well.

In numerical simulations, we demonstrated characteristic imprint of EFE on the distribution of aphelia of candidate OOC comets that migrated down to au. EFE is able to induce GCC in for the intermediate- bodies. The appearance of the simulated GCC qualitatively agrees with the observed one. The combined effect of EFE and tides leads to the distribution of aphelia resembling the observed one. As the external field rotates around the Sun with period of one Galactic year, the longitude of the GCC of aphelia rotates with the same rate.

We have found only one combination of Cassini-allowed interpolating function and Galactic parameters that leads to sufficiently strong EFE and as a consequence to the GCC in ; see the highlighted line in Table 1. Though, note that the list of considered combinations is far from extensive. These findings are in concordance with the work of Paučo (2017), where we found that similarly strong EFE is needed to explain the large perihelion distances of Sedna and 2012 VP.

As a future prospect, it remains to be demonstrated that the distribution of aphelia of OC comets delivered into the inner solar system can be explained also quantitatively by the combined effect of EFE, Galactic tides, and stellar encounters in the framework of MD. An important question is also whether the GCC induced by EFE matches the position of the observed one at the time of the actual injection into the inner solar system when the external field has the present-day orientation. To accomplish these tasks we would have to switch from the presented semi-analytical modeling to solving the QUMOND equations numerically.

###### Acknowledgements.

I am thankful to Jozef Világi and Leonard Kornoš for providing a computational facility. Simulations in this paper made use of the REBOUND\xspacecode which can be downloaded freely at http://github.com/hannorein/rebound.## References

- Banik & Zhao (2015) Banik, I. & Zhao, H. 2015, ArXiv e-prints [\eprint[arXiv]1509.08457]
- Bekenstein & Milgrom (1984) Bekenstein, J. & Milgrom, M. 1984, ApJ, 286, 7
- Biermann et al. (1983) Biermann, L., Huebner, W. F., & Lust, R. 1983, Proceedings of the National Academy of Science, 80, 5151
- Delsemme (1987) Delsemme, A. H. 1987, A&A, 187, 913
- Famaey & McGaugh (2012) Famaey, B. & McGaugh, S. S. 2012, Living Reviews in Relativity, 15, 10
- Feng & Bailer-Jones (2014) Feng, F. & Bailer-Jones, C. A. L. 2014, MNRAS, 442, 3653
- Feng & Bailer-Jones (2015) Feng, F. & Bailer-Jones, C. A. L. 2015, MNRAS, 454, 3267
- Fouchard (2004) Fouchard, M. 2004, MNRAS, 349, 347
- Hees et al. (2016) Hees, A., Famaey, B., Angus, G. W., & Gentile, G. 2016, MNRAS, 455, 449
- Hees et al. (2014) Hees, A., Folkner, W. M., Jacobson, R. A., & Park, R. S. 2014, Phys. Rev. D, 89, 102002
- Heisler & Tremaine (1986) Heisler, J. & Tremaine, S. 1986, Icarus, 65, 13
- Holmberg & Flynn (2004) Holmberg, J. & Flynn, C. 2004, MNRAS, 352, 440
- Horner & Evans (2002) Horner, J. & Evans, N. W. 2002, MNRAS, 335, 641
- Kroupa (2012) Kroupa, P. 2012, PASA, 29, 395
- Kroupa et al. (2010) Kroupa, P., Famaey, B., de Boer, K. S., et al. 2010, A&A, 523, A32
- Luhman (2014) Luhman, K. L. 2014, ApJ, 781, 4
- Marsden & Williams (2008) Marsden, J. J. & Williams, G. V. 2008, Catalogue of Cometary Orbits, 17th Ed. (Smithsonian Astrophysical Observatory, Cambridge, MA)
- Matese et al. (1998) Matese, J. J., Whitman, P. G., & Whitmire, D. P. 1998, Celestial Mechanics and Dynamical Astronomy, 69, 77
- Matese et al. (1999) Matese, J. J., Whitman, P. G., & Whitmire, D. P. 1999, Icarus, 141, 354
- Matese & Whitmire (2011) Matese, J. J. & Whitmire, D. P. 2011, Icarus, 211, 926
- McGaugh (2004) McGaugh, S. S. 2004, ApJ, 609, 652
- McGaugh et al. (2016) McGaugh, S. S., Lelli, F., & Schombert, J. M. 2016, Physical Review Letters, 117, 201101
- McGaugh et al. (2000) McGaugh, S. S., Schombert, J. M., Bothun, G. D., & de Blok, W. J. G. 2000, ApJ, 533, L99
- McMillan & Binney (2010) McMillan, P. J. & Binney, J. J. 2010, MNRAS, 402, 934
- Milgrom (1983) Milgrom, M. 1983, ApJ, 270, 365
- Milgrom (1986) Milgrom, M. 1986, ApJ, 302, 617
- Milgrom (2009) Milgrom, M. 2009, ApJ, 698, 1630
- Milgrom (2010) Milgrom, M. 2010, MNRAS, 403, 886
- Paučo (2017) Paučo, R. 2017, A&A in press, ArXiv e-prints [\eprint[arXiv]1703.06682]
- Paučo & Klačka (2016) Paučo, R. & Klačka, J. 2016, A&A, 589, A63
- Rein & Liu (2012) Rein, H. & Liu, S.-F. 2012, A&A, 537, A128
- Rein & Tamayo (2015) Rein, H. & Tamayo, D. 2015, MNRAS, 452, 376
- Rickman et al. (2008) Rickman, H., Fouchard, M., Froeschlé, C., & Valsecchi, G. B. 2008, Celestial Mechanics and Dynamical Astronomy, 102, 111
- Sanders (1990) Sanders, R. H. 1990, A&A Rev., 2, 1
- Torbett (1986) Torbett, M. V. 1986, MNRAS, 223, 885

## Appendix A Secular variation of perihelion distance due to EFE

We aim to find the secular evolution of perihelion distance of an OC body, i.e. , driven by EFE in the region where the external fields dominates. Let be components of the perturbing acceleration due to the potential in Eq. (4), in a way that is a component in the direction of the heliocentric position vector , lies in the plane of the body’s orbit and is perpendicular to , and completes the right-handed coordinate system. For simplicity, we consider only bodies orbiting in the Galactic midplane. Using the classical Gauss equations associated with the perturbed two-body problem Sun-OC body, we get the time variation of mutually tied-up orbital elements , , , where is perihelion longitude, as:

(12) |

where

(13) | |||||

For in Eq. (4), we can write on the basis of Eqs. (3) and (A) that

(14) |

Finally, by using

(15) |

we arrive at

(16) |

where

(17) |