Tidal migration of exoplanets around M-dwarfs: frequency-dependent tidal dissipation

The orbital architectures of short-period exoplanet systems are shaped by tidal dissipation in their host stars. For low-mass M-dwarfs whose dynamical tidal response comprises a dense spectrum of inertial modes at low frequencies, resolving the frequency dependence of tidal dissipation is crucial to capturing the effect of tides on planetary orbits throughout the evolutionary stages of the host star. We use non-perturbative spectral methods to calculate the normal mode oscillations of a fully-convective M-dwarf modeled using realistic stellar profiles from MESA. We compute the dissipative tidal response composed of contributions from each mode as well as non-adiabatic coupling between the modes, which we find to be an essential component of the dissipative calculations. Using our results for dissipation, we then compute of the evolution of circular, coplanar planetary orbits under the influence of tides in the host star. We find that orbital migration driven by resonance locking affects the orbits of Earth-mass planets at orbital periods $P_{\rm orb} \lesssim 1.5$ day and of Jupiter-mass planets at $P_{\rm orb} \lesssim 2.5$ day. Due to resonantly-driven orbital decay and outward migration, we predict a dearth of small planets closer than $P_{\rm orb} \sim 1$ day and similarly sparse numbers of more massive planets out to $P_{\rm orb} \sim 3$ day.


