# Social Phases in Human Movement

Yi Zhang  Department of Physics, University of Miami, Coral Gables, Florida
33146, USA  Debasish Sarker  Department of Physics, University of Miami, Coral
Gables, Florida 33146, USA  Samantha Mitsven  Department of Psychology,
University of Miami, Coral Gables, Florida 33146, USA  Lynn Perry  Department
of Psychology, University of Miami, Coral Gables, Florida 33146, USA  Daniel
Messinger  Department of Psychology, University of Miami, Coral Gables,
Florida 33146, USA  Chaoming Song<EMAIL_ADDRESS>Department of Physics,
University of Miami, Coral Gables, Florida 33146, USA

###### Abstract

Recent empirical studies found different thermodynamic phases for collective
motion in animals. However, such a thermodynamic description for human
movement remains unclear, mainly due to the limited resolution of existing
tracking technologies. In this Letter, we used a new ultra-wideband radio
frequency identification technology to collect high-resolution spatial-
temporal data for classroom movements. We measured the “effective
temperatures” for the underlying systems, and found two social phases, gas and
liquid-vapor coexistence, associated with different temperatures. We also
proposed a simple statistical physics model that can reproduce different
empirically observed phases.

Collective motion is a ubiquitous phenomenon in many biological and social
systems including bacteria [1], animals [2], and humans [3]. Over the past two
decades, substantial research has been devoted to the understanding of human
movement, thanks in part to the availability of large-scale human mobility
data collected by a broad range of new technologies including banknote
circulation [4], mobile-phone records [5, 6, 7], Bluetooth scanning [8] and
GPS tracking [9, 10]. Initial empirical observations suggest that human
movement at the population level satisfies a generalized continuous-time
random walk [4], while later work shows that individuals’ mobility patterns
are better captured by a more sophisticated exploration and preferential
return mechanism [6]. However, these empirical studies focus mainly on the
asymptotic behavior of moving individuals; little emphasis is placed on social
interaction among them.

In parallel, theoretical work on the social force model [11, 12, 13, 14]
suggested a potential liquid-vapor transition in human movements because of
social interaction [11, 15]. Moreover, recent empirical studies of animals use
synchronized cameras to reconstruct high-resolution moving trajectories for
swarming insects [16, 17], fish [18, 19], and birds [20]. This research
indicates the coexistence of a “condensed” and a dilute “vapor” phase in
equilibrium in insect swarms, calling attention to a thermodynamic description
of individual movements. However, to our best knowledge, an empirical
exploration of such macroscopic phases in human movement is missing in the
current literature, partly due to a lack of high-resolution records of
individuals’ trajectories during social events.

Existing tracking techniques such as GPS [9, 10] and mobile phone data [5, 6]
are too coarse for such a task because social interactions often occur within
1-2 meters [21], beyond the spatial resolution of these traditional
techniques. In addition, social events may take place indoors where GPS and
cell phone tracking are largely limited to track social contacts. Badge-based
technologies [16, 17, 8, 22] and Bluetooth scanning [8] allow detection of
social contacts when two individuals are within a certain distance. However,
badge-based technologies neither give an explicit description of the defined
social distance nor capture the exact locations where interactions take place.

In this Letter, we investigate human movement in classrooms using a commercial
tracking system based on ultra-wideband radio frequency identification (UWB-
RFID). The UWB-RFID allows us to measure the locations and orientations of
each individual with a high spatial-temporal resolution, and estimate the
“effective temperature” of classroom movement. By measuring local density
distribution, we discovered two different “social” phases at different
temperatures: (i) gas-like phase where individuals move freely without forming
small social groups, and (ii) liquid-vapor coexistence-like phase where
individuals form small social groups with free-moving individuals joining and
leaving these groups. We find that the effective temperature for the gas-like
social phase is significantly higher than the effective temperature of the
liquid-vapor coexistence phase. Finally, we propose a new social force model
inspired by the empirical data. Numerical simulations of our model
successfully reproduce the empirically observed social phases. Finally, we
show that physical quantities measured in our model agree qualitatively with
the experimental observations.

Data Collection. Using UWB-RFID technology, we collected human movement
datasets in classrooms in three schools. Receivers in the corners of each
classroom received UWB-RFID signals at $2\sim 4$ Hz, allowing continuous
measurement of individual’s locations with a $13$ cm spatial accuracy [23,
24]. D1) comprises records in $2018-2019$ academic year on the classroom
movements of two preschools in the US, each consisting of $10$ (D1A) and $17$
(D1B) children aged between $2.5\sim 4$ years, respectively. A total of $13$
three-hour observations in class D1A and $16$ observations in class D1B were
recorded. Vests with RFID tags over children’s hips were specially designed to
be worn by them during their regular classroom activities. Unifying the
spatial location information from tags on the left and right side of each
vest, the orientation information of each child was also extracted. D2)
comprises records of children aged around $5$ years in one kindergarten in
Germany in the fall $2014$. Children were equipped with a single tag on their
wrists, which allowed the UWB-RFID to trace their locations only. Three
observations were recorded in D2, each involving a one-hour free-play period
(see SM Section A for details).

Empirical measures. A preliminary visualization (Fig. 1a) of nine individuals
from class D1 during a continuous class time period shows that the
participants are not uniformly located in the classroom space, suggesting they
may form small social groups. The emergence of such groups is consistent over
different observations with individuals joining and leaving the groups. It
appears that the density of children varies both by physical location and over
time. Furthermore, individuals’ orientations, represented by their arrows,
appear to be correlated within groups.

Figure 1: Trajectories and distributions. (a) An illustration of 300 seconds
of the trajectories of nine individuals from the D1A, colors and arrows
identify the specific individuals and their orientations, respectively. (b)
Angular correlation function $C_{\theta}(r)$ for the dataset D1.
$C_{\theta}(r)>0$ and $<0$ correspond to ferromagnetic and antiferromagnetic
orientations, respectively. The Grey dashed line represents no correlation.
Note dataset D2 has no orientation measurement. (c) Radial distribution
function $g(r)$ for all three datasets D1-2.

To quantify orientational interaction among pairs, we compute the angular
correlation function [25]
$C_{\theta}(r)=\langle\cos{(\theta_{i}(0)-\theta_{j}(r))}\rangle$ that
measures the alignment of pair angles as a function of separation distance
$r$, averaged over all pairs. For instance, $C_{\theta}(r)=1$ indicates a
parallel alignment (side-by-side), i.e., $\theta_{i}=\theta_{j}$ whereas
$C_{\theta}(r)=-1$ represents an opposite (face-to-face) alignment as
$\theta_{i}=-\theta_{j}$. Figure 1b plots $C_{\theta}(r)$ for datasets D1A and
D1B, showing that $C_{\theta}(r)$ is positive when $r$ is small and becomes
negative at $r\approx 2\sim 4m$. When _r_ becomes large, $C_{\theta}(r)$
approaches zero as the correlation declines. This observation suggests that at
short distances individuals are more likely to be oriented parallel to each
other, like ferromagnets, due to the domination of side-by-side interactions.
On the other hand, at longer distances they tend to be oppositely oriented
like anti-ferromagnets due to the contribution of face-to-face interaction.

We next measure the radial distribution function (RDF), $g(r)$, the
probability of finding a pair of individuals at a separation distance $r$,
relative to random movements. $g(r)=1$ indicates no pair correlations, where
$g(r)>1$ and $g(r)<1$ correspond to attractive and repulsive pair interactions
respectively. Figure 1c plots RDF for all datasets, finding similar patterns
across different datasets. First, $g(r)$ is consistently peaked around
$r_{0}\approx 0.8m$, suggesting a typical social contact distance regardless
of experimental settings. Second, individuals are unlikely to appear at very
short distances around $r<0.3m$, due to physical constraints of body size.
Third, social interaction is persistent over a broad range of distances from
$r=0.3m$ through $2-3$$m$. On the other hand, the peak values of $g(r)$ vary
significantly for D1A-D1B, indicating variation in social interaction
strength. Not only are interactions weakest in D2, but its RDF also converges
to unity at long distances, a typical pattern for gas. In contrast, we find
$g(r)$ is lower than $1$ asymptotically for D1 datasets, a phenomenon often
observed in the liquid-vapor coexistence phase [16, 26, 27]. This finding
implies that D1 has a different macroscopic phase compared to the dataset D2,
evidenced by social clustering due to stronger interaction [17, 28].

Figure 2: Empirical dynamical measures. (a)-(c) Mean square displacement
$\langle\Delta r^{2}\rangle$. (d)-(f) Velocity autocorrelation function
$C_{v}(t)$. Dots represent experimental measures, whereas red curves represent
the least-squaring fitting of the theory (see SM Section B).

Effective temperature. To distinguish between different social phases, we
implement tools from statistical physics to explore a central thermodynamic
quantity, temperature $T$. The Stokes-Einstein relation $k_{B}T=D\zeta$
relates temperature in Brownian motion to the diffusion constant $D$ and
viscosity $\zeta$. For interactive systems, the Stokes-Einstein relation has
been generalized as the fluctuation-dissipation theorem [29, 30]. To find
diffusion constant $D$, we first measure the mean square displacement (MSD),
$\langle\Delta r^{2}\rangle$, that quantifies the variance of individual
movements over time $t$. For two-dimensional Brownian motion, the MSD scales
linearly with time as $\langle\Delta r^{2}\rangle=4Dt$, where the slope is
proportional to the diffusion constant $D$. Figure 2a-c plots averaged MSD for
datasets D1-D2, showing a linear dependency at a short time scale
$t\ll\tau_{s}$, where $\tau_{s}$ represents a typical time scale for social
interaction. We find that $\tau_{s}\approx 3s$ for datasets D1A-D1B, whereas
the MSD in D2 does not show this typical diffusion regime, due to experimental
limitation of temporal resolution (see SM Section B2). In the continuation of
$\langle\Delta r^{2}\rangle$, $t\gg\tau_{s}$; that is, individuals diffuse
sub-linearly, as $\langle\Delta r^{2}\rangle\sim t^{\alpha}$. The exponent
$\alpha\approx 0.3\sim 0.5$ (see SM Section B2) indicates a slowdown in
movement due to social interaction. We measure the diffusion constant $D$ for
D1-D2 (see SM Section B2), finding that the D2 dataset has the largest
diffusion constant.

In parallel, we also measure the velocity autocorrelation function (VACF),
$C_{v}(t)=\langle\vec{v}(t)\cdot\vec{v}(0)\rangle$ that quantifies the
temporal velocity correlation over a time interval $t$. For Brownian
particles, VACF satisfies an exponential decay [29, 31, 32],
$C_{v}(t)=e^{-t/\tau}$, where $\tau=m/\zeta$ represents the typical time scale
[29]. Figures 2d-f show VACF measured for all datasets. At short times, the
VACF follows exponential decay as Brownian particles. However, VACF becomes
negative for a longer time, implying attractive social interactions among
individuals. We find the empirical VACF can be fit by
$Ae^{-t/\tau}-Be^{-t/\tau^{\prime}}$, where former and later terms capture the
thermal diffusion and social interaction, respectively. Note that
$\tau^{\prime}$ diffs from $\tau_{s}$ measured in MSD but they are positively
correlated (see SM Section B1).

