Dynamical tide modified Roche limit in eccentric, asynchronous binaries

TL;DR

This paper develops a nonlinear hydrodynamical framework to refine the Roche limit in eccentric, asynchronous binaries, revealing that eccentricity and fluid dynamics significantly alter the threshold for mass transfer and instability.

Contribution

It introduces a comprehensive theoretical model for the Roche limit considering nonlinear hydrodynamics, eccentric orbits, and asynchronous rotation, extending beyond previous static and circular orbit assumptions.

Findings

Eccentricity increases the likelihood of fluid instability at pericenter.

The threshold separation for mass loss can be 30% higher in eccentric orbits compared to static predictions.

Chaotic fluid evolution can be triggered by pericenter passages even at moderate eccentricities.

Abstract

The Roche limit, or the threshold separation within which a celestial object (the donor) M cannot remain in a stable configuration due to a companion's tidal field, has been well established when M is in hydrostatic equilibrium and has synchronous rotation in a circular orbit. However, limited analyses exist considering corrections to the Roche limit due to hydrodynamical effects. We fill in the gap by providing a general theoretical framework involving nonlinear hydrodynamics. We consider both exact nonlinear equations derived from an affine model describing incompressible ellipsoids and series-expanded ones that can be calculated for realistic stars and planets. Our formulation addresses the Roche problem in generic orbits and synchronization levels of M, and fully accounts for the history-dependent hydrodynamical effects. We show that as the orbital eccentricity increases, fluid instability is more likely to develop at the pericenter due to the increased dynamical tide that accumulates over multiple orbits. When M moves in a highly eccentric orbit (with eccentricity around 0.9) and the damping of the fluid is small, the threshold pericenter separation at which mass loss from M can occur can be at least 30% higher than the value predicted for a circular orbit with hydrostatic equilibrium. If only a single passage is considered, however, the threshold separation is 20% smaller than the static limit. The nonlinear interaction at each pericenter passage can also trigger a chaotic fluid evolution inside M even with moderate eccentricities, complementing previous studies of chaotic tides caused by random propagation phases. Our work has broad implications for interacting binaries in eccentric orbits, including migrating gaseous exoplanets, repeated partial tidal disruption events, and more.

Equations

$$\ddot{\xi }={\omega }_f^2(\underset{\text{linear oscillator}}{\underbrace{-\xi }}+\underset{\text{linear tide}}{\underbrace{Iϵ}}+\underset{\text{combined nonlinear corrections}}{\underbrace{\kappa {\xi }^2}}),$$

$$\kappa {\xi }^2-\xi +Iϵ=0,$$

$$1-4\kappa Iϵ\ge 0\text{ or }ϵ\le \frac{1}{4\kappa I}\text{ for stability}.$$

$$U_{\mathrm{eff}}(\xi )=\underset{\text{internal energy + self-gravity}}{\underbrace{\frac{1}{2}{\omega }_f^2{\xi }^2-\frac{1}{3}\kappa {\omega }_f^2{\xi }^3}}\underset{\text{external interaction}}{\underbrace{-Iϵ\xi .}}$$

$${\omega }_f^2\kappa {\xi }^2-({\omega }_f^2-{\varpi }^2)\xi +{\omega }_f^{2}Iϵ=0,$$

$$({\omega }_f^2-{\varpi }^2{)}^2-4{\omega }_f^4\kappa Iϵ\ge 0\text{ or }ϵ\le \frac{({\omega }_f^2-{\varpi }^2{)}^2}{4{\omega }_f^4\kappa I}\text{ for stability}.$$

$$ϵ\le \frac{1-4({\varpi }^2/{\omega }_f^2)\kappa {\xi }_{\mathrm{dyn}}}{4\kappa I}\text{ for stability.}$$

$$\mathbf{\xi }(t,\mathbf{x})=\sum _m{q}_m(t){\mathbf{\xi }}_m(\mathbf{x}),$$

$${\mathbf{\xi }}_2=\sqrt{\frac{2\pi }{3{\omega }_f^2}}\nabla (x_r^2{Y}_{22})=\frac{5}{4}\left[(x+iy){\mathbf{e}}_x+(ix-y){\mathbf{e}}_{\mathbf{y}}\right],$$

$${\mathbf{\xi }}_0=\sqrt{\frac{2\pi }{3{\omega }_f^2}}\nabla (x_r^2{Y}_{20})=-\frac{5}{2\sqrt{6}}[x{\mathbf{e}}_x+y{\mathbf{e}}_{\mathbf{y}}-2z{\mathbf{e}}_z],$$

$${\mathbf{\xi }}_{-2}=\sqrt{\frac{2\pi }{3{\omega }_f^2}}\nabla (x_r^2{Y}_{2,-2})=\frac{5}{4}\left[(x-iy){\mathbf{e}}_x+(-ix-y){\mathbf{e}}_{\mathbf{y}}\right].$$

$${\ddot{q}}_m=-{\omega }_f^2{q}_m+{\omega }_f^2{K}_m,$$

$$K_2$$

$$K_0$$

$$r_i(t,\mathbf{x})=q_{ij}(t)x_j,$$

$$v_2^{(\mathrm{st})}=\frac{3}{4}ϵ>0,\text{ and }{v}_0^{(\mathrm{st})}=-\frac{\sqrt{6}}{4}ϵ<0,$$

$$q_m=c_{m,+}+c_{m,-}=c_{m,+}+c_{-m,+}^{\ast },$$

$${\dot{c}}_{m,+}=-i{\omega }_{f,+}{c}_{m,+}+\frac{i{\omega }_f{K}_m}{2},$$

$${\dot{c}}_{m,-}=i{\omega }_{f,-}{c}_{m,-}-\frac{i{\omega }_f{K}_m}{2}.$$

$${\omega }_{f,\pm }={\omega }_f\pm \frac{⟨{\mathbf{\xi }}_m,-i\mathbf{\omega }\times {\mathbf{\xi }}_m⟩}{⟨{\mathbf{\xi }}_m,{\mathbf{\xi }}_m⟩}={\omega }_f\mp \frac{m}{l}\omega .$$

$${\ddot{q}}_m=-{\omega }_f^2{q}_m+m{\omega }_f\omega (c_{m,+}-c_{m,-})+{\omega }_f^2{K}_m.$$

$$m{\omega }_f\omega (c_{m,+}-c_{m,-})\simeq \frac{(m\omega {)}^2}{l}{K}_m(\simeq \frac{3}{2}{\omega }^2ϵe^{-2i{\Phi }_C}\equiv C_2^{(\mathrm{st})},\text{for }m=2)$$

$${\ddot{q}}_m=-{\omega }_f^2{q}_m+{\omega }_f^2{K}_m+C_m^{(\mathrm{st})},$$

$$K_2$$

$$K_0$$

$${\ddot{q}}_m{e}^{im{\varphi }_q}=-{\omega }_f^2{v}_m+{\omega }_f^2{\bar{K}}_m+{\bar{C}}_m^{(\mathrm{st})},$$

$$a_1$$

$$a_2$$

$$a_3$$

$$v_0=\frac{\sqrt{6}[5a_{31}{v}_2+2(a_{31}-1)]}{5(2+a_{31})},$$

$$X_{22}{v}_2^2+X_{21}{v}_2+X_{20}=X_{02}{v}_2^2+X_{01}{v}_2+X_{00}=0.$$

$$\det [\mathbf{S}(X_{2n},X_{0n})]$$

$$\text{with }\mathbf{S}(X_{2n},X_{0n})$$

$$\sum _{n=0}^3{W}_{mn}{v}_2^n=0,$$

Full text

\uatRoche limit1404, \uatHydrodynamics1963, \uatExoplanet dynamics490, \uatExoplanet tides497, \uatTidal disruption1696, \uatX-ray transient sources1852

††software: NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2020), Matplotlib (Hunter, 2007), Numba (Lam et al., 2015), Mathematica (Wolfram Research, Inc., 2025).

I Introduction

The Roche model has long been the foundational framework for studying mass loss and transfer in interacting binary systems. Of particular importance is the Roche limit, a critical separation within which a donor object (with mass $M$ and radius $R$) can no longer hold itself but will lose mass due to the tidal force of the companion (treated as a point particle with mass $M^{\prime}$ here). Many applications of the Roche model rely on the fit by Eggleton (1983), which assumes that the donor is in hydrostatic equilibrium–a valid approximation for evolved circular binaries but one that may break down in eccentric systems, particularly when eccentricity is high and the donor has not reached synchronization.

Highly eccentric binaries can form in various environments. One such example lies in the realm of exoplanets. The formation of hot Jupiters (HJs) and hot Neptunes (HNs) at separations as close as $0.05\,{\rm AU}$ to their host stars is one of the longest outstanding puzzles in exoplanetary science (Mayor and Queloz, 1995; Dawson and Johnson, 2018). High-eccentricity migration of HJ and HN is a promising explanation (Fabrycky and Tremaine, 2007): A giant planet initially forms with a semimajor axis of $\mathcal{O}({\rm AU})$ and is dynamically excited onto a highly eccentric orbit. It then migrates inward due to tidal dissipation while the orbit circularizes. This picture naturally reproduces some key characteristics of HJs (Dawson and Johnson, 2018), including their highly inclined, sometimes retrograde orbits. While qualitatively successful, this channel still has large uncertainties on a quantitative level. For example, the pericenter separation required for inward migration can be dangerously close to tidal disruption of the planet, and a small variation in the disruption threshold can have a significant impact on the formation efficiency of HJs and HNs (Yu and Dai, 2024).

Another example is found in the nucleus of a galaxy that hosts a massive black hole (MBH). In such regions, a stellar-mass binary may be drawn close to the MBH and disrupted by its tidal forces through the Hills mechanism, which ejects one star while leaving the other bound to the MBH on an eccentric orbit (Hills, 1988). If the bound star approaches the MBH closely enough, it may itself be disrupted, resulting in a tidal disruption event (TDE; Rees 1988; Komossa 2015). At slightly greater separations, the star may survive a full disruption by the MBH and repeat mass transfer over multiple orbits, causing repeated partial TDEs (rpTDEs; Hayasaki et al. 2013; Coughlin and Nixon 2019; Cufari et al. 2022; Liu et al. 2024; Bandopadhyay et al. 2024, 2025; Yao and Quataert 2025). Recently, a new class of recurring X-ray bursts, known as quasi-periodic eruptions (QPEs), has been observed and is associated with low-mass MBHs (Miniutti et al., 2019). The mechanism driving QPEs remains uncertain, yet similar to the production of TDEs (Maguire et al., 2020), several proposed models involve a star or white dwarf (WD) in an eccentric orbit around an MBH, undergoing periodic mass transfer episodes that power the QPEs (King, 2020, 2022; Wang et al., 2022; Linial and Sari, 2023). The recent work of Lau and Yu (2025) demonstrates that both can be the result of a tidally-induced mass loss at different pericenter distances. For orbits with moderately high eccentricities ($0.90<e<0.99$) where the orbital period is short compared to the tidal dissipation timescale, the tide inside the WD builds up diffusively (i.e., as a random-walk process) over multiple orbits. The tidal energy eventually approaches the stellar binding energy and is released via a mass ejection as the extended tidal bulge breaks. Neither the WD nor the orbit will be disrupted by a single break, allowing the process to repeat several times. This may lead to rpTDEs with time intervals of a few years, consistent with the observations (e.g., Campana et al. 2015; Wevers et al. 2023; Miniutti et al. 2023; Somalwar et al. 2025). Meanwhile, the model also shows the possibility of QPEs being sourced by tidal stripping at closer pericenter separations, where the tidal field is strong enough such that a single passage is sufficient to induce mass transfer, though this mechanism is challenged by certain observed features in some QPEs, like the alternating long-short recurrence time and brightness of bursts (Miniutti et al., 2019; Giustini et al., 2020; Arcodia et al., 2022; Zhou et al., 2024). The proposed mechanism by Lau and Yu (2025) accounts for the dynamical nature of the tide over multiple orbits, improving previous studies of King (2022); Wang et al. (2022) based on the static Roche model that does not apply in eccentric orbits. Nonetheless, the threshold distance for the tidally induced mass transfer, either via multiple orbit buildup or a single passage, depends on a phenomenological parametrization that is not derived from fundamental hydrodynamics.

