1 Introduction↩︎

Neutron stars provide a unique arena to study matter under extreme conditions [1]. Their cores can reach densities exceeding those of atomic nuclei, making them natural laboratories to investigate the behavior of strongly interacting matter and its connection to nuclear and particle physics. Understanding their internal structure is therefore a central goal in modern astrophysics.

The macroscopic properties of neutron stars are governed by the equation of state (EoS), which encodes the relation between pressure, density, and temperature and depends sensitively on the underlying microphysics [1], [2]. The EoS determines observable quantities such as the mass-radius relation, the maximum mass, and the tidal deformability. The advent of gravitational-wave (GW) astronomy has opened a direct observational window onto these properties: signals from coalescing compact binaries carry imprints of the internal structure of neutron stars, allowing us to probe fundamental physics in regimes inaccessible to terrestrial experiments (see Refs. [3], [4] for reviews).

In a binary system, each neutron star experiences the gravitational field generated by its companion, which induces multipolar deformations [5]–[9]. These deformations are commonly characterized by the tidal Love numbers (TLNs), which relate the external tidal field to the induced multipole moments. Tidal interactions affect both the conservative dynamics of the binary and the emitted GW signal, leaving measurable imprints on the waveform [10], [11]. Since the first binary neutron-star observation GW170817 [12], tidal effects have become a key probe of dense-matter physics and compact objects (see Refs. [13], [14] for recent overviews).

Most studies to date have focused on the static (or adiabatic) tidal response, in which the induced multipole moments instantaneously track the external tidal field. This approximation is valid when the orbital timescale is much longer than the characteristic internal timescales of the star. As the binary inspirals toward merger, however, the orbital frequency increases and the tidal response acquires a nontrivial frequency dependence, giving rise to genuinely dynamical tidal effects. These include resonant excitations of stellar oscillation modes, dissipative phenomena, and phase lags between the tidal field and the induced deformation, all of which can modify the late-inspiral GW signal.

The physics of dynamical tides has a long history in Newtonian gravity, where the tidal response of a self-gravitating fluid body can be described in terms of the excitation of its normal modes [15]–[25]. In this framework, the frequency-dependent tidal deformability admits a mode-sum representation in which the dominant contribution typically arises from the fundamental fluid (\(f\)-) mode because of its relatively low frequency and large overlap with the external tidal field [15], [17], [19]. Even away from resonance, the coupling to the \(f\)-mode can produce sizable corrections to the static Love number, motivating effective descriptions in which the static tidal deformability is supplemented by more resonant contributions.

Although Newtonian models provide valuable physical intuition, a quantitatively accurate description of tidal interactions in compact binaries ultimately requires a relativistic treatment. Considerable progress in this direction has been achieved in recent years. A fully relativistic analysis of the dynamical tidal response of nonrotating stars was recently presented in Refs. [26]–[30]. The relativistic mode structure underlying the dynamical response has also been investigated in detail in Refs. [31], [32], where it was compared to the Newtonian overlap-integral formalism.

An important outcome of these studies is that the relativistic mode expansion differs qualitatively from its Newtonian counterpart. In particular, the operators defining the natural inner product for relativistic stellar perturbations are not positive definite, implying that the associated mode sums are not strictly convergent [31], [32]. This feature reflects the fact that, in general relativity, the dynamical tidal response cannot be fully decomposed in terms of a sum over normal-mode excitations. Closely related issues emerge when comparing the exact relativistic tidal response with its low-frequency expansion [31]: while the expansion reproduces the weak-field regime, sizable deviations can arise for sufficiently compact stars, signaling the increasing importance of genuinely relativistic effects.

The relativistic formulation of dynamical tides has also been developed from the perspective of worldline effective field theory (EFT) [33]–[35]. In this framework, the internal structure of a compact object is encoded in frequency-dependent response functions entering a point-particle effective action. Dynamical tidal effects in compact binaries have been investigated in EFT approaches in Refs. [36]–[42]. This formalism provides a coordinate-independent and systematically improvable description of finite-size effects and is particularly suited to connect stellar perturbation theory with post-Newtonian (PN) and waveform calculations. Furthermore, the EFT calculation clarifies the appearance of logarithmic terms characterizing the running of the dynamical Love numbers [36], [38], [39], [42]–[52] and highlighting the necessity of a rigorous matching procedure. As we discuss below, a complete matching has not yet been performed for the dynamical Love numbers of a neutron star and this would be one of our primary goals.

Dissipative aspects of the tidal response have also recently received considerable attention. Frequency-dependent dissipative tidal deformabilities, including the effects of viscosity and mode damping, were investigated in Refs. [27]–[30]. These works showed that dissipative effects are generally much smaller than the conservative dynamical response for realistic neutron-star models, although they remain sensitive to the internal composition and transport properties of dense matter. At the same time, they highlighted that the conservative dynamical tidal response can vary significantly with the EoS and stellar compactness.

From a PN perspective, dynamical tidal effects formally enter the GW phase at high order, comparable to other subleading finite-size contributions [26], [27], [36], [37], [39], [47], [50], [51], [53]–[55], including nonlinear tidal effects beyond linear perturbation theory [56]. Nevertheless, these contributions can be enhanced by resonant mode excitations, inverse compactness scaling, and by relativistic strong-field effects, making them potentially observable with future detectors. Accurate modeling of dynamical tides is therefore essential both to avoid systematic biases in parameter estimation and to fully exploit the scientific potential of third-generation GW observatories [57]–[60].

In this work we address the following question: How does the frequency dependence of the tidal response of relativistic neutron stars affect the dynamics of compact binaries and their GW emission? While previous studies established the qualitative importance of dynamical tides, a complete and systematic treatment connecting relativistic stellar perturbation theory, EFT response functions, binary dynamics, all the way down to GW observables is still lacking. This is particularly relevant in the presence of running terms,¹ as is the case for dynamical tidal effects in compact objects, and to avoid ambiguities that inevitably plague perturbation-theory calculations, if the latter are not properly matched to GW observables (see Refs. [13], [14] for recent reviews).

In this work, we complete this program. First, we compute the frequency-dependent tidal response of neutron stars by directly solving the relativistic perturbation equations, both perturbatively and nonperturbatively in frequency, without relying on a mode decomposition. Second, building on Ref. [50] (see also Ref. [49]), we match the resulting response to the worldline EFT at loop order, thereby obtaining for the first time the dynamical Love numbers including both the logarithmic running and the scheme-dependent finite contributions in dimensional regularization. Crucially, matching the EFT to the underlying general-relativistic neutron-star dynamics in order to obtain the complete Love-number couplings cannot be achieved at tree level within the EFT. Instead, it requires the evaluation of higher-loop worldline diagrams at the classical level. These diagrams are formally divergent and therefore require regularization and renormalization, including the introduction of a subtraction prescription and an associated renormalization scale—a structure that closely mirrors the renormalization procedure and renormalization-group running familiar from quantum field theory. Matching the results of stellar perturbation theory to the EFT couplings enables the construction of observable quantities that are independent of the choice of coordinates and perturbation variables.

We further compare the full relativistic response with resonant models based on the dominant \(f\)-mode, quantifying the regime of validity of these approximations in agreement with recent results [32]. Finally, after matching the stellar-perturbation theory response with the EFT action, we derive the complete leading-order dynamical tidal corrections to the conservative dynamics and GW signal of compact binaries. Since the tidal contributions to the PN waveform can themselves be derived from the EFT action [11], [68], fixing the split between the tidal field and the response through matching to the EFT guarantees that the coefficients we compute are precisely those entering the waveform. We then assess their detectability with current and future GW detectors for realistic neutron-star EoS. Interestingly, we find that—due to its strong enhancement at moderate compactness—the 8PN contribution coming from the dynamical Love number is relevant for parameter estimation with third-generation interferometers such as the Einstein Telescope (ET) [60].

Outline: The paper is organized as follows. In Sec. 2 we introduce the worldline EFT framework and define the static and dynamical Love numbers through the finite-size effective action. In Sec. 3 we describe the relativistic stellar perturbation theory, presenting both a perturbative expansion and a nonperturbative treatment in the frequency. In Sec. 4 we match the interior and exterior solutions, extract the tidal Wilson coefficients for a range of realistic EoS, and discuss universal relations and the connection to the \(f\)-mode resonance model. In Sec. 5 we derive the complete 8PN tidal waveform phase and assess the observability of the dynamical tide with LIGO-Virgo-KAGRA (LVK) O4 and ET [60]. We conclude in Sec. 6. Technical details are collected in four appendices: the EFT computation of the graviton one-point function (Appendix 7), details of stellar perturbation theory (Appendix 8), the exterior Zerilli matching (Appendix 9), and a comparison with the tidal constants used in the literature (Appendix 10).

Conventions: We adopt the mostly-plus metric signature, \((-,+,\cdots,+)\), and work in natural units with \(\hbar = c = 1\). The reduced Planck mass is defined as \(M_\mathrm{Pl}= (8 \pi G)^{-\frac{1}{D-2}}\), with \(G\) Newton’s gravitational constant and \(D\) the number of spacetime dimensions. Although our main interest is in \(D=4\), we keep \(D\) generic in certain sections where it is needed for dimensional regularization in the worldline effective theory. In Sec. 3, we will also set \(G=1\) to simplify the notation. Our curvature conventions is \({R^\rho}_{\sigma\mu\nu}=\partial_\mu\Gamma^\rho_{\nu\sigma}+\dots\) and \(R_{\mu\nu}={R^\rho}_{\mu\rho\nu}\). Fields are decomposed in spherical harmonics and Fourier modes as \(\psi(t,r,\theta,\varphi) = \sum_{\ell,m} \int \frac{\mathrm{d}\omega}{2\pi} \mathrm{e}^{-i\omega t} \, \psi^{\ell m}(\omega,r) \, Y_{\ell m}(\theta,\varphi)\) and, for brevity, we will often suppress the arguments of \(\psi^{\ell m}\) when context suffices. We use Greek indices, \(\mu, \nu, \rho, \dots\), to label spacetime coordinates, while Latin indices, \(i, j, k, \dots\), refer to spatial coordinates.

2 Worldline Effective Field Theory↩︎

We define the dynamical Love numbers of the neutron star within the point-particle EFT framework [33]–[35], which provides a robust and gauge-invariant description of an object’s tidal response. In this work, we focus on non-rotating, perfect-fluid stars, leaving the study of spin and dissipative effects to future work.

The EFT of a non-rotating particle takes the following form (see, e.g., Refs. [69]–[72] for some reviews): \[S=S_{\mathrm{pp}}+S_{\text{bulk }}+S_{\text{int}} \,. \label{eq:ppEFT0}\tag{1}\] Here, \(S_{\mathrm{pp}}\) is the point-particle action, given by the standard expression \[\begin{align} S_{\rm pp} = -M \int {\mathrm d}\tau = -M \int {\mathrm d}\sigma \sqrt{-g_{\mu\nu}(X)\frac{{\mathrm d}X^\mu}{{\mathrm d}\sigma}\frac{{\mathrm d}X^\nu}{{\mathrm d}\sigma}} \label{eq:Spp} \end{align}\tag{2}\] where \(\tau\) is the proper time measured along the particle’s worldline, \(M\) the particle’s mass, \(X^\mu(\sigma)\) its spacetime trajectory, and \(\sigma\) an affine parameter describing the worldline.