The Stokes-Einstein relation allows us to define an effective temperature
$k_{B}T\equiv mD/\tau$ for human movement. Although such an effective
temperature differs from the thermodynamic temperature of physical substances,
it quantifies the movement activity of classrooms. We measured the reduced
effective temperature, $\tilde{T}\equiv k_{B}T/m$ for both datasets, finding
$\tilde{T}=0.06J/kg$ and $0.10J/kg$ for D1A and D1B, respectively, while D2
has the highest effective temperature $\tilde{T}=1.47J/kg$, reflecting the
fact that D2 involves only free-play activities (see SM Table S2 for a
summary). This is also consistent with the RDF measurement in Figure 1c,
showing that dataset D2 has the weakest social interaction.

Local density distribution. The RDF function suggests that the higher
temperature dataset D2 falls into a gas-like social phase. On the other hand,
the lower temperature datasets D1A-B appear to be in the liquid-vapor
coexistence-like social phase. To distinguish these two different social
phases, we calculated the local density distribution (LDD) [33], $p(\rho)$,
that measures the spatial variation of the local number density of particles
(see SM Section B3). In the coexistence phase, $p(\rho)$ is bipeaked [15]
because the particles are inhomogeneously distributed in space where the low
and high-density peaks represent vapor and liquid phases, respectively. In
contrast, $p(\rho)$ is single-peaked around the average density, $\bar{\rho}$,
for the gas phase. Figure 3a-c, plot $p(\rho)$ for datasets D1-D2, showing
that $p(\rho)$ is mono-peaked for the high-temperature dataset (D2). In
contrast, we find that $p(\rho)$ in D1 agrees well with a bi-gamma
distribution, with a vapor density peak around $\rho_{v}=0.2\bar{\rho}\sim
0.4\bar{\rho}$, and a liquid density peak around $\rho_{l}=1.05\bar{\rho}\sim
1.20\bar{\rho}$. This result suggests a liquid-vapor coexistence phase in D1
and a gas phase in D2. While existing models suggest a similar liquid-vapor
phase transition [15], these transitions are often coupled with a change of
orientational order. We measured the magnetization order parameter
$m_{0}\equiv\frac{1}{N}\left|\left<\left(\sum_{i=1}^{N}\cos\theta_{i},\sum_{i=1}^{N}\sin\theta_{i}\right)\right>\right|$,
for the dataset D1, finding $m_{0}\approx 0.001$ and $0.002$ for D1A and D1B
respectively. This finding suggests that unlike the existing models [11, 15]
the liquid-vapor phase transition does not accompany by orientational order
changes, partly due to the mixture of ferro- and antiferromagnetic orders at
different distances.

Phase diagram. We illustrate the $\rho-T$ phase diagram for classroom movement
in Figure 3d. For temperature below the critical temperature $T_{c}$, the
system undergoes a vapor-liquid phase transition where different phases are
governed by the system density. The low-density regime ($\rho<\rho_{v}$)
corresponds to a vapor phase, whereas the high-density regime
($\rho>\rho_{l}$) represents a liquid. When density is between $\rho_{v}$ and
$\rho_{l}$, vapor and liquid coexist. The density difference
$\Delta\rho=\rho_{l}-\rho_{v}\sim(T_{c}-T)^{\beta}$ decreases as temperature
$T$ approaches $T_{c}$, where the critical exponent $\beta=1/2$ for mean field
theory [34] and $\beta\approx 0.35\pm 0.15$ experimentally for molecules [35,
33]. The distinction between vapor and liquid disappears when $T=T_{c}$, and
above $T_{c}$ the system is gas-like.

While social systems clearly differ from molecules, similar equations seem to
apply to our datasets. In particular, the bipeaked local density distributions
in school D1 suggest they exhibit a vapor-liquid coexistence phase.
Individuals in D1 are less mobile with a relatively low temperature. Most
individuals are trapped within small social groups, in analogy to droplets in
molecular systems. Occasionally, some active individuals escape from social
groups and move like a gas molecule. Social groups evolve dynamically with
individuals joining and leaving the group. To quantify the liquid-vapor phase
transition, we find $\Delta\rho_{D1}\approx 0.9$ for school D1A, and
$\Delta\rho_{D2}\approx 0.6$ for school D1B. Since our dataset is limited to
two data points in the coexistence phase, we can only roughly estimate the
critical temperature $\tilde{T}_{c}\approx 0.11\sim 0.12J/kg$ and
$\beta\approx 0.3\sim 0.6$. In contrast, school D2 falls into the gas phase,
where individuals move freely without forming social groups under a high-
temperature environment, in line with its mono-peaked local density
distribution.

Figure 3: Local density distribution measures of empirical data and phase
diagram. (a)-(c) Local density distribution for datasets D1A-B, respectively.
$x$ axis is the local density normalized by average density $\bar{\rho}$.
Scatter plots represent the empirical measures, and colored curves indicate
curve fitting. D1 shows a typical bipeak pattern, implying a vapor-liquid
coexistence phase. D2 demonstrates a typical mono-peaked local density
distribution, representing a single vapor phase. (d) Density-temperature phase
diagram. The grey dashed curves separate the different phase domains. The red
crosses are estimated by the local density distributions of D1.

Simple model. Unlike existing models [11, 15], our empirical observation
suggests that liquid-vapor phase transition is independent of the change of
orientational order. Inspired by the empirical observation, we model moving
individuals as small magnets with a $XY$-model like interaction [36],

$\displaystyle
U(\textbf{r},\theta)=\sum_{i<j}J(r_{ij})\cos{(\theta_{i}-\theta_{j})},$ (1)

where $r_{ij}\equiv|\textbf{r}_{i}-\textbf{r}_{j}|$, and $J(r)$ is the
coupling function that varies by distance. When $J(r)$ is the Heaviside step
function, Eq. (1) leads to a social force similar to the Vicsek model [11].
However, our empirical observation in Figure 1b suggests a more subtle
interaction with mixed ferromagnetic and antiferromagnetic effects. Therefore,
we use the Morse potential [37] for the coupling function $J(r)$, which
captures the empirical observation of ferromagnetic interaction at a short
extent and an antiferromagnetic interaction at longer distances (Fig. 1b). In
addition, we also used a soft sphere potential for each individual to take
their physical size into consideration (see SM Section C).

Equation (1) can apply both to passive and active particles [38]. To
investigate the impact of social interaction, we use a passive model for
simplicity, which allows us to apply Monte Carlo simulations. We tested our
model at different temperatures and density values. We found that the proposed
model reproduces radial distribution and angular correlation functions in
qualitative agreement with empirical measurements (see SM Section D1).
Moreover, for densities comparable to the observed empirical densities, we
found mono-peaked local density distributions for high-temperature values
corresponding to a gas phase (Fig. 4a), whereas bipeaked distributions for
lower temperatures corresponding to a liquid-vapor coexistence phase (Fig.
4b). The insets of Fig. 4a-b demonstrates snapshots of movement patterns for
these two phases.

Figure 4: Simulation results. (a)-(b) show the local density distributions for
relatively high and low temperatures with an average density ($0.25m^{-2}$),
respectively. (c) Maximum mean-absolute order parameter $|M|_{max}$ as a
function of effective temperature $\tilde{T}$ with a relatively high average
density ($20m^{-2}$). (d) Density-temperature phase diagram. Scatter points
are estimated by MC simulation. Black solid lines represent the boundaries
among phases. Insets in (a)-(c) are the snapshots of corresponding systems.

Similar to the empirical dataset, there are no orientational order changes
during the liquid-vapor phase transition in our model. It is interesting to
explore such a magnetization phase transition, similar to the celebrated
Berezinskii–Kosterlitz–Thouless transition observed in the classical
$XY$-model [39, 40] and existing social force models [11, 12, 13, 14]. Indeed,
we observe a similar phase transition that occurs at very high-density values
around $\bar{\rho}\gg 1$, about one order of magnitude larger than that
observed empirically. Nevertheless, the magnetization order in our model shows
a striped pattern instead of a ferromagnetic phase (see Inset of Fig. 4c).
This reflects the fact that our coupling function $J(r)$ is ferromagnetic in
short distances and antiferromagnetic at longer distances. Figure 4c plots the
stripped order parameter $m_{s}$ and susceptibility $\chi_{s}$ as functions of
temperature (see SM Section D3 for definitions), providing clear evidence of
such a phase transition. Note that this orientational phase transition occurs
independently of the vapor-liquid phase transition, implying that the
underlying physics of our model Eq. (1) is distinct from those of existing
models [11, 15].

Figure 4d summarizes the phase diagram for our model generated by numerical
simulations. For a relatively low density, the model reproduces the gas and
vapor-liquid coexistence phases similar to those observed in experiments (see
Fig. 3d). For a high density $\bar{\rho}\gg 1$ which is beyond the access of
the existing datasets, the model predicts an ordered phase consisted of a
striped pattern. Needless to say, this prediction requires further validation
from future experiments.

In conclusion, we utilized relatively high-resolution UWB-RFID technology to
track human movements. We measured RDFs and angular correlation functions for
two preschools and one kindergarten. Measurement of LDDs revealed two social
phases, i.e., a gas-like phase corresponding to a single-peaked LDD and a
liquid-vapor coexistence-like phase with a doubly peaked LDD. We measured the
MSDs and VACFs which allowed us to calculate the effective temperature of
classroom movement. Different social phases correspond to different effective
temperatures and population densities, suggesting a phase diagram of human
movements driven by social interaction. As noted by Ref [16], the effective
temperature measure described here could enrich our understanding of phase co-
existence in insects and other social animals.

Specifically, the liquid-vapor coexistence-like phase and the higher
temperature gas-like system reflect differences in the observations of
schools. The lower-temperature classrooms were acquired over a complete half-
day observation in which teachers interacted actively with children. This
appeared to reduce effective temperature allowing the formation of small
dynamically stable social groups that are similar to liquid droplets. Free-
moving individuals may join and leave these social groups with relatively high
mobility, similar to vapor in the coexistence phase. Similar coexistence has
been observed in swarm insects [16]. By contrast, the higher-temperature
classrooms recorded the free-play activities of kindergarten children in which
teachers were requested not to interact with children. Comparatively higher
velocities and lower RDF peak values observed in these classrooms suggested a
weaker social interaction, making the state of the system to be more gas-like.

In physical systems, individual properties change continuously and macroscopic
systems show emergent changes in their phases. Here we show that social
systems may exhibit parallel properties. To the best of our knowledge, our
work provides the first experimental evidence of the existence of multiple
phases in human movements, a hypothesis posited theoretically two decades ago
[15]. Future research is required to experimentally manipulate control
parameters such as temperature that are hypothesized to control these phase
transitions. Inspired by these empirical observations, we proposed a social
interactive model that captures the nature of human orientational interactions
at different distances. Our model not only reproduces the observed social
phases but also captures several physical quantities including radial
distribution and angular correlation functions, in qualitative agreement with
empirical measurements. The results shed light on the role of social
interactions in human movement, with potential implementations in fields
including behavioral science [41], biological science and epidemiology [42].

###### Acknowledgements.

This work was supported partially by the National Science Foundation under
Grants No. 2150830 and No. IBSS-1620294, the Institute of Education Sciences
under Grant No. R324A180203, and the National Institutes of Health under Grant
No. R01DC018542.

