# A Six Planet Resonance Chain in K2-138?

M. Cerioni  C. Beaugé  Instituto de Astronomía Teórica y Experimental (IATE),
Observatorio Astronómico, Universidad Nacional de Córdoba,
Laprida 854, (X5000BGR) Córdoba, Argentina

###### Abstract

The K2-138 system hosts six planets and presents an interesting case study due
to its distinctive dynamical structure. Its five inner planets are near a
chain of 3/2 two-body mean-motion resonances, while the outermost body (planet
g) is significantly detached, having a mean-motion ratio of $n_{f}/n_{g}\sim
3.3$ with its closest neighbor. We show that the orbit of $m_{g}$ is actually
consistent with the first-order three-planet resonance (3P-MMR) characterized
by the relation $2n_{e}-4n_{f}+3n_{g}=0$ and is the first time a pure first-
order 3P-MMR is found in a multi-planet system and tied to its current
dynamical structure. Adequate values for the masses allow to trace the
dynamical history of the system from an initial capture in a 6-planet chain
(with $n_{f}/n_{g}$ in a 3/1 resonance), up to its current configuration due
to tidal interactions over the age of the star. The increase in resonance
offset with semi-major axis, as well as its large value for $n_{f}/n_{g}$ can
be explained by the slopes of the pure three-planet resonances in the mean-
motion ratio plane. The triplets slide outward over these curves when the
innermost pair is pulled apart by tidal effects, in a pantograph-like manner.
The capture into the 3P-MMR is found to be surprisingly robust given similar
masses for $m_{g}$ and $m_{f}$, and it is possible that the same effect may
also be found in other compact planetary systems.

Exoplanets (498) — Orbital resonances (1181)

##  1 Introduction

One of the things we learnt from the Kepler Space Telescope is about the
significant lack of 2-planet mean-motion resonances (2P-MMRs) among low-mass
planets ($<20m_{\oplus}$) when compared to their high-mass counterparts
($>100m_{\oplus}$) (Lissauer et al., 2011; Winn & Fabrycky, 2015).
Nevertheless, we still find rich resonant dynamics in the former domain.

On one hand, it has been recently argued that, even when lacking in pair-wise
commensurabilities, their distribution of orbital periods shows a correlation
with the web of 3-planet mean-motion resonances (3P-MMRs) (Cerioni et al.,
2022), suggesting a more complex resonant evolution than previously thought.
On the other hand, planetary systems showing a high number of MMRs among its
constituents are predominantly low-mass (e.g Kepler-60, Kepler-80, Kepler-223,
TRAPPIST-1, TOI-178, K2-138, etc.), although some examples are also found in
systems harboring massive planets (e.g. GJ876 and HR8799). These multi-
resonant systems appear as ”resonance chains” with multiple links of adjacent
2P-MMRs and (consequently or not) zero-order 3P-MMRs.

We do not expect these commensurabilities to be the outcome of in-situ
formation with fortunate initial conditions. That scheme falls short because
it fails to account for the meaningful gravitational interactions between
planetesimals and a gaseous protoplanetary disk, mainly, disk-induced
migration. This mechanism is characterised by a smooth decrease in semi-major
axis, forcing a slow inward spiraling of the bodies towards the star.

Figure 1: Sketch of K2-138. The white circle on the left corresponds to the
central star. The solid lines correspond to their semi-major axes as listed in
Table 1. Black dot sizes are proportional to the measured radii.

Depending on their relative rates, migration can lead two planets into a
resonance. Although in the case of divergent migration (where planet
separation increases) planets are unlikely to become trapped in resonance
(e.g, Henrard & Lemaitre, 1983), the case of convergent migration (where
planet separation decreases) is different. Convergent migration can lead to
permanent capture into resonance (Goldreich & Tremaine, 1980; Lee & Peale,
2002), and is believed to be the cause of most observed MMRs, such as that of
Neptune and Pluto (Malhotra, 1991), the satellites of Jupiter and Saturn, as
well as several large-mass exoplanetary systems (e.g., Lee & Peale, 2002;
Beaugé et al., 2003; Ramos et al., 2017).

It is believed that resonance capture in stable solutions of first-order MMRs
is highly probable provided sufficiently slow smooth differential migration
and low initial eccentricities (see Henrard & Lemaitre, 1983; Beaugé et al.,
2006; Batygin, 2015). Resonance chains are therefore expected to be a frequent
outcome given any system of multiple migrating low-mass planetesimals,
although many may succumb under instabilities induced by the other planets and
disrupt after disk dispersal (e.g., Thommes et al., 2008; Lega et al., 2013;
Izidoro et al., 2017).

One resonant system that recently caught our attention is K2-138, the first K2
system discovered by citizen scientists through the Exoplanet Explorers
program on the Zooniverse platform (Christiansen et al., 2018).

K2-138 is a K-Dwarf type star with mass $m_{0}=0.93\pm 0.02$ $m_{\odot}$ and
radius $R_{0}=0.86\pm 0.03$ $R_{\odot}$, hosting one super-Earth and five sub-
Neptunes (Christiansen et al., 2018; Lopez et al., 2019; Hardegree-Ullman et
al., 2021). A sketch of the system’s architecture can be seen in Figure 1. In
addition to the Kepler lightcurves, HARPS radial velocity measurements and
Spitzer lightcurves have allowed to simultaneously constrain planetary orbits,
radii and (in some cases) masses. The complete list of K2 and HARPS derived
planetary parameters are presented in Table 1.

Table 1: K2-138 planetary parameters as derived and presented in Lopez et al. (2019). Planets | P | a | e | R | m
---|---|---|---|---|---
| (days) | (AU) |  | ($R_{\oplus}$) | ($m_{\oplus}$)
b | 2.35309 ${}^{+0.00022}_{-0.00022}$ | 0.03385 ${}^{+0.00023}_{-0.00029}$ | 0.048 ${}^{+0.054}_{-0.033}$ | 1.510 ${}^{+0.110}_{-0.084}$ | 3.1 ${}^{+1.1}_{-1.1}$ |  |
c | 3.56004 ${}^{+0.00012}_{-0.00011}$ | 0.04461 ${}^{+0.00030}_{-0.00038}$ | 0.045 ${}^{+0.051}_{-0.032}$ | 2.299 ${}^{+0.120}_{-0.087}$ | 6.3 ${}^{+1.1}_{-1.2}$ |  |
d | 5.40479 ${}^{+0.00021}_{-0.00021}$ | 0.05893 ${}^{+0.00040}_{-0.00050}$ | 0.043 ${}^{+0.041}_{-0.030}$ | 2.390 ${}^{+0.104}_{-0.084}$ | 7.9 ${}^{+1.4}_{-1.3}$ |  |
e | 8.26146 ${}^{+0.00021}_{-0.00022}$ | 0.07820 ${}^{+0.00053}_{-0.00066}$ | 0.077 ${}^{+0.048}_{-0.049}$ | 3.390 ${}^{+0.156}_{-0.110}$ | 13.0 ${}^{+2.0}_{-2.0}$ |  |
f | 12.75758 ${}^{+0.00050}_{-0.00048}$ | 0.10447 ${}^{+0.00070}_{-0.00088}$ | 0.062 ${}^{+0.064}_{-0.043}$ | 2.904 ${}^{+0.164}_{-0.111}$ |  $<$ 8.69 |  |
g | 41.96797 ${}^{+0.00843}_{-0.00725}$ | 0.23109 ${}^{+0.00154}_{-0.00196}$ | 0.059 ${}^{+0.063}_{-0.040}$ | 3.013 ${}^{+0.303}_{-0.251}$ |  $<$ 25.47 |  |

As we can see, its architecture is characterised by five tightly packed inner
planets separated by a significant gap from $m_{g}$. Relative separations, as
described by mean-motion ratios between adjacent planets, are
$n_{b}/n_{c}=1.513$, $n_{c}/n_{d}=1.518$, $n_{d}/n_{e}=1.529$,
$n_{e}/n_{f}=1.544$ and $n_{f}/n_{g}=3.290$. We use these values in Figure 2
to locate adjacent planetary triplets in the space of mean-motion ratios,
where we can study their proximity to different resonances (vertical and
horizontal dashed lines for 2P-MMRs ; diagonal solid lines for 3P-MMRs).

Figure 2: Mean-motion ratio plane for K2-138. Colored dots correspond to each
set of three adjacent planets (triplets). The $x$-axis corresponds to the
mean-motion ratio between the innermost and middle planet, while the $y$-axis
is for middle and outermost. Dashed vertical and horizontal lines correspond
to 2P-MMRs, while the diagonal solid line corresponds to the (2,-5,3) 3P-MMR.
Given the large separation of $n_{f}/n_{g}\sim 3.29$, the $efg$ triplet is
shown in the right frame.

A few different features immediately stand out. The five inner planets are
located near an intersection of 2P-MMRs, suggesting a 5-planet chain of
consecutive near first-order 3/2 MMRs and the longest 3/2 chain known to date.
The green dot, corresponding to the $efg$ triplet, is separated from the rest
because of the notoriously large gap between $m_{f}$ and $m_{g}$. In addition,
the inner triplets are located on top of the zero-order $(2,-5,3)$ 3P-MMR,
plotted as the solid black line. Lastly, we note that the offsets which
separate the planets from the exact 3/2 intersection interestingly increase
with the distance to the star.

The dynamical history which led K2-138 to its current configuration remains
uncertain from several angles. For instance: Could the inner planets truly be
in a 5-planet 3/2 resonance chain in spite of their sizable offsets? The
origin of the gap between $m_{f}$ and $m_{g}$ is also unclear. Why is the
sixth planet so isolated from the rest given that all underwent a similar
migration process? Is it really detached from the inner resonance chain?

Both Lopez et al. (2019) and Hardegree-Ullman et al. (2021) speculate about
two non-transiting planets at $\sim 19$ and $\sim 30$ days, which would
complete an 8-planet chain of 3/2 MMRs and help explain the current location
of $m_{g}$. However, this is not the only hypothesis.

Considering only the confirmed planets so far, it is not obvious whether
$m_{f}$ is detached from the inner chain or not. It is possible that during
the migration process $m_{g}$ was captured in a 3/1 MMR with $m_{f}$, and the
large observed resonance offset ($\sim$0.3) is either primordial or a
consequence of dynamical evolution after the dispersal of the protoplanetary
disk. In such a case, all six planets would have formed a complete resonance
chain, at least at some point of their history.

While this idea is intriguing, it is unclear how such a large offset may have
been achieved. There are two main processes which are frequently invoked to
explain deviations from an observed mean-motion ratio with respect to the
nominal resonant value: differential migration in a flared-disk (Ramos et al.,
2017) and tidal effects (e.g. Lithwick & Wu, 2012; Delisle & Laskar, 2014).
Although these effects are dependent on the masses, in both scenarios the
offset is expected to decrease with larger semi-major axis, contrary to what
is observed. Moreover, the observed mean-motion ratio between the two
outermost planets appears notoriously large, particularly for a second-order
resonance where the family of zero-amplitude solutions deviate only slightly
with respect to the exact resonance position.

