# On the contribution of the 40K geo-antineutrino to single Borexino events.
Version 2

L.B. Bezrukov, I.S. Karpikov, A.S. Kurlovich, A.K. Mezhokh,
S.V.Silaeva, V.V. Sinev and V.P. Zavarzina

###### Abstract

We propose to include in the analysis of Borexino single event energy spectrum
the scattering of 40K geo-antineutrinos by scintillator electrons. The
Hydridic Earth model predicts the concentration of potassium in modern Earth
from 1% to 4% of the Earth mass. We calculated contribution of 40K geo-
antineutrino interactions in single Borexino events for these concentrations.
This contribution is comparable to the contribution from the interaction of
CNO neutrinos. We discuss the reasons for using the Hydridic Earth model.

Institute for Nuclear Research of Russian academy of sciences, Moscow

##  1 Introduction.

Recently the Borexino collaboration [1] published the results of analysis of
the experimental energy spectrum of single events. In this work the authors
fitted the experimental spectrum with the theoretical one. As a result, the
interaction rates from various sources were obtained except for CNO solar
neutrinos interaction rate.

To obtain the values of the interaction rates from all sources, it is
necessary to have the correct theoretical information about the shape of the
measured energy spectrum from each of the possible sources. In the absence of
correct information, the used procedure for fitting of experimental data with
the model may not allow to obtain the interaction rates and their
uncertainties from the sources used in the model.

The results of [1] most likely demonstrate the lack of complete information on
possible sources of single events. The paper reports that there is no way to
get the interaction rate of CNO neutrinos. Only the upper bound of $8\div 10$
events per day in 100 tons of scintillator was obtained. From 3 to 5 events
per day are expected from CNO neutrino interactions depending on the model of
the Sun.

##  2 Energy spectrum of recoil electrons from scattering of 40K geo-
antineutrino.

The authors of the work [2] considered an additional source of single events
for the Borexino detector. This is the scattering of potassium geo-
antineutrinos by electrons. Such a source was not considered in [1] because
the Silicate Earth model was used in it. The concentration of potassium in
this model is very small (0.024% of the Earth mass) and the contribution to
single events of the Borexino detector from potassium geo-antineutrinos is
negligible.

Fig.1. Differential energy spectra of recoil electrons from neutrino (or
antineutrino) scattering on electrons for various sources (the oscillations
were not taken into account). The red line is the solar neutrino from 13N
decay in CNO cycle. The green is the solar neutrino from 15O decay in CNO
cycle. The blue is the solar neutrino from 17F decay in CNO cycle. The black
is 40K geo-antineutrino for the concentration of potassium in modern Earth 1%
of the Earth mass

Fig.2. Differential energy spectra of recoil electrons from solar CNO neutrino
scattering by scintillator electrons in 100 tons of scintillator (red curve)
and the sum of the recoil electron spectra from the scattering of solar CNO
neutrino and from the scattering of 40K geo-antineutrino for the concentration
of potassium in the modern Earth from 1% to 2% of the Earth mass. The
oscillations were taken into account.

We proceeded from the prediction of the Earth elemental composition by another
model of the Earth ”Hydridic Earth” or ”Hydrogen-rich Earth” [3, 4]. These
studies predict concentration of potassium in the original Earth up to 3.8% of
the Earth mass. In the modern Earth we can expect a larger proportion of
potassium due to the loss of hydrogen by the Earth during the time of it’s
existence. The predictions of the Hydridic Earth model have low accuracy, so
it is reasonable to search for the flux of potassium geo-antineutrinos in a
wide range of values.

The work [2] presented the results of calculation of the energy spectrum of
recoil electrons from geo-antineutrino scattering from 40K decay on
scintillator electrons and the counting rate of these events in 100 tons of
scintillator for the concentration of potassium in the modern Earth from 0.3%
to 3.8% of the Earth’s mass.

We present at figure 1 our calculation (using the approximations described in
[2] and the oscillations were not taken into account) of differential energy
spectra of recoil electrons produced by the scattering of neutrinos (or
antineutrinos) on electrons for different sources: red curve - solar neutrinos
from the decay of 13N in the CNO cycle, the green curve corresponds to the
solar neutrinos from the decay of 15O in the CNO cycle, blue curve - solar
neutrinos from the decay of 17F in the CNO cycle, black curve - geo-
antineutrinos from 40K decay for the concentration of potassium in the Earth
1% of weight of the Earth. We can see that the spectrum of 40K geo-
antineutrinos for the concentration of potassium in the Earth of 1% of the
Earth mass is similar to the spectrum of solar neutrinos from the decay of 13N
in the CNO cycle. It has an upper limit larger than that of the spectrum of
solar neutrinos from the decay of 13N. For our calculation we used spectra of
solar neutrinos from the CNO cycle [6] — original work of J. Bahcall, which
are normalized to the value of the flux for large metallicity of the Sun
($4.88\cdot 10^{8}$ cm-2 s-1).

Figure 2 shows the differential energy spectrum from CNO neutrino scattering
by electrons in 100 tons of scintillator (red line) and the total recoil
electron spectrum from the scattering of CNO neutrino and 40K geo-antineutrino
for the concentration of potassium in the modern Earth from 1% to 2% of the
Earth mass (green and blue curves respectively). The oscillations were taken
into account. It means that the contribution of muon and tau neutrinos (and
anty-neutrinos) scattering was taken into account.

The intensity of events from the CNO neutrino calculated for the curve from
figure 2 is 4.9 cpd/100 tons, and the additional intensities of events from
40K geo-antineutrinos is 1.6 to 3.4 cpd/100 tons. The Borexino collaboration
in [1] used a differential energy spectrum of recoil electrons from the
scattering of CNO neutrinos on electrons, similar to the one shown in figure 2
in red, and did not obtain expected value of interaction rate ($3\div 5$
cpd/100 tons).

We propose to perform the same study as in Borexino collaboration work [1]
using the total spectrum of recoil electrons from the scattering of CNO
neutrinos and from the scattering of 40K geo-antineutrino for the potassium
concentration in the modern Earth from 1% to 2% of the Earth’s mass. We expect
to find a certain value for the intensity of total events from the CNO
neutrinos and from the 40K geo-antineutrinos in the range of $6\div 9$ events
per day in 100 tons of scintillator. We also expect that, as a result of this
analysis, several events per day will be taken away from events attributed in
work [1] to solar 7Be neutrinos interactions and 210Bi decays. The errors of
the obtained intensities will become symmetrical.

We give the functions from figure 2 in table form (see table below).

The obtaining the intensity of total spectrum of recoil electrons from the
scattering of CNO neutrinos and 40K geo-antineutrinos as a result of fitting
analysis of Borexino single event data will be an important confirmation of
the validity of the Hydridic Earth model.

Table 1: Differential energy spectra of recoil electrons from solar CNO neutrino scattering by scintillator electrons in 100 tons of scintillator and the same with addition of 40K antineutrinos in units $MeV^{-1}year^{-1}(100tons)^{-1}$.  T, MeV | CNO |  CNO+1%(40K) |  CNO+1.5%(40K) |  CNO+2%(40K)
---|---|---|---|---
0.1 | 2997.64 | 5254.39 | 6334.59 | 7414.79
0.2 | 2736.43 | 4363.73 | 5142.65 | 5921.57
0.3 | 2447.10 | 3597.83 | 4148.63 | 4699.42
0.4 | 2129.82 | 2924.45 | 3304.78 | 3685.11
0.5 | 1787.78 | 2319.67 | 2574.24 | 2828.80
0.6 | 1439.26 | 1781.15 | 1944.77 | 2108.38
0.7 | 1100.80 | 1308.16 | 1407.39 | 1506.62
0.8 | 787.671 | 901.87 | 956.52 | 1011.16
0.9 | 526.847 | 579.58 | 604.81 | 630.04
1.0 | 355.736 | 371.46 | 378.98 | 386.5
1.1 | 249.9 | 249.9 | 249.9 | 249.9
1.2 | 156.71 | 156.71 | 156.71 | 156.71
1.3 | 81.17 | 81.17 | 81.17 | 81.17
1.4 | 26.75 | 26.75 | 26.75 | 26.75
1.5 | 0.303 | 0.303 | 0.303 | 0.303

##  3 Reasons to use the Hydridic Earth model.

The widespread belief in the fairness of Silicate Earth model and belief in
the validity of the results of work [5] that the heat flux from the Earth
interior is equal to 47$\pm$2 TW do not allow to include in analysis a new
source of single events in Borexino - the 40K geo-antineutrinos interaction.

The idea that the Hydridic Earth model is not correct is also widely accepted.
One of the objections is that the entire Earth could have melt due to
radiogenic heat and spent most part of its life in this state if the Earth
contains potassium more than 1% of the Earth weight. The current heat flux
from the Earth’s interior can be more than 200 TW in the frame of Hydridic
Earth model which contradicts to the result of the work [5].

We have learned how to respond to the objections to Hydridic Earth model. The
answers to the main objections are published in [7, 8, 9, 10]. Here we note
that the entire Earth could not have melt because the Hydridic Earth model
contains a subsurface cooling mechanism that is activated when the subsurface
is heated enough to decompose the metal hydrides. Therefore, in the Hydridic
Earth model, the subsurface temperature oscillates [3].

It is also noted in [7, 8, 9, 10] that thermal conductivity is not the main
mechanism of heat transfer in the Earth, but protons and hydrogen-containing
gases carry out the heat away. In these works the experimental evidences are
provided that the heat flux from Earth interior can reach the several hundreds
TW. These are the heating of the oceans, the temperature profile of ultra-deep
wells and non-direct evidence – Moon heat flux from interior.

