We first briefly review results from Newtonian simulations. We used a
reference simulation performed by Nesvorný et al. (2023), which provides
results equivalent to those discussed in Vokrouhlický et al. (2019). Both
references used a very similar methodology to this paper, the planetary
migration model in particular, but they are not identical. For instance,
Vokrouhlický et al. (2019) used the same initial trans-Neptunian disk model as
we use in the MONDian simulations described below, while Nesvorný et al.
(2023) had a broader initial disk of planetesimals extending from Jupiter’s
initial location to $30$ au heliocentric distance. Additionally, several
simulations in Nesvorný et al. (2023) modeled the dynamical effects from the
solar system’s birth cluster. The presence or absence of the birth cluster
does not affect the distribution of binding energies of currently observed
LPCs, at least in the simulations of Nesvorný et al. (2023).

Figure 15: The same as Fig. 13, but now for the simulation ${\cal S}1$ based
on the gradual transition function $\mu_{2}(x)$ and the nominal external
acceleration $g_{\rm e}=2.32\times 10^{-10}$ m s-2. Because of the strength of
the MONDian torques at heliocentric distances $\sim r_{\rm M}\simeq 7,000$ au
(Figs. 5 and 6), the distribution of new comets (class ${\cal C}_{1}$; blue
histogram) extends down to semimajor axes as small as $10^{3}$ au. With more
returns (classes ${\cal C}_{m},m\geq 2$; green and red histograms) the
distribution spreads to even smaller heliocentric distances. The observed Oort
peak of LPCs is not reproduced. This flaw is highlighted for the best-fit
distribution in the right panel. Figure 16: The same as Fig. 14, but now for
the simulation ${\cal S}1$ based on the gradual transition function
$\mu_{2}(x)$ and the nominal external acceleration $g_{\rm e}=2.32\times
10^{-10}$ m s-2. The Oort peak is too broad in the model. In addition, a
population of LPCs on (formally) hyperbolic orbits emerges in this model.

Figure 13 shows the simulated distribution of original semimajor axes $a$ for
LPCs having perihelion $q\leq 4$ au, compared with the same orbits in the
Marsden–Williams catalog. In the left panel, we show the $a$-distribution for
a set of models in which the observable comets are allowed to perform a
maximum number of perihelion passages with $q\leq 4$ au before they disappear
from the sample. This family starts with the class ${\cal C}_{1}$ of comets at
their first visit (“new” comets), continuing with the class ${\cal C}_{2}$ of
comets visiting the target zone at most twice, and so on. In each case we
normalized the simulated distribution to contain the same number of LPCs as
the data. We see from the figure that the classes with too few visits (${\cal
C}_{1}$, ${\cal C}_{2}$, or ${\cal C}_{3}$) over-exaggerate the strength of
the Oort peak ($a\geq 10^{4}$ au), while those for classes with many visible
returns (${\cal C}_{m}$ with $m\geq 15$, say) exaggerate the strength of the
hump of returning LPCs ($a\leq 10^{4}$ au). Clearly a compromise is needed, in
which a combination of results from the ${\cal C}_{m}$ families with a
weighting function depending on $m$ would provide a final model that fits the
data properly. This was known already to Oort (1950), and since then models
with various weighting schemes have been developed. Here we use the simple but
successful model of Whipple (1962), in which the probability that a comet
survives $m$ returns (i.e., $m$ perihelion passages with $q\leq 4$ au) is the
offset power law given by Eq. (38). The best-fit solution is shown by the red
curve on the right panel of Fig. 13 and has $\kappa=0.82$ and $c=1.14$ (this
solution is only slightly better than $\kappa\simeq 0.6$–$0.7$ and $c=0$, the
power-law model that was used, e.g., by Whipple 1962; Wiegert & Tremaine 1999;
Vokrouhlický et al. 2019). The normalized $\chi^{2}$ for this model is $1.5$,
somewhat larger than would be expected from Poisson statistics (a summary of
the best-fit solutions for all models is given in Table 2). Taken literally,
the likelihood of the best-fitting model is rather small, only $0.02$ (see the
$Q$ parameter in Ch. 15.1 of Press et al. 2007). However, the model is not
expected to fit the data perfectly because it does not account for non-
gravitational accelerations due to cometary activity (e.g., outgassing; see
Królikowska 2020), and any fading model based solely on the number of
perihelion passages is oversimplified. Despite these limitations, the model is
able to correctly reproduce the principal characteristics of the distribution
of observed LPC orbits, namely (i) the Oort peak composed of very weakly bound
cometary orbits ($a=10^{4}$–$10^{5}$ au); (ii) the hump of more strongly bound
returning comets ($a<10^{4}$ au), and (iii) the correct relative proportion of
these two classes.

Figure 14 provides a zoomed-in view of the LPC population with semimajor axes
$a>5,000$ au. Following the tradition since Oort (1950), we use the inverse of
the original semimajor axis $1/a$ as the abscissa, instead of $a$ directly.
The gray histogram shows $116$ comets from the Warsaw catalog with $q\leq 4$
au, and the color-coded histograms represent the sequence of ${\cal C}_{m}$
($m\geq 1$) models from our simulation, each normalized to the total number of
comets in the data. Here the new comets (${\cal C}_{1}$ class; blue histogram)
dominate, representing nearly $90$% of the signal in the Oort peak with
$1/a\leq 10^{-4}$ au-1. As a result, the properties of the Oort peak are only
weakly dependent on the uncertain modeling of the fading process. The match to
the data is satisfactory. The main discrepancy is that the simulated Oort peak
is slightly too narrow. Previous studies have shown that non-gravitational
accelerations (which we ignore in our model) are the most likely explanation
of the differences between the model and the data (e.g., Wiegert & Tremaine
1999; Królikowska 2020; Królikowska & Dones 2023). Non-gravitational forces
are likely important for all of the small number of comets on apparently
hyperbolic orbits ($1/a<0$).

Figure 17: The same as Fig. 15, but now for the simulation ${\cal S}2$ based
on the sharp transition function $\mu_{10}(x)$ and the nominal external
acceleration $g_{\rm e}=2.32\times 10^{-10}$ m s-2. The problems in matching
the distribution of the observed LPCs are the same as for simulation ${\cal
S}1$. Figure 18: The same as Fig. 16, but now for the simulation ${\cal S}2$
based on the sharp transition function $\mu_{10}(x)$ and the nominal external
acceleration $g_{\rm e}=2.32\times 10^{-10}$ m s-2. The signature of the Oort
peak is absent in the simulated data.

The Newtonian models we have described provide a benchmark for the success
that can be achieved in matching the data. The agreement with the data is
rather good although not perfect. If MOND is the correct theory of gravity,
then our MONDian simulations of the evolution of LPCs ought to match the data
at least as well and probably better.

Simulations with MOND and the nominal external acceleration (${\cal S}1$ and
${\cal S}2$).– Next, we discuss simulations assuming the nominal value of the
external acceleration, $g_{\rm e}=2.32\times 10^{-10}$ m s-2. We use the same
analysis as we did above for the Newtonian simulation.

Figure 19: The same as Fig. 17, but now for the simulation ${\cal S}3$ based
on the sharp transition function $\mu_{10}(x)$ and a larger external
acceleration, $g_{\rm e}=3.35\times 10^{-10}$ m s-2. Due to the stronger
external acceleration the MONDian perturbation is weaker and the signature of
the Oort peak emerges more clearly. However, even in the class ${\cal C}_{1}$
of new comets there is a large population of LPCs with semimajor axes as small
as $10^{3}$ au. The best-fit model using Whipple’s fading law (38), under-
represents the Oort peak. Figure 20: The same as Fig. 18, but now for the
simulation ${\cal S}3$ based on the sharp transition function $\mu_{10}(x)$
and a larger external acceleration, $g_{\rm e}=3.35\times 10^{-10}$ m s-2.
Here the correspondence between the model and the observations of the
distribution of binding energy $1/a$ for nearly parabolic LPCs appears
satisfactory (i.e., comparable to that in the Newtonian simulation, Fig. 14).
The problem with the simulation, namely the large population of new comets
with semimajor axes of $10^{3}$–$10^{4}$ au (Fig. 19), is beyond the range
displayed here.

Figure 15 compares the distribution of the original semimajor axes $a$ of the
observed LPCs with $q\leq 4$ au to the results of our simulation ${\cal S}1$.
The changes due to MOND are best illustrated by the distribution of new comets
(${\cal C}_{1}$ class; blue histogram). While previously the outer Oort cloud
($a\geq 15,000$ au) was the only source of new comets, the vigorous MONDian
torques at heliocentric distances of a few thousand au (roughly the MOND scale
$\sim r_{\rm M}$) provide an additional source of new comets. The ${\cal
C}_{1}$ class thus is dominated by much smaller semimajor axes. At subsequent
returns (classes ${\cal C}_{m}$, $m\geq 2$), the comets tend to bind more
strongly to the solar system, corresponding to even smaller semimajor axes. In
these classes there is no concentration of comets in the Oort peak between
$10^{4}$ and $10^{5}$ au. A combination of the classes ${\cal C}_{m}$,
weighted by Whipple’s fading law (38), may again be used to find the formal
best fit to the data (right panel on Fig. 15, and Table 2). The best fit, with
parameter values $\kappa=0.67$ and $c=0$, is poor, as the minimum reduced
$\chi^{2}$ is $\simeq 3.6$. The fit is not improved significantly with our
alternative fading law (Table 2).