## References

  * Zhang  _et al._ [2010] H.-P. Zhang, A. Be’er, E.-L. Florin, and H. L. Swinney, Collective motion and density fluctuations in bacterial colonies, Proceedings of the National Academy of Sciences 107, 13626 (2010).
  * Parrish and Edelstein-Keshet [1999] J. K. Parrish and L. Edelstein-Keshet, Complexity, pattern, and evolutionary trade-offs in animal aggregation, Science 284, 99 (1999).
  * Vicsek and Zafeiris [2012] T. Vicsek and A. Zafeiris, Collective motion, Physics reports 517, 71 (2012).
  * Brockmann  _et al._ [2006] D. Brockmann, L. Hufnagel, and T. Geisel, The scaling laws of human travel, Nature 439, 462 (2006).
  * Gonzalez  _et al._ [2008] M. C. Gonzalez, C. A. Hidalgo, and A.-L. Barabasi, Understanding individual human mobility patterns, nature 453, 779 (2008).
  * Song  _et al._ [2010] C. Song, T. Koren, P. Wang, and A.-L. Barabási, Modelling the scaling properties of human mobility, Nature physics 6, 818 (2010).
  * Simini  _et al._ [2012] F. Simini, M. C. González, A. Maritan, and A.-L. Barabási, A universal model for mobility and migration patterns, Nature 484, 96 (2012).
  * Sekara  _et al._ [2016] V. Sekara, A. Stopczynski, and S. Lehmann, Fundamental structures of dynamic social networks, Proceedings of the national academy of sciences 113, 9977 (2016).
  * Firth  _et al._ [2020] J. A. Firth, J. Hellewell, P. Klepac, S. Kissler, A. J. Kucharski, and L. G. Spurgin, Using a real-world network to model localized covid-19 control strategies, Nature medicine 26, 1616 (2020).
  * Nagy  _et al._ [2010] M. Nagy, Z. Ákos, D. Biro, and T. Vicsek, Hierarchical group dynamics in pigeon flocks, Nature 464, 890 (2010).
  * Vicsek  _et al._ [1995] T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Novel type of phase transition in a system of self-driven particles, Physical review letters 75, 1226 (1995).
  * Helbing and Molnar [1995] D. Helbing and P. Molnar, Social force model for pedestrian dynamics, Physical review E 51, 4282 (1995).
  * Helbing  _et al._ [2000] D. Helbing, I. Farkas, and T. Vicsek, Simulating dynamical features of escape panic, Nature 407, 487 (2000).
  * Helbing [2001] D. Helbing, Traffic and related self-driven many-particle systems, Reviews of modern physics 73, 1067 (2001).
  * Grégoire and Chaté [2004] G. Grégoire and H. Chaté, Onset of collective and cohesive motion, Physical review letters 92, 025702 (2004).
  * Sinhuber and Ouellette [2017] M. Sinhuber and N. T. Ouellette, Phase coexistence in insect swarms, Physical review letters 119, 178003 (2017).
  * Attanasi  _et al._ [2014a] A. Attanasi, A. Cavagna, L. Del Castello, I. Giardina, S. Melillo, L. Parisi, O. Pohl, B. Rossaro, E. Shen, E. Silvestri, _et al._ , Collective behaviour without collective order in wild swarms of midges, PLoS computational biology 10, e1003697 (2014a).
  * Katz  _et al._ [2011] Y. Katz, K. Tunstrøm, C. C. Ioannou, C. Huepe, and I. D. Couzin, Inferring the structure and dynamics of interactions in schooling fish, Proceedings of the National Academy of Sciences 108, 18720 (2011).
  * Herbert-Read  _et al._ [2011] J. E. Herbert-Read, A. Perna, R. P. Mann, T. M. Schaerf, D. J. Sumpter, and A. J. Ward, Inferring the rules of interaction of shoaling fish, Proceedings of the National Academy of Sciences 108, 18726 (2011).
  * Bialek  _et al._ [2012] W. Bialek, A. Cavagna, I. Giardina, T. Mora, E. Silvestri, M. Viale, and A. M. Walczak, Statistical mechanics for natural flocks of birds, Proceedings of the National Academy of Sciences 109, 4786 (2012).
  * Messinger  _et al._ [2019] D. S. Messinger, E. B. Prince, M. Zheng, K. Martin, S. G. Mitsven, S. Huang, T. Stölzel, N. Johnson, U. Rudolph, L. K. Perry, _et al._ , Continuous measurement of dynamic classroom social interactions, International Journal of Behavioral Development 43, 263 (2019).
  * Stehlé  _et al._ [2011] J. Stehlé, N. Voirin, A. Barrat, C. Cattuto, L. Isella, J.-F. Pinton, M. Quaggiotto, W. Van den Broeck, C. Régis, B. Lina, _et al._ , High-resolution measurements of face-to-face contact patterns in a primary school, PloS one 6, e23176 (2011).
  * Matakos  _et al._ [2017] A. Matakos, E. Terzi, and P. Tsaparas, Data Mining and Knowledge Discovery 31, 1480 (2017).
  * Irvin  _et al._ [2018] D. W. Irvin, S. A. Crutchfield, C. R. Greenwood, W. D. Kearns, and J. Buzhardt, An automated approach to measuring child movement and location in the early childhood classroom, Behavior research methods 50, 890 (2018).
  * Nho and Landau [2002] K. Nho and D. Landau, Spin-dynamics simulations of the triangular antiferromagnetic xy model, Physical Review B 66, 174403 (2002).
  * Doye and Wales [1996] J. P. Doye and D. J. Wales, The effect of the range of the potential on the structure and stability of simple liquids: from clusters to bulk, from sodium to, Journal of Physics B: Atomic, Molecular and Optical Physics 29, 4859 (1996).
  * Mochizuki and Koga [2015] K. Mochizuki and K. Koga, Solid- liquid critical behavior of water in nanopores, Proceedings of the National Academy of Sciences 112, 8221 (2015).
  * Attanasi  _et al._ [2014b] A. Attanasi, A. Cavagna, L. Del Castello, I. Giardina, S. Melillo, L. Parisi, O. Pohl, B. Rossaro, E. Shen, E. Silvestri, _et al._ , Finite-size scaling as a way to probe near-criticality in natural swarms, Physical review letters 113, 238102 (2014b).
  * Kubo [1966] R. Kubo, The fluctuation-dissipation theorem, Reports on progress in physics 29, 255 (1966).
  * Bruna  _et al._ [2017] M. Bruna, S. J. Chapman, and M. Robinson, Diffusion of particles with short-range interactions, SIAM Journal on Applied Mathematics 77, 2294 (2017).
  * Leegwater [1991] J. A. Leegwater, Velocity autocorrelation function of lennard-jones fluids, The Journal of chemical physics 94, 7402 (1991).
  * Li and Raizen [2013] T. Li and M. G. Raizen, Brownian motion at short time scales, Annalen der Physik 525, 281 (2013).
  * Hansen and McDonald [1988] J.-P. Hansen and I. R. McDonald, Theory of simple liquids, Physics Today 41, 89 (1988).
  * Fisher [1964] M. E. Fisher, Correlation functions and the critical region of simple fluids, Journal of Mathematical Physics 5, 944 (1964).
  * Goodwin [1970] R. D. Goodwin, Estimation of critical constants tc, $\rho$c from the $\rho$ (t) and t ($\rho$) relations at coexistence, Journal of Research of the National Bureau of Standards. Section A, Physics and Chemistry 74, 221 (1970).
  * Loft and DeGrand [1987] R. Loft and T. A. DeGrand, Numerical simulation of dynamics in the xy model, Physical Review B 35, 8528 (1987).
  * Kaplan [2003] I. Kaplan, Handbook of molecular physics and quantum chemistry, Wilson, S., Ed 3 (2003).
  * Bechinger  _et al._ [2016] C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, Active particles in complex and crowded environments, Reviews of Modern Physics 88, 045006 (2016).
  * Stanley [1971] H. E. Stanley,  _Phase transitions and critical phenomena_ , Vol. 7 (Clarendon Press, Oxford, 1971).
  * Jos [2013] J. V. Jos,  _40 Years of Berezinskii-Kosterlitz-Thouless Theory_ (World Scientific, 2013).
  * Saragosa-Harris  _et al._ [2022] N. M. Saragosa-Harris, A. O. Cohen, T. R. Reneau, W. J. Villano, A. S. Heller, and C. A. Hartley, Real-world exploration increases across adolescence and relates to affect, risk taking, and social connectivity, Psychological Science 33, 1664 (2022).
  * Zhang  _et al._ [2022] Y. Zhang, Y. Tao, M.-L. Shyu, L. K. Perry, P. R. Warde, D. S. Messinger, and C. Song, Simulating covid19 transmission from observed movement, Scientific reports 12, 1 (2022).


# Dissipative Capture of Planets Into First-Order Mean-Motion Resonances

Konstantin Batygin  Division of Geological and Planetary Sciences California
Institute of Technology, Pasadena, CA 91125, USA  Antoine C. Petit
Laboratoire Lagrange, Université Cote d’Azur, CNRS, Observatoire de la Cote
d’Azur, Nice, France

###### Abstract

The emergence of orbital resonances among planets is a natural consequence of
the early dynamical evolution of planetary systems. While it is well-
established that convergent migration is necessary for mean-motion
commensurabilities to emerge, recent numerical experiments have shown that the
existing adiabatic theory of resonant capture provides an incomplete
description of the relevant physics, leading to an erroneous mass scaling in
the regime of strong dissipation. In this work, we develop a new model for
resonance capture that self-consistently accounts for migration and
circularization of planetary orbits, and derive an analytic criterion based
upon stability analysis that describes the conditions necessary for the
formation of mean-motion resonances. We subsequently test our results against
numerical simulations and find satisfactory agreement. Our results elucidate
the critical role played by adiabaticity and resonant stability in shaping the
orbital architectures of planetary systems during the nebular epoch, and
provide a valuable tool for understanding their primordial dynamical
evolution.

Orbital dynamics, Perturbation theory

##  1 Introduction

Orbital resonances facilitate long-term exchange of energy and angular
momentum within planetary systems, thereby playing a critical role in their
long-term evolution. The preference for orbital commensurability – first
quantified in a statistically rigorous manner by Dermott (1968a, b) – is a
well-known attribute of the solar system’s architecture, which is particularly
pronounced among the satellites of Jupiter, Saturn, and Uranus (see Murray &
Dermott 1999 and the references therein). Beyond the realm of the solar
system, resonant configurations can be found in appreciable proportion within
the census of giant and sub-Jovian exoplanets alike (e.g., Goździewski et al.
2016; Mills et al. 2016; Luger et al. 2017; Petit et al. 2020; Nesvorný et al.
2022; Dai et al. 2023). Intriguingly, the importance of resonant dynamics
likely goes well beyond the population of planetary systems that are presently
entrained in mean-motion commensurabilities. That is to say, the evolutionary
role played by transient mean-motion resonances is almost certainly more
significant than a superficial examination of the data may indicate. To this
end, numerous lines of evidence suggest that the outer solar system itself
originated in a compact multi-resonant configuration before becoming
temporarily unstable and eventually settling in its current state by way of
dynamical friction (Batygin & Brown, 2010; Nesvorný & Morbidelli, 2012). Such
a sequence of events may in fact constitute a relatively typical post-nebular
evolutionary path of planetary systems, and recent modeling has shown that
both the period ratio distribution, as well as the degree of intra-system
uniformity of short-period super-Earths, can be satisfactorily reproduced if
the majority of systems originate as resonant chains that subsequently relax
towards more widely-spaced orbits through dynamical instabilities (Izidoro et
al., 2017, 2021; Goldberg & Batygin, 2022a; Batygin & Morbidelli, 2023).

Despite their extant and inferred prevalence, mean motion resonances do not
arise as an innate byproduct of the planet formation process itself. Instead,
they are established as a consequence of orbital convergence facilitated by
dissipative effects (Goldreich, 1965; Henrard, 1982). Within protoplanetary
nebulae, this occurs naturally due to planet-disk interactions (i.e., type-I
migration; Goldreich & Tremaine 1980; Ward 1997) – particularly in the inner
regions of disks, where magnetospheric cavities create bonafide traps for
planetary orbits (Masset et al., 2006).