INTRODUCTION
Through missions such as Kepler and TESS, thousands of exoplanets with a variety of sizes, masses, and orbital architectures have been discovered to date around stars of spectral types M through A. Due to the observational bias of both the transit and radial velocity methods, a large proportion of exoplanets are detected at short orbital periods, where their proximity to the host star can shape these planetary orbits through tidal effects.Tidal interactions may have modified the orbital inclinations, eccentricities, and semimajor axes across the history of these planetary systems.
In at least one case, tides are suspected to be acting in real time, as the orbit of exoplanet WASP-12b is measurably decaying from its current orbital period of P orb = 1.09 day according to transit-timing analysis over the past decade (Maciejewski et al. 2013(Maciejewski et al. , 2018;;Patra et al. 2017;Yee et al. 2020).Due to the short decay timescale of ∼3 Myr, it is thought that the influence of the so-called dynamical tide in the late-F star host is accelerating the infall (Weinberg et al. 2017;Bailey & Goodman 2019;Barker 2020).In addition to the equi-Corresponding author: Samantha Wu scwu@astro.caltech.edulibrium tide (Zahn 1966a), which represents the hydrostatic response of the star to its exoplanet perturber, the dynamical tide in the star is the oscillatory component of the response.The oscillation modes of the dynamical tide are damped due to local effects, such as turbulent viscosity in convective zones or radiative damping in stably-stratified stellar regions.This process leads to an exchange of energy and angular momentum between the orbit and the stellar interior.
As the majority of stars in the galaxy are low-mass, the characteristics of M-dwarfs as planetary hosts are important for understanding the effects of tides in the most common exoplanet systems.Though the low luminosity of M-dwarfs makes them intrinsically more difficult to observe, short-period exoplanet systems around M-dwarfs are of interest as they are likely to be located in the habitable zone of the M-dwarfs (e.g., Trappist-1, Gillon et al. 2017) and are relatively easy to detect.Statistical analysis of exoplanet detections indicates that close-in exoplanets around M-dwarfs are much more abundant than around other spectral types, and they are also more likely to be small with R ≲ R ⊕ (Mulders et al. 2015;Hsu et al. 2020;Dressing & Charbonneau 2013;Ment & Charbonneau 2023).With TESS already contributing dozens of exoplanets with measured masses and radii hosted by M-dwarfs and the Nancy Grace Roman Space Telescope anticipated to yield ∼ 10 3 small planets around early-and mid-M dwarfs (Tamburo et al. 2023), the plethora of observed systems will facilitate tests of theoretical predictions for how planetary architectures are shaped by effects such as tides.
Studies of tidal dissipation in low-mass stellar hosts have accounted for the influential contribution of the dynamical tide, which can provide orders of magnitude larger dissipation than the equilibrium tide (Ogilvie & Lin 2007;Bolmont & Mathis 2016).Many also follow how the strength of the tidal dissipation varies during the evolution of the star.In these low-mass stars (M ≲ 0.35 M ⊙ ) that are fully convective during the premain sequence and the main sequence for ∼ 10 11-12 yr, inertial modes (i-modes) excited in the stellar convective zone comprise the dynamical tide.Restored by the Coriolis force, i-modes propagate with rotating-frame frequencies −2Ω s < ω < 2Ω s , where Ω s is the stellar spin frequency, and may be damped by a convective turbulent viscosity.Stars containing radiative zones also experience dissipation due to radiative damping of internal gravity waves.At larger frequencies, the contribution from fundamental modes (f-modes), also described as surface gravity modes, constitutes the majority of the star's equilibrium tidal response.
Previous works have parameterized the strength of the dynamical tide due to i-modes as a function of the stellar properties in order to study the tidal dissipation across stellar evolution.In the simplest case of i-modes excited in a homogeneous envelope surrounding a rigid core with fractional radius α, averaging over the i-mode frequency range −2Ω s < ω < 2Ω s produces a dissipation rate that scales as α 5 and the dimensionless stellar rotation rate ϵ 2 (Goodman & Lackner 2009;Ogilvie 2013).With this formula for the dynamical tidal dissipation as a function of stellar properties, many studies have efficiently calculated the tidal dissipation coupled with stellar evolution (Mathis et al. 2016;Bolmont & Mathis 2016;Gallet et al. 2017;Bolmont et al. 2017).However, the formula that is most commonly used in these works is of limited use for fully convective stars with realistic density profiles, as it predicts negligible tidal dissipation as the core size α → 0. In contrast, tidally-excited inertial modes have actually been demonstrated to contribute significantly to the dynamical tide in rapidly-rotating, coreless isentropic bodies (Wu 2005a,b;Ogilvie 2013;Dewberry & Lai 2022;Dewberry 2023).
Furthermore, studies employing a frequency-averaged formalism do not account for the strong frequency dependence of the dynamical tidal response, which exhibits large peaks associated with resonances with the inertial modes (Ogilvie & Lin 2007;Ogilvie 2013).Ignoring the frequency dependence of tidal dissipation also entails the drawback that no resonance locking can be resolved in these systems.During resonance locking, the enhanced response resulting from tidal excitation of resonant modes persists over a significant period of the stellar lifetime.This mechanism has therefore been invoked to explain faster than expected orbital migration and circularization in stellar binaries and planet-moon systems (Witte & Savonije 1999;Fuller & Lai 2012;Fuller et al. 2016).In order to capture the amplified response near mode resonances, the tidal dissipation as a function of frequency must be calculated throughout the stellar evolution.For convective stars, resonance locking is of particular interest since the frequency range of i-modes coincides with the tidal forcing frequencies of most planetary companions.Thus, planetary orbits will have many opportunities to encounter resonances with the capacity to shape orbital architectures through strong tidal dissipation.
In this work, we perform frequency-dependent calculations of tidal dissipation, coupled with the stellar evolution of a M = 0.2 M ⊙ M-dwarf which is fully convective on the pre-main sequence (PMS) and main sequence (MS).This procedure is repeated for models with different initial rotation periods, allowing us to capture the diversity in outcomes that ensue by varying this initial condition.Using an expansion over the normal modes of stellar oscillation, we compute the dissipative response of the star, thereby resolving the resonant frequencies coinciding with inertial mode (i-mode) and fundamental mode (f-mode) frequencies as a function of stellar age.We use our results to calculate the orbital migration of Earth-mass and Jupiter-mass planets, initialized at a range of ages to represent a few possible formation pathways, and find that resonance locking dominates the migration of planets within ∼ 1-2 day orbits.Section 2 describes how we set up our calculations of the tidal dissipation, normal mode eigenfunctions, and orbital evolution; Section 3.1 demonstrates the dissipative response of the star across its evolution, noting how our approach compares to other methods; and Section 4 shows the associated orbital migration due to tides.We discuss the observational implications for exoplanet demographics and the uncertainties in our implementation of the tidal response in Section 5.

Tidal Formalism
We consider circular, coplanar orbits of a satellite of mass M ′ at a separation a from a body of mass M , which rotates at a spin frequency Ω s .The satellite orbital frequency is Ω o = [G(M + M ′ )/a 3 ] 1/2 .In a frame centered on and rotating with the body of mass M which expe-riences the tidal disturbance, the tidal potential of the satellite can be written as where From here onward, we refer to the perturbed body M as a star, and the satellite M ′ as a planet.The linearized momentum equation for Lagrangian displacement ξ in the star's rotating frame of reference is where C is a self-adjoint operator describing internal fluid forces.The normal modes of the star are ξα = ξ α (r) exp[−iω α t].For a rotating fluid, the normal modes satisfy the following conditions (in units of G = M = R = 1, where M and R are the stellar mass and radius respectively): where ⟨ , ⟩ is the inner product such that Assuming a phase space expansion of the tidal response in terms of normal modes of stellar oscillation, it can be shown (e.g.Schenk et al. 2001;Lai & Wu 2006) that the amplitude of oscillation mode α forced by the tidal potential in the absence of dissipation satisfies where is the tidal overlap coefficient for mode α.In this work, we will consider the effect of forcing by one tidal potential U ℓm at a time.Stationary solutions with ċα = −iω m c α for each mode α are then In the quasi-adiabatic formalism (e.g., Kumar et al. 1995;Burkart et al. 2012), we can write the individual mode damping rates γ α due to viscous damping as (e.g., Ipser & Lindblom 1991) γ α = I αα /(ω α ϵ α ), where where : indicates contraction along both indices of rank-2 tensors.Here µ = ρν is the dynamic viscosity given density ρ and kinematic viscosity ν, and is the shear tensor for mode α with velocity eigenfunction v α = −iω α ξ α .This normalizes the viscous dissipation rate I αα by the total energy of the mode α (Vick & Lai 2020).
In general, the time-averaged dissipation rate due to viscous dissipation of the tide is where the total shear tensor δS is given by Here the velocity fields v constitute the sum of contributions from the normal modes of stellar oscillation: (16) Above, we have modified the mode amplitudes ĉℓ α defined in Equation 11 to include the individual mode damping rates γ α .
An important quantifier of the frequency-dependent tidal response in the perturbed body is the potential Love number k n ℓm (ω m ), which is a dimensionless, complex number.In this notation, for given azimuthal wavenumber m, the Love number describes the effect of an isolated tidal potential of harmonic degree n on the gravitational response of harmonic degree ℓ, at each tidal frequency ω m .Note that in a spherically symmetric calculation, k n ℓm = 0 for ℓ ̸ = n.The imaginary part of the Love number where ℓ = n is related to the energy dissipation rate D by (Ogilvie 2013) in units of G = M = R = 1.We make use of this relation to calculate Im[k ℓ ℓm ], which enters into our orbital evolution equations in 2.2.Note that due to viscous coupling between modes (e.g., Braviner & Ogilvie 2015), values of Im[k ℓ ℓm ] computed from 14 and 17 deviate from results achieved by treating mode damping independently.We discuss this discrepancy further in Section 5.3, and in a companion paper (Dewberry & Wu 2023).

Orbital evolution due to tidal dissipation
Following Ogilvie (2014), the evolution of the orbital semimajor axis a and stellar spin Ω s for circular, copla-nar orbits under tidal dissipation are given by 1 a where the ratio of the orbital angular momentum to the spin angular momentum of the star is, generally, Here Ω o is the orbital frequency, I s is the moment of inertia of the star, and e is the orbital eccentricitythroughout this work e = 0.
We integrate these equations to study the tidal evolution of these quantities.In the calculations described in Section 4, we adopt the approach appropriate for Earthmass planets and ignore the tidal contribution to the spin from Equation 19; in Section 4.1, we discuss the validity of this assumption.
The timesteps required to accurately integrate Equations 18 and 19 are much smaller than the timescales of stellar evolution, and to recalculate the full frequency dependence of Im k 2 2,2 at each step in the orbital evolution is prohibitively expensive.As a result, for each model in a given set of snapshots of stellar evolution, we calculate the imaginary part of the Love number as a function of tidal forcing frequency, Im k 2 2,2 (ω m ).We discuss the details of this calculation and of our stellar models in Section 2.3.Combining the profiles of Im k 2 2,2 (ω m ) for all the stellar models provides a grid of values across tidal forcing frequency and stellar age, for a total of three dissipation grids corresponding to the initial rotation periods of 1, 5, and 10 days.With each pre-calculated grid in hand, at each timestep in the orbital evolution we interpolate over the grid to find the value of Im k 2 2,2 evaluated at the tidal forcing frequency ω m and age of that step.

Normal Mode Oscillations of Realistic Stellar Models
Using MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019)), we simulate the evolution of a star with zeroage main sequence (ZAMS) mass 0.2 M ⊙ , using three different initial rotation periods of 1, 5, and 10 days respectively.We prescribe the rotational evolution of the star due to magnetic braking using the model of Matt et al. (2015).In this formulation of magnetic braking, the torque on the star encounters two regimes: saturated magnetic torque for small Rossby number, which applies below a critical rotation period P sat ≃ 0.1P ⊙ τ cz /τ cz,⊙ ; and unsaturated magnetic torque for slowly rotating stars where magnetic activity strongly correlates with the Rossby number.Here, τ cz is the convective turnover time.The chosen initial rotation periods of 5 and 10 day pass through observed values for young clusters ∼ 5 Myr as well as older clusters at ∼ 500 Myr and a few Gyr (Matt et al. 2015).
In each simulation, we evolve the star from the PMS to the terminal-age MS.For this initial mass, the star remains fully convective until ≈ 70 Gyr, about 20 Gyr before the terminal age MS. Figure 1 shows the evolution of the stellar radius, dimensionless stellar rotation rate, and rotation period for each of the initial conditions.As the stellar radius contracts, the star spins up until ≈ 100 Myr, after which it begins to spin down to ∼100 day rotation periods by 10 Gyr.The dimensionless stellar rotation rate, which is the stellar spin in units of the dynamical frequency Ω dyn = (GM/R 3 ) 1/3 , drops below 10 −3 by ≈ 2.5 Gyr.
Throughout the PMS and MS evolution of our star, we use the 2D pseudospectral methods outlined in Dewberry et al. (2021); Dewberry & Lai (2022); Dewberry (2023) to compute normal modes of stellar oscillation in order to calculate the tidal response.These non-perturbative methods for finding stellar oscillation modes, which can capture the full effects of the Coriolis force by treating the linearized fluid equations as nonseparable in (r, θ) under rotation, assume rigid rotation and zero viscosity.The axisymmetry of the background state allows modes to be labeled by unique azimuthal wavenumbers m.Assuming spherical symmetry, the 1D, radial profiles of realistic stellar structure from MESA constitute the equilibrium profiles for the 2D pseudospectral methods, which we use to solve for adiabatic normal modes with time dependence exp [−iω α t] for rotating frame mode frequencies ω α .
In this work, we solve for modes with azimuthal wavenumber m = 2 and thus consider tidal forcing frequencies ω m = 2(Ω o − Ω s ).We keep track of the ℓ = 2, 4, 6 prograde and retrograde f-modes, where ℓ is the spherical harmonic degree, as well as 12 i-modes.Here, prograde refers to modes whose inertial-frame frequency ω α is positive, and retrograde modes have ω α < 0. For the i-modes, we characterize them by the number of nodes in directions that are roughly horizontal (n 1 ) and vertical (n 2 ), where horizontal describes the direction of the cylindrical radius R, and vertical is the cylindrical z direction (Wu 2005a).We consider the n 1 + n 2 = 1, n 1 + n 2 = 2, and n 1 + n 2 = 3 prograde and retrograde i-modes.To refer to the modes in the text, we follow the notation of Dewberry & Lai (2022) so that f ℓ,m,± denotes the f-modes and i m,n1,n2,± la-bels the i-modes, with + and − indicating prograde or retrograde modes respectively.Inertial modes are only labeled prograde or retrograde if the m, n 1 , n 2 values are duplicated in another mode.The chosen set of i-modes are the longest-wavelength modes; shorter-wavelength modes are unlikely to contribute significantly to tidal dissipation as they have smaller tidal overlap coefficients Q α ℓm and smaller amplitudes.
For each stellar profile, we calculate the dissipation D from Equation 14, with Equations 15 and 16 summed over the oscillation modes α computed for each stellar model.In each sum, we include the f-modes and i-modes listed above.Equation 14 involves a double sum over the mode velocity eigenfunctions that allows us to capture the dissipative coupling between different modes, e.g.I αβ (Equation 12).In most prior work that takes the modal decomposition approach, they instead perform a single sum over each mode to calculate tidal dissipation (using, e.g., Equation 30) that does not incorporate this coupling.However, we find that the terms I α,β can be quite large and contribute significantly to the dissipation away from mode resonances.Given D, the imaginary part of the Love number Im k 2 2,2 as a function of tidal forcing frequency ω m can be found from Equation 17.