Despite the importance of mass loss in eccentric binaries, illustrated by the examples above, theoretical advancements beyond the classic Roche model are still limited. Some pioneering work by Sepinsky et al. (2007) extended the Roche model to non-circular orbits with non-synchronized donors, identifying significant corrections to the classic fit of Eggleton (1983); however, their work was still based on a quasi-static approximation. In reality, the donor also exhibits a dynamical response, corresponding to oscillations of the eigenmodes (mainly the fundamental modes, or f-modes) excited at each pericenter passage (Press and Teukolsky, 1977). This component is known as the dynamical tide. The analysis of Diener et al. (1995), based on an affine model developed by Carter and Luminet (1983), accounted for the dynamical tide but focused on the hydrodynamical stability in a single passage (though they considered the effects due to varying the ellipsoid’s orientation). Unlike the equilibrium component, which depends solely on the instantaneous orbital configuration, the dynamical tide can accumulate over multiple orbits, allowing its amplitude to grow to a scale comparable to the unperturbed radius (Mardling, 1995; Ivanov and Papaloizou, 2004; Wu, 2018; Vick and Lai, 2018; Yu et al., 2021). Such an accumulation is expected to play a critical role in the dynamics of mass transfer.

The numerical study by Guillochon et al. (2011) examining the threshold separation at which a proto-HJ is disrupted by its host star strongly motivates the need to include DT in modeling mass transfer. According to the fit of Eggleton (1983), Roche lobe overflow should only occur when the planet’s pericenter separation $D_{p}$ reaches $D_{p}{\lesssim}2.0R_{t}$ with $R_{t}{=}(M^{\prime}/M)^{1/3}R$ the tidal radius. However, Guillochon et al. (2011) found that with high orbital eccentricity ($e{=}0.9$), a planet is disrupted at a much greater pericenter separation of $D_{p}\simeq 2.7R_{t}$, $35\%$ greater than the separation from Eggleton’s fit. The discrepancy is further amplified when the results of Sepinsky et al. (2007) are used, which predict a Roche lobe $\sim 10\%$ greater than Eggleton’s fit (mainly due to the lack of centrifugal force), requiring an even smaller separation for the overflow to occur. Similarly, the critical single-passage separation of $D_{p}\simeq 1.7R_{t}$ from Diener et al. (1995) also cannot explain the numerical results of Guillochon et al. (2011).

In this work, we will show that the dynamical tide is the needed ingredient to reconcile the discrepancy. For this, we will consider the nonlinear hydrodynamical equations following both the affine model of Diener et al. (1995) (see also Carter and Luminet 1983, 1985) and those derived by Van Hoolst (1994) following an order-by-order expansion (see also Weinberg et al. 2012; Yu et al. 2023). Our model is fully dynamical and makes no (quasi-)static assumptions. Consequently, it applies to arbitrary eccentricities of the orbit and asynchronicities of the donor (though we still restrict our analysis to aligned spins). On the other hand, when the perturbed fluid is in hydrostatic equilibrium, our model also reproduces the classic results obtained by Chandrasekhar (1963) (see also Chandrasekhar 1987 for a comprehensive summary) on homogeneous ellipsoids.

The rest of the paper is organized as follows. We will first describe a toy model in §I.1 that captures the essence of our model and illustrates the significance of the dynamical tides. The nonlinear hydrodynamical model is introduced in §II, after which we will demonstrate its consistency with the analysis of Chandrasekhar (1963) in §II.1 under the static limit. A theoretical investigation on instantaneous fluid stability at the pericenter follows in §II.2. Then, in §II.3, a numerical experiment is conducted to address the build-up of the instability throughout multiple orbits. How the model applies to different astrophysical settings is described in §III. Lastly, we conclude in §IV.

Throughout the paper, we will use natural units $G=M=R=1$.

I.1 A toy model

The appearance of the Roche limit can be schematically shown with a toy model:

[TABLE]

where $\xi\in\mathbb{R}$ denotes the fluid displacement due to a large-scale perturbation (i.e., the f-mode of the donor whose eigenfrequency is denoted by $\omega_{f}$). The first term on the right hand side (RHS) describes the linear restoration of the fluid to the equilibrium, the $I\epsilon$ term the linear tidal drive, with $I>0$ a coupling constant, $\epsilon=M^{\prime}/D^{3}$ the instantaneous tidal field exerted by a companion (the accretor) with mass $M^{\prime}$, and $D$ the instantaneous binary separation. Lastly, the $\kappa\xi^{2}$ term represents a combination of nonlinear effects in both the internal fluid response and the external tidal drive. Theoretical calculation predicts $\kappa>0$ when $I>0$ (see Appendix A).

For the linear problem (ignoring the $\kappa\xi^{2}$ term), a static solution ($\ddot{\xi}=0$) can always be arranged by setting $\xi=I\epsilon$. With the nonlinear $\kappa\xi^{2}$ term, the RHS of Eq. (1) must take the opposite (same) sign as $\xi$ for the system to be oscillatory (grow exponentially). We immediately see that with $\kappa\xi^{2}\geq 0$, the nonlinear term can lead to instability when $\xi>0$, and an increasing $\xi$ will make the system more vulnerable to instability. This is consistent with the standard Roche lobe overflow picture where mass loss occurs mainly at the L1 and L2 Lagrangian points, corresponding to locations of maximum stretch of the donor. If we assume the hydrostatic limit and write $\xi\sim e^{\sigma t}$, at the instability threshold ($\sigma=0$), we need

[TABLE]

which is formally a quadratic function of $\xi$, and a real root exists only if its discriminant is non-negative

[TABLE]

We thus obtain that $1/(4\kappa I)$ is the maximum perturbation allowed. This is the essence of the static Roche limit found in Chandrasekhar (1963). See later in §II.1. Alternatively, Eq. (1) can also be viewed as a particle with displacement $\xi$ in an effective potential $U_{\rm eff}(\xi)$, with $\ddot{\xi}=-dU_{\rm eff}/d\xi$ and

[TABLE]

As the donor is sufficiently deformed, the total energy of the internal fluid can decrease (see Appendix A), creating a local maximum in the effective potential. Going beyond the local maximum will then lead to a runaway in $\xi$ (Wu, 1998).111In numerical experiments to be presented in §II.3, we confirm the presence of the instability by solving the fully nonlinear equations derived in Appendix A. While previous studies (e.g., Wu 1998) suggested that the four-wave interaction (the $\propto\xi^{4}$ term) in the effective potential could prevent the instability, the stabilization effect is eventually taken over by even higher-order corrections. When considering the f-modes, the full stability threshold is similar to that obtained at the three-wave order, as presented in the toy model here.

To see the effect of the dynamical effects, we first consider $\xi\sim\cos(\varpi t)e^{\sigma t}$ with $\varpi^{2}\lesssim\omega_{f}^{2}$. Plugging it back to Eq. (1) and consider the threshold case, $\sigma=0$, we have

[TABLE]

which is again a quadratic function of $\xi$. It has a real solution only if

[TABLE]

Compared to the static case of Eq. (3), including the dynamical effect makes instability more likely to occur as the maximum $\epsilon$ is smaller by a factor of $(\omega_{f}^{2}-\varpi^{2})^{2}/\omega_{f}^{4}$. This model will be relevant for the discussions in §II.2 about finite-frequency corrections to the equilibrium tide, relevant mainly to asynchronous donors in mildly eccentric binaries.

Another way to incorporate the dynamical effect is to consider $\xi=\xi_{\rm dyn}\cos(\varpi t)+\eta e^{\sigma t}$ with $\xi_{\rm dyn}\in\mathbb{R}$ treated as a known quantity. Eq. (1) is now treated as a quadratic function of $\eta$, and a solution with $\sigma=0$ can be arranged, by again examining the discriminant, if

[TABLE]

The $\xi_{\rm dyn}$ term reduces (increases) the maximum $\epsilon$ a stable configuration can allow if $\xi_{\rm dyn}>0$ ($\xi_{\rm dyn}<0$). If a positive $\xi_{\rm dyn}$ has triggered some exponential growth in $\eta$, suddenly changing its sign will only make $\eta$ oscillatory but not damp away the growth. Therefore, for a sufficiently large $\xi_{\rm dyn}$ with a randomly fluctuating sign, the net effect will again make the instability more likely to occur. This model is relevant to the dynamical tide in a donor moving in a highly eccentric orbit in §II.3.

The goal of our study is to quantify the enhanced hydrodynamical instability caused by the dynamical tide through analytical and numerical investigations.

II Nonlinear hydrodynamics

We focus mainly on a homogeneous, incompressible fluid object with a polytropic index $n=0$ and an adiabatic index $\Gamma\to\infty$ for the theoretical discussion, though our numerical investigations later in §II.3 will be extended to touch upon $n=1$ polytropes. This choice is for simplicity and to assist a direct comparison with the classical results obtained by Chandrasekhar (1963, 1987). Denote the Lagrangian displacement of a mass element $dM$ as $\boldsymbol{\xi}=\boldsymbol{r}-\boldsymbol{x}$ with $\boldsymbol{r}$ and $\boldsymbol{x}$ the perturbed and original position vectors, respectively. We further decompose $\boldsymbol{\xi}$ into eigenmodes in the configuration space, as

[TABLE]

For the Roche limit calculation, it is sufficient to consider only the quadrupolar ($l=2$) fundamental modes (f-modes). We drop the $l\geq 3$ perturbations, restrict our analysis to aligned spins when rotation is incorporated, and ignore the inertial modes (Ho and Lai, 1999; Lai and Wu, 2006; Braviner and Ogilvie, 2015; Xu and Lai, 2017). Therefore, we label the modes only with their azimuthal quantum numbers $m$. The eigenfunctions $\boldsymbol{\xi}_{m}$ are computed from a non-rotating background as (e.g., Poisson and Will 2014)

[TABLE]

and we have used $\omega_{f}^{2}=4/5$ as the (non-rotating) eigenfrequency under the linear theory. The modes are normalized so that $\omega_{f}^{2}\langle\boldsymbol{\xi}_{m^{\prime}},\boldsymbol{\xi}_{m}\rangle=\omega_{f}^{2}\int\boldsymbol{\xi}_{m^{\prime}}^{\ast}\cdot\boldsymbol{\xi}_{m}\rho d^{3}x=\delta_{mm^{\prime}}$. A mode with $|q_{m}|=1$, therefore, has a potential energy of unity (also approximately the magnitude of the binding energy of $M$). Note that the mode amplitudes $q_{m}$ are complex and can be further decomposed as $q_{m}=v_{m}e^{im\phi_{q}}$ with $v_{m}\in{\mathbb{R}}$ and $m\phi_{q}$ describing the phase. The $m=\pm 1$ modes are not excited by the tide and are dropped from the discussion.

Using the modal expansion, $\boldsymbol{\xi}(t,\boldsymbol{x})$ is fully determined once $q_{m}(t)$ is obtained. First, we ignore rotation. The equations of motion of $q_{m}$ are given by

[TABLE]

where

[TABLE]

with $\epsilon=M^{\prime}/D^{3}$. We use $D$ and $\Phi$ to denote orbital separation and phase. The derivation follows a variational principle and is presented in detail in Appendix A. Note that in a homogeneous ellipsoid, the perturbation satisfies the affine model of Carter and Luminet (1983, 1985); Diener et al. (1995), where the perturbed position $\boldsymbol{r}$ is related to the unperturbed position $\boldsymbol{x}$ via

[TABLE]