In addition to the necessity of convergent migration, resonance capture
requires stability of the resonant equilibrium and adiabaticity. Crudely
speaking, this means that dissipative torques must not exceed gravitational
perturbations in magnitude, and that resonant dynamics must operate “faster”
than the timescale associated with extrinsic (that is, disk-driven) forcing of
the orbits. In this vein, Batygin (2015) proposed an analytic criterion for
adiabatic capture in the unrestricted 3-body problem, by equating the resonant
libration (bounded oscillation) period to the migratory resonance-crossing
time. While this criterion yields quantitatively adequate results in the
regime where orbital migration ensues in absence of other dissipative effects,
the recent simulation suite of Kajtazi et al. (2023) has shown that the
behavior of resonance capture is qualitatively different if convergent
migration is accompanied by strong eccentricity damping. Evidently, disk-
driven orbital circularization alters the efficiency of resonance capture in a
non-trivial manner.

The principle goal of this Letter is to understand the process of resonance
capture in presence of direct dissipation, from theoretical grounds. That is,
in this work, we employ perturbation theory to quantify conditions under which
stable resonant dynamics can be established, derive an analytic criterion for
such dissipative capture, and confirm our results with numerical experiments.
The remainder of the manuscript is organized as follows. In section 2, we
outline a simplified sketch of the resonance stability argument within the
context of the circular restricted 3-body problem. We generalize our
analytical framework to the unrestricted elliptic problem and compare our
results with numerical simulations in section 3. We summarize and discuss our
findings in section 4.

##  2 The Restricted Problem

As the simplest starting point for our analysis, we adopt the circular
restricted 3-body problem as a paradigm, wherein one of the secondary bodies
is taken to have negligible mass, while the other is assumed to reside on a
circular orbit. A similar approach has been undertaken in the recent study of
Huang & Ormel (2023). To be clear, we make the restricted approximation in
this section strictly for comprehensibility: the work of Sessin & Ferraz-Mello
(1984); Wisdom (1986), as well as a number of more recent studies (Batygin &
Morbidelli, 2013a; Petit et al., 2017; Hadden, 2019) have shown how the
perturbative treatment of first-order resonances within the restricted problem
can be generalized to the full 3-body problem111Within the context of the full
3-body problem, results depend predominantly on the sum of the planetary
masses, rather than their ratio (Deck et al., 2013; Deck & Batygin, 2015), and
we carry out this generalization in the next section.

###  2.1 Perturbation Theory

#### Model Hamiltonian.

Upon averaging over short-periodic terms and expanding the interaction
potential (i.e., the disturbing function) to leading order in eccentricity and
inclination, the governing Hamiltonian for a $k:k-1$ mean-motion resonance
takes the form:

$\displaystyle\mathcal{H}$
$\displaystyle=-\frac{\mathcal{G}\,M_{\star}}{2\,a^{\prime}}-\frac{\mathcal{G}\,m}{a^{\prime}}$
$\displaystyle\times\,f\,e^{\prime}\cos(k\,\lambda^{\prime}-(k-1)\,n\,t-\varpi^{\prime}).$
(1)

In the above expression, the Keplerian orbital elements have their usual
meanings, $f$ is constant of order unity222The coefficient $f$ depends only on
the semi-major axis ratio, and evaluates to $f\approx 1.2$ for $k=2$ and
$f\approx 0.8\,k$ for $k\geqslant 3$., and the primed variables refer to the
outer body, which we take to be massless. This choice circumvents any
consideration of over-stable librations, which can ensue if the massive
perturber resides on an exterior orbit (Deck & Batygin, 2015). We note
however, that the statistical analysis of Huang & Ormel (2023) indicates that
even under a reversed mass-ordering, over-stable librations are expected to be
rare in real protoplanetary disks. In addition, we have explicitly written the
mean longitude of the inner orbit explicitly as a product of its mean motion,
$n=\sqrt{\mathcal{G}\,M_{\star}/a^{3}}$, and time. In other words, the
physical setup of our problem has: $M\gg m\neq
0;m^{\prime}=0;e=0;\lambda=n\,t$.

To simplify the functional form of $\mathcal{H}$, we follow the well-
documented procedure of expanding the leading (Keplerian) term of equation
(2.1) in the vicinity of exact commensurability, to second order in
$\delta\mathcal{L}=\mathcal{L}-[\mathcal{L}]$, where
$\mathcal{L}=\sqrt{\mathcal{G}\,M_{\star}\,a’}$ and
$[\mathcal{L}]=\sqrt{\mathcal{G}\,M_{\star}\,[a’]}=\sqrt{\mathcal{G}\,M_{\star}\,(k/(k-1))^{2/3}\,a}$
represents the maximal specific angular momentum of the test-particle orbit,
evaluated at the nominal resonance semi-major axis, $[a^{\prime}]$. Switching
to the canonically conjugated action-angle variables (e.g., Peale 1986)

$\displaystyle\Phi=[\mathcal{L}]\,(e^{\prime})^{2}/2$
$\displaystyle\phi=k\,\lambda^{\prime}-(k-1)\,n\,t-\varpi^{\prime}$
$\displaystyle\Psi=\delta\mathcal{L}-k\,\Phi$
$\displaystyle\psi=\lambda^{\prime},$ (2)

Hamiltonian (2.1) takes the familiar form of the second fundamental model for
resonance (Henrard & Lemaitre, 1983):

$\displaystyle\mathcal{H}$
$\displaystyle=n^{\prime}(k\,\Phi+\Psi)-n\,(k-1)\,\Phi$
$\displaystyle-\frac{3}{2}\frac{n^{\prime}}{[\mathcal{L}]}(k\Phi+\Psi)^{2}-\mathcal{F}\,\sqrt{2\,\Phi}\cos(\phi),$
(3)

where $\mathcal{F}=(\mathcal{G}\,m/[a^{\prime}])\,(f/\sqrt{[\mathcal{L}]})$ is
a constant. Because the angle $\psi$ does not appear within $\mathcal{H}$, the
evolution of the conjugated action, $\Psi$, is dictated entirely by extrinsic
forces.

#### The case of pure migration.

Before proceeding to consider the full dissipative problem, let us pause and
recall some qualitative aspects of the well-studied instance of pure orbital
migration, where external (non-Hamiltonian) forces do not affect the
$(\Phi,\phi)$ degree of freedom directly. In this case, it is easy to see that
convergent migration will cause the action $\Psi$ to diminish without bound,
changing the topological structure of the phase-space portrait of
$\mathcal{H}$ in concert (Henrard & Lemaitre, 1983; Batygin, 2015). If the
evolution of $\Psi$ is slow compared to resonant dynamics, then an adiabatic
invariant – which corresponds to the phase-space area encircled by the orbit –
emerges as a quasi-integral of motion (Henrard, 1982; Neishtadt, 1984).

In the practically important case where orbits originate with zero
eccentricity far away from resonance, the initially occupied phase-space area
is null, meaning that as long as the adiabatic condition is satisfied333An
important additional caveat is that no encounters with the separatrix take
place., the trajectory must remain confined to the $\phi=\pi$ resonant
equilibrium point (since it encapsulates zero phase-space area).
Simultaneously, as $\Psi$ evolves to highly negative values, the eccentricity
grows perpetually as $e^{\prime}\sim\sqrt{-2\,\Psi/(k\,[\mathcal{L}])}$. In
this manner, forces that facilitate the convergent migration of orbits
translate to sustained eccentricity excitation, once the resonant lock is
established.

#### The case of concurrent migration and circularization.

The arguably more physically realistic scenario – wherein convergent migration
occurs together with efficient orbital damping – is different from the
aforementioned case of pure migration in a number of important ways. First and
foremost, the evolution of the action $\Psi$ is no longer unbounded, and
instead stabilizes at an equilibrium value, which in turn dictates the
equilibrium eccentricity. To quantify this, consider the following generic
parameterizations of semi-major axis decay and orbital circularization:

$\displaystyle\frac{1}{a^{\prime}}\frac{da^{\prime}}{dt}=-\frac{1}{\tau_{m}^{\prime}}-\frac{2\,e^{\prime
2}}{\tau_{m}^{\prime}/\mathcal{K}}$
$\displaystyle\frac{1}{e^{\prime}}\frac{de^{\prime}}{dt}=-\frac{\mathcal{K}}{\tau_{m}^{\prime}},$
(4)

where $\tau_{m}^{\prime}$ is the convergent migration timescale and
$\mathcal{K}=\tau_{m}^{\prime}/\tau_{e}^{\prime}$ is the ratio of semi-major
axis and eccentricity damping timescales.

It is worth noting that (Pichierri et al., 2022) have recently shown that the
eccentric contribution to semi-major axis damping in equations (4) arises
self-consistently within planet-disk interaction formulae that are routinely
implemented in $N-$body codes to mimic the effects of the gaseous nebula
(Papaloizou & Larwood, 2000). Moreover, for type-I migration, $\mathcal{K}$
has a well-defined dependence on the the geometric aspect ratio, $h/r$,
(Tanaka et al., 2002; Tanaka & Ward, 2004):

$\displaystyle\mathcal{K}=\bigg{(}\frac{1}{e^{\prime}}\frac{de^{\prime}}{dt}\bigg{)}\bigg{(}\frac{1}{a^{\prime}}\frac{da^{\prime}}{dt}\bigg{)}^{-1}\propto\bigg{(}\frac{h}{r}\bigg{)}^{-2},$
(5)

though the dimensionless pre-factor of this dependence is specified by the
disk’s particular structure. As an example, for a locally isothermal Mestel
(1963) type disk (where the surface density varies inversely with the semi-
major axis – i.e., $\Sigma\propto 1/a^{\prime}$), this pre-factor is
approximately $0.2$, such that $\mathcal{K}\sim\mathcal{O}(10^{2})$ for
$h/r\sim 0.05$ (see also Lee & Peale 2002). Nevertheless, values for
$\mathcal{K}$ that are substantially higher (and lower) are expected to arise
in realistic model nebulae; we will revisit the relevant scalings in the next
section.

With the relevant damping formulae defined, it is straightforward to derive
the equilibrium value for $e^{\prime}$. Recalling the definition of $\Psi$
from equations (2), let us set the time-derivative of $\Psi$ equal to zero:

$\displaystyle\frac{d\Psi}{dt}$
$\displaystyle=\frac{1}{2}\sqrt{\frac{\mathcal{G}\,M_{\star}}{a^{\prime}}}\frac{da^{\prime}}{dt}-2\,[\mathcal{L}]\,e^{\prime}\,\frac{de^{\prime}}{dt}$
$\displaystyle=-\frac{\mathcal{L}}{2}\frac{1+4\,\mathcal{K}\,\Phi/[\mathcal{L}]}{\tau_{m}^{\prime}}+2\,\mathcal{K}\,\frac{k\,\Phi}{\tau_{m}^{\prime}}=0.$
(6)

It is expected that equilibrium will be reached close to nominal
resonance444Notice that this is not the case for the case of dissipative
divergent migration, where the system can follow equilibrium loci far away
from nominal commensurability (Pichierri et al., 2019; Goldberg & Batygin,
2021). such that $\delta\mathcal{L}\approx 0$. Thus, replacing $\mathcal{L}$
by $[\mathcal{L}]$ in the above equation and recalling from equations (2) that
$\Phi=[\mathcal{L}]\,e^{\prime 2}/2$, we obtain:

$\displaystyle\big{(}e^{\prime}\big{)}_{\rm{eq}}\rightarrow\frac{1}{\sqrt{2\,(k-1)\,\mathcal{K}}}\sim\frac{h}{r}.$
(7)

Indeed, resonantly-excited planetary eccentricities within protoplanetary
nebulae are expected to be comparable to the disk aspect ratio (which is
itself equal to the inverse Mach number of the Keplerian flow), as the above
expression suggests (e.g., see also Pichierri et al. 2018).

The terminal step of the calculation is to evaluate the stability of the
resonant fixed point. Here, a second important distinction with the pure
migration case comes into view: in presence of explicit eccentricity damping,
the phase-space portrait rotates counter-clockwise, such that the equilibrium
value of the critical angle, $\phi$, shifts to $\phi_{\rm{eq}}=\pi+\epsilon$,
where $\epsilon$ is determined by the strength of dissipative effects (Batygin
& Morbidelli, 2013b). Importantly, criticality is achieved when
$\epsilon\rightarrow\pi/2$ and $\delta\mathcal{L}\rightarrow 0$ – a
configuration where the resonant torque is maximized. It thus follows that in
this state,

$\displaystyle\big{(}\Psi)_{\rm{eq}}\rightarrow-k\,\big{(}\Phi\big{)}_{\rm{eq}}=\frac{[\mathcal{L}]}{4\,(k-1)\,\mathcal{K}}.$
(8)

Figure 1: Numerical simulations of a $k=3$ resonant encounter within the
context of the circular restricted 3-body problem. The top and bottom panels
depict phase-space evolution and the time-series of the orbital frequency
ratio, respectively. Differently colored curves correspond to distinct
migration timescales, as labeled in the inset between the two panels. The
simulation setup is as follows: an exterior test particle, initialized on a
circular, planar $a^{\prime}=0.15\,$AU orbit around a $M_{\star}=1M_{\odot}$
star, migrates convergently towards a $m=1\times 10^{-5}\,M_{\odot}\approx
3\,M_{\oplus}$ planet residing at $a=0.1\,$AU, eventually encountering the 3:2
mean-motion commensurability. Eccentricity damping is applied with a
characteristic timescale that is a factor of $\mathcal{K}=1000$ shorter than
the migration time, $\tau_{m}^{\prime}$. Simultaneous migration and
circularization of the outer orbit is indicated with shaded arrows on the
diagram. As $\tau_{m}^{\prime}$ is reduced from $80\,$kyr toward the
analytically-predicted critical value of $5\,$kyr (equation 10), the
equilibrium value of the resonant angle, $\phi$, tends from $\pi$ to $3\pi/2$.
The equilibrium value of the eccentricity, on the other hand, stabilizes at
$e^{\prime}=1/\sqrt{4\,\mathcal{K}}\approx 0.016$ in all cases. The behavior
of this heavily damped system is fully consistent with expectations provided
by analytic theory.

Accounting for dissipative effects, the equilibrium equations take the form:

$\displaystyle\frac{d\phi}{dt}$
$\displaystyle=\frac{\partial\mathcal{H}}{\partial\Phi}=k\,n^{\prime}-(k-1)\,n-\frac{\mathcal{F}\,\cos((\phi)_{\rm{eq}})}{\sqrt{2\,(\Phi)_{\rm{eq}}}}$
$\displaystyle-\frac{3\,k\,n^{\prime}\,((\Psi)_{\rm{eq}}+k\,(\Phi)_{\rm{eq}})}{[\mathcal{L}]}=0$
$\displaystyle\frac{d\Phi}{dt}$
$\displaystyle=-\frac{\partial\mathcal{H}}{\partial\phi}-2\,\mathcal{K}\,\frac{(\Phi)_{\rm{eq}}}{\tau_{m}^{\prime}}=-\mathcal{F}\,\sqrt{2(\Phi)_{\rm{eq}}}\,\sin((\phi)_{\rm{eq}})$
$\displaystyle-2\,\mathcal{K}\,\frac{(\Phi)_{\rm{eq}}}{\tau_{m}^{\prime}}=0.$
(9)

The first $(\dot{\phi}=0)$ equation is trivially satisfied. The second
$(\dot{\Phi}=0)$ equation, on the other hand, yields the criterion for the
shortest migration timescale that allows for resonant capture:

$\displaystyle\tau_{m}^{\prime}$
$\displaystyle=\frac{M_{\star}}{f\,m}\sqrt{\frac{\mathcal{K}}{2\,(k-1)}\frac{a^{\prime
3}}{\mathcal{G}\,M_{\star}}}$
$\displaystyle\approx\frac{5}{4}\frac{M_{\star}}{m}\sqrt{\frac{\mathcal{K}}{2\,(k-1)^{3}}}\frac{1}{n},$
(10)

where we have used the compact $k\approx k-1$ approximation to evaluate
$f\approx 4\,k/5$. The expression agrees with the one recently obtained by
Huang & Ormel (2023), who arrived at it through a somewhat distinct – but
ultimately equivalent – approach.

###  2.2 Numerical Experiments

Whether the process of resonant capture is controlled by stability – and thus
follows the criterion given by equation (10) – or adiabaticity (as discussed
in e.g., Batygin 2015), depends on the efficiency of orbital circularization.
In the limit of weak damping (“small $\mathcal{K}$”), we may reasonably expect
that adiabaticity will serve as the more stringent constraint, while stability
will dominate in the regime of rapid (“large $\mathcal{K}$”) damping.
Moreover, the transitionary value of $\mathcal{K}$ is likely to significantly
exceed unity, since the eccentricity-damping timescale should be contrasted
against the resonance-crossing time – a quantity that is proportional to
$\tau_{m}^{\prime}$, but is much smaller in magnitude. In the remainder of
this section, we will use numerical experiments to verify the validity of
equation (10), as well as to explore the capture process in the heavily
($\mathcal{K}=1000$) and moderately ($\mathcal{K}=100$) damped regimes.

Our numerical experiments follow the conventional scheme of implementing
dissipative effects into an $N-$body framework using the formulae of
Papaloizou & Larwood (2000). For definitiveness, we adopted a physical setup
that is reminiscent of short-period extrasolar Super-Earth systems such as
Kepler-59 and Kepler-128: a $1\,M_{\odot}$ star encircled by a
$m=10^{-5}\,M_{\odot}\approx 3\,M_{\oplus}$ planet on a circular orbit at
$a=0.1\,$AU (Hadden & Lithwick, 2016; Saad-Olivera et al., 2020). In addition,
we initialized an exterior test-particle into the simulation on a
$e^{\prime}=0$, $a=0.15\,$AU orbit, such that the first commensurability
encountered by the system is the 3:2 period-ratio. To drive convergent
migration, non-gravitational forces were only applied to the test-particle.
The system of ODEs was integrated using the Bulirsch-Stoer algorithm with an
accuracy parameter of $\hat{\epsilon}=10^{-10}$.

In both the $\mathcal{K}=1000$ and $\mathcal{K}=100$ simulation suites, we
carried out five numerical experiments, first setting $\tau^{\prime}_{\rm{m}}$
equal to the value given by the stability criterion (10) and doubling it in
every run. Figures (1) and (2) depict the results of these simulations: the
bottom panels show the time-series of the orbital frequency ratio (equal to
the period ratio), and the top panels show the phase-space evolution of the
system.

Overall, the results of the numerical experiments with $\mathcal{K}=1000$
conform to the theoretical expectations outlined above. That is, as the test-
particle approaches the 3:2 resonance, $e^{\prime}$ stabilizes at the
equilibrium value given by equation (7) – shown on the phase-space plot with a
gray circle – independent of the adopted $\tau_{m}^{\prime}$. Conversely, the
stationary value of the critical angle ratchets up towards $\phi=3\,\pi/2$ as
$\tau_{m}^{\prime}$ approaches the critical value of $\sim 5000$ years, as
dictated by equation (10). Though all of the shown runs result in stable
capture within a 3:2 resonance, we have also confirmed that resonant locking
fails for shorter (e.g., $4500$ year) migration timescales, in agreement with
the predictions of the analytical theory.

Unlike their more heavily-damped counterparts, the $\mathcal{K}=100$
simulations deviate notably from our analytical stability arguments. While the
system equilibrates at the predicted state as long as the convergence time is
long, the process of resonance capture is accompanied by a growing libration
amplitude, as $\tau_{m}^{\prime}$ approaches the critical value of $\sim 1500$
years (recall that it scales as $\propto\sqrt{\mathcal{K}}$). For this reason,
the resonant equilibrium becomes compromised at a migration timescale that is
a factor of $\sim 2$ longer than that predicted by equation (10). Indeed, the
run with $\tau_{m}^{\prime}=1500$ years results in passage through the 3:2
commensurability and subsequent capture into the 5:4 resonance (a
configuration reminiscent of the Kepler-307 system; Jontof-Hutter et al.
2016). Evidently, for the given mass-ratio $m/M_{\star}$, a value of
$\mathcal{K}$ significantly in excess of $100$ is required to fully suppress
the growth of the phase-space area during orbital convergence, giving way to
adiabaticity as the process that largely determines the outcome of resonant
encounters.

Figure 2: Same as Figure 1 but with a reduced eccentricity damping factor of
$\mathcal{K}=100$. While resonant capture ensues for long migration
timescales, the analytically-predicted critical value of
$\tau_{m}^{\prime}=1.5\,$kyr leads to passage through the 3:2 commensurability
and capture into the 5:4 resonance instead. Failure of the stability criterion
(10) in this example indicates that the orbital circularization is
sufficiently slow that a different mechanism – namely, adiabaticity –
regulates resonance capture.

##  3 The Unrestricted Problem

###  3.1 Analytical Theory

With the qualitative picture outlined within the simplified framework of the
restricted problem above, generalization of our results to the full resonant
three-body problem is relatively undemanding. Retaining the nearly coplanar
and low-eccentricity approximations but putting no limits on the planetary
mass-ratio $m/m^{\prime}$, the governing Hamiltonian takes the form:

$\displaystyle\mathcal{H}$
$\displaystyle=-\frac{m^{3}}{2}\bigg{(}\frac{\mathcal{G}\,M_{\star}}{\Lambda}\bigg{)}^{2}-\frac{m^{\prime
3}}{2}\bigg{(}\frac{\mathcal{G}\,M_{\star}}{\Lambda^{\prime}}\bigg{)}^{2}$
$\displaystyle-\mathcal{A}\sqrt{2\,\Gamma}\,\cos(k\,\lambda^{\prime}-(k-1)\,\lambda+\gamma)$
$\displaystyle-\mathcal{B}\sqrt{2\,\Gamma^{\prime}}\,\cos(k\,\lambda^{\prime}-(k-1)\,\lambda+\gamma^{\prime}),$
(11)

where $\mathcal{A}=(\mathcal{G}^{2}\,M_{\star}\,m\,m^{\prime
3}/\Lambda^{\prime 2})\,(g/\sqrt{\Lambda})$ and
$\mathcal{B}=(\mathcal{G}^{2}\,M_{\star}\,m\,m^{\prime 3}/\Lambda^{\prime
2})\,(f/\sqrt{\Lambda^{\prime}})$ are pre-factors similar to $\mathcal{F}$
(Batygin & Morbidelli, 2013a; Hadden, 2019). Furthermore, in the above
expression, we have used the conventional system of Poincaré action-angle
variables, defined in the $m\ll M_{\star}$ and $e\ll 1$ limit as:

$\displaystyle\Lambda=m\,\sqrt{\mathcal{G}\,M_{\star}\,a}$
$\displaystyle\lambda=\mathcal{M}+\varpi$
$\displaystyle\Gamma=\Lambda\,e^{2}/2$ $\displaystyle\gamma=-\varpi,$ (12)

