2024 July 172024 PASAPASA

supernovae: general — stars: massive — binaries: close

#  Blue supergiants as a progenitor of intermediate-luminosity red transients

Takashi J. Moriya11affiliation: National Astronomical Observatory of Japan,
National Institutes of Natural Sciences, 2-21-1 Osawa, Mitaka, Tokyo 181-8588,
Japan 22affiliation: Graduate Institute for Advanced Studies, SOKENDAI, 2-21-1
Osawa, Mitaka, Tokyo 181-8588, Japan 33affiliation: School of Physics and
Astronomy, Monash University, Clayton, VIC 3800, Australia
<EMAIL_ADDRESS>Athira Menon44affiliation: Instituto de Astrofísica
de Canarias, Avenida Vía Láctea s/n, 38205 La Laguna, Tenerife, Spain
55affiliation: Universidad de La Laguna, Departamento de Astrofísica, Avenida
Astrofísico Francisco Sánchez s/n, 38206 La Laguna, Tenerife, Spain
<EMAIL_ADDRESS>
###### Abstract

The current perspective about the explosions of blue supergiants is that they
resemble SN 1987A. These so-called peculiar Type II supernovae, however, are
one of the rarest types of supernovae and may not hence be the fate of all
blue supergiants. In this work, we explore other explosion scenarios for blue
supergiants. We create synthetic light curves from the explosions of blue
supergiant models born from binary mergers, over a range of explosion energies
and ${}^{56}\mathrm{Ni}$ masses. We find that blue supergiant explosions may
also lead to intermediate-luminosity red transients. We thus identify two
categories of supernovae possible from blue supergiant explosions: those with
high ${}^{56}\mathrm{Ni}$ masses ($\gtrsim 0.01$ $\mathrm{M}_{\odot}$) result
in slow-rising, dome-shaped light curves like SN 1987A. Lower
${}^{56}\mathrm{Ni}$ masses result in low-luminosity, short-plateau light
curves resembling some intermediate-luminosity red transients and Type II
supernovae like SN 2008bp, which are possible from the explosions of compact
blue supergiants and not from the far more extended red supergiants. Our
results indicate that blue supergiant explosions are more diverse than SN
1987A-like events and may be hidden among different kinds of transients,
explaining the possible discrepancies between the expected fraction of blue
supergiants born from binary mergers and the observed fraction of SN
1987A-like supernovae.

cccccccccc Properties of BSG progenitor models at core-carbon depletion used
in this work, selected from Menon et al. (2024). $M_{\textrm{1,i}}^{(a)}$
$M_{\textrm{2,i}}^{(b)}$ $P_{\textrm{i}}^{(c)}$ $f_{\textrm{c}}^{(d)}$ log
$L_{\star}^{(e)}$ T${}_{\textrm{eff,$\star$}}^{(f)}$ $M_{\star}^{(g)}$
$M_{\star}^{\textrm{env}(h)}$ $M_{\star}^{\textrm{core}(i)}$ $R_{\star}^{(j)}$
( $\mathrm{M}_{\odot}$) ( $\mathrm{M}_{\odot}$) (d) (L⊙) (kK) (
$\mathrm{M}_{\odot}$) ( $\mathrm{M}_{\odot}$) ( $\mathrm{M}_{\odot}$) (
$\mathrm{R}_{\odot}$)
(a) Initial mass of the primary of the pre-merger binary.
(b) Initial mass of the secondary of the pre-merger binary.
(c) Initial orbital period of the pre-merger binary.
(d) Core dredge-up factor at the time of the merger as defined in Menon et al.
(2024).
(e) Luminosity of post-merger star at central-carbon depletion.
(f) Effective temperature of post-merger star at central-carbon depletion.
(g) Stellar mass of post-merger star at central-carbon depletion.
(h) Hydrogen-rich envelope mass of post-merger star at central-carbon
depletion.
(i) Hydrogen-free core mass of post-merger star at central-carbon depletion.
(j) Stellar radius of post-merger star at central-carbon depletion.
11.2 9.0 1413 0.06 5.0 30.6 16.9 13.2 3.7 10.8
20.0 15.0 3162 0.08 5.6 33.8 32.3 25.7 6.6 17.8
20.0 2.0 16 0.00 5.2 18.8 19.7 12.6 7.1 37.4
31.6 13.3 28 0.28 5.5 20.3 28.6 17.0 11.6 50.8
25.1 1.3 1000 0.00 5.4 15.0 24.6 14.7 9.9 76.0

##  1 Introduction

Blue supergiants (BSGs), with spectral classes between B0 and BI, are observed
in plenty in the local and distant Universe (Menon et al., 2024). Yet, despite
their abundance, the nature of their explosive fates has not been thoroughly
explored. SN 1987A was the explosion of a BSG, whose progenitor was imaged to
be a BIa supergiant (Sonneborn et al., 1987). Indeed, the progenitor
properties of SN 1987A can be well explained via stellar-merger scenario
(Podsiadlowski et al., 1992; Menon & Heger, 2017; Urushibata et al., 2018),
including the features of the supernova (SN) light curve as well, along with a
few other SNe with similar light-curve shapes (Menon et al., 2019).
Collectively these SN 1987A-like SNe are called “peculiar Type II SNe.”

Recently, Menon et al. (2024) demonstrated that the majority of BSGs observed
today have been born from binary mergers of post-main sequence Hertzsprung-gap
or red supergiants (RSGs) with their main-sequence companions. They further
predicted that a substantial fraction of the observed BSGs may also die as
BSGs. Stellar mergers are expected to be the outcome of about $\sim 30\%$ of
massive binary star systems (de Mink et al., 2014), of which the majority may
occur during the post-main sequence expansion of the primary star. SN
1987A-like events, however, are rare, comprising only $\lesssim 5$% of the
known core-collapse SNe (e.g., Kleiser et al. (2011); Pastorello et al.
(2012); Taddia et al. (2016); Sit et al. (2023)). This fraction is far fewer
than the observed fraction of BSGs among massive star populations, a large
number of which are expected to remain as BSGs when they undergo core-collapse
and explode, according to the binary-merger models of Menon et al. (2024).

The light curve of SN 1987A has a slow rise lasting for about 80 days (Arnett
et al., 1989), while typical Type II SNe have a luminosity plateau lasting for
around 100 days (Anderson et al., 2014). The explosions of BSGs have been
mostly associated with the characteristic slow-rising, dome-shaped light
curves of SN 1987A which are due to the compactness of BSGs
(R${}_{\star}\lesssim 70$ $\mathrm{R}_{\odot}$, e.g., Taddia et al. (2016);
Sit et al. (2023) but see also Kleiser et al. (2011); Tsuna et al. (2020)).
Due to the small radii of BSGs, the adiabatic expansion after the explosion is
efficient and, therefore, the luminosity declines quickly at first. Then, the
SN ejecta are heated by the nuclear decay energy of ${}^{56}\mathrm{Ni}$ and
this ${}^{56}\mathrm{Ni}$ heating makes the characteristic slow rise of SN
1987A (e.g., Woosley (1988); Shigeyama & Nomoto (1990); Blinnikov et al.
(2000); Dessart & Hillier (2019); Pumo et al. (2023)). In the case of SN
1987A, the amount of ${}^{56}\mathrm{Ni}$ synthesized at the explosion is
estimated to be around 0.07 $\mathrm{M}_{\odot}$ (Arnett et al., 1989).

Since the slow light-curve rise observed in SN 1987A-like SNe is powered by
${}^{56}\mathrm{Ni}$ decay heating, this characteristic slow light-curve rise
would become less significant if the amount of ${}^{56}\mathrm{Ni}$
synthesized during the explosions of BSGs is small. In such a case, BSG
explosions might not be observed as SN 1987A-like events and it is possible
that a population of BSG explosions do not become SN 1987A-like events.

Intermediate-luminosity red transients (ILRTs) are a class of “gap transients”
which have peak absolute magnitudes between those of classical novae and SNe
(Rau et al. (2009); Kasliwal (2012); see Pastorello & Fraser (2019) for a
review). They are characterized by a low-luminosity plateau lasting for a
couple of months (e.g., Cai et al. (2021)). The late-time light-curve decline
rate is consistent with those of the ${}^{56}\mathrm{Co}$ decay rate from
$\mathrm{M}_{{}^{56}\mathrm{Ni}}\sim 10^{-3}-10^{-4}~{}\,\mathrm{M}_{\odot}$,
indicating that they could be related to terminal SN explosions (e.g., Cai et
al. (2021)). Because of their faint nature as well as the low-mass progenitor
observed for SN 2008S (Prieto et al., 2008; Adams et al., 2016), ILRTs have
been related to low mass core-collapse SNe such as electron-capture SNe (e.g.,
Botticella et al. (2009); Cai et al. (2021), but see also Tominaga et al.
(2013); Moriya et al. (2014); Kozyreva et al. (2021)). In the case of a
potential ILRT AT 2019krl, a BSG was identified as a progenitor (Andrews et
al., 2021). Thus IRLTs can come from several different progenitors.

In this work, we explore the light-curve properties of BSG explosions over a
range of explosion energies ($E=0.001-1$ B, where $1~{}\mathrm{B}\equiv
10^{51}~{}\mathrm{erg}$) and ${}^{56}\mathrm{Ni}$ masses
($\mathrm{M}_{{}^{56}\mathrm{Ni}}=0.001-0.07\,\mathrm{M}_{\odot}$). The
progenitor models are taken from Menon et al. (2024) which are evolved from
detailed stellar-merger models until core carbon-depletion and have masses of
$16-40$ $\mathrm{M}_{\odot}$. We identify a dichotomy in SNe from BSG
explosions, which depends mainly on the ${}^{56}\mathrm{Ni}$ mass
($\mathrm{M}_{{}^{56}\mathrm{Ni}}$) ejected: for
$\mathrm{M}_{{}^{56}\mathrm{Ni}}\gtrsim 0.01~{}\,\mathrm{M}_{\odot}$ we get SN
1987A-like light curves and for $\mathrm{M}_{{}^{56}\mathrm{Ni}}\lesssim
0.01~{}\,\mathrm{M}_{\odot}$ we get low-luminosity, short-plateau light curves
that are consistent with certain ILRTs and Type II SNe like SN 2008bp.

The rest of this paper is organized as follows. We first introduce the BSG
progenitor models and our methods of SN light-curve calculations in Section 2.
We present the results of the light-curve calculations in Section 3 and show
that BSG explosions can result in low-luminosity, short-plateau Type II SNe
when the amount of ${}^{56}\mathrm{Ni}$ is sufficiently small. We compare the
low-luminosity, short-plateau light curves from BSGs with those of ILRTs in
Section 4 and show that they are indeed similar. We summarize our discoveries
of this paper in Section 5.