Figure 16 provides a zoomed-in view of the nearly parabolic LPCs. The
simulation shows only modest excursions from a nearly uniform distribution in
$1/a$, even for the new comets (class ${\cal C}_{1}$), and is inconsistent
with the data because it does not reproduce the width and height of the Oort
peak. Moreover, the simulated distribution now includes far more LPCs on
hyperbolic orbits ($1/a<0$) than the data. In the simulation, these comets are
bound to the solar system, but their orbits have been significantly affected
by MONDian effects. Because the original semimajor axis value is computed by a
Newtonian model (the same way as the data are processed), these objects appear
be approaching the inner solar system on unbound hyperbolic orbits. We have
verified, using smaller scale simulations, that this apparent hyperbolic
population of new LPCs becomes increasingly significant for smaller values of
the external acceleration value $g_{\rm e}$ and for more gradual MONDian
transition functions (see Fig. 8 for context).

The simulation ${\cal S}2$, with a sharp transition function, is shown in
Figs. 17 and 18. This simulation is also inconsistent with the data. The fit
using Whipple’s fading model (38) has parameter values $\kappa=0.85$ and
$c=7.8$, and the normalized $\chi^{2}\simeq 4.8$. The statistical likelihood
of a valid model to have this level of $\chi^{2}$ is now in the land of
despair $3\times 10^{-16}$. The principal mismatch between the model and
observations is in the height and width of the Oort peak. No plausible fading
model can build the Oort peak out of nothing.

The simulation ${\cal S}2$ used $\mu_{10}(x)$ as an example of a sharp
transition function. Using a smaller simulation we verified that the results
are similar for an even sharper transition function, $\mu_{20}(x)$. This is
not surprising since Fig. 1 shows that the difference between these two
transition functions is not very large.

Simulation with MOND and a larger external acceleration (${\cal S}3$).–
Finally, we performed a simulation ${\cal S}3$ using the sharp transition
function $\mu_{10}(x)$ and a 45% larger value of the external acceleration,
$g_{\rm e}=3.35\times 10^{-10}$ m s-2. This test was motivated by current
ideas about the migration of the solar system in the Galaxy: the solar system
probably formed closer to the Galactic center and subsequently evolved outward
to its current location. This would imply that the external acceleration
$g_{\rm e}$ was larger in the past. Ideally, we should conduct simulations in
which $g_{\rm e}$ has a time dependence determined by the Sun’s migration
history. Unfortunately, (i) the solar migration history is still poorly
constrained; (ii) such simulations would require the determination of a
continuous set of radial functions $\upsilon_{\ell}(r,g_{\rm e})$ (Eq. 23).
Given these difficulties, we decided to run a single additional simulation
(${\cal S}_{3}$) in which the external acceleration was constant, and set to
the largest value likely to have been experienced by the solar system during
its history. Recall that with larger values of $g_{\rm e}$ the MONDian effects
are generally reduced, which might help to remove the inconsistencies between
the simulated orbital distribution of the LPCs and the observational data.
Using the maximum possible value thus provides a strong test of whether MOND
is consistent with the orbits of LPCs.

Figures 19 and 20 show the simulated distribution of the original semimajor
axes of LPCs with $q\leq 4$ au, and the $1/a$ distribution for nearly
parabolic comets. While the latter distribution is in relatively good
agreement with the data, the problem with the distribution of semimajor axes
in the whole LPC population persists. Even the new comets (class ${\cal
C}_{1}$; blue histogram) still have a significant maximum below the $10^{4}$
au inner boundary of the Oort peak. For that reason, an attempt to fit the
observations using Whipple’s fading law (Fig. 19, right panel) still under-
represents the Oort peak. The normalized $\chi^{2}=3.6$, while still much
larger than in the Newtonian simulation ($\chi^{2}=1.5$), is lower than in
simulation ${\cal S}_{2}$, confirming that the predictions of MOND become less
inconsistent with the observations when the external field $g_{\rm e}$ is
pushed to higher values.

##  5 Discussion and conclusions

The basic conclusion from this work is that the predictions of the MOND theory
of gravity are in conflict with solar-system observations. However, there are
several layers in this statement, which deserve to be discussed in more
detail. This is because MOND has more degrees of freedom than conventional
theories of gravity (Newtonian gravity and general relativity).

The most severe conflicts with the observational data occur in the versions of
MOND with gradual transition functions. In this paper we have mostly studied a
particular example of a gradual transition function, $\mu_{2}(x)$ (Eq. 4).
MOND theories with this transition function disagree with solar-system data in
at least three ways: (i) They predict planetary ephemerides that disagree with
observations (Blanchet & Novak 2011; Hees et al. 2014; Fienga et al. 2018).
(ii) The torques from the EFE produce relatively broad distributions of
perihelion distances and orbital inclinations of distant TNOs, which are
inconsistent with DES observations. (iii) The distribution of binding energies
of the long-period comets with perihelia $q\leq 4$ au is inconsistent with the
data, because orbits with original semimajor axes $>10,000$ au (the Oort peak)
are under-represented relative to those with smaller semimajor axes. This
statement holds for any fading law that depends only on the number of
perihelion passages close to the Sun.

In spite of its multiple failures, variants of MOND with slow transition
functions are still popular in studies of galaxies because they fit galaxy
rotation curves more accurately (see, e.g., Famaey & McGaugh 2012; Candlish et
al. 2015; Lüghausen et al. 2015; Nagesh et al. 2021; Desmond et al. 2024).
This result may indicate that the transition function is asymmetric, i.e.,
$\mu(x)\to 1$ rapidly when $x\gg 1$, but $\mu(x)\to x$ gradually when $x\ll 1$
(see for example the transition function favored by Desmond et al. 2024). Such
a transition function would not affect our conclusions since we have shown
that MOND is incompatible with the LPC data even if we only consider the
regime $x\lesssim 1$, i.e., $r\gtrsim r_{\rm M}$.

The conflicts with the data are less severe for versions of MOND with sharp
transition functions. We illustrate this class using the transition function
$\mu_{10}(x)$. Adopting the external acceleration $g_{\rm e}=2.32\times
10^{-10}$ m s-2 from Gaia Collaboration et al. (2021) (an independent estimate
from McMillan 2017, is only 7% smaller), this model is consistent with
planetary ephemerides as well as with the currently available survey
constraints on the population of TNOs detached from the scattered disk.
However, the model still predicts a semimajor axis distribution of LPCs that
is inconsistent with the data between a few hundred au and $10^{5}$ au – the
Oort peak at $a>10,000$ au is still too small compared to the population of
comets on more tightly bound orbits (see Figs. 17 and 18).

We find that larger values of the external acceleration $g_{\rm e}$ reduce the
inconsistencies between MOND models with sharp transition functions and the
observational data. Our results imply that in order to be consistent with the
data, $g_{\rm e}$ must be substantially more than 50% larger than its current
value during much of the history of the solar system. It is unlikely that this
constraint can be satisfied in plausible models of the migration history of
the solar system; however, (i) the migration history is not well-understood
and (ii) simulations that would take into account a fully consistent model of
solar-system migration have not yet been done.

Our simulations have other limitations, although we do not believe that any of
these compromise our conclusions:

  * •

We have only tested MOND using a specific one-parameter family of models for
the transition function, $\mu_{n}(x)$ (Eq. 4). However, we believe that this
family is general enough that there are no plausible transition functions that
would be consistent with the observational data on the LPCs.

  * •

We only tested two fading laws (Eqs. 38 and 39), each with two free parameters
that we fit to the data. Both assume that the probability of fading depends
only on the number of perihelion passages close to the Sun. Other fading laws,
which might depend on parameters such as the size of the nucleus rather than
the number of perihelion passages (e.g., Jewitt 2022), might offer better fits
to the data. However, we note that the MONDian simulations of the distribution
of semimajor axes of LPCs already over-predict the number of comets with
$a<10^{4}$ au relative to the number in the Oort peak for comets on their
first perihelion passage (“new” comets) and any fading law that allows comets
to remain visible in subsequent passages makes the disagreement worse.

  * •

Perturbations due to non-gravitational forces have been neglected in our
analysis. These forces arise from outgassing, which typically begins at a
heliocentric distance of a few au, where water ice starts sublimating from the
cometary surface. While the inclusion of the effects of outgassing on the
orbit are important for modeling the details of the binding energy
distribution of LPCs, especially for the Oort-cloud component, their inclusion
is unlikely to cure the differences between the MONDian simulations and the
observational data. The work of Królikowska (2020) suggests the average effect
of non-gravitational forces on the orbits of LPCs is a slight shift toward a
more bound distribution in $1/a$ ($\Delta(1/a)\lesssim 10^{-5}$ au-1), which
is much too small and in the wrong direction to make the results in our ${\cal
S}1$ to ${\cal S}3$ simulations agree with the observed data.

  * •

All of our simulations based on MONDian dynamics ignore the possibility that
tidal fields and stellar encounters were much stronger during the brief period
when the Sun was still part of its birth cluster. However, the effects of the
birth cluster were examined in our Newtonian models (Fig. 9); the presence or
absence of the birth cluster did not affect the distribution of binding
energies of the currently observed LPCs in the Newtonian simulations and we
believe that it is unlikely to do so in MOND.

  * •

We have only examined a single value of the critical acceleration,
$a_{0}=1.2\times 10^{-10}$ m s-2. Fits to galactic rotation curves have led to
$a_{0}$ values that vary from the one we have used by up to $\pm 40\%$ (e.g.,
Hees et al. 2016). However, (i) the best-fit value of $a_{0}$ depends on the
transition function, and for the sharp transition functions that we have
argued are required by solar-system observations the allowed range of $a_{0}$
is smaller; (ii) the characteristic MOND radius $r_{\rm M}$ varies only as the
inverse square root of $a_{0}$, so the range of possible values of $r_{\rm M}$
remains much smaller than the range of semimajor axes in the Oort cloud.