where the mapping between $q_{ij}(t)$ and $q_{m}(t)$ is given in Eq. (A2). The affine model allows the exact nonlinear form of $K_{m}$ to be obtained (available in Appendix A, Eqs. A25-A30, and allows for arbitrary values of $n$ and $\Gamma$ with the caveat that Eq. 13 may not be accurate beyond $n=1/(\Gamma-1)=0$). Nonetheless, we expand the results to the three-wave order (that is, we count $q_{m}\sim\epsilon$ and include terms to $\mathcal{O}(\epsilon^{3})$ in the energy variations and $\mathcal{O}(\epsilon^{2})$ in the equation of motion). It is already sufficient to reproduce the Roche limit at this order. More importantly, such expanded equations can also be obtained following Van Hoolst (1994), which takes a general set of linear eigenmodes as input to compute their nonlinear couplings and applies to generic compressible fluids (realistic stars and planets). This is advantageous over the affine model, which only self-consistently describes idealized homogeneous ellipsoids. Therefore, our study bridges the two different formulations of nonlinear hydrodynamics. The exact agreement between the two formulations at the three-wave order for homogeneous ellipsoids is presented in Appendix B.

We pause here and briefly describe how the nonlinear term leads to instability. As the eigenfunctions $\boldsymbol{\xi}_{m}\propto e^{im\phi}$ (Eq. 9), we multiply both sides of Eq. (10) by $e^{im\phi}$ and set $\phi=\phi_{q}$. Effectively, this corresponds to evaluating the fluid displacement along the major axis of the perturbed ellipsoid. At the linear order and under the static limit $\ddot{q}_{m}=0$, we have

[TABLE]

with $v_{m}=(q_{m}e^{im\phi})_{\phi=\phi_{q}}$. Treating $(-v_{0})\sim v_{2}\sim\xi$, we see that Eq. (10) takes the same form as the toy model of Eq. (1), and the nonlinear terms in $K_{m}$ tend to destabilize the system when $\xi>0$. Thus, instability is most likely to occur when $\xi$ reaches its maximum along the major axis. In the equilibrium limit (without mode resonances; following the convention of Yu et al. 2024), where the tidal bulge follows the companion, we have $\phi_{q}=\Phi$, matching the intuition that mass transfer first occurs at the L1 Lagrangian point. At the quadrupolar order, the L2 point ($\phi=\Phi+\pi$) destabilizes at the same time.

Before solving for the instability threshold quantitatively, we now extend the equations of motion to include rotation. In this analysis, we focus on the case where the background spin is uniform and the spin axis is aligned with the orbital angular momentum (along $\boldsymbol{e}_{z}$). The equations of motion are now described in a frame that corotates with $M$ with three modifications.

The first is a simple Doppler effect, which corresponds to a phase shift $\Phi_{C}=\Phi-\omega t$. The mode amplitude changes accordingly, as $q_{mC}=q_{m}e^{im\omega t}$ while keeping $q_{m}e^{im\Phi}=q_{mC}e^{im\Phi_{C}}$. Quantities with a subscript $C$ are evaluated in the corotating frame. Since the analyses later will be in the corotating frame as the default, we will drop the subscript $C$ in the mode amplitude $q_{m}$ from this point onward. For the orbital phase, however, we will still make the distinction between the inertial frame $\Phi$ and the corotating frame $\Phi_{C}$.

The second correction comes from the centrifugal force whose detailed incorporation is presented in Appendix A. Since the force can be derived from a specific centrifugal potential $d\mathcal{R}/dM=-\omega^{2}(r^{2}-z^{2})/2$, its impact on the equations of motion can be obtained through a variational principle by considering the energy $\mathcal{R}$. When doing so, it is important to introduce an $l=0$ mode whose amplitude we denote as $q_{\rm tr}$ (a perturbation of the trace of $q_{ij}$ in Diener et al. 1995) and allow it to be coupled to the $q_{m}$ mode at the three-wave level to properly capture the correction at $\mathcal{O}(\epsilon^{2})$. The same result can also be obtained from the formulation of Van Hoolst (1994) by splitting $d\mathcal{R}/dM=(d\mathcal{R}/dM)_{l=0}+(d\mathcal{R}/dM)_{l=2}=-\omega^{2}r^{2}/3-\omega^{2}(x^{2}+y^{2}-2z^{2})/6$, and adding $(d\mathcal{R}/dM)_{l=0}$ to the background (monopolar) gravitational potential and $(d\mathcal{R}/dM)_{l=2}$ to the first-order, quadrupolar Eulerian potential perturbation (same as the tidal interaction potential). Note that our approach treats the centrifugal potential as a (nonlinear) drive of the modes defined by a non-rotating background, instead of redefining a new set of eigenmodes of a non-spherical background where the centrifugal perturbation is absorbed. In other words, the centrifugal potential is treated on an equal footing as the tidal potential.

Lastly, to incorporate the Coriolis force to the lowest order in $\omega$, we split the corotating frame mode amplitude $q_{m}$ as

[TABLE]

with

[TABLE]

This splitting is based on the phase-space expansion of Schenk et al. (2002), and the Coriolis force is absorbed as a frequency shift (Lai, 2021; Yu et al., 2024; Yu and Lau, 2025),

[TABLE]

Higher-order corrections from the Coriolis force are ignored (but see, e.g., Dewberry and Lai 2022). The corresponding second-order equation can be obtained by differentiating $\dot{c}_{m,\pm}$ one more time and eliminating the first-order derivatives with the equations of motion, leading to

[TABLE]

Compared to Eq. (10), there is an additional term $\propto(c_{m,+}-c_{m,-})$ due to the Coriolis effect. While this in principle represents a new degree of freedom (as prograde and retrograde modes are split due to rotation), under the static limit $c_{m,\pm}\simeq[1\pm(m/l)\omega/\omega_{f}](K_{m}/2)$, we have

[TABLE]

which is $\mathcal{O}(\omega^{2}\epsilon)\sim\mathcal{O}(\epsilon^{2})$. Since this is already as small as the smallest terms in Eqs. (10), we will approximate this term as $C_{2}^{\rm(st)}$ from Eq. (19) in our analytical discussion. When doing numerical experiments, we will solve directly the equations of $c_{m,\pm}$ (Eq. 16) to account for the additional degree of freedom.

Summarizing, the equations of motion we will solve in §II.1 and §II.2 are

[TABLE]

where

[TABLE]

and $C_{2}^{(\rm st)}$ is given by Eq. (19) while $C_{0}^{\rm(st)}=0$. We will also consider

[TABLE]

where $\bar{K}_{m}=K_{m}e^{im\phi_{q}}$ and $\bar{C}_{m}^{\rm(st)}=C_{m}^{\rm(st)}e^{im\phi_{q}}$. In the equilibrium limit (as considered in §II.1 and II.2) without free oscillations of the mode due to resonances, $\phi_{a}=\Phi_{C}$, and one just needs to replace $q_{m}$ with $v_{m}$ and drop the phases following $\epsilon$ in $K_{m}$ to obtain $\bar{K}_{m}$. Eq. (23) is real valued.

II.1 Consistency with previous analyses

First, we consider the case where $M$ moves in a circular orbit with synchronous rotation, so the left-hand sides of Eq. (20) or (23) vanish. Our goal is to establish the consistency between our formulation based on an order-by-order expansion with the classic analyses of Chandrasekhar (1963) under this static limit. We write the spin frequency $\omega$ as $\omega^{2}=p\epsilon$, with $p=1+M/M^{\prime}$ for synchronized rotation in a circular orbit. Note that our $p$ differs from the one in Chandrasekhar (1963) by 1, as $p=p^{\rm(Ch)}+1$. The superscript “(Ch)” will be used when we quote notations or results from Chandrasekhar (1963). In our analysis, we will focus on $M^{\prime}\gg M$ (as for planet migration or TDEs), so a synchronous donor has $p=1$.

Inspired by Chandrasekhar (1963), we note that the perturbed fluid can be described by an ellipsoid with $\bar{x}^{2}/a_{1}^{2}+\bar{y}^{2}/a_{2}^{2}+\bar{z}^{2}/a_{3}^{2}=1$ at the surface. The barred coordinates are aligned with the axes of the ellipsoid. In the limit with circular orbit and synchronous rotation, the barred frame coincides with the corotating frame. The axes are given by (Eqs. A5 and A31)

[TABLE]

For a given $a_{31}=a_{3}/a_{1}$, it defines a linear relation between $v_{0}$ and $v_{2}$,

[TABLE]

Using this relation to eliminate $v_{0}$, Eq. (23) now formally become two quadratic equations of $v_{2}$,

[TABLE]

where the coefficients $X_{mn}$ are functions of $(\epsilon,a_{31})$. The first subscript indicates the azimuthal number $m$ of the corresponding equation (of $\ddot{q}_{2}$ or $\ddot{q}_{0}$), while the second indicates the power in $v_{2}$. The expressions of the coefficients $X_{mn}$ are straightforward to obtain but lengthy, so we report their values in Appendix C.

If there exists a $v_{2}$ that simultaneously solves the above two polynomials, then the resultant of the two polynomials must vanish (Weisstein, 2025). That is,

[TABLE]

where $\boldsymbol{S}(X_{2n},X_{0n})$ is also known as the Sylvester matrix, and the resultant is its determinant. For given $a_{31}$, Eq. (27) can then be used to determine the value of $\epsilon$, thereby determining an equilibrium sequence.

Some examples of equilibrium sequences (plotted as $\epsilon$ versus $\cos^{-1}a_{31}$ or $e_{13}=\sqrt{1-a_{31}^{2}}$) are shown in Fig. 1. We use the color gray to represent a synchronized donor with $p=1$ (as $M\ll M^{\prime}$), corresponding to the Roche sequence, while yellow for the non-rotating case with $p=0$ (the Jeans sequence; note $\ddot{q}_{m}$ are ignored in this calculation). As the ellipsoid becomes more deformed (towards the right of the figure), the $\epsilon$ that leads to an equilibrium does not increase monotonically. Instead, it reaches a maximum, $\epsilon_{\rm max}$, at around $\cos^{-1}a_{31}\simeq 60^{\circ}-65^{\circ}$. This means for $\epsilon>\epsilon_{\rm max}$, no equilibrium configuration of the ellipsoid can be arranged, indicating hydrodynamical instability. We can also translate $\epsilon>\epsilon_{\rm max}$ to $(D/R_{t})<(1/\epsilon_{\rm max})^{1/3}$. In other words, when the binary separation $D$ measured in the tidal radius $R_{t}$ is less than $(1/\epsilon_{\rm max})^{1/3}$, the donor $M$ cannot remain in equilibrium, and hydrodynamic instability occurs, defining the Roche limit.

In our calculation, we find that it is sufficient to capture the classic results of Chandrasekhar (1963) with percent-level accuracy using only three-wave interactions (solid lines; Eqs. 21 and 22). For the synchronous ($p=1$) case, we find $\epsilon_{\rm max}\simeq 0.0747$ at $\cos^{-1}a_{31}\simeq 64.8^{\circ}$, meaning the threshold separation is at $D^{(\rm th)}=(1/\epsilon_{\rm max})^{1/3}R_{t}\simeq 2.38R_{t}$, matching the classic Roche limit (indicated by the horizontal dotted lines) of $2.44R_{t}$ to $2.5\%$. A similar level of accuracy is also obtained for the non-rotating $p=0$ case, with $\epsilon_{\rm max}\simeq 0.110$, $(1/\epsilon_{\rm max})^{1/3}\simeq 2.09$, achieved at $\cos^{-1}a_{31}\simeq 65.4^{\circ}$.

The accuracy can be improved further by including the next order, four-wave interactions. The equations of motion at this order are derived in Appendix A, Eqs. (A33a) and (A33b). Still using Eq. (25) to eliminate $v_{0}$ in terms of $v_{2}$ and $a_{31}$,222Strictly, at the four-wave order one should account for corrections to $a_{1,2,3}$ from the monopolar mode, $q_{\rm tr}$; Eqs. (A5) and (A31). This introduces nonlinear corrections to Eq. (25), which we ignore here.

we again obtain two polynomial equations of $v_{2}$,

[TABLE]

by setting $\ddot{q}_{m}=0$, with $m=2$ and [math]. The resultant (now functions of $W_{mn}$) needs to vanish as before for the two equations to hold simultaneously. We can thus find the equilibrium sequences at the four-wave order, shown in the dashed lines in Fig. 1. The result now matches those from Chandrasekhar (1963) almost exactly. For the $p=1$ case shown in gray, we have $\epsilon_{\rm max}\simeq 0.0687$ at $\cos^{-1}a_{31}=61.1^{\circ}$, corresponding to a Roche limit of $D^{(\rm th)}=(1/\epsilon_{\rm max})^{1/3}R_{t}=2.44R_{t}$.