Viscous dissipation in convective zones
To treat the viscous dissipation of the tides, we estimate the effective turbulent viscosity of convection throughout the fully convective models.The kinematic eddy viscosity ν of turbulent convective friction can be described using mixing length theory as (Zahn 1966b): where v con , l con , and L con are the convective velocity, mixing length, and convective luminosity respectively, H is the pressure scale height, and ρ is the density.This simple estimate ignores the effect of rotation on convective flows, as well as the attenuation of the dissipation efficiency when the tidal forcing frequency is much larger than the convective turnover frequency, ω m ≫ ω con ∼ v con /l con (see, e.g., Zahn 1966b;Goldreich & Keeley 1977;Duguid et al. 2020a).
Figure 2 shows the variation of l con = H, v con , and the kinematic eddy viscosity ν NR (solid lines in top panel) with radial coordinate, colored by the dimensionless stellar rotation rate for the P rot,0 = 5 day models.While the star is on the pre-MS, the scale height and convective velocity are elevated because the stellar radius and luminosity are much larger than on the MS, leading to higher viscosity on the pre-MS (solid yellow and orange curves in Figure 2).While on the MS, the stellar radius  26) and without (νNR, Equation 21) a rotational correction.As Equation 26is a function of Rossby number, νR varies with rotation and approaches the value of νNR towards the surface of the star where the Rossby number is large.On the main sequence, νNR, H, and vcon remain nearly constant so all lines lie on top of each other under the purple line.
remains nearly constant, so the scale height, convective velocity, and thus ν NR all remain similar throughout the rest of the star's lifetime as it spins down to small values of Ω s /Ω dyn (solid purple curves).
However, rapid rotation such as that experienced by our stellar models affects turbulent convective flows, in particular by inhibiting the efficiency of the turbulent friction.Mathis et al. (2016) address how the viscosity in Equation 21 should be modified for the rotating case by scaling the convective velocity and mixing length for rapid rotation to the non-rotating case (Stevenson 1979;Barker et al. 2014): Here Ro = ω con /Ω s is the Rossby number.Mathis et al. (2016) use the above scalings for the rapidly rotating regime of Ro ≲ 0.25, whereas in the slowly rotating regime of Ro ≳ 0.25 they employ the following: (24) With these scaling relations, the kinematic eddy viscosity in the rotating case may be estimated as .
As expected from the rotational evolution of the star (Figure 1), Ro is much smaller throughout the star during the pre-MS than the late MS.From the bottom panel of Figure 2, we see that Ro ≪ 1 in the stellar interior, but increases sharply towards the surface as H drops and v con rises.As a result, the viscosity transitions between the rapidly-rotating and slowly-rotating regimes within the star.The dotted lines in the top panel of Figure 2 illustrate that the effect of rapid rotation in our models is to lower the kinematic viscosity in the majority of the star, as Ro is small except towards the surface.In addition, the rotationally-corrected viscosity now varies on the MS due to the dependence on Ro, so that the slowly-rotating models towards the late MS ≳ 1 Gyr have 1-2 orders of magnitude larger ν R than at peak rotation around ∼ 10-100 Myr.
A second factor not implemented here is the frequencydependent attenuation of the efficiency for rapid tide where ω m ≫ ω con .In the literature, several scalings have been proposed for this reduction, including linear (Zahn 1966b;Penev et al. 2007Penev et al. , 2009) ) or quadratic (Goldreich & Keeley 1977;Goodman & Oh 1997;Ogilvie & Lesur 2012;Duguid et al. 2020a), with some arguing that no reduction is appropriate (Terquem 2023).For i-modes, whose frequencies scale with the stellar rotation rate, the rapid tide attenuation will effectively . Evolution of the inertial-frame mode frequencies σ in units of the dynamical frequency Ω dyn (top panel) and the absolute value of the ℓ = m = 2 tidal overlap coefficients Q ℓ,m (bottom panel) for the set of f-modes and i-modes we include during the pre-MS to MS evolution of the Prot,0 = 5 day (left) and 1 day (right) sets of stellar models.In the bottom panel, a subset of these modes are plotted.For the Prot,0 = 1 day models, i-modes that experience avoided crossings are plotted with a circle marker in the bottom panel.
scale as Ro or Ro 2 depending on which prescription is taken.In the former case, we expect similar effects on the mode damping compared to the considerations of rotating convection described above, but in the case of quadratic scaling, the suppression will be more dramatic than shown here.Furthermore, if both reductions act simultaneously, the viscosity would strongly scale with Ro and decrease even more.
We repeat our calculation of the dissipation as described in Section 2.3 for both assumptions of the kinematic viscosity ν R or ν NR as described above.The fiducial results shown throughout this work use the results from assuming ν R , and we address comparisons between the two assumptions in Sections 3.2.1 and 4.1.1.

Mode Evolution with Stellar Evolution
The top panel of Figure 3 shows the evolution of the inertial-frame mode frequency σ for the set of f-modes and i-modes tracked in this work, for initial conditions P rot,0 = 5 and 1 day.The description for P rot,0 = 5 day is representative of P rot,0 = 10 day as both are calculated for slowly-rotating stars.As the stellar spin increases and decreases (see Figure 1), the i-mode frequencies are generally proportional to the stellar spin and follow the same evolution.For the slower-spinning models, the f-mode frequencies also increase and decrease with stellar spin slightly, but remain fairly constant and well separated from the i-mode frequencies.In the P rot,0 = 1 day model on the right, the i-mode frequencies again increase and decrease with rotation, but here the f-mode frequencies also show a significant increase with rotation rate.As the star spins down to Ω s /Ω dyn ∼ 10 −3 at   2,2 as a function of tidal forcing frequency ωm (in units of the dynamical frequency Ω dyn ).The y-axis is scaled by ωm, and panels are labeled by the value of Prot,0 used to evolve the stellar model.In each panel, the red lines show the dissipation on the pre-MS for the stellar rotation rates Ωs (in units of Ω dyn ) listed in the legend, whereas the blue dotted lines show the same for a model on the MS with a very similar rotation rate.The x-axis limits in each panel are limited to [−2Ωs, 2Ωs] to show the inertial mode range.Each resonant peak is labeled by the mode responsible for the resonance.The dissipation is larger on the pre-main sequence due to larger turbulent convective viscosity.
∼ 2.5 Gyr, the i-mode frequencies similarly approach 0 in both models.
The bottom panels show the absolute value of the ℓ = m = 2 tidal overlap coefficients |Q 2,2 | (Equation 10), which affect the mode amplitudes; for all i-modes, these drop by several orders of magnitude once the star spins down to Ω s /Ω dyn ∼ 10 −3 at a few Gyr.Throughout the evolution, the prograde and retrograde ℓ = 2 f-modes still have much larger |Q 2,2 |, but the next largest values come from the longest-wavelength i modes, e.g.i 2,0,1 and i 2,1,0 , as well as the ℓ = 4, m = 2 f-modes.
In the P rot,0 = 1 day model, the prograde ℓ = 2, m = 2 f-mode frequency crosses those of the highest-frequency i-modes at large rotation rates Ω s ≈ 0.6-0.8Ω dyn .This leads to several instances of an "avoided crossing", during which the f-mode and i-mode mix in character as they approach one another in frequency (see Dewberry & Lai 2022, for a discussion of similar i-mode/f-mode mixing in isentropic polytropes).These avoided crossings occur between the f-mode f 2,2,+ and the i-modes i 2,1,1,+ and i 2,1,2,+ , which are distinguished with dots instead of crosses in the bottom right panel of Figure 3.The i-mode i 2,2,1,+ also approaches f 2,2,+ in frequency near the maximum stellar rotation rate (also shown as dots).These three i-modes show increases in |Q 2,2 | over several orders of magnitude as their mode frequencies approach or cross that of the f-mode.For the i 2,2,1,+ i-mode whose frequency nears, but does not cross, that of the f-mode, the enhancement is reduced yet still visible.Though none of these i-modes achieve |Q 2,2 | values that are comparable to the f-mode's tidal overlap integral, the effect of mode mixing allows these relatively short wavelength i-modes to contribute similar |Q 2,2 | coefficients compared to the longest-wavelength i-modes, whereas they would otherwise have negligibly small |Q 2,2 | in comparison.The consequence, as we discuss below in the context of Figure 4, is to enhance the amplitudes of these modes during these periods of rapid rotation, leading to greater resonant dissipation.