##  2 Methods

###  2.1 Progenitor models

The BSG models used in this paper are from Menon et al. (2024), which have
been evolved through a binary merger of a core hydrogen-depleted giant (on the
Hertzsprung gap or as a RSG) with its main-sequence companion. These post-
merger models evolved to core He-burning BSGs and matched the spectroscopic
properties of the majority of the sample of early B-type supergiants in the
Large Magellanic Cloud (LMC). Those models with masses $\approx 16-25$ M⊙
always exploded as BSGs, while those with higher masses (upto 40 M⊙) either
exploded as BSGs or as cooler stars with T${}_{\textrm{eff}}\leq 10$ kK,
depending on the dynamics of the merging process.

For this study, we select representative BSG progenitor models that span the
upper and lower mass range of B-type supergiants. These models are evolved
until core carbon-depletion and summarized in Table Blue supergiants as a
progenitor of intermediate-luminosity red transients. After the core carbon
depletion, the envelope structure determining light-curve properties is not
expected to change much until core collapse. The density structure of the BSG
models are presented in Fig. 1. In the same figure, we also show the 15.8
$\mathrm{M}_{\odot}$ (20.0 $\mathrm{M}_{\odot}$ at the zero-age main sequence
[ZAMS]) RSG SN progenitor model (1064 $\mathrm{R}_{\odot}$) from Sukhbold et
al. (2016) for comparison.

Figure 1:  Density structure of BSG progenitors adopted in this work. The top
panel shows the structure in the radius coordinate and the bottom panel shows
the structure in mass coordinate. The structure of a RSG progenitor model from
Sukhbold et al. (2016) is presented with the dashed line for comparison.

###  2.2 Explosion and light-curve modeling

We conduct explosion and subsequent light-curve modeling by using one-
dimensional radiation hydrodynamics code STELLA (Blinnikov et al., 1998, 2000,
2006). The central 1.4 $\mathrm{M}_{\odot}$ of the progenitor models is
removed as a compact remnant and the remaining structure is put as the initial
condition for the STELLA calculations. The explosion is triggered by inserting
thermal energy at the innermost layers of the progenitors above the mass cut
at 1.4 $\mathrm{M}_{\odot}$. STELLA does not treat explosive nucleosynthesis
and the amount of 56Ni is a free parameter. In this work, 56Ni is uniformly
mixed within the entire ejecta.

Evolution of spectral energy distributions (SEDs) is numerically followed in
each timestep in STELLA. We adopt 100 frequency bins between 1 Å and 50,000 Å
in a log scale. We can obtain light curves in any photometic bands by
convolving the SEDs with any filter response functions. We present bolometric
light curves as well as the V and R-band light curves in the Johnson-Cousins
filter system (Bessell, 1990) in this paper.

Figure 2:  Effects of ${}^{56}\mathrm{Ni}$ mass in shaping bolometric light
curves of BSG explosions. The top panel shows the $1~{}\mathrm{B}$ models and
the bottom panel shows the $0.3~{}\mathrm{B}$ models. We fix the progenitor
model ($R=37.4~{}\mathrm{R_{\odot}}$ and $M=19.7~{}\mathrm{M_{\odot}}$) in
this figure. When the ${}^{56}\mathrm{Ni}$ mass is large enough, we find the
luminosity peak powered by the ${}^{56}\mathrm{Ni}$ heating becomes apparent
and the SN 1987A-like SN light curves can be reproduced. When the
${}^{56}\mathrm{Ni}$ mass becomes too low, we stop finding the light curve
peak powered by the ${}^{56}\mathrm{Ni}$ heating and BSG explosions end up to
have low-luminosity, short-plateau light curves. The bolometric light curve of
SN 1987A (Hamuy et al., 1988) is presented as a reference.

Figure 3:  Bolometric light curves from the explosions of BSGs with various
radii ($10.8-76.0~{}\,\mathrm{R}_{\odot}$) with
$\mathrm{M}_{{}^{56}\mathrm{Ni}}=0.03,0.01,0.003$ and
$0.001~{}\,\mathrm{M}_{\odot}$ and $E=0.3~{}\mathrm{B}$. A light curve of a
RSG explosion ($R=414~{}\,\mathrm{R}_{\odot}$, $M=8.7~{}\,\mathrm{M}_{\odot}$
with the ZAMS mass of 9.0 $\mathrm{M}_{\odot}$ in Sukhbold et al. (2016)) with
$\mathrm{M}_{{}^{56}\mathrm{Ni}}=0.001~{}\,\mathrm{M}_{\odot}$ is presented
for comparison (dashed line). In one case, the numerical computation was
terminated before reaching 120 days because of numerical issues.

##  3 Results

Fig. 2 presents how ${}^{56}\mathrm{Ni}$ mass affects the light curves of BSG
explosions based on a progenitor model for SN 1987A. This progenitor has the
radius $R=37.4~{}\,\mathrm{R}_{\odot}$ and the mass
$M=19.7~{}\,\mathrm{M}_{\odot}$. This mass does not exclude the compact
remnant mass of 1.4 $\mathrm{M}_{\odot}$ and the ejecta mass in the
calculation is 18.3 $\mathrm{M}_{\odot}$. In Fig. 2a, we assume the explosion
energy $E=1~{}\mathrm{B}$. We can find that the bolometric light curve of the
explosion of this progenitor model with
$\mathrm{M}_{{}^{56}\mathrm{Ni}}=0.07~{}\,\mathrm{M}_{\odot}$ matches well to
that of SN 1987A.

As we decrease $\mathrm{M}_{{}^{56}\mathrm{Ni}}$ from this model, we can find
that the light-curve peak caused by the ${}^{56}\mathrm{Ni}$ heating become
earlier and fainter. Then, when the ${}^{56}\mathrm{Ni}$ mass becomes
sufficiently low ($\mathrm{M}_{{}^{56}\mathrm{Ni}}\lesssim
0.003~{}\,\mathrm{M}_{\odot}$ in the example in Fig. 2a), we no longer find
the luminosity increase due to the ${}^{56}\mathrm{Ni}$ heating. The previous
study by Pumo et al. (2023) also found that the effects of the
${}^{56}\mathrm{Ni}$ heating become negligible with a similar
${}^{56}\mathrm{Ni}$ mass. When the effects of the ${}^{56}\mathrm{Ni}$
heating become negligible, the luminosity is mainly powered by the thermal
energy provided by the initial shock wave triggered the explosion as in the
case of Type II SNe from more extended RSGs. Because BSGs are compact, the
adiabatic cooling becomes efficient and the plateau luminosity and duration
become low. In the example presented in Fig. 2, we can find that the plateau
luminosity is around $10^{41}~{}\mathrm{erg}$ with the plateau duration of
around 50 days. Such a short, low-luminosity plateau is difficult to be
realized by RSG explosions because low explosion energy required to realize a
faint plateau luminosity makes a plateau duration even longer (e.g., Kasen &
Woosley (2009)).

Because $\mathrm{M}_{{}^{56}\mathrm{Ni}}$ is required to be small not to have
a luminosity peak powered by the ${}^{56}\mathrm{Ni}$ heating, it would be
natural to investigate BSG explosions with a small explosion energy. In Fig.
2b, we present the bolometric light curves of BSG explosions with
$E=0.3~{}\mathrm{B}$ with various ${}^{56}\mathrm{Ni}$ masses. With the lower
explosion energy, we still find that the effect of the ${}^{56}\mathrm{Ni}$
heating disappears with $\mathrm{M}_{{}^{56}\mathrm{Ni}}\lesssim
0.003~{}\,\mathrm{M}_{\odot}$. The pleateau luminosity decreases and the
plateau duration increases with the lower explosion energy in a similar way as
in Type II SNe from RSGs.

Fig. 3 explores the effects of the ${}^{56}\mathrm{Ni}$ mass on the light
curves of the explosions of BSGs with various radii. The adiabatic cooling
depends on the progenitor radii, and they are expected to affect the
properties of BSG explosions with low $\mathrm{M}_{{}^{56}\mathrm{Ni}}$. All
the BSG explosion models present the ${}^{56}\mathrm{Ni}$-powered slow rise in
their light curves when
$\mathrm{M}_{{}^{56}\mathrm{Ni}}=0.03~{}\,\mathrm{M}_{\odot}$. When
$\mathrm{M}_{{}^{56}\mathrm{Ni}}=0.01~{}\,\mathrm{M}_{\odot}$, the slow rise
from the ${}^{56}\mathrm{Ni}$ heating starts to be insignificant in the models
with $R=37.4~{}\,\mathrm{R}_{\odot}$, 51.0 $\mathrm{R}_{\odot}$, and 76.0
$\mathrm{R}_{\odot}$, and they only show the plateau phase when
$\mathrm{M}_{{}^{56}\mathrm{Ni}}=0.003~{}\,\mathrm{M}_{\odot}$. The
$R=17.7~{}\,\mathrm{R}_{\odot}$ BSG explosion shows the low-luminosity, short
plateau with insignificant luminosity increase with
$\mathrm{M}_{{}^{56}\mathrm{Ni}}=0.001~{}\,\mathrm{M}_{\odot}$. Our most
compact BSG model ($R=10.8~{}\,\mathrm{R}_{\odot}$) may still have a small
effect of the ${}^{56}\mathrm{Ni}$ heating with
$\mathrm{M}_{{}^{56}\mathrm{Ni}}\simeq 0.001~{}\,\mathrm{M}_{\odot}$. We can
also find that low-luminosity, short-plateau light curves are difficult to be
realized by RSG explosions. Overall, BSG explosions need to have
$\mathrm{M}_{{}^{56}\mathrm{Ni}}\lesssim 0.01~{}\,\mathrm{M}_{\odot}$ in order
not to have a luminosity peak so that they can be observed as low-luminosity,
short-plateau Type II SNe. This critical mass is smaller than the average
${}^{56}\mathrm{Ni}$ mass estimated in Type II SNe ($\simeq
0.04~{}\,\mathrm{M}_{\odot}$, Anderson (2019)).

Because the plateau phase is caused by the recombination within the hydrogen-
rich envelope as in the case of RSG explosions, the photospheric temperature
during the plateau phase is found to be around 5000 K. The photospheric
velocities can be small when the explosion energy is small as presented in
Fig. 4.

##  4 Comparison with observations