While the main result of this Subsection is to reproduce the classic ones by Chandrasekhar (1963) on homogeneous ellipsoids, it is significant for two reasons. First, it validates our formulation based on which $\ddot{q}_{m}\neq 0$ corrections will be introduced. Second, it bridges the analyses of Chandrasekhar (1963) on idealized ellipsoids and Van Hoolst (1994) on general stars and planets (see also Weinberg et al. 2012; Weinberg 2016; Yu et al. 2023). In general, exact nonlinear expressions are challenging to find, but an order-by-order expansion as in Van Hoolst (1994) is straightforward to compute (though lengthy). Our result provides a framework for how one can start from Van Hoolst (1994) to determine the Roche limit of a realistic star/planet theoretically.

II.2 Finite-frequency corrections in asynchronous, low eccentricity systems

In this Subsection, we will extend the stability analysis to analytically consider finite-frequency corrections due to the equilibrium tide (following the convention of Yu et al. 2024), where the perturbation varies with the instantaneous tidal forcing. In other words, the tidal bulge still points towards the companion. The component that varies at a mode’s natural frequency and depends on the excitation history of the mode (what we call the dynamical tide) will be deferred to the numerical investigation in §II.3 and ignored in the analysis here. This means that the prediction derived here will be accurate for binaries with low eccentricity $e\lesssim 0.3$. For more eccentric binaries, the dynamical tides can exhibit complicated, often diffusive, evolutions that may trigger instabilities at separations greater than the prediction of this Subsection by tens of percent.

To proceed, we write $q_{m}=(1+z_{m})q_{m}^{(0)}$ where $q_{m}^{(0)}\simeq I_{m}\epsilon e^{-im\Phi_{C}}$ represents an approximate solution of the perturbation, with $I_{2}=3/4$ and $I_{0}=-\sqrt{6}/4$ (Eqs. 21 and 22). We treat $z_{m}$ as a low-frequency variable (whose time derivatives will be ignored) that we will use to determine the equilibrium structure of the donor. Note Eq. (23) still holds as long as we incorporate the time derivatives on the left-hand side due to $\ddot{q}_{m}\simeq\ddot{q}_{m}^{(0)}$. In essence, this model corresponds to the toy model described in Eq. (6).

The procedure for finding the Roche limit follows identically those described in §II.1. Choosing again a linear relation between $v_{0}$ and $v_{2}$ as in Eq. (25), Eq. (23) now becomes two quadratic equations of $z_{2}$:

[TABLE]

The coefficients, now depending on $(\epsilon,a_{31})$ as well as $q_{m}^{(0)}$ and their derivatives, are given in Appendix C. Requiring the resultant of the two polynomial equations in Eq. (30) to vanish allows one to find $\epsilon$ given $a_{31}$. For a given orbit and synchronization level, the maximum $\epsilon$ in the equilibrium sequence then determines the corresponding Roche limit. As the tidal deformation is the strongest at the pericenter, we will evaluate all the quantities entering $Z_{mn}$ there. Consequently, the analysis considers the instantaneous hydrodynamical stability at each pericenter passage.

To find $\ddot{q}_{m}^{(0)}$, we note

[TABLE]

at the pericenter (denoted with the subscript $p$) and $q=M/M^{\prime}(=0\text{ in this work by default})$. A simple approximation of $q_{m}^{(0)}\simeq I_{m}\epsilon e^{-im\Phi_{c}}$ with $\Phi_{c}=\Phi-\omega t$ leads to

[TABLE]

We have evaluated the result at the pericenter and assumed an unperturbed Keplerian orbit.

For practical usefulness, we show in Fig. 2 the threshold pericenter separation, $D_{p}^{(\rm th)}/R_{t}=(1/\epsilon_{\rm max})^{1/3}$, below which hydrodynamical instability can occur. Different colors represent different spin rates relative to the pseudo-synchronization values $\omega_{\rm ps}$ in Hut (1981),

[TABLE]

where $P_{\rm orb}$ is the orbital period. Also shown as a reference in the gray-dashed line is the threshold $D_{p}^{\rm(th)}$ under the static tide limit (obtained by setting $\ddot{q}_{m}=0$). At a fixed eccentricity, the nonrotating configuration has the largest $D_{p}^{(\rm th)}$ (most prone to instability) as it has the largest $|\ddot{q}_{m}^{(0)}|$ (Eq. 32) among the configurations we consider. As the rotational rate increases, $D_{p}^{(\rm th)}$ first decreases due to the reduction of the $|\ddot{q}_{m}^{(0)}|$ term, and then increases as the centrifugal potential dominates. The minimum $D^{\rm(th)}_{p}$ is achieved when $\omega/\omega_{\rm ps}\simeq 0.5$ for circular binaries. When the orbital eccentricity $e$ increases, $D^{\rm(th)}_{p}$ increases approximately linearly with $e$, especially for slowly rotating systems. For synchronous systems, there is a decrease in $D_{p}^{(\rm th)}$ when $e\lesssim 0.1$ due to the decrease in $\omega_{\rm ps}$, which reduces the centrifugal force.

We pause here to discuss further the meaning of Fig. 2. As the coefficients in Eq. (30) are evaluated at the pericenter, Fig. 2 should be interpreted as the instantaneous instability threshold at the pericenter. Therefore, a donor $M$ marginally unstable at the pericenter may become stable once it moves away from the pericenter and $\epsilon$ (defined by the instantaneous separation) drops below the instability threshold. As a result, the growth of hydrodynamical instability during a single pericenter passage may not disrupt $M$ sufficiently to cause mass loss. Nonetheless, as long as the fluid dissipation timescale is long compared to the orbital period (which is very well satisfied for the large-scale f-mode when $1-e>0.01$ under the linear theory assuming damping due to turbulent convection; Goldreich and Nicholson 1977; Wu 2005), the growth of instability at every pericenter passage will accumulate over multiple orbits, eventually disrupting $M$ and triggering mass loss (see also numerical examples to be presented in the next Subsection). For systems with low eccentricities $e\lesssim 0.3$, Fig. 2 provides a decent approximation of the onset of mass loss as a runaway in mode amplitude occurs within a few ($<10$) orbits. For systems with higher eccentricities, the amount of growth per passage is smaller, given the smaller fraction of time $M$ spends near the pericenter ($2\pi/(\dot{\Phi}_{p}P_{\rm orb})\propto\sqrt{(1-e)^{3}/(1+e)}$). This requires more orbits for the perturbation to build up until a runaway eventually occurs. However, the dynamical tide component can also be excited in those systems, significantly lowering the threshold $\epsilon$ for a stable configuration (see Eq. 6 and note the dynamical tide plays the role of $\xi_{\rm dyn}$ in the toy model). In the weak damping case, mass loss from $M$ can occur at $D_{p}$ significantly greater than the prediction of Fig. 2 for highly eccentric systems.

II.3 Dynamical tide induced instability and diffusive tide

For binaries with significant eccentricities, the truly dynamical component of a mode can be excited at a pericenter passage. Such a component oscillates at a mode’s natural frequency (hence phase incoherent with the orbit) and does not decrease in amplitude as the companion moves to the apocenter (assuming dissipation is small). In the next pericenter passage, the existing component of the dynamical tide can interfere with the new tidal kick it receives at the pericenter, hence making this component’s evolution history dependent (as opposed to the equilibrium tide that depends only on the instantaneous orbital configuration). This dynamical component formally modifies the stability threshold according to the toy model in Eq. (7), though its complicated evolution prohibits a simple analytical study as we have done in §II.1 and II.2. Nonetheless, the hydrodynamical equations are sufficiently generic for arbitrary rotations and orbital configurations. We will use a numerical approach to investigate the associated hydrodynamical stability.

In this Subsection, we will solve directly the phase-space expanded equations of $c_{m,\pm}$, Eqs. (16), so that the Coriolis effect can be accounted for naturally. For the driving terms ($K_{m}$), we will consider both the version expanded to the three-wave order, Eqs. (21) and (22), and the exact nonlinear version derived in Appendix A following the affine model of Diener et al. (1995) (see also Carter and Luminet 1983, 1985). The monopolar mode $q_{\rm tr}$ is treated fully dynamically, where a differential equation of $\ddot{q}_{\rm tr}$ is solved instead of reducing $q_{\rm tr}$ under the static limit with Eqs. (A31). Note that the affine model in Diener et al. (1995) allows for generic values of polytropic index $n$ and adiabatic index $\Gamma$. Thus, our numerical experiments will consider both a homogeneous ellipsoid with $(n,\ \Gamma)=(0,\ 20)$ as before333Numerically, $\Gamma\to\infty$ is hard to implement. We therefore use a finite $\Gamma=20\gg 1$ as an approximation. , and a polytrope with $(n,\ \Gamma)=(1,\ 2)$ that serves as a good approximation of a buoyancy-neutral Jovian planet. However, a caveat for the $n=1$ model is that the starting point of the affine model (Eq. 13) is not self-consistently satisfied as the Lagrangian displacement $\boldsymbol{\xi}$ is no longer linear in position. We fix the companion’s mass to be $M^{\prime}=10^{3}$ and treat it as a point particle. When computing the orbital motion, we ignore tidal back-reaction and simply model it as Keplerian, with

[TABLE]

The initial conditions are set at the apocenter, with the modes computed with their linear equilibrium tide solution. The set of differential equations is evolved using the solve_ivp function in SciPy (Virtanen et al., 2020) with an explicit Runge-Kutta method of order 8 (Hairer et al., 1993).

In the first experiment shown in Fig. 3, we validate the results of §II.2 by considering the numerical evolution of $|q_{2}|=|c_{2,+}+c_{-2,+}^{\ast}|$ for a non-rotating ($p=0$), homogeneous ($n=0$) star in a moderately eccentric orbit with $e=0.3$. We evolve the system using both the exact nonlinear equations (left panel) and the three-wave-expanded ones (right panel). Qualitatively, the two sets of equations yield similar results, with the exact version allowing the runaway to occur at a slightly further separation $D_{p}/R_{t}\simeq 2.72$. This is consistent with Fig. 1, which shows that the three-wave equations slightly underestimates $D_{p}^{(\rm th)}\propto(1/\epsilon_{\rm max})^{1/3}$.

The numerical result of Fig. 3 also matches the theoretical prediction of Fig. 2, which predicts $D_{p}^{(\rm th)}(e=0.3,p=0)\simeq 2.66R_{t}$ at the three-wave order. However, a subtle difference is worth emphasizing. Fig. 2 predicts the instantaneous stability at the pericenter. Fig. 3 shows that the instability may not directly lead to a runaway in the mode amplitude444In the following, we will use the words “instability” and “runaway” with different meanings. “Instability” is said to occur when a mode experiences exponential growth instantaneously (typically near the pericenter). If the pericenter passage timescale is much shorter than the growth timescale, the mode may become oscillatory, and the donor may remain bounded by its self-gravity once it moves away from the pericenter. In contrast, we use “runaway” when a mode’s amplitude grows significantly to exceed order unity within a single pericenter passage. A runaway in the mode amplitude will require the donor to significantly restructure itself to remove the energy in the oscillation, which likely results in mass loss. in a single passage. As the donor moves away from the pericenter and $\epsilon$ decreases, the mode stabilizes and becomes oscillatory. However, in the next pericenter passage, the growth resumes, eventually leading to a divergence in the mode amplitude, potentially triggering mass loss from the donor.