In this work, we study the cases of both a detached (non-resonant) and an
attached (resonant) outermost planet $m_{g}$ and test under which conditions
each may have occurred. Section 2 describes the migration framework as well as
the protoplanetary disk adopted for our N-body simulations of Type-I
migration. Section 3 estimates the necessary masses for the two outermost
planets leading to the detached or attached configurations. In Section 4 we
discuss the first possibility: a 5-planet resonance chain $m_{b}-m_{f}$ plus
the outermost planet $m_{g}$ as a non-resonant distant companion. Section 5
focuses on the second option and analyses how very large resonance offsets may
be achieved in multi-resonant systems. In both scenarios we analyze what would
be the necessary planetary masses and disk parameters to reproduce the
observed dynamical features and how robust these solutions are likely to be.
Finally, the pantographic effect leading to increasing offsets is explained in
detail in Section 6, while a general discussion and future prospects close the
paper in Section 7.

##  2 Planetary Migration

###  2.1 Migration Scheme

The gravitational interaction between a low-mass planet and a protoplanetary
disk results in orbital decay, usually known as Type-I planetary migration. An
analytical prescription was originally developed by Tanaka et al. (2002), who
considered a single fully formed planet in a circular orbit and a three-
dimensional isothermal gaseous disk with a surface density $\Sigma$ and aspect
ratio $H$ given by

$\Sigma(r)=\Sigma_{0}\left(\frac{r}{r_{0}}\right)^{-s},\quad
H(r)=H_{0}\left(\frac{r}{r_{0}}\right)^{f},$ (1)

where $s$ and the disk flare index $f$ characterise the surface density
distribution and shape of the disk, respectively. If $f=0$, the disk profile
is linear (flat), while $f>0$ describes a curved (flared) disk. Both
$(\Sigma_{0},H_{0})$ are the values at $r_{0}=1$ au.

The expression for the orbital decay $\tau_{a}$ was proposed by Tanaka et al.
(2002) for the case described. Later on, Tanaka & Ward (2004) extended the
model for planets in eccentric orbits, which would circularize when embedded
in a disk. The expressions for both the migration and circularization
timescales are given by

$\tau_{a}\equiv\left|\frac{a}{\dot{a}}\right|=\frac{Q_{a}t_{{\rm
w}}}{H(r)^{2}}\;\;\;;\;\;\;\tau_{e}\equiv\left|\frac{e}{\dot{e}}\right|=Q_{e}\frac{t_{{\rm
w}}}{0.78}$ (2)

where $Q_{a}^{-1}\simeq 2.7+1.1s$ and we assume the ad-hoc eccentricity
damping factor adopted by Cresswell & Nelson (2006) of $Q_{e}=0.1$. Lastly,
$t_{\rm{w}}$ is the wave timescale given by

$\quad t_{{\rm
w}}=\frac{m_{0}}{m_{p}}\frac{m_{0}}{\Sigma(r)r^{2}}\frac{H(r)^{4}}{\Omega_{K}(r)},$
(3)

where $m_{p}$ is the mass of the planet under consideration and $m_{0}$ the
mass of the star. The Keplerian angular velocity is given by
$\Omega_{K}(r)=\sqrt{G(m_{0}+m_{p})/r^{3}}$. Note that equations (2) and (3)
show that the migration velocity is directly proportional to the planetary
mass, as well as dependent on the disk parameters. Thus, the mass distribution
in a given system is largely responsible for the convergence or divergence
between two or more planets, as well as the possibility of resonance trapping
(see Beaugé & Cerioni, 2022).

In order to account for the coupling between semi-major axis and
eccentricities during the migration process, Goldreich & Schlichting (2014)
modified expressions (2) with a $\beta$ factor to quantify the conservation on
the angular momentum. The effective timescale for the radial drift is

$\frac{1}{\tau_{a_{{\rm
eff}}}}=\frac{1}{\tau_{a}}+2\beta\frac{e^{2}}{\tau_{e}}.$ (4)

When $\beta=0$ there is no coupling while $\beta=1$ implies no loss in the
orbital angular momentum. Goldreich & Schlichting (2014) suggest $\beta=0.3$
following estimations by Tanaka & Ward (2004).

###  2.2 Disk Parameters

For our N-body simulations and analytical estimates we assumed a thin laminar
flat disk ($f=0$) characterized by a surface density exponent of $s=1/2$ and
an aspect ratio of $H_{0}=0.05$ at $r_{0}=1$ au. We considered several values
of the surface density $\Sigma_{0}$ from $100$ up to $1700\mathrm{\ gr\
cm^{-2}}$. In most cases we observed no significant difference in the outcome.

These parameters were chosen for the sake of simplicity considering our lack
of information of the real disk. At the very least, we expect our results to
be robust with respect to reasonable variations of the disk parameters.
Equations (2) and (3) show that $\tau_{a}$ is directly proportional to
$H_{0}$, and inversely so to $s$, $f$, and $\Sigma_{0}$; we can therefore
predict the evolution timescale and necessary integration time of our
simulations as a function of these parameters.

While the general trend of the migrating system is expected to be fairly
robust to disk parameters, the conditions of resonance capture and the
subsequent orbital evolution of the resonance chain may be more sensitive. For
this reason, each new set of parameters led to a new calculation of capture
conditions with our analytical model and new N-body simulations. More details
will be given in the next section.

##  3 Mass estimates for detachment

As may be seen in Table 1, even if fairly reasonable estimations of planetary
masses are available for the first four planets, only upper bounds are given
for the two outermost bodies. Beaugé & Cerioni (2022) showed that the
formation of a resonance chain is particularly sensitive to the mass of the
last (and first) planet, so the estimation of ($m_{f}$,$m_{g}$) can provide
detachment regimes for the sixth planet.

In that spirit, we start by checking which combination of masses would produce
a complete 6-planet MMR chain. This calculation does not require N-body
simulations but may be undertaken using an analytical resonant capture
condition for $N$ planets developed by Beaugé & Cerioni (2022). Truth be said,
this method only checks whether a given (e.g. current) distribution of orbital
separations is an attractor under the effects of the exterior non-conservative
force and thus evolves towards more compact configuration, or whether the
evolution will tend to disperse the system.

As mentioned earlier, resonance capture conditions are sensitive to the
combination of planetary masses, so in Figure 3 we studied a grid of values
($m_{f}$,$m_{g}$), where each mass varied between zero and its upper bound.
For each pair of values we generated a set of $N_{\rm sys}=10^{3}$ fictitious
systems varying the mass of the other planets from a probability distribution
function (PDF) consistent with their observed values***For each planet $i$
with a measured mass of ${m_{i}}_{-\Delta_{L}}^{+\Delta_{R}}$, we produced a
Skew-Normal distribution with median in $m_{i}$ and such that its integral in
the interval [$m_{i}-\Delta_{L}$,$m_{i}+\Delta_{R}$] yields 0.683 (i.e the
68.3%-credibility interval) and used this PDF to draw random masses.. Then, we
checked what percentage of $N_{sys}$ led to the formation of a stable
resonance chain. The color assigned to each cell represents this ratio.

Figure 3: Percentage of fictitious 6-planet systems $m_{b}-m_{g}$ that lead to
global convergent migration and the formation of a complete resonance chain.
While $m_{f}$ and $m_{g}$ are varied in a regular grid, for each cell the
values of the other masses are drawn from a PDF consistent with Table 1.

As shown by the color code on the right, the percentage of positive values
(i.e. convergent infall of the complete 6-planet system) varies from zero for
low values of $m_{g}\lesssim 7\,m_{\oplus}$, to mixed values around $m_{g}\sim
8-12\,m_{\oplus}$, to $100\%$ for $m_{g}\gtrsim 13\,m_{\oplus}$. For the
migration of $m_{g}$ onto the inner system, these three regions respectively
represent: i) convergent migration leading to an attached $m_{g}$ and a
6-planet chain, ii) a mixed domain of slow divergence/convergence, and iii)
divergent migration leading to a detached $m_{g}$ and only an inner 5-planet
chain. These limits are only weakly dependent on the adopted value of $m_{f}$
and the rest of the planets.

We will consider these mass estimates in the following sections in order to
analyze two options: an $m_{g}$ planet which is detached from the inner
5-planet resonance chain, and $m_{g}$ planet attached to the inner system and
forming a complete 6-planet resonance chain.

##  4 Option 1: A 5-Planet Resonance Chain

###  4.1 A Non-Resonant yet Non-Hierarchical $m_{g}$

With the innermost planet currently located at $a_{b}\simeq 0.03$ au around a
$0.93m_{\odot}$ star, it is fair to say that the planets underwent a large-
scale migration from their original formation sites.

In this section we assume a non-resonant sixth planet, which means that for
$m_{g}$ we are avoiding the green region in Figure 3 and considering the red
and white areas instead. These are different mass regimes which imply
different migration rates. We will take it case by case.

To begin with, if the mass of the outer planet were significantly lower than
$7\,m_{\oplus}$ (red tones), then the infall rate of the outermost planet
would have been much smaller than that of the resonance chain, leading to a
marked hierarchical planetary system with a much larger relative separation
than observed. Recall that the mean-motion ratio between the two outer planets
is $n_{f}/n_{g}\simeq 3.3$ and not too far from the 3/1 MMR.

Therefore, in order to explain the non-resonant but also non-hierarchical
location of $m_{g}$, the only option is the white-toned region, corresponding
to $m_{g}\sim 8-12\,m_{\oplus}$. In this regime, the migration rate of $m_{g}$
can be slightly smaller or higher than that of the inner resonance chain,
leading to a very slow divergent or convergent migration, depending on the
particular combination of $m_{e},\,m_{f}$, and $m_{g}$. The scenario of slow
divergence would have seen an early formation of the $m_{b}-m_{f}$ chain
followed by a slow but persistent increase in the separation with $m_{g}$,
coming to a halt at the observed separation right at the time of dispersal of
the gas disk. The same would be true for the case of a slow convergent
migration of an initially farther out $m_{g}$ planet cut short before it could
reach the 3/1 resonance. While plausible, these alternatives are critically
sensitive to the combination of initial separation, disk lifetime and
differential migration rates. Nevertheless, given that this section only
pretends to study the formation of the inner chain with a non-resonant
$m_{g}$, we will take its value from this regime.

As a final aside, it is interesting to note that a non-resonant nature for
$m_{g}$ may be proposed as strong evidence against an inner disk edge, or, if
present, one that did not play an important role in the migration history of
the K2-138 system. If the resonance chain $m_{b}-m_{f}$ was stopped at such an
edge, then $m_{g}$ would have ultimately reached its inner companions
independently of its mass. Smaller values of $m_{g}$ would imply longer times,
but the outcome would have been the same.

###  4.2 Planetary Masses