Moreover, the effect was found that is predicted by the Hydridic Earth model
and not predicted by the Silicate Earth model. The Hydridic Earth model
predicts that the Earth’s crust is positively charged and contains a large
amount of positive charge in the form of protons and positive ions of various
hydrogen-containing gases. We tested the validity of this prediction
experimentally by detecting an excess of the concentration of positive air-
ions in underground rooms over the concentration of negative air-ions. In
order to make this prediction, a new model of terrestrial electricity was
developed based on the Hydridic Earth model [11]. The model named ”Hydridic
model of terrestrial electricity”. This model explains the origin of the
atmospheric electric field and all the observed effects of atmospheric
electricity in a single way. The model also explains the origin of telluric
currents. In [12], the Hydridic model of terrestrial electricity was
successfully used in the analysis of the reaction of telluric currents to
earthquakes.

We consider Hydridic Earth model as an adequate tool for the study of the
Earth.

##  4 Conclusion.

  1. 1.

We propose to include the scattering of 40K geo-antineutrinos by electrons in
the analysis of single event energy spectrum of Borexino detector.

  2. 2.

The Hydridic Earth model predicts the potassium concentrations up to 4% of the
Earth mass.

  3. 3.

It is expected that the total interaction rate of the solar CNO neutrino and
40K geo-antineutrino will be obtained in the range of $6\div 9$ events per day
in 100 tons of scintillator.

  4. 4.

We consider Hydridic Earth model as an adequate tool for the study of the
Earth.

##  5 Acknowledgment

The authors express their gratitude to M. D. Skorokhvatov for the opportunity
to discuss this note at a seminar at the NRC ”Kurchatov Institute”.

## References

  * [1] M. Agostini, K. Altenmuller, S. Appel et al., Simultaneous precision spectroscopy of pp, 7Be, and pep solar neutrinos with Borexino Phase-II, Phys. Rev. D 100 (2019) 082004, doi:10.1103/PhysRevD.100.082004; arXiv:1707.09279 [hep-ex].
  * [2] V.V. Sinev, L.B. Bezrukov, E.A. Litvinovich, I.N. Machulin, M.D. Skorokhvatov, Looking for Antineutrino Flux from 40K with Large Liquid Scintillator Detector, Physics of Particles and Nuclei 46 (2015) 186, doi:10.1134/S1063779615020173; arXiv:1405.3140 [physics.ins-det].
  * [3] V.N. Larin, Hydridic Earth: the New Geology of Our Primordially Hydrogen-Rich Planet. / Ed. C. Warren Hunt. Calgary, Alberta, Canada: Polar Publishing, 1993. 247p.
  * [4] Herve Toulhoat, Valerie Beaumont, Viacheslav Zgonnik, Nikolay Larin, Vladimir N. Larin, Chemical differentiation of planets: a core issue, arXiv:1208.2909 [astro-ph.EP].
  * [5] J.H. Davies and D.R. Davies, Earth’s surface heat flux, Solid Earth 1 (2010) 5, doi:10.5194/se-1-5-2010.
  * [6] J. N. Bahcall and R. K. Ulrich, Rev. Mod. Phys. 60, 297 (1988). J. N. Bahcall, https://www.sns.ias.edu/$\sim$ jnb/SNdata/cnospectra.html (2005).
  * [7] L.B. Bezrukov, A.S. Kurlovich, B.K. Lubsandorzhiev, A.K. Mezhokh, V.P. Morgalyuk, V.V. Sinev and V.P. Zavarzina, How Geoneutrinos can help in understanding of the Earth heat flux, J. Phys. Conf. Ser. 934 no.1 (2017) 012011, doi:10.1088/1742-6596/934/1/012011.
  * [8] L. B. Bezrukov, A. S. Kurlovich, B. K. Lubsandorzhiev, A. K. Mezhokh, V. P. Morgalyuk, V. V. Sinev, and V. P. Zavarzina, Geo-Neutrinos and the Earth’s Internal Heat Flux. Physics of Particles and Nuclei, 2018, Vol. 49, No. 4, pp. 674–677. Original Russian Text © L.B. Bezrukov, A.S. Kurlovich, B.K. Lubsandorzhiev, A.K. Mezhokh, V.P. Morgalyuk, V.V. Sinev, V.P. Zavarzina, 2018, published in Fizika Elementarnykh Chastits i Atomnogo Yadra, 2018, Vol. 49, No. 4, pp. https://doi.org/10.1134/S1063779618040135
  * [9] I.R. Barabanov, L.B. Bezrukov, V.P. Zavarzina, I.S. Karpikov, A.S. Kurlovich, B.K. Lubsandorzhiev, A.K. Mezhokh, V.P. Morgalyuk, V.V. Sinev, Study of Earth’s Heat Flux by Means of Geoneutrino Detection. Published in Phys.Atom.Nucl. 82 (2019) no.1, 8-12
DOI: 10.1134/S1063778819010034

  * [10] L. B. Bezrukov, A. S. Kurlovich, B. K. Lubsandorzhiev, V. V. Sinev, V. P. Zavarzina and V. P. Morgalyuk. Geo-neutrino, Earth heat flux, Earth electricity. EPJ Web of Conferences 191, 03005 (2018) QUARKS-2018.
DOI: https://doi.org/10.1051/epjconf/201819103005

  * [11] L.B. Bezrukov, V.P. Zavarzina, A.S. Kurlovich, B.K. Lubsandorzhiev, A.K. Mezhokh, V.P. Morgaluk, V.V. Sinev, On Negatively Charged Layer of the Earth’s electric Field //Doklady Physics, 2018, Vol. 63, No. 5, pp. 177–179. Original Russian Text © L.B. Bezrukov, V.P. Zavarzina, A.S. Kurlovich, B.K. Lubsandorzhiev, A.K. Mezhokh, V.P. Morgaluk, V.V. Sinev, Doklady Akademii Nauk, 2018, Vol. 480, No. 2, pp. 155–157. https://doi.org/10.1134/S1028335818050051
  * [12] L. B. Bezrukov, I. S. Karpikov, A. S. Kurlovich, B. K. Lubsandorzhiev, A. K. Mezokh, V. P. Morgaluk, V.V. Sinev and V. P. Zavarzina. INTERPRETATION OF RESULTS OF MEASUREMENTS OF VOLTAGE INTO BAYKAL LAKE. Geomagnetism and Aeronomy. 2019. V.59. No 5. pp 623-627. DOI: 10.1134/S0016793219040054 Original Russian Text © L.B. Bezrukov, V.P. Zavarzina, I.S.Karpikov, A.S. Kurlovich, B.K. Lubsandorzhiev, A.K. Mezhokh, V.P. Morgaluk, V.V. Sinev, Geomagnetizm I Aeromomiya, 2019, Vol. 59, No. 5, pp. 666-670.
DOI: 10.1134/S0016794019040059


# Can Cold Jupiters Sculpt the Edge-of-the-Multis?

Nicole Sobski  Department of Computer Science, Wellesley College, Wellesley,
MA 02481, USA  Kavli Institute for Astrophysics and Space Research,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA  Sarah C.
Millholland Kavli Institute for Astrophysics and Space Research, Massachusetts
Institute of Technology, Cambridge, MA 02139, USA  Department of Physics,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA
<EMAIL_ADDRESS>
###### Abstract

Compact systems of multiple close-in super-Earths/sub-Neptunes (“compact
multis”) are a ubiquitous outcome of planet formation. It was recently
discovered that the outer edges of compact multis are located at smaller
orbital periods than expected from geometric and detection biases alone,
suggesting some truncation or transition in the outer architectures. Here we
test whether this “edge-of-the-multis” might be explained in any part by
distant giant planets in the outer regions ($\gtrsim 1$ AU) of the systems. We
investigate the dynamical stability of observed compact multis in the presence
of hypothetical giant ($\gtrsim 0.5\ M_{\mathrm{Jup}}$) perturbing planets. We
identify what parameters would be required for hypothetical perturbing planets
if they were responsible for dynamically sculpting the outer edges of compact
multis. “Edge-sculpting” perturbers are generally in the range $P\sim 100-500$
days for the average compact multi, with most between $P\sim 200-300$ days.
Given the relatively close separation, we explore the detectability of the
hypothetical edge-sculpting perturbing planets, finding that they would be
readily detectable in transit and radial velocity data. We compare to
observational constraints and find it unlikely that dynamical sculpting from
distant giant planets contributes significantly to the edge-of-the-multis.
However, this conclusion could be strengthened in future work by a more
thorough analysis of the detection yields of the perturbing planets.

citnum

##  1 Introduction

Short-period super-Earths/sub-Neptunes are the most prevalent class of
observed extrasolar planets, orbiting around roughly half of Sun-like stars
(e.g. Petigura et al., 2013; Zhu et al., 2018; He et al., 2019). First
discovered by early Doppler surveys (e.g. Mayor et al., 2011) and then, to a
greater degree, by NASA’s Kepler Mission (Borucki et al., 2010), short-period
planets are often found in tightly-spaced configurations, with multiple
planets all having orbital periods in the range from days to months (Lissauer
et al., 2011, 2014; Rowe et al., 2014; Fabrycky et al., 2014). These so-called
“compact multiple-planet systems” (or “compact multis”) exhibit a remarkable
degree of structure and regularity in both their orbital and physical
properties. Their orbits are nearly coplanar and circular (e.g. Fang & Margot,
2012; Fabrycky et al., 2014; Van Eylen & Albrecht, 2015; Xie et al., 2016).
Moreover, their period ratios, radii, and masses exhibit a statistical
tendency towards uniformity within a given system (Weiss et al., 2018;
Millholland et al., 2017). This set of patterns is known collectively as the
“peas-in-a-pod patterns” or “intra-system uniformity” (for a review, see Weiss
et al. 2022).

