Tidal Synchronization Trapping in Stars and Planets with Convective Envelopes

Tidal torques can alter the spins of tidally interacting stars and planets, usually over shorter timescales than the tidal damping of orbital separations or eccentricities. Simple tidal models predict that, in eccentric binary or planetary systems, rotation periods will evolve toward a “pseudosynchronous” ratio with the orbital period. However, this prediction does not account for “inertial” waves that are present in stars or gaseous planets with (i) convective envelopes and (ii) even very slow rotation. We demonstrate that tidal driving of inertial oscillations in eccentric systems generically produces a network of stable “synchronization traps” at ratios of orbital to rotation period that are simple to predict but can deviate significantly from pseudosynchronization. The mechanism underlying spin synchronization trapping is similar to tidal resonance locking, involving a balance between torques that is maintained automatically by the scaling of inertial mode frequencies with the rotation rate. In contrast with many resonance locking scenarios, however, the torque balance required for synchronization trapping need not drive mode amplitudes to nonlinearity. Synchronization traps may provide an explanation for low-mass stars and hot Jupiters with observed rotation rates that deviate from pseudosynchronous or synchronous expectations.


INTRODUCTION
Tidal torques typically affect the spins of tidally interacting bodies more rapidly than their orbital motion, because the orbital angular momentum usually dwarfs constituent spin angular momenta.Unravelling the effects of tides on spin evolution therefore stands as a prerequisite to predicting their impact on orbital semimajor axes, eccentricities, or spin-orbit misalignments.Simple tidal models invoking a constant time lag in the tidally perturbed body's response to every component of a Fourier-decomposed tidal potential predict that in eccentric systems, tidally interacting planets or stars will evolve toward "pseudosynchronous" states in which the ratio of the spin frequency Ω to the orbital mean motion Ω o approaches a value that depends only on the eccentricity (Hut 1981).Pseudosynchronization or synchronization of orbital and rotation periods are therefore commonly assumed as initial conditions in both theoretical and observational studies of eccentricity or obliquity damping by tides (e.g., Goodman & Dickson 1998;Ogilvie & Lin 2007;Leconte et al. 2010;Vick et al. 2019;Barker 2022;Rice et al. 2022).However, both numerical calculations (Townsend & Sun 2023) and some observations (de Wit et al. 2016;Lurie et al. 2017) suggest that this assumption may not always be warranted.
We demonstrate that, in stars or gaseous planets with large convective regions, dynamical driving of oscilla-tions intrinsically linked to rotation (in particular, inertial modes) by eccentricity tides can hinder tidal evolution toward pseudosynchronization.Sufficiently large eccentricity aliases inertial mode resonances to tidal frequencies at which they can oppose and balance the torque from nonresonantly driven fundamental oscillations (the "equilibrium" tide), introducing a network of spin frequencies Ω at which the time derivative Ω vanishes.Because inertial modes have frequencies that are directly proportional to the rotation rate of a star or planet, a subset of these spin equilibria constitute stable fixed points in the parameter Ω/Ω o , attracting rotation rates toward particular ratios with the mean motion.We predict that many lower-mass stars and/or hot jupiters may be caught in these synchronization traps until eccentricities are damped away, with their rotation rates held at values that deviate significantly from pseudosynchronization.
This paper is structured as follows: Section 2 describes the basic mechanism underlying synchronization trapping, and Section 3 presents a quantitative demonstration of its efficacy.Section 4 then provides a discussion of astrophysical relevance and a comparison with observations, while Section 5 concludes.The appendices provide further details related to the derivation of the secular tidal evolution equations and to the numerical calculations described in Section 3.

SYNCHRONIZATION MECHANISM
Consider a rotating star with mass M , equatorial radius R, and equilibrium density ρ 0 = ρ 0 (r, θ), which is tidally perturbed by a secondary mass M 2 .We assume that the orbit of the two bodies is aligned with the spin of the primary tidally perturbed star, that the secondary is distant enough to be treated as a point mass, and that the tidal potential produced by this point mass can be approximated as purely quadrupolar.Under the additional assumption that the primary rotates rigidly, Appendix A reviews the derivation of secular equations for the time evolution of the rotation rate Ω, the orbital semimajor axis a, and the orbital eccentricity e.The spin evolution is governed by while the eccentricity and semimajor axis evolve according to Equation A14 and Equation A12, respectively.Here, I = V ρ 0 R 2 dV is the primary's moment of inertia, Ω o = [G(M + M 2 )/a 3 ]1/2 is the mean motion, X ℓmk are Hansen coefficients (Equation A4), and κ ℓmk = ℑ[k ℓm ] are the imaginary parts of tidal Love numbers k ℓm = k ℓm (kΩ o ) describing the response of the star in spherical harmonic degree ℓ, azimuthal wavenumber m, and tidal frequency kΩ o to a tidal potential with the form ℜ[A ℓmk (r/R) ℓ Y ℓm e −ikΩot ] (see Appendix A.1; Ogilvie 2014).Spin evolution takes place on a much shorter timescale than semimajor axis or eccentricity evolution, because the orbital angular momentum of the binary dwarfs the spin angular momentum of the primary.Quantitatively, the right-hand sides of Equations ( 1), (A12), and (A14) involve similar coefficients, but with prefactors in the dimensionless ratio where q = M 2 /M is the mass ratio, J = V ρ 0 R 2 ΩdV is the spin angular momentum, and P o = 2π/Ω o is the orbital period.For this reason, studies of tidal evolution in (for example) solar-mass binary star systems commonly assume that spin synchronization can be taken as a foregone conclusion prior to considering the effects of tides on a and e (e.g., Goodman & Dickson 1998;Barker 2022).We do not necessarily disagree.However, self-consistently including the effects of stellar rotation on the dissipative coefficients κ ℓmk can lead to equilibrium synchronization states (with Ω = 0) that deviate significantly from the pseudosynchronous rotation rate predicted by Hut (1981).