Assuming that the sixth planet will not come into resonance with the inner
system, the next question is about the formation of the inner 5-planet chain
on its own. In Figure 4 we perform a similar analysis to that of Section 3,
this time focusing on which values of $m_{e}$ and $m_{f}$ lead to the
formation of the inner chain. The values of $m_{f}$ were varied between zero
and the upper bound given in Table 1, while for its inner neighbor we
considered values in the range $m_{e}\in[13-2,13+2]\,m_{\oplus}$. As before,
for each cell of the rectangular grid we generated a series of fictitious
systems with the masses of the other planets drawn from a PDF consistent with
the observational data.

The formation of the 5-planet resonance chain seems very unlikely for the
estimated parameters of the system. As mentioned earlier, the convergence
condition is particularly sensitive to the values of the first and last planet
in the system (Beaugé & Cerioni, 2022). Since the best fit values imply
$m_{f}<m_{e}$ (even for the upper bound of $m_{f}$), convergent migration is
highly improbable, especially for $m_{e}$ larger than its reported median
value. The color bar shows some positive results, and with a high probability,
but only for $m_{e}\sim 11\,m_{\oplus}$ and $m_{f}\simeq 8.5\,m_{\oplus}$.
While the outer planet is still less massive than its inner neighbor, an early
formation of a resonance chain comprised of $m_{b}-m_{e}$ would slow the
infall of the inner system sufficiently to allow $m_{f}$ to catch up and join
its kin.

Figure 4: Analogous to Figure 3, but now considering only the 5-planet system
comprised of $m_{b}-m_{f}$. In both plots the color bars indicates the ratio
of systems which converge to a chain. Figure 5: Left Frame: N-body simulations
under mutual gravitational perturbations and migration of the 6-planet K2-138
system. Masses were taken equal to
$\\{m_{k}\\}=(3.1,6.3,7.9,10.0,8.7,8.0)\,m_{\oplus}$, while for the disk
parameters we adopted $(s,f,\Sigma_{0},H_{0})=(0.5,0,500\ {\rm
gr/cm}^{2},0.05)$. All planets were initially placed on circular coplanar
orbits with separation such that $n_{i}/n_{(i+1)}=1.55$. Center Frame: Same as
previous plot, but with initial separations equal to $n_{i}/n_{(i+1)}=1.70$.
Right Frame: Final value of $n_{e}/n_{f}$ as function of the initial
separation.

###  4.3 Initial Separations

A more subtle difficulty in the 6-planet system is related to their initial
separation leading to the migration process. Even considering the largest
value for $m_{f}$ and reasonably low values for $m_{e}$ (say,
$10\,m_{\oplus}$), we still have mass ratios $m_{e}/m_{f}\sim 8.69/10<1$. From
Ramos et al. (2017) we know that the mass ratio between adjacent planets
$m_{i}$ and $m_{i+1}$ leading to 2-planet convergent migration must satisfy
the condition

$\frac{m_{i+1}}{m_{i}}>\left(\frac{n_{i}}{n_{i+1}}\right)^{\frac{2}{3}(2f+s-1/2)}$
(5)

where $n_{i}/n_{i+1}$ is not the resonant mean-motion ratio we wish to achieve
but the initial value at the start of the migration. A flat disk ($f=0$) and a
surface density exponent $s=1/2$ give a minimum value, i.e. $m_{i+1}/m_{i}=1$
independent of the initial separation. Other disk parameters, especially those
with a steeper surface density, lead to more stringent conditions for the mass
ratio.

Since the estimated value of $m_{e}/m_{f}$ is smaller than unity, the inner
body of the pair migrates faster and $m_{f}$ lags behind increasing any
initial separation. Convergent migration between these two planets is only
expected to occur after the inner resonance chain $m_{b}-m_{e}$ is established
and the effective semi-major axis decay rate is slowed. This behavior is
expected for all disk properties and regardless of the initial conditions,
although the resulting MMR in which each pair is captured may (and will) vary.

Nevertheless, there is one hard condition on initial separations. Given that
the mass ratios of the first four planets is greater than unity, migration is
convergent and they will be trapped in a mean-motion commensurability that is
lower than their initial separations. Therefore, we set initial separations
greater than $n_{i}/n_{(i+1)}>1.5$.

###  4.4 N-Body Simulations

We performed N-body simulations employing a Bulirsch-Stoer integration scheme
with variable time-step and precision $ll=12$ (Bulirsch & Stoer, 1966). All
bodies begin with initial separations $n_{i}/n_{(i+1)}=1.55$ and migration
rates as outlined in 2.1. The left frame of Figure 5 shows results for the
aforementioned scenario, where we plot the mean-motion ratio of each adjacent
pair of planets as a function of time.

The mass ratio increases subsequently for the first three pairs, so migration
is convergent and they rapidly are trapped in a four-planet chain of MMRs.
However, since $m_{f}<m_{e}$ the fifth planet starts the simulation falling
behind and increasing the value of $n_{e}/n_{f}$. This process reverses once
the inner system is trapped. The net infall speed of the four-planet chain is
now lower than the decay rate of $m_{f}$ (see Beaugé & Cerioni, 2022) and this
planet ultimately catches up and joins the resonance chain. No such fate
awaits the outermost planet and $m_{g}$ continues to fall behind increasing
its distance with the rest of the pack.

The scenario described above appears consistent with the overall dynamical
characteristics of the system, and we were able to obtain similar results for
a wide range of disk properties. Variations in the masses also led to
analogous results as long as the mass ratios satisfied the
convergence/divergence conditions mentioned earlier and outlined in Figures 3
and 4. Different results, however, were found when changing the initial
separation between the planets. The central frame of Figure 5 shows a similar
N-body simulation as before, but now considering a wider initial system with
$n_{i}/n_{(i+1)}=1.70$.

Since the planets start the simulation further apart, the migration of the
inner pack towards the 3/2 MMRs takes longer, giving time for $m_{f}$ to
separate even more and encounter a larger number of commensurabilities in its
path. When the differential migration reverses, $m_{f}$ does not reach the 3/2
MMR with $m_{e}$ but becomes trapped in the 5/3, leading to a configuration
inconsistent with the one observed. The right frame shows how the final value
of the ratio $n_{e}/n_{f}$ varies according to the initial separation between
adjacent planets.

Unsurprisingly, initial separations close to the 3/2 MMR yield final systems
trapped in such a commensurability and thus analogous to the observed K2-138
6-planet model. However, even slightly wider separations lead to other
resonances and other configurations, with the 5/3 and 2/1 being the most
common, depending on the initial $n_{i}/n_{(i+1)}$. Although extremely high
migration rates will sometimes avoid all these alternate capture sites and
fall back to the 3/2, this was rarely the case. The capture condition for a
given resonance depends much more strongly on the mass ratio than on the
surface density of the disk, and planets with similar masses will suffer
similar semi-major axis decay rates even for high values of $\Sigma_{0}$. We
found results analogous to Figure 5 even for $\Sigma_{0}=1700$ gr/cm2, but we
did not test higher values.

These simulations seem to indicate that the current 3/2 MMR resonance chain
among the five inner planets of K2-138 is possible provided the mass of
$m_{f}$ is close to its upper bound, the mass of $m_{e}$ is significantly
lower than its median value (beyond the 1 or 2 $\sigma$-level), and the
initial separation between the planets was such that their mean-motion ratios
all lied between $1.5$ and $1.6$. This would require that the planets form
around the same time and with very similar mean-motion ratios as those
observed today. While all these conditions are, again, possible, they seem to
require very fine tuning and some of the parameters (such as the mass ratio
$m_{e}/m_{f}$) are very far from the expected values.

Figure 6: Observed position of the $(n_{e}/n_{f},n_{f}/n_{g})$ triplet (red
cross) in the mean-motion ratio space. Diagonal curves mark every low-integer
zero- and first-order 3P-MMR that passes through the domain.
Commensurabilities crossing the $(3/1,3/2)$ intersection are labeled and drawn
using thicker lines. Error bars are too small to be distinguishable in the
figure.

##  5 Option 2: A 6-Planet Resonance Chain

###  5.1 The (2,-4,3) 3-Planet Resonance

The idea that the outermost planet ($m_{g}$) may form part of the system’s
resonance chain is strongly suggested when we analyze the location of the
triplet $(n_{e}/n_{f},n_{f}/n_{g})$ in the mean-motion ratio plane. This is
shown with a red cross in Figure 6, together with the position of every zero-
and first-order $(k_{1},k_{2},k_{3})$ 3P-MMR in the region, with indexes
$|k_{i}|\leq 10$. The 2-planet commensurabilities $n_{e}/n_{f}=3/2$ and
$n_{f}/n_{g}=3/1$ are drawn in dashed lines, while 3P-MMRs are shown in
continuous lines. Of these, three resonances cross the $(3/1,3/2)$
intersection and are highlighted with thicker lines and labeled in the plot.

The location of the three outer planets in this plane is strikingly close to
the $(2,-4,3)$ first-order three-planet resonance, characterized by the
relation $2n_{e}-4n_{f}+3n_{g}=0$. Introducing the observed values for the
orbital periods from Table 1, as well as their uncertainties, we obtain

$2n_{e}-4n_{f}+3n_{g}=0.00020^{+0.00024}_{-0.00015}\ d^{-1}$ (6)

which comes extremely close to the exact nominal commensurability. Given the
limited number of two and three-planets resonances in this region, it seems
very unlikely that such a correlation could be a coincidence, even if the
3P-MMR is of first order and thus not the strongest dynamical feature in the
vicinity. Similarly to the inner triplets, Figure 6 suggests a capture of
$(n_{e}/n_{f},n_{f}/n_{g})$ in the $(3/2,3/1)$ intersection, with an excursion
of the triplet over the three-planet resonance that crosses that intersection.

To prove that the $(2,-4,3)$ resonance is not only relevant to the current
dynamics of the system but also tied to its migration history, we must solve
two independent issues: (i) test whether planetary masses consistent with the
outer K2-138 system can be trapped in this first-order 3P-MMR and survive the
dissipation of the gaseous disk without being expelled from the
commensurability, and (ii) explain the large separation of the current
configuration with respect to the $(3/2,3/1)$ intersection of two-planet
resonances.

###  5.2 Formation of the 6-Planet Chain

In principle, there are two possible routes that could lead to the capture of
the outer triplet in the $(2,-4,3)$ resonance; these are divergent and
convergent migration.

Divergent migration occurs in a scenario similar to that discussed in the
previous section, but with a mass of the outer planet in the range
$m_{g}\lesssim 8m_{\oplus}$ (Figure 3) and initial separations such that
$n_{f}/n_{g}<3$. In the previous section we assumed all planets initially
equidistant forming a compact formation; this is not strictly necessary and
the outer planet could have began its migration further out. However, we could
not find any set of parameters that led to the capture in the $(2,-4,3)$
three-planet commensurability during divergent migration, so this case is
ruled out. As a side note, some cases did result in the zero-order $(4,-7,3)$
3P-MMR, with stable orbits and a small but non-zero offset with respect to the
two-planet resonances.