Despite the abundance of information the Kepler Mission delivered about the
architectures of close-in super-Earths/sub-Neptunes, it was mostly limited to
the inner regions ($a\lesssim$ 1 AU) of the systems. However, various
observational and theoretical efforts are gradually revealing more about their
outer architectures. For instance, using radial velocity (RV) observations to
search for long-period companions, Zhu & Wu (2018) and Bryan et al. (2019)
showed evidence that distant giant planets (or “cold Jupiters”, with $a\gtrsim
1$ AU, $M_{p}\gtrsim 0.5\ M_{\mathrm{Jup}}$) are over-represented in systems
with inner super-Earths. They found that $\sim 30-40\%$ of systems with inner
super-Earths are found to have distant giant planets, relative to the $\sim
10\%$ distant giant occurrence for stars irrespective of small planet
presence. The super-Earth/distant giant correlation was also studied by
Rosenthal et al. (2022) in their analysis of the California Legacy Survey
(Rosenthal et al., 2021); they found tentative agreement with the earlier
studies, showing that stars with inner small planets may have an enhanced
occurrence of outer giants with $1.7\sigma$ significance.

A contradictory result was recently presented by Bonomo et al. (2023), who
studied RVs of 38 Kepler and K2 small-planet systems collected over nearly a
decade with the HARPS-N spectrograph, as well as publicly available
measurements collected with other facilities. They detected five cold Jupiters
in three systems and derived an occurrence rate of $9.3^{+7.7}_{-2.9}\%$ for
planets with $0.3-13\ M_{\mathrm{Jup}}$ and $1-10$ AU in systems with inner
small planets. Bonomo et al. (2023) found no evidence of the overabundance of
distant giant planets in super-Earth systems claimed by Zhu & Wu (2018) and
Bryan et al. (2019). However, Zhu (2023) recently suggested that the
discrepancy between the Bonomo et al. (2023) result and previous works can be
fully resolved by accounting for the metallicity dimension of the super
Earth–cold Jupiter relations. Further constraints may come from the Kepler
Giant Planet Survey (KGPS, Weiss et al. 2023), a decade-long survey of 63
Kepler systems using HIRES at the Keck Observatory designed to search for
long-period planets in Kepler systems. Weiss et al. (2023) presented RV-
detected companions to 20 stars in their sample, with occurrence rate analyses
forthcoming. It is clear that the relationship between inner super-Earths and
outer giants is still debated but gradually coming into better focus.

Distant giant planets in general often have moderate to high eccentricities,
$e\sim 0.1-0.8$ (e.g. Jones et al., 2006; Udry & Santos, 2007; Wright et al.,
2009; Dawson & Murray-Clay, 2013), and their orbits might be significantly
mutually-inclined with respect to the inner systems (e.g Gratia & Fabrycky,
2017). In some cases, an eccentric and/or inclined distant perturber can exert
a stronger gravitational influence on the inner planets than their mutual
gravitational coupling, and the inner system can gain significant
eccentricities and inclinations, potentially driving dynamical instabilities
(e.g. Lai & Pu, 2017; Hansen, 2017; Pu & Lai, 2018; Denham et al., 2019;
Spalding & Millholland, 2020). The dynamics of the inner system can thus be
tightly coupled to that of the outer system.

Given this dynamical interaction, if a system has both small inner planets and
a distant giant planet, it must experience an architectural transition point
in its overall layout. The transition point could be generated by the
perturbative influence of the distant giant planet. Specifically, each system
with an eccentric and possibly inclined distant giant planet must possess a
“zone of instability”, where any planets existing in that region would
experience eccentricity/inclination excitation and destabilization. It is
possible that the impacts of this are already detectable.

One potential link is the result of Millholland et al. (2022), who presented
statistical evidence for an outer truncation or transition (named the “edge-
of-the-multis”) in the observed high-multiplicity Kepler systems. They
explored the existence of hypothetical additional super-Earths/sub-Neptunes
orbiting beyond the outermost planets in Kepler compact multis, under the
assumption that the “peas-in-a-pod” patterns continue to larger separations
than observed. They then estimated the transit and detection probabilities of
these hypothetical planets and found that additional exterior planets should
have been detected in $\gtrsim 35\%$ of Kepler compact multis. Accordingly, to
$\gtrsim 7\sigma$ confidence, the Kepler compact multis truncate at smaller
orbital periods than expected from geometric and detection biases alone. These
results can be explained if compact multis experience an architectural
truncation (i.e. an occurrence rate decrease) or a transition (i.e. a
breakdown in the “peas-in-a-pod” patterns) at $\sim 100-300$ days.

Although various theories predict an “edge-of-the-multis” consistent with
Millholland et al. (2022)’s observation (e.g. Chatterjee & Tan, 2014; Zawadzki
et al., 2022; Batygin & Morbidelli, 2023), it is still unknown what physical
mechanism(s) are most responsible. In this paper, we explore the possibility
of a link between the outer edges of compact multis and distant giant planets.
That is, we test the hypothesis that dynamical perturbations from distant
giant planets are strong enough to sculpt the outer edges. We explore this
with the following steps: (1) We first consider the observed properties of
compact multis and identify where hypothetical distant giant planets with
various properties could exist without causing dynamical instabilities. (2)
Next, we determine the required properties of the distant giant planets if
they were sculpting the outer edges. (3) Finally, we compare the properties of
these “edge-sculpting” distant giants to observational constraints. It is
important to note that distant giants cannot be the only explanation of the
edge-of-the-multis, since there are not enough of them in compact multi
systems (even if they were as prevalent as suggested by Zhu & Wu 2018 and
Bryan et al. 2019). Our aim here is to address whether they could play any
significant role in the edge-of-the-multis.

This paper is organized as follows. We begin by describing our planet sample
and our methodology for assessing dynamical stability of compact multis in the
presence of distant giant planets (Section 2). We then present the results of
these calculations, both of individual systems and the whole population
(Section 3). This allows us to determine the properties of perturbing planets
that would be capable of sculpting the outer edges of the compact multis. We
explore the detectability of these perturbing planets and the resulting
implications for the edge-of-the-multis (Section 4). We discuss alternative
theories for the edge-of-the-multis and conclude in Section 5.

##  2 Methods

###  2.1 Planet sample

Figure 1: Planet sample. We display the architectures of the 64 Kepler compact systems with four or more transiting planets described in Section 2.1. The dot size is proportional to the planet size, as shown in the legend at top. Table 1: Parameter space scope. Parameters used in the dynamical stability calculations and their sampling scheme. Parameter | Inner planets | Outer perturber
---|---|---
Mass | $M_{p}\sim\mathrm{Normal}(\overline{M}_{p},\sigma_{M_{p}})$ | $M_{p,p}\sim\mathrm{Uniform}[0.5\ M_{\mathrm{Jup}},5\ M_{\mathrm{Jup}}]$
Semi-major axis |  $a=$ observed | $a_{p}\sim\mathrm{Uniform}[a(P_{p}=1.3P_{\mathrm{out}}),a(P_{p}=1000\ \mathrm{days})]$
Eccentricity | $e\sim\mathrm{Rayleigh}(0.04)$ | $e_{p}\sim\mathrm{Uniform}[0,0.7]$
Inclination | $i\sim\mathrm{Rayleigh}(1^{\circ})$ |  $i_{p}\sim\mathrm{Rayleigh}(10^{\circ})$ (Case 1); $i_{p}\sim\mathrm{Uniform}[0^{\circ},40^{\circ}]$ (Case 2)
Argument of periapse | $\omega\sim\mathrm{Uniform}[0^{\circ},360^{\circ}]$ | $\omega_{p}\sim\mathrm{Uniform}[0^{\circ},360^{\circ}]$
Longitude of ascending node | $\Omega\sim\mathrm{Uniform}[0^{\circ},360^{\circ}]$ | $\Omega_{p}\sim\mathrm{Uniform}[0^{\circ},360^{\circ}]$
Mean anomaly | $M\sim\mathrm{Uniform}[0^{\circ},360^{\circ}]$ | $M_{p}\sim\mathrm{Uniform}[0^{\circ},360^{\circ}]$

We begin by defining our sample of Kepler compact multiple-planet systems,
which is identical to the sample studied in Millholland et al. (2022). We
consider systems with four or more observed transiting planets only, since it
is a well-defined high multiplicity sample that is most suitable to
investigations of the outer edges. As our starting point, we use the Kepler
DR25 KOI catalog (Thompson et al., 2018; NASA Exoplanet Archive, 2022) and
consider all planets with “confirmed” and “candidate” dispositions. Where
possible, we replace the stellar parameters and planet radii in the DR25
catalog with parameters from the Gaia-Kepler Stellar Properties Catalog
(Berger et al., 2020a, b). We apply two quality cuts. First, we consider only
planets smaller than $16\ R_{\oplus}$ with fractional radius uncertainties
less than 100%. Second, we discard targets for which Furlan et al. (2017)
found a companion star that contributed more than 5% of the light in the
photometric aperture. For the purposes of this study, we consider only systems
with four or more observed transiting planets, which represent a high fidelity
sample of compact multis. We are left with 279 planets in 64 Kepler systems
with four or more observed transiting planets. The planet sample is displayed
in Figure 1.