The accumulation of the dynamical tide becomes more significant in systems with higher eccentricities, with one such example shown in Fig. 4 for a donor moving in a highly eccentric ($e=0.9$) orbit. We still consider a homogeneous donor but vary its rotation rates (left for the non-rotating case and right for the pseudo-synchronized case). Both cases are evolved using the three-wave equations of motion. As the hydrodynamical instability grows only momentarily at each pericenter passage, an eventual runaway in the mode amplitude requires multiple orbits to occur. On the other hand, this accumulation allows runaway to occur at pericenter separations much greater than the prediction of Fig. 2. For the non-rotating case, the left panel of Fig. 4 shows that the runaway can kick in at a separation as large as $D_{p}/R_{t}\simeq 3.20$, whereas without the dynamical tide, the separation needs to be $D_{p}\lesssim 2.93R_{t}$ based on Fig. 2. Rotation suppresses the dynamical tide slightly (Lai, 1997), requiring a smaller $D_{p}/R_{t}\simeq 3.05$ for the runaway to occur within 20 orbits. We also note a larger variation in $|q_{2}|$ away from the pericenter. This is due to the beat of the $c_{2,+}$ (prograde) and $c_{-2,+}$ (retrograde) modes entering $q_{2}$. Note from Eq. (17) that the two modes have different natural frequencies when the background is rotating.

Note the mode evolution largely resembles the diffusive tide (also known as the chaotic tide) studied in, e.g., Mardling (1995); Ivanov and Papaloizou (2004, 2007); Wu (2018); Vick and Lai (2018); Yu et al. (2021, 2022). We will elaborate on this point further in §III.1.

As our theoretical study is largely motivated by the numerical simulations of Guillochon et al. (2011), we consider in Fig. 5 the evolution of a donor with $n=1$ and $\Gamma=2$ (a decent approximation for gaseous giant planets) moving in an eccentric orbit with $e=0.9$. In the affine model, the linear tidal overlap ($\propto I_{m}$) of an $n>0$ model is smaller than that of the $n=0$ model by a factor of $(5-n)\mathcal{M}_{\ast}$ with $\mathcal{M}_{\ast}=(\int x^{2}dM)/3$; see Appendix A. For $n=1$, $(5-n)\mathcal{M}_{\ast}\simeq 0.523$. As a result, the $n=1$ model requires a $D_{p}/R_{t}$ smaller than the $n=0$ model by a factor of $\sim 1/[(5-n)\mathcal{M}_{\ast}]^{1/3}\simeq 1.24$ for the runaway to occur. Other than this overall scaling difference, the evolution of the $n=1$ model largely resembles the $n=0$ model, and the predictions of the exact nonlinear equations (left) can be largely captured by the three-wave equations (right). Remarkably, our results show nice consistency with the simulations of Guillochon et al. (2011). Although the donor is not significantly disrupted in a single passage, the perturbation accumulates diffusively, eventually leading to a runaway in the mode amplitude in $\mathcal{O}(10)$ orbits when $D_{p}/R_{t}\lesssim 2.7R_{t}$.

In Fig. 6, we show a survey of the maximum change in the mode energy $\Delta E_{\rm mode}$ reached in 22 orbits versus the orbital eccentricity $e$ and pericenter separation $D_{p}/R_{t}$. The mode energy is defined as $\Delta E_{\rm mode}=\Delta T+\Delta U+\Delta\Omega$, corresponding to the sum of the kinetic energy $\Delta T$ (Eq. A7), internal energy $\Delta U$ (Eq. A9), and self-gravity $\Delta\Omega$ (Eq. A18). At the leading order and for $n=0$, $\Delta E_{\rm mode}\simeq(1/2)\sum_{m}(\dot{q}_{m}^{2}/\omega_{f}^{2}+|q_{m}|^{2})$ in our normalization. The exact nonlinear equations from the affine model are used in the hydrodynamical evolution and the energy calculation (see Appendix A; Carter and Luminet 1985; Diener et al. 1995). The donor is assumed to be non-rotating. In the left and right panels, we use respectively $(n,\Gamma)=(0,20)$ and $(n,\Gamma)=(1,2)$. In the color scale, we saturate at values where runaways in the mode amplitudes are observed. We empirically determine these saturation values by plotting $\Delta E_{\rm mode}$ as a function of $|q_{2}|$, which leads to two clearly separate sectors. For smaller values $|q_{2}|$, we observe a nearly quadratic relation $\Delta E\propto|q_{2}|^{2}$, corresponding to bounded oscillations. The other sector has $|q_{2}|$ piled at unity (a criterion we use to terminate the integration), and at this fixed $|q_{2}|$, $\Delta E_{\rm mode}$ is multi-valued. The maximum $\Delta E_{\rm mode}$ from the first branch is taken to set the maximum in the color bar.

We note an approximately linear increase in the threshold $D_{p}/R_{t}$ with increasing orbital eccentricity $e$. This is the most critical modification of the earlier prediction of Sepinsky et al. (2007), which mainly considered the equipotential and found an expanded Roche lobe in eccentric orbits (harder for mass loss to occur; opposite of our results). Our study, therefore, highlights the significance of accounting for the hydrodynamical response of the fluid. Compared to the prediction assuming instantaneous equilibrium tide in Fig. 2, the accumulation of the dynamical tide over multiple orbits makes the runaways occur at even greater values of $D_{p}/R_{t}$. We also note that at fixed $e$, the energy $\Delta E_{\rm mode}$ does not vary monotonically from red to blue as $D_{p}$ increases. Instead, a random alternation of red and blue is observed near the threshold, indicating the stochastic nature of the evolution (see §III.1 for more discussions). We further emphasize that Fig. 6 serves only as a conservative estimation of the Roche limit, as the integration lasts only 22 orbits. Figs. 4 and 5 suggest that more systems may run away given more orbits, as the evolution is diffusive. However, integrating the system for longer would require a more careful treatment of the numerical scheme beyond brute-force integration of the differential equations to improve both efficiency and accuracy. This is left for future investigations.

So far, the discussions have ignored damping. If there exists some efficient damping mechanism that could remove all the energy in the dynamical tide within one orbit, we would need to consider the hydrodynamical stability in a single passage. Fig. 7 shows the energy change acquired after a single passage (still assuming non-rotating $M$). When plotting $\Delta E_{\rm mode}$ versus $|q_{2}|$, we do not see a clear split into two sectors, and $|q_{2}|$ does not grow exponentially. Instead, at locations where bounded oscillations terminate in the multi-orbit case, we note $\Delta E_{\rm mode}$ becomes multi-valued at a given $|q_{2}|$. Such branching suggests instability to additional perturbations (and runaways occur in the next orbit), so we still saturate the color bar as before. We’ve also checked with three-wave equations of motion (which predict runaways in a single passage) and found consistent thresholds in $D_{p}/R_{t}$ as in Fig. 7 (which also approximately match Diener et al. 1995).

This time, we note a decrease in $\Delta E_{\rm mode}$ with eccentricity. This is because the tidal interaction is significant only for a fraction of $2\pi/(\dot{\Phi}_{p}P_{\rm orb})\propto\sqrt{(1-e)^{3}/(1+e)}$ of the orbit. Although the mode becomes instantaneously unstable near the pericenter as long as $D_{p}/R_{t}$ is below the prediction of Fig. 2, the time for the instability to grow can be too short compared to the instability growth timescale, thereby preventing a runaway that could lead to mass loss. The threshold $D_{p}/R_{t}$ can differ by order unity when comparing Fig. 6 and 7 due to the dynamical tide. Using the linear damping calculation of Goldreich and Nicholson (1977) with a quality factor of $\mathcal{O}(10^{13})$, it requires $1-e\lesssim 10^{-8}$ to suppress the dynamical tide within a single orbit. Even if the quality factor is reduced to $10^{6}$ (Wu, 2005), damping is significant only when $1-e<10^{-3}$. For systems less eccentric, Fig. 6 should be a more accurate representation of the Roche limit. As a caveat, the discussion considers only linear damping. Kumar and Goodman (1996) suggested that nonlinear damping can be highly efficient with a quality factor of $\mathcal{O}(100)$. However, their calculation focused on main sequence stars with masses $0.5-0.75\,M_{\odot}$ and was based on order of magnitude estimations. Whether nonlinear damping can be as strong as claimed in Kumar and Goodman (1996) for donors of different types remain to be examined in details by future studies.

III Applications

After describing the general theoretical results in §II, we now describe the applications of the model in various astrophysical settings.

III.1 Relation to the diffusive tide

It is interesting to note that the numerical mode evolution studied in §II.3 largely resembles the diffusive (or chaotic) tide studied in, e.g., Mardling (1995); Ivanov and Papaloizou (2004, 2007); Wu (2018); Vick and Lai (2018); Yu et al. (2021, 2022); Lau and Yu (2025). In fact, our work presents a new mechanism to trigger the diffusive tide. In all previous analyses, the diffusive behavior ultimately originates from a random phase accumulation away from the pericenter. Such a random phase could be produced by fluctuations either in the orbital periods due to tidal back-reactions (Mardling, 1995; Ivanov and Papaloizou, 2004; Wu, 2018; Vick and Lai, 2018) or gravitational-wave radiation (Ivanov and Papaloizou, 2007; Lau and Yu, 2025), or in the mode frequencies due to nonlinear anharmonicities (Yu et al., 2021, 2022). A limitation of this mechanism is that the diffusive evolution can hardly be maintained once the orbital eccentricity drops below $e\simeq 0.9$ and the orbital period is not long enough for the random phase to accumulate. For example, for $M^{\prime}=10^{3}$, $e=0.7$, and $D_{p}/R_{t}=2.645$, we have $P_{\rm orb}=164.5$ (deliberately chosen to be far away from linear resonance) and $E_{\rm orb}=-MM^{\prime}(1-e)/(2D_{p})=-5.7$ (all in natural units). After a single passage, the mode acquires an energy $\Delta E_{\rm mode}\simeq 1.7\times 10^{-3}$ for the $(n,\Gamma)=(1,2)$ model. The one-pass perturbation of the mode phase from the tidal back-reaction is $2\Delta\phi_{q}=\omega_{f}\Delta P_{\rm orb}\simeq(-3/2)\omega_{f}P_{\rm orb}(\Delta E_{\rm mode}/E_{\rm orb})\simeq 0.091\,{\rm rad}$, much smaller than the threshold value of $2\Delta\phi_{q}\gtrsim 1\,{\rm rad}$ derived in Vick and Lai (2018). Moreover, back-reaction is turned off in our simulations. The anharmonicity-induced phase perturbation differs from the back-reaction phase by a factor of $\mathcal{O}(|E_{\rm orb}|)$ and typically has a greater numerical value (Yu et al., 2021), yet $2\Delta\phi_{q}=P_{\rm orb}\Delta\omega_{f}\sim\omega_{f}P_{\rm orb}\Delta E_{\rm mode}\simeq 0.34\,{\rm rad}$ is still smaller than unity to trigger the diffusive growth. Nonetheless, our Fig. 6 shows that this system indeed grows diffusively and eventually runs away.

Crucially, all studies referenced above used the linear theory of Press and Teukolsky (1977); Lai (1997) to compute the tidal kick each mode receives at the pericenter. In this work, we incorporate nonlinear hydrodynamics at each pericenter passage, which may allow the mode to grow momentarily due to instability. If the mode’s phase at pericenter is periodic, then the growth will repeat after a few orbits and accumulate. This allows the mode to reach an amplitude large enough to trigger the diffusive tide. The growth then continues stochastically until the mode eventually runs away.