We then turn to the second possibility, a convergent migration of $m_{g}$ from
an initial orbit beyond the 3/1 resonance with $m_{f}$ in the hope of trapping
the three outermost planets in the desired first-order 3P-MMR. As before we
fix the masses of the four inner planets to the values chosen in the previous
sections: $(m_{b},\ldots,m_{e})=(3.1,6.3,7.9,10.0)m_{\oplus}$ and with initial
separations among adjacent bodies between 1.5 and 1.6. The rest of the masses,
as well as their initial separations, and disk parameters were varied.

Figure 7:  N-body simulations under mutual gravitational perturbations and
migration of the K2-138 system in the mean-motion ratio plane. Colored lines
correspond to each planetary triplet. Inner triplets belong to Set#1, but
analogous results are obtained for Set#2. The evolution of of the $efg$
triplet in both sets are shown in green colors. Initial positions are shown
with white circles. Arrows indicate the direction of evolution. Bottom frame
is a zoom-in of the resting place of the last triplet.

Figure 8: Left Frame: Time evolution of resonant angles corresponding to the
simulation Set#1. The first four plots on the left show the resonant angles
$\sigma_{1}=3\lambda_{i+1}-2\lambda_{i}-\varpi_{i}$ (in blue) and
$\sigma_{2}=3\lambda_{i+1}-2\lambda_{i}-\varpi_{i+1}$ (in red) for the first
four inner planetary pairs (i.e. $i=1,2,3,4$). The lower center plot shows the
resonant angles of the 3/1 MMR between $m_{f}$ and $m_{g}$, where blue is used
for $\theta_{1}=3\lambda_{g}-\lambda_{f}-2\varpi_{f}$, red for
$\theta_{2}=3\lambda_{g}-\lambda_{f}-\varpi_{f}-\varpi_{g}$ and orange for
$\theta_{3}=3\lambda_{g}-\lambda_{f}-2\varpi_{g}$. The capture of the pair
into the two-planet resonance is only accompanied by a libration of
$\theta_{2}$ with the other angles circulating. The green data points in the
lower-right frame show the time evolution of the resonant angle
$\theta_{243}=2\lambda_{e}-4\lambda_{f}+3\lambda_{g}-\varpi_{g}$ corresponding
to the $(2,-4,3)$ 3P-MMR of the outer triplet. Right Frame: Eccentricities of
all planets for the same simulation. The disk is smoothly removed from the
simulation at $T=375$ kyrs (vertical black line). All subsequent orbital
variations occur in a gas-free scenario.

Figure 7 shows the results of two different simulations, labeled Set#1 and
Set#2. Their main characteristics are given below:

$\begin{split}{\rm
Set\\#1:}&\;\;m_{f}=9.3\;,\;m_{g}=11.0\;,\;\Sigma_{0}=100\;{\rm g/cm^{2}}\\\
{\rm Set\\#2:}&\;\;m_{f}=8.7\;,\;m_{g}=8.8\;,\;\Sigma_{0}=1000\;{\rm
g/cm^{2}}\\\ \end{split}$ (7)

where the masses are given in units of $m_{\oplus}$. In Set#1 we tested a
value of $m_{f}$ that is slightly over the reported upper bound as low values
of this parameter prevent the formation of the inner chain of multiple 3/2
resonances (Figure 4). Moreover, the mass of the outer planet in Set#2 was set
barely large enough to guarantee convergent migration once the inner system
was trapped in its chain. This means a very slow convergent migration w.r.t
its inner neighbors and thus its initial semi-major axis could be placed close
to the 3/1 MMR with $m_{f}$; we chose $n_{f}/n_{g}\simeq 3.4$, although
slightly smaller values lead to similar results. Conversely, the mass of
$m_{g}$ in Set#1 was significantly larger leading to a convergent migration of
the outer planet onto $m_{f}$ even before the resonance lock of the inner
system. In consequence the initial separation had to be larger indicative of a
more hierarchical primordial system. We opted for $n_{f}/n_{g}\simeq 3.8$,
although larger values would result in similar final configurations.

The top frame of Figure 7 shows the evolutionary tracks of the mean-motion
ratios of all triplets, each identified with a different color. The large
circles indicate the initial values, and arrows show the direction of the
evolution. The behavior of the inner system $m_{b}-m_{f}$ was similar in both
sets, resulting in the formation of the well-known 3/2 resonance chain. The
evolution of the outer triplet is shown in green, where a dark tone was chosen
for Set#1, and a lighter color for Set#2.

Both simulations show fairly similar characteristics. Since $m_{f}<m_{e}$,
initially $n_{e}/n_{f}$ increases signaling a divergent migration of this
pair. As soon as the inner resonance chain begins to assemble, their migration
rate slows and allows $m_{f}$ to reverse its differential migration. The
evolutionary track thus switches direction and the pair $m_{e}-m_{f}$ moves
towards the 3/2 MMR where it is finally trapped. The behavior of the outer
pair is simpler as the value of $m_{g}$ guarantees convergent migration in
both cases. However, while the mass chosen for Set#1 leads to a relatively
fast inward migration from the start of the run, in Set#2 the outer pair
approaches the rest of the system only after the resonance chain begins to
form. This explains the almost horizontal evolutionary track during the first
part of the simulation.

The bottom frame of Figure 7 zooms in on the vicinity of the final resting
place of the outer triplet. After the capture of $n_{e}/n_{f}$ in the 3/2 MMR
both runs are practically indistinguishable. As the mean-motion ratio of the
outer pair $n_{f}/n_{g}$ approaches the 3/1 MMR, it first passes through the
$(4,-7,3)$ 3P-MMR but is finally stopped at the $(2,-4,3)$ commensurability.
Even though this first-order resonance should have a lower capture probability
than its near neighbor, we found it surprisingly robust and most of our
simulations led to the same result.

Figure 8 shows the time evolution of several resonant angles and
eccentricities for Set#1. The first four plots on the left frame show
$\sigma_{1}=3\lambda_{i+1}-2\lambda_{i}-\varpi_{i}$ (in blue) and
$\sigma_{2}=3\lambda_{i+1}-2\lambda_{i}-\varpi_{i+1}$ (in red) for the first
four inner planetary pairs (i.e. $i=1,2,3,4$). The inlaid box in the upper-
right corners of each frame shows the planets involved in each case. All four
pairs evolve towards an ACR (see, e.g., Beaugé et al., 2003; Michtchenko et
al., 2006) characterized by equilibrium values
$(\sigma_{1},\sigma_{2})=(0,\pi)$. As each link of the resonance chain falls
into place the resonant angles and eccentricities (right frame in the figure)
suffer a temporary kick, but are rapidly restored to low-amplitude
oscillations around their stationary solutions.

The lower center plot shows the resonant angles of the 3/1 MMR between $m_{f}$
and $m_{g}$. The capture of the pair into the two-planet resonance is only
accompanied by a libration of
$\theta_{2}=3\lambda_{g}-\lambda_{f}-\varpi_{f}-\varpi_{g}$ while all other
possible resonant angles circulate. At the same time the outer pair enters the
3/1 commensurability, the outer triplet is captured in the $(2,-4,3)$ 3P-MMR.
The green data points in the lower-right plot show the time evolution of the
resonant angle
$\theta_{243}=2\lambda_{e}-4\lambda_{f}+3\lambda_{g}-\varpi_{g}$ displaying an
asymmetric libration during the time interval up to $\sim 375$ kyrs. All other
possible resonant angles circulate.

It is difficult to prove that the capture in the three-planet resonance is a
consequence of the trapping of the individual pairs in the 2-planet
commensurabilities; the fact that the system chose the first-order 3P-MMR over
its zero-order neighbor is intriguing; additional work is necessary before
establishing in detail how this stable configuration is reached.

In order to test the stability of such a multi-resonant configuration, we
introduced a smooth reduction in the surface density of the disk during a time
interval centered around $375$ kyrs and which lasted about ten thousand years.
The subsequent evolution of the system occurs without the effects of disk
damping. The librations of both the 3/1 MMR and the three-body resonance
angles turn towards symmetric solutions as the surface density of the disk
tends to zero, with $\theta_{2}\rightarrow 0$ and
$\theta_{243}\rightarrow\pi$. The orbital evolution of the full 6-planet
resonance chain remains stable throughout the simulation with no indication of
any incipient increase in amplitudes that could signal a future disruption of
the system.

###  5.3 Tidal Evolution of the 6-Planet Chain

While we have succeeded in trapping the outermost triplet in the $(2,-4,3)$
3P-MMR and thus constructing a full 6-planet resonance chain, the offsets with
respect to the nominal 2-planet resonances are much lower than observed. This
is particularly notorious in the case of $n_{f}/n_{g}$ whose current value
($\simeq 3.29$) is so distant from the exact resonance that it does not appear
related.

The rest of the planetary pairs are also significantly displaced from the
exact 3/2 resonance, with observed mean-motion ratios equal to $1.513$,
$1.518$, $1.529$ and $1.544$, respectively. None of our capture simulations
were able to obtain similar values, regardless of the chosen planetary masses
or disk parameters. Consequently, if our previous results are representative
of the primordial system, the current resonant offsets must have originated
after the dissipation of the primordial disk; tidal effects with the central
star appears as the only possible mechanism.

Several studies have analyzed how tidally-induced orbital variation generated
in an inner planet may propagate to outer companions where tidal direct
effects with the central star are negligible. In the case of secular (non-
resonant) systems, this synergy affects mainly the eccentricities, although
the finite equilibrium eccentricity induced by the mutual gravitational
interactions also affects the rate of semi-major axis decay (e.g., Rodríguez
et al., 2011; Greenberg & Van Laerhoven, 2011). The case of resonant
configurations, particularly involving 3-planet systems, is more intricate.
Even though tidal evolution generates divergent migration and, thus, escape
from 2-planet resonances, Papaloizou (2015) found that in the case 3 planets
in a double resonance the escape occurs along a 3P-MMR commensurability.
Similar results were reported by Charalambous et al. (2018) under a wide range
of initial conditions and migration rates.

Figure 9: Tidal evolution of the K2-138 system in the mean-motion ratio plane.
Colored lines correspond to the triplets of Set#1. Initial positions (white
circles) correspond to the final stage of the capture simulation discussed in
the previous section, with semi-major axes scaled to the observed system. The
current mean-motion ratios of the planets are shown in black circles.
Analogous results are obtained for Set#2.

Figure 10: Time evolution of the mean-motion ratios of adjacent planets due to
tidal interactions with the central star. As with the previous plot, initial
conditions correspond to the final stage of the capture simulation of Set#1.
Observed values are highlighted with horizontal dashed lines.

The offset observed in several planetary systems with respect to 2P-MMRs (e.g.
Kepler-80, TOI-178) appear consistent with such an evolutionary route (e.g.,
MacDonald et al., 2016; Leleu et al., 2021), perhaps indicating that the
resonant structure of small-mass planetary systems could have a stronger
correlation to 3P-MMRs than to two-planet commensurabilities (Cerioni et al.,
2022). However, in no case was a capture reported in a first-order three-
planet resonance, nor was such a configuration suggested as an evolutionary
route due to tidal effects.