Tidal dissipation coupled with stellar evolution
We compute the frequency-dependent dissipation D from Equation 14 in each stellar model, using the fmodes and i-modes shown in Figure 3 to construct the velocity fields of Equation 16.Models shown in this section use the rotation correction to viscosity of Equation 26 unless specified otherwise.
Figures 4-5 show the positive-definite quantity Im k 2 2,2 /ω m versus ω m at different snapshots in the stellar evolution.In each panel corresponding to different P rot,0 , the spikes that rise orders of magnitude above the baseline dissipation are peaks of resonance   with the stellar oscillation mode as labeled.Modes with shorter peaks are shorter-wavelength and have smaller maximum amplitudes due to damping (Equation 16) and smaller tidal overlap integrals (Equation 10).Note that Im k 2 2,2 /ω m is constant across ω m = 0 because of the scaling by ω m .The value of Im k 2 2,2 actually approaches zero at ω m = 0, which corresponds to the case of corotation (Ω s = Ω o ).
In red, the stellar model shown is on the pre-MS, whereas the blue dashed curve is for a model on the MS with a similar rotation rate (in units of Ω dyn ).Compared to the pre-MS model, the more compact star on the MS has a lower viscosity and therefore smaller offresonant dissipation even for the same dimensionless rotation rate.For the models shown here, the difference amounts to a factor of a few, but similar comparisons involving the more extended star on the early pre-MS can differ by up to an order of magnitude.
Figure 4 focuses on the region −2Ω s < ω m < 2Ω s , which is the range over which inertial modes are excited.For P rot,0 = 1 day at the rotation rate of Ω s /Ω dyn ∼ 0.8 shown, the peak for f 2,2,+ falls within this range and is very near in frequency to i 2,2,1+ and i 2,1,1,+ .In the vicinity of this rotation rate, the avoided crossings between the f-mode and i-modes as described in the context of Figure 3 amplify the i-mode dissipation.In contrast, the same two i-modes in slower-rotating models appear with some of the smallest resonant spikes, as they have quite short wavelengths and small |Q 2,2 | values.Throughout the evolution, all the modes we consider for P rot,0 = 1 day augment the dissipative response significantly near resonance, but fewer modes contribute significantly for the slower-rotating stars.This corresponds to overall smaller magnitudes of |Q 2,2 | exhibited by the slower-rotating models.
Figure 5 shows a larger range in tidal forcing frequency for the P rot,0 = 5 day model.In this region, we can see that the f-mode amplitudes dominates the dissipation away from resonant peaks, which can also be thought of as the equilibrium tidal response.As the initial rotation period increases, the dissipation away from resonance mildly increases due to the rotational correction to viscosity.The largest dissipation comes from ℓ = m = 2 f-modes and the longest-wavelength i-modes (i 2,1,0 and i 2,0,1 ).In older models with slower rotation rates, the strength of resonant i-mode dissipation fades as the star spins down greatly on the MS, tracking the rapid decline of Q 2,2 with age as seen in Figure 3.The ℓ = m = 2 f-mode peaks maintain similar heights throughout the evolution.
As the star spins up from a few-100 Myr, the viscosity grows more inefficient due to the decreasing Rossby number, and the dissipation away from resonance correspondingly decreases.Once the star spins down significantly, the Rossby number and consequently the viscosity grow.The non-resonant dissipation is similar on the pre-MS at ≲ 10 Myr and on the late MS at ∼ 2 Gyr, even though the star is spinning 10-100 times faster on the pre-MS.Since the fully convective star has a much larger radius on the pre-MS, the convective velocity and scale height are correspondingly larger, leading to larger viscosities on the pre-MS for a given rotation rate.

Comparison of different assumptions for the viscosity
We compare in Figure 6 the effects of including or excluding the correction to viscosity under rapid rotation.In red, the quantity Im k 2 2,2 /ω m is shown for the same pre-MS model as in Figures 4 and 5, for which we compute the dissipation assuming the kinematic viscosity to be ν R (Equation 26).In blue, the same model is shown for a calculation using ν NR (Equation 21) instead.As seen in Figure 2, ν NR ≫ ν R throughout the majority of the star, and this impacts the dissipation curves in two significant ways.The off-resonance dissipation is orders of magnitude larger when using the ν NR instead of ν R , 0.2 0.0 0.2 m ( dyn ) 10 9 10 7 10 5 10 3  2,2 ] scaled by ωm as a function of tidal forcing frequency ωm is shown, but now to compare the effect of viscosity.The red line shows the dissipation using the correction to viscosity due to rotation (Equation 26).The blue dotted lines show the same but using viscosity without the rotation correction (Equation 21).Including the effects of rotation lowers the viscosity, leading to sharper resonant peaks, but smaller off-resonant dissipation by 1-2 orders of magnitude.
which follows from the attenuation of ν R at small Ro by ≈ 2-3 orders of magnitude.Furthermore, the resonant spikes are wider and shorter for the higher viscosity, ν NR .As discussed in Section 4.1.1,we expect larger offresonance dissipation to cause more rapid orbital migration as a result of the equilibrium tide, whereas smaller peak heights decrease the likelihood of resonance locking for orbital evolutionary calculations that utilize the Ro-independent viscosity.We note that if the additional effect of frequency-dependent attenuation of the viscous efficiency were to be included such that the viscosity scaled even more strongly with Ro, then the profiles of Im k 2 2,2 /ω m would trend towards even smaller dissipation away from resonance and narrower, taller peaks.

