Dynamics of Relativistic Binaries in Structured and Stochastic Environments: A Lagrange-Fourier-Hansen Framework

Abstract: We develop a general framework to characterize non-vacuum perturbations to relativistic binaries in the gravitational-wave (GW) driven regime, for use in GW parameter estimation studies. The effect of smooth, structured and stochastic perturbations to the binary's motion is reduced to a resonant spectral projection defined on a rolling averaging window, with weights given by Hansen coefficients. This is combined with practical criteria for identifying and evaluating the corresponding dynamical response to perturbations, starting from either analytical models or numerical simulations of binaries in environments. The result is a set of coupled ODEs for the orbital elements that capture epi-cyclic, apsidal and nodal resonances, consistently incorporate feedback from radiation reaction and can be solved efficiently on a coarse time grid. We demonstrate the practical application of the framework in two representative astrophysical scenarios: a compact binary in a variable tidal field and an extreme-mass-ratio inspiral in an accretion disk. We propose the Lagrange-Fourier-Hansen framework as a unified tool for modeling environmental effects in GW templates for eccentric and precessing binary sources, and particularly for bridging the gap between phenomenological prescriptions and realistic models of binaries in environments.

• The paper presents the Lagrange–Fourier–Hansen (LFH) framework, a novel approach to model GW phase modulations from environmental forces in relativistic binary systems.
• The LFH framework incorporates a method to break down complex environmental forces into resonances with binary dynamics, enhancing the precision of GW waveform templates.
• In real-world scenarios with moderately eccentric binaries, the LFH method reduces the computation cost significantly over direct simulation methods, offering efficiency gains beneficial for analyzing parameter space in GW astronomy.

Motivation and scope

Gravitational-wave (GW) astronomy has reached a level of waveform fidelity in which the binary phase is mapped to orbital dynamics with high precision, but this precision has been achieved almost exclusively under the assumption of vacuum evolution. Astrophysical GW sources—stellar-mass binaries assembled dynamically, massive black hole binaries, and extreme-mass-ratio inspirals (EMRIs)—form in dense environments (galactic nuclei, stellar clusters, accretion discs) that can remain coupled to the binary throughout the observable inspiral. The dominant methodology for modeling such environmental effects (EEs) has been phenomenological: power-law corrections to the GW phase of the form $\delta\psi_{\rm GW}\sim f_{\rm GW}^{n}$. Zwick, Dyson, Seymour, Takátsy, and Samsing argue that this approach is fundamentally limited for two reasons: eccentricity and precession make relativistic binary dynamics intrinsically multi-frequency, and realistic environmental perturbations are structured, time-variable, or stochastic with broadband spectra rather than smooth secular torques. Their paper develops the Lagrange–Fourier–Hansen (LFH) framework (2606.27526), a general and computationally efficient map from an arbitrary perturbing force—including forces extracted directly from numerical simulations—to the orbital evolution and GW observables of an eccentric, precessing binary.

Structure of the framework

The LFH framework combines three classical ingredients: the Lagrange Planetary Equations (LPE) in Gauss form for perturbed osculating elements $\mathbf{Y}=(p,e,\omega,\Omega,\iota,M)$; windowed Fourier spectral analysis of the forcing; and Hansen coefficients $X_l^{jk}(e)$ to expand the orbital dependence of the LPE into harmonics of the mean anomaly. The carrier orbit is approximated as a rigidly precessing Keplerian conic section whose elements evolve through orbit-averaged radiation reaction and conservative precession. The derivation proceeds in four steps:

1. Adiabatic windowing. A local time window $\mathcal{I}_\tau$ of duration $\Delta t_W$ satisfying $T_{\rm orb}\ll \Delta t_W \ll T_{\rm rr}$ isolates the instantaneous response from the slow radiation-reaction-driven inspiral. Within each window, orbital phases evolve linearly in time.
2. Fourier–Hansen decomposition. The projected accelerations $(R,S,W)$ are Fourier transformed within each tapered window, and the LPE are expanded in Hansen coefficients, yielding phases $\Phi_l^{(q)}=\sigma t + lM_{\rm c}+q\omega_{\rm c}$ labeled by epicyclic index $l$ and apsidal sideband $q=0,\pm1$.
3. Resonance selection. Window averaging acts as a spectral sinc filter; in the long-window limit it becomes a delta function selecting resonant frequencies $\mathbf{Y}=(p,e,\omega,\Omega,\iota,M)$0. The residual phase $\mathbf{Y}=(p,e,\omega,\Omega,\iota,M)$1 tracks phase alignment between forcing and orbit across windows.
4. Stationary-phase integration. Resonant crossings over radiation-reaction timescales are evaluated via the stationary-phase approximation (SPA), where stationarity requires $\mathbf{Y}=(p,e,\omega,\Omega,\iota,M)$2, each stationary point contributing a kick proportional to the coherence time $\mathbf{Y}=(p,e,\omega,\Omega,\iota,M)$3.

The mean anomaly response is treated separately and decomposed into three physically distinct contributions: a Keplerian term accumulating historical semi-major-axis drift, a direct radial-forcing term modifying angular velocity at fixed geometry, and a geometric term from reorientation of the osculating ellipse. The final output is a closed system of coupled ODEs for the full ("renormalised") orbital elements, combining vacuum radiation reaction and precession with the resonant environmental terms, so that feedback between accumulated perturbations and the GW-driven inspiral is captured self-consistently. The framework also handles center-of-mass acceleration through a Rømer-delay timing correction, and converts time-domain phase shifts into frequency-domain dephasings per GW harmonic via $\mathbf{Y}=(p,e,\omega,\Omega,\iota,M)$4, enabling direct use in frequency-domain templates.