In order to check such a hypothesis for K2-138, we performed a series of
N-body simulations of Set#1 and Set#2 under the combined effects of mutual
gravitational perturbations and tides. The initial conditions were taken equal
to the end product of the capture simulations, scaling the semi-major axes to
better fit the observed system. Eccentricities and angles were untouched.

The tidal evolution was simulated following (Mignard, 1979). This model
considers tidal bulges raised on the planets as well as on the central star,
which will be slightly delayed from (or ahead of) their instantaneous
equipotential positions by a constant time-lag (CTL). We assumed a stellar
tidal parameter $Q^{\prime}_{0}=10^{6}$, while for the planets we chose
$Q^{\prime}_{i}=10^{2}$, regardless of the body’s mass or initial orbit.
Recall that in the CTL model these numerical values are defined at the onset
of the simulation and change during the orbital evolution of the system such
that $Q^{\prime}_{i}n_{i}=const$.

Figure 9 shows results for Set#1; analogous outcomes were obtained for Set#2.
Each plot shows the evolutionary tracks of mean-motion ratios associated to
adjacent planetary triplets. Initial values are depicted with white circles
while the evolution of the system in these planes is shown with colored lines;
arrows are drawn to help visualize the direction of the evolution. As
expected, tidal effects increase the 2-planet resonance offset while
preserving each triplet in its three-planet commensurability. This effect is
also noted in the case of the outer triplet where the guiding route is the
first-order 3P-MMR characterized by the index array $(2,-4,3)$. Finally, the
current mean-motion ratios of the planets are identified with black circles.

While not explicitly shown, semi-major axes decrease due to the tidal
interactions, and this effect is strongest for the innermost planet. This
differential migration pulls open the first planetary pair, which in turn
opens up every other pair in the resonance chain, much like the motion of a
pantograph-like mechanism. This is what is driving the evolution shown in
Figure 9. More details will be given in Section 6. We can thus see that the
observed values can in fact be obtained with tidal evolution, and that this
mechanism is effective even in the case of the outer pair $n_{f}/n_{g}$.

Additional proof is given in Figure 10 which shows the evolution of each mean-
motion ratio over time; current values are highlighted with horizontal dashed
lines. Not only is the observed orbital separation consistent with the tidal
evolution of the system, but the individual offsets of all planetary pairs are
attained at the same time.

In these simulations, and for the adopted values of the tidal parameters, the
current state of the system was reached at approximately $8$ Gyrs after the
dissipation of the primordial disk. Stellar age estimated from spectral
analysis yields $T_{\rm age}=2.3^{+0.44}_{-0.36}$ Gyrs (Lopez et al., 2019).
This requires a faster tidal evolution and thus a smaller value of
$Q^{\prime}_{i}$ for the planets; we find that tidal parameters of the order
of $Q^{\prime}_{i}\simeq 30$ would have led to the same configuration in
integration times consistent with $T_{\rm age}$. However, stellar ages are
notoriously difficult to evaluate with precision, and the actual value for
this system could be larger.

Figure 11: Sketch of two 3P-MMRs in mean-motion ratio space. Points correspond
to triplets locked in the respective three-body resonances. Right frames are
zoomed in with equal-scale axes, in order to better compare the slope of each
curve, which determines the ordering of the two-body resonant offsets:
$\Delta_{12}<\Delta_{23}$ for the blue curve and $\Delta_{12}>\Delta_{23}$ for
the red.

##  6 The Pantographic Effect

We repeated the previous simulations ignoring tidal interactions between the
star and all but the innermost planet ($m_{b}$); results were
indistinguishable, indicating that any change in the resonant offsets is
mainly driven by the gravitational interactions with the other bodies and not
by direct tidal forces. Thus, the current offsets of K2-138 may be used to
estimate the tidal parameter of the inner planet but not of the rest of the
system.

The observed increase of the resonance offsets with the semi-major axis is not
a consequence of tidal interactions but of interplay between planets in a
resonance chain. In order to understand the reason behind this effect, let us
concentrate on a generic triplet $m_{1}$, $m_{2}$, and $m_{3}$ trapped in a
double resonance but slightly away from the nominal commensurabilities:

$\frac{n_{1}}{n_{2}}=\frac{p_{1}+q_{1}}{p_{1}}+\Delta_{12}\hskip
11.38092pt;\hskip
11.38092pt\frac{n_{2}}{n_{3}}=\frac{p_{2}+q_{2}}{p_{2}}+\Delta_{23}$ (8)

where $\Delta_{12}$ and $\Delta_{23}$ are the offsets initially determined by
the capture process. Additionally, we assume that the three planets are also
trapped in a $(k_{1},k_{2},k_{3})$ three-body commensurability with no
noticeable offset, such that

$k_{1}n_{1}+k_{2}n_{2}+k_{3}n_{3}=0.$ (9)

where the order of the 3P-MMR is given by the sum of indexes
$k_{1}+k_{2}+k_{3}$. Dividing equation (9) by $n_{2}$ and replacing equation
(8), we get a condition for the offsets $\Delta_{12}$ and $\Delta_{23}$ that
must be satisfied in the face of 2P- and 3P-MMRs:

$\Delta_{23}=-k_{3}\ (k_{2}\ +\ k_{1}\ \frac{p_{1}+q_{1}}{p_{1}}\ +\ k_{1}\
\Delta_{12})^{-1}\ -\ \frac{p_{2}+q_{2}}{p_{2}}.$ (10)

Expression 10 shows that the offset of the outer pair, apart from being
precisely determined by the first, can also be either smaller or higher than
it, depending on the resonances involved. The same result can be deduced from
analyzing 3P-MMR slopes in the mean-motion ratio plane, as the offsets will be
the triplet’s coordinates with the origin placed at the 2P-MMR intersection.

Figure 11 shows two examples. A triplet trapped in the $(1,-4,3)$
commensurability close to the $n_{1}/n_{2}=2/1$ and $n_{2}/n_{3}=3/2$ will
have $\partial(n_{2}/n_{3})/\partial(n_{1}/n_{2})<1$ and thus the offset of
the outer pair will be smaller than that of the inner pair. The opposite
behavior is observed in the case of the three planets trapped in the
$(1,-3,2)$ 3P-MMR. In such a configuration the outer pair will suffer a larger
departure from the nominal two-planet resonance even if $m_{1}$ is the only
body driving the migration process.

Figure 12: Open black circles plot the observed two-planet resonance offsets
for the K2-138 system. Blue crosses show the predictions given by expression
(10) assuming the current $\Delta_{bc}$ as a starting value.

The same principle can be applied to any number of interacting triplets in a
long resonance chain, such as K2-138, where the $(2,-5,3)$ and $(2,-4,3)$
3P-MMRs have slopes $>1$ and the offsets increase with semi-major axis.
Furthermore, Expression 10 shows that only one given offset will determine the
rest, and that if any planetary pair is spread apart or brought together,
every other pair will simultaneously follow, while the triplets travel over
the curves of 3-planet resonances. This linked motion holds a resemblance to
the crossing points in a pantograph-like mechanism, and is what happens in
Figure 9 when tidal effects bring the innermost planet inwards, pulling open
the first planetary pair.

Figure 12 shows, in blue crosses, the observed offsets for each planetary
pair. As mentioned previously, the values grow with the distance to the star
with a surprisingly large jump in the case of $\Delta_{fg}$ due to the
pronounced slope of the (2,-4,3) resonance. The open black circles show the
values predicted from Equation (10) starting from the observed value of
$\Delta_{bc}$. The agreement with the real values is excellent in the case of
the first three links, while some discrepancy is noted in the case of the
outer pair. This difference is due to the accumulation of short-period
deviations in each successive pair of planets.

##  7 Discussion

In this paper we analyzed the current dynamical configuration of the complete
K2-138 system. While the resonance chain involving the five innermost planets
($m_{b}$ up to $m_{f}$) is an expected outcome of planetary migration, the
orbital characteristics of the outer planet $m_{g}$ are intriguing. Its large
distance from the rest and apparent non-resonant dynamics strike out in such
an otherwise compact system.

We analyzed two possible solutions to explain the dynamics of $m_{g}$. In the
first, a small-mass outermost planet could have migrated slower than the inner
system and become detached as the planets migrated closer to the star. This
requires a small mass for $m_{g}$ as well as fine-tuned initial separations
and disk parameters. While possible, this formation scenario appears unlikely.

A second, and more interesting, explanation arose when noting that the three
outer planets seem very close to the $(2,-4,3)$ first-order 3P-MMR. At first
glance this proximity did not appear relevant given that these resonances are
generally believed to be weak. Moreover, the mean-motion ratio
$n_{f}/n_{g}\simeq 3.3$ appears too far from any significant 2-planet
commensurability, thus making any correlation between the current separation
and a past resonance capture very unlikely. A closer look, however, showed
both beliefs unfounded. A more massive outer planet with $m_{g}\gtrsim
9m_{\oplus}$ initially located beyond the 3/1 MMR with $m_{f}$ could be
trapped in a double resonance $(n_{e}/n_{f}\ ,\ n_{f}/n_{g})=(3/2\ ,\ 3/1)$
leading to a libration of the critical angle of the 3P-MMR $(2,-4,3)$. This
outcome was found for a wide range of initial conditions and disk parameters.
The resulting configuration at the time of the dissipation of the primordial
disk was composed of all 6 planets in a stable resonance chain with small
offsets.

Tidal evolution of the system for times comparable with the system’s age
proved sufficient to convert the small resonance offsets to their present-day
values, at least for values of $Q^{\prime}$ of the order of
$Q^{\prime}_{b}\simeq 30$ for the innermost planet. Direct tidal interactions
between the other planets and the star proved negligible and thus their tidal
parameters could not be constrained. Surprisingly even the current offset of
the outer pair $n_{f}/n_{g}\simeq 3.3$ can be explained by tidal evolution
alone, indicating that perhaps other seemingly non-resonant systems may be
dynamically affected by two-planet or three-planet commensurabilities, even
more so than previously thought (e.g., Cerioni et al., 2022).

The propagation of resonance offsets among members of a resonance chain can
take several forms, depending on the 2P-MMRs and 3P-MMRs involved in each
link. As a result, offset-inducing mechanisms such as tidal effects will not
produce offsets that necessarily decrease with the distance from the star, as
previously thought. Offsets among chained planets will also be linked such
that the opening (contraction) of any single pair will simultaneously spread
apart (bring together) the rest, a behavior reminiscent of the
extension/contraction of a pantograph.

While at present it is not possible to prove that K2-138 forms a complete
6-planet resonant chain, including a first-order 3P-MMR, we believe the
dynamical evidence is sufficient for such a claim to be at least very
plausible.

The authors would like to express their gratitude to the computing facilities
of the Instituto de Astronomía Teórica y Experimental (IATE) and the
Universidad Nacional de Córdoba (UNC). This work was funded by research grants
from the Consejo Nacional de Investicaciones Científicas y Técnicas (CONICET)
and the Secretaría de Ciencia y Tecnología (SECYT/UNC).