with equivalent definitions for the primed quantities.

Figure 3: A dimensionless mass vs. migration time map of $k=3$ resonance
capture for an equal-mass ($m=m^{\prime}$) planetary system. The filled
circles indicate parameter combinations where numerical simulations yield
successful resonant capture, whereas empty circles indicate passage through
the 3:2 resonance. The left and right panels correspond to the heavily
($\mathcal{K}=1000$) and moderately ($\mathcal{K}=100$) damped regimes,
respectively. Both panels additionally show analytic stability and
adiabaticity capture criteria as dashed and solid lines. While the stability
criterion regulates capture in the heavily-dissipated case, adiabaticity
controls the outcome of resonant encounters in the moderately damped regime.

Notice that unlike equation (2.1), the Hamiltonian (11) now contains two
harmonic terms, and we have maintained the Keplerian contributions in their
unexpanded form (this will not affect our analysis). As in the preceding
section, we must now consider the stability of the global fixed point of the
resonant Hamiltonian, in presence of migration and dissipation. In fact,
similar analyses have previously been carried out in the literature (Batygin &
Morbidelli, 2013b; Terquem & Papaloizou, 2019), but typically in the limit of
weak friction.

Central to our calculation is the long-term evolution of the semi-major axes,
which – unlike the case of the restricted problem – can be directed inward,
outward, or can be stationary, depending on the signs and magnitudes of the
individual planetary migration timescales, $\tau_{m}$ and $\tau_{m}^{\prime}$.
Nevertheless, the maintenance of the resonant relationship between the orbits
necessitates that
$\dot{\Lambda}/\Lambda=\dot{\Lambda^{\prime}}/\Lambda^{\prime}$. This equality
can be computed directly from Hamilton’s equation:

$\displaystyle\frac{1}{\Lambda}\frac{d\Lambda}{dt}=-\frac{\partial\,\mathcal{H}}{\partial\,\lambda}-\frac{\Lambda\,\big{(}\tau_{e}+2\,\tau_{m}(2\,\Gamma/\Lambda)\big{)}}{2\,\tau_{m}\,\tau_{e}},$
(13)

where we have expressed the dissipative contribution (given by equation 4) in
canonical coordinates. Collecting the Hamiltonian terms on the LHS, we have:

$\displaystyle\frac{2\,\big{(}k\,(\Lambda+\Lambda^{\prime})-\Lambda^{\prime}\big{)}}{\Lambda\,\Lambda^{\prime}}\bigg{(}\mathcal{A}\,\sqrt{2\,\Gamma}\,\sin(\varphi)+\mathcal{B}\,\sqrt{2\,\Gamma^{\prime}}\,\sin(\phi)\bigg{)}$
$\displaystyle=\frac{1}{\tau_{m}}+\frac{4\,\Gamma}{\Lambda\,\tau_{e}}-\frac{1}{\tau_{m}^{\prime}}-\frac{4\,\Gamma^{\prime}}{\Lambda^{\prime}\,\tau_{e}^{\prime}}.$
(14)

Similarly to equation (2), in the above expression, we have used $\varphi$ and
$\phi$ to denote the resonant harmonics containing $\varpi$ and
$\varpi^{\prime}$, respectively.

The dependence of equation (14) on the resonant angles as well as the actions
$\Gamma$ and $\Gamma^{\prime}$ can be eliminated by considering the
equilibrium of the eccentricities. This equilibrium is given by
$\dot{\Gamma}=-\partial\mathcal{H}/\partial\gamma-2\,\Gamma/\tau_{e}=0$, with
an identical expression for $\Gamma^{\prime}$. In particular, we obtain:

$\displaystyle\mathcal{A}\,\sqrt{2\,\Gamma}\,\sin(\varphi)=-2\,\Gamma/\tau_{e}$
$\displaystyle\mathcal{B}\,\sqrt{2\,\Gamma^{\prime}}\,\sin(\phi)=-2\,\Gamma^{\prime}/\tau_{e}^{\prime}.$
(15)

From this expression, it is easy to see how the sinusoidal terms in equation
(14) can be eliminated. Moreover, in analogy with the preceding section,
criticality is attained as $\varphi\rightarrow\pi/2$, $\phi\rightarrow
3\pi/2$, which yields
$\Gamma=\mathcal{A}^{2}\,\tau_{e}^{2}/2;\Gamma^{\prime}=\mathcal{B}^{2}\,\tau_{e}^{\prime
2}/2$. Upon direct substitution into equation (14) and setting
$a=((k-1)/k)^{2/3}a^{\prime}=\alpha\,a^{\prime}$, we obtain the criterion for
resonant capture:

$\displaystyle\frac{1}{\tau_{m}^{\prime}}-\frac{1}{\tau_{m}}=\frac{2\,k\,\mathcal{G}\,M_{\star}}{a^{\prime
3}}\bigg{(}\frac{m}{M_{\star}}+\frac{m^{\prime}}{\sqrt{\alpha}\,M_{\star}}\bigg{)}$
$\displaystyle\times\bigg{(}\frac{g^{2}\,m^{\prime}}{\sqrt{\alpha}\,M_{\star}}\tau_{e}+\frac{k-1}{k}\frac{f^{2}\,m}{M_{\star}}\tau_{e}^{\prime}\bigg{)}$
$\displaystyle\approx\frac{32\,\mathcal{G}\,k^{3}\,\big{(}m+m^{\prime}\big{)}\,\big{(}m^{\prime}\,\tau_{e}+m\,\tau_{e}^{\prime}\big{)}}{25\,M_{\star}\,a^{\prime
3}}$ (16)

It is trivial to check that this criterion reproduces equation (10) in the
limit where $m^{\prime}\rightarrow 0$.

Although here we have derived criterion (16) directly from Hamiltonian (11),
an equivalent – albeit somewhat more mathematically involved – approach would
have been to first reduce $\mathcal{H}$ to an integrable form that only
contains a single resonant harmonic (see e.g., Batygin & Morbidelli 2013a),
and then analyze the stability of its equilibrium under dissipation. This
approach yields identical results to those delineated above and is reproduced
in the Appendix.

###  3.2 Numerical Experiments

As a quantitative test of the generalized criterion derived above, we have
repeated the numerical simulations described in the previous section, this
time setting $m^{\prime}=m$, and sampling a broad range of orbital convergence
times and planetary masses, in order to map the capture criterion for the 3:2
resonance on the ($m/M_{\star},\bar{\tau}_{m}/P$) plane, where $P$ is the
orbital period of the inner object. For definitiveness, in these simulations,
we assigned a common value of $\mathcal{K}$ to both planets (such that
$\tau_{e}=\tau_{e}^{\prime}=\bar{\tau}_{m}/\mathcal{K}$), but applied orbital
decay only to the outer body such that
$\bar{\tau}_{m}=\tau_{m}\,\tau_{m}^{\prime}/(\tau_{m}-\tau_{m}^{\prime})=\tau_{m}^{\prime}$.
Finally, to prevent the system from spiraling onto the central star without
careful modeling of the disk’s inner edge (e.g., Izidoro et al. 2017, 2021),
we simply rescaled both of the semi-major axes at every time-step555We have
simulated the assembly of the Galilean moons into the Laplace resonance using
an identical method in a previous study (Batygin & Morbidelli, 2020).,
maintaining the inner planet at $a=0.1\,$AU.

The results of these numerical experiments are shown in Figure (3). Instances
where capture into the 3:2 resonance was successful are shown with filled gray
points whereas runs that resulted in passage through the commensurability are
shown with empty circles. In both the $\mathcal{K}=1000$ (left panel) and
$\mathcal{K}=100$ (right panel) simulation suites, a clear power-law threshold
emerges on the diagrams, though the slope of this boundary is subtly distinct.
As already shown in the proceeding section, for the adopted mass-ratios, the
$\mathcal{K}=1000$ case is regulated by stability, whereas the
$\mathcal{K}=100$ case is controlled by adiabaticity. To confirm this
expectation, we have over-plotted the $k=3$, $m=m^{\prime}$ stability and
adiabaticity criteria (see Batygin 2015):

$\displaystyle\bigg{(}\frac{\bar{\tau}_{m}}{P}\bigg{)}_{\rm{stab}}=\frac{5}{32\,\pi}\frac{M_{\star}}{m}\sqrt{\frac{\mathcal{K}}{6}}$
$\displaystyle\bigg{(}\frac{\bar{\tau}_{m}}{P}\bigg{)}_{\rm{ad}}=\frac{5}{384}\bigg{(}\frac{125}{4}\bigg{)}^{1/9}\bigg{(}\frac{M_{\star}}{m}\bigg{)}^{4/3},$
(17)

with dashed and solid lines, respectively. These two criteria adequately
explain the numerical results and highlight the distinct regimes of resonance
capture that can ensue within moderately and heavily dissipative planet-
formation environments.

###  3.3 Scaling for Type-I Migration

As a final theme of this section, let us move away from parameterized
simulations considered above, and compare our results with more realistic
simulations of resonance capture driven by type-I disk migration. For a planet
of mass $m$ residing on an orbit with semi-major axis $a$, the characteristic
timescale associated with type-I migration is the spiral density wave
propagation time (Tanaka & Ward, 2004):

$\displaystyle\tau_{\rm{wave}}=\frac{1}{n}\frac{M_{\star}}{m}\frac{M_{\star}}{\Sigma\,a^{2}}\bigg{(}\frac{h}{r}\bigg{)}^{4},$
(18)

where $\Sigma$ is the surface density of the nebula, evaluated at the
planetary semi-major axis.

For nearly circular and planar orbits, the eccentricity damping timescale
differs from $\tau_{\rm{wave}}$ only by a numerical factor of order unity,
$f_{e}=0.78$, such that $\tau_{e}=\tau_{\rm{wave}}/f_{e}$. The semi-major axis
damping time, on the other hand, exceeds the wave propagation time by a large
margin: $\tau_{m}=\tau_{\rm{wave}}\,(r/h)^{2}/f_{a}$. Notably, the
proportionality constant depends on the index of the disk’s surface density
profile, $s$, in a linear manner: $f_{a}=2.7+1.1\,s$ (Tanaka et al., 2002).
Employing these dependencies, Kajtazi et al. (2023) used the surface density
of the disk as a physical proxy for migration speed, and formulated the
results of their numerical simulation suite in terms of a critical value of
$\Sigma$ required for resonance capture.

Carrying out their numerical experiments in a non-evolving, flared, $s=1$
nebula, Kajtazi et al. (2023) simulated the orbital convergence of a pair of
equal-mass planets, suppressing disk-driven migration of the inner object,
which was initialized at $0.1\,$AU. Intriguingly, they found that the critical
value of $\Sigma$ exhibits a clear dependence on the cube of the disk’s aspect
ratio but is independent of the planetary mass (i.e., critical
$\Sigma\propto\,m^{0}\,h^{3}$). Let us examine if these scalings can be
derived from our analytical criterion.

Plugging in the aforementioned type-I expressions for $\tau_{m}$ and
$\tau_{e}$ into equation (16), we obtain:

$\displaystyle\frac{f_{a}\,m^{\prime}\,n^{\prime}\,\Sigma\,a^{\prime
2}}{M_{\star}}\bigg{(}\frac{h}{r}\bigg{)}^{-2}=\frac{128\,\mathcal{G}\,k^{3}\,m^{\prime}\,M_{\star}}{25\,f_{e}\,n^{\prime}\,\Sigma\,a^{\prime
5}}\bigg{(}\frac{h}{r}\bigg{)}^{4},$ (19)