The action \(S_{\mathrm{bulk}}\) governs instead the dynamics of the gravitational field in the bulk spacetime and is given by the standard Einstein–Hilbert term, \[S_{\mathrm{bulk}} = \frac{M_{\mathrm{Pl}}^{2}}{2} \int \mathrm{d}^4 x \, \sqrt{-g} \, R \,.\] The object is a point particle only at leading order. To describe tidal responses, we introduce a series of higher-dimensional operators \(S_{\mathrm{int}}\) that account for the object’s finite size, evaluated along its worldline and coupled to the bulk fields [33], [36], [44], [50], [73], [74]: \[\begin{align} S_{\mathrm{int}} & = \int \mathrm{d}\tau \, Q_E^{ij}(\tau) \, E_{ij} \\ & \quad + \text{magnetic} + \text{higher orders} \,, \end{align} \label{SintE}\tag{3}\] where \(E_{ij}\) corresponds to the electric component of the Weyl tensor \(C_{\mu\nu\rho\sigma}\), \[E_{ij} \equiv C_{0i0j} \,.\] The quantity \(Q_E^{ij}(\tau)\) in Eq. 3 represents a composite operator, constructed from the internal degrees of freedom, and encodes all information about the microscopic interior dynamics, including tidal deformations, absorption, internal hydrodynamic modes, resonances, and any other finite-size effects.²

Note that in the second line of Eq. 3 we have omitted operators involving the magnetic component of the Weyl tensor. Because of the relative velocity suppression between the magnetic and electric tidal fields, magnetic tidal effects arise at higher PN order compared to their electric counterparts. For the same reason, additional higher-derivative terms involving gradients of \(E_{ij}\) have also been neglected in Eq. 3 , as they are subleading in the PN expansion. Since this work focuses on the dominant dynamical quadrupolar relativistic response of neutron stars, we will restrict our attention to the leading operator 3 , leaving a more general analysis, including subdominant terms, for future work.

The process we aim to compute is the following: we probe the point particle with an external gravito-electric tidal field and then determine the induced response field. Since a long-distance observer cannot resolve the detailed internal microdynamics of the object, we integrate out its internal degrees of freedom. Operationally, this amounts to solving for the composite operator \(Q_E^{ij}(\tau)\) within response theory. In the following, we will restrict ourselves to linear response,³ but retain frequency-dependent corrections. In addition, since our focus here is on perfect-fluid (non-rotating) neutron stars, for which viscous effects vanish, we will retain only the conservative part of Eq. 3 .⁴

In this limit, the most general quadratic finite-size action for the gravito-electric conservative sector, up to quadratic order in frequency, boils down to \[S_{\mathrm{int}} = \frac{1}{2} \int \mathrm{d}\tau \left( \frac{R_\star^5}{G}c_E E_{ij}E^{ij} + \frac{R_\star^8}{G^2M}c_{\dot{E}} \dot{E}_{ij}\dot{E}^{ij} \right)\,, \label{eq:Sintmain}\tag{4}\] where \(R_\star\) is the radius of the object, and \(c_{{E}}\) and \(c_{\dot{E}}\) are the (dimensionless) static and dynamical quadrupolar Love number couplings, respectively.

The action 1 can be used to compute any long-distance observable involving the point particle and the gravitational field. In this work, we focus on an off-shell quantity, namely the graviton one-point function induced by an external quadrupolar tidal field \(\bar{E}_{ij}\) coupled to the point particle. This response can be obtained by solving the equations of motion derived from the action 1 or, more generally, by using the Schwinger–Keldysh in-in formalism (which allows one to incorporate dissipative effects when present). The relevant formalism is reviewed in Appendix 7; see also Ref. [50]. The result for the graviton one-point function at tree-level (in \(D = 4\)) reads \[\begin{gather} \label{ErrD4} E_{rr} = -2\mathrm{e}^{-i\omega t}c_{\mathrm{ext}} Y_{2m} \\ \times\left[1 + 36 \, \left(\frac{R_\star}{r} \right)^{5} \left(c_E+ \omega^2 \frac{R_\star^3}{GM}c_{\dot{E}}\right)\right] \, , \end{gather}\tag{5}\] where \(\omega\) is the perturbation’s frequency and \(c_{\mathrm{ext}}\) characterizes the amplitude of the external (quadrupole) tidal field, located at a distance \(r\) from the object. Eq. 5 captures the leading distance behavior of the gravito-electric field in terms of the object’s tidal couplings.

To match it to the full relativistic solution, however, we must include the coupling to gravity, which generates subleading terms and induces the appearance of ultraviolet divergences, which must be regularized. As in standard quantum field theory, these divergences are absorbed by introducing suitable counterterms and adopting a renormalization scheme. Throughout this work, we employ dimensional regularization together with minimal subtraction [50]. As an intermediate step toward obtaining the renormalized quadrupolar response, we compute the Zerilli field, which satisfies the Zerilli equation [see Eq. 75 ] \[\left(\frac{{\mathrm d}^2}{{\mathrm d}r^2}-\frac{(2-\varepsilon)(3-\varepsilon)}{r^2} \right)\Psi_{\text{Z}}(r) = V_{\Psi_\text{Z}}(r)\Psi_\text{Z}(r)\,, \label{eq:Z95EWEmain}\tag{6}\] in the bulk, where \(\varepsilon=2-D/2\) and \(V_{\Psi_\text{Z}}\) is the effective potential capturing the nonlinearities of the Einstein–Hilbert theory. After resumming the gravitational nonlinearities through a Born-series expansion [49] and canceling the ultraviolet divergences with the appropriate counterterms, we obtain the renormalized Zerilli field [see Eq. 77 in Appendix 7, and Ref. [50] for details]: \[\begin{gather} \Psi_\text{Z}^{\mathrm{R}}=\, \bar{B}_\text{reg} \left\{r^3 \left[1-\bar{G}^2 \omega ^2 \left(\frac{214}{105} \log (\mu r)+\frac{2731}{9800}\right)\right] \right. \\ \left.+\frac{\bar{G}^7 \omega ^2}{r^2}\left(\frac{3011}{80} \log (\mu r)+\frac{1545209}{57600}\right)-\frac{189 \bar{G}^5}{32r^2}\right\}\\ +\frac{\bar{B}_\text{irr}}{r^2} \left[\bar{G}^2 \omega ^2 \left(\frac{214}{525} \log (\mu r)+\frac{111383}{441000}\right)+\frac{1}{5}\right]\,, \label{eq:psiZRappmain} \end{gather}\tag{7}\] where \(\bar G \equiv GM\), \(\mu\) is the renormalization scale, and \(\bar B_\text{reg}\) and \(\bar B_\text{irr}\) are the (frequency-dependent) renormalized integration constants of the flat-space Zerilli equation. Equation 7 can then be used to construct the gravito-electric component \(E_{ij}\) of the Weyl tensor [see Eq. 73 of Appendix 7]. Matching this result to the one-point function computed from the EFT action in Eq. 5 allows us to relate the renormalized tidal couplings \(c_E\) and \(c_{\dot{E}}\) to the coefficients \(\bar B_{\rm reg}\) and \(\bar B_{\rm irr}\) via [see Eqs. 64 and 74 ]:⁵ \[c_E + \omega^2 \frac{R_\star^3}{GM}c_{\dot{E}} = \frac{1}{180R_\star^5} \frac{\bar B_\text{irr}(\omega)}{ \bar B_\text{reg}(\omega)}\,, \label{eq:ccBB}\tag{8}\] where, consistently with our perturbative framework, the right-hand side should be expanded to quadratic order in \(\omega r_s\), where \(r_s=2GM\) is the Schwarzschild radius of the star. Once the Zerilli solution 7 is matched to the relativistic solution for \(\Psi_\text{Z}\), and the coefficients \(\bar B_\text{reg}\) and \(\bar B_\text{irr}\) are fixed in terms of the parameters of the full theory, the Love numbers \(c_E\) and \(c_{\dot{E}}\) are completely determined by the neutron star’s properties and EoS via Eq. 8 . This matching will be discussed in Sec. 4.1.

3 Stellar perturbation theory↩︎

In this section, we describe time-dependent, linear relativistic stellar perturbation theory. We employ two complementary approaches: one which is nonperturbative in the frequency and is valid to any order in \(M\omega\) (Sec. 3.2) and another which is perturbative in the frequency and valid up to quadratic order (Sec. 3.3). We set \(G=1\) throughout this section.

3.1 Background configuration↩︎

We consider as a background a static and spherically symmetric configuration. In Schwarzschild-like coordinates \((t,r,\theta,\varphi)\), the line element takes the form \[\begin{align} g_{\mu\nu} {\mathrm d}x^\mu {\mathrm d}x^\nu = -\mathrm{e}^\nu {\mathrm d}t^2+\mathrm{e}^\lambda {\mathrm d}r^2 + r^2\left({\mathrm d}\theta^2+\sin^2\theta \,{\mathrm d}\varphi^2\right),\label{eq:BGmetric} \end{align}\tag{9}\] where \(\nu=\nu(r)\) and \(\lambda=\lambda(r)\). We describe the matter configuration by a perfect-fluid model. The energy-momentum tensor of the unperturbed fluid is given by \[\begin{align} T_{\mu \nu}=\left(\epsilon+p\right)u_\mu u_\nu +pg_{\mu \nu}\,,\label{eq:BGfluid} \end{align}\tag{10}\] where \(\epsilon\) and \(p\) are the total energy density and pressure, respectively. The unperturbed fluid four-velocity \(u^\mu=(\mathrm{e}^{-\nu/2},0,0,0)\) is normalized by \(g_{\mu \nu} u^\mu u^\nu=-1\).