It is not our intention to review the many tests of MOND on the scale of
galaxies and clusters of galaxies (see for example Famaey & McGaugh 2012). We
note, however, that another potential test of Newtonian dynamics and gravity,
on comparable length and mass scales to the LPCs, is provided by wide binary
stars (Chae 2023, 2024; Hernandez 2023; Hernandez et al. 2023). These papers
argue for a breakdown of Newtonian dynamics in binary systems when the
semimajor axis exceeds a few thousand astronomical units, which for these
systems corresponds to a mutual acceleration $\lesssim a_{0}$. The analysis is
challenging, however, because (i) the masses of the stars must be estimated
from spectral fitting; (ii) only two of the three components of the spatial
separation of the two stars can be measured; (iii) only two of the three
relative velocity components can be measured (the referenced papers derive the
kinematics from Gaia measurements of the relative proper motion), and the
observational errors in the velocity are significant; (iv) the binary systems
may contain unseen third components, either bound or unbound (e.g., Clarke
2020). Indeed, other authors argue from similar data that the dynamics of this
class of binaries _is_ well described by Newtonian mechanics (e.g., Pittordis
& Sutherland 2023; Banik et al. 2024; Desmond et al. 2024, but see Hernandez &
Chae 2023). Chae (2023) reports a statistically significant acceleration
excess $g/g_{\rm N}-1$ of $\simeq(12\pm 3)$% at $g_{\rm N}=10^{-8.91}$ m s-2,
or $g_{\rm N}=10.25\,a_{0}$. If (i) this excess is due to MOND, (ii) the AQUAL
formulation is correct, and (iii) the transition function belongs to the
family $\\{\mu_{n}(x)\\}$ (Eq. 4), then the transition function must be
gradual, with $n\lesssim 1$ (see Fig. 1). As summarized in Sec. 2.2.2 such
gradual transition functions are excluded on solar-system scales from
planetary and comet dynamics, although it is easy to construct an _ad hoc_
transition function, slowly varying for $x\lesssim 10$ and rapidly approaching
unity for larger values of the argument, that would be compatible with both
sets of observations.

The heliocentric acceleration of comets in the Oort cloud is significantly
less than the critical MOND acceleration $a_{0}$, and comparable to the
acceleration measured in the outer parts of the rotation curves of galaxies.
However, the Oort cloud differs from galaxies in at least three relevant
aspects: the mass scale is much smaller ($1M_{\odot}$ vs. $10^{11}M_{\odot}$),
the spatial scale ($0.1\,$pc vs. $10^{4}\,$pc), and the strength of the
external field $g_{\rm e}$ ($1.9a_{0}$ vs. a few per cent of $a_{0}$). If MOND
is correct, one or more of these differences may explain why the theory works
on galaxies but fails in the Oort cloud.

Other theories of MOND.– AQUAL is only one of a number of theories of gravity
and dynamics that reproduce the most important features of MOND on galactic
scales: the breakdown of Newtonian dynamics and/or gravity at small
accelerations, the asymptotic flatness of galaxy rotation curves, etc. The
simplest alternative is QUMOND (Milgrom 2010), which we discussed briefly
around Eq. (16), and we believe that simulations based on QUMOND would give
very similar results. More elaborate alternatives include modified-inertia
theories (Milgrom 2011, 2023a), generalized QUMOND (Milgrom 2023b),
tripotential MOND (Milgrom 2023c), a MOND adaptation of Galileon-k-mouflage
theory (Babichev et al. 2011), etc.

Ideally, we should repeat our simulations for these other theories. However,
there are obstacles to doing so:

  * •

In some cases, the theory contains an adjustable screening length $R_{\rm sc}$
that suppresses MONDian effects on small scales. For example, Babichev et al.
(2011) (see also Milgrom 2023b and similar comments in the concluding section
of Banik et al. 2024) proposed a relativistic theory of gravitation of the
tensor-scalar class, for which MOND is the weak-field limit.101010In fact, the
work of Babichev et al. (2011) builds on a more general framework, in which
theorists aim at constructing gravity theories modified at large, usually
cosmological, distances, while preserving the validity of well-tested general
relativity on small spatial scales (e.g., Babichev et al. 2009; Ishak 2019).
The model has negligible deviations from general relativity on scales below
the screening length $R_{\rm sc}$. Although the screening phenomenon is
constructed only for this purpose, the resulting theory is self-consistent and
admissible. Generalizations of QUMOND that depend on higher derivatives of the
potential may also have screening lengths below which the MONDian effects are
strongly suppressed (Milgrom 2023b). In such theories our simulations would
not test MOND, but would only increase the lower limit to $R_{\rm sc}$ from
the size of the planetary system (a few au) to the size of the Oort cloud
(5–10 $\times 10^{4}$ au).

  * •

In some cases, the theory contains enough degrees of freedom that it is
impractical to test with a limited set of numerical simulations. For example,
in tripotential MOND (Milgrom 2023c), the transition function depends on three
variables, not one.

  * •

As discussed already by Milgrom (1983a), some MOND theories preserve Poisson’s
equation and instead modify the law of inertia. Modern variants of this
approach are developed in Milgrom (1994), Milgrom (2011) and Milgrom (2023a).
The more complicated modified-inertia theories bring a time-nonlocal (memory)
effect into the dynamics. For that reason, a thorough analysis of modified-
inertia models in specific astronomical problems has not been carried out so
far; in particular, standard $N$-body techniques could not be used for
simulations of orbit evolution.

The simulations in this paper appear to rule out AQUAL, one of the simplest
and most well-defined MOND theories. All existing MOND theories are heuristic,
and thus will eventually be superseded by a “fundamental” MOND theory (if one
exists). Whatever form this theory takes, it will be strongly constrained by
the requirement that it reproduce the orbital distributions of the long-period
comets and trans-Neptunian objects.

We thank Benoit Famaey and Mordehai Milgrom for thoughtful and constructive
comments. The simulations were performed on the NASA Pleiades Supercomputer.
We thank the NASA NAS computing division for continued support. DN was
supported by the NASA Emerging Worlds program. ST was supported in part by the
Natural Sciences and Engineering Research Council of Canada (NSERC), funding
reference number RGPIN-2020-03885.