Among Type II plateau SNe, low-luminosity ($\sim
10^{41}~{}\mathrm{erg~{}s^{-1}}$), short-duration plateaus of $\sim
50~{}\mathrm{days}$ are quite rare. Still, in the Type II SN light-curve
sample compiled by Anderson et al. (2014), one Type II SN, SN 2008bp, is found
to have such a property and its light curve indeed match well to our BSG
explosion model with a small amount of ${}^{56}\mathrm{Ni}$ having no slow
rise powered by the ${}^{56}\mathrm{Ni}$ heating (the top panel of Fig. 5).

The more common variety of low-luminosity, short-plateau light curves are
found among the sample of ILRTs. In the middle panel of Fig. 5, we compare our
low-luminosity, short-plateau BSG explosion models with some representative
ILRTs. The low-luminosity plateau phase of ILRTs lasting for a couple of
months is consistent with the predicted properties of low-energy, small
${}^{56}\mathrm{Ni}$ mass explosions of BSGs discussed in the previous
section. Indeed, their overall light-curve features are found to be similar to
ILRTs. Some differences such as the early luminosity bump in, e.g., SNhunt120,
may be caused by circumstellar interaction, for example. The low photospheric
velocities observed in IRLTs (e.g., Cai et al. (2021)) are also consistent
with our BSG models (Fig. 4). Therefore, we propose that a fraction ILRTs may
originate from faint BSG explosions with a small amount of
${}^{56}\mathrm{Ni}$.

The event rate of ILRTs is estimated to be $1-10$% of core-collapse SNe (Cai
et al., 2021; Karambelkar et al., 2023). If a fraction of ILRTs as well as SN
1987A-like SNe originate from BSGs, up to around 10% of core-collapse SNe may
originate from BSG explosions. We note again that some ILRTs are known to
originate from low-mass progenitors (e.g., SN 2008S) and not all of ILRTs
originate from BSGs.

The second group of light curves consist of various SN 1987A-like SNe, that
result from high $\mathrm{M}_{{}^{56}\mathrm{Ni}}$ explosions of the BSG
models (Fig. 5, bottom panel). These peculiar Type-II SNe are similar in their
light-curve shapes to the prototype SN 1987A, and can be well-explained by
high-energy explosions ($E\approx 1$ B) with
$\mathrm{M}_{{}^{56}\mathrm{Ni}}\geq 0.07$ $\mathrm{M}_{\odot}$. A more
precise reproduction of the observed light curves (such as in Menon et al.
(2019)), can be done by exploring the explosions of other BSG progenitor
models.

Figure 4:  Photospheric velocity evolution of the BSG explosion models with
$\mathrm{M}_{{}^{56}\mathrm{Ni}}=0.001~{}\,\mathrm{M}_{\odot}$ and
$E=0.3~{}\mathrm{B}$. Photosphere is set where the Rosseland-mean optical
depth becomes $2/3$.

Figure 5:  Top: V-band light curves of the BSG explosion models presented in
Fig. 2a. We can find that the low-luminosity, short-plateau Type II SN 2008bp
reported in Anderson et al. (2014) can be explained by the BSG explosion with
$\mathrm{M}_{{}^{56}\mathrm{Ni}}=0.003~{}\,\mathrm{M}_{\odot}$. Middle: The
R-band light curves of some BSG explosions with
$\mathrm{M}_{{}^{56}\mathrm{Ni}}=0.001~{}\,\mathrm{M}_{\odot}$ are compared to
those of ILRTs. Their overall light-curve features are similar in many cases.
In one case, the numerical computation was terminated before reaching 120 days
because of numerical issues. The light curves of AT 2013lb (Cai et al., 2021),
AT 2018aes (Cai et al., 2021), M85 OT2006-1 (Kulkarni et al., 2007; Pastorello
et al., 2007), NGC 300-2008OT1 (Berger et al., 2009; Bond et al., 2009; Prieto
et al., 2009; Kashi et al., 2010; Humphreys et al., 2011), and SNhunt120
(Stritzinger et al., 2020) are presented. Bottom: The bolometric light curves
of BSG explosions with $\mathrm{E=1~{}B}$ and
$\mathrm{M}_{{}^{56}\mathrm{Ni}}=0.07~{}\,\mathrm{M}_{\odot}$ compared with
those of SN 1987A (Hamuy et al., 1988) and other SN 1987A-like SNe: SN 1998A
(Pastorello et al., 2005), SN 2000cb (Kleiser et al., 2011; Utrobin & Chugai,
2011), SN 2009E (Pastorello et al., 2012), and SN 2009mw (Takáts et al.,
2016).

##  5 Summary

We explored explosion properties of BSGs especially when the amount of
${}^{56}\mathrm{Ni}$ produced during the explosions is small. When
$\mathrm{M}_{{}^{56}\mathrm{Ni}}\lesssim 0.01~{}\,\mathrm{M}_{\odot}$, the
characteristic slow rise in light curves observed in SN 1987A-like SN light
curves may disappear, and BSG explosions are observed as low-luminosity
(around $10^{40}-10^{41}~{}\mathrm{ergs^{-1}}$), short-plateau ($50-90$ days)
SNe. This is because of the compactness of the BSG progenitors that allows the
adiabatic cooling efficient. The exact threshold ${}^{56}\mathrm{Ni}$ mass for
BSG explosions to be observed as such Type II SNe depends on the radii of the
BSG progenitors. Low-energy, low-${}^{56}\mathrm{Ni}$ mass explosions of BSGs
can result in transients similar to those of some ILRTs. Therefore, we
conclude that a fraction of ILRTs may originate from BSG explosions. The
upcoming Legacy Survey of Space and Time (LSST) by the Vera C. Rubin
Observatory would allow us to explore the nature of these faint transients and
they may be able to probe diverse deaths of BSGs.

We thank Yongzhi Cai for sharing the photometric data of ILRTs, along with
Norbert Langer and Alexander Heger for vital discussions that have shaped this
paper. TJM is supported by the Grants-in-Aid for Scientific Research of the
Japan Society for the Promotion of Science (JP24K00682, JP24H01824,
JP21H04997, JP24H00002, JP24H00027, JP24K00668) and by the Australian Research
Council (ARC) through the ARC’s Discovery Projects funding scheme (project
DP240101786). AM acknowledges the support of the Juan de la Cierva-
Incorporación in conducting this research. AM also thanks the NAOJ Visiting
Joint Research program for hosting her as an invited researcher. Numerical
computations were in part carried out on PC cluster at Center for
Computational Astrophysics, National Astronomical Observatory of Japan. This
work was supported by the NAOJ Research Coordination Committee, NINS (NAOJ-
RCC-2401-0402).

## References

  * Adams et al. (2016) Adams, S. M., Kochanek, C. S., Prieto, J. L., et al. 2016, MNRAS, 460, 1645
  * Anderson (2019) Anderson, J. P. 2019, A&A, 628, A7
  * Anderson et al. (2014) Anderson, J. P., González-Gaitán, S., Hamuy, M., et al. 2014, ApJ, 786, 67
  * Andrews et al. (2021) Andrews, J. E., Jencson, J. E., Van Dyk, S. D., et al. 2021, ApJ, 917, 63
  * Arnett et al. (1989) Arnett, W. D., Bahcall, J. N., Kirshner, R. P., & Woosley, S. E. 1989, ARA&A, 27, 629
  * Berger et al. (2009) Berger, E., Soderberg, A. M., Chevalier, R. A., et al. 2009, ApJ, 699, 1850
  * Bessell (1990) Bessell, M. S. 1990, PASP, 102, 1181
  * Blinnikov et al. (2000) Blinnikov, S., Lundqvist, P., Bartunov, O., Nomoto, K., & Iwamoto, K. 2000, ApJ, 532, 1132
  * Blinnikov et al. (1998) Blinnikov, S. I., Eastman, R., Bartunov, O. S., Popolitov, V. A., & Woosley, S. E. 1998, ApJ, 496, 454
  * Blinnikov et al. (2006) Blinnikov, S. I., Röpke, F. K., Sorokina, E. I., et al. 2006, A&A, 453, 229
  * Bond et al. (2009) Bond, H. E., Bedin, L. R., Bonanos, A. Z., et al. 2009, ApJ, 695, L154
  * Botticella et al. (2009) Botticella, M. T., Pastorello, A., Smartt, S. J., et al. 2009, MNRAS, 398, 1041
  * Cai et al. (2021) Cai, Y. Z., Pastorello, A., Fraser, M., et al. 2021, A&A, 654, A157
  * de Mink et al. (2014) de Mink, S. E., Sana, H., Langer, N., Izzard, R. G., & Schneider, F. R. N. 2014, ApJ, 782, 7
  * Dessart & Hillier (2019) Dessart, L., & Hillier, D. J. 2019, A&A, 622, A70
  * Hamuy et al. (1988) Hamuy, M., Suntzeff, N. B., Gonzalez, R., & Martin, G. 1988, AJ, 95, 63
  * Humphreys et al. (2011) Humphreys, R. M., Bond, H. E., Bedin, L. R., et al. 2011, ApJ, 743, 118
  * Karambelkar et al. (2023) Karambelkar, V. R., Kasliwal, M. M., Blagorodnova, N., et al. 2023, ApJ, 948, 137
  * Kasen & Woosley (2009) Kasen, D., & Woosley, S. E. 2009, ApJ, 703, 2205
  * Kashi et al. (2010) Kashi, A., Frankowski, A., & Soker, N. 2010, ApJ, 709, L11
  * Kasliwal (2012) Kasliwal, M. M. 2012, PASA, 29, 482
  * Kleiser et al. (2011) Kleiser, I. K. W., Poznanski, D., Kasen, D., et al. 2011, MNRAS, 415, 372
  * Kozyreva et al. (2021) Kozyreva, A., Baklanov, P., Jones, S., Stockinger, G., & Janka, H.-T. 2021, MNRAS, 503, 797
  * Kulkarni et al. (2007) Kulkarni, S. R., Ofek, E. O., Rau, A., et al. 2007, Nature, 447, 458
  * Menon & Heger (2017) Menon, A., & Heger, A. 2017, MNRAS, 469, 4649
  * Menon et al. (2019) Menon, A., Utrobin, V., & Heger, A. 2019, MNRAS, 482, 438
  * Menon et al. (2024) Menon, A., Ercolino, A., Urbaneja, M. A., et al. 2024, ApJ, 963, L42
  * Moriya et al. (2014) Moriya, T. J., Tominaga, N., Langer, N., et al. 2014, A&A, 569, A57
  * Pastorello & Fraser (2019) Pastorello, A., & Fraser, M. 2019, Nature Astronomy, 3, 676
  * Pastorello et al. (2005) Pastorello, A., Baron, E., Branch, D., et al. 2005, MNRAS, 360, 950
  * Pastorello et al. (2007) Pastorello, A., Della Valle, M., Smartt, S. J., et al. 2007, Nature, 449, 1
  * Pastorello et al. (2012) Pastorello, A., Pumo, M. L., Navasardyan, H., et al. 2012, A&A, 537, A141
  * Podsiadlowski et al. (1992) Podsiadlowski, P., Joss, P. C., & Hsu, J. J. L. 1992, ApJ, 391, 246
  * Prieto et al. (2009) Prieto, J. L., Sellgren, K., Thompson, T. A., & Kochanek, C. S. 2009, ApJ, 705, 1425
  * Prieto et al. (2008) Prieto, J. L., Kistler, M. D., Thompson, T. A., et al. 2008, ApJ, 681, L9
  * Pumo et al. (2023) Pumo, M. L., Cosentino, S. P., Pastorello, A., et al. 2023, MNRAS, 521, 4801
  * Rau et al. (2009) Rau, A., Kulkarni, S. R., Law, N. M., et al. 2009, PASP, 121, 1334
  * Shigeyama & Nomoto (1990) Shigeyama, T., & Nomoto, K. 1990, ApJ, 360, 242
  * Sit et al. (2023) Sit, T., Kasliwal, M. M., Tzanidakis, A., et al. 2023, ApJ, 959, 142
  * Sonneborn et al. (1987) Sonneborn, G., Altner, B., & Kirshner, R. P. 1987, ApJ, 323, L35
  * Stritzinger et al. (2020) Stritzinger, M. D., Taddia, F., Fraser, M., et al. 2020, A&A, 639, A103
  * Sukhbold et al. (2016) Sukhbold, T., Ertl, T., Woosley, S. E., Brown, J. M., & Janka, H. T. 2016, ApJ, 821, 38
  * Taddia et al. (2016) Taddia, F., Sollerman, J., Fremling, C., et al. 2016, A&A, 588, A5
  * Takáts et al. (2016) Takáts, K., Pignata, G., Bersten, M., et al. 2016, MNRAS, 460, 3447
  * Tominaga et al. (2013) Tominaga, N., Blinnikov, S. I., & Nomoto, K. 2013, ApJ, 771, L12
  * Tsuna et al. (2020) Tsuna, D., Ishii, A., Kuriyama, N., Kashiyama, K., & Shigeyama, T. 2020, ApJ, 897, L44
  * Urushibata et al. (2018) Urushibata, T., Takahashi, K., Umeda, H., & Yoshida, T. 2018, MNRAS, 473, L101
  * Utrobin & Chugai (2011) Utrobin, V. P., & Chugai, N. N. 2011, A&A, 532, A100
  * Woosley (1988) Woosley, S. E. 1988, ApJ, 330, 218