Synchronization equilibria
To search for equilibrium rotation states with Ω = 0, first note that as e → 0, X 222 → 1, X 221 → −e/2, and X 223 → 7e/2, while X 22k ∼ O(e 2 ) for all other k.Then to second order in eccentricity, Strong nonmonotonic variation in the κ ℓmk can introduce additional rotation rates at which Ω = 0, though, specifically when tidal frequencies come into resonance with the natural frequencies of the seismic oscillation modes of the tidally perturbed body.Each κ ℓmk can be computed from the properties of oscillations via (e.g., Dewberry & Wu 2024) where α labels modes, and ∆ αk = ω α − ω k is the detuning between the tidal frequency ω k = kΩ o − mΩ and the frequency ω α of a given mode α (both considered in the rotating frame).Meanwhile ϵ α and Q α ℓm are constants characterizing the mode energy and spatial overlap with the tidal force (see Appendix C).Lastly, γ k = γ(kΩ o ) is the (tidal frequency dependent) damping rate experienced by all of the collectively driven oscillations due to internal processes. 1We refer to all quantities appearing in Equation 5 in a dimensionless form scaled by units with Away from resonances, Equation 5 is dominated by monotonic contributions from nonresonantly driven "fundamental" (f-)modes.In the absence of any other oscillation modes, nonresonant f-mode driving predicts 2.0 2.5 3.0 3.5 4.0 / o = P orb /P rot 10 4 10 5 10 6 10 7 10 8 10 9 0 10 9 10 8 10 7 10 6 10 5 10 4 Figure 1.Profiles of Ω/Ω (right-hand side of Equation 1) computed for an n = 1.5 model of a fully convective, rotating star (with rotation rate Ω ≃ 0.09Ω d , where Ω d is the star's dynamical frequency), which is tidally interacting with an equal mass body through an orbit with eccentricity e = 0.5, and a varying orbital mean motion Ωo (the y-axis switches from log to linear scale at 10 −9 ).We adopt a constant kinematic viscosity ν with Ekman number E k = ν/(R 2 Ω) = 10 −5 .The dark blue solid line plots a profile computed solely with f-modes, and roughly describes the "equilibrium tide" of weak frictional theories (Ogilvie 2014).The light blue line plots a profile of Ω/Ω computed with the inclusion of inertial modes, which introduce additional rotational equilibria (where Ω = 0).The arrows along the curves indicate the expected evolution toward (away from) the stable (unstable) fixed points indicated by the filled (unfilled) black points.
a single pseudosynchronous state in rough agreement with constant time-lag models (see the dark blue solid line in Fig. 1).On the other hand, small detunings ∆ αk (combined with weak damping and large ratios |Q α ℓm | 2 /ϵ α ) can introduce additional zero-crossings for Ω curves computed for a given set of orbital parameters, when the nonresonant f-mode torque is balanced by a near-resonant torque contribution of the opposite sign (see the light blue solid line in Fig. 1).
If the dependence of the mode's rotating frame frequency ω α on Ω is known, the resonance condition ω α −ω k ≃ 0 in principle predicts these additional equilibrium rotation states.The rotating-frame frequencies of, for example, gravito-inertial oscillations can be written as ω α ≃ ω α0 +β α Ω, where ω α0 is the mode's frequency in the nonrotating limit and β α is a constant that depends on internal structure.Balance with ω k would then imply the existence of equilibrium rotation states close to the ratio Equation 6 is less useful than Equation 4, though, since it predicts a rotation rate Ω that does not scale directly with the mean motion Ω o ; as the orbit changes, so too will the ratio Ω/Ω o .Both ω α0 and the coefficient β α are also subject to change with the internal structure of the tidally perturbed body as it evolves.
Several previous works have noted that eccentricity can introduce multiple possible spin equilibria (e.g., Storch & Lai 2014).The main point of this paper is to highlight the fact that the dynamics simplify when oscillations possess frequencies that vanish in the limit of zero rotation (ω α0 = 0).In particular, planets and stars with convective envelopes possess long-wavelength inertial modes (e.g., Wu 2005) or mode-like flows (Papaloizou & Ivanov 2010; Lin & Ogilvie 2021) that enhance dissipation at frequencies that are (at least in the limit of slow-to-moderate rotation) linearly related to Ω; inertial mode frequencies satisfy ω α ≈ β α Ω, where β α depends only weakly on internal structure (e.g., Dewberry & Lai 2022).We find that, in stellar models with a large convective envelope and for eccentricities ≳ 0.1, resonances with the pair of prograde and retrograde m = 2 inertial modes with the longest wavelengths (for which we write β α = β ± ) introduce numerous stable synchronization states close to the ratios The stable subset of these equilibria are indicated by the black filled points plotted at some of the light blue curve's zero-crossings in Figure 1.The key point is that, because inertial-wave frequencies scale directly with Ω, these synchronization states remain fixed near the ratio of Ω/Ω o given by Equation 7 as the rotation rate and orbital parameters evolve (see Figure 2), although they evaporate for sufficiently small eccentricity.Equation 7is even relatively robust to changes in internal structure (save for the complete disappearance of a convective envelope), since inertial waves' frequencies are determined primarily by the Coriolis force (this is also the case for Rossby modes; Papaloizou & Savonije 2023).

Mode amplitudes during torque balance
In Section 3 we demonstrate quantitatively that synchronization traps near the period ratios given by Equation 7 can impede progression toward pseudosynchronous rotation in stars and planets with substantial convective envelopes.First, though, in this section we explore predictions for the amplitudes achieved by modes involved in such a torque balance.Estimates of mode amplitudes are necessary to determine whether or not the linear approximation assumed throughout this paper will remain reasonable throughout tidal evolution.
Equation 7 provides only a rough estimate of the period ratios where we expect to find equilibrium rotation states; to estimate the mode amplitudes required to trap the primary's rotation rate, we search for states in which Ω = 0 due to a balance between the torque contribution associated with one mode (labeled by α ′ ) and torques attributed to other processes.Equation 1 and Equation 5then imply that during such a balance, the frequency detuning associated with the mode α ′ satisfies where is the contribution from nonresonantly driven f-modes (in principle, τ f could be modified to include torques from other processes such as magnetic braking), and j is the integer such that |∆ α ′ j | = |ω α ′ − ω j | ≪ 1 is smaller than the detuning associated with any other k-harmonic (for the mode α ′ ).For this quadratic equation to have real roots, ϵ α ′ τ f γ j must be both negative and small enough for the first term on the right-hand side to be larger than one.The dimensionless mode amplitude achieved during this balance can be estimated from (e.g., Dewberry & Wu 2024) where Ã22j = A 22j /(GM/R), and A 22j is given by Equation A2.Making use of Equation 8, we then find where Ĩ = I/(M R 2 ) and Ω = Ω/Ω d are the primary's scaled moment of inertia and rotation rate, Ω d = (GM/R 3 ) 1/2 is its dynamical frequency, tΩf = t Ωf Ω d , and gives the timescale for spin evolution due solely to nonresonant f-mode driving (i.e., equilibrium tidal torques).
The mode amplitude required for spin equilibrium therefore scales directly with the nonresonant tidal torque encapsulated in τ f ; in order to balance a stronger nonresonant torque, we might expect that a mode would need to be driven further into resonance.However, Equation 10 also involves an inverse dependence on the damping rate γ j , which balances the implicit dependence of τ f on internal dissipation (through Equation 5).Equivalently, smaller damping rates γ j due to weaker dissipation are compensated by longer timescales t Ωf .For spin evolution governed primarily by tides, we find that mode amplitudes in a synchronization trap remain small (see Figure 3) and essentially independent of the assumed viscosity (see Appendix E).

Interior models
As a proof of concept, we characterize fully selfconsistent tidal spin synchronization in polytropic stellar models with pressure and density related by P = Kρ 1+1/n , for n = 3/2.Together with a first adiabatic exponent Γ 1 = 5/3, this choice of index n produces models that are neutrally stratified (i.e., fully convective).Numerically computed polytropic models are characterized by dimensionless "eigenvalues" , where ρ c is the central density (e.g., Boyd 2011;Dewberry et al. 2022).To construct a sequence of rotating models representative of solarmass stars very early on the pre-main sequence, we first choose a primary mass M = 1M ⊙ and a nonrotating radius of R = 2R ⊙ .This fixes K, through λ.Assuming that this constant remains unchanged with changes in spin, the equatorial radius for a model with a given rotation rate can then be computed from where M = M/(ρ c R 3 ) is a dimensionless mass.Appendix B provides further details about models computed with rotation rates up to nearly the dynamical frequency Ω d = (GM/R 3 ) 1/2 , and compares the internal structures recovered to that of a 1M ⊙ , pre-main sequence stellar model computed with MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019;;Jermyn et al. 2023).