Most planets in these systems do not have measured masses due to observational
limitations. (For instance, the star is too faint for Doppler follow-up and/or
there are no detectable transit timing variations.) However, in cases where
they are available (60 of 279), we utilize observational mass estimates from
radial velocity or transit timing variation analyses, which we obtain from the
NASA Exoplanet Archive (2023). For the remaining planets, we obtain
approximate mass estimates from the Forecaster probabilistic mass-radius
prediction model by Chen & Kipping (2017). Forecaster is an open-source
package that predicts a missing planetary mass or radius based on a mass-
radius relationship established from a sample of 316 objects of various sizes
spanning small planets to late-type stars. Altogether, our approach provides
us with a heterogeneous collection of mass estimates, which may affect our
results at the detailed level. However, this approach is sufficient for our
purposes insofar as we are obtaining an average constraint across the compact
multis as a population.

###  2.2 Dynamical stability calculations without perturbers

We first assess the dynamical stability of the compact multi-planet systems
before introducing additional perturbing planets. These estimates will serve
as our basis of comparison to the later calculations. We perform our dynamical
stability analyses using SPOCK (Stability of Planetary Orbital Configuration
Klassifier; Tamayo et al. 2020), a machine learning model that predicts
dynamical stability over long timescales ($\sim 10^{9}$ orbits) by first
running a short integration ($10^{4}$ orbits), assessing summary statistics,
and classifying the stability using the XGBoost machine learning algorithm.
SPOCK provides dynamical stability estimates at a rate that is five orders of
magnitude faster than direct $N$-body integrations. Here, we will utilize the
SPOCK-derived probability of stability, which we denote $p_{\mathrm{stab}}$.

To estimate $p_{\mathrm{stab}}$, SPOCK requires the masses and orbital
elements as inputs. Our planet sample only has measured stellar masses,
planetary masses (most of which are estimated using Forecaster, as described
previously), and periods, so we must randomly sample the remaining parameters.
We take the orbital eccentricities and inclinations to be Rayleigh-distributed
with scale parameters consistent with current observational constraints (Xie
et al., 2016; Van Eylen et al., 2019; Mills et al., 2019). We take the
arguments of periastron, longitudes of the ascending node, and mean anomalies
to be uniformly distributed from $0^{\circ}$ to $360^{\circ}$. Because the
planet masses are subject to uncertainties (especially in the case where the
masses are predicted from the Forecaster model), we include random sampling of
the masses. We take the masses to be normally distributed around the mean
observed or estimated value, $\overline{M}_{p}$, with standard deviations
$\sigma_{M_{p}}$ equal to the average of the observational uncertainties
$\sigma_{M_{p}}=(\sigma_{M_{p},\mathrm{high}}+\sigma_{M_{p},\mathrm{low}})/2$
in the cases where the mass estimates are from RVs or TTVs, or equal to
$0.3\overline{M}_{p}$ in the cases where the mass estimates are from
Forecaster. The details of the parameter sampling are summarized in Table
LABEL:tab:_parameter_space_scope. This sampling is performed for each planet
in each system in our sample. New random variables are drawn for each
iteration of dynamical stability calculations, as will be described in further
detail in the next section.

For each system, we estimate $p_{\mathrm{stab}}$ 100 times, with each
iteration having a different set of random variables. We then find the mean
probability of stability, $\overline{p}_{\mathrm{stab}}$ and consider that to
be representative of the system. Although a majority of the compact multi-
planet systems are classified as stable, there are 20 systems with
$\overline{p}_{\mathrm{stab}}<0.5$, with 6 of these systems demonstrating
$\overline{p}_{\mathrm{stab}}<0.1$. These systems are preferentially near
mean-motion resonances.111Considering the subset of unstable systems, we found
that 45% of adjacent planet pairs are near a first-order mean-motion
resonance, compared to 15% of planet pairs within the subset of stable
systems. Moreover, 95% of the unstable systems contain at least one near-
resonant pair compared to 40% of the stable systems. This is unsurprising
given that we are randomly sampling the mean anomalies, and planets can become
destabilized near resonances if they not placed in specific parts of their
orbits. Additional unstable systems that are not near-resonance may have a
heightened sensitivity to mass variability (e.g. the real masses may be
significantly smaller than the estimated masses), and the combination of this
and the randomization of the other orbital elements could yield instability.
We disregard these systems in our subsequent analysis, but we believe our
analysis is still robust as a population average over many systems. We are
more interested in the influence of distant perturbing planets on the ensemble
of compact multis than on individual systems.

###  2.3 Dynamical stability calculations with perturbers

Having assessed the dynamical stability of the inner systems in isolation, we
are ready to test the limits of additional outer planets in each stable
system. We again opt to use SPOCK for our stability calculations. However, an
important caveat is that the system architectures with both inner multis and
distant giant planets are somewhat outside of SPOCK’s training models, which
only considered masses up to $0.1\ M_{\mathrm{Jup}}$ and inclinations
$\lesssim 10^{\circ}$ (Tamayo et al., 2020). Although SPOCK has been shown to
generalize reasonably well, it is still important to check its validity in
this context. We performed multiple checks using other numerical and
analytical stability estimates, which we describe in more detail in Appendix
A.

To assess the dynamical impact of the outer perturbing planet, we randomly
sample its physical and orbital properties within reasonable ranges. The outer
perturber is taken to be a massive gas giant planet on an eccentric and
moderately inclined orbit, in line with general observed properties of cold
Jupiters (e.g. Dawson & Murray-Clay, 2013). We sample its mass as
$M_{p,p}\sim\mathrm{Uniform}[0.5\ M_{\mathrm{Jup}},5\ M_{\mathrm{Jup}}]$. We
consider periods, $P_{p}$, between $1.3P_{\mathrm{out}}$ (where
$P_{\mathrm{out}}$ is the outermost observed planet in the inner system) and
$1000$ days. These limits were chosen because $1.3P_{\mathrm{out}}$ is
approximately the smallest possible stable period, and $1000$ days is
sufficiently distant to detect the transition between stable and unstable
orbits. The $1000$ day maximum period also extends beyond the maximum
detectable window via transit detection. We sample uniformly in values of the
perturber’s semi-major axis, $a_{p}$, according to this period range.

The outer perturber’s eccentricity is sampled as
${e_{p}\sim\mathrm{Uniform}[0,0.7]}$, which encompasses the range of
eccentricities of most observed distant giant planets (e.g. Rosenthal et al.,
2021). As for the perturber’s inclination, we consider two different cases. In
Case 1, we consider the perturber to be nearly coplanar with the inner system
and sample $i_{p}\sim\mathrm{Rayleigh}(10^{\circ})$. This is motivated by
Masuda et al. (2020), who suggested that distant giant planets orbiting
compact multis tend to be nearly coplanar with the inner systems. In Case 2,
we consider a wider range of inclinations for the perturber,
$i_{p}\sim\mathrm{Uniform}[0^{\circ},40^{\circ}]$. This range is consistent
with the expectations of energy equipartition between eccentricities and
inclinations. The remaining orbital elements (argument of periastron,
longitude of the ascending node, and mean anomaly) are sampled uniformly from
$0^{\circ}$ to $360^{\circ}$.

Figure 2: Dynamical stability calculations: effect of mass sampling
variations. The scatterplots depict the SPOCK probabilities of stability for
three example systems (KOI-812, KOI-408, and KOI-94) with the addition of an
outer perturbing planet. Each point corresponds to a single stability
calculation with an associated mass and semi-major axis for the outer
perturber. The left column illustrates stability calculations in which the
inner planet masses are fixed, while the right column has the inner planet
masses randomly sampled from a normal distribution as described in Section
2.2. The perturber’s orbital inclination is randomly sampled from a Rayleigh
distribution (Case 1 in Table LABEL:tab:_parameter_space_scope). The colorbar
indicates the SPOCK probability of stability. Blue horizontal lines are
plotted to indicate the semi-major axes of existing planets. The red line
bordered by a translucent gray belt defines the zone corresponding to
“metastability”, the transition from highly unstable to highly stable.

We use the set-up of the inner systems as described in Section 2.2 and use
SPOCK to assess the stability of the augmented systems, consisting of the
observed inner planets and the hypothetical outer perturbers. For each system,
we generate $10,000$ unique parameter combinations and calculate an associated
probability of stability for each. Note that the parameters of both the outer
perturber and inner planets are randomly sampled, as summarized in Table
LABEL:tab:_parameter_space_scope.

We can visualize the general trends resulting from varying the perturber’s
physical characteristics by constructing scatterplots of the stability
variations as a function of the perturber’s parameters. We first create
scatterplots of $a_{p}$ vs. $M_{p,p}$, with $p_{\mathrm{stab}}$ indicated with
the colors of the points. Examples are shown for the KOI-812, KOI-408, and
KOI-94 systems in Figure 2, which we will discuss at greater length in Section
3. There is a noticeable division in the $a_{p}$ vs. $M_{p,p}$ parameter
space, which separates regions of low and high probabilities of stability. In
all cases, the high probabilities max out at values near the probabilities of
stability from the simulations without the perturbers. However, the stable and
unstable regions are not starkly divided; we notice an upward sloping
transition zone with $p_{\mathrm{stab}}\sim 0.4-0.6$, a region that we will
refer to as “metastable”. This observation motivates the next section’s
analysis, which is to calculate exactly where this metastable region exists
for each system.

###  2.4 Metastable transition region