where we have set $\bar{\tau}_{m}=\tau_{m}^{\prime};\,m=m^{\prime}$ following
Kajtazi et al. (2023), and have assumed the compact approximation valid for
$k\gtrsim 3$ for simplicity. Solving for $\Sigma$ and relating the value to
the reference surface density, $\Sigma_{0}$, at $r_{0}=1\,$AU, we recover the
scalings found in simulations:

$\displaystyle\frac{\Sigma\,a^{\prime
2}}{M_{\star}}=\frac{\Sigma_{0}\,r_{0}\,a^{\prime}}{M_{\star}}=\sqrt{\frac{2\,k^{3}}{f_{a}\,f_{e}}}\bigg{(}\frac{h}{r}\bigg{)}^{3}.$
(20)

As a concluding step, to assess the degree of quantitative agreement between
theory and numerical experiments, in Figure (4), we plot our analytical
criterion for a range of resonance indexes together with the numerical results
of (Kajtazi et al., 2023, their Fig. 6). Though the agreement is not exact, we
find that the analytical result does not deviate from the numerical findings
by more than 20%.

Figure 4: Comparison between analytic and numeric determination of resonance
capture. Using the nebular surface density as a proxy for the critical
migration rate, the results of the Kajtazi et al. (2023) simulation suite
($\mathcal{K}\approx 570$) are shown with blue points for a range of resonant
indexes, $k$. Analytic stability and adiabaticity criteria are over-plotted
with gray and black points, respectively. Although not exact, the critical
value of the surface density given by the stability criterion (equation 16)
matches the numerical data to within $20\%$.

##  4 Discussion

Capture of planets into mean-motion resonances is an expected outcome of
evolution within protoplanetary nebulae, and in this work, we have derived an
analytic criterion for resonance locking in presence of migration and orbital
circularization. Our theoretical framework is based upon a perturbative
treatment of the unrestricted gravitational three-body problem (Peale 1986;
Murray & Dermott 1999 and the refs therein), and thus assumes small
eccentricities and inclinations, while placing no restrictions on the planet-
planet mass ratio. Fundamentally, our derivation rests upon the stability
analysis of the resonant equilibrium in presence of dissipation, and
complements the previously obtained criterion that stems from a consideration
of adiabaticity (Batygin, 2015).

We note that the stability argument considered here concerns the elementary
question of the existence of an equilibrium point for the resonant variables
and not the long-term post-capture evolution. To this end, it has been shown
by Goldreich & Schlichting (2014); Deck & Batygin (2015); Xu et al. (2018)
that depending on the ratio of the planetary masses and eccentricity damping
timescales, systems where resonant capture is initially successful can
eventually escape from resonance by way of over-stability. Although beyond the
scope of our study, these dynamics, along with stochastic growth of the
libration amplitude by nebular turbulence (Adams et al., 2008; Batygin &
Adams, 2017), synodic modulation of resonant angles (Pichierri & Morbidelli,
2020; Goldberg et al., 2022b), etc., give rise to additional constraints on
the kind of architectures that can be established within protoplanetary disks.

It is interesting to note that although the criteria for resonant capture
based upon stability and adiabaticity yield distinct scalings with the
planetary mass, for typical parameters pertinent to planetary migration within
circumstellar disks, they can give results that are quantitatively similar.
Ultimately, for resonant locking to ensue, both criteria – stability and
adiabaticity – must be satisfied, such that whichever one yields the longer
critical timescale for orbital convergence plays the controlling role.
Accordingly, the specific capture regime is dictated by the combination of the
ratio of the cumulative planetary mass to the mass of the star as well as the
ratio of the migration and circularization timescales.

For sub-Jovian planet pairs, the analysis carried out in section (3) indicates
that the outcomes of resonant encounters in systems with $\mathcal{K}=10^{2}$
and $\mathcal{K}=10^{3}$ are determined by adiabaticity and stability
respectively, such that the transitionary value of the timescale ratio lies in
between these two regimes. While the appropriate value of $\mathcal{K}$ is
bound to be case-specific, comparison of our results with the published
simulation suite of Kajtazi et al. (2023), indicates the that stability
criterion (16) reproduces the scalings seen in the numerical experiments.
Still, it remains important to keep in mind that within the context of the
flared disk model of Kajtazi et al. (2023), the disk aspect ratio at a
stellocentric distance of $0.1\,$AU (where the resonant encounter is set up to
take place) is $(h/r)=0.033\times(0.1)^{1/4}\approx 0.019$, corresponding to
$\mathcal{K}\approx 570$. Thus, it is reasonable to expect that if the
resonant encounters were to instead take place at $10\,$AU, where
$(h/r)\gtrsim 0.05$ and $\mathcal{K}\lesssim 100$, adiabaticity would have
played the deterministic role. In summary, the analytical framework developed
here adds to our understanding of the dynamics of resonant encounters, and
serves as a basis for further interpretation of the early evolution of
planetary systems.

#### Acknowledgments.

K. B. is grateful to Caltech, the David and Lucile Packard Foundation, and the
National Science Foundation (grant number: AST 2109276) for their generous
support. During the preparation of this paper, we have become aware that Huang
& Ormel (2023, submitted) arrived at similar arguments simultaneously and
independently.

## Derivation of the Stability Criterion From an Integrable Hamiltonian

In this appendix, we present an alternative derivation of the stability
criterion given by equation (16), based upon the reduction of the resonant
three-body problem to an integrable Hamiltonian. In particular, we follow the
formalism outlined in Deck & Batygin (2015), switching to their notation for
consistency. Because this reduction is well documented (_e.g._ Batygin &
Morbidelli, 2013a; Deck et al., 2013; Petit et al., 2017; Hadden, 2019), we
limit ourselves to recalling the relevant variables and refer the interested
reader to the cited works.

As in the main text, we consider the general problem of two planets of mass
$m_{1}$ and $m_{2}$ orbiting a star of mass $M_{\star}$ in the plane. We
assume that the planets experience a convergent migration and are close to the
crossing of the $k$:$k-1$ resonance. The integrable Hamiltonian has the form:

$\displaystyle\mathcal{H}=-\frac{1}{2}(\Phi^{\prime}-\Gamma^{\prime})^{2}-\sqrt{2\Phi^{\prime}}\cos(\phi),$
(1)

where $\Phi^{\prime}$ is the renormalized Sessin variable (Sessin & Ferraz-
Mello, 1984)

$\displaystyle\Phi^{\prime}$
$\displaystyle=\frac{1}{Q}\frac{k}{k-1}\frac{\alpha_{\rm{res}}}{\zeta+\alpha_{\rm{res}}}\frac{\zeta\sqrt{\alpha_{\rm{res}}}}{2(R^{2}+\zeta\sqrt{\alpha_{\rm{res}}})}\sigma^{2},$
(2) $\displaystyle\sigma^{2}$
$\displaystyle=R^{2}e_{1}^{2}+e_{2}^{2}-2Re_{1}e_{2}\cos(\Delta\varpi),$

where $R=|f_{k,27}(\alpha_{\rm{res}})|/f^{\prime}_{k,31}(\alpha_{\rm{res}})$,
$\zeta=m_{1}/m_{2}$, $\alpha_{\rm{res}}=(1-1/k)^{2/3}$ and the functions
$f_{27/31}$ are the resonant coefficients (Murray & Dermott, 1999). $Q$ is a
renormalization factor for the actions

$Q=\varepsilon_{\rm
p}^{2/3}\zeta\frac{\alpha_{\rm{res}}^{5/6}}{k-1}\left(\frac{f^{\prime
2}_{k,31}}{9k(1+\zeta)^{2}}\frac{R^{2}+\zeta\sqrt{\alpha_{\rm{res}}}}{(\alpha_{\rm{res}}+\zeta)^{5}}\right)^{1/3}$
(3)

defined by (Deck & Batygin, 2015, their Appendix A) and $\varepsilon_{\rm
p}=(m_{1}+m_{2})/M_{\star}$. The angle $\phi$ is conjugated with
$\Phi^{\prime}$ and is the generalized resonant angle. $\Gamma^{\prime}$, on
the other hand, is a parameter that is related to the system’s angular
momentum.

Maintaining the notation of Deck & Batygin (2015), we have the relative
migration timescale

$\frac{1}{\tau_{a}}=\frac{1}{\tau_{a,2}}-\frac{1}{\tau_{a,1}},$ (4)

and the eccentricity damping timescale

$\frac{1}{\tau_{e}}=\frac{1}{\tau_{e,1}}+\frac{\zeta}{\tau_{e,2}}.$ (5)

Note that without making any approximations in the computation of the
dissipation onto the action $\Phi^{\prime}$, one obtains the following
eccentricity damping timescale

$\frac{1}{\tau_{e}}=\frac{1}{\tau_{e,1}}+\frac{\sqrt{\alpha}}{R^{2}}\frac{\zeta}{\tau_{e,2}}$
(6)

This expression slightly differs from that given in Deck & Batygin (2015) by a
few percent for $k\geq 2$ but is more accurate for the 2:1 MMR. Deck & Batygin
(2015) also define the timescale

$\frac{1}{\tau_{a,e}}=\frac{1}{\tau_{e,1}}-\frac{{\alpha}}{R^{2}}\frac{\zeta^{2}}{\tau_{e,2}}.$
(7)

The damping on the action $\Phi^{\prime}$ is expressed as

$\left.\frac{\mathrm{d}\Phi^{\prime}}{\mathrm{d}t^{\prime}}\right|_{\rm
dis}=C_{0}\Phi^{\prime}$ (8)

where we assume $\tau_{a}\gg\tau_{e}$; neglecting higher order terms in
eccentricities, we have:

$C_{0}=-\frac{2\gamma_{1}}{\tau_{e}},$ (9)

where

$\gamma_{1}=\frac{R^{2}}{R^{2}+\zeta\sqrt{\alpha}},$ (10)

and for completness we define $\gamma_{2}=1-\gamma_{1}$. Damping on the
Hamiltonian parameter $\Gamma^{\prime}$ is expressed as

$\left.\frac{\mathrm{d}\Gamma^{\prime}}{\mathrm{d}t^{\prime}}\right|_{\rm
dis}=A_{0}+A_{1}\Phi^{\prime}$ (11)

with

$A_{0}=\frac{\zeta\sqrt{\alpha}}{2Qk\eta_{2}^{2}}\frac{1}{\tau_{a}}\quad\mathrm{and}\quad
A_{1}=\frac{-2p\gamma_{1}}{k\eta_{2}}\frac{1}{\tau_{a,e}},$ (12)

where

$\eta_{2}=\frac{k-1}{k}+\zeta\sqrt{\alpha}$ (13)

with $\eta_{1}=\eta_{2}/(\zeta\sqrt{\alpha})$.

The equations of motion666Expressed in terms of the renormalized time
$t^{\prime}$. of the dissipative problems are

$\displaystyle\frac{\mathrm{d}\Phi^{\prime}}{\mathrm{d}t^{\prime}}$
$\displaystyle=-\sqrt{2\Phi^{\prime}}\sin\phi+\frac{\eta_{2}^{3}C_{0}\Phi^{\prime}}{n_{2}Q|\mathcal{K}_{2}|}$
(14) $\displaystyle\frac{\mathrm{d}\phi}{\mathrm{d}t^{\prime}}$
$\displaystyle=-(\Phi^{\prime}-\Gamma^{\prime})-\frac{1}{\sqrt{2\Phi^{\prime}}}\cos\phi$
(15) $\displaystyle\frac{\mathrm{d}\Gamma^{\prime}}{\mathrm{d}t^{\prime}}$
$\displaystyle=\frac{\eta_{2}^{3}}{n_{2}Q|\mathcal{K}_{2}|}(A_{0}+(A_{1}+C_{0})\Phi^{\prime})$
(16)