## References

  * Adams (2010) Adams, F. C. 2010, ARA&A, 48, 47, doi: 10.1146/annurev-astro-081309-130830
  * Babichev et al. (2011) Babichev, E., Deffayet, C., & Esposito-Farèse, G. 2011, Phys. Rev. D, 84, 061502, doi: 10.1103/PhysRevD.84.061502
  * Babichev et al. (2009) Babichev, E., Deffayet, C., & Ziour, R. 2009, International Journal of Modern Physics D, 18, 2147, doi: 10.1142/S0218271809016107
  * Bailer-Jones (2018) Bailer-Jones, C. A. L. 2018, A&A, 609, A8, doi: 10.1051/0004-6361/201731453
  * Bailer-Jones (2022) —. 2022, ApJ, 935, L9, doi: 10.3847/2041-8213/ac816a
  * Bailer-Jones et al. (2018) Bailer-Jones, C. A. L., Rybizki, J., Andrae, R., & Fouesneau, M. 2018, A&A, 616, A37, doi: 10.1051/0004-6361/201833456
  * Banik et al. (2024) Banik, I., Pittordis, C., Sutherland, W., et al. 2024, MNRAS, 527, 4573, doi: 10.1093/mnras/stad3393
  * Batygin et al. (2019) Batygin, K., Adams, F. C., Brown, M. E., & Becker, J. C. 2019, Phys. Rep., 805, 1, doi: 10.1016/j.physrep.2019.01.009
  * Beaudoin et al. (2023) Beaudoin, M., Gladman, B., Huang, Y., et al. 2023, Planet. Sci. J., 4, 145, doi: 10.3847/PSJ/ace88d
  * Bekenstein & Magueijo (2006) Bekenstein, J., & Magueijo, J. 2006, Phys. Rev. D, 73, 103513, doi: 10.1103/PhysRevD.73.103513
  * Bekenstein & Milgrom (1984) Bekenstein, J., & Milgrom, M. 1984, ApJ, 286, 7, doi: 10.1086/162570
  * Bekenstein (2004) Bekenstein, J. D. 2004, Phys. Rev. D, 70, 083509, doi: 10.1103/PhysRevD.70.083509
  * Bernardinelli et al. (2022) Bernardinelli, P. H., Bernstein, G. M., Sako, M., et al. 2022, ApJS, 258, 41, doi: 10.3847/1538-4365/ac3914
  * Bernstein et al. (2004) Bernstein, G. M., Trilling, D. E., Allen, R. L., et al. 2004, AJ, 128, 1364, doi: 10.1086/422919
  * Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press)
  * Blanchet & Novak (2011) Blanchet, L., & Novak, J. 2011, MNRAS, 412, 2530, doi: 10.1111/j.1365-2966.2010.18076.x
  * Brasser et al. (2006) Brasser, R., Duncan, M. J., & Levison, H. F. 2006, Icarus, 184, 59, doi: 10.1016/j.icarus.2006.04.010
  * Brown & Mathur (2023) Brown, K., & Mathur, H. 2023, AJ, 166, 168, doi: 10.3847/1538-3881/acef1e
  * Brown et al. (2004) Brown, M. E., Trujillo, C., & Rabinowitz, D. 2004, ApJ, 617, 645, doi: 10.1086/422095
  * Candlish et al. (2015) Candlish, G. N., Smith, R., & Fellhauer, M. 2015, MNRAS, 446, 1060, doi: 10.1093/mnras/stu2158
  * Chae (2023) Chae, K.-H. 2023, ApJ, 952, 128, doi: 10.3847/1538-4357/ace101
  * Chae (2024) —. 2024, ApJ, 960, 114, doi: 10.3847/1538-4357/ad0ed5
  * Clarke (2020) Clarke, C. J. 2020, MNRAS, 491, L72, doi: 10.1093/mnrasl/slz161
  * Collins & Sari (2010) Collins, B. F., & Sari, R. 2010, AJ, 140, 1306, doi: 10.1088/0004-6256/140/5/1306
  * Crompvoets et al. (2022) Crompvoets, B. L., Lawler, S. M., Volk, K., et al. 2022, Planet. Sci. J., 3, 113, doi: 10.3847/PSJ/ac67e0
  * Delsemme (1987) Delsemme, A. H. 1987, A&A, 187, 913
  * Desmond et al. (2024) Desmond, H., Hees, A., & Famaey, B. 2024, arXiv e-prints, arXiv:2401.04796, doi: 10.48550/arXiv.2401.04796
  * Dones et al. (2015) Dones, L., Brasser, R., Kaib, N., & Rickman, H. 2015, Space Sci. Rev., 197, 191, doi: 10.1007/s11214-015-0223-2
  * Famaey & McGaugh (2012) Famaey, B., & McGaugh, S. S. 2012, Living Reviews in Relativity, 15, 10, doi: 10.12942/lrr-2012-10
  * Feast & Whitelock (1997) Feast, M., & Whitelock, P. 1997, MNRAS, 291, 683, doi: 10.1093/mnras/291.4.683
  * Feng & Bailer-Jones (2014) Feng, F., & Bailer-Jones, C. A. L. 2014, MNRAS, 442, 3653, doi: 10.1093/mnras/stu1128
  * Feng & Bailer-Jones (2015) —. 2015, MNRAS, 454, 3267, doi: 10.1093/mnras/stv2222
  * Fienga et al. (2018) Fienga, A., Laskar, J., Manche, H., & Gastineau, M. 2018, in Fourteenth Marcel Grossmann Meeting - MG14, ed. M. Bianchi, R. T. Jansen, & R. Ruffini, 3694–3695, doi: 10.1142/9789813226609_0482
  * Fienga & Minazzoli (2024) Fienga, A., & Minazzoli, O. 2024, Living Reviews in Relativity, 27, 1, doi: 10.1007/s41114-023-00047-0
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Klioner, S. A., Mignard, F., et al. 2021, A&A, 649, A9, doi: 10.1051/0004-6361/202039734
  * García-Sánchez et al. (2001) García-Sánchez, J., Weissman, P. R., Preston, R. A., et al. 2001, A&A, 379, 634, doi: 10.1051/0004-6361:20011330
  * Gladman & Chan (2006) Gladman, B., & Chan, C. 2006, ApJ, 643, L135, doi: 10.1086/505214
  * Gladman et al. (2008) Gladman, B., Marsden, B. G., & Vanlaerhoven, C. 2008, in The Solar System Beyond Neptune, ed. M. A. Barucci, H. Boehnhardt, D. P. Cruikshank, A. Morbidelli, & R. Dotson (University of Arizona Press), 43–57
  * Grandclément & Novak (2009) Grandclément, P., & Novak, J. 2009, Living Reviews in Relativity, 12, 1, doi: 10.12942/lrr-2009-1
  * Hees et al. (2016) Hees, A., Famaey, B., Angus, G. W., & Gentile, G. 2016, MNRAS, 455, 449, doi: 10.1093/mnras/stv2330
  * Hees et al. (2014) Hees, A., Folkner, W. M., Jacobson, R. A., & Park, R. S. 2014, Phys. Rev. D, 89, 102002, doi: 10.1103/PhysRevD.89.102002
  * Heisler & Tremaine (1986) Heisler, J., & Tremaine, S. 1986, Icarus, 65, 13, doi: 10.1016/0019-1035(86)90060-6
  * Hernandez (2023) Hernandez, X. 2023, MNRAS, 525, 1401, doi: 10.1093/mnras/stad2306
  * Hernandez & Chae (2023) Hernandez, X., & Chae, K.-H. 2023, arXiv e-prints, arXiv:2312.03162, doi: 10.48550/arXiv.2312.03162
  * Hernandez et al. (2023) Hernandez, X., Verteletskyi, V., Nasser, L., & Aguayo-Ortiz, A. 2023, MNRAS, in press, doi: 10.1093/mnras/stad3446
  * Ishak (2019) Ishak, M. 2019, Living Reviews in Relativity, 22, 1, doi: 10.1007/s41114-018-0017-4
  * Jewitt (2022) Jewitt, D. 2022, AJ, 164, 158, doi: 10.3847/1538-3881/ac886d
  * Kaib (2022) Kaib, N. A. 2022, Science Advances, 8, eabm9130, doi: 10.1126/sciadv.abm9130
  * Kaib & Quinn (2008) Kaib, N. A., & Quinn, T. 2008, Icarus, 197, 221, doi: 10.1016/j.icarus.2008.03.020
  * Kaib et al. (2011) Kaib, N. A., Roškar, R., & Quinn, T. 2011, Icarus, 215, 491, doi: 10.1016/j.icarus.2011.07.037
  * Kaib & Sheppard (2016) Kaib, N. A., & Sheppard, S. S. 2016, AJ, 152, 133, doi: 10.3847/0004-6256/152/5/133
  * Królikowska (2014) Królikowska, M. 2014, A&A, 567, A126, doi: 10.1051/0004-6361/201323263
  * Królikowska (2020) —. 2020, A&A, 633, A80, doi: 10.1051/0004-6361/201936316
  * Królikowska & Dones (2023) Królikowska, M., & Dones, L. 2023, A&A, 678, A113, doi: 10.1051/0004-6361/202347178
  * Królikowska & Dybczyński (2017) Królikowska, M., & Dybczyński, P. A. 2017, MNRAS, 472, 4634, doi: 10.1093/mnras/stx2157
  * Królikowska & Dybczyński (2020) —. 2020, A&A, 640, A97, doi: 10.1051/0004-6361/202038451
  * Królikowska et al. (2014) Królikowska, M., Sitarski, G., Pittich, E. M., et al. 2014, A&A, 571, A63, doi: 10.1051/0004-6361/201424329
  * Levison & Duncan (1994) Levison, H. F., & Duncan, M. J. 1994, Icarus, 108, 18, doi: 10.1006/icar.1994.1039
  * Lüghausen et al. (2015) Lüghausen, F., Famaey, B., & Kroupa, P. 2015, Canadian Journal of Physics, 93, 232, doi: 10.1139/cjp-2014-0168
  * Maquet & Pierret (2015) Maquet, L., & Pierret, F. 2015, Phys. Rev. D, 91, 084015, doi: 10.1103/PhysRevD.91.084015
  * Marsden & Williams (2008) Marsden, B. G., & Williams, G. V. 2008, Catalogue of Cometary Orbits 2008, 17th Edition (IAU Minor Planet Center/Central Bureau for Astronomical Telegrams, Cambridge, MA)
  * Matese et al. (1999) Matese, J. J., Whitman, P. G., & Whitmire, D. P. 1999, Icarus, 141, 354, doi: 10.1006/icar.1999.6177
  * Matese & Whitmire (2011) Matese, J. J., & Whitmire, D. P. 2011, Icarus, 211, 926, doi: 10.1016/j.icarus.2010.11.009
  * McMillan (2017) McMillan, P. J. 2017, MNRAS, 465, 76, doi: 10.1093/mnras/stw2759
  * Migaszewski (2023) Migaszewski, C. 2023, MNRAS, 525, 805, doi: 10.1093/mnras/stad2250
  * Milgrom (1983a) Milgrom, M. 1983a, ApJ, 270, 365, doi: 10.1086/161130
  * Milgrom (1983b) —. 1983b, ApJ, 270, 371, doi: 10.1086/161131
  * Milgrom (1983c) —. 1983c, ApJ, 270, 384, doi: 10.1086/161132
  * Milgrom (1986a) —. 1986a, ApJ, 302, 617, doi: 10.1086/164021
  * Milgrom (1986b) —. 1986b, ApJ, 306, 9, doi: 10.1086/164314
  * Milgrom (1994) —. 1994, Annals of Physics, 229, 384, doi: 10.1006/aphy.1994.1012
  * Milgrom (2009) —. 2009, MNRAS, 399, 474, doi: 10.1111/j.1365-2966.2009.15302.x
  * Milgrom (2010) —. 2010, MNRAS, 403, 886, doi: 10.1111/j.1365-2966.2009.16184.x
  * Milgrom (2011) —. 2011, Acta Physica Polonica B, 42, 2175, doi: 10.5506/APhysPolB.42.2175
  * Milgrom (2012) —. 2012, MNRAS, 426, 673, doi: 10.1111/j.1365-2966.2012.21719.x
  * Milgrom (2015) —. 2015, Canadian Journal of Physics, 93, 107, doi: 10.1139/cjp-2014-0211
  * Milgrom (2020) —. 2020, Studies in the History and Philosophy of Modern Physics, 71, 170, doi: 10.1016/j.shpsb.2020.02.004
  * Milgrom (2023a) —. 2023a, arXiv e-prints, arXiv:2310.14334, doi: 10.48550/arXiv.2310.14334
  * Milgrom (2023b) —. 2023b, Phys. Rev. D, 108, 084005, doi: 10.1103/PhysRevD.108.084005
  * Milgrom (2023c) —. 2023c, Phys. Rev. D, 108, 063009, doi: 10.1103/PhysRevD.108.063009
  * Morbidelli & Levison (2004) Morbidelli, A., & Levison, H. F. 2004, AJ, 128, 2564, doi: 10.1086/424617
  * Nagesh et al. (2021) Nagesh, S. T., Banik, I., Thies, I., et al. 2021, Canadian Journal of Physics, 99, 607, doi: 10.1139/cjp-2020-0624
  * Nesvorný (2015) Nesvorný, D. 2015, AJ, 150, 68, doi: 10.1088/0004-6256/150/3/68
  * Nesvorný (2018) —. 2018, ARA&A, 56, 137, doi: 10.1146/annurev-astro-081817-052028
  * Nesvorný et al. (2023) Nesvorný, D., Bernardinelli, P., Vokrouhlický, D., & Batygin, K. 2023, Icarus, 406, 115738, doi: 10.1016/j.icarus.2023.115738
  * Nesvorný & Morbidelli (2012) Nesvorný, D., & Morbidelli, A. 2012, AJ, 144, 117, doi: 10.1088/0004-6256/144/4/117
  * Nesvorný & Vokrouhlický (2016) Nesvorný, D., & Vokrouhlický, D. 2016, ApJ, 825, 94, doi: 10.3847/0004-637X/825/2/94
  * Nesvorný et al. (2017) Nesvorný, D., Vokrouhlický, D., Dones, L., et al. 2017, ApJ, 845, 27, doi: 10.3847/1538-4357/aa7cf6
  * Nesvorný et al. (2016) Nesvorný, D., Vokrouhlický, D., & Roig, F. 2016, ApJ, 827, L35, doi: 10.3847/2041-8205/827/2/L35
  * Nesvorný et al. (2020) Nesvorný, D., Vokrouhlický, D., Alexandersen, M., et al. 2020, AJ, 160, 46, doi: 10.3847/1538-3881/ab98fb
  * Oort (1950) Oort, J. H. 1950, Bull. Astron. Inst. Netherlands, 11, 91
  * Paučo (2017) Paučo, R. 2017, A&A, 603, A11, doi: 10.1051/0004-6361/201630335
  * Paučo & Klačka (2016) Paučo, R., & Klačka, J. 2016, A&A, 589, A63, doi: 10.1051/0004-6361/201527713
  * Paučo & Klačka (2017) —. 2017, arXiv e-prints, arXiv:1705.09273, doi: 10.48550/arXiv.1705.09273
  * Penner (2020a) Penner, A. R. 2020a, Ap&SS, 365, 154, doi: 10.1007/s10509-020-03870-x
  * Penner (2020b) —. 2020b, Ap&SS, 365, 65, doi: 10.1007/s10509-020-03776-8
  * Pittordis & Sutherland (2023) Pittordis, C., & Sutherland, W. 2023, The Open Journal of Astrophysics, 6, 4, doi: 10.21105/astro.2205.02846
  * Press et al. (2007) Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical Recipes: The Art of Scientific Computing (Cambridge University Press, Cambridge)
  * Rickman (2010) Rickman, H. 2010, in Lecture Notes in Physics, Berlin Springer Verlag, Vol. 790, Lecture Notes in Physics, Berlin Springer Verlag, ed. J. Souchay & R. Dvorak, 341–399, doi: 10.1007/978-3-642-04458-8_7
  * Rickman et al. (2008) Rickman, H., Fouchard, M., Froeschlé, C., & Valsecchi, G. B. 2008, Celestial Mechanics and Dynamical Astronomy, 102, 111, doi: 10.1007/s10569-008-9140-y
  * Sanders (2015) Sanders, R. H. 2015, Canadian Journal of Physics, 93, 126, doi: 10.1139/cjp-2014-0206
  * Silsbee & Tremaine (2018) Silsbee, K., & Tremaine, S. 2018, AJ, 155, 75, doi: 10.3847/1538-3881/aaa19b
  * Skordis (2009) Skordis, C. 2009, Classical and Quantum Gravity, 26, 143001, doi: 10.1088/0264-9381/26/14/143001
  * Vokrouhlický et al. (2019) Vokrouhlický, D., Nesvorný, D., & Dones, L. 2019, AJ, 157, 181, doi: 10.3847/1538-3881/ab13aa
  * Weissman (1979) Weissman, P. R. 1979, in Dynamics of the Solar System, ed. R. L. Duncombe, Vol. 81, 277
  * Weissman (1980) Weissman, P. R. 1980, A&A, 85, 191
  * Whipple (1962) Whipple, F. L. 1962, AJ, 67, 1, doi: 10.1086/108596
  * Wiegert & Tremaine (1999) Wiegert, P., & Tremaine, S. 1999, Icarus, 137, 84, doi: 10.1006/icar.1998.6040
  * Will (2014) Will, C. M. 2014, Living Reviews in Relativity, 17, 4, doi: 10.12942/lrr-2014-4
  * Wisdom & Holman (1991) Wisdom, J., & Holman, M. 1991, AJ, 102, 1528, doi: 10.1086/115978