We now describe how we utilize the collection of stability calculations to
extract the set of perturber parameters defining the metastable region. We
approximate the metastable region as a banded line that divides the stable and
unstable regions, and we identify its location using the following numerical
procedure. First, we divide the $a_{p}$ vs. $M_{p,p}$ parameter space into a
grid with 10 evenly-spaced segments in $M_{p,p}$ and 40 evenly-spaced segments
in $a_{p}$, and we calculate the mean probability of stability,
$\overline{p}_{\mathrm{stab}}$, in each grid cell. The number of grid cells
was chosen through trial and error; the spacing ensures that each grid cell is
large enough to include many stability values within it but small enough to
not smear out the relevant features. The results are not strongly sensitive to
these choices.

Next, we use the grid of $\overline{p}_{\mathrm{stab}}$ to find the best-fit
line of metastability. We do this by first iterating through columns in the
grid (each of which have a fixed $M_{p,p}$) and estimating the value of
$a_{p}$ that lies closest to the center of the metastable region. The
metastable region is observed to be a gradual transition between low and high
$\overline{p}_{\mathrm{stab}}$, and the curve of
$\overline{p}_{\mathrm{stab}}$ vs. $a_{p}$ for a fixed $M_{p,p}$ (i.e. a
single column of the grid) is well-fit by a piecewise linear function of the
following form:

$p(a_{p})=\begin{cases}\overline{p}_{\mathrm{stab},l}&a_{p}\leq a_{p,l}\\\
m_{\mathrm{trans}}(a_{p}-a_{p-l})+\overline{p}_{\mathrm{stab},l}&a_{p,l}<a_{p}<a_{p,h}\\\
\overline{p}_{\mathrm{stab},h}&a_{p}\geq a_{p,h}.\end{cases}$ (1)

Here, $a_{p,l}$ and $a_{p,h}$ are parameters of the function defining the low
and high boundaries of the metastable region, $\overline{p}_{\mathrm{stab},l}$
and $\overline{p}_{\mathrm{stab},h}$ are parameters defining the low and high
probabilities at which the curve levels out, and
$m_{\mathrm{trans}}=(\overline{p}_{\mathrm{stab},h}-\overline{p}_{\mathrm{stab},l})/(a_{p,h}-a_{p,l})$
is the slope within the metastable region.

We fit the piecewise linear function to the curve of
$\overline{p}_{\mathrm{stab}}$ vs. $a_{p}$ for each column of fixed $M_{p,p}$.
We then find the value of $a_{p}$ for which $p(a_{p})=0.5$. This value, which
we denote $a_{p}(p=0.5)$, defines the center of the metastable region within
the column. In principle, the edges of the metastable region are defined by
$a_{p,l}$ and $a_{p,h}$. However, we find that $a_{p,l}$ is very near to
$\min(a_{p})$ in most cases, so we define the lower edge of the metastable
region to be $a_{p}(p=0.3)$. We define the upper edge of the metastable region
to be $a_{p}(p=0.7)$. By repeating this process for each column, we obtain a
set of ($M_{p,p},a_{p}$) values defining three lines: the center, lower edge,
and upper edge of the metastable region. Finally, we perform linear
regressions on each curve so as to obtain a smooth demarcation of the
metastable region. The centers and edges of the metastable regions are
indicated with red lines and gray bands for the examples in Figure 2.

##  3 Results

We now present the results of our dynamical stability calculations. The first
thing we explore is the impact of different sampling choices (Section 3.1).
Many parameters are sampled in the exploration of these systems, so it is
important to verify that these choices do not strongly impact the results. We
will then examine the averaged results for all systems (Section 3.2).

###  3.1 Effect of sampling variations

Figure 3: Dynamical stability calculations: effect of inclination sampling
variations. The plotting scheme and example systems are identical to Figure 2,
except we now explore different methods of sampling the outer perturber’s
inclination. The left column indicate uniform sampling
$i_{p}\sim\mathrm{Uniform}[0^{\circ},40^{\circ}]$, and the right indicates
Rayleigh-distributed sampling ${i_{p}\sim\mathrm{Rayleigh}(10^{\circ})}$. In
all cases, the inner planet masses are randomly sampled from a normal
distribution as in the right column of Figure 2.

Figure 4: Dynamical stability calculations: dependence on perturber semi-major
axis vs. periastron distance. The plotting scheme and example systems are
identical to Figure 2 and Figure 3, but now the left column shows the semi-
major axis vs. mass, while the right column shows the perturber periastron
distance vs. mass. In all cases, the inner planet masses are randomly sampled
from a normal distribution as in the right column of Figure 2 and the orbital
inclination of the outer perturber is Rayleigh-distributed.

Figure 2 shows the impact of sampling the inner planet masses, where the left
and right columns show the results when the masses are held fixed and varied,
respectively. Sampling the inner planet masses makes some systems a bit more
stable and other systems a bit more unstable. However, for most systems, there
is very little overall change in the location and width of the metastable
region, suggesting that the results are not strongly sensitive to the inner
planet mass sampling scheme.

Figure 3 shows the impact of the different methods of sampling the outer
perturber’s inclination. The left and right columns show the results of the
uniformly-distributed and Rayleigh-distributed sampling (recall the two cases
in Table LABEL:tab:_parameter_space_scope). As in Figure 2, there are some
slight differences between the two versions, but overall the general features
and the locations and widths of the metastable regions are the same regardless
of the sampling method.

Finally, Figure 4 displays the results for the same three example systems in a
different parameter space. The left column is the same as the right column of
Figure 3, displaying perturber semi-major axis vs. mass. (Here we use
Rayleigh-distributed perturber inclinations.) The right column replaces the
y-axis with the perturber periastron distance, $a_{p}(1-e_{p})$. It is
immediately clear that plotting the results in this space yields a sharper
transition from the unstable to stable region and thus a smaller metastable
region. This makes physical sense because the inner planets are more sensitive
to the perturber’s close approach distance than its semi-major axis.

###  3.2 Combined constraints

Figure 5: Aggregate results of stability calculations. We display the
“metastable regions” (grey belts in Figures 2-4) for all systems in our
sample. The top, middle, and bottom panels show the results in perturber semi-
major axis, periastron distance, and period vs. mass. The depicted metastable
regions are taken from the results in which the inner planet masses are
randomly sampled from a normal distribution around the mean and the orbital
inclination of the outer perturber is Rayleigh-distributed. The darkest blue
indicates where there is greatest overlap among the various systems. In the
bottom plot, we include contours of the outer perturber’s RV semi-amplitude
(curved lines) and, if it was transiting, the probability that it would
transit three or more times in Kepler data (horizontal lines).

We now combine the results for all systems in order to explore the influence
of distant perturbing planets on the ensemble of compact multi-planet systems.
Figure 5 shows the extracted metastable regions for all systems. The darkest
color indicates where there is the greatest overlap among the various systems.
The metastable region is concentrated in the period range $P\sim 100-500$ days
($a\sim 0.4-1.2$ AU for $M_{\star}=M_{\odot}$).

We also display information on the detectability of the perturbing planet
using RV and Kepler transit data. As for RV detection, we plot contours of
constant RV semi-amplitude assuming $M_{\star}=M_{\odot}$ and $e=0$. For the
vast majority of the parameter space, the perturber would be easily detected
in RV data of a sufficiently long baseline, since the semi-amplitude is much
greater than typical RV uncertainties.

As for transit detection, we calculate the conditional probability that,
provided the perturber was transiting, it would transit its host star at least
three times in Kepler data. We obtain this estimate using a data product
provided by Kepler DR25 (Thompson et al., 2018; NASA Exoplanet Archive, 2022)
called the window function, which is the fraction of unique transit ephemeris
epochs that permit three or more transits to be observed as a function of
orbital period (Burke & Catanzarite, 2017). This quantity is complex due to
the data gaps within the limited span of Kepler photometry. Kepler DR25
provided tabulations of window functions for each star as a function of
orbital period and transit duration. We use the downsampled tables provided by
Hsu et al. (2019) to extract the window function data for each of the host
stars in our sample and obtain the values for periods in the range $P=1-800$
days, picking the transit duration such that it aligns with expectations of
circular orbits. We then average over the different stars in our sample.

Using the averaged window function data, we calculate the orbital periods at
which the window function is equal to 0.1, 0.3, 0.5, 0.7, and 0.9. The periods
are indicated with the horizontal lines in Figure 5. Since the perturbing
planets would have large radii (and thus large transit signal-to-noise
ratios), the main factor in their detection by the Kepler pipeline would
simply be whether they have three or more transits. Accordingly, the window
function effectively provides the probability that the perturbing planets
would be detected if they were on transiting orbits. We observe that roughly
half of the combined metastable region would have $p(\geq 3\
\mathrm{transits})=0.7$, indicating that many of them would be detected in
Kepler data if they existed and were transiting.

##  4 Discussion

###  4.1 Implications for the edge-of-the-multis

The results of Figure 5 indicate that perturbing planets with parameters
falling in the average metastable region would be readily detectable with RV
data (of a sufficiently long baseline) and transit data. What are the
implications of this? We argue that it indicates that distant giant planets
likely cannot explain the edge-of-the-multis, according to the following
logic. First, the metastable region delineates two regimes: (1) the unstable
regime in which the distant giant planet is strongly coupled to the inner
planets, and they cannot maintain dynamical stability, and (2) the stable
regime in which the distant giant planet’s perturbative influence is weak
relative to the coupling among the inner planets, and they maintain stability.
Accordingly, the metastable region (or, more accurately, the region just
slightly wide of it) is an approximate indicator of the distant giant “edge-
sculpting” properties, that is, the range of possible properties that a
perturbing planet could have if it is responsible for sculpting the outer edge
of the inner system. To see this, consider a system where the perturbing
planet has a period larger than the metastable location. With such properties,
the inner system could potentially have another planet added to its outer edge
and still be stable, so the giant planet would not be consistent with
sculpting the edge of the observed system.

