# Planet-Planet Scattering and ZLK Migration – The Dynamical History of
HAT-P-11

Tiger Lu Department of Astronomy, Yale University, 52 Hillhouse, New Haven, CT
06511, USA  Qier An Department of Physics, University of California, Santa
Barbara, Santa Barbara, CA 93106, USA  Gongjie Li Center for Relativistic
Astrophysics, School of Physics, Georgia Institute of Technology, Atlanta GA
30332, USA  Sarah C. Millholland MIT Kavli Institute for Astrophysics and
Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139,
USA  G. Mirek Brandt Department of Physics, University of California, Santa
Barbara, Santa Barbara, CA 93106, USA  Timothy D. Brandt Space Telescope
Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA  Department
of Physics, University of California, Santa Barbara, Santa Barbara, CA 93106,
USA

###### Abstract

The two planets of the HAT-P-11 system represent fascinating dynamical puzzles
due to their significant eccentricities and orbital misalignments. In
particular, HAT-P-11 b is on a close-in orbit that tides should have
circularized well within the age of the system. Here we propose a two-step
dynamical process that can reproduce all intriguing aspects of the system. We
first invoke planet-planet scattering to generate significant eccentricities
and mutual inclinations between the planets. We then propose that this
misalignment initiated von-Zeipel-Lidov-Kozai cycles and high-eccentricity
migration that ultimately brought HAT-P-11 b to its present-day orbit. We find
that this scenario is fully consistent only when significant tidally-driven
radius inflation is accounted for during the tidal migration. We present a
suite of N-body simulations exploring each phase of evolution and show that
this scenario is consistent with all observational posteriors and the reported
age of the system.

celestial mechanics, planet-star interactions, numerical methods