Tidal calculations
We compute the dissipative coefficients κ 22k for our sequence of rotating stellar models using an aggressively truncated expansions in m = 2 modes that consist solely of the prograde and retrograde sectoral (ℓ ≈ m = 2) fmodes, and the two (prograde and retrograde) m = 2 inertial modes with the longest wavelengths (see Figure 1 in Dewberry & Lai 2022).Appendix C describes the numerically computed mode properties that enter into the calculation, and compares profiles of κ 222 = ℑ[k 22 (2Ω o )] computed for our sequence of rotating polytropes (via Equation 5) to direct tidal calculations for a MESA model of a 1M ⊙ star on the pre-main sequence.

Secular evolution integration
Adopting the stellar and tidal models described in Section 3.1 and Section 3.2, we integrate Equation 1 while assuming a fixed orbital eccentricity and semimajor axis.Equation 1 can be written in dimensionless form as where again t = tΩ d , Ω = Ω/Ω d , Ĩ = I/(M R 2 ), and ã = a/R.Note that Ĩ and ã are both functions of Ω (even for a fixed physical semimajor axis), since the equatorial radius R = R(Ω).We solve Equation 13using the Gragg-Bulirsch-Stoer integrator included in REBOUND (Rein & Liu 2012;Rein & Spiegel 2015), self-consistently computing Ĩ, ã, and κ 22k for all k such that |X 22k | 2 > 10 −6 (decreasing this tolerance to 10 −10 had no effect on our results); at each time-step, we interpolate from the sequence of model and mode data described in Appendix B and Appendix C (respectively) in order to generate the coefficients appearing on the right-hand side of Equation 13.

