Astronomy & Astrophysics manuscript no.

aanda ©ESO 2025

September 9, 2025

Dynamical Constraints on a Population of Massive Interstellar

Objects
Oem Trivedil1 and Abraham Loeb2,

1
Department of Physics and Astronomy, Vanderbilt University, Nashville, TN 37235, USA
e-mail: oem.trivedi@vanderbilt.edu
2
Astronomy Department, Harvard University, 60 Garden St., Cambridge, MA 02138, USA
e-mail: aloeb@cfa.harvard.edu
arXiv:2509.06300v1 [astro-ph.EP] 8 Sep 2025

September 9, 2025

ABSTRACT

Context. The detection of kilometer-scale interstellar objects (ISOs) places strong dynamical constraints on their underlying popula-
tion and distribution.
Aims. We aim to assess whether the observed encounter rates of large ISOs can be explained by their expected velocity distribution,
or whether anomalous velocity anisotropies are required.
Methods. We first derive the cumulative encounter rate scaling with size from brightness-limited detection models. We then employ
the Eddington inversion method to link power-law density profiles ρ(r) ∝ r−k with phase-space distributions, illustrating how steep
slopes imply strong velocity-space focusing. Finally, we develop a Liouville mapping formalism based on energy and angular momen-
tum conservation, propagating the interstellar distribution inward to the detection region while explicitly incorporating gravitational
focusing and anisotropy.
Results. We show that encounter rates of large ISOs require flux enhancements far beyond expectations, that steep density profiles
produce the necessary inward velocity bias, and that Liouville mapping provides a physically self-consistent way to reproduce the ob-
served size dependent detection rates. The main results are framed in the context of the parameters for 3I/ATLAS, but the implications
are general and go on to sharpen the distinction between natural dynamical mechanisms and potential artificial origins for ISOs.
Key words. Interstellar Objects – 3I/ATLAS – Non-Natural Trajectories

1. Introduction velocity, and specific trajectory presents significant dynamical

challenges. Unlike smaller ISOs which may be ejected more
The discovery of the first two interstellar objects (ISOs), 1I/’ readily from planetary systems during early dynamical instabil-
Oumuamua Meech et al. (2017) in 2017 and 2I/Borisov in 2019 ities, the presence of massive interstellar bodies at observable
Guzik et al. (2020) marked a turning point in planetary science, number densities like this strains standard models of ejection
providing the first direct probes of small bodies formed outside efficiency and population statistics. As with 1I/‘Oumuamua
the Solar System. Despite their common interstellar provenance and 2I/Borisov, we are hence led to consider the detection
these objects displayed strikingly divergent properties, with of 3I/ATLAS in a broader framework. Much work has hence
1I/’Oumuamua being photometrically inactive yet exhibited followed to understand various properties of this object Loeb
statistically significant non-gravitational acceleration, while (2025); Hibberd et al. (2025); Loeb et al. (2025); De La
2I/Borisov showing unmistakable cometary activity. These Fuente Marcos et al. (2025); Hopkins et al. (2025b), as was
contrasting characteristics ignited a broad range of theoretical the case for 1I/’Ouamuamua Loeb (2022); Forbes & Loeb
explanations. (2019); Siraj & Loeb (2022, 2019a); Bialy & Loeb (2018).One
of the central questions is whether the appearance of a large
The recent discovery of 3I/ATLAS Seligman et al. (2025), object within the inner solar system is consistent with expected
with an orbital eccentricity e ∼ 6.1, an inclination of 175◦ , and dynamical pathways, or does it point to the need for a revised
a heliocentric excess velocity of V∞ ∼ 58 km s−1 , adds a third understanding of either the ISO size distribution or their delivery
case to this growing sample of ISOs. Early observations reveal mechanisms?
that 3I/ATLAS hosts a faint but detectable coma, yet exhibits
minimal rotational light curve variation and a moderately red In this work, we build upon the population level argu-
spectral slope of ∼ 17%/100 nm. Its absolute magnitude near ments articulated in the rarity of 3I/ATLAS observation in
HV ≈ 12.4 suggests a nucleus radius in the range below 2.8 Loeb (2025) and quantify the dynamical constraints implied
km Jewitt et al. (2025) and up to 23 km Lisse et al. (2025), by the detection of 3I/ATLAS. Specifically, we evaluate the
depending on the relative contribution of the coma and the phase-space conditions and angular momentum distributions
nucleus to the brightness. Altogether, 3I/ATLAS represents the that would allow an object of its inferred size to be delivered
largest interstellar object observed to date. through the observed trajectory, and we calculate the likelihood
of such events under both natural and artificial population
The combination of its large inferred size, high incoming hypotheses. By considering both the small radius and large