# Confirming Resonance in Three Transiting Systems

Tyler Quinn Department of Astronomy & Astrophysics, Center for Exoplanets and
Habitable Worlds, The Pennsylvania State University, University Park, PA
16802, USA  Mariah G. MacDonald Department of Astronomy & Astrophysics, Center
for Exoplanets and Habitable Worlds, The Pennsylvania State University,
University Park, PA 16802, USA  Department of Physics, The College of New
Jersey, 2000 Pennington Road, Ewing, NJ 08628, USA

###### Abstract

Although resonant planets have orbital periods near commensurability,
resonance is also dictated by other factors, such as the planets’
eccentricities and masses, and therefore must be confirmed through a study of
the system’s dynamics. Here, we perform such a study for five multi-planet
systems: Kepler-226, Kepler-254, Kepler-363, Kepler-1542, and K2-32. For each
system, we run a suite of N-body simulations that span the full parameter-
space that is consistent with the constrained orbital and planetary
properties. We study the stability of each system and look for resonances
based on the libration of the critical resonant angles. We find strong
evidence for a two-body resonance in each system; we confirm a 3:2 resonance
between Kepler-226c and Kepler-226d, confirm a 3:2 resonance between
Kepler-254c and Kepler-254d, and confirm a three-body 1:2:3 resonant chain
between the three planets of Kepler-363. We explore the dynamical history of
two of these systems and find that these resonances most likely formed without
migration. Migration leads to the libration of the three-body resonant angle,
but these angles circulate in both Kepler-254 and Kepler-363. Applying our
methods to additional near-resonant systems could help us identify which
systems are truly resonant or non-resonant and which systems require
additional follow-up analysis.

Exoplanet dynamics (490), Exoplanet migration (2205), Exoplanet structure
(495)

##  1 Introduction

While in operation, the Kepler space telescope discovered over 4,500 planet
candidates during both the Kepler and K2 missions. Today, many of these
candidates have been confirmed, and Kepler-era exoplanets have contributed to
the growth of the confirmed exoplanet catalog to over 5,000 and the catalog of
candidate planets to over 8,500. This large sample size has led to many
investigations into planetary composition, formation, dynamics, and evolution
through astrobiological studies.

One intrigue raised by these studies is mean-motion resonance (MMR). MMR
occurs when two or more orbiting bodies periodically exert gravitational
perturbations on each other, leading to a repeated exchange of energy and
angular momentum. We can predict MMR by observing the orbital frequency of
neighboring planets. If in resonance, the ratio of neighboring planets’
periods will reduce to a ratio of small integers, such as 2:1 or 12:5.
However, determining resonance requires a deeper study into the system’s
dynamics since a period ratio of small integers does not necessarily mean the
system is in resonance. Such in-depth studies have been conducted and
confirmed resonance in a handful of Kepler systems such as Kepler-80
(MacDonald et al., 2016), Kepler-223 (Mills et al., 2016), and K2-138
(MacDonald et al., 2022).

Mean-motion resonance can form in systems with two or more orbiting bodies.
The simplest form of MMR is the two-body resonance. Mathematically, this is
defined as the oscillation or libration of the two-body critical angle:

$\Theta_{b,c}=j_{1}\lambda_{b}+j_{2}\lambda_{c}+j_{3}\omega_{b}+j_{4}\omega_{c}+j_{5}\Omega_{b}+j_{6}\Omega_{c}$
(1)

where $\lambda_{p}$ is the mean longitude of planet $p$, $\omega_{p}$ is the
argument of periapsis, $\Omega_{p}$ is the longitude of the ascending node,
$j_{i}$ are coefficients which sum to zero, and planet $b$ orbits closer to
the host star than planet $c$.

In systems with three or more orbiting bodies, numerous bodies may be in
resonance, either in a chain of two-body resonances or in a three or more body
resonance. A zeroth-order three-body MMR is defined by the difference of the
two-body resonant angles:

$\phi_{b,c,d}=\Theta_{c,d}-\Theta_{b,c}=m\lambda_{d}-(m+n)\lambda_{c}+n\lambda_{b}$
(2)

where $\lambda_{p}$ is the mean longitude of planet $p$, and $m$ and $n$ are
integers. This angle is independent of all longitudes of periapsis
($\bar{\omega}=\Omega+\omega$), making it ideal for resonant study in systems
with poorly constrained orbital angles and eccentricity.

Traditionally, such resonances are confirmed if all solutions to the system’s
RV or TTV forward modeling lead to librating angles. Unfortunately, few
systems produce large enough perturbations that could be detected with a
typical survey cadence (30 minute cadence from photometry and $\sim$few day
cadence from radial velocities). Due to a lack of high-precision measurements
of these systems, we must model all solutions to a system—across all potential
parameters that are consistent with the data—to confirm resonance. In the case
that all solutions result in the planets locked in MMRs, we are able to
confirm resonance in the system.

MacDonald et al. (2022) were the first to confirm a resonance without forward
modeling either transit times or the radial velocity signal of the planets.
They found that three of the planets of K2-138 are locked in a resonant chain
in 99% of $N$-body simulations that spanned the entirety of parameter space
that was previously constrained by both photometry and radial velocity
measurements, providing a method of MMR confirmation in the absence of high-
cadence, high-precision data.

Such a method, if applied on a larger scale to more systems, would enable us
to confirm more resonances. Since resonances allow for the constraint of
planetary properties, the system’s formation history, and the planets’ long
term stability, a significant number of confirmed resonant systems would allow
us to start leveraging these dynamics to better understand planet formation
and evolution.

Here, we perform such an analysis on five systems: Kepler-226, Kepler-254,
Kepler-363, Kepler-1542, and K2-32. Each of these systems was suggested to be
a “broken,” full-system 3:2 resonant chain, where the discovery of an
additional planet would complete the chain (Christiansen et al., 2018).
However, the period ratios of adjacent known planets suggest the presence of
resonant chains. Very few known systems with similar architecture exist
(Livingston et al., 2018), and confirmation of such a chain can provide
valuable insight into the dynamics, history, and composition of systems of
this architecture.

In Section 2, we briefly describe the five systems we study and discuss the
initial conditions and parameters of our N-body simulations. We then present
our results and analyze the resonant configurations of each system in Section
3. For two of the systems in which we confirm resonance, we use the resonances
to constrain the planetary masses and orbital periods and discuss forming the
chain in Section 4 before summarizing and concluding our work in Section 5.

##  2 Methods

Kepler-226 is a G-type star hosting a super-Earth and two Earth-sized planets
with orbital periods between 4 and 8 days. These three planets could be locked
in a 2:3:4 resonant chain. Since their initial confirmation (Rowe et al.,
2014), the anti-correlated TTVs of planets b and c constrained their masses to
$M_{b}=24.0^{+11.8}_{-10.1}~{}M_{\oplus}$ and
$M_{c}=45.2^{+22.5}_{-19.1}~{}M_{\oplus}$, although the radii of these two
planets ($R_{b}=1.64~{}R_{\oplus}$ and $R_{c}=2.47~{}R_{\oplus}$, Berger et
al., 2018) suggest these values to be overestimates. Although the TTVs and
period ratios of the system suggest this chain of resonances, the specific
dynamics of the system have yet to be explore.

Kepler-254 is a relatively dim ($V=16.012$) G-type star, hosting three
confirmed exoplanets with orbital periods ranging from 5.8 days to 18.7 days.
The period ratios of adjacent planets suggest the system could be locked in a
1:2:3 resonant chain. Jontof-Hutter et al. (2021) suggest that Kepler-254d and
Kepler-254c could be locked in a 3:2 resonance. However, the orbital dynamics
of Kepler-254 have yet to be included in an in-depth study to confirm MMRs.