Results
Figure 2 illustrates the spin evolution of a star in an equal-mass system with (i) an eccentricity of e = 0.5, (ii) a constant Ekman number of E k = ν/(R 2 Ω) = 10 −5 (here ν is a kinematic viscosity), and (iii) rotation rates (indicated by the colorbar) ranging from Ω ≈ 0.02Ω d to 0.23Ω d (rotation periods ranging from P rot ≈ 20.5 days to 1.5 days).The curves in the top panel plot the inverse of rotational evolutionary timescales t Ω defined by t −1 Ω = Ω/Ω (in units of inverse years) as a function of Ω/Ω o = P orb /P rot , with darker (lighter) colors indicating slower (faster) rotation rates Ω.Zero-crossings of these curves delineate equilibrium states with Ω = 0.The red dot in Figure 2 (top) indicates the pseudosynchronous equilibrium state predicted by Hut (1981).
As discussed in Section 2.1, calculations including only f-modes (see the dark blue line in Figure 1) predict vanishing torques at only a single period ratio that generally falls close to this prediction, except when rotation becomes fast enough and/or eccentricity large enough for f-modes to be driven resonantly.The inclusion of the two longest-wavelength m = 2 inertial modes leads to more interesting predictions for spin evolution, producing numerous zero-crossings in Ω close to the ratios of Equation 7(indicated by orange dashed lines for the prograde inertial mode, and blue dotted lines for the retrograde mode).The black points indicate predictions for the zero-crossings close to these ratios, computed using the frequency detuning estimate given by Equation 8. Small deviations between these predictions and the actual zero-crossings of Ω occur when the f-mode torque is balanced by more than the driving of an inertial mode by a single Fourier component of the tidal potential.
Importantly, the scaling of inertial mode frequencies with rotation rate ensures that these zero-crossings occur at ratios of orbital to rotation period that do not shift with changing rotation period (i.e., the light and dark colored lines cross zero at the same values of Ω/Ω o , regardless of the individual values of Ω or Ω o ).The consequence of this commonality is that every ratio Ω/Ω o = k/(2 + β ± ) involving an inertial mode driven strongly enough by the k'th tidal component to balance the nonresonant (f-mode) torque introduces a pair of zero-crossings in Ω that constitute fixed points in the spin evolution.
The fixed points involving a negative slope in Ω/Ω with respect to Ω/Ω o are stable (since a positive Ω will push the spin evolution toward larger values of Ω; see the arrows in Figure 1), and the bottom panel of Figure 2 illustrates the efficacy of these stable fixed points in "trapping" the primary body's rotation rate.The curves in this panel plot the same ratio Ω/Ω o (x-axis), computed as a function of time (y-axis) by integrating Equation 13while assuming a fixed orbital period of 10 days.These curves indicate first of all that the torques plotted in the top panel act to drive the star's spin toward the pseudosynchronous state predicted by Hut (1981).Along the way, though, the inertial mode resonances invariably halt progression toward pseudosynchronization at ratios of Ω/Ω o just above (below) k/(2 + β ± ) if the nonresonant f-mode torque is negative (positive).
The time taken for rotation rates to become stuck in one of these synchronization traps decreases with increasing eccentricity (see Appendix D), and with increasing viscosity (see Appendix E).Because we have for simplicity assumed a constant Ekman number E k = ν/(R 2 Ω), the latter point explains why Figure 2 indi-  (1981).From dark to light, the colormap indicates increasing rotation rate Ω.Because inertial-wave frequencies scale directly with Ω, the values of Ω/Ωo where Ω vanishes do not vary with changes in rotation rate.Bottom: ratios Ω/Ωo (x-axis) as a function of time (y-axis) computed by directly integrating Equation 13 for a fixed P orb = 10 days, and initial rotation rates Ω indicated by the color-scale.The period ratio-invariant zero-crossings in Ω due to inertial modes (top panel) impede the star's progression toward the pseudosynchronization.
cates more rapid spin evolution in stars that are initially spinning more rapidly.The ability of an inertial mode resonance to act as a synchronization trap is relatively insensitive to the nature and amplitude of viscous dissipation, though, requiring only that the near-resonant response be large enough to balance the (typically small) nonresonant f-mode torque.Indeed, we find that the period ratio where Ω vanishes (as determined from Equation 8, or from the numerically computed curves) is essentially invariant with Ekman number.This invariance (demonstrated by Figure 4) derives from the fact that although tidal resonances generally become narrower with decreasing viscosity (see Figure 9), a more weakly damped inertial mode only has to balance the nonresonant f-mode torque (which also weakens with decreasing viscosity).
Because the torque balance required for a synchronization trap does not drive inertial modes deeply into resonance, mode amplitudes need not be large.Figure 3 (left) plots the amplitudes predicted by Equation 10, computed from values of τ f and R evaluated at the end of each of the integrations shown in Figure 2 (bottom).The right-hand panel then compares these predictions to the mode amplitudes actually computed during the integrations, demonstrating agreement within ∼ 10%.For the long-wavelength inertial modes considered here, the dimensionless amplitudes of |c α | ∼ 10 −4 to 10 −3 shown in Figure 3 8in the lowest-viscosity case.These curves demonstrate that the period ratios where the total torque vanishes do not vary with viscosity, except when that viscosity is very large.
In this respect synchronization trapping differs from resonance locking (e.g., Fuller 2017), which involves continuous resonant tidal wave excitation that can often drive modes to nonlinear amplitudes (Ma & Fuller 2021); our spin equilibria are determined by a balance that causes the total torque to vanish, rather than a con-tinuous transfer of energy into mode driving.Spin regulation by synchronization trapping also promises to be more or less automatic (barring extreme changes in internal structure), since it does not rely on the coincident evolution of tidal and mode frequencies (just the continued existence of a large convective region).