Article number, page 1 of 6

 A&A proofs: manuscript no. aanda

radius interpretations of 3I/ATLAS, we demonstrate how the Assuming spatial uniformity, the number of objects per unit vol-
rarity of its trajectory places stringent constraints on the abun- ume in a size range [R, R + dR] is
dance of massive interstellar objects. This analysis provides a
dN
new perspective on the emerging census of ISOs and directly n(R)dR = , (8)
addresses the question of whether the detection of 3I/ATLAS dV
is dynamically compatible with a natural origin, or whether and so
alternative explanations such as artificial or directional delivery
need to be considered. n(R)dR ∝ R−q dR, (9)
which implies the differential number density scales as
2. Encounter rate analysis
n(R) ∝ R−q . (10)
Since direct detections of ISOs remain rare, the statistical rate of
encounters provides a powerful indirect diagnostic of the under- Now, rather than compute the differential detection rate, we
lying ISO population. In particular, encounter rate formulations derive the cumulative detection rate which is the rate of en-
allow us to quantify the scaling of detection probabilities with countering all ISOs with radius greater than or equal to R. To
object size, velocity distribution, and survey sensitivity. This is do so, we integrate the differential rate from R to ∞. Since
especially interesting in the case of an object like 3I/ATLAS, as Γ1 (R′ ) ∝ n(R′ ) · σ(R′ ), and we have just shown that n(R′ ) ∝ R′−q
it could highlight the enhancements in flux required to recon- and σ(R′ ) ∝ R′ , the integrand becomes R′−(q−1) . The cumulative
cile the observed ISO detections with natural population models. encounter rate is
Z ∞
We begin with the standard expression for the detection Γ1 (R) ∝ R′−(q−1) dR′ . (11)
rate of ISOs, written as an encounter rate R

Γ1 (R) = n(R) · ⟨v⟩ · σ(R), (1) This evaluates (for q > 2) to

where n(R) is the number density of objects within a narrow bin Γ1 (R) ∝ R−(q−2) . (12)
dR of radius R, ⟨v⟩ is the average encounter speed of interstel- Finally, including the average encounter speed ⟨v⟩, which is
lar objects relative to the Sun, and σ(R) is the effective detection computed from the shifted Maxwellian distribution as
cross section for an object of radius R. Rather than modeling the
cross section geometrically, we instead adopt a brightness-based 2
!
v2⊙
!
detection model, following a similar formulation as in Seligman ⟨v⟩ = √ σ s 1 + , (13)
et al. (2025). In this approach, the detection cross section is set π 6σ2s
by the maximum distance at which an object of size R is de- with σ = 25 km/s v⊙ = 17.4 km/s. We have taken the velocity
tectable, given a fixed limiting apparent magnitude mlim . Using dispersion σ s to be of this value, given recent results showing
standard photometric relations for reflected sunlight, the appar- that the dispersion in the range 20 − 30 km/s regime is reason-
ent magnitude of an object at heliocentric distance r ≫ 1 au, able. For example, Kohandel et al. (2020) showed a velocity dis-
geocentric distance ∆, and radius R, for a fixed geometric albedo persion in the region 23 − 28km/s region, Eubanks et al. (2021)
p, is given approximately by discussed a value of 26.14 for local stellar population as well,
while other recent works have also also considered this regime
m ≈ H + 5 log10 (r∆), (2)
Hopkins et al. (2025b,a). The final expression for the cumulative
where the absolute magnitude H itself depends on the object ra- detection rate becomes
dius as
 p  Γ1 (R) = k · R−(q−2) · ⟨v⟩. (14)