If the metastable region thus approximates the parameters of edge-sculpting
perturbing planets, and if most of the metastable region would be readily
detectable with sufficient radial velocity and transit data, then edge-
sculpting perturbing planets would generally be detectable. However, such
perturbing planets are not widely detected in this region of parameter space.
Here we will briefly review the observational constraints.

As for transits, all but one of the compact multis in our sample are devoid of
a transiting giant planet with a period in the range $P=100-500$ days
confirmed by the Kepler search pipeline. One exception is the Kepler-90
(KOI-351) system, which contains a sun-like star and eight confirmed
transiting planets, the outermost of which has period $P=332$ days and mass
$203\pm 5\ M_{\oplus}$ (Liang et al., 2021). This outer planet may be
partially sculpting the edge of the inner seven planets. Another possible
exception is Kepler-87 (KOI-1574), which contains two planets beyond 100 days,
one planet at 115 days with mass $324\pm 9\ M_{\oplus}$ and one planet at 191
days with mass $6.4\pm 0.8\ M_{\oplus}$ (Ofir et al., 2014). The planet at 115
days could be responsible for dynamical perturbations, but it does not
strictly fit our definition of edge-sculpting giant planets, since the system
contains only two planets interior to it and a small exterior planet.

Beyond detections from the Kepler pipeline, Foreman-Mackey et al. (2016)
performed a systematic search for transiting long-period giant planets in
Kepler data, requiring only one or two transits (rather than the usual $\geq
3$). They found 16 planet candidates, one of which is in a system in our
sample. Kepler-154 (KOI-435) contains five inner super-Earths/sub-Neptunes
with periods ranging from $3-62$ days and an outer transiting giant with
$P=4.3^{+4.7}_{-1.3}$ yr and $R_{p}=0.83^{+0.12}_{-0.11}\ R_{\mathrm{Jup}}$
(Foreman-Mackey et al., 2016). However, this outer giant is too far from the
inner system to be responsible for sculpting the edge. Using injection-
recovery experiments, Foreman-Mackey et al. (2016) found that their search
method was $\sim 70\%$ complete to planets with periods $P\lesssim 3$ yr and
$R_{p}\gtrsim 0.5\ R_{\mathrm{Jup}}$. Roughly speaking, this indicates that
they would be $\sim 70\%$ complete to any transiting edge-sculpting giant
planets orbiting the stars in our sample, of which none were detected other
than in Kepler-154.

The lack of many transit detections is suggestive but inconclusive, since the
intrinsic transit probability is very small at long periods, and it depends on
the unknown orbital inclinations. However, we can still obtain a rough
estimate of the number of edge-sculpting giant planets we would expect to be
transiting in our sample. We use our calculated metastable region for each
system222For systems without SPOCK calculations (see Section 2.2), we randomly
pick a metastable region from the rest of the sample. and randomly sample an
orbital period for a giant planet within it. We then sample its inclination as
$i\sim\mathrm{Uniform}[80^{\circ},100^{\circ}]$. This assumes that the inner
system is approximately edge-on ($i\approx 90^{\circ}$), and it roughly agrees
with observed constraints on the mutual inclinations between inner super-
Earths and outer cold Jupiters. Masuda et al. (2020) found that these mutual
inclinations tend to be small, especially for inner systems with higher
transit multiplicity. Using the sampled period and inclination, we determine
whether or not the planet would be transiting. We then sum up the number of
transiting planets across the full sample of 64 systems. Finally, we repeat
this process 1,000 times to create a distribution of the number of transiting
planets. This yields an expectation of $2.0\pm 1.3$ transiting planets across
the 64 systems. The percentages of cases with 0, 1, 2, 3, 4, and 5 transiting
planets are 18%, 32%, 25%, 16%, 7%, and 2%. This indicates that it would be
fairly likely to observe zero transiting giant planets, even if all of the
systems hosted edge-sculpting giant planets. Constraints from transits are
thus not very powerful; RVs are likely to be more informative.

As for RV detections, several constraints have recently become available
through long term HARPS-N and Keck/HIRES surveys presented by Bonomo et al.
(2023) and Weiss et al. (2023), respectively. As for Bonomo et al. (2023),
their sample of 38 Kepler and K2 systems of small planets has only three
systems that overlap with our sample. None of those have new detected planets
or long-term trends. The Kepler Giant Planet Survey (Weiss et al., 2023)
observed 63 systems, 13 of which are common to our sample. Of those systems,
nine have no new detections or long-term trends, and the data places strong
constraints on the existence of additional long-period planets. In most cases,
the RV coverage places a $3\sigma$ upper limit of $M\sin i\lesssim 0.5\
M_{\mathrm{Jup}}$ at 5 AU (although there are variations from system to
system). The existence of edge-sculpting perturbers can thus be ruled out in
these systems.

The remaining four systems have RV-detected non-transiting companions:
Kepler-20 (KOI-70), Kepler-106 (KOI-116), Kepler-1130 (KOI-2169), and
Kepler-444 (KOI-3158). In Kepler-20 (KOI-70), there is a detected non-
transiting planet with $P=34.9$ days and $M\sin i=21.0\pm 3.4\ M_{\oplus}$ in
between the planets with $P=19.6$ days and $P=77.6$ days. This planet was
initially discovered by Buchhave et al. (2016), although it remains possible
that the signal is a false positive due to stellar activity (Nava et al.,
2020). Kepler-106 (KOI-116) has four small transiting planets and reveals
peaks in the RV residual periodogram near P = 90, 180, and 365 days. The
veracity of these signals is uncertain because they could be driven by the
seasonal window function. However, if there is a planet at 90 days, it would
have $M\sin i=46\pm 6\ M_{\oplus}$. Kepler-1130 (KOI-2169) contains a stellar-
mass companion with $P=14000\pm 1800$ days ($\sim 12$ AU), $M\sin i=218\pm 3\
M_{\mathrm{Jup}}$, and $e=0.65\pm 0.024$. The periapse distance is $\sim 4$
AU. Finally, Kepler-444 (KOI-3158) also contains an M-dwarf binary companion
(components B and C) first characterized by Dupuy et al. (2016). The orbital
solution in Weiss et al. (2023) is consistent with a recent analysis by Zhang
et al. (2023) and finds $P=87000\pm 3000$ days ($\sim 52$ AU), $M\sin i=629\pm
21\ M_{\mathrm{Jup}}$, and $e=0.55\pm 0.05$. The periapse distance is $\sim
23$ AU.

Among the systems just described, the Kepler-1130 system contains the only
perturber with properties that might be consistent with sculpting the outer
edge of the inner system, although the perturber is a star, not a giant planet
as considered in this work. It is worth noting that the two systems with
stellar-mass companions, Kepler-1130 (KOI-2169) and Kepler-444 (KOI-3158),
contain some of the smallest inner planets in our sample. The planets all have
smaller radii than Earth, and Mills & Fabrycky (2017) constrained the masses
of the Kepler-444 planets to be similar to that of Mars. It is thought that
the stellar companions limited the extent of the protoplanetary disk in which
the inner planets formed, thus limiting the solid material for planet
formation (e.g. Zhang et al., 2023). However, this is a different process than
the dynamical sculpting we have been considering in this paper.

Altogether, the existing constraints from both transits and radial velocities
suggest a lack of perturbing planets that are close enough to the inner
compact multis to sculpt their outer edges. The observational constraints are
incomplete, but they are sufficient to conclude that distant perturbers likely
contribute very little to the observed edge-of-the-multis. It is important to
clarify that the preceding discussion is not intended to be exhaustive, as it
neglected several complications. For instance, we neglected the possibility
that a perturbing planet that is not close enough to an inner system to render
it dynamically unstable could still potentially increase the mutual
inclinations among that inner system and reduce the transit multiplicity. This
would be a different form of edge-sculpting. Moreover, our estimation of the
expected number of transiting edge-sculpting giant planets used several
simplifying assumptions. It is beyond the scope of this work to perform a
thorough estimate of how many giant planets should have been detected in
existing transit and radial velocity observations if they are the primary
cause for the edge-of-the-multis. We leave a more thorough accounting for
these detection statistics to future work.

###  4.2 Alternative theories for the edge-of-the-multis

If distant giant planets are not the explanation of the edge-of-the-multis, we
must turn towards alternative hypotheses. Several theories were outlined in
Millholland et al. (2022), and additional theories have since been proposed.
We briefly recap them here. First, the formation process of compact multis may
naturally confine a system of $\sim 4-6$ inner planets to the $\lesssim 0.5-1$
AU region, with a truncation beyond them. For instance, the theory of “inside-
out planet formation” (Chatterjee & Tan, 2014, 2015) posits that planets may
form sequentially from successive gravitational instabilities of rings of
pebbles. The rings form via pebble drift due to gas drag and build up at the
pressure maximum associated with the dead-zone boundary. After the planets
form from the pebble ring, they migrate inwards, the dead-zone boundary and
pressure maximum retreat, and the process repeats. A similar theory of
successive formation of super-Earths from narrow rings of solid material was
recently proposed by Batygin & Morbidelli (2023). In this framework,
planetesimals form rapidly at the silicate sublimation line and grow through
pairwise collisions, until they achieve terminal proto-planet masses regulated
by isolation and migrate inwards. Systems that host distant giant planets have
also been shown to potentially drive the formation of high-density rings of
planetesimals due to the sweeping of secular resonances that transport
planetesimals into the inner disk (Best et al., 2023). All of these theories
naturally predict an outer edge of compact multis near the observed edge.