## References

  * Batygin (2015) Batygin, K. 2015, MNRAS, 451, 2589, doi: 10.1093/mnras/stv1063
  * Beaugé & Cerioni (2022) Beaugé, C., & Cerioni, M. 2022, Celestial Mechanics and Dynamical Astronomy, 134, 57, doi: 10.1007/s10569-022-10113-4
  * Beaugé et al. (2003) Beaugé, C., Ferraz-Mello, S., & Michtchenko, T. A. 2003, ApJ, 593, 1124, doi: 10.1086/376568
  * Beaugé et al. (2006) Beaugé, C., Michtchenko, T. A., & Ferraz-Mello, S. 2006, MNRAS, 365, 1160, doi: 10.1111/j.1365-2966.2005.09779.x
  * Bulirsch & Stoer (1966) Bulirsch, R., & Stoer, J. 1966, Numerische Mathematik, 8, 1, doi: 10.1007/bf02165234
  * Cerioni et al. (2022) Cerioni, M., Beaugé, C., & Gallardo, T. 2022, MNRAS, 513, 541, doi: 10.1093/mnras/stac876
  * Charalambous et al. (2018) Charalambous, C., Martí, J. G., Beaugé, C., & Ramos, X. S. 2018, MNRAS, 477, 1414, doi: 10.1093/mnras/sty676
  * 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
  * Cresswell & Nelson (2006) Cresswell, P., & Nelson, R. P. 2006, A&A, 450, 833, doi: 10.1051/0004-6361:20054551
  * Delisle & Laskar (2014) Delisle, J. B., & Laskar, J. 2014, A&A, 570, L7, doi: 10.1051/0004-6361/201424227
  * Goldreich & Schlichting (2014) Goldreich, P., & Schlichting, H. E. 2014, AJ, 147, 32, doi: 10.1088/0004-6256/147/2/32
  * Goldreich & Tremaine (1980) Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425, doi: 10.1086/158356
  * Greenberg & Van Laerhoven (2011) Greenberg, R., & Van Laerhoven, C. 2011, ApJ, 733, 8, doi: 10.1088/0004-637X/733/1/8
  * Hardegree-Ullman et al. (2021) Hardegree-Ullman, K. K., Christiansen, J. L., Ciardi, D. R., et al. 2021, AJ, 161, 219, doi: 10.3847/1538-3881/abeab0
  * Henrard & Lemaitre (1983) Henrard, J., & Lemaitre, A. 1983, Celestial Mechanics, 30, 197, doi: 10.1007/BF01234306
  * Izidoro et al. (2017) Izidoro, A., Ogihara, M., Raymond, S. N., et al. 2017, MNRAS, 470, 1750, doi: 10.1093/mnras/stx1232
  * Lee & Peale (2002) Lee, M. H., & Peale, S. J. 2002, ApJ, 567, 596, doi: 10.1086/338504
  * Lega et al. (2013) Lega, E., Morbidelli, A., & Nesvorný, D. 2013, MNRAS, 431, 3494, doi: 10.1093/mnras/stt431
  * Leleu et al. (2021) Leleu, A., Alibert, Y., Hara, N. C., et al. 2021, A&A, 649, A26, doi: 10.1051/0004-6361/202039767
  * 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
  * Lithwick & Wu (2012) Lithwick, Y., & Wu, Y. 2012, ApJ, 756, L11, doi: 10.1088/2041-8205/756/1/L11
  * Lopez et al. (2019) Lopez, T. A., Barros, S. C. C., Santerne, A., et al. 2019, A&A, 631, A90, doi: 10.1051/0004-6361/201936267
  * MacDonald et al. (2016) MacDonald, M. G., Ragozzine, D., Fabrycky, D. C., et al. 2016, AJ, 152, 105, doi: 10.3847/0004-6256/152/4/105
  * Malhotra (1991) Malhotra, R. 1991, Icarus, 94, 399, doi: 10.1016/0019-1035(91)90237-N
  * Michtchenko et al. (2006) Michtchenko, T. A., Beaugé, C., & Ferraz-Mello, S. 2006, Celestial Mechanics and Dynamical Astronomy, 94, 411, doi: 10.1007/s10569-006-9009-x
  * Mignard (1979) Mignard, F. 1979, Moon and Planets, 20, 301, doi: 10.1007/BF00907581
  * Papaloizou (2015) Papaloizou, J. C. B. 2015, International Journal of Astrobiology, 14, 291, doi: 10.1017/S1473550414000147
  * Ramos et al. (2017) Ramos, X. S., Charalambous, C., Benítez-Llambay, P., & Beaugé, C. 2017, A&A, 602, A101, doi: 10.1051/0004-6361/201629642
  * Rodríguez et al. (2011) Rodríguez, A., Ferraz-Mello, S., Michtchenko, T. A., Beaugé, C., & Miloni, O. 2011, MNRAS, 415, 2349, doi: 10.1111/j.1365-2966.2011.18861.x
  * Tanaka et al. (2002) Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257, doi: 10.1086/324713
  * Tanaka & Ward (2004) Tanaka, H., & Ward, W. R. 2004, ApJ, 602, 388, doi: 10.1086/380992
  * Thommes et al. (2008) Thommes, E. W., Bryden, G., Wu, Y., & Rasio, F. A. 2008, ApJ, 675, 1538, doi: 10.1086/525244
  * Winn & Fabrycky (2015) Winn, J. N., & Fabrycky, D. C. 2015, ARA&A, 53, 409, doi: 10.1146/annurev-astro-082214-122246


# Collective states of active matter with stochastic reversals: emergent
chiral states and spontaneous current switching

Kristian Stølevik Olsen  Nordita, Royal Institute of Technology and Stockholm
University, Hannes Alfvéns väg 12, 23, SE-106 91 Stockholm, Sweden
Luiza Angheluta  The Njord Center, Department of Physics, University of Oslo,
Blindern, 0316 Oslo, Norway
Eirik Grude Flekkøy  PoreLab, the Njord Center, Department of Physics,
University of Oslo, Blindern, 0316 Oslo, Norway

###### Abstract

We study collective dynamical behavior of active particles with topological
interactions and directional reversals. Surprising phenomena are shown to
emerge as the interaction relaxation time is varied relative to the reversal
rate, such as spontaneous formation of collective chiral states due to phase
synchronization, and collective directional reversals in the presence of
confinement. The results have a direct relevance to modelling and
understanding collective reversals and synchronization phenomena in active
matter.

Active matter; Collective phenomena; Kuramoto model;

###### pacs:

Valid PACS appear here

††preprint: APS/123-QED

Collective behavior is ubiquitous in biological systems, ranging from
collective migration of cells and bacteria to flocks of birds [1].
Multistability is often present in realistic swarming behavior, where the
system undergoes a series of transitions between different collective states,
triggered by internal or external perturbations. Examples include the co-
existence of translational flocking and milling states observed for example in
fish [2] or bi-stable clockwise and anti-clockwise states of swarming of
locust in circular confinements or annuli [3, 4].

Understanding the role disorder plays in collective phenomena has become a
central question in modern active matter research. The case of external
quenched disorder has been studied rather thoroughly both for interacting and
non-interacting active particles [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
17, 18, 19, 20, 21, 22]. More interesting perhaps is the role played by
internal types of disorder and complexity, originating not from external
factors but rather from the dynamics of the particles themselves. In the case
of mobility of animals or microorganisms such intrinsic disorder corresponds
for example to behavioral variability in a population. Recent studies include
inhomogeneous chirality disorder throughout a population [23], speed
inhomogeneity [24] or the inclusion of dissenters [25] or particles with
contrarian tendencies [26].

Active matter in confinement has also gained more interest because it is
representative in biological systems and has potential microfluidic
applications. Surprising behaviors emergence under confinement, including
accumulation phenomena near boundaries [27, 28, 29, 30], non-trivial escape
dynamics [31, 32] and oscillating collective motion [33, 34]. In such bounded
spaces, the collective states typically conforms to the symmetry of the
confinement.

In this paper, we report on how topological alignment interactions and
individual reversals in systems of non-chiral active particles produce
striking collective behaviors, such as global chiral states in open spaces and
collective directional reversals under confinement. For this, we propose a
minimal model with two dimensionless parameters and study its mean-field limit
to predict the transition point from disordered to partially-ordered states.
We use extensive numerical simulations to explore the bistable dynamics of
partially-ordered states in both open space and under a confinement that
breaks rotational symmetry.

Figure 1: Sketch of the system under consideration. Self-propelled particles
align through a Kuramoto interaction and undergo Poissonian stochastic
reversals of their direction of motion (eg. blue particle). Interactions are
topological, and each particle interacts with a fixed number $N_{\text{int}}$
of randomly chosen particles ($N_{\text{int}}=4$ in sketch). The case of open
unbounded space (A) and confinement in the form of a channel (B) is
considered.

Crucial to the collective reversal phenomena observed in our model is the
competition between the time-scale associated with single-particle reversals
(SPR) and a relaxation time associated with the alignment interactions. Each
active particle has an orientation $\varphi_{i}$ which evolves according to a
Kuromoto equation, i.e.
$\dot{\varphi}_{i}=-\frac{1}{\tau_{a}}\left\langle\sin\left(\varphi_{i}-\varphi_{j}\right)\right\rangle_{\mathcal{I}_{i}}$
where $i=1,2,...,N$ is the particle index. This is similar to the phase
evolution of coupled oscillators, e.g. [35], with a homogeneous coupling
strength given here by the inverse of the relaxation time $\tau_{a}$. However,
there is fundamental difference here in that our ”oscillators” are moving in
space in the direction of their phase and with a constant speed. Alignment is
modelled using a simple topological interaction rule whereby particle $i$ is
given a random set $\mathcal{I}_{i}$ of $N_{\textbf{int}}<N$ interaction
partners. Such interactions correspond, for example, to non-local interactions
mediated by sensory input in animal swarms [36, 37]. In addition, we induce
single-particle reversals (SPR) as a stochastic switching modelled by a
Poisson process, where the probability of reversal in a time window $\Delta t$
is $p_{\text{rev}}=e^{-\gamma\Delta t}\gamma\Delta t$, with rate parameter
$\gamma$. Active matter systems with instantaneous Poissonian reversals have
gained more attention in recent times both for interacting and non-interacting
cases [38, 39, 40, 41]. Such abrupt reversals are also observed in Nature, for
example in some species of bacteria, including P. Putida and M. Xanthus [42,
43]. The reversal rate $\gamma$ sets the other timescale competing with
alignment timescale. The ratio between them provides us with a dimensionless
parameter $\mathcal{D}=(\gamma\tau_{a})^{-1}$, which we will refer to as the
delay number. Thus, in its dimensionless form, the evolution of the
orientation angle $\varphi_{i}(t)$ reads as

$\dot{\varphi}_{i}=-\mathcal{D}\left\langle\sin\left(\varphi_{i}-\varphi_{j}\right)\right\rangle_{\mathcal{I}_{i}}+\pi\sum_{\alpha}\delta(t-t_{\alpha}^{(i)}).$
(1)

For simplicity, we neglect the rotational diffusion term in the limit where
the timescale associated with the rotational meanderings is much larger than
the other relevant timescales.