Astrophysical relevance
The astrophysical relevance of the physical mechanism described in this paper depends on the timescales over which tidal synchronization takes place, which in turn depend on the highly uncertain interaction of tidal and convective flows.By adopting a constant Ekman number in our calculations, we omit a reduction in the efficiency of convective damping that is expected with increasing tidal frequency (Zahn 1989;Goodman & Oh 1997;Duguid et al. 2020), as well as the question of whether convective turbulence should even be expected to act as an effective viscosity when tidal frequencies outstrip convective turnover times (Terquem 2021;Barker & Astoul 2021;Terquem & Martin 2021).Similarly, the effect of tides on spin rate may be overcome by additional processes (like magnetic braking; Matt et al. 2015) in many systems, although additional weak torques may only modify the balance between the mode resonance and the nonresonant f-mode torque (i.e., change the location of zero-crossings in Ω).Additionally, our assumption that the tidally perturbed body spins as a rigid body is oversimplified (Goldreich & Nicholson 1989;Astoul & Barker 2022).
However, we contend that-in the event that tidal dissipation does act in eccentric systems to regulate the spins of stars or planets with considerable convective envelopes-synchronization traps produced by inertial modes will hinder progression toward pseudosynchronization.Rossby modes (Papaloizou & Savonije 2023) might also play a similar role in bodies with significant stable stratification, although their ability to regulate spin evolution will likely be outweighed by resonance locks involving gravito-inertial oscillations with frequencies that do not scale directly with the rotation rate, and hence do not attract rotation rates to period ratios that are as simple to predict (Witte & Savonije 1999, 2001).
We consider fully convective models as an approximation to solar-mass stars on the pre-main sequence.The dynamical tidal response will of course vary as the radiative core in such stars grows, but Lin & Ogilvie (2021) demonstrated that the dominant dissipation in models with convective envelopes and rigid cores is still associated with long-wavelength flows closely resembling the inertial modes of fully isentropic bodies (in both spatial structure and frequency).We therefore expect the synchronization traps associated with the two longestwavelength inertial modes to remain relevant.However, both radiative and solid cores generically produce enhanced inertial wave dissipation at many more frequencies within the range [−2Ω, 2Ω] than are relevant in fully convective models (Ogilvie 2009(Ogilvie , 2013;;Lin 2023;Dewberry 2023;Pontin et al. 2023).The presence of a core may therefore lead to synchronization trapping within an even more granular set of spin equilibria.
How spins evolve within the network of rotation equilibria described in this paper will in turn affect the dominant channels for tidal damping of semimajor axes, eccentricities, and spin-orbit misalignments, by altering the tidal frequencies felt in the corotating frame of the fluid (Ogilvie & Lin 2007).Frequency-dependent spin regulation by inertial and/or Rossby waves may therefore be important to consider in attempting to explain the circularization of solar-type binaries (Meibom & Mathieu 2005;Zanazzi 2022).

Comparison with observations
Observations offer some evidence that pseudosynchronization may not operate unimpeded in nature.Figure 5 plots the period ratios of eclipsing binaries obtained by Lurie et al. (2017) against orbital periods (top panel) and eccentricity estimates (bottom panel).We computed these eccentricities by cross-referencing with the Villanova Kepler Eclipsing Binary Catalog (Conroy et al. 2014), and using the relation e cos ϖ = (π/2)(s − 0.5), where ϖ is the pericentre angle relative  to line-of-sight, and s = (t s − t p )/P orb is the difference between primary and secondary eclipses (Zanazzi 2022).Taking these estimates of |e cos ϖ| as an approximation of (in fact a lower bound for) e, point colors in the top panel indicate deviations from the pseudosynchronization expected for a given eccentricity.Although the data cluster around the pseudosynchronous prediction given by Equation 4(red line in the bottom panel), they also exhibit significant departures from this prediction.These deviations could simply originate in binaries that have not had enough time to synchronize.However, the physical mechanism explored in this paper also predicts that trajectories toward pseudosynchronization in eccentric systems would be impeded by the numerous synchronization traps indicated by the vertical lines in the bottom panel.
In stars hosting exoplanets, stellar rotation rates are broadly distributed over a range of asynchronous pe- / d = 0.016 Figure 6.Inverse spin evolution timescales (right-hand side of Equation 1) computed for parameters roughly appropriate to the hot Jupiter HD 80606 b (we assume a rotation rate Ω/Ω d ≈ 0.016, mass ratio M * /Mp = M⊙/(3.94MJ) ≈ 266, a planetary radius Rp = RJ , an eccentricity e = 0.93366, and an Ekman number E k = 10 −5 ).As in Figure 1, the light and dark blue curves plot profiles of Ω/Ω computed both with and without the dynamical tidal effects of inertial modes (respectively).Due to the high eccentricity, we plot these curves as a function of (1 − e) 3/2 Ω/Ωo rather than the period ratio.The red point indicates the pseudosynchronous prediction of Hut (1981), while the vertical black line and gray shading indicate the period ratio and uncertainties reported by de Wit et al. (2016).As shown in the inset, this asynchronous observation could be explained by numerous synchronization traps produced where the retrograde inertial mode causes Ω/Ω to cross zero.
riod ratios (Canto Martins et al. 2023).Synchronization trapping likely is not needed to explain this asynchronicity, though, since many exoplanets are not massive enough to efficiently regulate the spins of their host stars.The exoplanets themselves, on the other hand, ought to be spun up or down more efficiently.The hot Jupiter HD 80606 b (de Wit et al. 2016) is highly eccentric, and possesses a spin period that deviates significantly from the expected pseudosynchronous rotation rate.Assuming that HD 80606 b is largely convective and possesses inertial modes with frequencies that scale similarly to the models considered here, the observed period ratio of P orb /P rot ∼ 29 could be explained by a synchronization trap involving the retrograde inertial mode and k ∼ 25 − 30. Figure 6 plots spin evolution timescales computed for a polytropic model with parameters scaled to match those of HD 80606 b.The spin evolution profiles shown in Figure 6 involve a dense forest of resonances, owing to the system's high eccentricity.The retrograde inertial mode generically produces synchronization traps that coincide with the observed period ratio constraints.

SUMMARY
We have demonstrated that in stars and planets with large convective regions, tidal driving of inertial oscillations by an eccentric orbit can introduce numerous spin equilibria beyond the pseudosynchronous state associated with constant time-lag or equilibrium tidal models (Hut 1981).A direct consequence of the fact that inertial-wave frequencies scale with the rotation rate is that a subset of these additional equilibria act as stable fixed points at particular values (given by Equation 7) of the period ratio P orb /P rot .They can therefore inhibit progression toward pseudosynchronization (a process commonly assumed to take place rapidly in tidally interacting systems), "trapping" spins at values that are closer to initial rotation rates.
The mechanism underlying this process of synchronization trapping is similar to resonance locking, except that near-resonances are maintained automatically by the scaling of inertial waves' frequencies with the rotation rate.Additionally, we find that the balance of torques does not drive inertial mode amplitudes to nonlinearity, and that the mode amplitudes and frequency detunings are insensitive to the effective viscosity ascribed to convective turbulence (at least in situations where tides dominate spin evolution).
The calculations presented in this paper are simplified, but find motivation in observations of both stellar binaries and exoplanet systems with rotation rates that do not conform to pseudosynchronous expectations.We also expect synchronization trapping to affect subsequent damping of orbital eccentricities, semimajor axes, and spin-orbit misalignments, by altering the frequencies of tidal driving.Future studies should investigate the interaction of the spin regulation considered here with the secular evolution of the other orbital elements.
Acknowledgements: I thank the reviewer for constructive comments that improved the quality of this paper.I also thank Yanqin Wu, Jim Fuller, and J. J. Zanazzi for helpful conversations.This work was supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) [funding reference #CITA 490888-16].
Software: numpy (Harris et al. 2020), matplotlib Hunter (2007), REBOUND (Rein & Liu 2012), MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019;;Jermyn et al. 2023), CMasher (van der Velden 2020) APPENDIX A. TIDAL EVOLUTION EQUATIONS A.1.Tidal and response potentials Working with spherical coordinates (r, θ, ϕ) in an inertial frame with the origin fixed at the center of a tidally disturbed fluid body that is axisymmetric in isolation (with mass M and equatorial radius R), consider the gravitational potential U = U (r, t) of a tidal perturber with separation d(t), mass M 2 , and negligible spin-orbit misalignment.Assuming that this perturber is sufficiently distant that it can be treated as a point-mass, the tidal potential can be expanded in spherical harmonics as Here Y nm are spherical harmonics, Ω o = [G(M + M 2 )/a 3 ] 1/2 is the mean motion, a is the orbital semimajor axis, and we omit the n = 0, 1 terms responsible for basic Keplerian orbital motion (e.g., Ogilvie 2014).The coefficients A nmk are given by where Here e is the orbital eccentricity, and are Hansen coefficients.We truncate the tidal potential at the quadrupolar spherical harmonic degree n = 2.The tidal force −∇U induces fluid motions, density variations, and consequently an Eulerian perturbation δΦ(r, t) to the equilibrium potential Φ 0 (r) of the tidally perturbed body.The gravitational perturbation can similarly be expanded as We assume that the coefficients Φ ′ ℓmk and A nmk of a given k and m can be linearly related via (complex-valued) potential Love numbers k n ℓm = k n ℓm (kΩ o ): For our purely quadrupolar tidal potential, this linear relation reduces to Φ ′ ℓmk = k 2 ℓm A 2mk .