ORBITAL EVOLUTION DUE TO STELLAR TIDES
To study the orbital dynamics caused by the frequencydependent tidal dissipation coupled with stellar evolu-tion presented in the previous section, we integrate planetary orbits beginning at stellar ages of 10 Myr, 100 Myr, 0.5 Gyr, and 1 Gyr.This in initial ages spans scenarios from in-situ or disk-driven formation of close-in planets, in which planets develop ≲ 10 day orbits while within the proto-planetary disk which persists until ∼ 10 Myr, to dynamically-driven migration, in which planets may reach short-period orbits throughout the stellar lifetime.At each stellar age, we initialize planets at a range of orbital periods from < 1 day to a few day.We repeat this for an Earth-mass companion (M p = M e ) and a Jupiter-mass companion (M p = M j ), and we also complete the entire procedure for each dissipation grid associated with initial stellar rotation rates P rot,0 = 1, 5, and 10 day.For each orbital integration, we assume co-planar, circular orbits and integrate Equations 18-19 from the respective initial stellar age to a final age of 2.5 Gyr.The stellar lifetime is of course much longer than this endpoint, but the stellar spin has slowed so drastically by 2.5 Gyr that our i-mode frequencies approach zero, so beyond this age we expect the f-mode dissipation (i.e., equilibrium tides) to govern the system.
Figure 7 shows a sample of orbital integrations for M p = M e for each initial stellar rotation rate.The top panels show the semi-major axis evolution, while bottom panels show the evolution of the tidal forcing frequency ω m and the rotating-frame mode frequencies ω α .For all values of the initial stellar rotation period P rot,0 , orbits with initial orbital period P orb,0 ≳ 1.5 day experience negligible orbital migration.For these planets that are further away such that the ratio (R/a) is small, ȧ, which depends on (R/a) 5 (Equation 18), remains tiny, and thus the rate of change of the orbital frequency Ωo does not appreciably affect the evolution of the tidal forcing frequency ωm .Nevertheless, the tidal forcing frequency still changes over time for these orbits due to the stellar spin evolution, but it usually traverses a path where the dissipation away from resonance is too low to affect the orbits.
For the shortest initial orbital periods shown of P orb,0 = 0.5 and 1 day, we observe that resonance locking between the tidal forcing frequency and stellar oscillation modes causes significant orbital migration.We first consider the Earth-mass planetary orbits shown in Figure 7.For P rot,0 = 1 day, the star spins faster than all planetary orbits at P orb,0 ≳ 0.5 day so that the planets migrate outward, no matter at what age the integration begins.The planet at P orb,0 = 0.5 day evolves into a brief resonance with the i-mode i 2,0,2 , which is indicated in the bottom panel by the dark blue line overlapping with the dotted pink line (shown more i 2, 0, 2 i 2, 0, 1 i 2, 1, 0 10 7 10 8 10 9 Age (yr) Figure 7. Top row: Evolution of the semimajor axis or orbital period of a Mp = Me Earth-mass companion for the different initial orbital periods shown in the legend in the upper left.Each column depicts results from each stellar evolutionary model with initial rotation periods Prot,0 as labeled, and the dissipation is calculated using the correction to viscosity under rotation (Equation 26).The dashed black line shows the radius of the star.Bottom row: The tidal forcing frequency ωm (in units of the dynamical frequency) for the same models.Dotted lines show the rotating-frame mode frequencies ωα for the i-modes listed in the legend in the center panel.Significant migration occurs when the companion's tidal forcing frequency enters a resonance lock with a stellar oscillation mode, which is visualized as a solid line overlapping with a dotted line in the bottom panels.
clearly in the inset panel).After it breaks out of resonance with i 2,0,2 having migrated slightly outward, it enters another resonance with the longest-wavelength retrograde i-mode, i 2,0,1 , and remains in that resonance until ∼ 2 Gyr.The planet at P orb,0 = 1 day experiences the resonance with i 2,0,1 for a brief interval as well, and its orbit then converges with that of the P orb,0 = 0.5 day planet.By 2.5 Gyr, both of these planets arrive at P orb > 1 day orbits, regardless of the time of formation.
For dissipation in a more slowly rotating star with P rot,0 = 5 day, most orbits begin retrograde (Ω o < Ω s ); however, planets initialized with P orb,0 = 0.5 day at early times ≤ 10 Myr are instead prograde (Ω o > Ω s ).These earlier prograde orbits enter into a resonance lock with the longest-wavelength prograde i-mode i 2,1,0 that lasts from ≈ 10-100 Myr, causing significant orbital decay to P orb ∼ 0.1 day.Eventually, between 700-800 Myr, the planet plunges into the star because of the equilibrium tide, which can be understood equivalently as dissipation due to the ℓ = m = 2 f-mode.If these planets arrive at P orb,0 = 0.5 day between 50-400 Myr, they will instead migrate outward in a resonance lock with i 2,0,1 that breaks at ∼1 Gyr; however, if they land at P orb,0 = 0.5 day from ≈ 500 Myr onward, they are unable to lock into resonance and remain stationary.At P orb,0 = 1 day, orbits remain unaffected by dissipation until ∼ 700 Myr, at which point they will also lock with i 2,0,1 and converge with the P orb,0 = 0.5 day orbits to migrate slightly outward to P orb > 1 day.However, planets arriving at P orb,0 = 1 day after ∼1 Gyr do not undergo migration.
For our slowest set of stellar rotation rates with P rot,0 = 10 day, planets closer than P orb ≲ 1.5 day will migrate inward until ∼100 Myr due to resonance locking  with i 2,1,0 .All these orbits converge to P orb ∼ 0.2 day by 2.5 Gyr due to equilibrium tidal dissipation.Only planets at P orb,0 = 1 day initialized between 50-500 Myr migrate outward.Overall, if Earth-mass planets form short-period orbits of P orb,0 ≲ 1 day at around 10 Myr, then tides tend to pull the planets from this initial position by a few Gyr.However, the final fate (migrating outward, orbital decay, or plunge-in) depends on the initial rotation rate of the star.
Figure 8 depicts the orbital evolution for M p = M j companions.Jupiter-mass planets suffer stronger tidal dissipation because ȧ is proportional to the mass ratio (Equation 18).Thus, even when orbits are not at resonant frequencies, more examples of migration are visible in each panel for P orb,0 < 2 day.In general, the results for Jupiter mass planets are similar to those for Earth mass planets as they lock into resonance with similar modes at similar times, but with more drastic changes in the semimajor axis.Of note are the Jupiter-mass planets that form at P orb,0 = 0.5 day, which mostly plunge towards the star in the P rot,0 = 5 and 10 day examples.One exception is a planet for the P rot,0 = 5 day stellar model that migrates outward in a resonance lock with i 2,1,1,− at ≈ 500 Myr.Planets that migrate outward in resonance locks with i 2,0,2 and i 2,0,1 tend to approach similar 3 day orbits by 2.5 Gyr, independent of the stellar rotation rate.
As noted in Section 2.2, we calculate the change in stellar spin without incorporating this feedback into our orbital integrations, since including this effect is computationally prohibitive at this stage.This assumption holds well for Earth-mass planets as long as the planet is not plunging in, as Figure 9 shows.Even when the planet's tidal forcing frequency passes through resonance with mode frequencies, the ratio |dΩ s,tide /dΩ s,evol | is small such that the change in stellar spin is dominated by the stellar evolution.In contrast, Jupiter-mass planets exert a stronger tidal torque on the star that significantly affects the stellar spin as the planetary orbit passes through resonances with various modes of the star.This introduces appreciable uncertainty in the ability of certain resonance locks with Jupiter-mass planets to hold for the duration presented in this work, which is discussed further in the next  19, Ωs,tide , compared to that due to stellar evolution (magnetic braking, contraction, or expansion), Ωs,evol .This is shown for an Earth-mass planet in blue and and a Jupiter-mass planet in red for the models listed in the legend.Note that the nearly vertical red and blue lines at 250 Myr and 800 Myr respectively are due to plunge-in of the planet (see corresponding times for the same models in top right panel of Figure 8 and top middle panel of Figure 7).section.Moreover, both Jupiter-mass and Earth-mass planets during plunge-in traverse nearly vertical lines in Figure 9 and attain large values of |dΩ s,tide /dΩ s,evol |.This trajectory indicates that both types of planets can greatly increase the stellar spin as their orbits decay into the star, potentially leading to tidal synchronization of the spin and orbit.

Resonance locking
In our orbital evolution calculations, we observe resonance locks that tend to occur between a few-100 Myr for inwardly-migrating orbits, as well as later at a few hundred Myr-1 Gyr outwardly-migrating orbits.To understand why resonance locks occur, we first note the conditions for resonance locking.The tidal forcing frequency must be nearly equal to the frequency of a mode: ω m ≃ ω α .For the lock to be maintained, the time derivatives must remain equal as well: ωm ≃ ωα .Once the resonance lock begins, the companion orbit will remain at a stable fixed point where this condition is true, unless something occurs to break the lock.The orbital migration timescale during the resonance lock follows the evolutionary timescale of the stellar oscillation mode: as the mode frequency changes over time, the tidal forcing frequency follows closely (Fuller et al. 2016).
Recalling that ω m = m(Ω o −Ω s ) ⇒ ωm = m( Ωo − Ωs ), we see that the ability to resonance lock depends on the interplay between the stellar spin evolution and the orbital evolution.In this work, we do not incorporate the (2 + c)M ′ a 2 /3Is = 2. Above this value, resonance locks will break due to tidal torque on the stellar spin, which was not included in these models.
feedback of the tidal dissipation onto the stellar spin, so Ωs is given by the slope of the curves in Figure 1.Tidal dissipation drives the change in semimajor axis ȧ, and consequently Ωo , through Equation 18.Given that the planet is close enough so that (R/a) 5 is non-negligible, then the term Ωo will begin to contribute when the dissipation from Im k 2 2,2 becomes large enough.This will occur near normal mode frequencies where the dissipation rises sharply into peaks.
We can understand migration driven by resonance locking from the condition ωα ≃ ωm .Inertial modes have ω α ≃ cΩ s , where −2 ≤ c ≤ 2 is nearly constant over time.Hence a resonance lock with an inertial mode requires In a resonance lock with an inertial mode, the planet can only migrate inwards if the star is spinning up, and it can only migrate outwards if the star is spinning down, such that Ωo and Ωs have the same sign.We can further decompose the spin evolution into a component from stellar evolution, Ωs,evol , and a component due to tidal torques, Ωs,tide .The latter can be related to the orbital evolution and Ωs,evol via equations 18 and 19.The orbital evolution during the resonance lock can then be written as Ωs,evol .2,2 (positive for prograde orbits ωm > 0, negative for retrograde orbits ωm < 0).The bottom panels show tidal frequency ωm and the rotating-frame mode frequency ωα for the i-mode in the legend.The left column shows an inwardly migrating resonance lock a prograde mode ωα > 0, and the right exemplifies an outwardly migrating resonance lock with a retrograde mode ωα < 0.
The second term in brackets accounts for the backreaction of the tidal torque on the star's spin, which we presently ignore in our orbital evolution calculations.It becomes important when the planet's orbital angular momentum is comparable to that of the star.In particular, as long as (2+c)M ′ a 2 /3I s < 2, then Ωo and Ωs,evol have the same sign and resonance locking therefore continues to be possible.
In Figure 10, the second term in brackets in Equation 28 is shown for some orbits that experience resonance locking.Equation 28 is only valid throughout the duration of the resonance lock, which lasts for ∼100 Myr in each orbit shown (dashed lines).At the start of each resonance lock, Ωo and Ωs are required to have the same sign in order to enter the resonance lock and satisfy Equations 27 and 28.Comparing the right hand sides of Equations 27 and 28 shows that this is equivalent to Ωo and Ωs,evol having the same sign in our models, as (2 + c)M ′ a 2 /3I s < 2 at the start of each resonance lock.For the majority of the models shown, the term (2 + c)M ′ a 2 /3I s never exceeds 2, so the resonance lock may continue.
For the outwardly migrating Jupiter-mass planet whose orbital evolution begins at 50 Myr, the term (2 + c)M ′ a 2 /3I s crosses 2 (dotted gray line) during the resonance lock at 700 Myr, when the planet migrates past P orb ≈ 2 day.Once (2 + c)M ′ a 2 /3I s > 2, the tidal torque on the star will become important and cause the right-hand side of Equation 28 to switch sign.This implies that the star would have to spin up due to stellar evolution, Ωs,evol > 0, in order to maintain the resonance lock condition.Since the star is still spinning down at this point in our stellar models due to magnetic braking, we know the resonance lock condition should actually be broken at this point.As a result, we predict that including the back-reaction of the tidal torque on the stellar rotation will cause resonance locks with outwardly-migrating Jupiter-mass planets to break once they reach P orb ≈ 2 day orbits.In fact, the tidal torque could synchronize the stellar spin with the orbits of these outwardly-migrating Jupiter-mass planets.For Earth mass planets whose resonance locking is treated accurately in our models, Figure 11 shows some examples of the dissipation and frequency evolution during the lock.The left panel exemplifies a resonance lock with a prograde i-mode.The bottom panel shows that the tidal forcing frequency approaches the value of ω α , which is also changing with time.The relevant mode frequency increases as the star spins up until ∼ 100 Myr, then decreases with the stellar spin-down thereafter.Once ω m ≃ ω α , the value of Im k 2 2,2 jumps by 3 orders of magnitude as the planetary orbit encounters the peak in dissipation due to resonance with i 2,1,0 .On the right, Figure 11 depicts a resonance lock for retrograde orbits with ω m < 0 that occurs around a few 100 Myr-1 Gyr for a subset of initial conditions.These planets lock into resonance with i 2,0,1 .
In both cases, the resonance locks break after ∼ 100s of Myr, which occurs when the stable fixed point of the resonance lock can no longer be maintained.For the prograde orbits on the left, the resonance lock is initially maintained as the mode frequency and tidal forcing fre-quency are both increasing.However, once the stellar rotation switches from spinning up to spinning down, the i-mode frequency begins to decrease and the resonant mode peak moves to decreasing forcing frequencies.The orbit thus loses the stable fixed point.
In the case of the retrograde resonance lock, the star is continuously spinning down by the time the lock begins at a few 100 Myr, so the lock breaks for a different reason.Figure 12 depicts the evolution of the relevant resonant mode peak for retrograde orbits shown in Figure 11, in terms of the value of Im k 2 2,2 (top) and the rate of change of the semimajor axis (bottom).The black dashed line shows the value of each quantity achieved by the planetary orbit, which in the regime of heightened dissipation can be thought of as the necessary value to maintain the fixed point.Since the peaks in both panels become shorter as the stellar rotation rate decreases, the maximum possible dissipation diminishes significantly over time.Furthermore, the value of Im k 2 2,2 required to maintain the fixed point actually increases with time, since more dissipation is required to sustain the necessary ȧ as the planet migrates away.Thus, the maximum dissipation available from the resonant peak is eventually unable to sustain the torque required for the fixed point.
For the P orb,0 = 0.5-2 day orbits of Earth-mass planets shown in Figure 7, the only resonance locks occur for i 2,0,1 , i 2,1,0 , and i 2,0,2 .On the other hand, for much closer orbits ≲ 0.5 day, we find that nearly every stellar oscillation mode is able to sustain resonance locks, given amenable initial conditions such that the tidal forcing frequency of the orbit can intersect with the mode frequency.