##  Appendix A Testing the accuracy of the numerical solution for the MOND
potential with analytical results

The radial functions $\upsilon_{\ell}(r)$ of the MOND potential $u({\bf r})$
(23) are determined numerically in terms of Chebyshev polynomials, segmented
in a number of radial zones between the origin and spatial infinity (see Sec.
2.2.1 and Appendix B). We seek tests of the accuracy of our solutions using
comparisons to previously published numerical results and analytic asymptotic
behaviors at $r\rightarrow\infty$ (tests of the solution near the origin $r=0$
have already been discussed in Sec. 2.2.2).

Testing the far zone: $\propto 1/r$ asymptotic behavior.– At spatial infinity,
Blanchet & Novak (2011) re-derive in their Appendix A an earlier result of
Bekenstein & Milgrom (1984), namely an asymptotic behavior of $u$

$u\left(r,\theta\right)=\frac{GM}{r\mu_{e}\sqrt{1+\lambda_{e}\sin^{2}\theta}}+O\left(\frac{1}{r^{2}}\right)\;,$
(A1)

where $\mu_{e}=\mu(x_{e})$ and $\lambda_{e}=[d\log\mu(x)/d\log x]_{x_{\rm
e}}$, with $x_{e}=g_{\rm e}/a_{0}$. To connect this asymptotic approximation
of $u$ with its multipole series representation (23), we rewrite this in terms
of Legendre polynomials (recall $\tau=\cos\theta$)

$u\left(r,\tau\right)=\frac{GM}{r}\sum_{\ell\geq
0}w_{2\ell}\,P_{2\ell}(\tau)+O\left(\frac{1}{r^{2}}\right)\;.$ (A2)

Straightforward algebra yields

$w_{2\ell}=\frac{\kappa_{2\ell}}{\mu_{e}\alpha_{e}^{\ell-1}\sqrt{\alpha_{e}\lambda_{e}}}\left[\mathbb{A}_{2\ell}\left(\alpha_{e}\right)\frac{{\rm
arcsin}\sqrt{\alpha_{e}}}{\sqrt{\alpha_{e}}}-\mathbb{B}_{2\ell}\left(\alpha_{e}\right)\sqrt{1-\alpha_{e}}\right]\;,$
(A3)

where $\alpha_{e}=\lambda_{e}/(1+\lambda_{e})$, $\kappa_{2\ell}$ are numerical
coefficients, and $\mathbb{A}_{2\ell}\left(\alpha_{e}\right)$ and
$\mathbb{B}_{2\ell}\left(\alpha_{e}\right)$ are polynomials of order $\ell$
and $\ell-1$. For the first few terms111111Milgrom (2010), Eq. (65), derives
within the QUMOND variant of MOND a result equivalent to our monopole term
given in (A4) with (A3). In particular, he shows that the effective dynamical
mass of the system is equal to $M\,w_{0}$ (in our notation), whose excess over
$M$ arises due to the MONDian phantom mass. However, this interpretation must
be taken with some caution, since the asymptotic acceleration is not
spherically symmetric (see Eq. A1). we obtained $\kappa_{0}=1$,
$\kappa_{2}=15/4$, $\kappa_{4}=945/64$, $\kappa_{6}=15015/256$,
$\kappa_{8}=3828825/16384$, $\kappa_{10}=61108047/65536$, and the $\mathbb{A}$
and $\mathbb{B}$ polynomials read

$\displaystyle\mathbb{A}_{0}\left(\alpha_{e}\right)$ $\displaystyle=$
$\displaystyle 1\;,\;\;\,\mathbb{B}_{0}\left(\alpha_{e}\right)=0\;,$ (A4)
$\displaystyle\mathbb{A}_{2}\left(\alpha_{e}\right)$ $\displaystyle=$
$\displaystyle
1-\frac{2\alpha_{e}}{3}\;,\;\;\,\mathbb{B}_{2}\left(\alpha_{e}\right)=1\;,$
(A5) $\displaystyle\mathbb{A}_{4}\left(\alpha_{e}\right)$ $\displaystyle=$
$\displaystyle
1-\frac{8\alpha_{e}}{7}+\frac{8\alpha_{e}^{2}}{35}\;,\;\;\,\mathbb{B}_{4}\left(\alpha_{e}\right)=1-\frac{10\alpha_{e}}{21}\;,$
(A6) $\displaystyle\mathbb{A}_{6}\left(\alpha_{e}\right)$ $\displaystyle=$
$\displaystyle
1-\frac{126\alpha_{e}}{77}+\frac{8\alpha_{e}^{2}}{11}-\frac{16\alpha_{e}^{3}}{231}\;,\;\;\,\mathbb{B}_{6}\left(\alpha_{e}\right)=1-\frac{32\alpha_{e}}{33}+\frac{28\alpha_{e}^{2}}{165}\;,$
(A7) $\displaystyle\mathbb{A}_{8}\left(\alpha_{e}\right)$ $\displaystyle=$
$\displaystyle
1-\frac{32\alpha_{e}}{15}+\frac{96\alpha_{e}^{2}}{65}-\frac{256\alpha_{e}^{3}}{715}+\frac{128\alpha_{e}^{4}}{6435}\;,\;\;\,\mathbb{B}_{8}\left(\alpha_{e}\right)=1-\frac{22\alpha_{e}}{15}+\frac{344\alpha_{e}^{2}}{585}-\frac{12176\alpha_{e}^{3}}{225225}\;,$
(A8) $\displaystyle\mathbb{A}_{10}\left(\alpha_{e}\right)$ $\displaystyle=$
$\displaystyle
1-\frac{50\alpha_{e}}{19}+\frac{800\alpha_{e}^{2}}{323}-\frac{320\alpha_{e}^{3}}{323}+\frac{640\alpha_{e}^{4}}{4199}-\frac{256\alpha_{e}^{5}}{46189}\;,\;\;\,\mathbb{B}_{10}\left(\alpha_{e}\right)=1-\frac{112\alpha_{e}}{57}+\frac{2028\alpha_{e}^{2}}{1615}-\frac{3232\alpha_{e}^{3}}{11305}+\frac{21472\alpha_{e}^{4}}{1322685}\;.$