A.2. Tidal power and torque
The rate of energy transfer in our inertial frame is given by (Ogilvie 2013) where the surface integral is carried out on any spherical surface of radius s > R that contains the primary and excludes the secondary.Inserting the expansions in Equation A1 and Equation A5, the time-steady rate of energy transfer (tidal power) in the inertial frame can be found as With our drastic truncation of the tidal potential at n = 2, the tidal power simplifies to where gives the imaginary part of the Love number k ℓm ≡ k ℓ ℓm describing the response in the same harmonic as the driving.
Similarly, for aligned spin and orbital angular momenta the tidal torque imposed on the perturbed body is given by where δρ is the Eulerian density perturbation.The steady rate of angular momentum transfer due to a quadrupolar tidal potential is then

A.3. Secular equations
The orbital energy is Meanwhile, assuming that the tidally perturbed body rotates rigidly, and equating the rate of angular momentum transfer with the time derivative of the tidally perturbed body's spin angular momentum J = IΩ (here Here we have ignored İ, which is technically nonzero for nonzero Ω (since changes in Ω lead to changes in the equilibrium density distribution ρ 0 ).Finally, writing T = − Lo , where

B. ROTATING MODEL SEQUENCE
Figure 7 summarizes the properties of the sequence of rotating models considered in this study.The lefthand panel plots pressure vs. density (in cgs units) for a nonrotating n = 1.5 polytrope with mass M = 1M ⊙ and radius R = 2R ⊙ (black line), and for snapshots of a MESA model of a 1M ⊙ star during the pre-main sequence (colored lines).The middle panel in Figure 7 plots equatorial radii (blue), rotation periods (orange), and dimensionless moments of inertia as a function of rotation rate for the sequence of n = 1.5 polytropes described in Section 3.1.Lastly, the meridional cuts in the right-hand panel show the surface radius for each rotating model in the sequence.The sequence of polytropic models demonstrated in Figure 7 was computed using a Newton-Kantorovich iteration (as described in Boyd 2011; Dewberry et al. 2022).2This approach yields results in agreement with classic calculations of rapidly rotating polytropic models (e.g., James 1964;Ostriker & Mark 1968;Hachisu 1986).The MESA model was computed with version r15140, using the 1M pre ms to wd test suite example.