Viscosity and non-linearity
As discussed in Section 3.2.1, the off-resonant dissipation is orders of magnitude larger when calculating tidal dissipation using the viscosity independent of Rossby number ν NR (Equation 21).Furthermore, the resonant peaks associated with stellar oscillation modes are both wider and shorter.For the resulting orbital evolution, we find that the effect of using the larger viscosity ν NR is an increased prevalence of planets plunging towards the star due to non-resonant dissipation.In addition, there are fewer initial conditions P orb,0 for which resonance locking occurs.Compared to the fiducial models shown in Figure 7, we find that the prograde orbits decay on a shorter timescale, and the P orb,0 ≳ 1 day orbits do not catch into resonance.Since the resonant mode peaks are shorter for this higher viscosity, P orb,0 =1 day orbits are now too far away for the resonant peak in the dissipation to be able to sustain the requisite torque for a resonance lock.
If the convective viscosity is much larger (or if the interaction of convective and tidal flows produces enhanced dissipation that is not well-described with a convective viscosity; Terquem 2021), then equilibrium tidal dissipation would likely dominate the evolution rather than locks with inertial modes.However, if convective viscosity is smaller via the Goldreich & Keeley (1977) prescription as demonstrated by recent work (Duguid et al. 2020a,b), then equilibrium tidal dissipation will be less important.Tidal migration in that case is likely to be dominated by resonance locks, even for fairly massive planets.We leave the consideration of different prescriptions for the viscosity in the regime of rapid tides to future work.
Our study assumes that modes are linear, which we find holds true even during the resonance locks shown in Figures 7-8.Though resonance locking drives modes to larger amplitudes, typically the amplitude (Equation 11) remains small.In the case of Earth-mass planets, the mass ratio is small enough to ensure small mode amplitudes throughout resonance locks.Even for Jupitermass planets, in many cases the tidal overlap integral Q α ℓ,m is small because the resonance lock occurs while the star is rotating slowly; during other resonance locks, either the frequency separation ω α − ω m or the mode damping rate γ α is relatively large.Overall, we estimate mode amplitudes ξ on the order of ∼ 10 −5 -10 −4 , and wavenumbers k ∼ n 2 1 + n 2 2 + m 2 ≲ 4, so that the nonlinearity measure kξ < 1.Whether weakly non-linear dissipation (e.g., Weinberg et al. 2012) can increase tidal dissipation or prevent resonance locking is unclear for fully convective stars, and should be studied in future work.

Comparison with Observed Systems
The occurrence rate of Earth-like planets around Mdwarfs in the Kepler sample is less than a few percent for short orbital periods P orb,0 ≲ 1 day (Mulders et al. 2015;Hsu et al. 2020), where our models predict migration due to resonance locking.The story appears to be similar for TESS, as recent analyses find only a few planets with Earth-like radii at ≲ 1 day orbits out of the dozens of planets orbiting M-dwarfs in their sample (Rodríguez Martínez et al. 2023).The vast majority of planets lie at a few-10 day orbits, where negligible migration occurs in our fully convective stellar models.
Our results indicate that a dearth of Earth-mass planets at P orb ≲ 1 day around fully-convective M-dwarfs is expected due to migration induced by resonance locks.For Earth-mass planets at P orb,0 ≲ 1 day, planets orbiting slower than the stellar spin tend to migrate out to similar ∼1.5 day orbits, whereas planetary orbits rotating faster than the stellar spin tend to decay towards the star.For Jupiter-mass planets, we expect few planets out to ∼3 days as a result of migration due to resonance locking.
By initializing our orbital evolutionary calculations at different times, we are able to explore how our predicted outcomes vary for the timescales associated with different theories of the origin of short-period planets.For insitu planet formation or disk migration theories, changes to the orbital separation are dominated by the protoplanetary disk until the disk dissipates at a few-10 Myr (Dawson & Johnson 2018).High-eccentricity tidal migration may operate on a large range of timescales, as it is highly sensitive to the initial conditions of the underlying planet-planet interactions.For instance, the secular interactions and the ensuing tidal migration to short orbital periods may span anywhere from ∼Myr to tens of Gyr (Dawson & Johnson 2018).
For the fastest-spinning P rot,0 = 1 day model, the Earth-mass planets will experience negligible migration until the resonance lock at ∼1-2 Gyr causes planets at P orb,0 < 1 day to migrate outward.This is true during the first Gyr of the stellar evolution, no matter what time the orbital integration begins.As a result, we predict the same behavior regardless of when the planet reached its short orbital period, as long as it occurs before ∼ 1 Gyr.
In our slower-rotating stellar models, an assortment of behaviors may ensue depending on the formation timescale of short-period planets.For instance, in the P rot,0 = 5 day set of integrations, formation theories placing planets at P orb,0 ≲ 1.5 day at stellar ages ≲ 50 Myr are susceptible to resonance locking, either immediately beginning inward migration or experiencing delayed outward migration at a few 100 Myr-1 Gyr.However, if the Earth-mass planet is driven to short orbital periods at ≳ 50 Myr, the orbit may not be affected.For example, evolving a P orb,0 = 0.5 day planet from 500 Myr in the P rot,0 = 5 day model is too late for a resonance lock.In the P rot,0 = 10 day model, the P orb,0 = 0.5 day evolution initialized at 100 Myr is unable to lock into resonance, though the same evolution initialized at 10 Myr migrates inward due to resonance locking.Similar variations in behavior for different initial conditions can be observed for the P orb,0 = 1-1.5 day curves in Figures 7.
In general, the dependence of evolutionary outcomes on the different timescales of short-period planet formation can be attributed to the variation in initial tidal forcing frequency ω m as the stellar rotation rate changes over time.This aspect of the stellar evolution may prevent the planet from achieving the same resonance lock and enhanced migration that occurs for another initial condition.Figure 8 presents an analogous range in evolutionary outcomes that depend on the potential formation timescales for short-period Jupiter-mass planets, if they are found around low-mass M-dwarfs.Nevertheless, the caveat discussed in Section 4.1 for outwardlymigrating Jupiter-mass planets remains relevant here.
The majority of the migration illustrated in this work occurs before a few Gyr, so observing evidence of migration in the orbital architectures of exoplanets would be facilitated by the detection of young M-dwarf systems.In our predictions, planets with P orb ≲ 1 day that move outward to larger P orb often begin experiencing migration at ∼ 1 Gyr.For older M-dwarf systems found at a few Gyr, these planets will have already reached larger separations and vacated the region of short-period orbits.If detected in future observations, e.g. from TESS or the Roman Space Telescope, shortperiod planets around M-dwarfs younger than ≲ 1 Gyr likely have not yet migrated outward.As a result, we predict short-period planets to be more common around young M-dwarf systems than older M-dwarf systems.
Looking at the NASA Exoplanet Archive, we note that there are not many young, low-mass stars hosting short-period planets due to the challenges in observing planets around active young stars.Nevertheless, some young, rapidly rotating stars hosting planets at short periods such as TOI 540b (Ment et al. 2021) and K2-25b (Thao et al. 2020;Stefansson et al. 2020) are interesting candidates to undergo tidal evolution in the context of our models.In the future, we will be able to better test our model's predictions given more observations of planets around young stars.Nevertheless, to make a definitive prediction for evidence of tidal migration in the period distribution of observed exoplanet systems would require more extensive modeling than done in this work.The period distribution is shaped by combination of the physics of planet formation and the effect of tides, so it will be necessary to model both of these aspects in tandem (e.g., Lee & Chiang 2017) order to make more robust predictions in the future.
In exoplanet systems with more massive host stars, the star is not fully convective.Thus, the dynamical tide likely comprises gravito-inertial waves instead of the pure inertial modes examined in this work.For young, rapidly-rotating solar type stars with planetary companions at ≲ 5 day, for instance the host stars TOI 942 (Wirth et al. 2021) and WASP 25 (Bonomo et al. 2017), we do expect similar orbital dynamics due to gravito-inertial modes relative to what we discuss in this work.However, the denser spectrum of gravito-inertial modes at small frequencies will complicate the picture by involving contributions to the tidal dissipation from a larger number of normal modes.We anticipate that the rich spectrum of mode peaks will warrant detailed exploration of resonance locking effects in these systems.