††journal: ApJ††software:  numpy (Oliphant, 2006; van der Walt et al., 2011),
REBOUND (Rein & Liu, 2012), REBOUNDx (Tamayo et al., 2020), Jupyter
(https://jupyter.org/).

tablenum SIXtablenum

##  1 Introduction

The mid-K dwarf HAT-P-11 hosts an intriguing system that has been the subject
of much interest in recent years. The system has two known planets. The first,
HAT-P-11 b, is a close-in ($a_{\mathrm{b}}=0.525$ AU) eccentric
($e_{\mathrm{b}}=0.218$) super-Neptune ($m_{\mathrm{b}}=23.4\>\mathrm{M}_{E}$)
first identified by Bakos et al. (2010). Follow-up Rossiter-McLaughlin
(Rossiter, 1924; McLaughlin, 1924) analysis of the system by Hirano et al.
(2011) and Sanchis-Ojeda & Winn (2011) revealed that the orbit of HAT-P-11 b
is polar ($\psi_{\mathrm{Ab}}=106^{\circ}$), which makes it a member of the
potential population of perpendicular planets (Albrecht et al., 2021; Dong &
Foreman-Mackey, 2023; Siegel et al., 2023).

HAT-P-11 b presents an interesting dynamical puzzle. Planets as close-in as
HAT-P-11 b are expected to be in circular orbits due to tidal forces. The fact
that HAT-P-11 b is significantly eccentric, coupled with its unusual
perpendicular orbit, points to a dynamically hot history. Yee et al. (2018)
significantly advanced our understanding of the dynamical history of the
system with their discovery of the second planet in the system, HAT-P-11 c, an
eccentric ($e_{\mathrm{c}}\approx 0.6$) super-Jupiter ($m_{\mathrm{c}}\sin
i\approx 1.6\>\mathrm{M}_{J}$) on an $a_{\mathrm{c}}\approx 4.1$ AU orbit.
Their dynamical analysis of the system concluded that with this additional
companion, both the eccentricity and stellar obliquity of HAT-P-11 b could be
explained on the condition that HAT-P-11 c was also misaligned. At the time,
the orbital inclination of HAT-P-11 c was unconstrained.

The final piece of the puzzle was provided by An et. al (in prep, hereafter
Paper 1), who combined radial velocities (RVs), astrometric data from
Hipparcos and Gaia, and projected obliquity constraints to confirm that the
orbits of HAT-P-11 b and HAT-P-11 c are indeed misaligned, with a bimodal
inclination posterior for HAT-P-11 c at values of $i_{\mathrm{c}}\approx
34^{\circ}$ or $i_{\mathrm{c}}\approx 144^{\circ}$. The mutual inclination
between the orbits of HAT-P-11 c and HAT-P-11 b, and the stellar obliquity of
HAT-P-11 c, were both constrained to be $>$30∘ at $3\sigma$-equivalent
significance, i.e., with significant misalignments measured for both. Given
the constraints provided by Paper 1, all relevant orbital quantities to the
dynamics of the system are well-constrained and are presented in Table 1. A
schematic view of the orbits is shown in Figure 1. We are now in a position to
postulate the complete dynamical history of HAT-P-11.

Table 1: Parameters of the HAT-P-11 System  Parameter | Value | Reference
---|---|---
Host Star
Mass (M⊙) | $0.81^{+0.02}_{-0.03}$ | Bakos et al. (2010)
Radius (R⊙) | $0.75^{+0.02}_{-0.02}$ | Bakos et al. (2010)
Age (Gyr) | $6.5^{+5.9}_{-4.1}$ | Yee et al. (2018)
HAT-P-11 baaHirano et al. (2011) and Huber et al. (2017) give different values
for stellar obliquity and eccentricity, respectively.
Mass (M⊕) | $23.4^{+1.5}_{-1.5}$ | Yee et al. (2018)
Radius (R⊕) | $4.36^{+0.06}_{-0.06}$ | Huber et al. (2017)
a (AU) | $0.05254^{+0.00064}_{-0.00066}$ | Yee et al. (2018)
e | $0.218^{+0.034}_{-0.031}$ | Yee et al. (2018)
$\psi_{\mathrm{Ab}}(^{\circ})$ | $106^{+15}_{-14}$ | Sanchis-Ojeda & Winn (2011)
HAT-P-11 c
Mass (MJ) | $3.06^{+0.42}_{-0.42}$ | Paper 1
a (AU) | $4.192^{+0.07}_{-0.07}$ | Paper 1
e | $0.56^{+0.036}_{-0.035}$ | Paper 1
$\psi_{\mathrm{Ac}}(^{\circ})$ | 43 to 138bbRange enclosing 95% of the (highly non-Gaussian) posterior | Paper 1
$\psi_{\mathrm{bc}}(^{\circ})$ | 50 to 130bbRange enclosing 95% of the (highly non-Gaussian) posterior | Paper 1
Figure 1: Schematic view of one possible orbital architecture of the HAT-P-11
system, based on posteriors from Table 1. We label the eccentricities and
misalignments between each orbit and the stellar spin axis.

Both planets in the HAT-P-11 system represent dynamical puzzles, and the
dynamical history of the system has been explored by Petrovich et al. (2020);
Pu & Lai (2021) among others. In this work we propose a two-step formation
mechanism: we posit that both HAT-P-11 b and HAT-P-11 c initially formed in
circular aligned orbits, along with two other large planets. Through dynamical
instability and subsequent planet-planet scattering, both additional planets
were ejected from the system. As a result, HAT-P-11 b and HAT-P-11 c were left
on eccentric and misaligned orbits, which triggered von-Zeipel-Lidov-Kozai
(ZLK) migration in HAT-P-11 b. Through these mechanisms, all aspects of the
present-day configuration of the system are reproduced.

The paper is organized as follows. In Section 2 we show that planet-planet
scattering results in significant eccentricities and orbital misalignments for
both planets, and we reproduce the present-day orbit of HAT-P-11 c. In Section
3 we similarly reproduce the present close-in state of HAT-P-11 b using a
series of N-body ZLK simulations accounting for both tidal and thermally
driven radius evolution. We discuss implications of our study and conclude in
Section 4.

##  2 Planet-Planet Scattering

The high eccentricity of HAT-P-11 c, coupled with the misaligned stellar
obliquity reported in Paper 1, is the first dynamical puzzle. Naively, one
would expect planets to form in circular, aligned orbits, thanks to strong
eccentricity damping within the protoplanetary disk. Any deviation from this
paradigm is believed to be the imprint of a dynamically active history. A
number of well-understood mechanisms are known to excite the eccentricities
and inclinations of planets, including but not limited to the von-Zeipel-
Lidov-Kozai (ZLK) effect (von Zeipel, 1910; Lidov, 1962; Kozai, 1962) and
secular chaos (Wu & Lithwick, 2011). Both of these mechanisms require an
additional undetected companion in the system, something for which there is no
direct evidence. The existing radial velocity data are satisfactorily
explained by two planets with no residual trend, while the measured
astrometric acceleration is similarly consistent with the two known planets
(Paper 1). To evade these constraints, an additional planet would have to be
low-mass and/or very widely separated and hence be dynamically insignificant.
Planet-disk interactions during planetary migration can excite eccentricities
without the need for a companion, but this mechanism has been shown to be
limited to relatively low eccentricities, of order $e\sim 0.1$ (Duffell &
Chiang, 2015), and it also cannot generate stellar obliquity misalignment.

In the absence of a more distant perturber, the most promising pathway is
planet-planet scattering, where close encounters between pairs of planets lead
to strong gravitational interactions that lead to large eccentricities and
obliquities. While additional planets would be required for this mechanism,
violent scattering is prone to ejecting planets or scattering them to wide
orbits where they would be difficult to detect (e.g., Chatterjee et al., 2008;
Carrera et al., 2019; Frelikh et al., 2019). In this section, we explore the
generation of HAT-P-11 c’s significant eccentricity and orbital misalignment
through N-body simulations of planet-planet scattering.

###  2.1 Theoretical Background

A wealth of literature exists on the subject of planet-planet scattering.
Analytic treatment of this problem is difficult, so most analyses have been
performed with numerical simulations. For a comprehensive review, see Davies
et al. (2014). Here we summarize some of the most relevant key results for our
study.

The stability of systems involving only two planets is well-understood
analytically (Gladman, 1993). Studies such as Ford et al. (2001) and Petrovich
et al. (2014) have shown that two-planet systems can attain a wide range of
eccentricities but cannot reach high mutual inclinations. Three-body
scattering, on the other hand, can both excite eccentricity to arbitrarily
large values and raise inclination up to $90^{\circ}$ (Chatterjee et al.,
2008; Carrera et al., 2019). In the stellar obliquity posterior derived in
Paper 1, the $95\%$ confidence interval is given by $None$ to $None$; thus the
$\sim$40∘ lower limit requires three-body scattering to be invoked.

The stability of three-body systems cannot be obtained analytically. From
numerical studies, system stability is both highly chaotic and dependent on
mean-motion resonances (Marzari, 2014; Rath et al., 2022, e.g.,). However, to
first-order the time to dynamical instability grows logarithmically with
mutual separation (Marzari & Weidenschilling, 2002). The mutual separation is
commonly parameterized by the mutual Hill radius:

$R_{H}=\frac{a_{i}+a_{i-1}}{2}\left(\frac{m_{i}+m_{i-1}}{3M_{*}}\right)^{1/3},$
(1)

and $\Delta$ denotes the separation in unit of the mutual Hill radius,

$\Delta=\frac{a_{i}-a_{i-1}}{R_{H}}.$ (2)

The stability timescale grows as the separation between the planets increases
(e.g., Chambers et al., 1996; Pu & Wu, 2015; Tamayo et al., 2020).

###  2.2 N-Body Simulations

Scattering simulations are computationally expensive, and there is a large
parameter space that the primordial system could inhabit. To work around this,
we adopt a sequential approach.

We first use three-planet simulations, with HAT-P-11 c and two additional
bodies but without HAT-P-11 b, to show that planet-planet scattering can
exactly reproduce HAT-P-11 c’s present-day orbit with the other bodies either
ejected from the system or scattered to very wide orbits. This step assumes
that we may neglect HAT-P-11 b for the purposes of reproducing HAT-P-11 c’s
observed parameters. We test this assumption by running an additional suite of
scattering simulations with HAT-P-11 b present on an initial 0.3 AU orbit
(roughly five times its present-day semimajor axis). In this suite of
simulations with HAT-P-11 b, we do not seek to exactly reproduce HAT-P-11 c’s
orbit, given the extremely tight constraints on its semimajor axis and the
large chaotic parameter space introduced by adding another surviving planet to
the system. Rather, we content ourselves with showing that these simulations
qualitatively match the architecture expected from the three-planet scattering
simulations.

Our tiered approach to N-body simulations allows us to construct a statistical
picture of the likely configurations of HAT-P-11 system after violent
scattering from an initial four-planet configuration. We can then estimate the
fraction of the available parameter space that can produce the high mutual
inclinations and significant eccentricities that are observed for HAT-P-11 c,
along with a range of plausible initial conditions for HAT-P-11 b’s subsequent
orbital evolution. This evolution under the influence of ZLK oscillations and
tides is discussed in Section 3.

####  2.2.1 Three-Planet Simulations

We perform 1000 three-planet simulations using the IAS15 integrator (Rein &
Spiegel, 2014) in the REBOUND N-body integrator package (Rein & Liu, 2012). We
use the adaptive timestep criterion described in Pham et al. (2024). The setup
of our simulations is as follows: we consider three identical planets with
$m_{1}=m_{2}=m_{3}=3.06\>M_{\mathrm{J}}$, the best-fit value for HAT-P-11 c’s
mass. Three-body scattering tends to result in a single surviving planet, or
additional survivors on highly eccentric/marginally unbound orbits (e.g
Carrera et al., 2019). In both cases, the orbital energy of the non-innermost
planets is nearly zero. Conservation of energy thus demands

$\frac{1}{a_{f}}=\frac{1}{a_{1,i}}+\frac{1}{a_{2,i}}+\frac{1}{a_{3,i}}.$ (3)

We set $a_{f}=4.192$ AU, the present-day semi-major axis of HAT-P-11 c. In the
interest of minimizing computation time, we set $\Delta=3$ to induce rapid
instability and scattering in our systems. This is unrealistically compact for
primordial planetary systems, where separations would be expected to be around
5-10 mutual Hill radii (Raymond & Morbidelli, 2022). Changing the mutual
separation of the planets does not qualitatively affect the scattering
outcome, only the timescale until dynamical instability (Chatterjee et al.,
2008). The above criterion hence informs the initial semi-major axes of the
three planets: $a_{1,i}=8.838$ AU, $a_{2,i}=13.281$ AU, $a_{3,i}=19.959$ AU.
In each simulation, the other relevant orbital parameters for each planet are
chosen at random from the following uniform distributions:

$\begin{split}e&\sim\mathcal{U}(0.01,0.05)\\\
i&\sim\mathcal{U}(0^{\circ},2^{\circ})\\\
f,\omega,\Omega&\sim\mathcal{U}(0,2\pi)\\\ \end{split}$ (4)

Note that we do not track the motion of the stellar spin axis in our
simulations. Rather, we assume that the star’s initial spin axis points along
the $z$-axis and does not evolve significantly over the course of our
simulations. With this assumption, the orbital inclination $i$ can be used as
a proxy for the stellar obliquity $\psi_{\mathrm{Ac}}$. We account for
collisions via merging planets (conserving mass, volume and momentum but not
energy) using the reb_collision_resolve_merge collision module. We also remove
particles if they exceed a distance of $10^{4}$ AU from the origin, using the
exit_max_distance feature.

Figure 2 shows an example of one of our simulations. The system rapidly enters
a regime of dynamical instability, as expected from the initialization of the
planets at small mutual Hill radii. This epoch is characterized by sudden
perturbations in semimajor axis, eccentricity and inclination triggered by
close encounters. At around $5$ Myr, the first planet is completely ejected
from the system. At this point, both remaining planets in the system gain
significant eccentricities and inclinations. The system remains stable over
the remaining duration of the simulation. The final orbital elements of
HAT-P-11 c are $a_{\mathrm{c}}=4.24$ AU, $e_{\rm c}=0.58$, and
$\psi_{\mathrm{Ac}}=68.5^{\circ}$, while the orbital elements of the surviving
outer planet are $a_{\mathrm{outer}}=240$ AU, $e_{\mathrm{outer}}=0.73$, and
$\psi_{\mathrm{outer}}=$$$.

Figure 2: Example three-body scattering results for the origin of HAT-P-11 c.
The thick maroon line and the two thin blue lines represent the orbital
elements of HAT-P-11 c and the two ejected planets (EP), respectively. The
subplots show the semimajor axis, eccentricity, and stellar obliquity
evolution of the three bodies, from top to bottom. When a planet is fully
ejected from the system (heliocentric distance $>10^{4}$ AU) or becomes
unbound $(a<0)$, the line associated with the evolution of that planet is
terminated. The black dot on the right-hand side of each subplot denotes our
derived values for the present-day system, with the error bar representing the
$1\sigma$ posteriors for semimajor axis and eccentricity and the $2\sigma$
posterior for stellar obliquity.

In Figure 3 we report statistics from the entire ensemble of simulations. We
plot histograms of the distributions in $a_{\mathrm{c}}$, $e_{\rm c}$ and
$\psi_{\rm Ac}$ of the innermost surviving planet in each simulation, as well
as the present-day $1\sigma$ constraints ($2\sigma$ for $\psi_{\rm Ac}$) on
the orbit of HAT-P-11 c from Paper 1. The simulation results are as expected:
$a_{\mathrm{c}}$ is strongly peaked near the present-day value, and both
$e_{\rm c}$ and $i_{\mathrm{c}}$ show a wide range of values that easily
encompass the constraints. This shows that our system setup is capable of
consistently producing orbits similar to that of HAT-P-11 c.

Figure 3: Results from our ensemble of three-body scattering simulations. We
plot distributions in $a$, $e$ and $\psi_{\rm Ac}$ for the innermost planet in
our three-planet simulations. The black lines correspond to the best-fit
present-day values of HAT-P-11 c from Table 1, with the dotted lines
representing the $1\sigma$ error bars for $a_{\mathrm{c}}$ and $e_{\rm c}$ and
$2\sigma$ for $\psi_{\rm Ac}$. All three orbital element ranges are well-
represented in our sample, showing that this system can easily be produced
through planet-planet scattering without fine-tuning.

Of the $1000$ simulations we ran, the number with $1$, $2$ and $3$ surviving
planets were $100,878$ and $22$, respectively. The most common outcome by far
is a system with two surviving planets, which would imply an external
undiscovered planet in the HAT-P-11 system. This is not in tension with
observations, because the extraneous surviving planets are generally very far
out, and hence they would both be undetectable and play a negligible role in
the future dynamics of the system. In our simulations there are a total of
$908$ extraneous planets that survive. The $1\sigma$ semimajor axis
distribution of these planets is $a_{\mathrm{extra}}=80^{+150}_{-37}$ AU,
while the perihelion distribution of these
$a_{p,\mathrm{extra}}=32^{+20}_{-12}$ AU. There is only one simulation where
there is an orbit crossing between any of the extra planets and the surviving
HAT-P-11 c. We do not analyze the stability of the surviving systems in-depth;
an analytic criteria such as the one introduced by Hadden & Lithwick (2018)
could be used to this end. We conclude that planet-planet scattering in a
three planet system well-reproduces the unusual stellar obliquity and
eccentricity of HAT-P-11 c.

####  2.2.2 Four-Planet Simulations

We now consider the effect of planet-planet scattering on the orbit of
HAT-P-11 b by running a suite of four-planet scattering simulations. The setup
of these simulations is as follows. We perform 10,000 simulations using the
new hybrid integrator TRACE (Lu et al., 2024) in REBOUND. We initialize the
outer three bodies as in the previous section. However, we now initialize
HAT-P-11 b on an $a_{\mathrm{b}}=0.3$ AU orbit (roughly $5\times$ its present-
day orbit), with all other orbital elements randomized in the same fashion as
the three outer planets. The semimajor axis of HAT-P-11 b does not change
significantly for most of our simulations. We adopt a timestep of $0.022$
years, which is roughly $1/9$th of HAT-P-11 b’s initialized orbit. The mean
fractional energy error over our 10,000 simulations is $10^{-1.35}$.

Figure 4 shows statistics of orbital elements from our simulations. We
restrict our discussion to simulations in which HAT-P-11 b survives on an
orbit within $a_{\mathrm{b}}<4$ AU; $2449$ simulations satisfy this criterion.
Of these simulations, the vast majority ($2111$, or $86\%$) result in three-
planet systems comprised of HAT-P-11 b, HAT-P-11 c, and a very distant outer
perturber. This is consistent with our results from the previous subsection,
which shows that HAT-P-11 b does not play a significant role in the scattering
dynamics of the giant planets. $116$ simultions resulted in two-planet systems
and $222$ simulations resulted in four-planet systems. We plot the
eccentricities and inclinations of all planets corresponding to HAT-P-11 b,
HAT-P-11 c and extraneous outer planets in their respective simulations as a
function of semimajor axis. The semimajor axis $1\sigma$ distributions are
given by $a_{\mathrm{b}}=0.38_{-0.11}^{+0.15}$ AU,
$a_{\mathrm{c}}=3.28_{-0.20}^{+0.24}$ AU and
$a_{\mathrm{outers}}=44_{-19}^{+110}$ AU. We also see that HAT-P-11 b can be
scattered onto a wide range of inclinations and eccentricities. These
$1\sigma$ distributions are given by $e_{\mathrm{b}}=0.21^{+0.42}_{-0.17}$ and
$\psi_{\mathrm{Ab}}=34_{-21}^{+30}$ degrees.

Figure 4: Results from an ensemble of 10,000 four-planet scattering
simulations. We plot the eccentricities and inclinations of the planets
representing HAT-P-11 b (green), HAT-P-11 c (maroon) and external planets
(blue) as a function of semimajor axis. The distributions of HAT-P-11 c and
the outer planets are qualitatively consistent with our results from the
three-planet scattering simulations. We see that HAT-P-11 b may be scattered
to arbitrarily high eccentricities and very high inclinations.

We conclude that the planet-planet scattering which was likely responsible for
the eccentric and misaligned orbit of HAT-P-11 c also likely resulted in
HAT-P-11 b gaining an eccentric and misaligned orbit. Crucially, this creates
significant misalignment between the orbits of the two planets.

##  3 ZLK Migration

The orbit of HAT-P-11 b is an even bigger puzzle than planet c. Planets as
close-in as HAT-P-11 b would naively be expected to have perfectly circular
orbits, due to significant tidal effects acting quickly to circularize the
orbit. Indeed, the vast majority of hot Jupiters observed do have circular
orbits. In addition, its near-polar orbit with $\psi_{Ab}\sim 106^{\circ}$
(Stassun et al., 2017) is also unusual. In this section, we explore ZLK
oscillations (driven by HAT-P-11 c) and subsequent tidal circularization and
migration as an avenue to explain both of these aspects of HAT-P-11 b’s orbit.
We use initial conditions consistent with our four-body scattering
simulations.

###  3.1 ZLK Oscillations

Lidov (1962) and Kozai (1962) independently discovered a fascinating behavior
characteristic of hierarchical three-body systems. Consider such a system, and
denote the three bodies as the star, b, and c. If the mutual inclination
between the orbits of b and c are significantly misaligned
($39.2^{\circ}<\psi_{\mathrm{bc}}<140.8^{\circ}$), then the orbit of b will
experience high-amplitude coupled oscillations in eccentricity and inclination
in the absence of short-range forces such as general relativity. Commonly
referred to as the Kozai-Lidov effect, the initial discovery of this mechanism
by von Zeipel (1910) has been recently brought to light by Ito & Ohtsuka
(2019). We hence adopt the name von Zeipel-Lidov-Kozai, or ZLK effect.

For an in-depth review of the ZLK effect as well as its many applications, see
Naoz (2016). We provide in this paper a brief overview of the analytic
understanding of the ZLK effect necessary to interpret our results. The
simplest way to investigate the ZLK effect analytically is in the hierarchical
secular approximation, which averages over the mean motion of both b and c
(Ford et al., 2000). The validity of this approximation holds as long as
$a_{\mathrm{c}}\gg a_{\mathrm{b}}$. The initial parameter space of the
HAT-P-11 system in our simulations does not always satisfy this requirement
($a_{\rm b}/a_{\rm c}<0.1$), so we will use direct N-body simulations to
explore our system to obtain more accurate results. However, good first-order
understanding can still be gained from the secular approximation. In this
approximation and assuming b is a test particle, c’s orbit remains unchanged
throughout the system’s evolution, and b’s semimajor axis remains fixed. The
maximum amplitude of the eccentricity oscillations depends on the initial
mutual inclination: assuming the inner begins on a circular orbit,
$e_{\mathrm{b,max}}=\sqrt{1-(5/3)\cos^{2}\psi_{\mathrm{bc}}}$ (in the
quadrupole limit). The period of the ZLK oscillations is of order (e.g.,
Fabrycky & Tremaine, 2007):

$\tau_{\rm ZLK}=\frac{2P^{2}_{\rm c}}{3\pi P_{\rm
b}}\frac{m_{\mathrm{star}}+m_{\rm b}+m_{\rm c}}{m_{\rm c}}(1-e^{2}_{\rm
c})^{3/2}$ (5)

It turns out that the genesis of these effects is the apsidal precession of
b’s orbit, $\dot{\omega}_{\rm b}$, generated by the weak perturbations of the
outer. In the secular approximation, the apsidal precession rate can be
explicitly calculated. The orbit-averaged Hamiltonian of the system is (Ford
et al., 2000; Fabrycky & Tremaine, 2007):

$\begin{split}\langle\mathcal{H}\rangle=&-\frac{Gm_{\mathrm{star}}m_{\rm
b}}{2a_{\mathrm{b}}}-\frac{G(m_{\mathrm{star}}+m_{\rm b})m_{\rm
c}}{2a_{\mathrm{c}}}\\\ &-\frac{Gm_{\mathrm{star}}m_{\rm b}m_{\rm
c}}{m_{\mathrm{star}}+m_{\rm
b}}\frac{a_{\mathrm{b}}^{2}}{8a_{\mathrm{c}}^{3}(1-e_{\rm c}^{2})^{3/2}}\\\
&\times\big{[}2+3e_{\mathrm{b}}^{2}-\sin^{2}\psi_{\mathrm{bc}}(3+12e_{\mathrm{b}}^{2}-15e_{\mathrm{b}}^{2}\cos^{2}\omega_{\mathrm{b}})\big{]}\end{split}$
(6)

This equation can be solved for the apsidal precession rate; see Naoz et al.
(2013); Naoz (2016) for an in-depth derivation:

$\begin{split}\dot{\omega}_{\rm b,\mathrm{ZLK}}=6C_{2}\bigg{[}&\frac{1}{G_{\rm
b}}\big{(}4\cos^{2}\psi_{\mathrm{bc}}+(5\cos(2\omega_{\mathrm{b}})-1)\\\
&\times(1-e_{\mathrm{b}}^{2}-\cos^{2}\psi_{\mathrm{bc}})\big{)}\\\
&+\frac{\cos\psi_{\mathrm{bc}}}{G_{\rm
c}}\big{(}2+e_{\mathrm{b}}^{2}(3-5\cos(2\omega_{\mathrm{b}})\big{)}\bigg{]}\end{split}$
(7)

where

$C_{2}\equiv\frac{G^{4}}{16}\frac{(m_{\rm s}+m_{\rm b})^{7}}{(m_{\rm s}+m_{\rm
b}+m_{\rm c})^{3}}\frac{m_{\rm c}^{7}}{(m_{\rm s}m_{\rm b})^{3}}\frac{L_{\rm
b}^{4}}{L_{\rm c}^{3}G_{\rm c}^{3}}$ (8)

Note that this level of approximation is only strictly valid in the case of a
circular outer orbit and a test particle inner, neither of which are good
assumptions in our system. The next level of approximation, the octupole-
level, would correctly account for both of these complications. However, the
dynamics are chaotic in the octupole order (Li et al., 2014). As we aim only
to draw first-order analytical understanding from the secular model (to
eventually be confirmed with full N-body simulations), we work for now in the
quadrupole approximation for greater clarity into the most impactful
quantities in this problem.

ZLK oscillations will persist so long as the dominant source of precession in
$\omega_{\mathrm{b}}$ is due to perturbations from c. As $a_{\mathrm{b}}$
becomes smaller due to tidal migration, a host of other effects become
significant and may suppress the precession from c, and thus the ZLK
oscillations themselves. For instance, the precession due to general
relativity can be quantified (Eggleton & Kiseleva-Eggleton, 2001; Fabrycky &
Tremaine, 2007):

$\dot{\omega}_{\rm b,\mathrm{GR}}=\frac{3G^{3/2}(m_{\mathrm{star}}+m_{\rm
b})^{3/2}}{a_{\mathrm{b}}^{5/2}c^{2}(1-e_{\mathrm{b}}^{2})}$ (9)

where $c$ here is the speed of light, not to be confused with a quantity
associated with body c. Yee et al. (2018) ruled out the possibility of ZLK
oscillations being active in the present-day HAT-P-11 system, on the basis of
general relativity suppressing ZLK precession at HAT-P-11 b’s present-day
semimajor axis. Plugging in values associated with the present-day system, we
see that ZLK oscillations are expected to be quenched by GR precession, with
$\dot{\omega}_{\rm b,\mathrm{ZLK}}/\dot{\omega}_{\rm b,\mathrm{GR}}\sim 0.6$.
Hence, HAT-P-11 b’s large eccentricity cannot be attributed to it currently
being in a high-eccentricity phase of a ZLK oscillation. However, this does
not rule out the possibility of ZLK oscillations having occurred in the past.
In order for this to have been the case, HAT-P-11 b must have originated at a
larger semimajor axis (consistent with our understanding that sub-Neptune-like
planets are challenging to form in-situ at the present-day semimajor axis of
HAT-P-11 b).

###  3.2 Tidal Friction

We turn to the process of ZLK migration, which has been invoked to explain the
anomalously high eccentricity and tight orbits of a number of hot Jupiters
(see Wu & Murray (2003), Beust et al. (2012) among others). The idea behind
ZLK migration is to combine the concepts of ZLK oscillations and tidal
evolution. During each high-eccentricity phase of the ZLK oscillation, tidal
friction becomes significant due to b’s close approach with the host star.

In this work, we will consider the prescription of equilibrium tides only
(Hut, 1981; Eggleton et al., 1998; Mardling & Lin, 2002). In this framework,
the strength of tidal friction is parameterized by the constant time lag
$\tau$, which represents the time lag between the tidal bulge of a body and
the line of centers connecting to the tidal perturber. $\tau$ may be related
(with caution, see Leconte et al. (2010); Lu et al. (2023)) to the commonly
used tidal quality factor $Q$ via

$\tau=(2nQ)^{-1}$ (10)

with $n$ the mean motion of the tidal perturber. We note that more complex and
nuanced tidal models are available. In particular, for the high-eccentricity
epochs, ZLK oscillations are characterized by dynamic tides may be a more
accurate prescription (see Vick et al. (2019) among others). However, no self-
consistent N-body package exists with these more complex prescriptions. Hence,
we restrict ourselves to the simpler equilibrium tide model in the present
work, though we encourage future works to explore the effect of more complex
tidal models.

In the equilibrium tide model, bodies are no longer point particles but rather
are endowed with structure. The two improvements on the point particle model
are the effect of a misaligned tidal bulge and the rotational flattening of
each body. Both are significant to the dynamics of ZLK migration, not only
because they act to tighten the orbit of b over time, but also because they
both introduce additional precession terms. For our system, tidal precession
is expected to be more relevant, so we will focus on that aspect. The
precession rate associated with tides is given by (Fabrycky & Tremaine, 2007)

$\begin{split}\dot{\omega}_{\mathrm{Tide}}=&\frac{15\sqrt{G(m_{\mathrm{star}}+m_{\rm
b})}}{8a_{\mathrm{b}}^{13/2}}\frac{8+12e_{\mathrm{b}}^{2}+e_{\mathrm{b}}^{4}}{(1-e_{\mathrm{b}}^{2})^{9/2}}\\\
&\times\frac{1}{2}\bigg{[}\frac{m_{\rm
b}}{m_{\mathrm{star}}}k_{\mathrm{star}}R_{\mathrm{star}}^{5}+\frac{m_{\mathrm{star}}}{m_{\rm
b}}k_{\rm b}R_{\rm b}^{5}\bigg{]}.\end{split}$ (11)

Figure 5: Comparison of precessions from the ZLK effect to general relativity
and tides in semimajor axis-eccentricity space, using HAT-P-11’s present-day
configuration. The maroon line represents where the precessions due to general
relativity and ZLK are exactly equal, while the gold lines represent tides
with HAT-P-11 b’s present day radius (solid) and tides with twice this radius
(dashed). Lines of constant angular momentum are shown in gray. HAT-P-11 b’s
present-day position is shown as a black cross. A planet experiencing ZLK-
induced high-eccentricity migration will move leftwards in phase space until
it encounters one of the colored curves, at which point ZLK oscillations are
quenched and it follows one of the contours of constant angular momentum to
zero eccentricity.

Figure 5 demonstrates the significance of including the effects of tidal
friction. In this figure we consider a system with a central object consistent
with HAT-P-11A and an outer perturber consistent with HAT-P-11 c. The
semimajor axis and eccentricity of HAT-P-11 b are varied, and we analyze where
in semimajor axis-eccentricity space ZLK oscillations are allowed. The
present-day orbit of HAT-P-11 b is also marked. The maroon and gold lines mark
where $|\dot{\omega}_{\mathrm{ZLK}}/\dot{\omega}_{\mathrm{GR}}|=1$ and
$|\dot{\omega}_{\mathrm{ZLK}}/\dot{\omega}_{\mathrm{Tide}}|=1$, respectively.
Towards the upper left (lower semimajor axis and higher eccentricity) of each
line, ZLK oscillations are quenched by the respective additional effect, while
towards the bottom right ZLK oscillations are allowed to be active. From
Equation (11), we see that $\dot{\omega}_{\mathrm{Tide}}$ scales strongly with
the eccentricity of b. As a result, while tides are less influential than GR
at low eccentricities, the opposite is true at higher eccentricities. This is
highly impactful for ZLK oscillations: a system where b is at a sufficiently
high semimajor axis will be able to sustain ZLK oscillations of extremely high
amplitude if only GR is considered, but if tides are accounted for as well
during a high-eccentricity epoch the ZLK effect will be quenched.

Consider now the evolution of a body under the influence of tidal effects
experiencing ZLK oscillations. If octopule level effects are accounted for
(Naoz et al., 2013; Naoz, 2016), the maximum eccentricity of the ZLK
oscillations is not constant; they change slowly over the octupole timescale.
In addition, during each high-eccentricity epoch tides will shrink the orbit.
The result is that over the course of ZLK oscillations, b’s orbit will move
towards the upper left corner of Figure 5. Inevitably, at some point there
will be a crossing (mostly likely with the gold line), and ZLK oscillations
will be quenched. Then b is marooned on a high-eccentricity orbit. This state
is only temporary, however, as tides continue to be significant. Now, tidal
dissipation acts efficiently to circularize and shrink the planet’s orbit,
until b settles into a close-in, circular orbit. Tidal circularization
conserves angular momentum $L_{\rm b}$, so b’s evolution through phase space
will follow the contour of constant $L_{\rm b}$ (plotted in grey) intersecting
the gold line until the planet reaches zero eccentricity.

Taking this together, the basic idea of our hypothesis is as follows. HAT-P-11
b initially formed at a significantly greater semimajor axis than we see
today. The violent scattering events described in the previous section
resulted in a significant orbital misalignment between HAT-P-11 b and HAT-P-11
c, which initiated ZLK oscillations in HAT-P-11 b’s orbit. Tidal dissipation
gradually tightened HAT-P-11 b’s orbit until tidal precession broke the ZLK
oscillations, at which point the highly eccentric HAT-P-11 b began the process
of tidal circularization. We propose that HAT-P-11 b is in this final state,
having not yet fully circularized. This is similar to the story proposed by
Beust et al. (2012) to explain the orbit of the similarly eccentric close-in
sub-Neptune GJ 436, but in their case the outer perturber is not well
constrained.

###  3.3 Radius Evolution

It turns out that story we just elucidated is incomplete, since simple ZLK
migration alone cannot explain HAT-P-11 b’s orbit. This is due to how
relatively close-in HAT-P-11 c is. From Equation (7), the precession rate
associated with the ZLK effect scales very strongly with $a_{\mathrm{c}}$.
Since $a_{\mathrm{c}}$ is small, the precession associated with ZLK
oscillations is strong, and thus it requires very high eccentricities for
tides to break the ZLK cycles (e.g., in comparison to GJ 436). Inspection of
Figure 5 provides a clearer picture. The contour that HAT-P-11 b presently
lies on does not intersect the solid golden line in the phase space shown.
While it does intersect the dark red line denoting suppression of ZLK
oscillations due to GR precession, this mechanism would require suppression of
ZLK oscillations at an eccentricity of just $\approx$0.5 at a semimajor axis
where tidal damping is slow, and is not something that we see in our
simulations. In order for HAT-P-11 b to have a constant $L$ path from tidal
quenching of ZLK oscillations to its current orbit, we would have to extend
Figure 5 well to the right, at which point HAT-P-11 b would participate in the
scattering events discussed in Section 2. In other words, the quenching of ZLK
by tidal effects in a reasonable initial parameter space cannot bring HAT-P-11
b to its present-day orbit via tidal circularization. We see that more
realistic contours of constant $L$ terminate significantly closer to the star,
sometimes by up to an order of magnitude. This tendency to overshoot HAT-P-11
b’s present-day orbit is insensitive to initial conditions and is confirmed
with N-body simulations. Thus, the present-day configuration cannot be
explained with traditional ZLK migration with fixed planet radius and
equilibrium tides alone.

However, the planet’s radius can be inflated due to tidal heating during ZLK
migration. Radius evolution of planets due to tidal heating is a degree of
freedom that is typically unexplored in studies of the ZLK effect. During ZLK
oscillations, HAT-P-11 b will reach very high eccentricities and hence
experience many close pericenter approaches with the host star. During these
close approaches, tidal heating becomes very significant and acts to inflate
the planet’s envelope. Radius inflation from tidal heating is a well-explored
concept: it has been invoked to explain the anomalously large radii of hot
Jupiters (Bodenheimer et al., 2001) and sub-Saturns (Millholland, 2019;
Millholland et al., 2020). These works show that the radius of a sub-Neptune
or sub-Saturn can reasonably be inflated over a factor of two via this
mechanism.

This effect may play a significant role in the dynamics of the system, as from
Equation (11) tidal precession scales as $\propto R_{\rm b}^{5}$. Indeed, we
see from the dotted golden line in Figure 5 that increasing the radius of
HAT-P-11 b by a factor of two breaks the ZLK oscillations significantly
earlier, and now there does exist a reasonable contour of constant $L_{\rm b}$
that can deliver HAT-P-11 b to its present-day orbit. Coupled dynamical and
thermal evolution of hot Jupiters were explored in the context of planet-
planet scattering by Rozner et al. (2022); Glanz et al. (2022), who found that
accounting for radius inflation significantly accelerated the formation and
destruction of hot Jupiters and enhanced the production of warm Jupiters.
Fully self-consistent coupled dynamical and thermal evolution in the context
of ZLK migration has not been explored in detail. Given the uncertainties in
our system, we employ a simplistic model for the radius inflation of HAT-P-11
b in this work, and we defer a more nuanced prescription for self-consistent
thermal and orbital evolution to future work.

Millholland et al. (2020) generated 10,000 MESA models, which modify and build
upon the publicly available MESA models for sub-Saturns developed by Chen &
Rogers (2016). These models spanned four principle parameters: mass $M$,
envelope mass fraction $f_{\mathrm{env}}$, flux from the host star $F$, and
internal luminosity $\mathcal{L}$111Note that Millholland et al. (2020)
parameterized the tidal luminosity with the factor
$\Gamma\equiv\log_{10}\left[\frac{Q^{\prime}(1+\cos^{2}\epsilon)}{\sin^{2}\epsilon}\right]$,
where $Q^{\prime}\equiv 3Q/2k_{2}$ is the reduced tidal quality factor and
$\epsilon$ is the planetary obliquity. This is because both Millholland et al.
(2020) and their progenitor study Millholland (2019) focused on obliquity
tides. However, we are agnostic to the specific source of tides in this study,
and as such require only the total luminosity deposited at the core-envelope
interface.. Each model was then evolved forward in time for $10$ Gyr, after
which a final radius $R$ value was reported. We fit a subset of the
Millholland et al. (2020) models appropriate to HAT-P-11 b in order to
construct a luminosity-radius relationship using a fourth-order
polylogarithmic function:

$\begin{split}\frac{R}{R_{\mathrm{E}}}=&A\left[\log_{10}\left(\frac{\mathcal{L}}{\mathcal{L}_{\odot}}\right)\right]^{4}+B\left[\log_{10}\left(\frac{\mathcal{L}}{\mathcal{L}_{\odot}}\right)\right]^{3}\\\
&+C\left[\log_{10}\left(\frac{\mathcal{L}}{\mathcal{L}_{\odot}}\right)\right]^{2}+D\left[\log_{10}\left(\frac{\mathcal{L}}{\mathcal{L}_{\odot}}\right)\right]\\\
&+E\end{split}$ (12)

where $A=5.23\times 10^{-4}$, $B=3.37\times 10^{-2}$, $C=0.801$, $D=8.357$,
and $E=3.596\times 10^{1}$. Here the orbit-averaged tidal luminosity of the
planet is given (Mardling & Lin, 2002):

$\begin{split}\langle\mathcal{L}\rangle=&-\mu_{\rm b*}a_{\mathrm{b}}^{2}n_{\rm
b}\left(\frac{m_{*}}{m_{\rm b}}\right)\left(\frac{R_{\rm
b}}{a}\right)^{5}\left(\frac{3k_{2,b}}{2Q_{\rm b}}\right)\\\
&\times\bigg{[}\frac{1}{2}\left(\Omega_{e}^{2}h_{1}(e)+\Omega_{q}^{2}h_{2}(e)\right)+\Omega_{h}^{2}h_{3}(e)\\\
&-2n_{\rm b}\Omega_{h}h_{4}(e)+n^{2}h_{5}(e)\bigg{]},\end{split}$ (13)

where $\Omega_{e},\Omega_{h},\Omega_{q}$ are functions of the planetary spin
axis and present orbital configuration, while $h_{1},h_{2},h_{3},h_{4},h_{5}$
are functions of eccentricity. The exact expressions are given in the
Appendix. Note the dependence on both orbital parameters and the spin axis.
Since our code self-consistently tracks the spin and orbital evolution, this
expression correctly accounts for both eccentricity and obliquity tides.

We select from the models generated by Millholland et al. (2020) as follows.
Out of the 10,000 models they generated, we take a slice in mass
($20\>\mathrm{M}_{\mathrm{E}}<\mathrm{M}<25\>\mathrm{M}_{\mathrm{E}}$)
appropriate for the observational constraints on HAT-P-11 b. The value of
$f_{\mathrm{env}}$ is not at all constrained, but the present-day radius and
orbital parameters of HAT-P-11 b are well-constrained (Yee et al., 2018). It
is therefore vital to match these parameters at the end of our simulations.
There is a degeneracy between $f_{\mathrm{env}}$ and the value of the tidal
quality factor $Q$ in reproducing HAT-P-11 b’s radius at its present-day
orbit, since the observed radius could be equally well-described by a small
envelope mass with efficient dissipation (low $Q$) or by a large envelope mass
with weaker dissipation (high $Q$). We assume a fiducial value of $Q=10^{5}$,
and then use Equations (12) and (13) to match the present-day radius of
HAT-P-11 b given its present-day orbital parameters. There is some degree of
uncertainty in the precise values of many of the values in Equation (13) at
the end of our simulations, but we can make a number of simplifying
assumptions to obtain a ballpark estimate. Most of the uncertainties lie in
the values of $\Omega_{e}$, $\Omega_{h}$ and $\Omega_{q}$, which are dot
products of the spin vector with the Runge-Lenz vector, orbit normal, and the
cross product of the two former, respectively. If we assume zero planetary
obliquity as expected given significant tidal dissipation (Su & Lai, 2022),
the spin axis will be aligned with the orbital angular momentum, and
$\Omega_{e},\Omega_{q}=0$. If we further assume that the planet rotates at the
pseudo-synchronous spin rate, $\Omega_{h}$ can be calculated. Plugging these
values and orbital parameters associated with HAT-P-11 b’s present-day orbit,
along with $Q=10^{5}$, gives the present-day tidal luminosity of HAT-P-11 b.
We select a range of $f_{\mathrm{env}}$ such that there is a good match with
this value. We ultimately choose $0.05\><f_{\mathrm{env}}<0.1$. With these
slices we consider 158 models. The models we selected, as well as our best-fit
line, are plotted in luminosity-radius space in Figure 6.

Figure 6: Radius-luminosity relationship derived from the models of
Millholland et al. (2020). In blue are 158 of their MESA model results using
the cuts described in the main text. The red line represents our best-fit to
the data, using a fourth-order polylogarithmmic function given in Equation
(12). The gold vertical line shows the tidal heating at HAT-P-11 b’s present-
day orbital configuration given $Q=10^{5}$. The black horizontal lines show
the observational constraints on HAT-P-11 b’s radius. Note that the gold and
red lines intersect with the black line, meaning HAT-P-11 b’s present-day
radius is reproduced.

We now use our fit to the MESA models to evolve HAT-P-11 b’s radius during the
simulations. In the interest of computation time and as our expression is
orbit-averaged, we do not recompute the planet’s radius every timestep.
Rather, radius evolution is computed once per orbit of HAT-P-11 b. Once per
orbit, the orbit-averaged luminosity is calculated and fed into our
luminosity-radius relation Eq. (12), and the radius of HAT-P-11 b is
correspondingly updated. We again emphasize that our prescription for the
thermal radius evolution of the planet is a first-order approximation, where
we assume the planet’s radius responds effectively instantaneously to the
tidal luminosity. A more realistic prescription would include a time-delayed
response to the tidal heating, and a fully self-consistent numerical
prescription could involve codes that couple the dynamical evolution of the
system to the planet’s thermal evolution (e.g. see Glanz et al. (2022), who
coupled outputs from REBOUND and MESA). While we do not use these nuanced
prescriptions in the present work, we argue that our first-order approximation
qualitatively achieves similar results. As our orbital code tracks tidal
luminosity self-consistently, our planet will ultimately be inflated to the
same degree, and the ZLK cycles will be broken at similar times. The
improvement which would likely produce the largest discrepancy between our
model is our choice of equilibrium tides, but we defer treatment of the more
realistic dynamical tide model to future work.

###  3.4 N-body simulations

We investigate the present-day orbit of HAT-P-11 b via ZLK migration through a
suite of N-body simulations using REBOUND. The setup of our numerical
simulations are as follows. We initialize HAT-P-11A and c in their present-day
orbital configurations. We select the best-fit values for the mass, semimajor
axis and eccentricity of HAT-P-11 c from Table 1. The stellar obliquity of
HAT-P-11 c is less constrained; we draw from a uniform distribution
$\psi_{\mathrm{c,i}}\in\\{$$,$$\\}$, where the lower bound is selected from
the $3\sigma$ confidence interval from the observational posteriors, and the
upper bound from the $2\sigma$ confidence interval of our four-body scattering
simulations. HAT-P-11 b is initialized with orbital elements randomly drawn
from uniform distributions informed from the $1\sigma$ distributions obtained
in our four-body scattering simulations: $a_{\mathrm{b}}=\\{0.27\
\mathrm{AU},0.52\ \mathrm{AU}\\}$, $e_{\mathrm{b}}=\\{0.04,0.627\\}$ and
$\psi_{\mathrm{b}}=\\{12.9^{\circ},64.7^{\circ}\\}$. We initialize the mass
with the best-fit value from Table 1, and the radius with Equation (12).

We also model the effects of general relativity and tides. We use REBOUNDx
(Tamayo et al., 2020) to model both. For GR precession, we use the gr
prescription (Anderson et al., 1975), appropriate for systems with a dominant
central mass. We consider self-consistent spin, tidal and orbital evolution
using the tides_spin effect (Lu et al., 2023), which provides the ability to
evolve the evolution of the stellar spin axis for consistent stellar obliquity
tracking as well. For planet b, we use fiducial parameters of $k_{2}=0.5$,
$I=0.25$ appropriate for a Uranus-like planet (Yoder, 1995; Lainey, 2016), and
we initialize its spin state as both pseudo-synchronized and aligned with the
orbit. Planet c is considered a point particle, and the effects of structure
are not accounted for.

Tides on the host star and planet c are not considered, but we set the value
of planet b’s tidal parameter to $\tau_{\rm b}=10^{-4}$ in the default REBOUND
units, where $G=1$. This is an unrealistically high value of $\tau$, chosen to
minimize computation time: direct N-body simulations with a realistic $\tau$
value proved to be computationally infeasible. This is a common practice in
numerical simulations (e.g Bolmont et al., 2015; Becker et al., 2020). Scaling
$\tau$ does not change the qualitative behavior of the system. Neither the ZLK
timescale nor the tidal precession rate, Equations (5) and (11) respectively,
depend on the value of the tidal parameter. The only dependence on $\tau$ is
the rate of orbital energy drained from tides. Thus, our scaled-up value of
$\tau$ serves only to drain more energy from the planet’s orbit every ZLK
cycle (therefore decreasing the semimajor axis more) and increase the rate at
which HAT-P-11 b circularizes once the ZLK oscillations are ultimately
quenched. These both only serve to linearly modify the timescale in which the
planet reaches the final stage of evolution.

We now briefly discuss what realistic values of $\tau$ may be. While the
values of tidal parameters in exoplanets are very poorly constrained,
estimates of Neptune’s $\tau\sim 10^{-7}$ are based on the orbit of its
satellites Proteus and Larissa (Zhang & Hamilton, 2007). Therefore, the tidal
timescale in the N-body simulations is a factor $\sim$103 longer than the
simulation timescale.

###  3.5 Simulation Results

Figure 7: Results from a population synthesis of ZLK simulations. We consider
$68$ simulations with radius evolution (orange) and $64$ with fixed radius
(blue). We compare distributions of semimajor axis, stellar obliquity and
mutual inclination. We see that the fixed radius model cannot generate
HAT-P-11 b’s present-day semi-major axis, but the model with tidally-driven
radius inflation can.

Figure 7 shows distributions of stellar obliquity and eccentricity for a
population synthesis of $100$ systems. We initialize each simulation as
described in the previous section and run until the ZLK cycles are quenched
and the orbit has circularized to the present-day eccentricity
$e_{\mathrm{b}}=0.218$. At this point, the simulation is halted. To illustrate
the effects of radius inflation, we also run the same initial condition but
fix the radius of HAT-P-11 b to its present-day value throughout the course of
the simulation. We allow each simulation to run for 10 Myr, which corresponds
to 10 Gyr in when time is rescaled appropriately. $68$ simulations with radius
inflation are able to quench ZLK oscillations and circularize in this
timeframe, with the remaining $32$ either never initiating ZLK oscillations or
being unable to quench them in time. For the simulations with fixed radius,
$64$ simulations quench and circularize in time. For the population synthesis
that considers radius inflation, we see that while the best-fit observational
values are not the most favored there is nonetheless a high proportion of
simulations that match the constraints, and we do not require fine-tuning to
reproduce our posteriors. The population synthesis without radius inflation
does not exhibit significant differences in the stellar obliquity or mutual
inclination distributions. However, as expected this model greatly
underpredicts the value the semimajor axis of HAT-P-11 b: all simulations
predict an orbit with $a_{\mathrm{b}}<0.045$ AU (consistent with the
implication of Figure 5).

Figure 8: Fiducial ZLK migration of HAT-P-11 b. From top to bottom, HAT-P-11
b’s semimajor axis, eccentricity, stellar obliquity and radius are shown with
the colored lines. The black dots represent HAT-P-11 b’s measured orbital
parameters, with the error bar indicating the $1\sigma$ confidence interval.
We posit that HAT-P-11 b is still in the tidal circularization stage after the
suppression of ZLK oscillations. There is excellent agreement between our
simulation results and observed values.

Figure 8 shows one of our simulations in detail. Table 2 lists the initial
orbital parameters of the simulation, as well as the final values of of
interest. During the first phase of the simulation, ZLK oscillations are
active and the planet experiences the characteristic coupled eccentricity and
inclination oscillations. At each high-eccentricity epoch, tides act on the
planet to lower its semimajor axis. The radius evolution of the planet is also
plotted, and we see that during each high-eccentricity epoch the planet
inflates to over twice its present-day radius. Note that the minimum
pericenter distance (${\sim}0.02\ \mathrm{AU}$) is well outside twice of the
tidal disruption radius ($r_{t}=R_{p}(M_{*}/M_{p})^{1/3}{\sim}0.008\
\mathrm{AU}$), so we do not include mass loss due to tidal stripping in our
model, as the effects are weak (Guillochon et al., 2011). Later, the ZLK
oscillations are damped during a high eccentricity epoch. Tides continue to
act to circularize the planet’s orbit. In the meantime, the orbit normal of
HAT-P-11 b precesses about the invariant plane of the system (essentially
HAT-P-11 c’s orbit normal), as expected, resulting in periodic stellar
obliquity oscillations. As expected, the stellar spin axis and the orbit of
HAT-P-11 c do not change significantly over the course of the integration.

We comment again on the interpretation of timescales, given our choice of
$\tau$ in order to speed up computations. Please refer to Section 3.4 for
justification regarding the choice of an enhanced dissipation rate. The
enhanced dissipation rate affects the phase when tidal dissipation dominates,
and thus we scale the timescales on the right-hand side of the panel to the
realistic values accordingly. Note that with a more realistic tidal time-lag,
we would expect to see more stellar obliquity cycles on the right-hand side of
the plot when tides dominate (since the orbital precession frequency does not
depend on tidal dissipation and shouldn’t be scaled). In addition, we would
expect to have more ZLK cycles on the left-hand side, since each cycle leads
to a weaker semi-major axis decay. However, the qualitative behavior of the
evolution is the same.

The age of the HAT-P-11 system is not well-constrained: Yee et al. (2018)
propose an age of $6.5^{+5.9}_{-4.1}$ Gyr, while Morton et al. (2016) present
an estimate of $2.69^{+2.88}_{-1.24}$ Gyr. Our results are consistent with the
system age. Note that the timescale of ZLK migration varied by around an order
of magnitude, depending on initial conditions.

Table 2: Relevant REBOUND simulation parameters

Parameter | Initial b | Final b | Initial c | Final c
---|---|---|---|---
$\mathrm{R}$ ($\mathrm{R_{Earth}}$) | 3.60 | 4.38 | N/A | N/A
$a$ (AU) | 0.496 | 0.0525 | 4.192 | 4.194
$e$ | 0.31 | 0.21 | 0.56 | 0.56
$\psi$ $(^{\circ})$ | 49.0 | 133.2 | 55.6 | 55.0

Note. — Relevant input parameters and output values of the REBOUND simulation
shown in Figure 8. The final values are the values at the time corresponding
to the black dot in Figure 8, which is where $a_{\mathrm{b}}$ is equal to the
best fit value. The orbit of HAT-P-11 c does not change significantly over the
course of the simulation.

The present-day observational constraints on HAT-P-11 b’s orbit are also shown
on Figure 8. Our simulation is able to match all values very well. We propose
that HAT-P-11 b is currently in this final stage of evolution and is still in
the process of tidal circularization. We conclude that ZLK migration coupled
with thermally driven radius evolution is capable of explaining the present-
day orbit of HAT-P-11 b remarkably well. Again, we emphasize that we did not
perform an exhaustive parameter space analysis of this system. Rather, our aim
is to qualitatively show that such a process is capable of reproducing
HAT-P-11 b’s orbit.

##  4 Conclusions

We have proposed a two-step process of planet-planet scattering followed by
subsequent ZLK migration to explain the unusual architecture of the HAT-P-11
system. A violent scattering history in which two former planets are ejected
from the system is a viable explanation for the high eccentricity and stellar
obliquity of HAT-P-11 c. Then, the ZLK mechanism follows naturally with large
mutual inclination between HAT-P-11 b and c. We find that traditional ZLK
migration with fixed planetary radius cannot reproduce the orbit of HAT-P-11
b, but accounting for thermally-driven radius inflation can. Using N-body
simulations, we have verified that this scenario is consistent with all
observational constraints as well as the estimated age of the system. The
HAT-P-11 system is thus an excellent case study of planet-planet scattering
and ZLK migration, both of which are common effects believed to sculpt the
architectures of many planetary systems.

HAT-P-11 b adds to a growing census of planets on polar orbits around their
host stars (Albrecht et al., 2021) – the exact nature of this peak in the
stellar obliquity distribution is not fully understood on a statistical level
(Siegel et al., 2023; Dong & Foreman-Mackey, 2023). ZLK migration triggered by
the presence of an external planetary perturber has been shown to be
theoretically capable of generating the peak of perpendicular planets
(Petrovich & Tremaine, 2016), and Vick et al. (2023) demonstrated that if the
inner planet starts on an initially misaligned orbit, then perpendicular
planets are preferentially produced. They attributed this primordial
misalignment to an inclined binary companion torquing the protoplanetary disk
(Spalding & Batygin, 2014; Zanazzi & Lai, 2018). We show that planet-planet
scattering can lead to similar ZLK initial conditions for single star systems,
and hence we demonstrate that the interplay between planet-planet scattering
and ZLK migration can easily produce perpendicular planets as well. While
stellar obliquity measurements of close-in planets are relatively common, the
same cannot be said of the longer-period giants which would constitute said
external perturbers (Rice et al., 2021). The possibility of more stellar
obliquity measurements of these longer-period planets in the near future will
allow more comprehensive analysis of systems characterized by planet-triggered
ZLK migration.

Our study demonstrates the importance of considering the evolution of the
physical characteristics of a planet during ZLK migration, in addition to the
orbital parameters. In the high-eccentricity epochs of ZLK migration, the
planet’s envelope can inflate due to tidal heating in the interior. This is
highly significant, affecting not only the timescale of ZLK migration but also
the equilibrium semimajor axis the planet settles at. This is particularly
important in the close companion ZLK cases, since the strength of the ZLK
perturbations is relatively much stronger. We show that radius inflation
unlocks areas of parameter space which otherwise would be impossible to access
through standard ZLK migration with fixed radius. HAT-P-11 is one of the first
candidates with two confirmed planets in which ZLK migration is hypothesized
to have occurred. A few other systems for which planet-planet ZLK may be
possible have begun to emerge (Beust et al., 2012; Petrovich et al., 2018;
Bardalez Gagliuffi et al., 2021), and as more such systems are discovered it
will be crucial to consider the evolution of the physical parameters of the
planet.

While our study highlights the importance of considering physical evolution of
the planet, our prescription is a first-order approximation. Given the
significant uncertainties in quantities such as internal composition and tidal
quality factor present in all exoplanetary systems we have not explored more
sophisticated models, and argue that our first-order approximation provides
sufficient qualitative insight. However, future work could more deeply
consider avenues such as 1) a more sophisticated tidal model – calculating
tidal heating and orbital evolution based on dynamical tides, which are more
accurate for high eccentricity orbits, 2) coupling tidal heating and radius
inflation during orbital evolution (e.g., implementing a lag time such that
the radius does not respond instantaneously to tidal forcing, and considering
the change in planetary structure when calculating tidal effects (Ogilvie,
2009)) and 3) considering the effects of atmospheric mass loss (Vissapragada
et al., 2022).