C. OSCILLATION MODE PROPERTIES
Figure 8 plots the mode properties required to compute the dissipative coefficients κ 22k (see Equation 5), calculated (as described in Dewberry et al. 2022) as a function of rotation rate for the polytropic model sequence described in Section 3.1 and Appendix B. In particular, ω α is the mode frequency in the corotating frame, and where ξ α is the Lagrangian displacement associated with a given mode α, v α is its velocity perturbation, and ⟨ξ α , ξ β ⟩ = V ρ 0 ξ * α • ξ β dV defines an inner product.The modes of rotating stars are not orthogonal under this inner product, but we normalize each so that ⟨ξ α , ξ α ⟩ = M R 2 .
Figure 9 plots frequency-dependent profiles of ℑ[k 22 ](2Ω o ) = κ 222 computed using these mode properties and a rotation rate of Ω = 0.1Ω d , for different values of Ekman number (indicated by the colormap).These profiles can be compared with direct tidal dissipation rates computed (using the method described in Dewberry 2023) from snapshots of the pre-main sequence MESA model described in Appendix B (filled points).For the latter calculations, we adopt a spatially variable kinematic viscosity ν = ℓ c v c , where ℓ c and v c are the convective mixing length and velocity calculated by MESA.In order to remove discontinuities in the derivative ∂ν/∂r, we smooth the transition of ν to zero inside the radiative core with an exponential decay function.
Figure 9 demonstrates that our reduced polytropic model captures the essential frequency dependence exhibited by the tidal response computed directly from MESA models on the pre-main sequence.The curves for E k = 10 −4 − 10 −3 match the direct calculations most closely.We adopt a more conservative E k = 10 −5 for most of our calculations due to the expectation that short tidal periods should diminish the efficiency of damping by convective turbulence (Zahn 1989;Goodman & Oh 1997;Duguid et al. 2020), and because of uncertainties surrounding the interaction of tidal and convective flows (Terquem 2021;Barker & Astoul 2021;Terquem & Martin 2021).