Comparison to frequency-averaged dissipation
Models of similar star-planet systems in the past have relied upon a prescription for the frequency-averaged dissipation due to the dynamical tide.For instance, Gallet et al. (2017) modeled dynamical tides in stars between 0.3-1.4M ⊙ based on a prescription for frequencyaveraged dissipation from Ogilvie (2013): where ϵ = Ω s /Ω dyn and α = R c /R is the fractional size of the radiative core.In their models, no tidal dissipation occurred during periods of the evolution where the star was fully convective because of their prescription's dependence on α.As a result, for our models where α = 0, equation 29 underestimates the amount of tidal dissipation.
We can also compare the how dissipation scales with ϵ. Figure 13 shows the value of ⟨Im k 2 2,2 (ω) ⟩ ω in our models, with the fiducial models using ν R on the left, and the models with ν NR on the right.We see the same scaling of ϵ 2 in our models up to values of ϵ ≳ 0.2, where scatter begins to accumulate around the relation.This is not surprising, as the scaling was derived for slowly rotating bodies where the f-modes are well separated from the i-modes in frequency, assuming a constant dynamic viscosity (Ogilvie 2013).The P rot,0 = 1 day models form a group of outliers as they encounter avoided crossings between f-modes and i-modes for ϵ ≳ 0.6, such that the frequency range of [−2Ω s , 2Ω s ] includes the large fmode dissipation.This strongly amplifies the value of the frequency-averaged dissipation.We also note that the models with ν NR follow the trend more closely than models using ν R because the variation of viscosity with Rossby number within the star introduces more scatter.
Compared to Figure 4 of Gallet et al. (2017), our frequency-averaged dissipation values, which represent α = 0 in the context of their work, span a similar orderof-magnitude range of ⟨Im k 2 2,2 (ω) ⟩ ω ≈ 10 −8 -10 −4 for ϵ ≲ 0.1.In their higher-mass models, stars achieve this level of dissipation with core sizes α ≈ 0.2-0.8throughout the evolution.Their most similar model in terms  2013).The left panel shows the result for dissipation using rotationally modified viscosity νR (Equation 26), and the right shows the result with the unmodified viscosity νNR (Equation 21).For comparison, the scaling in black matches the prescription in Ogilvie (2013) of stellar mass is an 0.3 M ⊙ model which only hosts a radiative core for a few tens of Myr and correspondingly only exhibits dissipation during that period.In contrast, for our low-mass M-dwarfs that are fully convective, we predict frequency-averaged dissipation at the same level as their solar-mass models throughout the stellar lifetime.Furthermore, the goal of this work was to instead use a frequency-dependent tidal dissipation in our calculations, which reveal a variety of scenarios involving resonance locking that are not captured by an average over frequency.

Total dissipation vs. individual mode damping
Multiple approaches to calculating dissipation from the dynamical tide exist in the literature, primarily by direct calculation or by modal decomposition.The direct calculation refers to numerically solving the linearized fluid dynamics equations to find the perturbed response of a fluid under the influence of a tidal potential with forcing frequency ω m (see, e.g., Dewberry 2023; Ogilvie 2013).To recover the spectrum of Im k 2 2,2 across forcing frequency, the calculation must be repeated at each value of ω m .Other works use an expansion over oscillation modes, which constitutes the theoretical framework of this work (e.g., Press & Teukolsky 1977;Kumar et al. 1995;Lai 1997;Fuller & Lai 2012;Fuller 2017).This popular approach allows for a less numerically expensive analysis that readily captures resonances with normal modes of the star.
However, Townsend & Sun (2023) find a discrepancy between the two methods in the context of modeling secular tidal torques in binary star systems.They attribute the issue to a subtlety about expanding over normal modes: it is insufficient to consider each mode as 0.3 0.2 0.1 0.0 0.1 0.2 0.3 m ( dyn ) 10 0 10 1 10 2 10 3 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 10 3 10 2 10 1 10 0 10 1 Equation 17 Direct calculation Equation 30Figure 14.Comparison of different methods for computing the Love number, shown for a Prot,0 = 5 day, Ωs/Ω dyn = 0.150 model.The red solid curve shows the result of using Equations 14-17 to infer Im k 2 2,2 , which includes dissipative mode coupling terms.The blue dots show the results from a direct forced calculation, which are nearly indistinguishable from the red line.In contrast, the dotted blue line shows the result from a sum over normal modes without dissipative coupling terms (Equation 30), which deviates significantly from the other two methods away from resonance.
contributing independent coefficients characterized by individual mode damping rates γ α .For instance, works using modal decomposition typically find the Love number from the following sum (e.g., Dewberry 2023): This method employs Equation 16 to compute the mode amplitudes.
Though we also use Equation 16to calculate the mode amplitudes, we obtain coefficients in the damping integral (Equation 14) that couple different oscillation modes together (more details are discussed in Dewberry & Wu 2023).The discrepancy in the dissipation away from resonances as noted by Townsend & Sun (2023) motivates us to take the approach described in Section 2 of inverting Equation 17to find the Love number instead of using Equation 30.
In agreement with the inconsistency put forth by Townsend & Sun (2023), we demonstrate in Figure 14 that the spectrum of Im k 2 2,2 from Equation 30 diverges from the results of Equations 14-17 significantly away from resonance.In both calculations, we sum over the subset of f-modes and i-modes utilized in this work.However, Equation 30 (blue dotted line) can yield unphysical results.For instance, it does not cross zero at ω m = 0 with this set of modes, and the magnitude of the off-resonance dissipation is lower than what we find from our approach using the same viscosity.Figure 14 also verifies that our approach (red line) agrees very well with the results of a direct calculation (brown dots).Thus, the mode decomposition method can yield accurate results, but only if mode interaction cross terms due to dissipation (equation 14 in this work) are included when summing over modes.