We emphasize that the instability at the pericenter enters at the three-wave order (as we have demonstrated in previous sections both theoretically and numerically), a lower-order effect compared to the nonlinear anharmonicity of Yu et al. (2021), which fundamentally enters at the four-wave order555While Yu et al. (2021) also considered three-wave interaction between f- and p-modes, the p-modes are not directly sourced by the tide. Therefore, to produce the anharmonicity, it requires nonlinearly sourcing the p-modes first and then letting the p-mode back-react on the f-mode; see Eq. (36). Since the three-wave coupling is used twice, it effectively becomes a four-wave interaction. The direct four-wave interaction was incorrectly ignored in Yu et al. (2021). There is also a numerical error in Yu et al. (2021) that overestimates the coupling strength between f- and p-modes.. It is also distinct from the nonlinear kick of Arras et al. (2023). Importantly, the correction to the tidal kick in this work is nonlinear in the fluid displacement $\boldsymbol{\xi}$, which enters the modes’ equations of motion at $\mathcal{O}(|q_{m}|^{2})\sim\mathcal{O}(M^{\prime 2}/D_{p}^{6})$; see, e.g., Eq. (21). In Arras et al. (2023), however, the correction is nonlinear in $(\Delta D_{p}/D_{p})$, where $\Delta D_{p}$ is the change in the pericenter separation due to tidal back-reaction. In other words, the tidal kick in Arras et al. (2023) is evaluated along a tidally perturbed orbit instead of the background Keplerian one. While ignored in this study, one can show that $(\Delta D_{p}/D_{p})\propto(D_{p}/M^{\prime})|q_{m}|^{2}\sim(M^{\prime}/D_{p}^{5})$ (see eq. 13 of Arras et al. 2023; note that the second term is smaller than the first one by $\mathcal{O}[(R_{t}/D_{p})^{3/2}]$). When evaluating $\epsilon$ at $(D_{p}+\Delta D_{p})$ instead of $D_{p}$, the correction to the tidal forcing is $\Delta\epsilon\propto(M^{\prime}/D_{p}^{3})(\Delta D_{p}/D_{p})\sim(M^{\prime 2}/D_{p}^{8})$, which is smaller than the nonlinear-in-$\boldsymbol{\xi}$ correction by $(1/D_{p}^{2})$. Since the new correction considered here is of lower order than those studied previously, the diffusive (chaotic) tide it leads to may occur over an expanded parameter space (especially for low eccentricities) than predicted by linear tidal kicks at the pericenter.

In turn, for highly eccentric systems, the diffusive tide can make the threshold pericenter separation greater than the prediction of Fig. 6. Note that in Fig. 6, the modes are initiated by their leading-order equilibrium tide solution at the apocenter. For a fixed $D_{p}$, a large eccentricity essentially suppresses the initial amplitude by a factor of $[(1+e)/(1-e)]^{3}$. Over merely 22 orbits, the diffusive growth is limited as, on average, the mode energy increases linearly with the number of passages (Wu, 2018; Vick and Lai, 2018). One can nonetheless capture the long-term growth history of the modes by giving them different initial conditions. One such example is shown in Fig. 8. We consider in this example the non-rotating $n=1$ model with $e=0.9$ and $D_{p}/R_{t}=2.9$. Starting with the equilibrium amplitude, the mode amplitude remains oscillatory within 22 orbits (the bluest curve; also Fig. 6). However, if a significant initial amplitude is used (indicated by curves with different colors), the mode may run away at this large separation. As the diffusive growth is almost unbounded (until reaching near equipartition with the orbit), the threshold $D_{p}/R_{t}$ for the mode’s runaway can be as large as that for the diffusive tide to occur.

In fact, the runaways in our simulation correspond to the catastrophic wave breakings during the diffusive evolution postulated in Wu (2018). When a mode’s amplitude runs away, it cannot be supported by the background and has to break somehow. However, our model does not describe the consequence of such a break. One possibility, as suggested in MacLeod et al. (2022), would be shock-heating the surface layer with some potential mass ejection from the donor. It further suggests that energy dissipation in highly eccentric systems can be very efficient due to not a continuous dissipation as in the linear theory but discrete epochs of diffusive growth followed by wave breakings as pointed out by Wu (2018). Meanwhile, mass loss is likely nonconservative based on Guillochon et al. (2011). More detailed analysis (theoretical or numerical) of the saturation of the instability is left for future investigations.

III.2 High eccentricity migration of exoplanets

Related to specific astrophysical settings, the first application we consider is the formation of HJs via the high-eccentricity migration channel (see, e.g., Dawson and Johnson 2018 for a review). In particular, our model provides a theoretical explanation for why Guillochon et al. (2011) observed mass loss from a migrating Jovian planet at a pericenter separation as large as $D_{p}\simeq 2.7R_{t}$, $35\%$ greater than the Roche-lobe-overfilling separation of $D_{p}\simeq 2.0R_{t}$ using the fit by Eggleton (1983). Extrapolating the low-eccentricity Roche lobe calculation of Sepinsky et al. (2007) will increase the discrepancy even further. We highlight that our Fig. 5 shows remarkable similarities to the simulations of Guillochon et al. (2011); see, e.g., their fig. 3. Specifically, there should be no significant mass loss from $M$ (the planet) in a single passage. Instead, the perturbation accumulates over orbits, and significant mass loss from $M$ eventually occurs after $\mathcal{O}(10)$ orbits. This consistency validates our analysis.

More importantly, our theoretical model (in terms of simple ordinary differential equations; Eq. 20 or 16) allows the more expensive three-dimensional hydrodynamical simulations of Guillochon et al. (2011) to be faithfully extrapolated to different systems based on each system’s unique configuration. In particular, the commonly adopted disruption threshold of $2.7R_{t}$ (e.g., Anderson et al. 2016; Vick et al. 2019) is only true for the specific system studied in Guillochon et al. (2011) over $\mathcal{O}(10)$ passages. More orbits of integration may increase the threshold further due to the diffusive build-up of the mode (captured by a non-equilibrium initial condition shown in Fig. 8). Recent studies by Yu and Dai (2024); Weldon et al. (2025) adopted a (potentially more realistic) soft disruption prescription where the planetary mass loss per orbit is parameterized as an exponentially decaying function of $D_{p}/R_{t}$, yet the parameters are still calibrated to Guillochon et al. (2011) and held constant throughout the evolution.

If either the orbital eccentricity or the planet’s rotational rate varies, a different disruption threshold (or a different set of mass loss prescriptions) should apply; see Figs. 6 and 4. The dependence of $D_{p}^{\rm(th)}$ on these properties may indicate interesting and unanticipated correlations between the current HJ distribution and progenitor configurations.

Our Fig. 8 further suggests that the Roche limit for highly eccentric systems can be as large as the threshold separation to trigger the diffusive tide at about $3.8R_{t}$ (Wu, 2018) for a typical proto-HJ. Interestingly, this distance coincides with the location where a pile-up of HJs is observed (with a current period of about 3 days; Wu 2018; Dawson and Johnson 2018). On the other hand, HJ systems with separations shorter than this could imply a high spin rate of the proto-HJ (young system), as the thresholds for both diffusive tide (e.g., fig. 4 of Yu et al. 2021) and Roche limit (Fig. 4) are suppressed by rotation.

Our Fig. 6 may also have implications for formation channels other than the standard high-eccentricity migration. For example, a planet may have disk-migrated before experiencing a high-eccentricity migration. In such a hybrid scenario, the eccentricity excitation required in the second stage is only moderate, with $e<0.9$ (e.g., figs. 5 and 10 of Yu and Dai 2024). This two-stage migration is suggested by Guillochon et al. (2011) based on observations available at the time and supported by the recent investigation of Yu and Dai (2024) on the formation of the WASP-107 system (their section 4.3). Both studies found that the initial apocenter distance of the planet of the high-eccentricity migration stage is within the snow line. Our Fig. 6 suggests that the HJ/HN formation efficiency for such a system with $e<0.9$ may be higher than originally expected based on a fixed $D_{p}^{(\rm th)}=2.7R_{t}$ derived from $e=0.9$. An increase in the predicted HJ formation rate will improve the agreement between theory and observation (see, e.g., table 2 of Dawson and Johnson 2018).

As discussed above in §III.1, the nonlinear tidal interactions leading to the Roche limit can also help trigger and maintain the diffusive tide. For highly eccentric systems with $e>0.9$, an expanded parameter space for the diffusive tide is crucial for explaining the lack of super-eccentric proto-HJs Socrates et al. (2012); Dawson et al. (2015). The efficient energy dissipation due to the diffusive tide naturally explains the paucity, yet the original study of Wu (2018) applies only to those with current periods less than 3 days, a small fraction of the missing population (Dawson and Johnson, 2018). Expanding the threshold for the diffusive tide is therefore crucial for reducing the tension with observations. For those with $e<0.9$, the propagation phase of a mode itself is not sufficiently random, and the nonlinear instability at the pericenter is likely what maintains the diffusive tide and drives the orbital evolution. It is important to note that there is no consensus on the dissipation mechanism in moderately eccentric systems after the quench of the diffusive tide. The predicted dissipation rate of f-modes from first-principle calculations (Goldreich and Nicholson, 1977) is orders of magnitude below what is needed for a planet to circularize without the diffusive tide (but see, e.g., Wu 2005; Ogilvie 2014).

We plan to follow up on the points discussed above with a population synthesis study utilizing a disruption threshold described in this work instead of fixing it to $D_{p}^{(\rm th)}=2.7R_{t}$, to explore new features and correlations in the current distribution.

III.3 rpTDEs and QPEs

Repeating X-ray bursts at galactic centers provide valuable information about the properties of the central MBH and its surrounding environments. A rpTDE corresponds to a star partially disrupted by the tidal force of the MBH, with the remnant continuing to orbit and getting disrupted again in the future. Compared to a full TDE, rpTDEs have a larger $D_{p}/R_{t}$, and a different mass fallback rate Coughlin and Nixon (2019). Through a hydrodynamic simulation, Ryu et al. (2020) finds that a main sequence star undergoes a partially disrupting instead of a full TDE for $2\gtrsim D_{p}/R_{t}\gtrsim 0.5$ for a single passage in the Newtonian limit. Rosswog et al. (2009b) performs a simulation of the tidal disruption of a WD orbiting a $10^{3}~M_{\odot}$ intermediate mass black hole, where mass loss occurs at $1.1R_{t}$ in a single passage. Unlike full TDEs that are typically assumed to happen on nearly parabolic orbits, rpTDEs are likely to have more moderate but still significant eccentricities with $1-e\sim 0.1-0.01$ (Hayasaki et al., 2013; Guillochon and Ramirez-Ruiz, 2013; Cufari et al., 2022; Liu et al., 2024; Bandopadhyay et al., 2024, 2025). The long orbital period limits previous numerical simulations to binaries with $D_{p}/R_{t}\lesssim 2$ in order for mass loss to occur in a single pericenter passage. However, the orbital period is not long enough to damp out the dynamical tide (unless the quality factor is significantly lower than $10^{6}$ in which case the study of Yao and Quataert 2025 applies; see the discussions at the end of §II.3), allowing it to build up over multiple orbits and potentially change the single-passage predication significantly. Similar to the case with Jovian planets, our study provides a new, theoretical model to describe the disruption, accounting for the accumulation history of the dynamical tide over multiple orbits. In our picture, the dynamical tide builds up diffusively until it runs away and ejects mass to power a burst. The stochastic growth can be a promising mechanism to produce rpTDEs with irregular recurrence times (Lau and Yu, 2025, 2026), which has been observed in, e.g., J0456-20 (Liu et al., 2024) and HLX-1 (Godet et al., 2014; Webb et al., 2023) that previous models requiring mass loss at every pericenter passage struggle to reproduce. For example, the simulations in Liu et al. (2024) require a fractional mass loss of almost 90% to reproduce the $\sim 30\%$ fluctuations in the recurrence times of J0456-20 from $\sim 300\,{\rm days}$ to $\sim 200\,{\rm days}$. Furthermore, this amount of mass loss is also in tension with the inferred eruption luminosity.

Moreover, the threshold pericenter separation $D_{p}^{(\rm th)}$ for the disruption to marginally occur may increase by $\gtrsim 60\%$ with the buildup compared to the single-passage predictions (e.g., Guillochon and Ramirez-Ruiz 2013; Ryu et al. 2020; comparing Figs. 6 and 7 at $e=0.9$). Since $D_{p}^{\rm(th)}$ sets the loss cone radius, an increase in $D_{p}^{(\rm th)}$ increases the rpTDE event rate (Bortolas et al., 2023). If we assume an isotropic distribution of a star’s velocity vector at a fixed distance for simplicity, the increase in loss cone event rates is proportional to the increase in its radius (Merritt, 2013).