D. VARYING ECCENTRICITY
The panels in Figure 10 show the same torques and spin evolution as Figure 2, but computed for a 10 day orbit with eccentricities e = 0.1 (top) and 0.3 (bottom).For these lower eccentricities, smaller values of high-k Hansen coefficients produce fewer spin equilibria at which aliased resonances with the inertial modes balance the nonresonant f-mode torque.However, the synchronization traps that do remain still act to inhibit pseudosynchronization.

E. VARYING VISCOSITY
The panels in Figure 11 also replicate those in Figure 2, but for Ekman numbers E k = 10 −4 (top) and 10 −6 (bottom).Fig. Figure 12 plots the same amplitude calculations as Figure 3, but for E k = 10 −4 .The differences in spin evolution primarily come from larger and smaller values for the nonresonant f-mode torques: stronger (weaker) equilibrium tides drive more (less) rapid evolution toward the synchronization traps produced by the inertial oscillations.A qualitative change results from viscosities so large that inertial mode resonances are erased (see the lightest curve in Figure 9, and the synchronization trap near Ω/Ω o ∼ 1 in Figure 11, top).However, as demonstrated in section 3 (see Figs. 3 and 4), for the astrophysically relevant case of weak damping the efficacy of a synchronization trap is effectively independent of the viscosity ascribed to convective turbulence.

Figure 2 .
Figure2.Top: similar to Figure1, the solid curves plot inverse spin evolution timescales for an n = 1.5 polytropic, 1M⊙ star (with a radius R = 2R⊙ in the limit of zero rotation, and an Ekman number of E k = 10 −5 ), which is tidally interacting with an equal-mass companion with an eccentricity of e = 0.5 (the y-axis switches from log to linear scale at 10 −12 ).The black points now show the zero-crossings in Ω predicted by Equation8, and the red point indicates the pseudosynchronous prediction ofHut (1981).From dark to light, the colormap indicates increasing rotation rate Ω.Because inertial-wave frequencies scale directly with Ω, the values of Ω/Ωo where Ω vanishes do not vary with changes in rotation rate.Bottom: ratios Ω/Ωo (x-axis) as a function of time (y-axis) computed by directly integrating Equation 13 for a fixed P orb = 10 days, and initial rotation rates Ω indicated by the color-scale.The period ratio-invariant zero-crossings in Ω due to inertial modes (top panel) impede the star's progression toward the pseudosynchronization.
correspond to radial derivatives of the radial Lagrangian displacement of |∂ξ r /∂r| ≲ 10 −3 (|∂ξ r /∂r| ≳ 1 is typically taken as a criterion for wave breaking, since it indicates fluid velocities in excess of the wave's phase speed; Wienkers & Ogilvie 2018).

Figure 3 .
Figure 3. Left: dimensionless amplitudes predicted by Equation 10 for the end states of all of the integrations shown in Figure 2 (bottom).Circles (squares) indicate spins trapped by torque balances due to the retrograde (prograde) inertial mode (because the integrations all end up with trapped spins, many points lie on top of one another).Right: mode amplitudes computed as a function of time for the numerical integrations shown in Figure 2 (bottom), normalized by the end-state predictions shown in the left panel.The curves demonstrate that (i) Equation 10 accurately predicts the amplitudes of modes involved in a synchronization trap, and (ii) these amplitudes remain linear.

Figure 4 .
Figure 4. Inverse spin synchronization timescales computed for a model with fixed Ω ≈ 0.046Ω d and varying Ekman number (indicated by color-scale), plotted as a function of varying period ratios close to the k = 10 synchronization trap involving the prograde inertial mode.The black points again plot the predictions of Equation8in the lowest-viscosity case.These curves demonstrate that the period ratios where the total torque vanishes do not vary with viscosity, except when that viscosity is very large.

Figure 5 .
Figure 5. Orbital periods (top) and eccentricity estimates (bottom), plotted as a function of period ratio for the sample of eclipsing binaries considered by Lurie et al. (2017).In the top panel, point colors indicate departures from pseudosynchronization (assuming eccentricities e ∼ |e cos ϖ| from the bottom panel).Point colors in the bottom panel indicate orbital period, and in both panels black outlines indicate systems with |e cos ϖ| ≳ 0.1.Discrepancies from the pseudosynchronous prediction of Hut (1981) could be explained by the synchronization traps indicated by the vertical lines in the bottom panel.

Figure 7 .
Figure 7. Left: Pressure vs. density, plotted for a nonrotating n = 1.5 polytropic star with M = 1M⊙ and R = 2R⊙ (black line), and snapshots of a 1M⊙ MESA model on the pre-main sequence (colored lines).Middle: Equatorial radius (blue), rotation period (orange), and dimensionless moment of inertia (green), plotted as a function of Ω/Ω d for the 1M⊙, n = 1.5 polytropic sequence of rotating models incorporated into the integration of Equation 13.Right: meridiodonal cuts showing the surface radii of the polytropic sequence, from Ω/Ω d ≈ 0.016 (dark blue) to Ω/Ω d ≈ 0.95 (bright yellow).