Figure 2: A. Time-series of the order parameter $|Z|(t)$ plateaus with small
fluctuations around an asymptotic constant value $|Z_{0}|$. B) Asymptotic
values of $|Z|$ as a function of the delay number $\mathcal{D}$, shows both
for the mean-field limit $N_{\text{int}}=N-1$ and away from the mean-field
limit. Dashed line indicates the predicted (mean-field) critical value of the
delay number where order starts to form. Figures C-E shows center-of-mass
trajectories for $N=200$ , $N_{\text{int}}=0.3N$ (corresponding to the green
curve in B.), with delay numbers $\mathcal{D}=1,6,15$ respectively. Relative
scale is shown by the bar in the bottom of figures C-E). $\text{Pe}_{s}=1$ in
open space.  Figure 3: A,B) Behavior of the model for $\mathcal{D}=30$ and
$N_{\text{int}}=N/2$, $N=200$. Oscillations in the two (normalized )order
parameters $(\text{Re}(Z)/|Z|,\text{Im}(Z)/|Z|)$ signaling a non-zero
collective angular velocity of $\psi(t)$, similar to the results of Bonilla
_et al_ [26]. C) Time-series of the global direction of motion $\psi(t)$ for
$\mathcal{D}=7$, $N=200$, $N_{\text{int}}=0.3N$. Stochastic sign changes in
the global angular velocity $\dot{\psi}$ are observed.

Active particles have an intrinsic mobility, and move at a constant speed
$v_{0}$ in the direction dictated by the orientation angle $\varphi_{i}(t)$
which changes in time as in Eq. (1). In dimensionless units set by a scale
$\ell$, we have that $\dot{\mathbf{x}}_{i}=\text{Pe}_{s}\hat{e}_{i}(t)$ where
$\text{Pe}_{s}=v_{0}/(\gamma\ell)$ is the Péclet number of the microswimmer,
and $\hat{e}_{i}(t)=[\cos\phi_{i},\sin\phi_{i}]\equiv[e_{x}^{i},e_{y}^{i}]$.
Note that we have two model parameters $(\mathcal{D}$ and $\text{Pe}_{s})$
that control internal and external properties of the system. The delay number
says something about how quickly particles synchronize their orientation angle
compared to how quickly they undergo SPRs. The Péclet number measures the
persistence length in terms of some reference lengthscale. In open space, this
scale can be taken to be the persistence length itself, and the Péclet number
is just unity. Under confinement, a natural choice for $\ell$ is the system
size, in which case the Péclet number becomes a tunable parameter measuring
the strength of confinement relative to intrinsic dynamics. Our model is
formally equivalent to a Kuramoto model of coupled oscillators on a random
network driven by a Poissonian noise with a coupling to spatial dynamics. In
the absence of any boundaries or obstacles that alter the particle’s direction
of motion, this is a one-way coupling; the angular dynamics feeds into the
spatial dynamics, while the converse is not true. In the presence of
confinement, the behavior is more complex, since collisions with boundaries
may induce (approximately) instantaneous changes in the direction of motion.
We investigate the role of SPRs on the collective dynamics both under
confinement (fig. (1 A)) and in open space (fig. (1 B)).

To quantify the ordered flocking states of the interacting active particles,
we define two real order parameters $(r,\psi)$ or equivalently a complex order
parameter $Z$ analogous to that of the phase synchronization in the Kuramoto
model,

$Z(t)=r(t)e^{-i\psi(t)}=\frac{1}{n}\sum_{i=1}^{n}e^{-i\varphi_{i}(t)},$ (2)

which may alternatively be expressed as
$r^{2}=\overline{e_{x}}^{2}+\overline{e_{y}}^{2}$ and
$\psi=\tan^{-1}(\overline{e_{y}}/\overline{e_{x}})$, with the bar notation
denoting averages over all particles. We see that $r$ is the mean-squared
normalized velocity of the particles, while $\psi$ is their average direction
of motion. In the disordered state where $r=0$, Eq. (2) implies that the
orientation angles are uniformly distributed, while highly ordered states
$r\approx 1$ happens when $\varphi_{i}\approx\psi$.

Valuable insights can be gained in the mean-field approximation of our model,
corresponding to $N_{\text{int}}=n$. In the absence of confinement, the
orientation $\phi$ has a probability density $\rho(\phi,t)$ that follows the
Fokker-Plank equation

$\partial_{t}\rho(\phi)=\partial_{\phi}\left[r\mathcal{D}\rho(\phi)\sin(\phi-\psi)\right]-\rho(\phi)+\rho(\phi+\pi),$
(3)

where the first term is the same as in the mean-field Kuramoto model [44],
while the last two terms describe loss and gain terms correspond to Poissonian
reversals. The order parameter from Eq. (2) reduces to
$Z=\langle\exp(-i\phi)\rangle$, with brackets indicating average over the
probability distribution. Note that $Z$ is nothing but the first Fourier mode
of the probability density.

The transition from disordered to ordered states takes place at a critical
value of the delay number $\mathcal{D}$. In the mean-field limit, this
critical number is determined from the linear stability analysis of Eq. (3)
[44]. By linearizing Eq. (3) around the disordered state, i.e. $r=0$ and
$\rho(\phi)=1/(2\pi)$, and using the fact that the complex order parameter is
the first Fourier mode of the density, one can derive an evolution equation
for the order parameter as $\dot{r}/r=\mathcal{D}/2-2$. Hence for the order
parameter to grow, we need $\mathcal{D}>4\equiv\mathcal{D}_{c}$. In terms of
dimension-full quantities the transition takes place when the mean waiting
time between single-particle reversals is larger than four times the
relaxation time of the interactions. Fig. (2 A-B) shows the order parameter
$Z(t)$ as function of time (panel A) and the dependence of its asymptotic
value on the delay number both in the mean-field limit ($N_{\text{int}}=N$)
and for a more sparse interaction network ($N_{\text{int}}\ll N$). We see a
clear agreement with the predicted $\mathcal{D}_{c}=4$ even beyond the mean
field approximation for sufficiently large number of particles. As expected
with decreasing the number of particles, the transition is less abrupt, but
remains consistent with predictions.

Fig. (2 C-E) shows representative center-of-mass trajectories for different
$\mathcal{D}$, which below $\mathcal{D}_{c}$ results in random meanderings
(panel C). Above this critical value (panels D and E) chiral collective motion
develops. At intermediate values of $\mathcal{D}$, the chirality switches
handedness spontaneously, as can also be seen from the time series of $\psi$
in Fig. (3 C). See movie 1 in the Supplemental Information for an example of
this behavior. Fig. (3 A-B) shows the behavior for large $\mathcal{D}$, where
the system locks into a state of fixed chirality, corresponding to a flocking
state where the flocks direction of motion rotates with a fixed handedness. If
we interpret the collective flock as a quasiparticle, the dynamics of this is
reminiscent of that of a chiral active Brownian particle with stochastic
chirality switching [41]. It is important to note that the individual
particles do not have an imposed chirality, and that the collective chiral
behavior in the (partially) ordered state is an emergent property of the
system originating in the Poissonian noise and finite-time interactions.
Similar emergence of chiral macroscopic states have been observed in
Myxococcus xanthus bacterial colonies [45].

Figure 4: A) Channel order parameter $\cos\psi$ as a function of time, showing
clear signs of reversals. Here $N=500,N_{\text{int}}=0.2N,Pe=1,\mathcal{D}=7$,
resulting in an order parameter $r\approx 0.32$. B) Collective waiting time
density associated with a long timeseries as in A. C) Collective reversal rate
as a function of the dimensionless number $\mathcal{D}$. As $\mathcal{D}$
increases, corresponding to increasing interaction strengths, the collective
reversal rate decreases. For smaller values of $\mathcal{D}$ (and smaller
values of the order parameter $r$) the reversals are indistinguishable from
noise. D) In the limit of strong confinement ($\text{Pe}_{s}$ increasing) the
collective reversal rate decreases as a consequence of the boundary
interactions promoting ordered states.

The mean-field approximation again provides us a qualitative understanding how
the Poissonian noise and the alignment interaction influence the stability of
the ordered state. We now linearize at the short-time dynamics of the mean-
field orientation $\phi$ away from the perfectly ordered state, i.e.
$\phi=\psi+\delta\phi$. The evolution of its average is
$\partial_{t}\langle\phi\rangle=-rD\langle\sin(\phi-\psi)\rangle-\pi$, where
the last term comes from the Poissonian noise, such that the perturbation away
from the order state evolves as

$\partial_{t}\langle\delta\phi\rangle\approx-r\mathcal{D}\langle\delta\phi\rangle-\pi$
(4)

We see that the first term has a stabilizing contribution, while the second
term corresponding to Poissonian reversals however is destabilizing the
ordered state. If instead, we consider the case relevant when a single
particle reverses its direction of motion and moves opposite to the flock
$\phi=\psi+\pi+\delta\phi$, the linearized dynamics reads

$\partial_{t}\langle\delta\phi\rangle\approx
r\mathcal{D}\langle\delta\phi\rangle-\pi$ (5)

Small deviations are now destabilized both by alignment interactions and
Poissonian reversals.

This destabilizing effect is crucial to understanding the collective motion in
the system. If a particle randomly reverses its direction of motion, any small
perturbation causes either a clockwise or anti-clockwise drift towards the
global direction of motion. This would also affect the evolution of the global
orientation $\psi$. When several SPR can take place over the same time window
and strong correlations are present, one particle’s random angular drift in
orientation will induce quick adjustments of the others, resulting in a biased
macroscopic angular drift. In this regime, we expect $\psi$ to evolve similar
to that of dichotomous diffusion, where the angular velocity $\dot{\psi}$
behaves like a telegraph noise. This is indeed what is observed in Fig. (3 C).
In the limit of high $\mathcal{D}$, the system gets locked into a chiral state
with a fixed handedness (see Fig. (2 E) or Fig. (3 A)), and displays little
spatial dynamics. Interestingly, the reversals on the microscopic scale do not
lead to sudden collective reversals of the global direction of motion as one
might expect. We attribute this to the lack of a symmetry breaking in the
system; the system preserves its a rotational invariance. Rather than
undergoing spontaneous macroscopic reversals $\psi\to\psi+\pi$ the system
explores all the equivalent macrostates which can be related through the
rotational symmetry.