Within the present setup, Einstein’s equations reduce to the standard set of Tolman–Oppenheimer–Volkov differential equations: \[\begin{align} m'&=4\pi r^2 \epsilon\,,\\ \nu'&=\frac{2m+8\pi r^3 p}{r\left(r-2m\right)}\,,\\ p'&=-\left(\epsilon+p\right)\frac{m+4\pi r^3p}{r\left(r-2m\right)}\,, \end{align}\] where the prime denotes a derivative with respect to \(r\) and \(m(r)=(r/2)(1-\mathrm{e}^{-\lambda})\) is the gravitational mass in a sphere of radius \(r\). The system is closed by imposing the EoS of the fluid composing the star, which we assume to be barotropic, i.e. \(p=p(\epsilon)\). The radius of the star \(r=R_\star\) is defined by \(p(R_\star)=0\). In the exterior \(r>R_\star\), \(\epsilon=p=0\) and the geometry is described by the Schwarzschild metric, i.e. \(m=M\), \(\mathrm{e}^\nu = \mathrm{e}^{-\lambda}=1-2M/r\).

3.2 Dynamical tidal perturbations: Nonperturbative treatment in frequency↩︎

3.2.1 Lindblom-Detweiler equations↩︎

We consider non-radial, time-dependent perturbations of the background 9 –10 with polar parity. Following Lindblowm and Detweiler (LD) [80], [81], we choose the Regge–Wheeler gauge [82] and expand the metric perturbations as: \[\begin{gather} \delta g_{\mu\nu}{\mathrm d}x^\mu {\mathrm d}x^\nu=- r^\ell \left[\mathrm{e}^\nu H^{\ell m}_0 {\mathrm d}t^2-2i \omega r H^{\ell m}_1 {\mathrm d}t{\mathrm d}r\right.\\ \left.+\mathrm{e}^\lambda H_0^{\ell m} {\mathrm d}r^2 +r^2 K^{\ell m}({\mathrm d}\theta^2+\sin^2\theta {\mathrm d}\varphi) \right]Y_{\ell m}\mathrm{e}^{-i\omega t}\,,\label{eq:PerturbedMetric} \end{gather}\tag{11}\] where \(Y_{\ell m}=Y_{\ell m}(\theta,\varphi)\) are the spherical harmonics, \(H_0^{\ell m}=H_0^{\ell m}(r)\), \(H^{\ell m}_1=H^{\ell m}_1(r)\), and \(K^{\ell m}=K^{\ell m}(r)\) are radial functions describing the perturbations, and we leave implicit an integration in \({\mathrm d}\omega\). Note the opposite sign convention for the Fourier transform to that of LD [80], [81].

In the perturbed configuration, the fluid variables are displaced as \(\epsilon\to \epsilon+\delta \epsilon\), \(\rho\to \rho+\delta \rho\), \(u^\mu \to u^\mu+\delta u^\mu\). Here, \(\delta\) denotes the Eulerian perturbation and is related to the Lagrangian perturbation \(\Delta\) through \(\Delta=\delta +{\cal L}_\xi\), where \({\cal L}_\xi\) denotes the Lie derivative along the displacement vector field \(\xi^\mu\). The latter is parametrized by two radial functions \(W^{\ell m}(r)\) and \(V^{\ell m}(r)\) as [80], [81] \[\begin{align} \xi^r& =r^{\ell-1}\mathrm{e}^{-\lambda/2}W^{\ell m}Y_{\ell m}\mathrm{e}^{-i\omega t}\,,\tag{12} \\ \xi_A& =-r^\ell V^{\ell m} {\cal Y}_A^{\ell m} \mathrm{e}^{-i\omega t}\,,\tag{13} \end{align}\] where \({\cal Y}_A^{\ell m}=(\partial_\theta Y_{\ell m},\partial_\varphi Y_{\ell m})\) are the even-parity vector spherical harmonics. Finally, we expand the Lagrangian perturbation of the pressure in terms of a further radial function, \(X^{\ell m}\), as: \[\begin{align} \Delta p=-\mathrm{e}^{\nu/2}r^{-\ell} X^{\ell m}Y_{\ell m} \mathrm{e}^{-i \omega t}\,. \end{align}\] Henceforth, we omit the superscript \(\ell m\) of the perturbation radial variables.

For a static background, the space components of the Eulerian four-velocity perturbations are \(\delta u^i=u^0 \partial_t \xi^i =-i\omega \mathrm{e}^{-\nu/2}\xi^i\) and can thus be expanded in terms of \(W,V\) from Eqs. 12 and 13 . By enforcing the normalization \[(g_{\mu\nu}+\delta g_{\mu\nu})(u^\mu+\delta u^\mu) (u^\nu+\delta u^\nu)=-1\,,\] we find \(\delta u^0=-\frac{1}{2}\mathrm{e}^{-\nu/2}r^\ell H_0Y_{\ell m}\mathrm{e}^{-i\omega t}\).

Then, the conservation laws of the baryon number current and the energy-momentum tensor for the perturbed configuration allow us to express \(\Delta \epsilon\) as \[\begin{gather} \Delta \epsilon=\left(\epsilon+p\right)r^\ell Y_{\ell m} \mathrm{e}^{-i\omega t}\\ \times \left[\frac{H_0}{2}+K-\frac{\ell\left(\ell+1\right)}{r^2}V-\frac{\left(\ell+1\right)W+r W'}{\mathrm{e}^{\lambda/2} r^2}\right]. \end{gather}\] Assuming that the perturbed configuration remains barotropic, \(\Delta p\) is related to \(\Delta \epsilon\) via \[\Delta p =c_s^2 \Delta \epsilon\,,\] where \(c_s:=(\partial p/\partial \epsilon)^{1/2}\) is the adiabatic speed of sound.

Then, the linearized Einstein equations and the energy-momentum conservation reduce to the LD equations [80], [81], a set of four first-order differential equations in the variables \((H_1,K,W,X)\):

\[\begin{align} &\omega^2(\epsilon+p)\mathrm{e}^{-\nu/2}V=X+\frac{p'}{r}\mathrm{e}^{(\nu-\lambda)/2}W-\frac{1}{2}(\epsilon+p)\mathrm{e}^{\nu/2}H_0 \\ &\left(3M+\frac{1}{2}(l+2)(l-1)r+4\pi pr^3\right)H_0=\nonumber\\ &8\pi r^3\mathrm{e}^{-\nu/2}X-\left(\frac{1}{2}l(l+1)(M+4\pi pr^3)-\omega^2r^3\mathrm{e}^{-\nu+\lambda}\right)H_1\nonumber\\ &+\left(\frac{1}{2}(l+2)(l-1)r-\omega^2r^3\mathrm{e}^{-\nu}-\frac{\mathrm{e}^{\lambda}}{r} (M+4\pi pr^3)\right.\times\nonumber\\ &\times (3M-r+4\pi pr^3)\Big)K \end{align}\] and \[\begin{align} H_1'&=\dots\nonumber\\ K'&=\dots\nonumber\\ W'&=re^{(\lambda-\nu)/2}\frac{c_s^{-2}}{(\epsilon+p)}X+\dots\nonumber\\ X'&=\dots \end{align}\] To obtain physical acceptable solutions, we integrate Eqs. 14 –17 imposing regularity at the center of the star. Given the two central amplitudes \((K,W)|_{r=0}\), one can determine the central values of the other variables by [81] \[\begin{align} H_1\big|_{r=0}&=\frac{1}{\ell(\ell+1)}\left[2\ell K+16\pi(\epsilon+p)W\right]\bigg|_{r=0}\,,\\ X\big|_{r=0}&=(\epsilon+p) \mathrm{e}^{\nu/2} \nonumber\\ & \times \left[\left(\frac{4\pi (\epsilon+3p)}{3}-\frac{\omega^2\mathrm{e}^{-\nu}}{\ell}\right)W+\frac{K}{2}\right]\bigg|_{r=0}\,. \end{align}\] Moreover, we impose that the Lagrangian pressure perturbation vanishes at the surface of the star, i.e. \(X(R_\star)=0\). This fixes the ratio \(W/K|_{r=0}\), and we are left with one free constant at the center of the star, which controls the overall normalization of the perturbations and hence does not affect the extraction of the TLNs. We set it to unity without loss of generality, thereby constructing a unique regular perturbed configuration for a given stellar model.

In Sec. 4, we will match the interior solution obtained above with an exterior solution at \(r=R_\star\). To this end, we construct the Zerilli function \(\Psi_\mathrm{Z}\) and its derivative at \(r=R_\star\) through [83] \[\begin{align} \Psi_\mathrm{Z}(r)=&\, \frac{r^{\ell+2}}{n r+3M}\left(K-\mathrm{e}^\nu H_1\right)\,,\nonumber\\ \frac{{\mathrm d}\Psi_\mathrm{Z}}{{\mathrm d}r_*}(r)=&\, \frac{r^\ell}{\left(n r+3M\right)^2}\left[-(nr^2-3M(n r+M))K\right.\label{eq:LogDeriveOfZerilli}\\ &\left.+\mathrm{e}^\nu\left(n (n +1)r^2+3M(n r+2M)\right)H_1\right]\,,\nonumber \end{align}\tag{20}\] where \(r_*=r+2M\log (r/2M-1)\) is the Schwarzschild tortoise coordinate, and \(n\equiv(\ell+2)(\ell-1)/2\). The matching to the exterior solution, outlined in Sec. 4, will be performed by requiring the logarithmic derivative of \(\Psi_\mathrm{Z}\) to be continuous at the stellar surface. To this end, we define the interior logarithmic derivative as \[\begin{align} \left.y^{\rm int}(\omega)\equiv \frac{2M}{\Psi_\mathrm{Z}}\frac{{\mathrm d}\Psi_\mathrm{Z}}{{\mathrm d}r_*}\right|_{r=R_\star}\,,\label{eq:yint} \end{align}\tag{21}\] where the quantities on the right-hand side are computed from Eq. 20 .

We compute \(y^{\rm int}(\omega)\) in the low-frequency regime, \(\omega r_s\ll1\). Although the equations become numerically unstable as \(\omega\to0\), they yield reliable results for \(\omega r_s\gtrsim0.01\), where Eq. 21 is well approximated by a quadratic fit \(y^{\rm int}\simeq y^{\rm int}_0+(\omega r_s)^2y^{\rm int}_2\). The leading-order coefficient agrees with the result for static perturbations (see e.g. Ref. [6]) within percent accuracy.

3.2.2 Alternative formulation↩︎

As an independent analysis, we solve the stellar perturbation problem within a different formulation, equivalent to that of LD. The LD equations are reformulated as three differential equations [84]: \[\begin{align} {\cal H}''& =\alpha_{\cal H, H'}{\cal H}'+\alpha_{\cal H,H}H+\alpha_{{\cal H},W} W+\alpha_{{\cal H},V}V\,,~\tag{22}\\ W'& =\alpha_{W,{\cal H}'}{\cal H}'+\alpha_{W,{\cal H}}{\cal H}+\alpha_{W,W}W+\alpha_{W,V}V\,,\tag{23}\\ V'& =\alpha_{V,{\cal H}}{\cal H}'+\alpha_{V,{\cal H}}{\cal H}+\alpha_{V,W}W+\alpha_{V,V}V\,,\tag{24} \end{align}\] where the radial function \({\cal H}\) is defined through \({\cal H}\equiv -r^\ell H_0\) in Eq. 11 . The coefficients depend on the background quantities and \(\omega\). The explicit forms are provided in Ref. [85]. As in the LD equations, there are two unconstrained parameters \(({\cal H},W)|_{r=0}\). Their ratio is fixed by imposing \(\Delta p=0\) at \(r=R_\star\) and the remaining central amplitude is set to unity without loss of generality. We obtain a regular interior solution for any given stellar model.

To match the interior solution with the exterior at linear order in the perturbation, we compute \(y^{\rm int}\) given in Eq. 21 by using Eq. 20 . The functions \(H_1\) and \(K\) are obtained from \({\cal H}\), \(W\), and \(V\) by using the following relations: \[\begin{align} -i \omega r^{\ell+1}H_1& =\alpha_{H_1,{\cal H}'} {\cal H}'+\alpha_{H_1,{\cal H}} {\cal H}+\alpha_{H_1,V} V+\alpha_{H_1,W}W\,,\tag{25}\\ - r^\ell K& =\alpha_{K,{\cal H}'} {\cal H}'+\alpha_{K,{\cal H}} {\cal H}+\alpha_{K,V} V+\alpha_{K,W} W\,,\tag{26} \end{align}\] where the coefficients can be found in Ref. [85]. We have verified that the values computed in this approach agree with those in the LD approach.

3.3 Dynamical tidal perturbations: perturbative treatment in frequency↩︎

To match the perturbation theory result with the EFT expansion 4 , it is useful to obtain the perturbative solutions order-by-order in \(\omega r_s\ll1\), up to quadratic order.⁶

First, we expand \({\cal H}, V, W\) as \[\begin{align} {\cal H}& ={\cal H}_{0}+\left(\omega r_s\right)^2 {\cal H}_{2}\,,\nonumber\\ V& =V_{0}+\left(\omega r_s\right)^2 V_{2}\,,\label{eq:PerturbativeHVW}\\ W& =W_{0}+\left(\omega r_s\right)^2 W_{2}\,.\nonumber \end{align}\tag{27}\] Then, Eqs. 22 –24 reduce to \[\begin{align} {\cal L}_{\cal H}\left[ {\cal H}_0\right]& =0\,,\tag{28}\\ {\cal L}_V\left[V_0,W_0,{\cal H}_0\right]& =0\,,\tag{29}\\ {\cal L}_W\left[W_0,V_0,{\cal H}_0\right]& =0\tag{30}\,, \end{align}\] at the leading order, and \[\begin{align} {\cal L}_{\cal H}\left[{\cal H}_2\right]&={\cal S}_{\cal H}\,,\tag{31}\\\ {\cal L}_V\left[V_2,W_2, {\cal H}_2\right]&={\cal S}_V\,,\\ {\cal L}_W\left[W_2,V_2,{\cal H}_2\right]&={\cal S}_W\,,\tag{32} \end{align}\] at the quadratic order. Here, we introduce the differential operators \[\begin{gather} {\cal L}_{\cal H}=\frac{{\mathrm d}^2}{{\mathrm d}r^2}+\left[\frac{2}{r}+\mathrm{e}^\lambda\left(\frac{2m}{r^2}+ 4 \pi r \left(p-\epsilon\right) \right)\right]\frac{{\mathrm d}}{{\mathrm d}r}\\ -\mathrm{e}^\lambda \left[ \frac{\ell \left(\ell+1\right)}{r^2}-4 \pi \frac{p+\epsilon}{c_s^2} - 4\pi \left(5 \epsilon+9p\right) \right]-\nu'^2\,, \end{gather}\] and \[\begin{align} {\cal L}_V\left[V_i,W_i,{\cal H}_i\right]=& \, V_i'+\beta_{V,{\cal H'}}{\cal H}_i'+ \beta_{V,{\cal H}}{\cal H}_i \tag{33}\\ &+\beta_{V,V}V_i+\beta_{V,W}W_i\,,\nonumber\\ {\cal L}_W\left[W_i,V_i,{\cal H}_i\right]=& \, W_i'+\beta_{W,{\cal H'}}{\cal H}_i'+ \beta_{W,{\cal H}}{\cal H}_i\tag{34}\\ &+\beta_{W,V}V_i+\beta_{W,W}W_i\,,\nonumber \end{align}\] where \(i=0,2\), and the source terms \[\begin{align} {\cal S}_{\cal H} & = S_{\cal H, H'} {\cal H}_0'+ S_{\cal H, H} {\cal H}_0 + S_{{\cal H},V} V_0 + S_{{\cal H},W} W_0\,,\tag{35}\\ {\cal S}_{V} & = S_{V, \cal H'} {\cal H}_0'+ S_{V, \cal H} {\cal H}_0 + S_{V,V} V_0 + S_{V,W} W_0\,,\tag{36}\\ {\cal S}_W & = S_{W, \cal H'} {H}_0'+ S_{W, \cal H} {\cal H}_0 + S_{W,V} V_0 + S_{W,W} W_0\,,\tag{37} \end{align}\] where the explicit forms of the coefficients are provided in Appendix 8.

We now briefly describe the perturbative construction of the interior solution. At zeroth order in \(\omega\), Eq. 28 reduces to the standard static tidal perturbation equation, coinciding, for \(\ell=2\), with Eq. (15) of Ref. [6]. After obtaining the regular solution for \({\cal H}_0\) at the stellar center, we solve Eqs. 29 and 30 simultaneously to determine the corresponding regular solutions for \(V_0\) and \(W_0\).

At quadratic order in \(\omega\), the reduced equations 31 –32 are sourced by the leading-order solutions. Solving them with the same regularity conditions yields \({\cal H}_2\), \(V_2\), and \(W_2\). As in the nonperturbative analysis, we impose \(\Delta p=0\) at \(r=R_\star\) up to \({\cal O}(\omega^2r_s^2)\), thereby selecting the physical perturbative solution.

Substituting the perturbative solutions for \({\cal H}\), \(V\), and \(W\) into the low-frequency expansions of Eqs. 25 and 26 , we compute \(H_1\) and \(K\) perturbatively and construct the logarithmic derivative of the Zerilli function through Eq. 20 . The resulting perturbative expression for \(y^{\rm int}\) in Eq. 21 agrees with the fully numerical calculation in the low-frequency regime, at zeroth and second order in \(\omega r_s\), within few percent.

Note that the inhomogeneous equations 31 –32 are defined modulo the addition of the homogeneous solution which is regular at the center (with an arbitrary amplitude). However, as expected, \(y^{\rm int}\) is independent of the particular choice of the homogeneous regular solution for \({\cal H}_2\), \(V_2\), and \(W_2\), since this ambiguity only affects the internal decomposition of the regular solution, whose physical content is already fixed (up to normalization) by regularity and the boundary condition \(\Delta p=0\).

4 Dynamical tidal response of a neutron star↩︎

In this section, we first describe the procedure to extract the dynamical tidal response of a neutron star. We then present our main results and discuss the qualitative and quantitative behavior of the dynamical TLN \(c_{\dot{E}}\) computed for several stellar models.

4.1 Extraction of the tidal response↩︎

The tidal perturbation in the exterior is described by the Zerilli function \(\Psi_\mathrm{Z}\), which is obtained by solving the Zerilli equation \[\begin{align} &f(r)\partial_r\left[f(r) \partial_r \Psi_{\mathrm{Z}}\right]\\ &+\left[\omega^2-f(r)\left(\frac{r_s}{r^3}+\frac{2 n}{3 r^2}+\frac{8 n^2(2 n+3)}{3\left(2 n r+3 r_s\right)^2}\right)\right] \Psi_{\mathrm{Z}}=0\,,\nonumber \end{align}\] where \(f(r)= 1-r_s/r\). In the low-frequency expansion \(\omega r_s\ll 1\), the general solution for \(\Psi_\mathrm{Z}\) can be written schematically, up to quadratic order, as \[\begin{align} \Psi_\mathrm{Z}(r)=\Psi_{\mathrm{Z},0}(r) +\left(\omega r_s\right)^2 \Psi_{\mathrm{Z},2}(r)\,. \end{align}\] Following the procedure outlined in Appendix 9, we obtain the general solution for \(\Psi_{\mathrm{Z},0}\) and \(\Psi_{\mathrm{Z},2}\). The large-distance expansion of \(\Psi_\mathrm{Z}\) for the quadrupole mode, \(\ell=2\), then reads:⁷

where \(a_0,b_0,a_2,b_2\) are integration constants. The full solution, including subleading terms in the \(1/r\) expansion, is provided in the Mathematica code available at the GitHub repository [86]. We match this expression to the renormalized Zerilli function \(\Psi_\mathrm{Z}^{\rm R}\) given in Eq. 7 , thereby relating \(a_0,b_0,a_2\), and \(b_2\) to the renormalized integration constants of \(\Psi_\mathrm{Z}^{\rm R}\), namely \(\bar{B}_{\rm reg}\) and \(\bar{B}_{\rm irr}\), via⁸

Thus, once \(a_0\), \(b_0\), \(a_2\), and \(b_2\) are determined by matching to the interior solution, the static and dynamical TLNs, \(c_E\) and \(c_{\dot{E}}\), can be obtained from Eq. 8 , as shown in Appendix 9: \[\begin{align} c_E & = \frac{8G^5M^5}{45 R_\star^5} \frac{b_0}{a_0}\,,\tag{39} \\ c_{\dot{E}} & = \frac{32G^8M^8}{45 R_\star^8} \bigg[\!- \left(1+\frac{107}{105}\frac{b_0}{a_0}\right)\log (\mu r_s) \nonumber \\ &\qquad\qquad\quad\,\, -\frac{67981 }{176400}\frac{ b_0}{ a_0} -\frac{b_0a_2}{a_0^2} \tag{40} \\ &\qquad\qquad\quad\,\, + \frac{945 a_2}{1024 a_0}+\frac{5 b_2}{a_0}-\frac{1015283}{1032192} \bigg].\nonumber \end{align}\] Before proceeding, it is worth commenting on the structure of the result 40 . Note that the coefficient of \(\log(\mu r_s)\) in \(c_{\dot{E}}\) contains two types of contributions. One is a universal term, independent of the integration constants, which arises solely from the gravitational nonlinearities of the Einstein–Hilbert action and is therefore common to any object, in particular black holes (see, e.g., Refs. [36], [44], [50], [51], [55]). The other contribution is \(107 b_0/(105 a_0)\), which is proportional to the static Love number \(c_E\) and is thus unambiguously fixed by the ratio \(b_0/a_0\). This should be contrasted with the remaining (scheme-dependent) finite terms on the second and third lines of Eq. 40 , which depend explicitly on the free parameters \(a_2\) and \(b_2\), and therefore require a full matching at second order in \(\omega\) to be determined.

The exterior solution for \(\Psi_\mathrm{Z}\) is then matched to the interior solution, obtained in the previous section, by requiring continuity of the logarithmic derivative of the Zerilli function at the stellar surface. For the exterior solution, we introduce \[\begin{align} y^{\rm ext}(\omega) \equiv &\left.\frac{2M}{\Psi_\mathrm{Z}}\frac{{\mathrm d}\Psi_\mathrm{Z}}{{\mathrm d}r_*}\right|_{r=R_\star}.\label{eq:yext} \end{align}\tag{41}\] By equating the exterior logarithmic derivative \(y^{\rm ext}\) to the interior one, \(y^{\rm int}\), given in Eq. 21 , we match the interior and exterior solutions for \(\Psi_Z\) to second order in \(\omega r_s\). This determines the ratios between \(b_0\), \(a_2\), and \(b_2\) with respect to \(a_0\) for a given stellar model, from which the TLNs \(c_E\) and \(c_{\dot{E}}\) are obtained.

4.2 Dynamical Love numbers↩︎

In what follows, we present several properties of the dynamical tidal response of neutron stars. We emphasize that the TLNs computed through Eqs. 39 and 40 are defined within the point-particle EFT framework and admit a direct connection to the gravitational waveform discussed in Sec. 5.1. For the reader’s convenience, in Appendix 10 we summarize the connections with the tidal constants introduced in alternative approaches based on the PN framework [26], [27], [54].

First of all, we construct several stellar configurations using realistic, tabulated EoS models released online [87]. Specifically, we consider SLy4 [88], AP4 [89], FPS [90], MS1 [91], and ENG [92]. Figure 1 presents the mass-radius relations for the stellar models studied in this work.

The results of the matching are shown in Fig. 2. The upper-left and lower-left panels display the static (\(c_E\)) and dynamical (\(c_{\dot{E}}\)) TLNs, respectively, as functions of the neutron-star compactness for various EoS. Accounting for the different conventions adopted in the literature—namely \(c_E=k_0/3\), where \(k_0\) is the Love number defined in Refs. [6], [9]—the upper panel reproduces the results for the static TLN reported in previous works, see e.g. [10], [93].

The bottom panel shows that, for the EoS considered in this work, \(c_{\dot{E}}\) exhibits trends qualitatively similar to those of \(c_E\): (i) it reaches a maximum around \(GM/R_\star \simeq 0.12\); (ii) it decreases toward both the low- and high-compactness regimes; and (iii) stiffer EoS generally yield larger peak values. A distinctive feature of \(c_{\dot{E}}\) is its explicit scale dependence through the \(\log (\mu r_s)\) term in Eq. 40 . We find that this dependence is relatively mild at low compactness, whereas smaller values of \(\mu\) enhance the peak of \(c_{\dot{E}}\) and produce a slower decay in the high-compactness regime. This is consistent with the fact that the running term is a relativistic effect.

It is worth noting that \(\mu\) is associated with the characteristic length scale of the binary system and is therefore expected to scale as \(\mu\sim 1/r\), where \(r\) is the orbital separation. This implies that \(\mu r_s\) is much smaller than unity at the early stage of the late-inspiral phase and grows to \({\cal O}(10^{-1})\) as the binary evolves. Equation 40 then indicates that, because the logarithmic term enters with a minus sign, \(c_{\dot{E}}\) may decrease during the binary evolution. In this respect, the behavior is reminiscent of an asymptotically free coupling: as the binary shrinks, the renormalization scale \(\mu \sim 1/r\) grows and the dynamical Love number coupling decreases logarithmically, much like the running of the strong coupling constant in quantum chromodynamics toward higher energies. In both cases, the running is driven by the self-interactions of the mediator: in quantum chromodynamics, gluons carry color charge and interact among themselves, antishielding the source and yielding a negative beta function; in gravity, gravitons carry energy-momentum and self-interact through the nonlinearities of general relativity, generating the analogous logarithmic flow of the tidal coupling. However, the decrease in \(c_{\dot{E}}\) does not necessarily imply that the tidal response is suppressed in the late inspiral. Indeed, Fig. 3 shows that the combination \(c_E+\omega^2 R_\star^3/(G M) c_{\dot{E}}\) of Eq. 8 increases as the orbital length scale decreases. This behavior results from the competition between the power-law growth in \(\omega^2 (r)\) and the logarithmic running. In the range shown, the former dominates over the latter.

In the right panels of Fig. 2, we show the dimensionless tidal coefficients \(\Lambda_0\) and \(\Lambda_2\) as functions of the stellar compactness, \(C \equiv GM/R_\star\), defined by \[\Lambda_0 = 2\left(\frac{GM}{R_\star}\right)^{-5} c_E, \qquad \Lambda_2 =2 \left(\frac{GM}{R_\star}\right)^{-8} c_{\dot{E}}\,. \label{Lambdas}\tag{42}\] The top panel reproduces the well-known behavior of the static tidal deformability; see, e.g., Ref. [94]. The bottom panel shows that \(\Lambda_2\) is typically about three orders of magnitude larger than \(\Lambda_0\). This difference is primarily driven by the stronger compactness dependence of \(\Lambda_2\), which scales as \(C^{-8}\) compared to the \(C^{-5}\) scaling of \(\Lambda_0\). As we discuss in Sec. 5, \(\Lambda_0\) and \(\Lambda_2\) are precisely the combinations that enter the gravitational waveform and therefore provide the relevant parameters for describing dynamical tidal effects in GW observations.

4.3 Approximate universal relations↩︎

Neutron-star observables are known to exhibit approximately EoS-insensitive relations [94]–[96]. It is therefore reasonable to expect that an approximate universal relation holds between \(\Lambda_0\) and \(\Lambda_2\). Indeed, such a relation was found in Ref. [31] for the dynamical Love numbers defined in that work. Here, we establish the approximate universality for our TLNs, which are defined unambiguously and directly connected to the gravitational waveform. We stress that this approximate universality becomes manifest only when using the waveform couplings \(\Lambda_0\) and \(\Lambda_2\), while it remains hidden for other combinations of the tidal parameters [97]. Figure 4 shows that an approximate universal relation emerges only for fixed values of \(\mu r_s\). We fit the numerical data using the following fitting formulae: \[\begin{align} Y|_{\mu r_s=0.01}&= 3.9171 +1.0707 X+0.023100 X^2\,,\tag{43}\\ Y|_{\mu r_s=1}&= 3.3840 +1.1621 X+0.019180 X^2\,,\tag{44} \end{align}\] where \((X,Y)=(\log \Lambda_0, \log \Lambda_2)\). The typical relative errors between the data and these empirical fits are smaller than \(10\%\), consistent with the results of Ref. [31]. Unlike the standard scale-independent I-Love-Q relations [94], [95], the relation between \(\Lambda_0\) and \(\Lambda_2\) inherits the explicit scale dependence of the dynamical Love number through \(\mu\), so that the fit coefficients in Eqs. 43 and 44 change with \(\mu r_s\).

As we shall discuss in Sec. 5, the running of the dynamical TLNs has only a subdominant impact on GW observations, with the overall tidal effect being primarily controlled by the scale-independent component. For this reason, it will be convenient to adopt the choice \(\mu r_s=1\), which allows using the approximate universal relation in Eq. 44 .

4.4 Comparison with mode-sum representation↩︎

As in the Newtonian case [15]–[24], [26], one may expect the dynamical tidal response to be closely related to the characteristic oscillation spectrum of the star, namely its quasi-normal modes. These are the complex eigenfreqencies for which the homogeneous perturbation equations discussed in Sec. 3.2 admit solutions that are regular at the stellar center and purely outgoing at infinity [98].

For neutron stars, the oscillation spectrum is typically dominated by the fundamental fluid (\(f\)-) mode, whose frequency is roughly given by [98] \[\omega_f \simeq c_f\sqrt{\frac{GM}{R_\star^3}} \approx \mathcal{O}(1)\,\mathrm{kHz}\,,\] where \(c_f\) is an \(\mathcal{O}(1)\) coefficient that depends on the EoS (see, for example, Ref. [99] for an explicit evaluation for polytropic EoS). At the Newtonian level, the dynamical tidal deformability admits a mode-sum representation in which the \(f\)-mode contribution plays a dominant role because of its relatively low frequency [15], [17], [19], [26]. Consistent with this picture, \(\omega_f\) typically agrees with an effective frequency \(\omega_*\) within \(0.1\%\), where \(\omega_*\) is defined by [26]⁹ \[\begin{align} \omega_* \equiv \sqrt{\frac{GM}{R_\star^3}\frac{c_E}{c_{\dot{E}}}}\,.\label{eq:effectiveomega} \end{align}\tag{45}\] Therefore, in the Newtonian regime the dynamical TLN can be very well approximated in terms of the \(f\)-mode frequency and the static TLN, \[\begin{align} c_{\dot{E}}\approx\frac{GM}{R_\star}\frac{c_{E}}{\omega_f^2}\,. \end{align}\] In the relativistic regime, however, the relation between tidal response and quasi-normal modes becomes more subtle. In particular, a mode-based description may not be complete, and additional contributions can cause \(\omega_*\) to deviate from \(\omega_f\) [26], [31], [32], and hence the dynamical TLN cannot be simply obtained from the quasi-normal modes.

To quantify this discrepancy, in Fig. 5 we show the ratio \(\omega_*/\omega_f\) as a function of compactness for several values of \(\mu\). For sufficiently small values of \(\mu\), the deviation from unity decreases with compactness, systematically lowering the ratio \(\omega_*/\omega_f\). Remarkably, throughout the most interesting interval \(0.1\lesssim \mu r_s\lesssim 1\), the ratio remains within approximately \(5\%\) of unity over the entire compactness range considered. This suggests that, despite relativistic corrections, the dynamical tidal response of realistic neutron stars during the late inspiral is still largely controlled by the \(f\)-mode frequency. Furthermore, we stress that the ratio is compatible with the corresponding prediction obtained from the tidal constants introduced in Refs. [26], [54], once the convention difference in Eq. 86 is accounted for (see Appendix 10 for their definition and a comparison with our dynamical TLN, as well as with other approaches [27], [31]).

Our findings are qualitatively consistent with those of Ref. [32], while providing an independent validation of this behavior within the EFT framework developed here, including the effects of renormalization-group running.

5 Measurability with GW detections↩︎

Having computed the dynamical TLNs of neutron stars, we now assess their impact on the GW signal emitted during the inspiral of a compact binary and their measurability with current and future detectors. In Sec. 5.1 we derive the PN waveform phase, including both the static and dynamical tidal contributions and the distinctive logarithmic running of the latter, following Ref. [51]. As shown in that work, the dynamical Love number enters the waveform at the 8PN order. In Sec. 5.2 we quantify the observability of these effects through a Fisher-matrix analysis and a mismatch calculation, considering the LVK O4 network and the ET at its nominal sensitivity. Together, these two analyses provide complementary views of the same physics: the Fisher matrix quantifies the statistical precision with which the dynamical Love number can be extracted, while the mismatch measures the systematic error incurred by neglecting it.

5.1 PN waveform derivation↩︎

In this section, we use the EFT Lagrangian 4 , with the renormalized induced dynamical tidal response couplings obtained in the previous sections, to study a binary inspiral. Although the renormalized couplings were derived for an isolated neutron star, no additional divergences are expected when the same Lagrangian is applied to the two-body problem. As in any EFT, once all counterterms consistent with the symmetries have been included up to a given order in the derivative expansion, the divergences are removed once and for all. The theory becomes fully predictive at that order and can be used to compute any process consistently, with finite results. In particular, the renormalization scale \(\mu\) drops out of physical observables—such as cross sections, scattering angles, or waveform phases—which depend only on the renormalized couplings evaluated at a chosen scale together with the corresponding logarithmic dependence on the relevant physical scale (e.g., energy or distance).

The procedure for deriving the interaction Hamiltonian, including the counterterms, is discussed in detail in e.g. Ref. [39].¹⁰ For the purposes of this work, however, we follow a simpler approach [51]. Since the renormalized dynamical Love number \(c_{\dot{E}}\) depends on the renormalization scale \(\mu\), and the natural scale in the binary problem is set by the inverse orbital separation, we identify \(\mu \sim \gamma/r\), with \(\gamma\) some constant number. Because our goal is to estimate the potential observability of dynamical Love number effects, the precise relation between \(\mu\) and \(r\) will not be important. As the inspiral proceeds and the orbital separation decreases, the tidal response evolves logarithmically, as illustrated for example by Eq. 38 . Importantly, this running cannot in general be absorbed into a redefinition of the scale-independent part of \(c_{\dot{E}}\). It therefore leaves a characteristic imprint on the GW phase in the form of a frequency-dependent 8PN correction, qualitatively distinct from the standard frequency-independent tidal contributions. In practice, however, this effect is suppressed in the small-compactness regime and, as we shall see explicitly below, does not significantly affect our estimates.

To connect with the waveform calculation, recall that the dynamical couplings \(c_E\) and \(c_{\dot{E}}\) enter the effective two-body Lagrangian through Eq. 4 . Their effects on the GW phase are conveniently encoded in the dimensionless tidal parameters \(\Lambda_0\) and \(\Lambda_2\), defined in Eq. 42 . It is therefore useful to express the renormalization-group running directly in terms of these quantities. Following Ref. [51], we write the running of the dynamical tidal parameter as \[\Lambda_2(\mu) \;\to\; \Lambda_2(r) \equiv \bar{\Lambda}_2 - \beta_2\,\log\left(\frac{\gamma\,r_s}{r}\right), \label{eq:Lambda295running}\tag{46}\] where \(\bar{\Lambda}_2\) captures the value of the dynamical TLN at the UV matching scale \(r_s\), given by the neutron-star Schwarzschild radius, while \(\beta_2\) encodes the logarithmic running associated with the renormalization procedure. A similar approach has been followed in Ref. [39], where an effective radial-dependent Love number has been included in the binary Lagrangian. The relation between \((\bar{\Lambda}_2,\,\beta_2)\) and the EFT coupling \(c_{\dot{E}}\) of Sec. 2 is given by \[\begin{align} \bar{\Lambda}_2 &= 2\left(\frac{GM}{R_\star}\right)^{-8} c_{\dot{E}} \, (\mu = r_s^{-1})\,, \nonumber \\ \quad \beta_2 &= 2\left(\frac{GM}{R_\star}\right)^{-8} c_{\dot{E}}^{(\log)} = - \frac{64}{45} \left(1 + \frac{321}{224} \Lambda_0 \right)\,, \label{eq:Lambda295dict} \end{align}\tag{47}\] where \(c_{\dot{E}}^{(\log)}\) is the coefficient of \(\log(\mu r_s)\) in Eq. 40 , while the static counterpart is \(\Lambda_0 = 2(GM/R_\star)^{-5}c_E\) as in Eq. 42 . The key observation is that, since \(\Lambda_2\propto C^{-8}\) while \(\Lambda_0\propto C^{-5}\), the dynamical parameter is enhanced by three additional inverse powers of the compactness \(C = GM/R_\star\) relative to the static one, an amplification that, as shown in Fig. 2, reaches several orders of magnitude at typical neutron-star compactnesses and is the principal reason why the 8PN dynamical tide may be observable despite its formally high PN order (see Ref. [56] for a similar enhancement for the quadratic TLNs).

We now derive the GW phase following Ref. [51]. Consider a binary of two non-rotating neutron stars with masses \(m_1\), \(m_2\), total mass \(M_\text{\tiny tot} = m_1+m_2\), symmetric mass ratio \(\eta = m_1 m_2/M_\text{\tiny tot}^2\), and orbital angular frequency \(\omega\). The effective binary Lagrangian for circular orbits, including the leading tidal interactions of both bodies described in Eq. 4 , reads \[\begin{align} L &= \frac{1}{2}\eta M_\text{\tiny tot} \dot{r}^2 + \frac{1}{2}\eta M_\text{\tiny tot} r^2\dot{\varphi}^2 + \frac{G \eta M_\text{\tiny tot}^2}{r} \nonumber \\ &- \frac{m_2 (Gm_2)^4}{4} \left[\Lambda^{(2)}_0\,E_{ij}E^{ij} + (G m_2)^2\,\Lambda^{(2)}_2(r)\,\dot{E}_{ij}\dot{E}^{ij}\right] \nonumber \\ & + (1 \leftrightarrow 2)\,, \label{eq:Lagrangian} \end{align}\tag{48}\] where \({E}_{ij}\) describes the tidal field sourced by the companion neutron star [11]. Time derivatives are denoted by a dot, and the superscript \((i)\) labels the tidal response of body \(i = 1,2\). Evaluating the Euler-Lagrange equations on circular orbits, one obtains the expression for the binding energy [51] \[\begin{align} \mathcal{E}(\omega) & = -\frac{1}{2}\eta M_\text{\tiny tot} (G M_\text{\tiny tot} \omega)^{2/3} \bigg\{1 + \frac{9\,(G m_1 \omega)^4}{(GM_\text{\tiny tot} \omega)^{2/3}} \frac{m_2}{M_\text{\tiny tot}} \Lambda^{(1)}_0 \nonumber \\ & + \frac{9\, (G m_1 \omega)^6}{(GM_\text{\tiny tot} \omega)^{2/3}} \frac{m_2}{M_\text{\tiny tot}} \left[ 5 \bar{\Lambda}^{(1)}_2 - \frac{2}{3} \beta^{(1)}_2 \right. \nonumber \\ & \left. + 5 \beta^{(1)}_2\log\left(\frac{(GM_\text{\tiny tot} \omega)^{1/3}}{\omega \, \gamma \, r_s^{(1)}} \right) \right] + (1\!\leftrightarrow\!2) \bigg\}\,, \label{eq:Ebind} \end{align}\tag{49}\] and the GW flux \[\begin{align} \mathcal{F}(\omega) &= \frac{32}{5} \eta^2 (GM_\text{\tiny tot} \omega)^{10/3} \bigg\{1 \nonumber \\ & - \frac{ 6 \, (Gm_1 \omega)^4 }{(GM_\text{\tiny tot} \omega)^{2/3}} \left(2\frac{m_2}{M_\text{\tiny tot}}+1\right)\Lambda^{(1)}_0 \nonumber \\ & + \frac{6 \, (Gm_1 \omega)^6}{(GM_\text{\tiny tot} \omega)^{2/3}} \left[ \frac{m_2}{M_\text{\tiny tot}}\beta^{(1)}_2 - 2\left(3\frac{m_2}{M_\text{\tiny tot}}+2\right) \right. \nonumber \\ & \left. \times \left(\bar{\Lambda}^{(1)}_2 + \beta^{(1)}_2 \log\left(\frac{ (GM_\text{\tiny tot} \omega)^{1/3}}{\omega \, \gamma \, r_s^{(1)}}\right) \right) \right] + (1\!\leftrightarrow\!2) \bigg\}\,. \label{eq:flux} \end{align}\tag{50}\] Note that, in deriving these expressions, the orbital separation \(r\) is eliminated in favor of the orbital frequency \(\omega\) via the condition for circular orbits, \(\partial L/\partial r = 0\); crucially, the resulting relation \(r(\omega)\) receives leading-order dynamical tidal corrections [51], which must be consistently included when expanding the energy and flux in powers of the frequency.

In the stationary-phase approximation, \({\mathrm d}^2\psi/{\mathrm d}\omega^2 = -(2/\mathcal{F})\,{\mathrm d}\mathcal{E}/{\mathrm d}\omega\) [11], and integrating twice in frequency yields the waveform phase as a function of the dimensionless PN parameter \(x = (GM_\text{\tiny tot}\omega)^{2/3}\) [51]¹¹ \[\begin{align} \psi_\text{\tiny tidal}(x) & = \frac{3}{128\,\eta\,x^{5/2}} \left[ 1 + \frac{39}{2}\,\Lambda_\text{\tiny 5PN}\,x^5 + \frac{15}{11} \Lambda_\text{\tiny 8PN}\,x^8 \right. \nonumber \\ & \left. + \frac{15}{11}\sum_i B^{(i)}_\text{\tiny 8PN}\,x^8\log\left(\frac{ GM_\text{\tiny tot}}{x\, \gamma \, r_s^{(i)}} \right) \right]\,, \label{eq:phase} \end{align}\tag{51}\] where we have defined the combinations \[\begin{align} \Lambda_\text{\tiny 5PN} &= \frac{16}{13} \sum_i \frac{m_i^3(m_i+11\eta M_\text{\tiny tot})}{M_\text{\tiny tot}^4} \Lambda^{(i)}_0\,, \nonumber \\ \Lambda_\text{\tiny 8PN} &= \sum_i \frac{m_i^5}{M_\text{\tiny tot}^6} \left[(8m_i+147\eta M_\text{\tiny tot})\bar{\Lambda}^{(i)}_2 \right. \nonumber \\ & \left. \quad +\left(\frac{38}{11}m_i+\frac{1253}{44}\eta M_\text{\tiny tot}\right)\beta^{(i)}_2\right]\,, \nonumber \\ B^{(i)}_\text{\tiny 8PN} &= \frac{m_i^5(8m_i+147\eta M_\text{\tiny tot})}{M_\text{\tiny tot}^6}\,\beta^{(i)}_2\,. \label{eq:B295def} \end{align}\tag{52}\] Equation 51 is the central result of this section. It contains three qualitatively distinct contributions. The first is the well-known 5PN static tidal term [6], [7], proportional to \(\Lambda_0\). The second and third arise from the dynamical tidal coupling \(c_{\dot{E}}\) computed in Sec. 4, and enter the waveform at 8PN order. The term proportional to \(\Lambda_\text{\tiny 8PN}\) captures the finite, EoS-dependent contribution of the dynamical TLN, while the logarithmic terms \(B^{(i)}_\text{\tiny 8PN}\) encode its renormalization-group running. Crucially, because \(\bar{\Lambda}_2 \propto C^{-8}\), the 8PN dynamical tide is enhanced by three additional inverse powers of the compactness compared to the 5PN static term \(\Lambda_0 \propto C^{-5}\). This makes the dynamical tide potentially observable despite its formally high PN order and also makes it much larger than (currently unknown) 8PN point-particle terms [56]. Furthermore, the logarithmic term in Eq. 51 provides a qualitative handle to, in principle, disentangle the dynamical contribution from higher-order point-particle corrections entering at the same PN order, since the latter do not carry the same frequency dependence. However, as we will comment below, this logarithmic running appears to be hardly measurable even with next-generation GW detectors: it enters at subleading order relative to the bare \(\Lambda_\text{\tiny 8PN}\) term, and the numerical coefficient \(\gamma\) that relates the renormalization scale to the orbital separation does not significantly affect the waveform. The dominant and observationally relevant quantity is therefore the EoS-dependent constant \(\Lambda_\text{\tiny 8PN}\).

This hierarchy is manifest when looking at Fig. 6, which shows the absolute contributions to the tidal GW phase \(\psi_\text{\tiny tidal}\) from the static (solid), finite dynamical (dashed and dot-dashed), and logarithmic (dotted) TLNs, computed from Eqs. 51 and 52 , as a function of the dimensionless binary frequency \(x\), taken to span the frequency range between the ET minimum frequency and the neutron star contact frequency (see Eq. 57 below). The hierarchy among these contributions is immediately apparent: the static Love numbers, entering at 5PN order, overwhelmingly dominate the tidal phase throughout the inspiral. At 8PN order, the finite dynamical contribution \(\bar{\Lambda}_2\) exceeds the logarithmic term \(\beta_2\) associated with the renormalization-group running of the Love numbers. Despite the logarithmic coefficient being proportional to the static Love number, see Eq. 47 , it is suppressed relative to the finite term due to the additional factors of compactness \(C^{-3}\). This clear separation of scales motivates our choice, in Sec. 5.2, to focus the data-analysis study on the bare 8PN coefficient \(\Lambda_\text{\tiny 8PN}\). Furthermore, as discussed in Sec. 4.3, this hierarchy suggests that the approximate universal relation between the 8PN and 5PN tidal parameters (neglecting the running) could be leveraged to reduce the waveform parameter space, potentially tightening constraints on the latter; we leave a detailed exploration of this avenue to future work.

5.2 Data analysis↩︎

We now assess the observability of the 8PN dynamical tidal corrections of Eq. 51 , following the approach of Ref. [51]. We consider two complementary diagnostics: a Fisher-matrix analysis, which quantifies statistical measurement precision, and a mismatch calculation, which directly probes the systematic bias incurred when the dynamical tide is omitted from the template.

We model the GW signal using the IMRPhenomD waveform family [100], [101], which captures the full inspiral evolution. The tidal phase of Eq. 51 is appended to the point-particle phase, including the 5PN static term \(\Lambda_\text{\tiny 5PN}\) and the 8PN dynamical term \(\Lambda_\text{\tiny 8PN}\). In the frequency domain the template reads \[\tilde{h}(f;\,\boldsymbol{\theta}) = \mathcal{C}_\Omega\,\mathcal{A}(f;\,\boldsymbol{\theta})\, \mathrm{e}^{\,i\left[\psi_\text{\tiny PP}(f;\,\boldsymbol{\theta}) +\psi_\text{\tiny tidal}(f;\,\boldsymbol{\theta})\right]}\,, \label{eq:waveform}\tag{53}\] where \(\psi_\text{\tiny PP}\) describes the standard point-particle phase, while the leading-order amplitude is \[\mathcal{A}(f;\,\boldsymbol{\theta}) = \sqrt{\frac{5}{24}}\, \frac{G^{5/6} \mathcal{M}_c^{5/6}}{\pi^{2/3}\,d_\text{\tiny L}\,f^{7/6}}\,, \label{eq:amplitude}\tag{54}\] in terms of the chirp mass \(\mathcal{M}_c = \eta^{3/5}M_\text{\tiny tot}\) and luminosity distance \(d_\text{\tiny L}\). The geometric factor \(\mathcal{C}_\Omega\) captures the detector’s response and depends on the binary’s inclination angle and antenna pattern functions (which in turn vary with the source’s sky location and polarization angle). For simplicity, we assume optimally oriented binaries, thereby neglecting these dependencies in our analysis. We consider equal-mass (\(m_1 = m_2 = m_\text{\tiny NS}\)), non-spinning binary neutron-star systems at \(d_\text{\tiny L} = 100\,\text{Mpc}\) and three representative EoS: SLy4 [88], AP4 [89], and MS1 [91].

To forecast measurement uncertainties, we use the Fisher information matrix formalism [102], which is reliable in the limit of large signal-to-noise ratio (SNR). The Fisher matrix is \[\Gamma_{ij} = \left\langle\frac{\partial h}{\partial\theta^i} \,\bigg|\, \frac{\partial h}{\partial\theta^j}\right\rangle \bigg|_{\boldsymbol{\theta}=\hat{\boldsymbol{\theta}}}\,, \label{eq:fisher}\tag{55}\] where the noise-weighted inner product is \[\langle h_1 | h_2 \rangle = 4\,\operatorname{Re}\int_{f_\text{\tiny min}}^{f_\text{\tiny max}} \frac{\tilde{h}_1^*(f)\,\tilde{h}_2(f)}{S_n(f)}\,{\rm d}f\,, \label{eq:inner95product}\tag{56}\] with \(S_n(f)\) the one-sided noise power spectral density. The \(1\sigma\) uncertainty on the \(i\)-th parameter is \(\sigma_i = \sqrt{(\Gamma^{-1})_{ii}}\) and the SNR is \(\rho = \sqrt{\langle h|h\rangle}\). The full parameter vector is \(\boldsymbol{\theta} = \{\mathcal{M}_c,\,\eta,\,\chi_s,\,\chi_a,\,t_c,\,\phi_c,\,\Lambda_\text{\tiny 5PN},\,\Lambda_\text{\tiny 8PN}\}\) with injected coalescence time and phase \((t_c,\phi_c)=(0,0)\), and vanishing injected spin terms, \(\chi_s = \chi_a = 0\). We consider two detector configurations: the LVK network at O4 sensitivity with \((f_\text{\tiny min},f_\text{\tiny max}) = (10, 4 \cdot 10^3)\,{\rm Hz}\), and the ET in its ET-D configuration [60], [103] with \((f_\text{\tiny min},f_\text{\tiny max}) = (2, 10^4)\,{\rm Hz}\). The inspiral sequence proceeds until the neutron stars come in contact or get tidally disrupted by their companion, which happens around the Roche frequency [104] \[\label{eq:Roche} f_\text{\tiny Roche} = 12 \sqrt{6} \, C^{3/2} f_\text{\tiny ISCO}\,,\tag{57}\] in terms of their compactness and frequency at the innermost circular orbit (ISCO) \[f_\text{\tiny ISCO} = 4.4 \, {\rm kHz} \left(\frac{M_\odot}{M_\text{\tiny tot}} \right) \,.\] As a complementary and more conservative diagnostic, we compute the mismatch between a signal template that includes the dynamical tide and one that neglects it: \[\mathcal{M} = 1 - \underset{\{\phi_c, t_c\}}{\rm max} \frac{\langle h_1 | h_2 \rangle}{\sqrt{\langle h_1|h_1\rangle\,\langle h_2|h_2\rangle}}\,, \label{eq:mismatch}\tag{58}\] where \(h_1\) includes \(\Lambda_\text{\tiny 8PN}\neq 0\) and \(h_2\) is computed with \(\Lambda_\text{\tiny 8PN}=0\). A mismatch \(\mathcal{M} \lesssim 1/(2\rho^2)\) guarantees that two waveforms are equivalent within statistical errors at a given SNR, while exceeding this limit might indicate a systematic bias in parameter estimation [105].

Figure 7 shows the relative \(1\sigma\) uncertainties on \(\Lambda_\text{\tiny 5PN}\) (left panel) and \(\Lambda_\text{\tiny 8PN}\) (right panel) as functions of the individual neutron star mass \(m_\text{\tiny NS}\), for the three EoS at \(d_\text{\tiny L}=100\,\text{Mpc}\); the horizontal dashed line marks a relative error of 100%. For the static Love number, ET (solid lines) lies well below this threshold across essentially the entire mass range, reaching \(\sigma_{\Lambda_\text{\tiny 5PN}}/\Lambda_\text{\tiny 5PN} \lesssim 10\%\) for \(m_{\text{\tiny NS}} \lesssim 1.4\,M_\odot\), and approaching the 100% line only for the most massive, softest configurations. By contrast, the LVK O4 curves (dashed lines) sit close to or above the 100% line, so that \(\Lambda_\text{\tiny 5PN}\) is at best marginally constrained at current sensitivity.

Two features of the curves are noteworthy. First, the relative error grows steeply with mass. This happens because more massive neutron stars are more compact, so that both tidal deformabilities fall rapidly with compactness, shrinking the tidal signal; in addition, the cutoff of the inspiral scales inversely with the total mass through \(f_\text{\tiny ISCO} \propto M_\text{\tiny tot}^{-1}\), so that the contact/Roche frequency of Eq. 57 moves to lower values and fewer wave cycles accumulate in the tidal-sensitive band. Both effects reduce the measurable tidal imprint and inflate the statistical uncertainty. Second, at fixed mass the relative error is smallest for the stiffest EoS (MS1) and largest for the softer ones (AP4, SLy). Since the tidal terms enter the phase linearly, the absolute Fisher uncertainty \(\sigma_{\Lambda_\text{\tiny 5PN}}\) is set essentially by the detector, the masses, and the available frequency band, and depends only weakly on the EoS; the relative error \(\sigma_{\Lambda_\text{\tiny 5PN}}/\Lambda_\text{\tiny 5PN}\) is therefore controlled primarily by the magnitude of \(\Lambda_\text{\tiny 5PN}\) itself. Stiffer EoS produce larger radii and hence lower compactness (see Fig. 1), which through the steep inverse-compactness scaling translates into substantially larger \(\Lambda_0\) and \(\Lambda_2\), and thus a stronger and more easily measured tidal imprint. The lower compactness of stiffer stars also lowers their contact/Roche frequency, an effect that by itself would reduce the in-band tidal information.

For the dynamical Love number the picture is more nuanced: ET reaches \(\sigma_{\Lambda_\text{\tiny 8PN}}/\Lambda_\text{\tiny 8PN} \sim 10\%\) only for \(m_{\text{\tiny NS}} \lesssim 1.0\,M_\odot\) and the stiffest EoS (MS1), the only case that drops appreciably below the 100% line; for the softer EoS (SLy4, AP4) the relative error stays above \(100\%\) for masses \(m_{\text{\tiny NS}} \gtrsim 1.2\,M_\odot\), where \(\Lambda_\text{\tiny 8PN}\) is not individually measurable on its own. LVK O4 yields uncertainties orders of magnitude above the 100 % for all configurations. Clearly, the presented relative errors have a simple linear scaling with the source distance.

The key factor enabling ET to approach the dynamical tide is the large enhancement of \(\Lambda_2\) relative to \(\Lambda_0\) visible in Fig. 2: since \(\Lambda_2 \propto C^{-8}\) while \(\Lambda_0 \propto C^{-5}\), the dynamical parameter is larger by roughly three orders of magnitude at typical neutron-star compactnesses, compensating the formal \(x^3\) suppression of the 8PN term relative to the 5PN one. The strong EoS dependence of \(\Lambda_\text{\tiny 8PN}\) implies that a detection would simultaneously constrain the stiffness of dense matter at supranuclear densities. This situation is qualitatively different from the black-hole case studied in Ref. [51], where the dynamical tidal correction is unobservable even with future-generation detectors because the analogous \((r_s/M)^8\) amplification factor reduces to \(\mathcal{O}(1)\). We also note that our Fisher analysis treats \(\Lambda_\text{\tiny 8PN}\) as the primary free parameter; the logarithmic correction \(B^{(i)}_\text{\tiny 8PN}\) is not independently varied, as it is numerically subleading relative to \(\Lambda_\text{\tiny 8PN}\) for the orbital frequencies relevant to the late inspiral, and the sensitivity to the free coefficient \(\gamma\) is correspondingly negligible. Finally, employing the approximate universality relation between \(\Lambda_0\) and \(\Lambda_2\) established above would allow one to express \(\Lambda_\text{\tiny 8PN}\) in terms of \(\Lambda_\text{\tiny 5PN}\), thereby reducing the dimensionality of the parameter space and improving the measurability of the static coefficient.

Figure 8 shows the mismatch \(\mathcal{M}\) as a function of \(m_{\text{\tiny NS}}\) for the same configurations, with the horizontal dashed lines marking the systematic-bias threshold \(1/(2\rho^2)\) for LVK O4 (upper line) and ET (lower line). A curve lying above its detector threshold signals that the dephasing induced by omitting the dynamical tide is larger than the statistical error, so that neglecting \(\Lambda_\text{\tiny 8PN}\) may bias the recovery of the remaining parameters. Because the dynamical tide is strongest for less compact stars, the mismatch is largest at low masses and for the stiffer EoS, and decreases steeply as \(m_\text{\tiny NS}\) grows.

For ET and a binary neutron star at \(d_\text{\tiny L}=100\,\text{Mpc}\), the mismatch exceeds the threshold for \(m_{\text{\tiny NS}} \lesssim 1.2\,M_\odot\) and the softest EoS (AP4, SLy), while being always above it for the stiffer EoS (MS1), whereas for LVK O4 it stays below the (much higher) O4 threshold for all configurations, in agreement with the Fisher-matrix result that the dynamical tide is inaccessible at current sensitivity.

The two diagnostics are complementary, and their joint reading is the central message of this section. The Fisher analysis quantifies whether \(\Lambda_\text{\tiny 8PN}\) can be measured as a parameter, whereas the mismatch quantifies the bias incurred by leaving it out of the template altogether. The two need not coincide: AP4, for instance, crosses the ET mismatch threshold even though its dynamical TLN is not individually measurable at the \(10 \%\) level according to Fig. 7. The reason is that the mismatch is sensitive to the total dephasing accumulated across the inspiral band, and can flag a systematic bias even when the corresponding parameter uncertainty is large. In practice this means that, for ET sources at \(d_\text{\tiny L} \sim 100 \, {\rm Mpc}\) (\({\rm SNR} \gtrsim 100\)), omitting the dynamical tidal phase will bias the inference of the static Love number—and hence the extracted nuclear EoS—even in regimes where the dynamical tide itself cannot be pinned down. Including the full phase of Eq. 51 in the template is therefore necessary to avoid such biases, irrespective of whether \(\Lambda_\text{\tiny 8PN}\) is independently resolvable.

In summary, our analysis shows that the dynamical Love number is potentially measurable with third-generation detectors for a range of neutron star masses and EoS, and that neglecting it can introduce systematic biases in the inference of nuclear-matter properties. The degree of observability depends sensitively on the EoS, with stiffer EoS yielding both a larger \(\Lambda_\text{\tiny 8PN}\) and a larger mismatch. These findings motivate including the full dynamical tidal phase of Eq. 51 in next-generation GW data-analysis pipelines, and suggest that joint inference of \(\Lambda_\text{\tiny 5PN}\) and \(\Lambda_\text{\tiny 8PN}\) with ET could provide a new window into the frequency-dependent tidal response of dense nuclear matter.

6 Discussion↩︎

In this work we have carried out a systematic and fully relativistic study of the dynamical tidal response of neutron stars by combining worldline EFT with stellar perturbation theory. We computed the dynamical quadrupolar Love number \(c_{\dot{E}}\) for a broad set of realistic EoS. The corresponding dimensionless parameter \(\Lambda_2 \propto C^{-8}\) is enhanced by three additional inverse powers of the compactness relative to the static counterpart \(\Lambda_0 \propto C^{-5}\), being significantly larger than its static counterpart for typical neutron-star compactnesses. The dynamical response coupling also exhibits a universal logarithmic running, in addition to a set of scheme-dependent finite terms, which are unambiguously fixed by the EFT matching and which we determine here for the first time using dimensional regularization. We have shown that, for a fixed running parameter, approximate EoS-insensitive relations exist between \(\Lambda_0\) and \(\Lambda_2\), analogous to the I-Love-Q universality [94], [95], and that, taking motivated values of the running into account, the \(f\)-mode still provides an accurate single-mode approximation to the dynamical response throughout the late inspiral, as in the Newtonian case. Using these results, we derived the complete leading-order 8PN tidal phase correction to the gravitational waveform of a binary neutron-star inspiral, including both the EoS-dependent constant term and the logarithmic running. While the logarithmic term is qualitatively distinctive—in principle allowing one to disentangle the dynamical tide from unknown point-particle contributions at the same PN order—it is numerically subleading during the late inspiral.

The waveform contribution associated with the dynamical Love number is also quantitatively comparable to that arising from the quadratic Love numbers [56], [106]. Both effects enter the GW phase at 8PN order and are enhanced in the low-compactness regime, scaling as \(\sim C^{-8}\). Comparing our results with those of Ref. [56], we find that the two contributions are generally of the same order of magnitude. This suggests that a consistent description of the leading 8PN tidal correction requires both dynamical and quadratic Love numbers to be included in the waveform model.

Our Fisher-matrix analysis shows that ET could measure the 8PN tidal contribution below the \(100\%\) level for \(m_\text{\tiny NS}\lesssim 1.3\,M_\odot\) and \(m_\text{\tiny NS}\lesssim 1.8\,M_\odot\) for the softest and stiffest EoS, respectively. For \(m_\text{\tiny NS}\lesssim \,M_\odot\), the accuracy improves to \(10\%\) and percent level, respectively. In contrast, LVK O4 cannot constrain the dynamical tide for any of the configurations considered. The mismatch analysis reveals a complementary and practically important finding: omitting the 8PN dynamical tide from the waveform template introduces systematic biases in the inference of the static Love number for \(m_\text{\tiny NS}\lesssim 1.2\,M_\odot\) and soft EoS, and larger masses for the stiffer ones, at ET sensitivity. Since the latter is the primary probe of the nuclear EoS in current and planned GW observations, this finding motivates including the full dynamical tidal phase in next-generation waveform models.

Our work can be extended along various directions. Incorporating stellar rotation would introduce gravito-magnetic dynamical tides and mode-mixing relevant for spinning neutron stars. Extending the EFT matching to include tidal corrections to the GW amplitude, beyond the phase computed here, would complete the leading-order dynamical tidal waveform. Finally, combining the dynamical TLN measurement with independent nuclear-physics constraints could significantly sharpen EoS inference with third-generation detectors, making the framework developed here a concrete ingredient for high-precision GW science.

Note added: While this work was being completed, Ref. [107] appeared, containing an independent discussion of the calculation of dynamical TLNs of neutron stars through EFT matching. We note that, differently from the approach of Ref. [107], we work at the level of off-shell quantities rather than on-shell ones, and perform the matching at loop level, resumming gravitational nonlinearities within the EFT.

Acknowledgments↩︎

We acknowledge interesting discussions with Nils Andersson, Abhishek Hegade, Tanja Hinderer, Soumodeep Mitra, Raj Patil, Hector O. Silva, and Nicolas Yunes. V.D.L. and L.S. are grateful to the Institute for Fundamental Physics of the Universe (IFPU) in Trieste for hosting the workshop “Perspectives on Gravity: From Theory to Observation”, where this work was completed. V.D.L. is supported by NSF Grants No. AST-2307146, No. PHY-2513337, No. PHY-090003, and No. PHY-20043, by NASA Grant No. 21-ATP21-0010, by John Templeton Foundation Grant No. 62840, by the Simons Foundation [MPS-SIP-00001698, Emanuele Berti], by the Simons Foundation International [SFI-MPS-BH-00012593-02], and by Italian Ministry of Foreign Affairs and International Cooperation Grant No. PGR01167. This work was carried out at the Advanced Research Computing at Hopkins (ARCH) core facility (https://www.arch.jhu.edu/), which is supported by the NSF Grant No. OAC-1920103. We acknowledge support by the MUR FIS2 Advanced Grant ET-NOW (CUP: B53C25001080001), and by the INFN TEONGRAV initiative. The research of L.S. has been funded, in part, by the French National Research Agency (ANR) under project ANR-24-CE31-1097-01. This work has received support under the program “Investissement d’Avenir” launched by the French Government and implemented by ANR, with the reference ANR-18-IdEx-0001 as part of its program “Emergence”.

References↩︎

---

1. Other examples where static Love number couplings can run include black holes in higher dimensions or in theories beyond general relativity, see e.g. [61]–[67].↩︎

2. This operator coincides with the quadrupolar tensor defined in Ref. [7], up to a normalization factor of \(1/2\).↩︎

3. See Refs. [56], [75]–[79] for a generalization to nonlinear response in the point-particle EFT.↩︎

4. Although we match the EFT to a non-dissipative relativistic star in this work, Appendix 7 presents the general treatment using the Schwinger–Keldysh formalism (see also Ref. [50] for details). This is useful for discussing general properties of the response that hold for any object, independently of the details of the UV matching.↩︎

5. We stress that the couplings on the left-hand side of Eq. 8 are the renormalized ones. For ease of notation, we drop the superscript R that distinguishes them from the bare couplings. The distinction should be clear from the context and should not be a source of confusion.↩︎

6. Note that, although the action 4 is formally expressed as an expansion in the number of time derivatives, one could instead keep the EFT nonperturbative in frequency (while still expanding in spatial gradients) and match the relativistic result directly to the kernel \(K_{+-}^{(E)}\left(\omega\right)\) defined in Appendix 7.↩︎

7. Since we will connect the exterior solution to the interior solution obtained from a non-dissipative stellar model, we have set \(a_1=b_1=0\) in Eq. 80 .↩︎

8. See Eqs. 81 –82 with \(a_1=b_1=0\).↩︎

9. Notice that the definition of \(\omega_*\) in Eq. 45 differs slightly from the quantity introduced in Eq. (1.14) of Ref. [26], as the latter is defined in terms of the tidal constants adopted therein.↩︎

10. Although the renormalization scheme adopted in Ref. [39] differs from ours, a nontrivial consistency check is that both approaches yield the same counterterms and renormalization-group running of the couplings.↩︎

11. Notice that the expression for the 8PN phase differs from the one of Ref. [51] by a typo in their last derivation.↩︎

12. Note that what we refer to as the “near zone” here corresponds to the “intermediate zone” of Ref. [50]. In that work, a finer distinction between near and intermediate zones was necessary because \(r\) extends down to \(r_s\). Here, since \(R_\star\) is never parametrically close to \(r_s\)—namely \(R_\star \gtrsim \mathcal{O}(1)r_s\)—this distinction is not required. We also stress that, despite its name, the near zone covers a wide range of distances, reaching distances parametrically far from \(R_\star\) (\(r\gg R_\star\)), while remaining smaller than \(1/\omega\).↩︎

13. We refer to Ref. [62] for the precise relation between the Zerilli field \(\Psi_{\text{Z}}\) and the fluctuation of the metric tensor in standard perturbation theory of Sec. 3.↩︎