H = 17.1 − 2.5 log10 − 5 log10 (R). (3) To fix the normalization, we impose the condition for 3I/ATLAS,
0.04
Assuming constant albedo and phase behavior, we can express Γ1 (0.6) = 0.2 yr−1 (15)
the apparent magnitude as
which gives
m ∝ −5 log10 (R) + 5 log10 (r∆). (4)
0.2
Fixing the apparent magnitude to the detection threshold mlim , k= . (16)
(0.6)q−2 · ⟨v⟩
we can then invert the above expression to solve for the maxi-
mum distance at which an object of size R can be detected. This Thus, the full expression becomes
leads to  R −(q−2)
0.2
rmax (R) ∝ R 1/2
, (5) Γ1 (R) = · (0.6)q−2 · R−(q−2) · ⟨v⟩ = 0.2 · . (17)
⟨v⟩ 0.6
and hence, since the detection cross section scales with the ob- This allows us to define the flux enhancement factor as
servable area (a circle of radius rmax ), we have
Γ2  R q−2
σ(R) ∝ rmax
2
∝ R. (6) ≈ , (18)
Γ1 0.6
To determine the number density n(R), we assume that the ISOs which gives the required enhancement in flux beyond the natural
follow a differential size distribution of the form population in order to explain the detection rate of objects larger
dN than a given radius R. This formulation is based on cumulative
∝ R−q . (7) detection rate considerations and incorporates both the size
dR
Article number, page 2 of 6
 Oem Trivedil and Abraham Loeb: Dynamical Constraints on a Population of Massive Interstellar Objects

under the Sun’s gravitational potential Ψ(r) = GM⊙ /r. This

method self-consistently accounts for gravitational focusing and
allows us to understand what kind of velocity distribution is
needed to reproduce the enhanced detection rate of ISOs like
3I/ATLAS near the Sun.

We begin by assuming a power-law density profile for the

incoming ISO population Napier et al. (2021); Siraj & Loeb
(2019b); Peñarrubia (2023) in the gravitational potential of the
Sun
GM⊙ dΨ GM⊙
Ψ(r) = so that =− 2 . (19)
r dr r
We now use the Eddington inversion formula
Fig. 1. Cumulative encounter rate enhancement factor as a function of
d ρ dΨ
"Z E 2 #
minimum ISO radius R for different power-law slopes q, showing the 1
steep increase in required enhancement for larger objects. The adopted f (E) = √ √ , (20)
parameters reflect representative values for 3I/ATLAS. 8π2 0 dΨ2 E − Ψ
where E = Ψ(r) − 12 v2 is the relative energy per unit mass. We
distribution and the average velocity derived from the velocity now invert ρ(Ψ) ∝ Ψk by re-parametrizing our assumed ρ(r) in
distribution. terms of Ψ, then perform the differentiation and integration to
compute f (E).
To determine the proportionality constant k, we use the re-
sult from Loeb (2025) for 3I/ATLAS and require that the We assume a power-law density profile
cumulative detection rate for all objects of radius greater
GM GM
than R0 = 0.6 km matches the observed detection rate of ρ(r) ∝ r−k , with Ψ(r) = ⇒r= Ψ (21)
Γ1 (0.6) = 0.2 yr−1 . r .
Figure 1 shows the cumulative encounter rate enhancement Substitute into ρ(r)
demonstrating how the expected detection frequency of interstel-
lar objects depends on their radius and on the slope of the size  GM −k
distribution q. Note that we have considered a wide range of pos- ρ(Ψ) ∝ ∝ Ψk , (22)
Ψ
sible radii for 3I/ATLAS, from ∼ 2.8km Jewitt et al. (2025) to
23 km Lisse et al. (2025). While the natural component of the and so
encounter rate sets the baseline expectation, the curves reveal
dρ d2 ρ
that for larger radii the required enhancement factor rises steeply ∝ kΨk−1 , ∝ k(k − 1)Ψk−2 . (23)
with increasing q and in particular for steep size distributions dΨ dΨ2
(q ≳ 4), objects in the multi-kilometer regime demand encounter We this plug this into the Eddington formula to get
rates that exceed natural expectations by orders of magnitude.
This might suggest an artificial or non-Maxwellian component Z E
dΨ
in the detection statistics of large ISOs such as 3I/ATLAS Selig- f (E) ∝ Ψk−2 √ (24)
man et al. (2025). 0 E−Ψ
Let Ψ = Eu ⇒ dΨ = Edu, and E − Ψ = E(1 − u) and then
3. Eddington-Inversion Analysis in Phase-Space Z 1
f (E) ∝ Ek−3/2
uk−2 (1 − u)−1/2 du (25)
The Eddington inversion method provides a powerful framework 0
to connect spatial density profiles with the underlying velocity
distribution functions in a self-consistent manner Binney & This integral yields a Beta function solution as
Tremaine (2011); Lacroix et al. (2018); Suárez et al. (2010).
Γ(a)Γ(b)
Z 1
For ISOs this approach is particularly useful as it naturally ua−1 (1 − u)b−1 du = B(a, b) = . (26)
incorporates gravitational focusing by the Sun and reveals 0 Γ(a + b)
how an observed radial overdensity can emerge from specific The final result comes out to be (as the Beta function just acts as
velocity-space structures. By inverting an assumed power-law a constant)
density profile, one directly obtains the phase-space distribution
function f (E), which determines the energy dependence of ISO f (E) ∝ Ek−3/2 . (27)
orbits near the Sun.
This yields a direct analytic expression for the phase-space dis-
This not only allows us to identify the kinds of velocity tribution function as a function of the relative energy E, parame-
biases required to reproduce enhanced detection rates, such as terized by k. This f (E) gives the required energy distribution for
those observed for 3I/ATLAS, but also provides a means of the ISO population that results in a radial density enhancement
distinguishing between shallow, weakly focused populations ρ(r) ∝ r−k and since this formulation assumes spherical symme-
and steep, radially concentrated ones. We use the Eddington try, the enhancement we get is entirely due to velocity-space bias
inversion method to determine the phase-space distribution towards the Sun. Figure 2 shows the Eddington-inverted phase-
function f (E) that would yield a radial density profile ρ(r) ∝ r−k space density f (v = 30 km/s, r) as a function of log10 (r/AU) for
Article number, page 3 of 6
 A&A proofs: manuscript no. aanda