Our results also complement some previous models on QPEs. In King (2022); Wang et al. (2022), it is proposed that a highly eccentric WD-MBH binary can simultaneously produce TDEs and QPEs. The predicted tidal stripping based on Eggleton’s formula is strictly periodic with $D_{p}<2R_{t}$. The QPEs produced via this mechanism would therefore be highly regular, which contradicts the observed fluctuations in the recurrence time of bursts. The stochastic nature of the mode provides a natural explanation for these fluctuations. The diffusive evolution also means that the tidal disruption does not repeat every orbit, as the mean mode energy builds up diffusively over time. Its corresponding impact on the orbital evolution and mass transfer has been studied in Lau and Yu (2025). A significant shortcoming within this work is that although the parametrized mass loss formalism is flexible enough to model the case with dynamical tide, the range of the parameters is prescribed to match the observations and previous hydrodynamic simulation results in Guillochon et al. (2011). In other words, it still lacks a direct connection to the first-principles calculation considering the hydrodynamics under an external tidal field. This paper provides a theoretical description of the development from mode oscillations into a runaway process that can potentially lead to mass loss.

For a deep encounter ($D_{p}/R_{t}<0.5$) of a WD with the MBH, a thermonuclear runaway can be triggered by the tidal compression Luminet and Pichon (1989); Rosswog et al. (2009a); MacLeod et al. (2016); Tanikawa (2018), depending on its composition. Based on the model by Carter and Luminet (1983), the star undergoes a phase in which it becomes flattened in the orbital plane and the interior is strongly compressed. Luminet and Pichon (1989) shows that it can trigger nucleosynthesis processes similar to that at the core of a main-sequence star. The energy release can be comparable to that of Type 1a supernovae Rosswog et al. (2009a). This scenario is favored for lower mass MBHs, of $\sim 10^{3}M_{\odot}$, where the WD orbit can penetrate deeper without getting directly swallowed. Whether a WD ignites or loses mass and powers rpTDE depends on how $\Delta E_{\rm mode}$ from a runaway is released. If $\Delta E_{\rm mode}$ is mostly used to eject mass, then rpTDE is the outcome. If, however, it is mostly converted into heat in the core, then the WD may be ignited in a thermonuclear explosion. Determining the outcome is beyond the scope of the current model, but an interesting question to follow up on in the future.

IV Conclusion and discussion

In this work, we studied the nonlinear hydrodynamical stability of a donor $M$ in a compact orbit around a point-particle companion $M^{\prime}$. Our study considers both the exact nonlinear equations derived from the affine model (Diener et al., 1995), which self-consistently apply to homogeneous, incompressible ellipsoids, and the expanded equations following Van Hoolst (1994), which can be readily extended to realistic stars and planets. Our study extends the classic Roche limit derived by Chandrasekhar (1963) of a donor in hydrostatic equilibrium with circular orbit and synchronous rotation to the fully dynamical regime, allowing for arbitrary orbital eccentricities and donor rotation rates (though restricted to aligned spins). Schematically, the finite-frequency corrections to the equilibrium tide (accurate for low eccentricities; §II.2) can be captured by the toy model in Eq. (6). The dynamical tide excited in more eccentric systems (§II.3) corresponds to the toy model of Eq. (7). More specifically, the Roche limit for moderately eccentric systems ($e\lesssim 0.3$) is described by our Fig. 2. While a small asynchronicity tends to decrease $D_{p}^{(\rm th)}/R_{t}$ (i.e., making the donor more tightly bound) by reducing the centrifugal force, eccentricity can dominate and make $D_{p}^{(\rm th)}/R_{t}$ greater (donor more subject to mass loss) than the prediction of Chandrasekhar (1963) by the finite-frequency corrections. For a more eccentric binary with $e\sim 0.9$, the fluid perturbation can grow diffusively due to nonlinear interactions at the pericenter. This eventually leads to a runaway in mode amplitude at a pericenter separation $\gtrsim 30\%$ greater than the Roche limit of circular orbits with synchronous rotations (Figs. 6 vs. 1). Our theoretical result explains the numerical simulation of Guillochon et al. (2011) on the high-eccentricity migration of gaseous planets, and allows the result to be physically extrapolated across possible parameter space and to other similar systems, e.g., rpTDEs. This study also provides a new mechanism to trigger and maintain the diffusive (chaotic) tide that is of lower order in $(1/D_{p})$ than all previously considered mechanisms (§III.1).

We mainly focused on a homogeneous, incompressible ellipsoid to assist a direct comparison with Chandrasekhar (1963). This simple model is also sufficient to illustrate the crucial role of dynamical effects in the Roche problem and its corrections to the commonly used fit of Eggleton (1983). We established that predictions based on the exact nonlinear equations of an affine ellipsoid can be largely reproduced by the expanded equations to the three-wave order, or $\mathcal{O}(\xi^{3})$ in $\Delta E_{\rm mode}$. Including the next order improves the accuracy, but does not modify the qualitative features (Fig. 1). While the affine model of Diener et al. (1995) allows for different values of $n$ and $\Gamma$, the starting point, Eq. (13), is not self-consistently satisfied by a general fluid. Nonetheless, the expanded equations of Van Hoolst (1994) will apply generically; see also Weinberg et al. (2012); Weinberg (2016). This allows the theoretical study of the Roche limit to be extended to realistic astrophysical objects. It also allows other modes, such as g-modes (present if $\Gamma>1+1/n$), or those with $l\geq 3$ (important for binaries with $M^{\prime}\sim 1$) to be included in the problem. Based on the affine model, an increasing $n$ reduces the single-passage $D_{p}^{(\rm th)}/R_{t}$ (Diener et al., 1995). It would be interesting to revisit this prediction quantitatively with generic shapes of $\boldsymbol{\xi}$ (non-affine) under (Van Hoolst, 1994) and consider the multi-orbit effects.

In fact, the set of equations required to derive the Roche limit has been used in Yu et al. (2023) to study the frequency shift of the f-mode due to the companion’s tidal field in coalescing binary neutron stars (see also Pitre and Poisson 2025 for its relativistic extension). The two effects share the same physical origin. In particular, our analysis also predicts an effective frequency shift of the $m=2$ f-mode given by (cf. Eq. 21 but allowing for arbitrary $n$ using results from the Appendix A)

[TABLE]

where the linear, static limit of $q_{0}\simeq-\sqrt{6}(5-n)\mathcal{M}_{\ast}/4$ is used. This agrees with eq. 16 of Yu et al. (2023) in both the functional form and the numerical value, validating the affine approximation.

The analysis of Yu et al. (2023) also provides further support for why we can truncate the analysis at the three-wave order for the Roche limit. Specifically, in this work, the interacting $m=2$ and $m=0$ modes are both directly linearly driven by the tide, allowing them to directly interact at the three-wave order. This is of the same order as the frequency shift of Yu et al. (2023) and lower order than the anharmonicity and pg-instability (Weinberg et al., 2013; Weinberg, 2016; Venumadhav et al., 2014), which fundamentally originate from four-wave interaction ($\mathcal{O}(\xi^{3})$ in the equation of motion). In the pg-instability, the p-mode is not linearly forced by the tide, but instead nonlinearly driven through the three-wave channel by an f-mode and a g-mode, schematically written as $q_{p}\simeq\kappa q_{f}q_{g}$ with the subscript indicating the mode (see appendix A of Weinberg et al. 2013). The p-mode then back-reacts on the g-mode, as

[TABLE]

causing the instability when $\omega_{g,{\rm eff}}^{2}<0$ or $\kappa q_{f}q_{f}^{\ast}>1$. Because the three-wave interaction is used twice with the p-mode acting only as a mediator, the pg-three-wave interaction comes in the same order as a direct four-wave interaction between two f-modes and two g-modes, corresponding to a forcing term of the form $\zeta q_{f}^{\ast}q_{f}q_{g}$ in the equation above. See also similar discussions about the nonlinear anharmonicity studied in Yu et al. (2021) and Kwon et al. (2024, 2025) caused by effective four-wave interactions. In this work, however, the three-wave interaction is only used once because the $m=0$ mode has a linearly forced component ($q_{0}\sim q_{2}$ here instead of $q_{0}\sim\epsilon q_{2}$ as in the pg-instability), making it a lower-order effect than the anharmonicity. This can also be seen from the left panel of fig. B1 in Yu et al. (2023); compare specifically the gray (equivalent to the instability of this work) and the yellow (equivalent to the pg-instability/anharmonicity) curves. This theoretical argument is reflected in Fig. 1, and is numerically justified by Figs. 3 and 5. Including the four-wave interaction (Eqs. A33a and A33b) modifies the threshold $D_{p}^{\rm(th)}/R_{t}$ by only a few percent (Fig. 1).

Similar to Chandrasekhar (1963), our study does not directly address the evolution of a fluid element after a mode’s amplitude runs away. Two crucial questions remain to be addressed. The first is how much mass is removed from $M$ and at what rate. Whereas in circular or low-eccentricity cases, a steady stream of mass flow from the donor to the accretor may form (Shu and Lubow, 1981), in the highly eccentric cases, the mass loss should occur as discrete bursts whenever the wave breaks (after the mode amplitude runs away). An upper bound on the mass loss $\Delta M$ can be set based on the energy released, $\Delta E_{\rm mode}\sim 0.1$. If most of the energy is converted to heat (MacLeod et al., 2022), then $\Delta M$ from the wave breaking may be much smaller than $\Delta E_{\rm mode}$. Nonetheless, a large amount of heat may inflate the donor, making it more likely to lose mass in the subsequent evolution, which quickly disrupts the donor (Yu and Dai, 2024). On the other hand, if the donor’s density increases as it loses mass, it may survive the mass loss (Liu et al., 2013). One such example could be planets in the HN desert (Armstrong et al., 2020). The second question is how a mass element lost from the donor affects the orbit. The numerical simulation of Guillochon et al. (2011) suggests that the total mass and angular momentum of the system are likely to be non-conserved when the mass loss happens with a highly eccentric orbit. One may study the effects on the orbital dynamics following, e.g., Huang (1963); Sepinsky et al. (2009); Hamers and Dosopoulou (2019); Lau and Yu (2025). Alternatively, since the exterior potential of a perturbed, oscillatory star is simply $\propto Y_{2m}/D^{3}$, it is straightforward to find a test particle’s trajectory under a Roche potential modified by the quadrupolar potential of the donor (and potentially the accretor). The trajectory will determine the fate of the ejected mass and its impact on the orbit. Of particular interest are the changes in the binary’s total mass and orbital angular momentum. We plan to study these in follow-up investigations.

We thank Eric Coughlin, Yanqin Wu, Brad Hansen, Smandar Naoz, Grant Weldon, Fei Dai and Zhen Pan for useful discussions during the preparation of this manuscript. This work is supported by Montana NASA EPSCoR Research Infrastructure Development under award No. 80NSSC22M0042 and NSF grant No. PHY-2308415.

Appendix A Equations of motion of the affine model

Write

[TABLE]

where $\boldsymbol{x}$ and $\boldsymbol{r}$ are the original and perturbed locations of a mass element. The Lagrangian displacement is denoted by $\boldsymbol{\xi}$, which is further decomposed into eigenmodes with eigenfunctions $\boldsymbol{\xi}_{m}$ according to Eqs. (8) and (9). As argued in the main text, inertial modes and $m=\pm 1$ f-modes are both ignored in this work. Because $\xi_{m}$ is linear in the original position, perturbations of a homogeneous ellipsoid are therefore affine as in Carter and Luminet (1983, 1985); Diener et al. (1995). We additionally introduce a trace term $\boldsymbol{\xi}_{\rm tr}=q_{\rm tr}\boldsymbol{x}$, corresponding to a monopolar (radial with $l=0$) perturbation. While a homogeneous ellipsoid has no radial modes in the linear problem, we find $q_{\rm tr}$ is needed to account for rotation and nonlinear interactions starting at the four-wave order. Note that the trace ($l=0$) mode is normalized differently from the $l=2$ ones; $q_{\rm tr}=1$ corresponds to a surface displacement of unity. This allows us to rewrite

[TABLE]

which satisfies the form of an affine model $r_{i}(t,\boldsymbol{x})=q_{ij}(t)x_{j}$, with Einstein summation assumed and the indices standing for coordinates in Cartesian. We further use $\boldsymbol{q}$ to denote the matrix whose components are $q_{ij}$. The determinant of the matrix will be useful later and is explicitly given

[TABLE]

Note that $\boldsymbol{q}$ is diagonalized as

[TABLE]