Finally, we note that while it may be tempting to place constraints on
quantities such as the tidal quality factor with this work, this must be done
with caution. The final state of HAT-P-11 b depends on the magnitude of tidal
forcing compared to the magnitude of the ZLK effect. In practice, this means
that it depends on a complex interplay of the initial semimajor axis, tidal
quality factor, tidal love number, initial mutual inclination, and internal
composition, none of which are well-constrained. To explore such a large
parameter space would be a herculean task, and we again stress that this study
does not claim to have done so. Rather, our work should be viewed more
qualitatively as a proof-of-concept that thermally-driven radius inflation is
highly significant towards the final products of planet-planet ZLK migration.

We thank Doug Lin, Yubo Su, Yasmeen Asali, Chloe Neufeld and Yurou (Nina) Liu
for insightful discussions.

## References

  * Albrecht et al. (2021) Albrecht, S. H., Marcussen, M. L., Winn, J. N., Dawson, R. I., & Knudstrup, E. 2021, ApJ, 916, L1, doi: 10.3847/2041-8213/ac0f03
  * Alexander (1973) Alexander, M. 1973, Astrophysics and Space Science, 23, 459
  * Anderson et al. (1975) Anderson, J. D., Esposito, P. B., Martin, W., Thornton, C. L., & Muhleman, D. O. 1975, ApJ, 200, 221, doi: 10.1086/153779
  * Bakos et al. (2010) Bakos, G. Á., Torres, G., Pál, A., et al. 2010, ApJ, 710, 1724, doi: 10.1088/0004-637X/710/2/1724
  * Bardalez Gagliuffi et al. (2021) Bardalez Gagliuffi, D. C., Faherty, J. K., Li, Y., et al. 2021, ApJ, 922, L43, doi: 10.3847/2041-8213/ac382c
  * Becker et al. (2020) Becker, J., Batygin, K., Fabrycky, D., et al. 2020, AJ, 160, 254, doi: 10.3847/1538-3881/abbad3
  * Beust et al. (2012) Beust, H., Bonfils, X., Montagnier, G., Delfosse, X., & Forveille, T. 2012, A&A, 545, A88, doi: 10.1051/0004-6361/201219183
  * Bodenheimer et al. (2001) Bodenheimer, P., Lin, D. N. C., & Mardling, R. A. 2001, ApJ, 548, 466, doi: 10.1086/318667
  * Bolmont et al. (2015) Bolmont, E., Raymond, S. N., Leconte, J., Hersant, F., & Correia, A. C. M. 2015, A&A, 583, A116, doi: 10.1051/0004-6361/201525909
  * Carrera et al. (2019) Carrera, D., Raymond, S. N., & Davies, M. B. 2019, Astronomy & Astrophysics, 629, L7
  * Chambers et al. (1996) Chambers, J. E., Wetherill, G. W., & Boss, A. P. 1996, Icarus, 119, 261, doi: 10.1006/icar.1996.0019
  * Chatterjee et al. (2008) Chatterjee, S., Ford, E. B., Matsumura, S., & Rasio, F. A. 2008, The Astrophysical Journal, 686, 580
  * Chen & Rogers (2016) Chen, H., & Rogers, L. A. 2016, ApJ, 831, 180, doi: 10.3847/0004-637X/831/2/180
  * Davies et al. (2014) Davies, M. B., Adams, F. C., Armitage, P., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 787–808, doi: 10.2458/azu_uapress_9780816531240-ch034
  * Dong & Foreman-Mackey (2023) Dong, J., & Foreman-Mackey, D. 2023, AJ, 166, 112, doi: 10.3847/1538-3881/ace105
  * Duffell & Chiang (2015) Duffell, P. C., & Chiang, E. 2015, ApJ, 812, 94, doi: 10.1088/0004-637X/812/2/94
  * Eggleton et al. (1998) Eggleton, P. P., Kiseleva, L. G., & Hut, P. 1998, The Astrophysical Journal, 499, 853
  * Eggleton & Kiseleva-Eggleton (2001) Eggleton, P. P., & Kiseleva-Eggleton, L. 2001, ApJ, 562, 1012, doi: 10.1086/323843
  * Fabrycky & Tremaine (2007) Fabrycky, D., & Tremaine, S. 2007, ApJ, 669, 1298, doi: 10.1086/521702
  * Ford et al. (2001) Ford, E. B., Havlickova, M., & Rasio, F. A. 2001, Icarus, 150, 303, doi: 10.1006/icar.2001.6588
  * Ford et al. (2000) Ford, E. B., Kozinsky, B., & Rasio, F. A. 2000, ApJ, 535, 385, doi: 10.1086/308815
  * Frelikh et al. (2019) Frelikh, R., Jang, H., Murray-Clay, R. A., & Petrovich, C. 2019, ApJ, 884, L47, doi: 10.3847/2041-8213/ab4a7b
  * Gladman (1993) Gladman, B. 1993, Icarus, 106, 247, doi: 10.1006/icar.1993.1169
  * Glanz et al. (2022) Glanz, H., Rozner, M., Perets, H. B., & Grishin, E. 2022, ApJ, 931, 11, doi: 10.3847/1538-4357/ac6807
  * Guillochon et al. (2011) Guillochon, J., Ramirez-Ruiz, E., & Lin, D. 2011, ApJ, 732, 74, doi: 10.1088/0004-637X/732/2/74
  * Hadden & Lithwick (2018) Hadden, S., & Lithwick, Y. 2018, AJ, 156, 95, doi: 10.3847/1538-3881/aad32c
  * Hirano et al. (2011) Hirano, T., Narita, N., Shporer, A., et al. 2011, PASJ, 63, 531, doi: 10.1093/pasj/63.sp2.S531
  * Huber et al. (2017) Huber, K. F., Czesla, S., & Schmitt, J. H. M. M. 2017, A&A, 597, A113, doi: 10.1051/0004-6361/201629699
  * Hut (1981) Hut, P. 1981, A&A, 99, 126
  * Ito & Ohtsuka (2019) Ito, T., & Ohtsuka, K. 2019, Monographs on Environment, Earth and Planets, 7, 1, doi: 10.5047/meep.2019.00701.0001
  * Kozai (1962) Kozai, Y. 1962, AJ, 67, 591, doi: 10.1086/108790
  * Lainey (2016) Lainey, V. 2016, Celestial Mechanics and Dynamical Astronomy, 126, 145, doi: 10.1007/s10569-016-9695-y
  * Leconte et al. (2010) Leconte, J., Chabrier, G., Baraffe, I., & Levrard, B. 2010, Astronomy & Astrophysics, 516, A64
  * Li et al. (2014) Li, G., Naoz, S., Holman, M., & Loeb, A. 2014, ApJ, 791, 86, doi: 10.1088/0004-637X/791/2/86
  * Lidov (1962) Lidov, M. L. 1962, Planet. Space Sci., 9, 719, doi: 10.1016/0032-0633(62)90129-0
  * Lu et al. (2024) Lu, T., Hernandez, D. M., & Rein, H. 2024, arXiv e-prints, arXiv:2405.03800, doi: 10.48550/arXiv.2405.03800
  * Lu et al. (2023) Lu, T., Rein, H., Tamayo, D., et al. 2023, The Astrophysical Journal, 948, 41, doi: 10.3847/1538-4357/acc06d
  * Mardling & Lin (2002) Mardling, R. A., & Lin, D. N. 2002, The Astrophysical Journal, 573, 829
  * Marzari (2014) Marzari, F. 2014, MNRAS, 442, 1110, doi: 10.1093/mnras/stu929
  * Marzari & Weidenschilling (2002) Marzari, F., & Weidenschilling, S. J. 2002, Icarus, 156, 570, doi: 10.1006/icar.2001.6786
  * McLaughlin (1924) McLaughlin, D. B. 1924, ApJ, 60, 22, doi: 10.1086/142826
  * Millholland (2019) Millholland, S. 2019, ApJ, 886, 72, doi: 10.3847/1538-4357/ab4c3f
  * Millholland et al. (2020) Millholland, S., Petigura, E., & Batygin, K. 2020, ApJ, 897, 7, doi: 10.3847/1538-4357/ab959c
  * 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
  * Naoz (2016) Naoz, S. 2016, ARA&A, 54, 441, doi: 10.1146/annurev-astro-081915-023315
  * Naoz et al. (2013) Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2013, MNRAS, 431, 2155, doi: 10.1093/mnras/stt302
  * Ogilvie (2009) Ogilvie, G. I. 2009, MNRAS, 396, 794, doi: 10.1111/j.1365-2966.2009.14814.x
  * Oliphant (2006) Oliphant, T. 2006, NumPy: A guide to NumPy, USA: Trelgol Publishing.  http://www.numpy.org/
  * Petrovich et al. (2020) Petrovich, C., Muñoz, D. J., Kratter, K. M., & Malhotra, R. 2020, ApJ, 902, L5, doi: 10.3847/2041-8213/abb952
  * Petrovich & Tremaine (2016) Petrovich, C., & Tremaine, S. 2016, ApJ, 829, 132, doi: 10.3847/0004-637X/829/2/132
  * Petrovich et al. (2014) Petrovich, C., Tremaine, S., & Rafikov, R. 2014, ApJ, 786, 101, doi: 10.1088/0004-637X/786/2/101
  * Petrovich et al. (2018) Petrovich, C., Wu, Y., & Ali-Dib, M. 2018, The Astronomical Journal, 157, 5, doi: 10.3847/1538-3881/aaeed9
  * Pham et al. (2024) Pham, D., Rein, H., & Spiegel, D. S. 2024, The Open Journal of Astrophysics, 7, 1, doi: 10.21105/astro.2401.02849
  * Pu & Lai (2021) Pu, B., & Lai, D. 2021, MNRAS, 508, 597, doi: 10.1093/mnras/stab2504
  * Pu & Wu (2015) Pu, B., & Wu, Y. 2015, ApJ, 807, 44, doi: 10.1088/0004-637X/807/1/44
  * Rath et al. (2022) Rath, J., Hadden, S., & Lithwick, Y. 2022, ApJ, 932, 61, doi: 10.3847/1538-4357/ac5f57
  * Raymond & Morbidelli (2022) Raymond, S. N., & Morbidelli, A. 2022, Demographics of Exoplanetary Systems: Lecture Notes of the 3rd Advanced School on Exoplanetary Science, 3
  * Rein & Liu (2012) Rein, H., & Liu, S.-F. 2012, A&A, 537, A128, doi: 10.1051/0004-6361/201118085
  * Rein & Liu (2012) Rein, H., & Liu, S. F. 2012, A&A, 537, A128, doi: 10.1051/0004-6361/201118085
  * Rein & Spiegel (2014) Rein, H., & Spiegel, D. S. 2014, Monthly Notices of the Royal Astronomical Society, 446, 1424, doi: 10.1093/mnras/stu2164
  * Rice et al. (2021) Rice, M., Wang, S., Howard, A. W., et al. 2021, AJ, 162, 182, doi: 10.3847/1538-3881/ac1f8f
  * Rossiter (1924) Rossiter, R. A. 1924, ApJ, 60, 15, doi: 10.1086/142825
  * Rozner et al. (2022) Rozner, M., Glanz, H., Perets, H. B., & Grishin, E. 2022, ApJ, 931, 10, doi: 10.3847/1538-4357/ac6808
  * Sanchis-Ojeda & Winn (2011) Sanchis-Ojeda, R., & Winn, J. N. 2011, ApJ, 743, 61, doi: 10.1088/0004-637X/743/1/61
  * Siegel et al. (2023) Siegel, J. C., Winn, J. N., & Albrecht, S. H. 2023, ApJ, 950, L2, doi: 10.3847/2041-8213/acd62f
  * Spalding & Batygin (2014) Spalding, C., & Batygin, K. 2014, ApJ, 790, 42, doi: 10.1088/0004-637X/790/1/42
  * Stassun et al. (2017) Stassun, K. G., Collins, K. A., & Gaudi, B. S. 2017, AJ, 153, 136, doi: 10.3847/1538-3881/aa5df3
  * Su & Lai (2022) Su, Y., & Lai, D. 2022, MNRAS, 509, 3301, doi: 10.1093/mnras/stab3172
  * Tamayo et al. (2020) Tamayo, D., Rein, H., Shi, P., & Hernandez, D. M. 2020, Monthly Notices of the Royal Astronomical Society, 491, 2885
  * Tamayo et al. (2020) Tamayo, D., Cranmer, M., Hadden, S., et al. 2020, Proceedings of the National Academy of Science, 117, 18194, doi: 10.1073/pnas.2001258117
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Vick et al. (2019) Vick, M., Lai, D., & Anderson, K. R. 2019, MNRAS, 484, 5645, doi: 10.1093/mnras/stz354
  * Vick et al. (2023) Vick, M., Su, Y., & Lai, D. 2023, ApJ, 943, L13, doi: 10.3847/2041-8213/acaea6
  * Vissapragada et al. (2022) Vissapragada, S., Knutson, H. A., dos Santos, L. A., Wang, L., & Dai, F. 2022, ApJ, 927, 96, doi: 10.3847/1538-4357/ac4e8a
  * von Zeipel (1910) von Zeipel, H. 1910, Astronomische Nachrichten, 183, 345, doi: 10.1002/asna.19091832202
  * Wu & Lithwick (2011) Wu, Y., & Lithwick, Y. 2011, ApJ, 735, 109, doi: 10.1088/0004-637X/735/2/109
  * Wu & Murray (2003) Wu, Y., & Murray, N. 2003, ApJ, 589, 605, doi: 10.1086/374598
  * Yee et al. (2018) Yee, S. W., Petigura, E. A., Fulton, B. J., et al. 2018, AJ, 155, 255, doi: 10.3847/1538-3881/aabfec
  * Yoder (1995) Yoder, C. F. 1995, Global Earth Physics: A Handbook of Physical Constants, 1, 1
  * Zanazzi & Lai (2018) Zanazzi, J. J., & Lai, D. 2018, MNRAS, 478, 835, doi: 10.1093/mnras/sty1075
  * Zhang & Hamilton (2007) Zhang, K., & Hamilton, D. P. 2007, Icarus, 188, 386, doi: 10.1016/j.icarus.2006.12.002