Practical workflow and worked examples

The paper demonstrates the pipeline on two representative systems: a $\mathbf{Y}=(p,e,\omega,\Omega,\iota,M)$5 binary embedded in a rotating tidal field (with $\mathbf{Y}=(p,e,\omega,\Omega,\iota,M)$6 chosen to exhibit an apsidal co-rotation resonance), and a $\mathbf{Y}=(p,e,\omega,\Omega,\iota,M)$7 EMRI around a $\mathbf{Y}=(p,e,\omega,\Omega,\iota,M)$8 black hole interacting with an accretion disc, using torque time series extracted from the hydrodynamic simulations of Derdzinski et al. Mode truncation is controlled by the Hansen scaling $\mathbf{Y}=(p,e,\omega,\Omega,\iota,M)$9, so near-circular binaries ($X_l^{jk}(e)$0 for the EMRI example) require only the lowest harmonics while moderately eccentric systems ($X_l^{jk}(e)$1) retain modes up to $X_l^{jk}(e)$2.

Spectrograms of the windowed force reveal which modes carry power as the inspiral progresses; in the tidal-field example, a component of the forcing sweeps through zero frequency ("DC"), at which point the $X_l^{jk}(e)$3 mode responds strongly, producing a resonant kick near orbit 1520. For the fully analytic tidal case, an additional Hansen expansion of the perturbation itself yields closed-form commensurability conditions—for instance, the $X_l^{jk}(e)$4 mode resonates precisely when $X_l^{jk}(e)$5 (apsidal co-rotation)—and the relative phase reduces to $X_l^{jk}(e)$6, confirming that the Fourier phase cancels for perturbations depending only on the instantaneous binary state.

The most consequential practical result concerns computational cost: because the ODE system evolves on the coarse $X_l^{jk}(e)$7 grid rather than the orbital timescale, the resonance kick is accurately reproduced with an integration step every roughly 80 orbits, whereas direct integration of the equations of motion requires several steps per orbit. Alternatively, SPA kicks implemented as smooth Gaussian ramps reproduce the fine-grid evolution analytically. This efficiency gain is what makes the framework viable for parameter-inference studies, where many waveform evaluations are required. The authors note, however, that systematic benchmarking of window-length choices, grid resolution, and stationary-point truncation criteria is deferred to a forthcoming release of their Python package, AstroWaveforms.

Application landscape

The discussion section surveys where structured or stochastic forcing is expected to matter. Tidal perturbations from tertiary companions produce line-of-sight accelerations, secular dephasing, Kozai–Lidov oscillations, and—beyond orbit-averaged treatments—time-dependent resonances from outer-orbit harmonics, with richer structure for eccentric outer orbits and finite-size horizon effects. Gas-embedded binaries are emphasized as a case where simple migration-torque or drag prescriptions fail: spectral analyses of hydrodynamic simulations show substantial torque power at multiples of the orbital frequency, and in thinner discs the binary response can be dominated by resonant coupling to fluctuations rather than by the mean torque. For EMRIs, estimates suggest a substantial fraction of LISA-detectable sources may form in AGN discs, and the quasi-periodic-eruption class has been interpreted as a precursor stage to gas-embedded EMRIs; long-lived phase coherences of $X_l^{jk}(e)$8–$X_l^{jk}(e)$9 orbital cycles observed in the simulated torque imply numerous stationary points and hence frequent resonant kicks. Related applications include dark-matter and ultralight-boson environments, where recent work shows multiple resonances significantly alter eccentric EMRI evolution relative to circular orbits, and modified-gravity scenarios with intrinsically time-dependent effective forces.

Limitations and open questions

The framework's principal structural limitation is its orbital parametrization. At 1PN order the quasi-Keplerian (QK) description is exactly a rigidly precessing conic section, so the resonant formalism carries over directly once QK parameters are adopted; however, beginning at 2PN order additional oscillatory terms ($\mathcal{I}_\tau$0, $\mathcal{I}_\tau$1, etc.) generate harmonics absent from the present treatment, which could in principle participate in resonances. Whether these higher-order corrections are observationally relevant depends on resonance strength and source signal-to-noise ratio, and remains unquantified here. For EMRIs, a fully relativistic extension via osculating geodesics and Mino-time frequencies $\mathcal{I}_\tau$2 is outlined but not yet developed; connecting generic environmental forcing to locally integrated torques in Kerr spacetime is work in progress by some of the authors. Additional caveats stated explicitly include: the EMRI torque series was computed for circular orbits and assumed representative for $\mathcal{I}_\tau$3; the tidal example considered only a circular outer orbit; and optimal windowing, truncation, and SPA-kick automation await the AstroWaveforms benchmarking study. Finally, the authors note that applying the results to GW observables requires re-expanding the vacuum radiation-reaction terms in terms of waveform phase to account for conservative frequency shifts induced by the environment—an implementation detail not carried out in this paper.

Conclusion

The LFH framework provides a systematic reduction of arbitrary environmental forcing—analytic, structured, or stochastic, including simulation output—to a sparse set of resonantly coupled orbital modes, expressed as efficient coupled ODEs on a coarse time grid with optional stationary-phase treatment of sharp resonances. Its central contribution is to replace single-parameter secular dephasing prescriptions with a mode-resolved treatment appropriate to eccentric and precessing binaries, thereby enabling waveform templates that retain the full dynamical interaction between relativistic binaries and their surroundings. The immediate open problems are the extension to higher-order quasi-Keplerian and self-force-based parametrizations, and the quantitative benchmarking needed before the method can be deployed in production-level parameter estimation.