In the limit $x_{e}\rightarrow 0$ the exterior field becomes negligible, and
we obtain $\lambda_{e}\rightarrow 1$, $\alpha_{e}\rightarrow 1/2$ and
$\mu_{e}\rightarrow 0$. In this case, however, the limit is discontinuous, as
the assumption $u(r)\propto 1/r$ from (A1) does not hold – recall from Sec. 2
that the MONDian potential of an isolated spherical system has a logarithmic
divergence at infinity (Eq. 11). In the opposite limit
$x_{e}\rightarrow\infty$, we have $\lambda_{e}\rightarrow 0$,
$\alpha_{e}\rightarrow 0$ and $\mu_{e}\rightarrow 1$ – in this situation the
MONDian effects are removed, $w_{0}\rightarrow 1$, $w_{2\ell}\rightarrow 0$
for $\ell\geq 1$; and $u(r,\tau)=GM/r$, the Newtonian result.

We used both the complete asymptotic expression (A1), and the coefficients of
its lowest multipoles in (A2) computed using (A3) to (LABEL:w10asym), to
verify the accuracy of our numerical evaluation of $u(r,\theta)$. First, we
computed the potential $u$ at the outermost radial domain $r_{\rm
max}^{I}=R^{\prime}$, thus $u(R^{\prime},\theta)$, and we compared it with
(A1). The values matched with fractional accuracy of $\simeq
10^{-4}$–$10^{-5}$ for all our solutions. Second, we used the behavior of the
radial functions $\upsilon_{\ell}(r)$ at the outermost radial zone to
determine the coefficients $w_{2\ell}$ numerically. Their correspondence to
the analytic form in (A3) was again very accurate (see Fig. 8).

Testing the far zone: $\propto 1/r^{2}$ asymptotic behavior.– As indicated in
Eq. (A2), the above analysis permits us to test our numerical evaluation of
the even-degree multipoles of the MOND potential $u(r,\theta)$, which
asymptotically vary $\propto 1/r$. The odd-degree multipoles decay faster
($\propto 1/r^{2}$). We would like to check their asymptotic behavior, in
particular to verify the dipole and octupole terms which are important in the
formulation of the MONDian potential. We generalized Eq. (A1) using the
expression

$u\left(r,\theta\right)=\frac{GM}{r\mu_{e}\sqrt{1+\lambda_{e}\sin^{2}\theta}}+\frac{K_{2}}{r^{2}}\frac{g_{2}\left(\theta\right)}{1+\lambda_{e}\sin^{2}\theta}+O\left(\frac{1}{r^{3}}\right)\;,$
(A10)

where $g_{2}(\theta)$ is an unknown function at this moment. Following the
approach in the Appendix of Blanchet & Novak (2011), we seek $g_{2}(\theta)$
such that ${\bf g}=\nabla U$ satisfies $\nabla\cdot(\mu{\bf g})=0$ far from
mass sources (Eq. 12)-. Switching to a new independent variable
$\tau=\cos\theta$, introducing $G_{2}(\tau)=g_{2}(\theta)$, and collecting
terms straightforward algebra results in a linear, second-order differential
equation in the Sturm–Liouville form

$\frac{d}{d\tau}\left[\left(1-\tau^{2}\right)\phi\,\frac{dG_{2}}{d\tau}\right]+\frac{2\alpha_{e}}{\lambda_{e}}\frac{G_{2}}{\phi^{3}}=0\;,$
(A11)

with $\phi=\phi(\tau)=\sqrt{1-\alpha_{e}\tau^{2}}$. Further simplification
follows from a transformation
$\tau\rightarrow\zeta=\sqrt{1-\alpha_{e}}\,\tau/\phi$, which brings Eq. (A11)
to a Legendre differential equation for eigenvalue equal to $2$. Choosing a
dipole solution which is bound on the interval $[-1,1]$, explicitly
$G_{2}(\zeta)=\zeta$, we thus finally obtain the remarkably simple form

$u\left(r,\theta\right)=\frac{GM}{r\mu_{e}\sqrt{1+\lambda_{e}\sin^{2}\theta}}\left(1+\frac{K^{\prime}_{2}}{r}\frac{\cos\theta}{1+\lambda_{e}\sin^{2}\theta}\right)+O\left(\frac{1}{r^{3}}\right)\;.$
(A12)

All terms asymptotically behaving as $\propto 1/r^{2}$ are antisymmetric under
the $\tau\rightarrow-\tau$ transformation. They become the leading-order odd-
degree terms in Legendre series representation of $u(r,\theta)$, namely

$u\left(r,\tau\right)=\frac{GM}{r}\sum_{\ell\geq
0}w_{2\ell}\,P_{2\ell}(\tau)+\frac{K_{2}}{r^{2}}\sum_{\ell\geq
0}w_{2\ell+1}\,P_{2\ell+1}(\tau)+O\left(\frac{1}{r^{3}}\right)\;,$ (A13)

with

$w_{2\ell+1}=\frac{\kappa_{2\ell+1}}{\alpha_{e}^{\ell-2}\left(\alpha_{e}\lambda_{e}\right)^{3/2}}\left[\mathbb{C}_{2\ell+1}\left(\alpha_{e}\right)\frac{1}{\sqrt{1-\alpha_{e}}}-\mathbb{D}_{2\ell+1}\left(\alpha_{e}\right)\frac{{\rm
arcsin}\sqrt{\alpha_{e}}}{\sqrt{\alpha_{e}}}\right]\;.$ (A14)

Here again $\kappa_{2\ell+1}$ are numerical coefficients,
while$\mathbb{C}_{2\ell+1}\left(\alpha_{e}\right)$ and
$\mathbb{D}_{2\ell+1}\left(\alpha_{e}\right)$ are polynomials of the order
$\ell$. For the lowest terms we obtained $\kappa_{1}=3$, $\kappa_{3}=105/4$,
$\kappa_{5}=10395/64$, $\kappa_{7}=225225/256$, and the $\mathbb{C}$ and
$\mathbb{D}$ polynomials read

$\displaystyle\mathbb{C}_{1}\left(\alpha_{e}\right)$ $\displaystyle=$
$\displaystyle 1\;,\;\;\,\mathbb{D}_{1}\left(\alpha_{e}\right)=1\;,$ (A15)
$\displaystyle\mathbb{C}_{3}\left(\alpha_{e}\right)$ $\displaystyle=$
$\displaystyle
1-\frac{11\alpha_{e}}{15}\;,\;\;\,\mathbb{D}_{3}\left(\alpha_{e}\right)=1-\frac{2\alpha_{e}}{5}\;,$
(A16) $\displaystyle\mathbb{C}_{5}\left(\alpha_{e}\right)$ $\displaystyle=$
$\displaystyle
1-\frac{11\alpha_{e}}{9}+\frac{274\alpha_{e}^{2}}{945}\;,\;\;\,\mathbb{D}_{5}\left(\alpha_{e}\right)=1-\frac{8\alpha_{e}}{9}+\frac{8\alpha_{e}^{2}}{63}\;,$
(A17) $\displaystyle\mathbb{C}_{7}\left(\alpha_{e}\right)$ $\displaystyle=$
$\displaystyle
1-\frac{67\alpha_{e}}{39}+\frac{1784\alpha_{e}^{2}}{2145}-\frac{484\alpha_{e}^{3}}{5005}\;,\;\;\,\mathbb{D}_{7}\left(\alpha_{e}\right)=1-\frac{18\alpha_{e}}{13}+\frac{72\alpha_{e}^{2}}{143}-\frac{16\alpha_{e}^{3}}{429}\;.$
(A18)

We did not find a constraint on the scaling constant $K_{2}$ in (A13), because
Eq. (A23) is at the $1/r^{3}$ level in $\mu{\bf g}$ automatically satisfied by
virtue of the antisymmetry of the odd-degree multipoles. As a result, our
procedure can only check the relative (rather than absolute) values of the
multipole coefficients $w_{2\ell+1}$. Nevertheless, this still provides a
stringent test of our numerical solution for the potential $u(r,\theta)$. For
instance, we easily find that $w_{1}:w_{3}:w_{5}:w_{7}\simeq
1:8.73:115.03:1754.99$ for the transition function $\mu_{2}(x)$ and the
nominal external acceleration $g_{\rm e}=2.32\times 10^{-10}$ m s-2, and we
verified these ratios hold exactly in the full-scale solution of $u(r,\theta)$
shown in Fig. 4.