##  Appendix A Self-Consistent Spin, Tidal and Dynamical Equations of Motion

The ZLK Migration results in Section 3 account for spin and tidal evolution in
the equilibrium tide framework. The equations of motion used are reproduced
from Eggleton et al. (1998), but see also Alexander (1973); Hut (1981);
Mardling & Lin (2002) for a deeper review of equilibrium tide theory. Also see
Lu et al. (2023) for details on the specific implementation.

In addition to point-particle gravity, we consider the acceleration due to the
quadrupole distortion of body 1, which accounts for both spin distortion and
tidal perturbation from another body 2:

$\begin{split}\textbf{f}_{\text{QD}}^{(1,2)}=r_{1}^{5}k_{L,1}\left(1+\frac{m_{2}}{m_{1}}\right)\cdot\bigg{[}&\frac{5(\mathbf{\Omega}_{1}\cdot\mathbf{d})^{2}\mathbf{d}}{2d^{7}}-\frac{\Omega_{1}^{2}\mathbf{d}}{2d^{5}}-\frac{(\mathbf{\Omega}_{1}\cdot\mathbf{d})\mathbf{\Omega}_{1}}{d^{5}}-\frac{6Gm_{2}\mathbf{d}}{d^{8}}\bigg{]},\end{split}$
(A1)