Even if the planets in compact multis do not form at well-defined discrete
locations in the inner disk, an edge-of-the-multis is also predicted by
theories of formation at larger orbital distances ($\gtrsim 1$ AU) followed by
inward migration. Specifically, Zawadzki et al. (2022) showed that migration
traps (generated by regions where the outward corotation torque dominates over
the inward Lindblad torque) yield planetary systems that are bifurcated into
two groups at small and larger periods with a gap at $\sim 100-300$ days.

To summarize, these theories predict that the formation of compact multis is a
balance between accretion and orbital migration that results in a compact
system of planets confined to the inner $\sim 0.5-1$ AU, with a truncation
just outside. We reviewed these hypotheses here for the sake of promoting
additional and probable arguments based on our findings, but further
investigation is beyond the scope of this paper. Future work on both planet
formation modeling and observational characterization of the $\gtrsim 1$ AU
regions of these systems should help determine which, if any, of these
mechanisms are dominant in sculpting the edge-of-the-multis.

##  5 Conclusion

This work was motivated by the recent finding that the outer edges of compact
multi-planet systems appear to truncate at smaller periods than expected from
geometric and detection biases alone (Millholland et al., 2022). The “edge-of-
the-multis” suggests the existence of some truncation (i.e. occurrence rate
fall-off) or transition (i.e. to smaller and/or more widely-spaced planets) in
the outer regions ($P\gtrsim 100$ days) of compact multis. Here we
investigated the question of whether the edge-of-the-multis could be caused by
the dynamical influence of distant giant planets, which are thought to be
common in systems with inner super-Earths/sub-Neptunes (Zhu & Wu, 2018; Bryan
et al., 2019).

We considered a sample of Kepler compact multi-planet systems with four or
more observed planets (Figure 1). We explored the dynamical stability of the
observed systems in the presence of hypothetical exterior giant planets. The
SPOCK machine learning stability classifier (Tamayo et al., 2020) was used to
compute the probability of dynamical stability of a given orbital
configuration. We randomly sampled the perturber period and mass in the range
$P_{p}\lesssim 1000$ days and $M_{p,p}\sim 0.5-5\ M_{\mathrm{Jup}}$. We tested
multiple sampling schemes for other parameters (e.g. the perturber
inclination) and found minimal sensitivity to these choices (Figures 2 and 3).

Given a set of perturber configurations for each system, we identified the
“metastable region”, the region in perturber mass/semi-major axis space that
divides the stable and unstable regimes (e.g. Figure 2). The metastable region
associated with an observed inner system approximately defines the parameter
regime in which a perturber would dynamically sculpt the edge of the system.
That is, a perturber with a wider orbit than the metastable region could leave
room for additional (undetected) inner planets and thus would not be an “edge-
sculpting” perturber. We found the metastable region to be in the range $P\sim
100-500$ days, with the strongest concentration around $P\sim 200$ days for
$M_{p,p}\sim 0.5\ M_{\mathrm{Jup}}$ and $P\sim 400$ days for $M_{p,p}=5\
M_{\mathrm{Jup}}$.

We explored the detectability of perturbing planets with “edge-sculpting”
properties, finding that they would generally be fairly easy to detect. If
such perturbers were transiting in Kepler data, many of them would transit at
least three times. Roughly half would have a $>0.7$ probability of three or
more transits. However, transit probabilities are small for long-period
planets, and they depend on the unknown orbital inclinations. A thorough model
of predicted transit yields would be necessary for a complete characterization
of the transit detectability. As for RV detection, the vast majority of
parameter space of edge-sculpting perturbers would have RV semi-amplitudes
well above 50 m/s.

Despite their apparent detectability, “edge-sculpting” perturbing planets have
not been identified in many compact multi systems. We reviewed the current
observational constraints, finding one system (Kepler-90) with a transiting
outer giant planet that may be consistent with edge sculpting. As for radial
velocity constraints, about 14 (perhaps more) of the 64 systems in our sample
have sufficient RV coverage to rule out perturbing planets with edge-sculpting
properties. The observational evidence is incomplete, but it clearly does not
point towards a large number of edge-sculpting perturbers. Future observations
will strengthen this result.

Taken as a whole, our results indicate that distant giant planets are unlikely
to play a significant role in sculpting the outer edges of compact multis. Our
finding narrows the search for the true cause of this phenomenon. The most
likely hypothesis is that the edge-of-the-multis is a signature of the
formation process, which may be relatively insensitive to the presence or
absence of distant perturbing planets. Future theoretical and observational
efforts may test this interpretation. In particular, we expect the upcoming
PLATO mission (Rauer et al., 2014) will yield further constraints on the outer
architectures of compact multis, especially if it revisits the Kepler field
for a long baseline.

##  6 Acknowledgements

We thank Dan Tamayo and the anonymous reviewer for their helpful comments and
suggestions. We thank the MIT Undergraduate Research Opportunities Program
(UROP) for their support in making this research possible. This research has
made use of the NASA Exoplanet Archive, which is operated by the California
Institute of Technology, under contract with the National Aeronautics and
Space Administration under the Exoplanet Exploration Program. Some of the
calculations presented in this paper were performed on the MIT Engaging
cluster at the Massachusetts Green High Performance Computing Center (MGHPCC)
facility. We gratefully MGHPCC for their computational resources and support.

##  Appendix A Testing the validity of the SPOCK calculations

Here we test the validity of the SPOCK calculations for our systems of inner
super-Earths and outer giant planets by verifying that the metastable region
agrees with the transition between stable and unstable regions identified by
alternative metrics. We first use the Mean Exponential Growth factor of Nearby
Orbits (MEGNO) chaos indicator, which quantifies the degree of divergence of
initially closely-separated trajectories in phase space (Cincotta & Simó,
2000; Cincotta et al., 2003). The time-averaged MEGNO value distinguishes
between stable/regular and unstable/chaotic orbital motion and thus provides
an independent calculation of the stability regions. We considered one of our
systems (KOI-812) as a case study and generated a similar $a_{p}$ vs.
$M_{p,p}$ map with 1,000 unique parameter combinations. Overall, we find good
agreement between the SPOCK and MEGNO maps, which increases confidence in the
validity of our SPOCK calculations.

We also check our results using an analytic stability criterion. Various
criteria have been developed for compact multis with distant perturbers (e.g.
Pu & Lai, 2018; Denham et al., 2019; Tamayo et al., 2021). We use the results
from Pu & Lai (2018) and compute the averaged coupling parameter,
$\bar{\epsilon}$ (their equation 60), which quantifies the extent to which the
inner planets are coupled to the giant planet ($\bar{\epsilon}\gg 1$) versus
each other ($\bar{\epsilon}\ll 1$). The coupling parameter is an approximate
stability criterion, since the degree of eccentricity and inclination
excitation in the inner system is much higher when they are tightly-coupled to
the giant planet. We compute $\bar{\epsilon}$ for the KOI-812 system and find
that the $\bar{\epsilon}=1$ transition line agrees well with the SPOCK
results. Given these numerical and analytic checks, we are confident in using
SPOCK for our stability calculations.