Testing the far zone: $\propto 1/r^{3}$ asymptotic behavior.– The procedure
outlined above may be pushed to further orders in the power $1/r$, but it
becomes more algebraically complicated. Here we briefly comment on the
$\propto 1/r^{3}$ asymptotics. Introducing the corresponding term in the non-
Newtonian potential

$\delta
u\left(r,\theta\right)=\frac{K_{3}}{r^{3}}\frac{g_{3}\left(\theta\right)}{(1+\lambda_{e}\sin^{2}\theta)^{3/2}}\;,$
(A19)

we seek again $g_{3}(\theta)$ to satisfy the source-free field equation
$\nabla\cdot(\mu{\bf g})=0$. Choosing $\tau=\cos\theta$ as the independent
variable, and $G_{3}(\tau)=g_{3}(\theta)$, we find

$\frac{d}{d\tau}\left[\left(1-\tau^{2}\right)\phi\,\frac{dG_{3}}{d\tau}\right]+\frac{6\alpha_{e}}{\lambda_{e}}\frac{G_{3}}{\phi^{3}}=S_{3}\left(\tau\right)\;,$
(A20)

where again $\phi=\sqrt{1-\alpha_{e}\tau^{2}}$. The important novelty is the
source term $S_{3}(\tau)$ (compare with Eq. A11), which is a non-linear
function of the potential asymptotics behaving as $\propto 1/r$ from Eq. (A1),
known at this moment. We do not develop it in full detail here, just note that
one can easily show that $S_{3}(\tau)$ is an antisymmetric function and may be
developed in odd-degree Legendre polynomials. The inhomogeneity of (A20)
implies that a general solution is composed of a free part, from the
homogeneous reduction, and a forced part, specific to the source on the right-
hand side, therefore $G_{3}=G_{3}^{\rm free}+G_{3}^{\rm forced}$. Upon a final
transformation to $\zeta(\tau)=\sqrt{1-\alpha_{e}}\,\tau/\phi$, as above, one
easily finds that (A20) transforms to an inhomogeneous Legendre equation with
eigenvalue $6$, unsurprisingly corresponding now to the quadrupole solution.
As a result

$\displaystyle G_{3}^{\rm free}\left(\tau\right)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(3\zeta^{2}-1\right)\;,$ (A21) $\displaystyle
G_{3}^{\rm forced}\left(\tau\right)$ $\displaystyle=$
$\displaystyle\sum_{\ell\geq 0}\gamma_{2\ell+1}\,P_{2\ell+1}(\tau)\;.$ (A22)

The free part can be multiplied by an arbitrary constant, and the amplitudes
$\gamma_{2\ell+1}$ of the forced part can be in principle computed from the
decomposition of the source term $S_{3}(\tau)$. We did not carry out this
calculation here, but we have verified that the odd-multipole terms
$\upsilon_{2\ell+1}(r)$, comprising the general MONDian potential
$u(r,\theta)$ in (23), have the corresponding $\propto 1/r^{3}$ asymptotic
behavior after the $\propto 1/r$ part has been eliminated. The free solution
$G_{3}^{\rm free}$ projects onto the even-multipole terms in $u(r,\theta)$.
However, it is easy to see that the contribution to the monopole is zero.
Therefore, the monopole radial function $\upsilon_{0}(r)$ from (23) must not
contain either the $\propto 1/r^{2}$ or the $\propto 1/r^{3}$ asymptotics. We
have verified that this conclusion holds in in all our numerical solutions.

More tests.– Other tests follow from integrating the field equation (12) over
a sphere of arbitrary radius $r$ (see also Eq. 31 in Blanchet & Novak 2011):
the right-hand side is equal to the source mass $M$ multiplied by a constant
factor $-4\pi G$, while the left-hand side may be transformed to an integral
over the surface of the sphere using Gauss’s theorem. If ${\cal S}_{r}$ is
used to denote the surface of a sphere with radius $r$, we have

$\int_{{\cal S}_{r}}\mu\left(g/a_{0}\right)\,{\bf g}\cdot d{\bf
S}=-4\pi\,GM\;,$ (A23)

or

$r^{2}\int_{-1}^{1}d\left(\cos\theta\right)\,\mu\left(g/a_{0}\right)\,\left[\sum_{\ell=0}^{L}\upsilon^{\prime}_{\ell}\left(r\right)\,P_{\ell}\left(\cos\theta\right)+g_{\rm
e}\cos\theta\right]=-2\,GM\;.$ (A24)

We may compute the integral in the left-hand side of Eq. (A24) at the
boundaries of all the radial zones used in our solution (typically 17 of
them), and each time verify how accurately Eq. (A24) is satisfied; we find a
fractional accuracy better than $10^{-5}$.

##  Appendix B Technical details about the numerical solution of the MOND
equation

We first focus on the determination of the source term $\sigma(u,\rho)$ in
Eqs. (19) and (20). Following Blanchet & Novak (2011) we numerically solve for
the MOND potential $u$ outside the Sun only, therefore in the region where the
mass density is zero ($\rho=0$). With ${\bf g}=\nabla u+{\bf g}_{\rm e}$,
orientation of the external field along the $z$-axis, and axisymmetry of $u$,
we have ($\tau=\cos\theta$)

$g^{2}=g_{\rm e}^{2}+\left(\partial_{r}u+2\tau g_{\rm
e}\right)\partial_{r}u+\frac{1-\tau^{2}}{r^{2}}\left(\partial_{\tau}u+2rg_{\rm
e}\right)\partial_{\tau}u\;,$ (B1)

and similarly

$g^{i}\partial_{i}g^{2}=\left(\partial_{r}u+\tau g_{\rm
e}\right)\partial_{r}g^{2}+\frac{1-\tau^{2}}{r^{2}}\left(\partial_{\tau}u+rg_{\rm
e}\right)\partial_{\tau}g^{2}\;.$ (B2)

The source term on the right-hand side of (20), given $u=u_{n}$ at some stage
of the iterative process, reads

$\sigma\left(u,0\right)=-\frac{\mu^{\prime}}{a_{0}\mu}\frac{g^{i}\partial_{i}g^{2}}{2g}\;,$
(B3)

with both $\mu$ and $\mu^{\prime}$ functions of $x=g/a_{0}$. Choosing $\mu(x)$
from one of the family $\\{\mu_{n}(x)\\}$ in Eq. (4), we finally obtain

$\sigma\left(u,0\right)=-\frac{1}{1+(g/a_{0})^{n}}\frac{g^{i}\partial_{i}g^{2}}{2g^{2}}\;.$
(B4)

Thus if we know $u(r,\tau)$ we can use Eqs. (B1), (B2) and (B4) to compute
$\sigma(r,\tau)$ on a large grid of $r$ and $\tau$. All necessary partial
derivatives are determined analytically using the Legendre and Chebyshev
representations of the latitudinal and radial dependence.

The iterative solution of the potential function $u$ (Eq. 19, Sec. 2.2.1)
starts with the Newtonian model ($u=u_{\rm N}=GM/r$) of a mass monopole at the
origin immersed in an asymptotically homogeneous field ${\bf g}_{\rm e}$ along
the $z$-axis of the coordinate system (contributing to $U_{\rm N}$ by a term
${\bf g}_{\rm e}\cdot{\bf r}$: $U_{\rm N}=u_{\rm N}+{\bf g}_{\rm e}\cdot{\bf
r}$). This simple approximation defines a unique point $P$ on the positive
$z$-axis at which the acceleration vanishes. Its distance from the origin is
$r_{\rm e}=r_{\rm M}/\sqrt{\eta}=\sqrt{GM/g_{\rm e}}$, with $r_{\rm M}$ from
Eq. (3) and $\eta=g_{\rm e}/a_{0}$ as in Sec. 2.2.2. Define auxiliary
functions ${\cal Z}(r,\tau)=(r_{\rm e}/r)^{2}-\tau$,
$\zeta_{1}(r,\tau)=1-\tau^{2}+{\cal Z}^{2}$ and
$\zeta_{2}(r,\tau)=1-\tau^{2}-2{\cal Z}^{2}$. Then the source term $\sigma$
from Eq. (B4) with $u=u_{\rm N}$ reads

$\sigma\left(u_{\rm N},0\right)=\frac{g_{\rm e}}{r}\left(\frac{r_{\rm
e}}{r}\right)^{2}\frac{1}{1+\left(\eta\sqrt{\zeta_{1}}\right)^{n}}\,\frac{\zeta_{2}}{\zeta_{1}}\;.$
(B5)

The right-hand side of (B5) – the ratio $\zeta_{2}/\zeta_{1}$, in particular –
is well-behaved except at the critical point $P$, where it exhibits a
singularity. As the iterations proceed, a critical phenomenon around $P$
develops, resulting eventually in an integrable singularity (compare with
similar discussion in Milgrom 1986b; Penner 2020b; Brown & Mathur 2023). This
phenomenon causes some loss of accuracy in our numerical solution for the
potential, especially for sharp transition functions such as $\mu_{10}(x)$;
see Fig. 5. Nevertheless, we find this loss of accuracy is limited to a
relatively small region near $P$, and we do not believe that its affects our
main conclusions.

The solution of the Poisson equation (20) is performed using spectral
decomposition, which proceeds in two steps. First we represent the angular
part of the potential in a series of Legendre polynomials, Eq. (22). Next, we
develop all radial functions (the potential and source) using series of
Chebyshev polynomials, such as

$\upsilon^{(n)}_{\ell}(r)=\sum_{0}^{N}\upsilon^{(n)}_{\ell k}T_{k}(\xi)$ (B6)

and

$S^{(n)}_{\ell}(r)=\sum_{0}^{N}s^{(n)}_{\ell k}T_{k}(\xi)\;,$ (B7)