with gravitational focusing effects. In contrast, the Eddington

inversion presumes an isotropic, steady-state bound population,
making it best suited for illustrative explorations of how dif-
ferent slopes k bias the velocity distribution. However, for the
boundary-driven ISO problem, Liouville mapping provides the
more self-consistent and physically motivated route.

We formulate the phase space analysis in terms of the constants

of motion for the Kepler problem which are the specific orbital
energy E and the specific angular momentum J. This change
of coordinates greatly improves numerical stability because
the relevant flux of low angular momentum trajectories can
be represented explicitly. The starting point is the distribution
function at infinity, which we model as a Maxwellian shifted
by the solar motion and for simplicity we adopt an isotropic
Maxwellian with one-dimensional velocity dispersion σ, so that
the velocity distribution at infinity is
Fig. 2. Comparison of Eddington-inverted phase-space densities for k =
3, 3.5, 4 with a shifted Maxwellian background. v2
!
f∞ (v) ∝ v2 exp − 2 .
2σ
different power-law density profiles alongside the value expected Transforming variables to the specific orbital energy
from a shifted Maxwellian background
E = 12 v2∞ , (29)
v2
!
f (v) ∝ v exp − 2 .
2
(28)
2σ we obtain an effective one dimensional energy density
√
!
Adopting for Sun’s velocity relative to the Local Standard E
of Rest v = 30 km/s and σ = 25 km/s yields f (v) ∼ ϕ(E) ∝ E exp − 2 . (30)
σ
167.56 (km/s)−1 . Note that this is just the background expected
ISO flux at that velocity, assuming an isotropic, unperturbed dis- The second conserved quantity is the specific angular momen-
tribution. Figure 2 illustrates that the curves for k = 3, 3.5, 4 all tum J = rv⊥ . At the outer boundary of the system, rout , the maxi-
exhibit a strong enhancement in f (v, r) near the Sun (small r) mum allowed angular momentum for a given E is Jmax ≃ rout v∞ .
and fall off rapidly at large distances. This behavior indicates Closer to the Sun, at a given radius r, the kinematic constraint
a phase-space overdensity consistent with a population of ISOs from v2r ≥ 0 imposes
concentrated on radial orbits toward the inner solar system. Such p
steep k-values are physically motivated when trying to explain J ≤ Jcap (E, r) = r 2(E + GM⊙ /r). (31)
the detection of large ISOs. By contrast, shallower profiles like
Thus, the allowed domain in (E, J) space is compact and
k = 1 or 2 imply only weak gravitational focusing and fail to
bounded by Jmax and Jcap .
generate the necessary enhancement in f (v, r), underestimating
the encounter rate at small heliocentric distances. The relation
To represent the anisotropy required to enhance low angu-
f (E) ∝ Ek−3/2 means that larger values of k require a stronger
lar momentum trajectories we introduced a weighting function
concentration of low-angular-momentum orbits. Therefore, the
gR (J) which suppresses high-J orbits relative to low-J ones. The
steep Eddington profiles better reflect the degree of inward di-
family of functions we employed has the form
rectional bias required to explain the observational data.
!α(R,a)
J
gR (J) ∝ 1 − for J ≤ J0 (R), (32)
4. Directional Bias J0 (R)
While the Eddington inversion method provided a useful and gR (J) = 0 otherwise and here J0 (R) sets the cutoff scale and
framework for exploring the connection between idealized is parameterized as
density profiles ρ(r) ∝ r−k and their corresponding velocity !−s
distributions, it is less naturally suited for problems where the R
J0 (R) = J0,base , (33)
boundary distribution is known explicitly at large radii. For the R0
present analysis, where we aim to propagate the interstellar  
Maxwellian distribution inward from rout ∼ 105 AU to the de- while the exponent is taken as α(R, a) = a RR0 − 1 , where
tection region near rmax = 4 AU, a Liouville mapping approach a is a tunable anisotropy amplitude calibrated to match the
is more appropriate. Liouville’s theorem guarantees that the observational constraints. For α = 0, gR (J) is flat, while for
phase space density is conserved along collisionless trajectories α > 0, the function favors orbits with lower J.
and this allows us to start from the physically meaningful
boundary condition f (v, rout ) = fISM (v) and map each velocity Given this phase space distribution, the density of parti-
vector inward using the constants of motion, which in this cles crossing a spherical shell at radius r for a population of
case are orbital energy E and angular momentum J Miller objects of physical radius R is
(1988); Lindblad (1959). This not only makes the boundary Z Emax Z Jmax (E,r)
condition explicit, but also directly incorporates the dynamics J
of unbound, hyperbolic orbits characteristic of ISOs, together ρ(r; R, a) ∝ dE ϕ(E) dJ gR (J) , (34)
0 0 vr (r, E, J)
Article number, page 4 of 6
 Oem Trivedil and Abraham Loeb: Dynamical Constraints on a Population of Massive Interstellar Objects