Kepler-363 is a relatively bright ($V=13.472$) G-type star, hosting three
confirmed exoplanets. These planets orbit their stair fairly rapidly, with
orbital periods ranging from 3.6 days to 11.9 days. The period ratios of
adjacent planets suggest the system could be locked into a 1:2:3 resonant
chain. The orbital dynamics of Kepler-363 have yet to be included in any in-
depth study to confirm resonance in the system.

Kepler-1542 is a G-type star that hosts four transiting planets and one
planetary candidate, all smaller than Earth and orbiting within 8 days. The
orbital periods of the planets suggest a chain of resonances of 4:3, 5:4, 7:6,
and 6:5 if we include the candidate. Validated by Morton et al. (2016), the
four planets have never been included in an in-depth study of the system.

K2-32 is a G-type star in a binary system, hosting four transiting planets.
The innermost planet K2-32e was most recently discovered and validated by
Heller et al. (2019), suggesting that these four planets are in a 1:2:5:7
chain of mean motion resonances. Although the orbital periods suggest this
resonance, as do many follow-up studies (e.g., Lillo-Box et al., 2020), the
dynamics of this system have yet to be explored.

Following the methods of MacDonald et al. (2022), we seek to understand the
dynamics of these systems by running N-body simulations using the python
module REBOUND (Rein & Liu, 2012). We run a suite of 1000 simulations, drawing
initial values for planetary masses, inclinations, and orbital periods from
independent, normal distributions that are centered on values constrained by
current photometry. For Kepler-226, Kepler-254, and Kepler-363, we use the
results from Thompson et al. (2018) for all parameters except planetary radii,
for which we use the updated stellar, and therefore planetary, radii from
Berger et al. (2018). For Kepler-1542, we use parameters from Morton et al.
(2016), and for K2-32 we use the values from Heller et al. (2019). For planets
without mass constraints, we draw masses from the mass-radius relationship
described in Weiss & Marcy (2014)111We explore a large range of masses for
each planet and use the resulting resonances to constrain the planet masses.
We therefore are not sensitive to any specific mass-radius relationship.. Each
simulation therefore initializes with a set of parameters that is unique from
other simulations but consistent with current data. Using the WHFast
integrator (Rein & Tamayo, 2015), we integrate the modeled systems for 10 Myr
with a timestep of 5% the innermost planet’s period. We summarize the
simulation initial conditions for our simulations in Table 1.

Table 1:  Planetary Properties for Determining Resonance

Kepler-226 | b | c | d |
---|---|---|---|---
$P$ [d] | $3.940997\pm 0.000020$ | $5.349555\pm 0.000014$ | $8.109044\pm 0.000094$ |
t0 [d] | $69.09337$ | $104.80599$ | $65.80333$ |
$i$ [∘] | 88.88$\pm$0.2 | 89.62$\pm$0.2 | 89.92$\pm$0.2 |
$M_{p}$ [$M_{\oplus}$] | $4.271^{+1.933}_{-1.825}{}^{*}$ | $6.237^{+2.071}_{-1.952}{}^{*}$ | $2.440^{+1.985}_{-1.244}{}^{*}$ |
Kepler-254 | b | c | d |
$P$ [d] | $5.82666\pm 0.00001$ | $12.41218\pm 0.00008$ | $18.7464\pm 0.0001$ |
t0 [d] | $106.01$ | $75.54$ | $80.13$ |
$i$ [∘] | 89.88$\pm$0.2 | 89.95$\pm$0.2 | 89.11$\pm$0.2 |
$M_{p}$ [$M_{\oplus}$] | 8.84${}^{+2.02}_{-1.94}{}^{*}$ | 5.75${}^{+1.99}_{-2.05}{}^{*}$ | 6.72${}^{+2.03}_{-1.98}{}^{*}$ |
Kepler-363 | b | c | d |
$P$ [d] | $3.61460279\pm 0.00003$ | $7.54235832\pm 0.00004$ | $11.93205399\pm 0.00005$ |
t0 [d] | $67.695$ | $245965.961$ | $245975.106$ |
$i$ [∘] | 86.02$\pm$0.2 | 88.44$\pm$0.2 | 89.52$\pm$0.2 |
$M_{p}$ [$M_{\oplus}$] | 3.05${}^{+1.83}_{-1.65}{}^{*}$ | 4.67${}^{+2.12}_{-1.90}{}^{*}$ | 5.34${}^{+2.06}_{-1.94}{}^{*}$ |
Kepler-1542 | c | b | e | d
$P$ [d] | $2.8922302\pm 1.472e-05$ | $3.95116882\pm 1.633e-05$ | $5.10115756\pm 2.409e-05$ | $5.99273738\pm 2.26e-05$
t0 [d] | $65.86465$ | $67.22178$ | $65.42378$ | $64.74864$
$i$ [∘] | 89.89$\pm$0.2 | 88.05$\pm$0.2 | 89.68$\pm$0.2 | 88.08$\pm$0.2
$M_{p}$ [$M_{\oplus}$] | $0.429^{+0.386}_{-0.228}{}^{*}$ | $0.803^{+0.823}_{-0.420}{}^{*}$ | $0.850^{+0.801}_{-0.445}{}^{*}$ | $1.083^{+0.979}_{-0.570}{}^{*}$
K2-32 | e | b | c | d
$P$ [d] | $4.34882^{+0.00069}_{-0.00075}$ | $8.991828^{0.000083+}_{-0.000084}$ | $20.66186^{+0.00102}_{-0.00098}$ | $31.7142^{+0.0011}_{-0.0010}$
t0 [d] | $1998.886$ | $2000.92713$ | $1999.42271$ | $2003.7913$
$i$ [∘] | 90.0∗∗ | $89.1\pm 0.7$ | $89.3\pm 0.9$ | $89.3\pm 0.9$
$M_{p}$ [$M_{\oplus}$] | $1.095^{+2.248}_{-0.625}{}^{*}$ | $16.5^{+2.7}_{-2.7}$ | $<12.1$ | $10.3^{+4.8}_{-4.3}$

Note. — Initial conditions used for the simulations, including orbital period
$P$, mid-transit time $t_{0}$, sky-plane inclination $i$, and planetary mass
$M_{p}$. We initialize all planets on circular orbits. We use the values
published by Rowe et al. (2014) for all parameters of Kepler-226, Kepler-254,
and Kepler-363, except for planetary radii, where we use the updated stellar
and therefore planetary radii from Berger et al. (2018). For Kepler-1542, we
use parameters from Morton et al. (2016), and for K2-32 we use the values from
Heller et al. (2019). We assume stellar masses of $0.831~{}M_{\odot}$
(Thompson et al., 2018), $0.943~{}M_{\odot}$ (Berger et al., 2018),
$1.173~{}M_{\odot}$ (Thompson et al., 2018), $0.933~{}M_{\odot}$ (Thompson et
al., 2018), and $0.856~{}M_{\odot}$ (Heller et al., 2019) for the stars as
ordered in the table. All parameters were drawn from independent, normal
distributions, centered on the nominal values with widths equal to the value’s
uncertainty; for parameters with unequal upper and lower uncertainties, we
take the larger uncertainty as the width.

∗ planetary masses were drawn from the mass-radius relation Weiss & Marcy
(2014).

∗∗ At the time of this work, no estimate existed for this value, so we fix the
parameter and do not draw it from a normal distribution.

##  3 Results

For each of our five systems of interest, we run a suite of 1000 N-body
simulations for 10 Myr and analyze the results of each suite for two-body and
three-body resonances. We stop integrations when any planet experiences a
close encounter, defined by a distance of less than three Hill radii. To
confirm a chain of resonances, we search for simulations where the three-body
angle is librating or where both of the two-body angles are librating.

We find it unlikely that Kepler-1542 and K2-32 contain any resonant chains;
for each of these systems, no three-body angle librated in our simulations,
regardless of planetary mass. In Kepler-1542, the resonant angle
$\Theta_{e,d}=7\lambda_{d}-6\lambda_{e}-\omega_{e}$ librated in 82% of
simulations, and in K2-32 the resonant angles
$\Theta_{c,d}=3\lambda_{d}-2\lambda_{c}-\omega_{c}$ and
$\Theta_{e,b}=2\lambda_{b}-1\lambda_{e}-\omega_{e}$ librated in 70% and 68% of
simulations, respectively. Because not all solutions to our current data lead
to these angles librating, we cannot claim the planets are in resonance.

In Kepler-226, we find that the two-body angle
$\Theta^{\prime}_{c,d}=3\lambda_{d}-2\lambda_{c}-\omega_{d}$ librates about
180∘ in 99.8% of our simulations, but with large libration amplitudes of
90.5${}^{+23.19}_{15.22}$. The two-body angle
$4\lambda_{c}-3\lambda_{b}-\omega_{c}$ librates in 42% of our simulations, and
the three-body angle circulates in all simulations. While we are therefore
able to confirm the 3:2 resonance between Kepler-226c and Kepler-226d, we are
not able to confirm a resonant chain.

We focus the rest of this work on the two remaining systems, Kepler-254 and
Kepler-363. We summarize the results of the resonance analysis for all systems
in Table 2.

Table 2:  Resonance Results