(the relation between $\xi$ and $r$ is described in Sec. 2.2.1); $N$ is the
adopted maximum degree of the development. For radial zones of stretching
factor $3$ (ratio between outer and inner boundary), we found $N\simeq 99$ is
sufficient. The coefficients $s^{(n)}_{\ell k}$ were computed using Gaussian
quadrature.

Using the spectral decomposition, the differential operator on the left-hand
side of Eq. (22) is represented with an upper triangular matrix ${\bf
M}_{\ell}$, and the whole differential equation becomes a matrix/algebraic
relation between the coefficients $\upsilon^{(n)}_{\ell k}$ and $s^{(n)}_{\ell
k}$ (useful details can be found in Grandclément & Novak 2009). The matrices
${\bf M}_{\ell}$ are singular, due to the additional degrees of freedom
expressed by an arbitrary solution of the homogeneous equation (22). This
singularity is constrained by the boundary conditions at $\xi=\pm 1$ of each
domain where we impose continuity of the potential and its radial derivative.
In the innermost radial domain, i.e. inside the Sun, we require $\partial
u/\partial r=0$ at the origin $r=0$ (Eq. 17), and in the outermost radial zone
we require $u=0$ at $r\rightarrow\infty$. In order to implement the set of all
these boundary conditions, we use the tau method discussed by Grandclément &
Novak (2009, Sec. 2.5.2). After this operation is done, the matrices ${\bf
M}_{\ell}$ are no longer singular, and are easily invertible using singular
value decomposition (e.g., Press et al. 2007).

In order to allow for the large dynamic range of the numerical values of the
variables manipulated during the solution, and to ensure that rounding errors
do not disturb the method, we conservatively use quadruple precision for
representing all real quantities.

Figure 21: The same as in Fig. 17, but now for the simulation ${\cal
S}2^{\prime}$, which includes (i) a static configuration of the giant planets,
the same as seen at present, and (ii) a limited population of $50,000$ test
particles, initially located just exterior to Neptune. Here again the Oort
peak, consisting of LPCs with semimajor axes $\geq 10,000$ au, is too small
compared to the population of LPCs with semimajor axes $\leq 10,000$ au.
Figure 22: The same as Fig. 18, but now for the simulation ${\cal
S}2^{\prime}$.

##  Appendix C Simulation with MOND and no planetary migration

The core set of simulations discussed in Sec. 4 was designed to test MOND
within a scenario that conforms to our best understanding of the evolution of
the solar system. This includes the assumptions that the source of the LPCs
was a planetesimal disk initially exterior to the giant planets, which was
dispersed when Uranus and Neptune migrated outward. While many lines of
evidence indicate that planetary migration must have happened early in the
history of the solar system (e.g. Nesvorný 2018), the details of the migration
history and the architecture of the planetesimal disk are still debated.
Therefore, one may ask whether our conclusions about MOND depend on a
particular set of initial conditions for the planetesimal disk and a
particular model of planetary migration. Here we test this issue. Instead of
considering variants of the planetary migration history we simply assume that
all of the giant planets are fixed on their present orbits, and place the
comets initially in a narrow range of semimajor axes just outside Neptune.
This is our model ${\cal S}2^{\prime}$.

The giant-planet heliocentric orbits were taken from JPL ephemerides at epoch
J2000.0. Since the terrestrial planets were not included, we performed a
barycentric correction and added their composite mass to the Sun. We used
$50,000$ massless particles (a.k.a. comets or planetesimals) placed initially
in a compact and dynamically cold disk just exterior to Neptune – this is much
smaller than the population of $10^{6}$ particles used in our core
simulations, but large enough for our purposes here. Their semimajor axes were
uniformly spread between $31$ and $34$ au, the eccentricities and inclinations
were drawn from a Rayleigh distribution with scale parameters $0.05$ and
$2^{\circ}$, and all other angular variables were distributed uniformly random
between $0^{\circ}$ and $360^{\circ}$. We then followed the system for $4.5$
Gyr using the same methods as for our other simulations (Sec. 3). The MOND
theory was the same as in simulation ${\cal S}2$ , that is, we used the sharp
transition function $\mu_{10}(x)$ and a external Galactic field with
acceleration $g_{\rm e}=2.32\times 10^{-10}$ m s-2.

We found that the results are nearly identical to those of simulation ${\cal
S}2$ (Sec. 4) indicating that the details of the initial distribution of
planetesimals and the planetary migration history do not strongly affect the
currently observed population of LPCs. We illustrate this conclusion by
showing the distribution of LPC binding energies for comets with perihelia
$q\leq 4$ au (Figs. 21 and 22). These figures are to be compared to Figs. 17
and 18. While the details are slightly different, and the noise in simulation
${\cal S}2^{\prime}$ is larger because of the smaller population of test
particles, the principal message is the same. We conclude that the mismatch
between the orbital parameters of the observed LPCs and the predictions of
MOND cannot be ascribed to particular assumptions about the initial spatial
distribution of the planetesimal disk or the planetary migration history.

---
Figure 23: The aphelion directions of nearly parabolic LPCs in our simulations
${\cal S}0$ (top), ${\cal S}2$ (middle), and ${\cal S}1$ (bottom) are shown by
red symbols in the Galactic coordinate system. This is a Mollweide projection,
in which the direction to the Galactic center is at the origin and the north
Galactic pole is at the top. The plot shows LPCs from our simulations with:
(i) perihelion distance $q\leq 4$ au, (ii) inverse value of the original
semimajor axis $a$ in the range $0.1\leq 10^{4}/a\leq 0.7$ au-1, (iii) first
visit to the solar system (“new comet”). The data are collected from the last
$500$ Myr of our simulations. The black triangles show the aphelion direction
for 81 class-1 LPCs from the Warsaw catalog using similar selection criteria
(see text). The histograms on both sides of the figures show the latitude
distribution of the aphelia, for observed LPCs on the left, in black, and for
the simulated data on the right, in gray. The histogram at the top shows the
longitude distribution of the observed LPCs. The bars at the bin center
indicate a formal Poissonian uncertainty of the value.

##  Appendix D Distribution of aphelion directions of nearly parabolic LPCs

Finally, we briefly analyze the distribution of aphelion directions for the
original121212As usual, the “original” orbital elements are the barycentric
orbital elements that the comets had at a heliocentric distance of $250$ au,
before the orbits are influenced by planetary perturbations. Note that the
aphelion direction is much less sensitive to planetary perturbations than the
semimajor axis. orbits of nearly parabolic LPCs. We use this subclass of LPCs
as their extremely large aphelion distances make them vulnerable to exterior
perturbations. While encounters with individual stars, if strong enough (i.e.,
high-mass star and/or small impact parameter and/or slow relative velocity),
naturally leave their imprint in the LPC population over a few Myr interval
(e.g., Feng & Bailer-Jones 2014, 2015), some authors argued that there might
also be persistent anisotropies. For instance, Matese et al. (1999) and Matese
& Whitmire (2011) claimed to identify an anomalous concentration of LPC
aphelia along a great circle located at Galactic longitudes $-45^{\circ}$ and
$135^{\circ}$, and that this concentration was evidence for a distant
planetary companion of a few Jupiter masses, residing in the Oort cloud. While
it is not our intention to discuss in detail such scenarios, we bring up his
topic here because Paučo & Klačka (2017) suggested that a similar effect may
be produced by MOND itself.

Here we use the simulations analyzed in Sec. 4 and describe the aphelion
distributions they predict. Similar to the study of the semimajor axis
distribution, we compare our results with data in the Warsaw catalog of nearly
parabolic LPCs (Królikowska & Dybczyński 2020). In order to minimize
contamination by uncertain orbits and to focus on the largest orbits, we
select: (i) class 1 solutions, (ii) with perihelia $q\leq 4$ au, and (iii)
with binding energies $1/a$ in the range $0.1\leq 10^{4}/a\leq 0.7$ au-1
(these represent the bulk of the Oort peak). We find 81 comets that satisfy
these criteria in the Warsaw catalog. The distribution of aphelion directions
in Galactic coordinates is shown by black triangles in Fig. 23. Without
conducting a statistical analysis, we believe that the evidence for a
concentration about the great circle discussed in Matese et al. (1999) (a
polar circle along longitudes $135^{\circ}$ and $315^{\circ}$) is rather weak:
the maxima in the corresponding longitude bins, shown in the histogram at the
top of the figure, are not more than 1$\sigma$ above the level set by a
uniform distribution. However, here we focus more on comparison of our
simulations.

The upper panel on Fig. 23 shows results from the reference model ${\cal S}0$
which contains only Newtonian perturbations (red circles). The longitude
distribution is fairly uniform as expected. There is, however, a clear
latitudinal variation with maxima at $\pm 45^{\circ}$ galactic latitude (see
the gray histogram on the right). A similar variation in seen in the
observational data (black histogram on the left), although the signal is
biased by the uneven geographical distribution of comet discoveries and the
difficulty of finding comets near the Galactic plane. This is evidence for the
substantial role of the smooth Galactic tide in injecting new comets from the
Oort cloud to the inner solar system (e.g., Delsemme 1987; Rickman 2010). The
remaining panels of Fig. 23 show the results from simulations ${\cal S}2$
(middle) and ${\cal S}1$ (bottom). These are based on MONDian dynamics, with
two forms of the transition function $\mu(x)$. While there are some
longitudinal variations in the distribution, especially in the case of the
${\cal S}1$ simulation that uses the gradual $\mu_{2}(x)$ transition function,
neither simulation shows the level of aphelion anisotropy that was described
in Paučo & Klačka (2017). We believe the methodology in that paper was too
simplified to realistically predict the angular distribution of the observed
comets.