with radial velocity that it is possible to match the observational size distribution
r slope by tuning the anisotropy in J, there is still no known

J2 astrophysical mechanism that would actually generate such a
vr (r, E, J) = 2 E + GM⊙ /r − 2 . (35) distribution of low angular momentum interstellar objects in
r
the first place. This gap leaves room for speculation that the
To extract the effective slope k(R) that measures the relative required anisotropy could be the result of artificial intervention,
enhancement of encounter rates with object size, we compute like from a directed release of objects with engineered orbits.
ρ(r; R, a) in a narrow logarithmic window around the detection
radius rmax = 4 AU and fit a local power law
5. Conclusions
ρ(r; R, a) ∝ r−k(R) . (36) In this work we have developed a dynamical framework to con-
This procedure gives an estimate of k(R) by comparing ρ at radii strain the population of large ISOs using three complementary
slightly above and below rmax and since only relative values approaches. Encounter rate analysis demonstrated that detection
matter, we normalize such that k(R0 ) = 0 at R0 = 0.6 km, in of ISO with radii of tens of kilometers cannot be reconciled
accordance to Loeb (2025). with a purely Maxwellian background, implying the need for
significant flux enhancements. By applying the Eddington
We then adjust the anisotropy amplitude a to be in agree- inversion method, we established how steep radial density
ment with the required enhancement implied by the cumulative profiles naturally correspond to phase-space distributions biased
detection rates and once a is fixed, the function k(R) follows toward radial orbits, providing a clear dynamical mechanism for
uniquely from the dynamics. The resulting plot in Figure 3 of gravitational focusing near the Sun. Extending this discussion,
we then employed a Liouville mapping in (E, J) space to prop-
agate the interstellar velocity distribution inward under explicit
conservation laws, showing that the required size dependent
enhancement of encounter rates can arise from anisotropies that
favor low angular momentum trajectories.

These results suggest that the detection of ISOs with R ≳ 10 km