## References

  * Batygin & Morbidelli (2023) Batygin, K., & Morbidelli, A. 2023, Nature Astronomy, 7, 330, doi: 10.1038/s41550-022-01850-5
  * Berger et al. (2020a) Berger, T. A., Huber, D., Gaidos, E., van Saders, J. L., & Weiss, L. M. 2020a, AJ, 160, 108, doi: 10.3847/1538-3881/aba18a
  * Berger et al. (2020b) Berger, T. A., Huber, D., van Saders, J. L., et al. 2020b, AJ, 159, 280, doi: 10.3847/1538-3881/159/6/280
  * Best et al. (2023) Best, S., Sefilian, A. A., & Petrovich, C. 2023, arXiv e-prints, arXiv:2304.02045, doi: 10.48550/arXiv.2304.02045
  * Bonomo et al. (2023) Bonomo, A. S., Dumusque, X., Massa, A., et al. 2023, arXiv e-prints, arXiv:2304.05773, doi: 10.48550/arXiv.2304.05773
  * Borucki et al. (2010) Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977, doi: 10.1126/science.1185402
  * Bryan et al. (2019) Bryan, M. L., Knutson, H. A., Lee, E. J., et al. 2019, AJ, 157, 52, doi: 10.3847/1538-3881/aaf57f
  * Buchhave et al. (2016) Buchhave, L. A., Dressing, C. D., Dumusque, X., et al. 2016, AJ, 152, 160, doi: 10.3847/0004-6256/152/6/160
  * Burke & Catanzarite (2017) Burke, C. J., & Catanzarite, J. 2017, Planet Detection Metrics: Window and One-Sigma Depth Functions for Data Release 25, Kepler Science Document KSCI-19101-002, id. 14. Edited by Michael R. Haas and Natalie M. Batalha
  * Chatterjee & Tan (2014) Chatterjee, S., & Tan, J. C. 2014, ApJ, 780, 53, doi: 10.1088/0004-637X/780/1/53
  * Chatterjee & Tan (2015) —. 2015, ApJ, 798, L32, doi: 10.1088/2041-8205/798/2/L32
  * Chen & Kipping (2017) Chen, J., & Kipping, D. 2017, ApJ, 834, 17, doi: 10.3847/1538-4357/834/1/17
  * Cincotta et al. (2003) Cincotta, P. M., Giordano, C. M., & Simó, C. 2003, Physica D Nonlinear Phenomena, 182, 151, doi: 10.1016/S0167-2789(03)00103-9
  * Cincotta & Simó (2000) Cincotta, P. M., & Simó, C. 2000, A&AS, 147, 205, doi: 10.1051/aas:2000108
  * Dawson & Murray-Clay (2013) Dawson, R. I., & Murray-Clay, R. A. 2013, ApJ, 767, L24, doi: 10.1088/2041-8205/767/2/L24
  * Denham et al. (2019) Denham, P., Naoz, S., Hoang, B.-M., Stephan, A. P., & Farr, W. M. 2019, MNRAS, 482, 4146, doi: 10.1093/mnras/sty2830
  * Dupuy et al. (2016) Dupuy, T. J., Kratter, K. M., Kraus, A. L., et al. 2016, ApJ, 817, 80, doi: 10.3847/0004-637X/817/1/80
  * Fabrycky et al. (2014) Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ, 790, 146, doi: 10.1088/0004-637X/790/2/146
  * Fang & Margot (2012) Fang, J., & Margot, J.-L. 2012, ApJ, 761, 92, doi: 10.1088/0004-637X/761/2/92
  * Foreman-Mackey et al. (2016) Foreman-Mackey, D., Morton, T. D., Hogg, D. W., Agol, E., & Schölkopf, B. 2016, AJ, 152, 206, doi: 10.3847/0004-6256/152/6/206
  * Furlan et al. (2017) Furlan, E., Ciardi, D. R., Everett, M. E., et al. 2017, AJ, 153, 71, doi: 10.3847/1538-3881/153/2/71
  * Gratia & Fabrycky (2017) Gratia, P., & Fabrycky, D. 2017, MNRAS, 464, 1709, doi: 10.1093/mnras/stw2180
  * Hansen (2017) Hansen, B. M. S. 2017, MNRAS, 467, 1531, doi: 10.1093/mnras/stx182
  * He et al. (2019) He, M. Y., Ford, E. B., & Ragozzine, D. 2019, MNRAS, 490, 4575, doi: 10.1093/mnras/stz2869
  * Hsu et al. (2019) Hsu, D. C., Ford, E. B., Ragozzine, D., & Ashby, K. 2019, AJ, 158, 109, doi: 10.3847/1538-3881/ab31ab
  * Jones et al. (2006) Jones, H. R. A., Butler, R. P., Tinney, C. G., et al. 2006, MNRAS, 369, 249, doi: 10.1111/j.1365-2966.2006.10298.x
  * Lai & Pu (2017) Lai, D., & Pu, B. 2017, AJ, 153, 42, doi: 10.3847/1538-3881/153/1/42
  * Liang et al. (2021) Liang, Y., Robnik, J., & Seljak, U. 2021, AJ, 161, 202, doi: 10.3847/1538-3881/abe6a7
  * Lissauer et al. (2011) Lissauer, J. J., Ragozzine, D., Fabrycky, D. C., et al. 2011, ApJS, 197, 8, doi: 10.1088/0067-0049/197/1/8
  * Lissauer et al. (2014) Lissauer, J. J., Marcy, G. W., Bryson, S. T., et al. 2014, ApJ, 784, 44, doi: 10.1088/0004-637X/784/1/44
  * Masuda et al. (2020) Masuda, K., Winn, J. N., & Kawahara, H. 2020, AJ, 159, 38, doi: 10.3847/1538-3881/ab5c1d
  * Mayor et al. (2011) Mayor, M., Marmier, M., Lovis, C., et al. 2011, arXiv e-prints, arXiv:1109.2497.  https://arxiv.org/abs/1109.2497
  * Millholland et al. (2017) Millholland, S., Wang, S., & Laughlin, G. 2017, ApJ, 849, L33, doi: 10.3847/2041-8213/aa9714
  * Millholland et al. (2022) Millholland, S. C., He, M. Y., & Zink, J. K. 2022, AJ, 164, 72, doi: 10.3847/1538-3881/ac7c67
  * Mills & Fabrycky (2017) Mills, S. M., & Fabrycky, D. C. 2017, ApJ, 838, L11, doi: 10.3847/2041-8213/aa6543
  * Mills et al. (2019) Mills, S. M., Howard, A. W., Petigura, E. A., et al. 2019, AJ, 157, 198, doi: 10.3847/1538-3881/ab1009
  * NASA Exoplanet Archive (2022) NASA Exoplanet Archive. 2022, Kepler Objects of Interest DR25, Version: 2022-02-21 15:06, NExScI-Caltech/IPAC, doi: 10.26133/NEA5
  * NASA Exoplanet Archive (2023) —. 2023, Planetary Systems Composite Parameters, Version: 2023-02-09 08:12, NExScI-Caltech/IPAC, doi: 10.26133/NEA13
  * Nava et al. (2020) Nava, C., López-Morales, M., Haywood, R. D., & Giles, H. A. C. 2020, AJ, 159, 23, doi: 10.3847/1538-3881/ab53ec
  * Ofir et al. (2014) Ofir, A., Dreizler, S., Zechmeister, M., & Husser, T.-O. 2014, A&A, 561, A103, doi: 10.1051/0004-6361/201220935
  * Petigura et al. (2013) Petigura, E. A., Howard, A. W., & Marcy, G. W. 2013, Proceedings of the National Academy of Science, 110, 19273, doi: 10.1073/pnas.1319909110
  * Pu & Lai (2018) Pu, B., & Lai, D. 2018, MNRAS, 478, 197, doi: 10.1093/mnras/sty1098
  * Rauer et al. (2014) Rauer, H., Catala, C., Aerts, C., et al. 2014, Experimental Astronomy, 38, 249, doi: 10.1007/s10686-014-9383-4
  * Rosenthal et al. (2021) Rosenthal, L. J., Fulton, B. J., Hirsch, L. A., et al. 2021, ApJS, 255, 8, doi: 10.3847/1538-4365/abe23c
  * Rosenthal et al. (2022) Rosenthal, L. J., Knutson, H. A., Chachan, Y., et al. 2022, ApJS, 262, 1, doi: 10.3847/1538-4365/ac7230
  * Rowe et al. (2014) Rowe, J. F., Bryson, S. T., Marcy, G. W., et al. 2014, ApJ, 784, 45, doi: 10.1088/0004-637X/784/1/45
  * Spalding & Millholland (2020) Spalding, C., & Millholland, S. C. 2020, AJ, 160, 105, doi: 10.3847/1538-3881/aba629
  * Tamayo et al. (2021) Tamayo, D., Murray, N., Tremaine, S., & Winn, J. 2021, AJ, 162, 220, doi: 10.3847/1538-3881/ac1c6a
  * Tamayo et al. (2020) Tamayo, D., Cranmer, M., Hadden, S., et al. 2020, Proceedings of the National Academy of Science, 117, 18194, doi: 10.1073/pnas.2001258117
  * Thompson et al. (2018) Thompson, S. E., Coughlin, J. L., Hoffman, K., et al. 2018, ApJS, 235, 38, doi: 10.3847/1538-4365/aab4f9
  * Udry & Santos (2007) Udry, S., & Santos, N. C. 2007, ARA&A, 45, 397, doi: 10.1146/annurev.astro.45.051806.110529
  * Van Eylen & Albrecht (2015) Van Eylen, V., & Albrecht, S. 2015, ApJ, 808, 126, doi: 10.1088/0004-637X/808/2/126
  * Van Eylen et al. (2019) Van Eylen, V., Albrecht, S., Huang, X., et al. 2019, AJ, 157, 61, doi: 10.3847/1538-3881/aaf22f
  * Weiss et al. (2022) Weiss, L. M., Millholland, S. C., Petigura, E. A., et al. 2022, arXiv e-prints, arXiv:2203.10076.  https://arxiv.org/abs/2203.10076
  * Weiss et al. (2018) Weiss, L. M., Marcy, G. W., Petigura, E. A., et al. 2018, AJ, 155, 48, doi: 10.3847/1538-3881/aa9ff6
  * Weiss et al. (2023) Weiss, L. M., Isaacson, H., Marcy, G. W., et al. 2023, arXiv e-prints, arXiv:2304.00071, doi: 10.48550/arXiv.2304.00071
  * Wright et al. (2009) Wright, J. T., Upadhyay, S., Marcy, G. W., et al. 2009, ApJ, 693, 1084, doi: 10.1088/0004-637X/693/2/1084
  * Xie et al. (2016) Xie, J.-W., Dong, S., Zhu, Z., et al. 2016, Proceedings of the National Academy of Science, 113, 11431, doi: 10.1073/pnas.1604692113
  * Zawadzki et al. (2022) Zawadzki, B., Carrera, D., & Ford, E. B. 2022, ApJ, 937, 53, doi: 10.3847/1538-4357/ac8b04
  * Zhang et al. (2023) Zhang, Z., Bowler, B. P., Dupuy, T. J., et al. 2023, AJ, 165, 73, doi: 10.3847/1538-3881/aca88c
  * Zhu (2023) Zhu, W. 2023, arXiv e-prints, arXiv:2306.16691, doi: 10.48550/arXiv.2306.16691
  * Zhu et al. (2018) Zhu, W., Petrovich, C., Wu, Y., Dong, S., & Xie, J. 2018, ApJ, 860, 101, doi: 10.3847/1538-4357/aac6d5
  * Zhu & Wu (2018) Zhu, W., & Wu, Y. 2018, AJ, 156, 92, doi: 10.3847/1538-3881/aad22a