where the angle $m\phi_{q}$ describes the phase of the mode $q_{m}$ and

[TABLE]

It shows the affine model describes an ellipsoid with semiaxes $a_{1}$, $a_{2}$, $a_{3}$. Note that $a_{1}>a_{2}$ and $a_{1}>a_{3}$. When the dynamical tide dominates, we further have $a_{2}\geq a_{3}$ in general, though it is also possible to have $a_{2}<a_{3}$ occasionally. We further use an upper bar to denote quantities evaluated in a frame co-rotating with $q_{2}$, e.g., displacement $\bar{\boldsymbol{r}}=\boldsymbol{R}\boldsymbol{r}$ and quadrupole moment

[TABLE]

The equation of motion is obtained from variations of the energies. While the affine model is self-consistently satisfied only for an incompressible, homogeneous fluid with a polytropic index $n=0$ and adiabatic index $\Gamma=\infty$, we nonetheless follow (Diener et al., 1995) (see also Carter and Luminet 1983, 1985) and allow the energies to be described by a generic pair of $(n,\Gamma)$, without requiring buoyancy neutral (i.e., $\Gamma$ can differ from $(1+1/n)$), while still assuming the perturbed fluid is described by Eq. (A2). The error due to this inconsistency will be assessed in future studies through direct comparisons with Van Hoolst (1994).

The kinetic energy can be computed from

[TABLE]

where $\mathcal{M}_{\ast}=\int x_{i}x_{i}dM/3$. Numerically, $\mathcal{M}_{\ast}=1/5$ for $n=0$ and $\mathcal{M}_{\ast}=0.1306910$ for $n=1$.

The gas energy is $U=\int(\epsilon/\rho)dM$, where $\epsilon$ and $\rho$ are respectively the internal energy density and mass density. We assume that the gas follows a polytropic equation of state, where

[TABLE]

with $P$ the pressure and $\Gamma=(\partial\log P/\partial\log\rho)_{\rm entropy}$ is the adiabatic index. Using further $\rho(\boldsymbol{r})=\rho(\boldsymbol{x})/\det[q]$, we have

[TABLE]

The explicit form of $U$ is readily obtained by plugging in $\det[\boldsymbol{q}]$ from Eq. (A3), and its derivative with respect to $q_{ij}$ is $\partial U/\partial q_{ij}=-\Pi(\boldsymbol{q}^{-1})_{ji}$. This will be used to derive the equation of motion in terms of $q_{ij}$.

To assist comparisons with Van Hoolst (1994), we also provide a series expansion in terms of $q_{m}$’s, as

[TABLE]

We have verified that the series matches exactly eq. (19) of Van Hoolst (1994) for $n=0$, in which case $\boldsymbol{\xi}_{m}$ are given by Eq. (9).

We then consider the self-gravity. For this, we introduce the self-gravitational potential energy tensor

[TABLE]

From the definition, it is easy to see that the total gravitational binding energy $\Omega$ is given by the trace of $\Omega_{ij}$, $\Omega={\rm Tr}[\Omega_{ij}]$. Because the equations of motion depends on $\partial\Omega/\partial q_{ij}=-\Omega_{ik}(\boldsymbol{q}^{-1})_{jk}$, we compute the full $\Omega_{ij}$ matrix here. While Carter and Luminet (1983) gave a formula for $\Omega_{ij}$ in their eq. (6.15), it involves solving a complicated integral in their eq. (6.4) whose explicit solution is non-trivial to obtain. We instead follow Chandrasekhar and Lebovitz (1962) to obtain a closed-form solution. For this, it is convenient to first consider the potentials in the barred frame obtained by rotating the system azimuthally by an angle $\phi_{q}$ via the $\boldsymbol{R}$ matrix. In this system, the potential energy tensor is diagonal and given by

[TABLE]

where $\bar{I}_{ij}$ are still given by Eq. (A6), and we explicitly include the $5/(5-n)$ factor to generalize the result to arbitrary values of $n$ (with the caveat that for $n>0$, $\boldsymbol{r}$ is not strictly affine). The one-index symbols $A_{i}$ are defined in eq. (5) of Chandrasekhar and Lebovitz (1962) and computed explicitly in terms of $a_{i}$ in their eqs. (15)-(17). For completeness, we reproduce them below. Define

[TABLE]

and denote incomplete elliptic integrals of the first and second kinds, respectively, as

[TABLE]

we can then write

[TABLE]

Note that the results are obtained by swapping the subscripts $(2,3)$ in the corresponding equations in Chandrasekhar and Lebovitz (1962). This is because in our coordinate system, $a_{1}>a_{3}>a_{2}$ most of the time, though occasionally, it is possible to have $a_{2}>a_{3}$. In this case, Eq. (A14) still holds as long as one lets $\sin^{2}\theta_{e}>1$ and uses the reciprocal modulus transformation of elliptical integrals to handle $m>1$ (Byrd and Friedman, 1954)

[TABLE]

In the limiting case where $a_{2}=a_{3}$ (achieved in the adiabatic tide limit), we have

[TABLE]

where $a_{21}=a_{2}/a_{1}$.

Once $\bar{\Omega}_{ij}$ is obtained, $\Omega_{ij}$ can be readily computed by rotating back to the original, unbarred frame, as

[TABLE]

and the total gravitational energy is

[TABLE]

where we again expanded the exact solution in powers of $q_{m}$ in the last equality to assist with a direct comparison with Van Hoolst (1994). For a homogeneous ellipsoid, discontinuities in the derivatives of the gravitational potentials at the surface are non-trivial to evaluate, and we present an explicit evaluation of those terms in Appendix B.

We consider both the tide and rotation as external perturbations to the system. For the tide, the interaction energy truncated to the quadrupolar order can be written as

[TABLE]

where

[TABLE]

Note that we work in the frame corotating with the donor and use $\Phi_{C}=\Phi-\omega t$ to account for the Doppler effect. In comparison to, e.g., eq. (B1) of Diener et al. (1995), we additionally included the centrifugal force due to a rotation $\omega$ (assumed to be along the $z$ axis, also the direction of the orbital angular momentum), giving rise to an energy

[TABLE]

By varying the energies, the equations of motion are

[TABLE]

By further contracting Eq. (A22) with the symmetric-tracefree tensor \mathcal{Y}_{lm}^{ij}{\color[rgb]{0,0,0}=(\mathcal{Y}_{l,-m}^{ij})^{\ast}} (Poisson and Will, 2014), with

[TABLE]

we obtain the equations of motion for each spherical harmonic mode.

[TABLE]

where on the right-hand side, the driving terms are due to variations of, respectively, the internal energy, the gravitational energy, the tidal interaction, and the centrifugal energy. Explicit expressions are provided as follows. By varying the internal energy, we have

[TABLE]

And from the gravity,

[TABLE]

Note that the equations above are exact up to the affine approximation, Eq. (A2). We also present expanded forms of the driving terms that are complete to the four-wave interaction level. We will further write $q_{\rm tr}=q_{\rm tr}^{(1)}+q_{\rm tr}^{(2)}$ where $q_{\rm tr}^{(1)}\sim\omega^{2}\sim\epsilon$ will be sourced by the ($l=0$ part of the ) centrifugal potential, while $q_{\rm tr}^{(2)}\sim\epsilon^{2}$ from nonlinear mode interactions. From the internal pressure,

[TABLE]

And from self-gravity,

[TABLE]

We now turn to external perturbations. For the tidal interaction, we have

[TABLE]

Similarly, the centrifugal drive is decomposed as

[TABLE]

Since the affine model is self-consistent only for homogeneous ellipsoids with $(n=0,\Gamma=\infty)$, we discuss the equations of motion for this case with further details. First note that the natural frequency of the $l=0$ mode is \omega_{\rm tr}^{2}=(3\Gamma-4)/(5-n){\color[rgb]{0,0,0}/\mathcal{M}_{\ast}}, which approaches infinity as $\Gamma\to\infty$. Therefore, we can approximate $\ddot{q}_{\rm tr}=0$, allowing us to solve for the trace part as

[TABLE]

We have expanded the results in $1/\Gamma$ and kept only the leading-order terms. Substituting them back into the equations of motion for quadrupolar modes, we obtain the complete description of the dynamics for $(n,\Gamma)=(0,\infty)$ at the four-wave order up to the ignorance of the Coriolis force (which can be included at the leading order following Eqs. 15-17; Yu and Lau 2025) which reads

[TABLE]

where,

[TABLE]

Note that terms proportional to $\Gamma$ cancel out exactly in $\Pi_{m}^{(3)}$, making the result regular when $\Gamma\to\infty$.

Appendix B Surface terms in the third order perturbations of the gravitational potential energy

At the formulation level, it is clear that the affine model agrees with that in Van Hoolst (1994) when applied to incompressible stars. However, it is not trivial to show the agreement, especially for terms of $\mathcal{O}(\xi^{3})$ and above. This is because these terms in Van Hoolst (1994) often involve high-order derivatives of the quantities from the background or the lower perturbation levels, which produce surface terms that require integrating the Dirac delta functions. In this appendix, we explicitly calculate the $\mathcal{O}(\xi^{3})$ terms in the gravitational potential energy and show their agreement with the affine model under careful treatment of the surface terms.

We first write the $\mathcal{O}(\xi^{3})$ terms as

[TABLE]

where

[TABLE]

where the integral is carried within the unperturbed volume. The Lagrangian displacement is expressed in spherical coordinates

[TABLE]

The unperturbed potential $\Omega_{0}$ in all space is given by

[TABLE]

where $\Theta(x)$ is the Heaviside step function666The one-half convention of the Heaviside step function is the result from treating both $\delta(x)$ and $\Theta(x)$ as finite distribution functions (with even and odd parity respectively) and taking the zero width limit for both of them. This gives $\int_{-\infty}^{\infty}dx\delta(x)\Theta(x)=1/2$. with $\Theta(0)=1/2$. We also replace $\rho$ by $\rho\Theta(R-r)$ since the surface terms depend on the boundary value. The triple derivative contributes only to a Dirac delta term within the domain of integration:

[TABLE]

The Lagrangian displacement expanded in terms of the quadrupolar f-mode eigenfunctions is given by Eq. (8). The integral of three spherical harmonics can be evaluated using the Wigner 3-j symbols (see, e.g., Weinberg et al. (2012)):

[TABLE]

The non-zero coefficients for $l=2$, $m=-2,0,2$ are

[TABLE]

Focusing on the $l=2$, $m=-2,0,2$ modes, we have

[TABLE]

The integral in $\Delta\Omega^{(3b)}$ contains both a surface integral term from the Dirac delta function and a bulk term.

[TABLE]

where the subscripts s and b denote the surface term and bulk terms, respectively. To see the delta function arising from the derivatives of the Eulerian potential, we write the quadrupolar f-mode $\Omega^{\prime}$ as

[TABLE]

Hence, $\partial_{r}^{2}\Omega^{\prime}$ gives a term with the Dirac delta in the integrand.

[TABLE]

For $l=2$, $m=-2,0,2$ modes, and using Eq. (B7), we have

[TABLE]

Note that the delta term within $\partial_{r}^{2}\Omega^{\prime}$ is equivalent to $-\xi^{r}\partial_{r}^{3}\Omega_{0}$. Therefore, the factor of three in Eq. (B13) comes from the ratio of the prefactors $1/6$ and $1/2$ of the potential terms.

The remaining terms without the Dirac delta function can be integrated to give

[TABLE]

Collecting Eqs. (B9), (B13) and (B14), we obtain

[TABLE]

which agrees with Eq. (A18).

Appendix C Coefficients of polynomials

We present here the explicit expressions for the numerical coefficients of Eqs. (26), (29), and (30).

For Eq. (26), we have

[TABLE]

where $a_{31}=a_{3}/a_{1}$, and

[TABLE]

Note Eq. (19) has been used to write $C_{2}^{(\rm st)}$.

For Eq. (29), we have

[TABLE]

and

[TABLE]

Lastly, for Eq. (30), the coefficients are

[TABLE]

and

[TABLE]

Again, Eq. (19) has been used to eliminate the Coriolis term, $C_{2}^{(\rm st)}$.