Angle | % librating | Center [∘] | Amplitude [∘]
---|---|---|---
K2-32 | stable = 984 | resonant = 664 |
$\Theta_{e,b}=2\lambda_{b}-\lambda_{e}-\omega_{e}$ | 67.58% | -0.005 ${}^{+0.349}_{-0.315}$ | 48.4 ${}^{+23.8}_{-20.2}$
$\Theta_{b,c}=2\lambda_{c}-\lambda_{b}-\omega_{b}$ | 14.43% | 0.036 ${}^{+0.501}_{-0.513}$ | 58.3 ${}^{+14.0}_{-31.8}$
$\Theta_{c,d}=3\lambda_{d}-2\lambda_{c}-\omega_{c}$ | 69.92% | 0.015 ${}^{+2.073}_{-2.065}$ | 64.5 ${}^{+9.7}_{-18.9}$
$\Theta^{\prime}_{e,b}=2\lambda_{b}-\lambda_{e}-\omega_{b}$ | 0.90% | -5.38${}^{+12.75}_{-14.09}$ | 134.64${}^{+6.23}_{-4.81}$
$\Theta^{\prime}_{b,c}=2\lambda_{c}-\lambda_{b}-\omega_{c}$ | 0.00% |  $\cdots$ |  $\cdots$
$\Theta^{\prime}_{c,d}=3\lambda_{d}-2\lambda_{c}-\omega_{d}$ | 6.80% | -0.06${}^{+15.61}_{-12.71}$ |  $\cdots$
Kepler-226 | stable = 998 | resonant = 457 |
$\Theta_{b,c}=\textbf{4}\lambda_{c}-\textbf{3}\lambda_{b}-\omega_{c}$ | 42.00% | -0.052 ${}^{+0.571}_{-0.504}$ | 119.1 ${}^{+22.2}_{-20.9}$
$\Theta_{c,d}=3\lambda_{d}-2\lambda_{c}-\omega_{c}$ | 45.80% | -0.05 ${}^{+0.91}_{-0.88}$ | 135.03 ${}^{+11.8}_{-30.3}$
$\Theta^{\prime}_{b,c}=\textbf{4}\lambda_{c}-\textbf{3}\lambda_{b}-\omega_{b}$ | 41.60% | 179.9 ${}^{+0.577}_{-0.428}$ | 137.25 ${}^{+9.21}_{-11.85}$
$\Theta^{\prime}_{c,d}=3\lambda_{d}-2\lambda_{c}-\omega_{d}$ | 99.8% | 179.9 ${}^{+1.38}_{-1.11}$ | 90.5 ${}^{+23.19}_{-15.22}$
Kepler-254 | stable = 996 | resonant = 422 |
$\Theta_{b,c}=2\lambda_{c}-\lambda_{b}-\omega_{c}$ | 42.40% | 0.021 ${}^{+0.35}_{-0.39}$ | 118.6 ${}^{+21.48}_{-49.9}$
$\Theta_{c,d}=3\lambda_{d}-2\lambda_{c}-\omega_{c}$ | 99.20% | -0.15 ${}^{+2.41}_{-2.22}$ | $65.1^{+4.6}_{-5.0}$
$\Theta^{\prime}_{b,c}=2\lambda_{c}-\lambda_{b}-\omega_{b}$ | 0.00% |  $\cdots$ |  $\cdots$
$\Theta^{\prime}_{c,d}=3\lambda_{d}-2\lambda_{c}-\omega_{d}$ | 100.0% | 179.98 ${}^{+1.79}_{-1.81}$ | 87.1 ${}^{+12.34}_{-14.21}$
Kepler-363 | stable = 998 | resonant = 924 |
$\Theta_{b,c}=2\lambda_{c}-\lambda_{b}-\omega_{c}$ | 99.2% | 0.0029 ${}^{+0.224}_{-0.243}$ | 35.1${}^{+30.0}_{-17.8}$
$\Theta_{c,d}=3\lambda_{d}-2\lambda_{c}-\omega_{c}$ | 92.59% | -0.02 ${}^{+0.54}_{-0.44}$ | 55.1 ${}^{+13.9}_{-13.7}$
$\Theta^{\prime}_{b,c}=2\lambda_{c}-\lambda_{b}-\omega_{b}$ | 0.0% |  $\cdots$ |  $\cdots$
$\Theta^{\prime}_{c,d}=3\lambda_{d}-2\lambda_{c}-\omega_{d}$ | 98.8% | 179.97 ${}^{+0.41}_{-0.37}$ | 96.98${}^{+34.82}_{-35.54}$
Kepler-1542 | stable = 897 | resonant = 0 |
$\Theta_{c,b}=4\lambda_{b}-3\lambda_{c}-\omega_{b}$ | 9.81% | 0.16 ${}^{+0.66}_{-0.67}$ | 64.5 ${}^{+8.7}_{-24.3}$
$\Theta_{b,e}=5\lambda_{e}-4\lambda_{b}-\omega_{b}$ | 5.13% | -0.17 ${}^{+1.30}_{-0.88}$ | 74.1 ${}^{+2.9}_{-7.4}$
$\Theta_{e,d}=7\lambda_{d}-6\lambda_{e}-\omega_{e}$ | 81.94% | -0.05 ${}^{+0.88}_{-0.81}$ | 61.2 ${}^{+10.8}_{-16.5}$
$\Theta^{\prime}_{c,b}=4\lambda_{b}-3\lambda_{c}-\omega_{c}$ | 0.50% | -5.19${}^{+3.04}_{-0.85}$ | 132.50${}^{+3.34}_{-1.50}$
$\Theta^{\prime}_{b,e}=5\lambda_{e}-4\lambda_{b}-\omega_{e}$ | 2.80% | -2.47${}^{+12.97}_{-8.74}$ | 127.21${}^{+3.25}_{-2.12}$
$\Theta^{\prime}_{e,d}=7\lambda_{d}-6\lambda_{e}-\omega_{d}$ | 29.2% | 1.08${}^{+11.11}_{-13.93}$ | 131.72${}^{+7.99}_{-4.75}$

Note. — For each system, the number of simulations out of 1000 that survived
10 Myr, the number of simulations where all planets participate in the chain,
then, for each angle, the percentage of simulations where the angle librates
and the center and amplitude of the libration. For each system, all three-body
angles were circulating.

###  3.1 Kepler-254

Through our analysis, we find that nearly all (99.6%) simulations of
Kepler-254 remained stable during the 10 Myr integrations, i.e. no planets
experienced a close encounter or were ejected, regardless of initial parameter
values. Of these simulations, 42.4% result in a 1:2:3 three-body resonant
chain. The two-body angle $\Theta_{b,c}=2\lambda_{c}-\lambda_{b}-\omega_{c}$
librates in 42.4% of the simulations, and the two-body angle
$\Theta_{c,d}=3\lambda_{d}-2\lambda_{c}-\omega_{c}$ librates in 99.2% of the
simulations. The three-body angle
$\phi_{1}=3\lambda_{d}-4\lambda_{c}+\lambda_{b}$ circulated in all of the
simulations. We show the evolution of one the N-body simulations in Figure 1.
Given these results, we are therefore able to confirm a two-body resonance
between Kepler-254c and Kepler-254d where the angle $\Theta_{c,d}$ librates
around 0∘ with an amplitude of $65.1^{+4.6}_{-5.0}$. A three-body resonant
chain is probable but requires further analysis and more precise orbits to
confirm. The system could therefore benefit from follow-up observation and
analysis.

###  3.2 Kepler-363

Regardless of the initial parameters, nearly all 1000 simulations of
Kepler-363 remained stable for the 10 Myr integration. We find the 2:1
resonant angle $\Theta_{b,c}=2\lambda_{c}-\lambda_{b}-\omega_{c}$ librates in
99.2% of simulations, and the 3:2 resonant angle
$\Theta_{c,d}=3\lambda_{d}-2\lambda_{c}-\omega_{c}$ librates in 92.6% of
simulations. Of all 1000 simulations, 92.4% result in a three-body 1:2:3
resonant chain.

The two-body angles $\Theta_{b,c}$ and $\Theta_{c,d}$ librate about 0∘ with
moderate amplitudes of 35.1${}^{+30.0}_{-17.8}$ and 55.1${}^{+13.9}_{-13.7}$,
respectively, and the two-body angle and $\Theta^{\prime}_{c,d}$ librates
around 180∘ with large amplitudes of 96.98${}^{+34.82}_{-35.54}$. Curious
enough, the three-body angle $\phi=3\lambda_{d}-4\lambda_{c}+\lambda_{b}$ does
not librate in any of our simulations. We discuss the implications of this
circulating angle in more detail in Section 4. We show the evolution of one
the N-body simulations in Figure 2.

Figure 1: Example evolution of the orbital periods, eccentricities,
inclinations, all four two-body resonant angles, and the three-body resonant
angle of the three planets of Kepler-254. We find that the two-body angle
$\Theta_{c,d}$ librates in nearly all of our simulations, the two-body angle
$\Theta_{b,c}$ only librates in approximately 40%, and the corresponding
three-body angle circulates in each one. The initial values for this
simulation were drawn from independent, normal distributions, as described in
Section 2 and summarized in Table 1. We integrate this simulation beyond 10
Myr for visualization purposes.

Figure 2: Example evolution of the orbital periods, eccentricities,
inclinations, all four two-body resonant angles, and the three-body resonant
angle of the three planets of Kepler-363. We find that the two-body angles
$\Theta_{b,c}$, $\Theta_{c,d}$, and $\Theta^{\prime}_{c,d}$, librate in nearly
all of our simulations, but the corresponding three-body angle circulates in
each one. The initial values for this simulation were drawn from independent,
normal distributions, as described in Section 2 and summarized in Table 1. We
integrate this simulation beyond 10 Myr for visualization purposes.

##  4 Discussion

With the confirmation of resonance, we are able to study additional
information about a system and its planets. In particular, resonances allow us
to constrain planetary masses and orbits and to explore the formation and
subsequent dynamical history of the planets.

###  4.1 Using Resonance to Constrain Masses and Orbits

We explore the differences in planetary parameters between simulations that
resulted in resonance and those that did not. We perform a two-sample
Kolmogorov-Smirnov test, exploring the null hypotheses that the masses,
eccentricities, and orbital periods of the planets in resonance and the
planets not in resonance are drawn from the same distribution. As an example,
we take the distribution of masses of Kepler-363b from simulations where
$\Theta_{b,c}$ librates as one sample for the K-S test, and the distribution
of that planet’s mass from simulations where the same angle circulates as the
second sample.

For all parameters except the eccentricity of Kepler-363c, we recover large
p-values (p$~{}>~{}$0.05) and fail to reject the null hypothesis. For
Kepler-363c’s eccentricity, we recover a p-value of 0.018, suggesting that the
two distributions are statistically different. We find that the resulting
eccentricity for simulations with a librating $\Theta_{c,d}$ is smaller than
for those with a circulating $\Theta_{c,d}$ (2.3${}^{+1.8}_{-1.4}\times
10^{-4}$ and 3.0${}^{+1.5}_{-1.7}\times 10^{-4}$, respectively). Although we
are not able to use the system’s resonances to constrain the planets’ masses,
we do find that this system’s resonant state is not very dependant on the
planetary masses, confirming that more precise mass measurements are not
necessary to confirm these resonances.

###  4.2 Constraining dynamical history

With confirmed resonances, we are now able to study each system’s formation
and evolution. Although resonant chains are typically seen as the hallmark of
disk-driven migration, two additional pathways exist to form resonant chains
that are each consistent with in situ formation (MacDonald & Dawson, 2018).
Following the prescription of MacDonald & Dawson (2018), the three chain
formation pathways are long-scale migration (LM; hypothesizes the planets were
formed both further from their star and each other when compared to current
observations), short-scale migration (SM; planets formed near current
observations, just outside of resonance, where small shifts in the planets’
semi-major axes will lead to resonance), and eccentricity dampening (ECC;
planets formed near current observations, just outside of resonance, where
damping to the planets’ eccentricities will lead to resonance).