One way to break the rotational symmetry is to include a non-rotationally
invariant confinement. When such confinement is introduced, the story becomes
more subtle. We consider an infinite channel geometry that is periodic in $x$
and has a finite extent in the $y$-direction. The channel’s width is used to
set the relevant length scale $\ell$ which allows us to tune the Peclet number
away from unity. First, let us consider the expected ordered states. Because
of the confinement’s geometry, there can be no average motion in the vertical
direction, and hence $\overline{e_{y}}\approx 0$. This is of course not true
exactly at every time $t$, but, on a slightly coarser time scale, the
fluctuations average to zero. This immediately implies from Eq. (2) that
$r\sin\psi=0$. For disordered states ($r=0$), this is trivially satisfied,
while for states with some degree of order ($r>0$), this can only hold if the
global direction of motion satisfies $\psi=n\pi$, implying that
$\cos\psi=\overline{e_{x}}/r=\pm 1$. Hence in a channel confinement, it will
be useful to quantify the order in the system by the order parameters
$(r,\cos\psi)$, where the latter dictates left $(-1)$ or right-moving $(+1)$
collective motion. The rotational symmetry in open space where $\psi$ takes
any values $\psi\in(0,2\pi)$ is now broken, and the remaining ordered states
are now connected instead by a discrete symmetry. The system spontaneously
transitions between the left- and right-moving states, similar to what has
been observed in some systems of collective animal behavior confined to an
annulus [46, 3, 4]. Fig. (4 A) shows the order parameter $\cos\psi$ undergoing
reversals like a telegraphic noise. To gain further insights into the
statistical nature of these reversals, we study the statistics of the inter-
reversal time, or the waiting time between collective reversals. Fig. (4 B)
shows a clear exponential distribution of the inter-reversal times, which
points to the fact there is the collective reversals are not correlated with
one another. Furthermore, Fig. (4 C) shows that the characteristic reversal
rate decreases with increasing the delay number $\mathcal{D}$, hence this
macroscopic reversal behavior tends to fade away into noise for sufficiently
frequent SPRs. At the opposite end of the spectrum, when sufficiently strong
interactions make it difficult to destabilize the ordered state, such that
macroscopic reversals become extremely rare events. A nontrivial dependence is
also observed with $\text{Pe}_{s}$. For small values of the Peclet number the
system becomes effectively unconfined and it is not possible to identify
collective reversals between the two macroscopic states (left- and right-
moving flocks). As the Peclet number increases the particles more strongly
interact with the solid boundaries. Such interactions promote ordered states
with $\cos\psi=\pm 1$ and hence the flocks will be locked into their
collective direction of motion. As with the delay number, the strongest
reversal phenomena are observed at intermediate values (See Movie 2 in the
Supplemental Information for an example of collective directional reversal in
confinement).

In summary, we studied a minimal model of active particles with topological
alignment interactions and random single particle reversals and found
surprising emergent behavior both in open spaces and under symmetry-breaking
confinement. In open space, where the rotational symmetry is preserved, we
predict that there is a critical delay number $\mathcal{D}_{c}=4$ above which
the system transitions spontaneously from a disordered state to emergent
collective chiral states, resulting from the combined effect of stochastic
SPR’s and a finite timescale of interactions. When the active particles are
confined to a channel geometry, collective spontaneous switching between left-
and right-going currents are observed. The collective directional switching
resemble a telegraphic noise and has a typical rate that decays quickly with
both $\mathcal{D}$ and $\text{Pe}_{s}$. Hence we expect that such reversal
phenomena are not as easily observed in the traditional Vicsek model, since in
that case the alignment timescale is zero and hence $\mathcal{D}\to\infty$.
Similar phenomena are however expected to be present in biological systems
where interactions have a non-zero timescale, for example in systems with
sensorial delay or rotational inertia.

###### Acknowledgements.

K.S.O acknowledges support from the Nordita Fellowship program. L.A and E.G.F
acknowledge support from the Research Council of Norway through the Center of
Excellence funding scheme, Project No. 262644 (PoreLab).

## References

  * Marchetti  _et al._ [2013] M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Reviews of Modern Physics 85, 1143 (2013).
  * Tunstrøm  _et al._ [2013] K. Tunstrøm, Y. Katz, C. C. Ioannou, C. Huepe, M. J. Lutz, and I. D. Couzin, PLoS computational biology 9, e1002915 (2013).
  * Yates  _et al._ [2009] C. A. Yates, R. Erban, C. Escudero, I. D. Couzin, J. Buhl, I. G. Kevrekidis, P. K. Maini, and D. J. Sumpter, Proceedings of the National Academy of Sciences 106, 5464 (2009).
  * Buhl  _et al._ [2006] J. Buhl, D. J. Sumpter, I. D. Couzin, J. J. Hale, E. Despland, E. R. Miller, and S. J. Simpson, Science 312, 1402 (2006).
  * Martínez-Calvo  _et al._ [2021] A. Martínez-Calvo, C. Trenado-Yuste, and S. S. Datta, arXiv preprint arXiv:2108.07011 (2021).
  * Martinez  _et al._ [2018] R. Martinez, F. Alarcon, D. R. Rodriguez, J. L. Aragones, and C. Valeriani, The European Physical Journal E 41, 1 (2018).
  * Duan  _et al._ [2021] Y. Duan, B. Mahault, Y.-q. Ma, X.-q. Shi, and H. Chaté, Physical Review Letters 126, 178001 (2021).
  * Chardac  _et al._ [2021] A. Chardac, S. Shankar, M. C. Marchetti, and D. Bartolo, Proceedings of the National Academy of Sciences 118 (2021).
  * Reichhardt and Reichhardt [2021a] C. Reichhardt and C. Reichhardt, Physical Review E 103, 022602 (2021a).
  * Reichhardt and Reichhardt [2021b] C. Reichhardt and C. Reichhardt, Physical Review E 103, 062603 (2021b).
  * Martinez  _et al._ [2020] R. Martinez, F. Alarcon, J. L. Aragones, and C. Valeriani, Soft matter 16, 4739 (2020).
  * Morin  _et al._ [2017] A. Morin, N. Desreumaux, J.-B. Caussin, and D. Bartolo, Nature Physics 13, 63 (2017).
  * Olsen  _et al._ [2021] K. S. Olsen, L. Angheluta, and E. G. Flekkøy, Soft Matter 17, 2151 (2021).
  * Alonso-Matilla  _et al._ [2019] R. Alonso-Matilla, B. Chakrabarti, and D. Saintillan, Physical Review Fluids 4, 043101 (2019).
  * Peruani and Aranson [2018] F. Peruani and I. S. Aranson, Physical review letters 120, 238101 (2018).
  * Sándor  _et al._ [2017] C. Sándor, A. Libal, C. Reichhardt, and C. O. Reichhardt, Physical Review E 95, 032606 (2017).
  * Borba  _et al._ [2020] A. Borba, J. L. Domingos, E. Moraes, F. Potiguar, and W. Ferreira, Physical Review E 101, 022601 (2020).
  * Adhikary and Santra [2021] S. Adhikary and S. B. Santra, EPL (Europhysics Letters) (2021).
  * Chepizhko and Peruani [2015] O. Chepizhko and F. Peruani, The European Physical Journal Special Topics 224, 1287 (2015).
  * Chepizhko  _et al._ [2013] O. Chepizhko, E. G. Altmann, and F. Peruani, Physical review letters 110, 238101 (2013).
  * Pattanayak  _et al._ [2019] S. Pattanayak, R. Das, M. Kumar, and S. Mishra, The European Physical Journal E 42, 1 (2019).
  * Rahmani  _et al._ [2021] P. Rahmani, F. Peruani, and P. Romanczuk, Communications Physics 4, 1 (2021).
  * Ventejou  _et al._ [2021] B. Ventejou, H. Chaté, R. Montagne, and X.-q. Shi, arXiv preprint arXiv:2107.14106 (2021).
  * Pattanayak  _et al._ [2020] S. Pattanayak, J. P. Singh, M. Kumar, and S. Mishra, Physical Review E 101, 052602 (2020).
  * Yllanes  _et al._ [2017] D. Yllanes, M. Leoni, and M. Marchetti, New Journal of Physics 19, 103026 (2017).
  * Bonilla and Trenado [2019] L. Bonilla and C. Trenado, Physical Review E 99, 012612 (2019).
  * Fily  _et al._ [2014] Y. Fily, A. Baskaran, and M. F. Hagan, Soft matter 10, 5609 (2014).
  * Fazli and Naji [2021] Z. Fazli and A. Naji, Physical Review E 103, 022601 (2021).
  * Leoni  _et al._ [2020] M. Leoni, M. Paoluzzi, S. Eldeen, A. Estrada, L. Nguyen, M. Alexandrescu, K. Sherb, and W. W. Ahmed, Physical Review Research 2, 043299 (2020).
  * Caprini and Marconi [2018] L. Caprini and U. M. B. Marconi, Soft matter 14, 9044 (2018).
  * Olsen  _et al._ [2020] K. S. Olsen, L. Angheluta, and E. G. Flekkøy, Phys. Rev. Research 2, 043314 (2020).
  * Paoluzzi  _et al._ [2020] M. Paoluzzi, L. Angelani, and A. Puglisi, Phys. Rev. E 102, 042617 (2020).
  * Zhang  _et al._ [2020] B. Zhang, B. Hilton, C. Short, A. Souslov, and A. Snezhko, Physical Review Research 2, 043225 (2020).
  * Liu  _et al._ [2020] P. Liu, H. Zhu, Y. Zeng, G. Du, L. Ning, D. Wang, K. Chen, Y. Lu, N. Zheng, F. Ye, _et al._ , Proceedings of the National Academy of Sciences 117, 11901 (2020).
  * Strogatz [2000] S. H. Strogatz, Physica D: Nonlinear Phenomena 143, 1 (2000).
  * Ballerini  _et al._ [2008] M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, V. Lecomte, A. Orlandi, G. Parisi, A. Procaccini, _et al._ , Proceedings of the national academy of sciences 105, 1232 (2008).
  * Zumaya  _et al._ [2018] M. Zumaya, H. Larralde, and M. Aldana, Scientific reports 8, 1 (2018).
  * Abbaspour  _et al._ [2021] L. Abbaspour, A. Malek, S. Karpitschka, and S. Klumpp, Effects of direction reversals on patterns of active filaments (2021), arXiv:2112.09188 [cond-mat.soft] .
  * Mahault  _et al._ [2018] B. Mahault, X.-c. Jiang, E. Bertin, Y.-q. Ma, A. Patelli, X.-q. Shi, and H. Chaté, Physical review letters 120, 258002 (2018).
  * Santra  _et al._ [2021a] I. Santra, U. Basu, and S. Sabhapandit, arXiv preprint arXiv:2101.11327 (2021a).
  * Olsen [2021] K. S. Olsen, Phys. Rev. E 103, 052608 (2021).
  * Pohl  _et al._ [2017] O. Pohl, M. Hintsche, Z. Alirezaeizanjani, M. Seyrich, C. Beta, and H. Stark, PLoS computational biology 13, e1005329 (2017).
  * Santra  _et al._ [2021b] I. Santra, U. Basu, and S. Sabhapandit, Phys. Rev. E 104, L012601 (2021b).
  * Strogatz and Mirollo [1991] S. H. Strogatz and R. E. Mirollo, Journal of Statistical Physics 63, 613 (1991).
  * Balagam  _et al._ [2021] R. Balagam, P. Cao, G. P. Sah, Z. A. Zhang, D. Wall, and O. A. Igoshin, bioRxiv (2021).
  * Chen  _et al._ [2021] D. Chen, Y. Sun, G. Shao, W. Yu, H.-T. Zhang, and W. Lin, Royal Society open science 8, 210649 (2021).