is dynamically consistent with strong gravitational focusing, but
only under specific velocity-space biases whose natural astro-
physical origin remains unknown. This tension suggests two
possible explanations, one being that new physical mechanisms
like stellar perturbations or local structure in the Galactic envi-
ronment produce the required velocity anisotropies. The other
possibility is an orbital fine-tuning for ISOs of technological
origin. Future ISO discoveries by the Vera C. Rubin Observatory
Thomas et al. (2020); Blum et al. (2022); Sebag et al. (2020)
will provide the crucial statistics needed to test these dynamical
Fig. 3. k(R) as a function of ISO radius R, taking into account Liouville characteristics.
mapping formalism
Acknowledgements. OT was supported in part by the Vanderbilt Discovery Al-
liance Fellowship. AL was supported in part by the Black Hole Initiative, which
k(R) shows a smooth, monotonic rise from zero at R = 0.6 km is funded by GBMF anf JTF.
to k(10) ≈ 4 at R = 10 km. We can extend it to R = 23km
for the very large size mentioned in Lisse et al. (2025), but the
trend remains the same. Physically, this curve quantifies how
a bias toward low angular momentum trajectories, consistent References
with Liouville’s theorem, can reproduce the cumulative size Bialy, S. & Loeb, A. 2018, The Astrophysical Journal Letters, 868, L1
dependence of the observed interstellar object detection rate. Binney, J. & Tremaine, S. 2011, Galactic dynamics (Princeton university press)
The fact that the curve is regular and monotonic indicates Blum, B., Digel, S. W., Drlica-Wagner, A., et al. 2022, arXiv preprint
arXiv:2203.07220
that the physics is internally consistent, and that the earlier De La Fuente Marcos, R., Alarcon, M., Licandro, J., et al. 2025, Astronomy &
numerical difficulties in direct velocity space integration were Astrophysics, 700, L9
purely artifacts of coordinate choice. Eubanks, T. M., Hein, A. M., Lingam, M., et al. 2021, arXiv preprint
arXiv:2103.03289
Forbes, J. C. & Loeb, A. 2019, The Astrophysical Journal Letters, 875, L23
This curve makes an important point for the artificial ori- Guzik, P., Drahus, M., Rusek, K., et al. 2020, Nature Astronomy, 4, 53
gin discussion as it shows that within a physically consistent Hibberd, A., Crowl, A., & Loeb, A. 2025, arXiv preprint arXiv:2507.12213
Liouville mapping in E, J space, one can in fact generate the Hopkins, M. J., Bannister, M. T., & Lintott, C. 2025a, The Astronomical Journal,
necessary enhancement in the encounter rate with increasing 169, 78
object size purely by favoring low-angular momentum orbits. Hopkins, M. J., Dorsey, R. C., Forbes, J. C., et al. 2025b, arXiv preprint
arXiv:2507.05318
The bias needed to explain why large interstellar objects are Jewitt, D., Hui, M.-T., Mutchler, M., Kim, Y., & Agarwal, J. 2025, The Astro-
detected at the observed rate can arise naturally from dynamical physical Journal Letters, 990, L2
phase–space effects, without having to invoke any special Kohandel, M., Pallottini, A., Ferrara, A., et al. 2020, Monthly Notices of the
mechanisms. The monotonic rise of k(R) to values ∼ 4 is a Royal Astronomical Society, 499, 1250
Lacroix, T., Stref, M., & Lavalle, J. 2018, Journal of Cosmology and Astroparti-
signature that gravitational dynamics can in principle account cle Physics, 2018, 040
for the observed scaling. But at the same time, the artificial Lindblad, B. 1959, in Astrophysik IV: Sternsysteme/Astrophysics IV: Stellar
scenario remains compelling as while the mapping shows Systems (Springer), 21–99

Article number, page 5 of 6

 A&A proofs: manuscript no. aanda

Lisse, C., Bach, Y., Bryan, S., et al. 2025, arXiv preprint arXiv:2508.15469
Loeb, A. 2022, Astrobiology, 22, 1392
Loeb, A. 2025, Research Notes of the AAS, 9, 178
Loeb, A., Hibberd, A., & Crowl, A. 2025, arXiv preprint arXiv:2507.21402
Meech, K. J., Weryk, R., Micheli, M., et al. 2017, Nature, 552, 378
Miller, R. 1988, Celestial mechanics, 45, 19
Napier, K. J., Adams, F. C., & Batygin, K. 2021, The Planetary Science Journal,
2, 53
Peñarrubia, J. 2023, Monthly Notices of the Royal Astronomical Society, 519,
1955
Sebag, J., Claver, C. F., Thomas, S. J., et al. 2020, in Ground-based and Airborne
Telescopes VIII, Vol. 11445, SPIE, 396–415
Seligman, D. Z., Micheli, M., Farnocchia, D., et al. 2025, The Astrophysical
Journal Letters, 989, L36
Siraj, A. & Loeb, A. 2019a, arXiv preprint arXiv:1904.07224
Siraj, A. & Loeb, A. 2019b, The Astrophysical Journal Letters, 872, L10
Siraj, A. & Loeb, A. 2022, New Astronomy, 92, 101730
Suárez, D. O., Rubio, L. B., Vögler, A., & Del Toro Iniesta, J. 2010, Astronomy
& Astrophysics, 518, A2
Thomas, S. J., Barr, J., Callahan, S., et al. 2020, in Ground-based and Airborne
Telescopes VIII, Vol. 11445, SPIE, 68–82

Article number, page 6 of 6