where $n_{j}$ is the mean motion and
$\mathcal{K}_{2}=-3(k-1)^{2}(\alpha_{\rm{res}}+\zeta)^{5}/(\zeta\alpha_{\rm{res}})$
is the second order derivative of the Keplerian Hamiltonian at the exact
resonance. Looking for an equilibrium, we have

$\Phi_{\rm eq}^{\prime}=\frac{-A_{0}}{C_{0}+A_{1}}$ (17)

and,

$\sin\phi_{\rm eq}=\frac{\eta_{2}^{3}\sqrt{\Phi_{\rm
eq}^{\prime}}C_{0}}{\sqrt{2}Q|\mathcal{K}_{2}|n_{2}}.$ (18)

The equilibrium exists if the absolute value of the right hand side of the
above expression is smaller than 1. Accordingly, one can derive a condition on
the migration speed $\tau_{a}$ that is, without further approximation:

$\frac{1}{\tau_{a}n_{2}}<\sqrt{2}\varepsilon_{\rm
p}\sqrt{\frac{\tau_{e}}{\tau_{a}}}G(k,\zeta,p_{a,e})$ (19)

where $G$ is a function solely of the resonant index $k$, the planet mass
ratio $\zeta$ and the ratio of the eccentricity damping timescales through the
parameter $p_{a,e}$

$G(k,\zeta,p_{a,e})=\frac{\zeta\sqrt{\alpha_{\rm{res}}}+R^{2}}{1+\zeta}\sqrt{\frac{\frac{k-1}{k}+\zeta\sqrt{\alpha_{\rm{res}}}}{R^{2}\alpha_{\rm{res}}}}f^{\prime}_{k,31}\sqrt{k}\sqrt{1+p_{a,e}},$
(20)

where

$p_{a,e}=\frac{p}{k-1+k\zeta\sqrt{\alpha_{\rm{res}}}}\frac{\tau_{e}}{\tau_{a,e}}.$
(21)

The expression (19) is equivalent to Eq. (16) for $p=0$.

In the case of an outer test particle, $\zeta\to+\infty$, using expression Eq.
(6) for $\tau_{e}$, the criterion becomes

$\frac{1}{\tau_{a}n_{2}}<\sqrt{2}\varepsilon_{\rm
p}\sqrt{\frac{\tau_{e,2}}{\tau_{a}}}f^{\prime}_{k,31}\sqrt{k}\sqrt{1-\frac{p}{k}}.$
(22)

Taking $p=0$ yields the expression given by equation (10) derived from the
restricted problem. Similarly $p=1$ corresponds to the result obtained by
Huang & Ormel (2023).

## References

  * Adams et al. (2008) Adams, F. C., Laughlin, G., & Bloch, A. M. 2008, ApJ, 683, 1117. doi:10.1086/589986
  * Batygin & Brown (2010) Batygin, K. & Brown, M. E. 2010, ApJ, 716, 1323. doi:10.1088/0004-637X/716/2/1323
  * Batygin & Morbidelli (2013a) Batygin, K. & Morbidelli, A. 2013a, A&A, 556, A28. doi:10.1051/0004-6361/201220907
  * Batygin & Morbidelli (2013b) Batygin, K. & Morbidelli, A. 2013b, AJ, 145, 1. doi:10.1088/0004-6256/145/1/1
  * Batygin (2015) Batygin, K. 2015, MNRAS, 451, 2589. doi:10.1093/mnras/stv1063
  * Batygin & Adams (2017) Batygin, K. & Adams, F. C. 2017, AJ, 153, 120. doi:10.3847/1538-3881/153/3/120
  * Batygin & Morbidelli (2020) Batygin, K. & Morbidelli, A. 2020, ApJ, 894, 143. doi:10.3847/1538-4357/ab8937
  * Batygin & Morbidelli (2023) Batygin, K. & Morbidelli, A. 2023, Nature Astronomy. doi:10.1038/s41550-022-01850-5
  * Cresswell & Nelson (2006) Cresswell, P. & Nelson, R. P. 2006, A&A, 450, A833. doi:10.1051/0004-6361:20054551
  * Deck et al. (2013) Deck, K. M., Payne, M., & Holman, M. J. 2013, ApJ, 774, 129. doi:10.1088/0004-637X/774/2/129
  * Deck & Batygin (2015) Deck, K. M. & Batygin, K. 2015, ApJ, 810, 119. doi:10.1088/0004-637X/810/2/119
  * Dai et al. (2023) Dai, F., Masuda, K., Beard, C., et al. 2023, AJ, 165, 33. doi:10.3847/1538-3881/aca327
  * Dermott (1968a) Dermott, S. F. 1968a, MNRAS, 141, 349. doi:10.1093/mnras/141.3.349
  * Dermott (1968b) Dermott, S. F. 1968b, MNRAS, 141, 363. doi:10.1093/mnras/141.3.363
  * Goldberg & Batygin (2021) Goldberg, M. & Batygin, K. 2021, AJ, 162, 16. doi:10.3847/1538-3881/abfb78
  * Goldberg & Batygin (2022a) Goldberg, M. & Batygin, K. 2022, AJ, 163, 201. doi:10.3847/1538-3881/ac5961
  * Goldberg et al. (2022b) Goldberg, M., Batygin, K., & Morbidelli, A. 2022, Icarus, 388, 115206. doi:10.1016/j.icarus.2022.115206
  * Goldreich (1965) Goldreich, P. 1965, MNRAS, 130, 159. doi:10.1093/mnras/130.3.159
  * Goldreich & Tremaine (1980) Goldreich, P. & Tremaine, S. 1980, ApJ, 241, 425. doi:10.1086/158356
  * Goldreich & Schlichting (2014) Goldreich, P. & Schlichting, H. E. 2014, AJ, 147, 32. doi:10.1088/0004-6256/147/2/32
  * Goździewski et al. (2016) Goździewski, K., Migaszewski, C., Panichi, F., et al. 2016, MNRAS, 455, L104. doi:10.1093/mnrasl/slv156
  * Hadden & Lithwick (2016) Hadden, S. & Lithwick, Y. 2016, ApJ, 828, 44. doi:10.3847/0004-637X/828/1/44
  * Hadden (2019) Hadden, S. 2019, AJ, 158, 238. doi:10.3847/1538-3881/ab5287
  * Henrard (1982) Henrard, J. 1982, Celestial Mechanics, 27, 3. doi:10.1007/BF01228946
  * Henrard & Lemaitre (1983) Henrard, J. & Lemaitre, A. 1983, Celestial Mechanics, 30, 197. doi:10.1007/BF01234306
  * Huang & Ormel (2023) Huang, S. & Ormel, C. 2023, arXiv:2302.03070. doi:10.48550/arXiv.2302.03070
  * Izidoro et al. (2017) Izidoro, A., Ogihara, M., Raymond, S. N., et al. 2017, MNRAS, 470, 1750. doi:10.1093/mnras/stx1232
  * Izidoro et al. (2021) Izidoro, A., Bitsch, B., Raymond, S. N., et al. 2021, A&A, 650, A152. doi:10.1051/0004-6361/201935336
  * Jontof-Hutter et al. (2016) Jontof-Hutter, D., Ford, E. B., Rowe, J. F., et al. 2016, ApJ, 820, 39. doi:10.3847/0004-637X/820/1/39
  * Kajtazi et al. (2023) Kajtazi, K., Petit, A. C., & Johansen, A. 2023, A&A, 669, A44. doi:10.1051/0004-6361/202244460
  * Lee & Peale (2002) Lee, M. H. & Peale, S. J. 2002, ApJ, 567, 596. doi:10.1086/338504
  * Luger et al. (2017) Luger, R., Sestovic, M., Kruse, E., et al. 2017, Nature Astronomy, 1, 0129. doi:10.1038/s41550-017-0129
  * Masset et al. (2006) Masset, F. S., Morbidelli, A., Crida, A., et al. 2006, ApJ, 642, 478. doi:10.1086/500967
  * Mestel (1963) Mestel, L. 1963, MNRAS, 126, 553. doi:10.1093/mnras/126.6.553
  * Mills et al. (2016) Mills, S. M., Fabrycky, D. C., Migaszewski, C., et al. 2016, Nature, 533, 509. doi:10.1038/nature17445
  * Murray & Dermott (1999) Murray, C. D. & Dermott, S. F. 1999, Solar system dynamics. Cambridge, UK: Cambridge University Press, ISBN 0-521-57295-9
  * Neishtadt (1984) Neishtadt, A. I. 1984, Prikladnaia Matematika i Mekhanika, 48, 197
  * Nesvorný & Morbidelli (2012) Nesvorný, D. & Morbidelli, A. 2012, AJ, 144, 117. doi:10.1088/0004-6256/144/4/117
  * Nesvorný et al. (2022) Nesvorný, D., Chrenko, O., & Flock, M. 2022, ApJ, 925, 38. doi:10.3847/1538-4357/ac36cd
  * Papaloizou & Larwood (2000) Papaloizou, J. C. B. & Larwood, J. D. 2000, MNRAS, 315, 823. doi:10.1046/j.1365-8711.2000.03466.x
  * Peale (1986) Peale, S. J. 1986, IAU Colloq. 77: Some Background about Satellites, 159
  * Petit et al. (2017) Petit, A. C., Laskar, J., & Boué, G. 2017, A&A, 607, A35. doi:10.1051/0004-6361/201731196
  * Petit et al. (2020) Petit, A. C., Petigura, E. A., Davies, M. B., et al. 2020, MNRAS, 496, 3101. doi:10.1093/mnras/staa1736
  * Pichierri et al. (2018) Pichierri, G., Morbidelli, A., & Crida, A. 2018, Celestial Mechanics and Dynamical Astronomy, 130, 54. doi:10.1007/s10569-018-9848-2
  * Pichierri et al. (2019) Pichierri, G., Batygin, K., & Morbidelli, A. 2019, A&A, 625, A7. doi:10.1051/0004-6361/201935259
  * Pichierri & Morbidelli (2020) Pichierri, G. & Morbidelli, A. 2020, MNRAS, 494, 4950. doi:10.1093/mnras/staa1102
  * Pichierri et al. (2022) Pichierri, G., Bitsch, B., & Lega, E. 2022, arXiv:2212.03608. doi:10.48550/arXiv.2212.03608
  * Saad-Olivera et al. (2020) Saad-Olivera, X., Martinez, C. F., Costa de Souza, A., et al. 2020, MNRAS, 491, 5238. doi:10.1093/mnras/stz3369
  * Sessin & Ferraz-Mello (1984) Sessin, W. & Ferraz-Mello, S. 1984, Celestial Mechanics, 32, 307. doi:10.1007/BF01229087
  * 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
  * Terquem & Papaloizou (2019) Terquem, C. & Papaloizou, J. C. B. 2019, MNRAS, 482, 530. doi:10.1093/mnras/sty2693
  * Ward (1997) Ward, W. R. 1997, Icarus, 126, 261. doi:10.1006/icar.1996.5647
  * Wisdom (1986) Wisdom, J. 1986, Celestial Mechanics, 38, 175. doi:10.1007/BF01230429
  * Xu et al. (2018) Xu, W., Lai, D., & Morbidelli, A. 2018, MNRAS, 481, 1538. doi:10.1093/mnras/sty2406