with an analagous expression for the quadrupole distortion of body 2 due to
body 1. We also consider the acceleration due to the tidal damping of body 1:

$\begin{split}\textbf{f}_{\text{TF}}^{(1,2)}=&-\frac{9\sigma_{1}k_{L,1}^{2}r_{1}^{10}}{2d^{10}}\left(m_{2}+\frac{m_{2}^{2}}{m_{1}}\right)\cdot\left[3\mathbf{d}(\mathbf{d}\cdot\dot{\mathbf{d}})+(\mathbf{d}\times\dot{\mathbf{d}}-\mathbf{\Omega}_{1}d^{2})\times\mathbf{d}\right],\end{split}$
(A2)

again, with an analogous expression for the tidal damping of body 2.

##  Appendix B Orbital Parameter Definitions

In Equation (13) several terms relating to the precise orientation of the
orbit are mentioned but not defined. These are provided here. First, we define
the spin vector $\mathbf{\Omega}$, which parameterizes both the direction and
magnitude of the planet’s rotation. Next, we define two vectors
$\mathbf{e},\mathbf{h}$. These are the Runge-Lenz vector (pointing in the
direction of periastron, with magnitude equal to the eccentricity) and the
orbital angular momentum vector. We define a third vector
$\mathbf{q}\equiv\mathbf{e}\times\mathbf{h}$. We can thus define the following
projections of the spin axis:

$\Omega_{e}=\mathbf{\Omega}\cdot\mathbf{\hat{e}},\Omega_{h}=\mathbf{\Omega}\cdot\mathbf{\hat{h}},\Omega_{q}=\mathbf{\Omega}\cdot\mathbf{\hat{q}}$
(B1)

The functions of eccentricity $h_{i}(e)$ are also given:

$h_{1}(e)=\frac{1+(3/2)e^{2}+(1/8)e^{4}}{(1-e^{2})^{9/2}}$ (B2)

$h_{2}(e)=\frac{1+(9/2)e^{2}+(5/8)e^{4}}{(1-e^{2})^{9/2}}$ (B3)

$h_{3}(e)=\frac{1+3e^{2}+(3/8)e^{4}}{(1-e^{2})^{9/2}}$ (B4)

$h_{4}(e)=\frac{1+(15/2)e^{2}+(45/8)e^{4}+(5/16)e^{6}}{(1-e^{2})^{6}}$ (B5)

$h_{5}(e)=\frac{1+(31/2)e^{2}+(255/8)e^{4}+(185/16)e^{6}+(25/64)e^{8}}{(1-e^{2})^{15/2}}$
(B6)