CONCLUSION
We have modeled the frequency-dependent tidal response throughout the pre-main sequence and main sequence evolution of fully convective stars using realistic stellar models.Using a non-perturbative spectral method that fully resolves the effects of the Coriolis force to calculate fundamental and inertial mode oscillations of the star, we then compute the dissipation in terms of contributions from these normal modes.Our formalism accounts for the importance of viscous coupling between different modes in producing a reliable dissipation spectrum across tidal forcing frequency.
We performed our analysis for a fully-convective 0.2 M ⊙ M-dwarf at three different initial rotation periods, one of which (P rot,0 = 1 day) rotates very quickly relative to observed cluster stars, and two of which (P rot,0 = 5 and 10 day) span a more common range of stellar rotation periods.As the star traverses the premain sequence and main sequence, we prescribe its spin evolution using the saturated magnetic braking model of Matt et al. (2015).The tidal response within this fully convective star is expected to be damped by a turbulent convective viscosity, which can be estimated from mixing-length theory as independent of Rossby number Ro (Equation 21) or modified to account for the reduced efficiency of turbulent convective viscosity at high rotation rates (Equation 26).Including the effects of rotation on viscosity decreases the magnitude of the viscosity in the stellar interior, leading to smaller equilibrium tidal dissipation but larger dissipation at resonances with stellar oscillation modes.
Using the formalism described in Section 2, we compute the spectrum of dynamical tidal dissipation across tidal forcing frequency ω m throughout the stellar evolution.Peaks in the dissipation spectrum form at frequencies which are resonant with the star's inertial modes (imodes) and fundamental modes (f-modes).For rapidlyrotating models Ω s /Ω dyn ≳ 0.6, avoided crossings between f-modes and i-modes greatly enhance the dissipation at tidal forcing frequencies |ω m | ≤ 2Ω s .As the star evolves, the mode frequencies change, allowing them to sweep through resonances with orbiting planets.
We integrate the semimajor axis evolution of planets causing tidal dissipation in the host star.We find that both inward and outward orbital migration can occur for Earth-mass planets at P orb,0 ≲ 1.5 day and Jupitermass planets at P orb,0 ≲ 2.5 day.In general, significant migration is induced by resonance locking, whereby the planetary orbit's tidal forcing frequency resonantly excites a stellar oscillation mode and continues to closely follow this mode frequency over time.This leads to enhanced dissipation during the resonance lock by several orders of magnitude.Due to orbital migration during the first few Gyr of the star's lifetime, we predict that the region of P orb ≲ 1 day is likely to be cleared of small planets, and giant planets are made scarce out to P orb ≲ 3 day around fully-convective M-dwarfs.
In the future, we hope to incorporate the feedback of the tidal torque on the stellar spin in our calculations, which we find to be significant for Jupiter-mass planets.This will contextualize our expectations for resonance locking of Jupiter-mass planets, as the spin-up or spin-down of the host star will affect the mode frequencies as a function of time.Another promising next step will be to apply the formalism established in this work to higher-mass, partially-convective stars.These stars will host gravito-inertial modes, which we anticipate will lead to a denser spectrum of dynamical tidal dissipation that can also drive resonance locking in planetary systems.These solar-type stars are representative of the majority of known exoplanet hosts and in particular the host star in the only planetary system whose orbit is known to decay, WASP-12b (Maciejewski et al. 2013(Maciejewski et al. , 2018;;Patra et al. 2017;Yee et al. 2020).
In addition, observations of exoplanet hosts in this mass range indicate that cooler stars with convective envelopes, which are often observed with spin that is aligned with their planetary orbits, may experience greater tidal dissipation than hotter stars with radiative envelopes, whose spin-orbit alignments trend larger (Albrecht et al. 2012(Albrecht et al. , 2021;;Winn et al. 2010;Winn & Fabrycky 2015;Spalding & Winn 2022).While we only considered aligned spin and orbits, our methods can easily be extended to planets with misaligned orbits to compute the tidal damping of obliquities.Studying the dynamical tide in low and high-mass stars will elucidate the puzzle of obliquity damping in these systems.

Figure 1 .
Figure1.Evolution of the stellar radius (top panel), stellar spin frequency in units of the dynamical frequency Ω dyn = (GM/R 3 ) 1/3 (middle panel), and stellar rotation period (bottom panel) with stellar age for the 0.2 M⊙ stellar model used in this work.Different initial rotation periods Prot,0 of 1, 5, and 10 day are shown with the linestyles indicated in the legend; the radius evolution is the same for all three initial conditions.Blue stars indicate the stellar profiles used to calculate normal modes of oscillation and construct the grid of dissipation for each Prot,0.

Figure 2 .
Figure 2. The kinematic viscosity (first panel), scale height (second panel), convective velocity (third panel), and Rossby number (fourth panel) versus radial coordinate (scaled by the total stellar radius) for a set of profiles during the evolution of a Prot,0 = 5 day model.The top panel shows the turbulent convective viscosity with (νR, Equation26) and without (νNR, Equation21) a rotational correction.As Equation26is a function of Rossby number, νR varies with rotation and approaches the value of νNR towards the surface of the star where the Rossby number is large.On the main sequence, νNR, H, and vcon remain nearly constant so all lines lie on top of each other under the purple line.

Figure 4 .
Figure 4. Imaginary parts of k 22,2 as a function of tidal forcing frequency ωm (in units of the dynamical frequency Ω dyn ).The y-axis is scaled by ωm, and panels are labeled by the value of Prot,0 used to evolve the stellar model.In each panel, the red lines show the dissipation on the pre-MS for the stellar rotation rates Ωs (in units of Ω dyn ) listed in the legend, whereas the blue dotted lines show the same for a model on the MS with a very similar rotation rate.The x-axis limits in each panel are limited to[−2Ωs, 2Ωs]  to show the inertial mode range.Each resonant peak is labeled by the mode responsible for the resonance.The dissipation is larger on the pre-main sequence due to larger turbulent convective viscosity.

Figure 5 .
Figure 5. Same as Figure 4, but zoomed out to show the fmodes for the Prot,0 = 5 day models.The f-mode and largest i-mode resonances are labeled by the mode responsible for each resonance.

Figure 6 .
Figure 6.Similar to the center panel of Figure 4 (Prot,0 = 5 day), Im[k 22,2 ] scaled by ωm as a function of tidal forcing frequency ωm is shown, but now to compare the effect of viscosity.The red line shows the dissipation using the correction to viscosity due to rotation (Equation26).The blue dotted lines show the same but using viscosity without the rotation correction (Equation21).Including the effects of rotation lowers the viscosity, leading to sharper resonant peaks, but smaller off-resonant dissipation by 1-2 orders of magnitude.

Figure 8 .
Figure 8. Same as Figure 8, but for a Mp = Mj Jupiter mass companion.

PFigure 9 .
Figure9.The rate of change of stellar spin from tides due to Equation 19, Ωs,tide , compared to that due to stellar evolution (magnetic braking, contraction, or expansion), Ωs,evol .This is shown for an Earth-mass planet in blue and and a Jupiter-mass planet in red for the models listed in the legend.Note that the nearly vertical red and blue lines at 250 Myr and 800 Myr respectively are due to plunge-in of the planet (see corresponding times for the same models in top right panel of Figure8and top middle panel of Figure7).

Figure 10 .
Figure10.The value of the term (2 + c)M ′ a 2 /3Is in Equation 28 during resonance locking shown for the models listed in the legend.Orbits evolving from 10 Myr that migrate inward and from 100 Myr that migrate outward, as well as both Earth-mass (blue) and Jupiter-mass (red) planets, are shown.The resonance lock duration is shown as dashed sections in each curve.The dotted gray line denotes where (2 + c)M ′ a 2 /3Is = 2. Above this value, resonance locks will break due to tidal torque on the stellar spin, which was not included in these models.

Figure 11 .
Figure11.Examples of resonance locking for Earth-mass exoplanets from the Prot,0 = 10 day initial condition, for initial orbital periods 0.5 (left) and 1 day (right).The top panels show the time evolution of the imaginary part of the Love number k 2 2,2 (positive for prograde orbits ωm > 0, negative for retrograde orbits ωm < 0).The bottom panels show tidal frequency ωm and the rotating-frame mode frequency ωα for the i-mode in the legend.The left column shows an inwardly migrating resonance lock a prograde mode ωα > 0, and the right exemplifies an outwardly migrating resonance lock with a retrograde mode ωα < 0.

Figure 12 .
Figure 12.The imaginary part of the Love number (top panel) and the rate of change of the semimajor axis (bottom panel) versus tidal forcing frequency ωm.For this example of a retrograde resonance lock, Im[k 2 2,2 ] < 0 so −Im[k 2 2,2 ] is plotted.Black dashed lines show the evolution of each quantity for the planetary orbit.The peaked lines show the shape of the dissipation curve for a subset of times in the orbital evolution, shaded by the stellar age as shown in the colorbar.The dissipation is sharply peaked due to the resonant mode i2,0,1, and the peak center moves to the right as the mode frequency ωα evolves in time.The peak height decreases with time as the star spins down and the mode dissipation weakens.The resonance lock breaks when the dashed black line exceeds the height of the resonance peak.