To study the formation of the resonances in Kepler-254 and Kepler-363, we
follow the methods of MacDonald & Dawson (2018) which we briefly describe
here. For each formation pathway, we run a suite of 500 N-body simulations
with the same initial conditions shown in Table 1 except with inflated orbital
periods. We use the modify_orbits_forces routine in the REBOUNDx library
(Tamayo et al., 2020) and the WHFast integrator (Rein & Tamayo, 2015). For the
LM simulations, we initialize the inner planet at 1 au from its host star and
start the other planets just wide of the observed resonances222We round each
planet’s semimajor axis up to the nearest 0.1 au to ensure the planets start
out of resonance.. For the SM and ECC simulations, we initialize the planets a
small percentage wide of their observed orbits, where we draw this percentage
for each planet and each simulation from a normal distribution of $N[5,3]$%.
All simulations start with the planets out of resonance. We then form the
resonant chains by damping the semi-major axes and/or eccentricities of the
planets, following the prescription in Papaloizou & Larwood (2000). For the LM
and SM simulations, we damp only the outer planet’s eccentricity and semi-
major axis333Both the direction and rate of migration for each planet will
depend on conditions of the disk and are therefore unknown. By simulating the
migration of only the outer planet, we implicitly assume that the migration
timescale of the inner planets is much longer., and for the ECC simulations,
we damp the eccentricity of all planets. We draw the timescales for the semi-
major axis damping ($\tau_{a}$) and eccentricity damping ($\tau_{e}$) from
independent, log-uniform distributions of log $\tau_{a}$ = U[7, 9] yr, log
$\tau_{e}$ = U[4, 6] yr; log $\tau_{a}$ = U[6, 9] yr, log $\tau_{e}$ = U[4, 7]
yr; and log $\tau_{e}$ = U[5, 7] yr for the LM, SM, and ECC suites,
respectively. We explore a wide range of damping timescales, representing a
wide range of disk conditions, to avoid fine-tuning our simulations.

We integrate each system forward with a timestep of 5% the innermost planet’s
observed orbital period. After 5$\times 10^{6}$ years, we “turn off” the
damping effects and integrate for another 0.25 Myr to ensure stability after
the gas disk would dissipate. We then study each resulting simulation for
librating two- and three-body resonant angles.

We find we are able to produce a full three-body resonant chain in systems
like Kepler-254 and Kepler-363 through all three formation pathways. However,
each formation pathway yields unique results, which we discuss in turn below.
We summarize the centers and amplitudes of librating angles resulting from
each formation pathway in Table 3, and we compare examples from each of these
formation pathways in Figures 3 and 4.

Short-scale migration: For both systems, short-scale migration (SM suite)
results in the three-body angle $\phi=\Theta_{b,c}-\Theta_{c,d}$ librating in
some of the simulations (34% and 25% for Kepler-254 and Kepler-363,
respectively), and librating about 180∘, $\sim$285∘, and a third center with
moderate amplitudes ($\sim 10-20^{\circ}$).

Long-scale migration: Since very few of the LM simulations for Kepler-254
remained stable for the full integration time, and with only one simulation in
resonance, we are unable to perform any meaningful statistical analysis on
this suite. The long-scale migration for Kepler-363 resulted in very few
simulations where $\phi$ librates and only 27% of the stable simulations with
a three-body resonant chain.

Eccentricity-damping: We find that eccentricity-damping results in the
libration of the two-body angles $\Theta_{b,c}$, $\Theta_{c,d}$, and
$\Theta^{\prime}_{c,d}$ for both Kepler-254 and Kepler-363 in about half of
the simulations, but very rarely results in the libration of
$\Theta^{\prime}_{b,c}=2\lambda_{c}-\lambda_{b}-\omega_{b}$ or of the three-
body angle $\phi$. For Kepler-254, $\Theta_{b,c}$ and $\Theta_{c,d}$ each
librate about 0∘ with small amplitudes of $5.96^{+5.23}_{-0.62}$ and
$5.76^{+9.26}_{-1.21}$, respectively, similar to the centers we recover in
Section 3.1 but with significantly smaller amplitudes. For Kepler-363,
$\Theta_{b,c}$ and $\Theta_{c,d}$ each librate about 0∘ with amplitudes of
$4.22^{+2.45}_{-0.47}$ and $28.32^{+4.05}_{-2.13}$, respectively, similar to
the centers we recover in Section 3.2 but, again, with significantly smaller
amplitudes.

In Section 3.1, we confirmed the two-body resonance between Kepler-254c and
Kepler-254d, but we were unable to confirm a resonance between the inner
planet pair in the system. Since each of the formation pathways resulted in
the libration of this angle and therefore each pathway is possible given our
current data, we cannot select one pathway over another as more probable. We
find it likely that the resonant chain of Kepler-363 formed through
eccentricity- damping, which we discuss in more detail below in Section 4.3.

Figure 3: Example evolution of systems like Kepler-363, forming the resonant
chain through three formation pathways: eccentricity damping only, short-scale
migration, and long-scale migration. The period ratio marked as black dots is
the ratio between planets b and c, the period ratio marked as green dots is
the ratio between planets c and d, and the vertical red line indicates when we
“turn-off” the damping effects. Although each pathway is able to lock the
planets into both two-body resonances, both short-scale migration and long-
scale migration result in the libration of the three-body angle
$3\lambda_{d}-4\lambda_{c}-\lambda_{b}$ which we find to be circulating.
Figure 4: Example evolution of systems like Kepler-254, forming the resonant
chain through three formation pathways: eccentricity damping only, short-scale
migration, and long-scale migration. The period ratio marked as black dots is
the ratio between planets b and c, the period ratio marked as green dots is
the ratio between planets c and d, and the vertical red line indicates when we
“turn-off” the damping effects. Both eccentricity dampening and short-scale
migration are lock the planets into both two-body resonances while only short-
scale migration results in the libration of the three-body angle
$3\lambda_{d}-4\lambda_{c}-\lambda_{b}$ which we find to be circulating. Long-
scale migration did not lead to enough simulations remaining stable to yield
statistically significant results. Table 3:  Formation Pathway Results

Angle | % librating | Center [∘] | Amplitude [∘] | Angle | % librating | Center [∘] | Amplitude [∘]
---|---|---|---|---|---|---|---
Kepler-254 | SM | stable = 481/500 | res = 368/500 | Kepler-363 | SM | stable = 474/500 | res = 356/500
$\phi_{1}$ | 8.11 | $180.32^{+2.43}_{-2.20}$ | $15.33^{+18.7}_{-7.98}$ | $\phi_{1}$ | 11.6 | $90.99^{+19.23}_{-17.78}$ | $13.64^{+17.38}_{-9.18}$
| 17.0 | $81.75^{+28.74}_{-11.25}$ | $14.07^{+19.00}_{-8.12}$ |  | 6.3 | $180.37^{+10.54}_{-2.55}$ | $20.86^{+30.08}_{-10.24}$
| 9.1 | $287.56^{+7.94}_{-6.50}$ | $22.48^{+19.64}_{-12.91}$ |  | 7.4 | $285.27^{+5.85}_{-15.72}$ | $20.93^{+12.12}_{-12.03}$
$\Theta_{b,c}$ | 49.9 | $0.11^{+3.38}_{-0.85}$ | $7.95^{+37.17}_{-4.04}$ | $\Theta_{b,c}$ | 8.6 | $-48.88^{+14.25}_{-8.56}$ | $14.05^{+9.50}_{-7.08}$
| 10.8 | $-48.62^{+11.05}_{-8.02}$ | $10.64^{+8.53}_{-5.50}$ |  | 68.6 | $0.17^{+16.90}_{-1.06}$ | $12.09^{+28.53}_{-7.99}$
| 16.2 | $45.53^{+9.70}_{-10.71}$ | $11.48^{+8.50}_{-6.27}$ | $\Theta_{c,d}$ | 82.1 | $0.02^{+2.63}_{-1.71}$ | $15.72^{+14.21}_{-12.57}$
$\Theta_{c,d}$ | 84.0 | $0.04^{+3.40}_{-1.09}$ | $7.44^{+20.09}_{-5.34}$ | $\Theta^{\prime}_{b,c}$ | 13.1 | $279.32^{+16.34}_{-26.93}$ | $8.53^{+9.41}_{-4.82}$
$\Theta^{\prime}_{b,c}$ | 21.2 | $287.45^{+6.09}_{-33.69}$ | $7.24^{+8.55}_{-3.98}$ |  | 8.4 | $65.87^{+10.33}_{-6.30}$ | $11.62^{+6.35}_{-5.09}$
| 11.6 | $66.59^{+8.83}_{-4.17}$ | $10.12^{+8.16}_{-5.96}$ |  | 7.2 | $179.04^{+3.88}_{-21.02}$ | $31.02^{+43.83}_{-21.05}$
| 6.7 | $179.70^{+2.20}_{-2.93}$ | $15.55^{+36.56}_{-8.64}$ | $\Theta^{\prime}_{c,d}$ | 82.9 | $179.97^{+1.42}_{-1.58}$ | $11.30^{+18.78}_{-8.21}$
$\Theta^{\prime}_{c,d}$ | 81.3 | $179.99^{+14.76}_{-1.02}$ | $7.78^{+18.20}_{-6.24}$ |  |  |  |
Kepler-254 | ECC | stable = 500/500 | res = 238/500 | Kepler-363 | ECC | stable = 500/500 | res = 229/500
$\phi_{1}$ | 0.0 |  $\cdots$ |  $\cdots$ | $\phi_{1}$ | 0.6 |  $\cdots$ |  $\cdots$
$\Theta_{b,c}$ | 47.6 | $-0.017^{+0.587}_{-0.541}$ | $5.96^{5.23}_{-0.62}$ | $\Theta_{b,c}$ | 46.8 | $0.019^{+0.392}_{-0.412}$ | $4.22^{+2.45}_{-0.47}$
$\Theta_{c,d}$ | 62.0 | $-0.013^{+0.478}_{-0.522}$ | $5.76^{+9.26}_{-1.21}$ | $\Theta_{c,d}$ | 47.0 | $-0.031^{+2.242}_{-1.735}$ | $28.32^{+4.05}_{-2.13}$
$\Theta^{\prime}_{b,c}$ | 0.0 |  $\cdots$ |  $\cdots$ | $\Theta^{\prime}_{b,c}$ | 1.4 |  $\cdots$ |  $\cdots$
$\Theta^{\prime}_{c,d}$ | 62.2 | $180.00^{+0.31}_{-0.22}$ | $3.69^{+10.61}_{-0.83}$ | $\Theta^{\prime}_{c,d}$ | 47.2 | $179.99^{+0.84}_{-0.79}$ | $15.57^{+1.46}_{-0.89}$
Kepler-254 | LM | stable = 13/500 | res = 1/500 | Kepler-363 | LM | stable = 73/500 | res = 21/500
$\phi_{1}$ | 0.2 |  $\cdots$ |  $\cdots$ | $\phi_{1}$ | 1.2 |  $\cdots$ |  $\cdots$
$\Theta_{b,c}$ | 0.2 |  $\cdots$ |  $\cdots$ | $\Theta_{b,c}$ | 6.4 | $0.28^{+6.81}_{-2.38}$ | $9.16^{+50.69}_{-5.05}$
$\Theta_{c,d}$ | 0.2 |  $\cdots$ |  $\cdots$ | $\Theta_{c,d}$ | 6.8 | $0.19^{+5.94}_{-1.98}$ | $40.42^{+38.66}_{-38.84}$
$\Theta^{\prime}_{b,c}$ | 0.2 |  $\cdots$ |  $\cdots$ | $\Theta^{\prime}_{b,c}$ | 1.2 |  $\cdots$ |  $\cdots$
$\Theta^{\prime}_{c,d}$ | 0.2 |  $\cdots$ |  $\cdots$ | $\Theta^{\prime}_{c,d}$ | 8.8 | $179.99^{+2.97}_{-3.84}$ | $26.93^{+20.19}_{-23.70}$

Note. — For each system, the number of simulations that survived the full
integration, the number of simulations where all planets participate in a
1:2:3 chain, then, for each angle, the percentage of simulations where the
angle librates and the center and amplitude of the libration. We do not
include center or amplitude data for angles librating in fewer than 5% of
simulations.

###  4.3 Unique Dynamical Configuration

The three planets of Kepler-363 are locked in a three-body resonance, where
both two-body angles librate and the three-body angle
$\phi=\Theta_{c,d}-\Theta_{b,c}=3\lambda_{d}-4\lambda_{c}+\lambda_{b}$
circulates; the three-body angle even circulates in most of our chain-
formation simulations (see Table 2). Typically, the three-body angle will
librate if the associated two-body angles librate444Although the opposite is
not true in the case of pure three-body resonance, and so we must ask: how
could this resonant chain form without the libration of this three-body angle?
We also find that the angle
$\Theta^{\prime}_{b,c}=2\lambda_{c}-\lambda_{b}-\omega_{b}$ always circulates
in our simulations and the angle
$\Theta^{\prime}_{c,d}=2\lambda_{d}-\lambda_{c}-\omega_{d}$ always librates in
our simulations. We can use all five resonant angles ($\Theta_{b,c}$,
$\Theta_{c,d}$, $\Theta^{\prime}_{b,c}$, $\Theta^{\prime}_{c,d}$, and $\phi$)
to study the possible formation history of Kepler-363; a likely formation
pathway would result in systems with dynamics similar to those we observe:
$\Theta_{b,c}$, $\Theta_{c,d}$, and $\Theta^{\prime}_{c,d}$ are librating but
$\Theta^{\prime}_{b,c}$ and $\phi$ are circulating.

Short-scale migration: The angle $\Theta^{\prime}_{b,c}$ librates in 38% of
our SM simulations, and $\Theta^{\prime}_{c,d}$ librates in 82.9% of our SM
simulations. In addition, the three-body angle $\phi$ librates in 33% of our
simulations. If $\Theta^{\prime}_{b,c}$ and $\phi$ are indeed circulating, we
find it unlikely that the resonant chain formed through short-scale migration.

Long-scale migration: As discussed above, it is challenging to form this chain
through long-scale migration as the system becomes unstable without large
eccentricity damping. However, we still find numerous sets of initial
parameters that result in $\phi$ librating. It is therefore possible that that
this resonant chain formed through long-scale migration but requires more
fine-tuning of parameters.

Eccentricity damping: From our 500 simulations, only seven (1.4%) result in
the libration of $\Theta^{\prime}_{b,c}$, and only three (0.6%) result in the
libration of $\phi$. Of the seven simulations resulting in the libration of
$\Theta^{\prime}_{b,c}$, one simulation has only this angle librating and all
other angles circulating, one simulation does not result in $\Theta_{c,d}$
librating, one simulation results in all angles librating, and the remaining
four simulations result in all two-body angles librating. The angle $\phi$
librates in one simulation where all angles librate and in two simulations
where all other angles circulate. We therefore find that it is challenging for
$\Theta^{\prime}_{b,c}$ and $\phi$ to librate if this chain was formed without
any change in the planets’ semi-major axes.

Since we are only able to simulate the formation of resonant chains in systems
similar to Kepler-363, we caution against claims of one formation mechanism;
however, we find that the angles $\Theta^{\prime}_{b,c}$ and $\phi$ do not
librate in chains formed with eccentricity-damping when the angles
$\Theta_{b,c}$, $\Theta_{c,d}$, and $\Theta^{\prime}_{c,d}$ do librate,
resulting in the dynamics we observe. Resonant chains formed through short-
scale and long-scale migration both result in the libration of
$\Theta^{\prime}_{b,c}$ and $\phi$ in the majority of simulations where the
other angles librate.

##  5 Conclusion

Planets in mean motion resonance with one another periodically exchange energy
and angular momentum, enabling us to constrain the formation history of
individual systems and identify indicators of formation history in other
systems. Because the confirmation of resonance requires an in-depth study of a
system’s dynamics, most resonances have not been confirmed. Here, we perform
such a dynamical study of five multi-planet systems whose period ratios
suggest they could be in resonance.

For each system, we run a suite of N-body simulations, exploring the full
range of possible planetary and orbital parameters as constrained by available
data. We confirm that two planets are in resonance if their critical resonant
angle librates in at least 90% of our simulations. Kepler-1542 and K2-32 each
contain at least one planet pair that is likely in resonance, but the
uncertainties on the planet masses and orbits prohibit us from confirming
these resonances. We confirm the 3:2 resonance between Kepler-226c and
Kepler-226d, confirm the 3:2 resonance between Kepler-254c and Kepler-254d,
and confirm the 1:2:3 resonant chain between the three planets of Kepler-363.
For each of these systems, we find that the three-body critical angle
$\phi=\Theta_{c,d}-\Theta_{b,c}=3\lambda_{d}-4\lambda_{c}+\lambda_{b}$
circulates in all of our simulations, even when both $\Theta_{c,d}$ and
$\Theta_{b,c}$ librate. All five of these systems could benefit from
additional data and certainly additional analysis, as their proximity to
resonance likely results in measurable TTVs.

We explore the dynamical history of Kepler-254 and Kepler-363, integrating the
systems through three potential resonant chain formation pathways: long-scale
migration, short-scale migration, and only eccentricity damping. Under our
simple migration model, both migration pathways lead to the libration of the
three-body angle, suggesting that the resonances in these two systems are more
likely to have formed in the absence of migration.

Our methods to confirm or constrain resonances within systems in the absence
of high-precision data can be applied to other systems with near-resonant
planets and would provide a list of potential new resonances that require
further analysis. With the confirmation of new resonances and particularly new
resonant chains, we are able to fully leverage the benefits of resonances and
constrain the formation history of exoplanetary systems.

We thank the anonymous referee for the constructive review that improved this
work. The authors acknowledge use of the ELSA high performance computing
cluster at The College of New Jersey for conducting the research reported in
this paper. This cluster is funded in part by the National Science Foundation
under grant numbers OAC-1826915 and OAC-1828163.

## References

  * Berger et al. (2018) Berger, T. A., Huber, D., Gaidos, E., & van Saders, J. L. 2018, ApJ, 866, 99, doi: 10.3847/1538-4357/aada83
  * Christiansen et al. (2018) Christiansen, J. L., Crossfield, I. J. M., Barentsen, G., et al. 2018, AJ, 155, 57, doi: 10.3847/1538-3881/aa9be0
  * Heller et al. (2019) Heller, R., Rodenbeck, K., & Hippke, M. 2019, A&A, 625, A31, doi: 10.1051/0004-6361/201935276
  * Jontof-Hutter et al. (2021) Jontof-Hutter, D., Wolfgang, A., Ford, E. B., et al. 2021, The Astronomical Journal, 161, 246, doi: 10.3847/1538-3881/abd93f
  * Lillo-Box et al. (2020) Lillo-Box, J., Lopez, T. A., Santerne, A., et al. 2020, A&A, 640, A48, doi: 10.1051/0004-6361/202037896
  * Livingston et al. (2018) Livingston, J. H., Dai, F., Hirano, T., et al. 2018, The Astronomical Journal, 155, 115, doi: 10.3847/1538-3881/aaa841
  * MacDonald & Dawson (2018) MacDonald, M. G., & Dawson, R. I. 2018, The Astronomical Journal, 156, 228
  * MacDonald et al. (2022) MacDonald, M. G., Feil, L., Quinn, T., & Rice, D. 2022, The Astronomical Journal, 163, 162, doi: 10.3847/1538-3881/ac524c
  * MacDonald et al. (2016) MacDonald, M. G., Ragozzine, D., Fabrycky, D. C., et al. 2016, The Astronomical Journal, 152, 105, doi: 10.3847/0004-6256/152/4/105
  * Mills et al. (2016) Mills, S. M., Fabrycky, D. C., Migaszewski, C., et al. 2016, Nature, 533, 509
  * Morton et al. (2016) Morton, T. D., Bryson, S. T., Coughlin, J. L., et al. 2016, ApJ, 822, 86, doi: 10.3847/0004-637X/822/2/86
  * Papaloizou & Larwood (2000) Papaloizou, J. C., & Larwood, J. D. 2000, Monthly Notices of the Royal Astronomical Society, 315, 823
  * Rein & Liu (2012) Rein, H., & Liu, S.-F. 2012, Astronomy & Astrophysics, 537, A128
  * Rein & Tamayo (2015) Rein, H., & Tamayo, D. 2015, Monthly Notices of the Royal Astronomical Society, 452, 376
  * 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
  * Tamayo et al. (2020) Tamayo, D., Rein, H., Shi, P., & Hernandez, D. M. 2020, MNRAS, 491, 2885, doi: 10.1093/mnras/stz2870
  * Thompson et al. (2018) Thompson, S. E., Coughlin, J. L., Hoffman, K., et al. 2018, ApJS, 235, 38, doi: 10.3847/1538-4365/aab4f9
  * Weiss & Marcy (2014) Weiss, L., & Marcy, G. 2014, The Astrophysical Journal Letters, 783, L6, doi: 10.1088/2041-8205/783/1/L6