Nonradial oscillations of stratified neutron stars with solid crusts: Mode characterization and tidal resonances in coalescing binaries

TL;DR

This paper models the tidal response of stratified neutron stars with solid crusts, analyzing mode spectra and resonances, and explores implications for gravitational wave signals and electromagnetic precursors in binary mergers.

Contribution

It provides a detailed mode analysis including stratification effects and identifies how various modes influence tidal interactions and crust breaking in neutron star binaries.

Findings

Crust-core interface mode vanishes in stratified stars, replaced by compositional gravity modes.

Fundamental mode dominates energy transfer despite high frequency.

Crustal shear and gravity modes can trigger crust breaking and influence electromagnetic precursors.

Abstract

Dynamical tides of neutron stars in the late stages of binary inspirals provide a viable probe into dense matter through gravitational waves, and potentially trigger electromagnetic precursors. We model the tidal response as a set of driven harmonic oscillators, where the natural frequencies are given by the quasinormal modes of a nonrotating neutron star. These modes are calculated in general relativity by applying linear perturbation theory to stellar models that include a solid crust and compositional stratification. For the mode spectrum, we find that the canonical interface mode associated with the crust-core boundary vanishes in stratified neutron stars and is replaced by compositional gravity modes with mixed gravity-interfacial character, driven primarily by strong buoyancy in the outer core. We also find that fluid modes such as the core gravity mode and the fundamental mode can penetrate the crust, and we establish a criterion for such penetration. Regarding the tidal interaction, we find that transfer of binding energy to oscillations is dominated by the fundamental mode despite its frequency being too high to resonate with the tidal forcing. In general, we find that lower-frequency modes induce gravitational-wave phase shifts smaller than $\sim 10^{-3},\rm rad$ for the equation of state we consider. We discover that nonresonant fundamental and crustal shear modes can trigger crust breaking already near the first gravity-mode resonance, while gravity-mode resonance concentrates strain at the base of the crust and may marginally crack it. These results suggest that both resonant and nonresonant excitations can overstress the crust and may channel energy into the magnetosphere prior to merger, potentially powering electromagnetic precursors. Our work represents an important step toward realistic modeling of dynamical tides of neutron stars in multimessenger observations.

Figures

Tables

Table 1. Table 1: Parameters for the canonical model ( M = 1.35 ​ M ⊙ M=1.35\,M_{\odot} ) and a higher-mass model ( M = 1.8 ​ M ⊙ M=1.8\,M_{\odot} ). The columns list the mass M M , central rest-mass density ρ c \rho_{\rm c} , stellar radius R R , crust–core interface radius R cc R_{\rm cc} , and neutron-drip radius R nd R_{\rm nd} .

$M\left[M_{\odot }\right]$	${\rho }_c\left[{10}^{14}g/{\mathrm{cm}}^3\right]$	$R\left[\mathrm{km}\right]$	$R_{\mathrm{cc}}\left[\mathrm{km}\right]$	$R_{\mathrm{nd}}\left[\mathrm{km}\right]$
1.35	$7.59$	$12.29$	$11.30$	11.86
1.80	$10.5$	$11.89$	$11.31$	11.65

Table 2. Table 2: The spectral and resonance parameters of the f f , s 1 s_{1} – s 2 s_{2} , g 1 g_{1} – g 2 g_{2} , and i i ( Γ = Γ 0 \Gamma=\Gamma_{0} ) modes for a canonical NS with mass M = 1.35 ​ M ⊙ M=1.35\,M_{\odot} in an equal-mass binary system. The columns list the mode frequency f α f_{\alpha} , the tidal overlap integral | Q n ​ ℓ | |Q_{n\ell}| , the maximum mode amplitude | q max α | |q^{\alpha}_{\rm max}| attained after resonance by integrating Eq.   49 , and its approximation from Eq.   53 . Also included are the mode energy after resonance E mode max E_{\rm mode}^{\rm max} , the resulting GW dephasing Δ ​ Φ \Delta\Phi , the energy deposited in the crust E crust E_{\rm crust} , and the post-resonance elastic energy E ela E_{\rm ela} .

mode	$f_{\alpha }\left[\mathrm{Hz}\right]$	$\left|Q_{n2}\right|$	$\left|q_{\max }^{\alpha }\right|\left[\text{<a><span>Eq.</span> <span>49</span></a>}\right]$	$\left|q_{\max }^{\alpha }\right|\left[\text{<a><span>Eq.</span> <span>53</span></a>}\right]$	$E_{\mathrm{mode}}^{\max }\left[\mathrm{erg}\right]$	$\Delta \Phi$	$E_{\mathrm{crust}}\left[\mathrm{erg}\right]$	$E_{\mathrm{ela}}\left[\mathrm{erg}\right]$
$f$	$1785.4$	0.394	$0.102\text{<sup>a</sup>}$	$0.337$	$4.52\times {10}^{51}\text{<sup>a</sup>}$	$-0.04\text{<sup>a</sup>}$	$7.58\times {10}^{49}\text{<sup>a</sup>}$	$1.99\times {10}^{47}\text{<sup>a</sup>}$
$s_2$	$1179.5$	$2.48\times {10}^{-4}$	$3.12\times {10}^{-4}$	$3.00\times {10}^{-4}$	$2.96\times {10}^{46}$	$-1.16\times {10}^{-5}$	$2.76\times {10}^{46}$	$2.67\times {10}^{45}$
$s_1$	$709.8$	$2.81\times {10}^{-4}$	$5.32\times {10}^{-4}$	$5.17\times {10}^{-4}$	$2.65\times {10}^{46}$	$-4.10\times {10}^{-5}$	$2.21\times {10}^{46}$	$2.24\times {10}^{45}$
$g_1$	$187.2$	$2.21\times {10}^{-4}$	$1.25\times {10}^{-3}$	$1.23\times {10}^{-3}$	$8.42\times {10}^{45}$	$-3.64\times {10}^{-4}$	$1.83\times {10}^{45}$	$6.70\times {10}^{42}$
$g_2$	$106.9$	$5.12\times {10}^{-5}$	$4.55\times {10}^{-4}$	$4.56\times {10}^{-4}$	$3.98\times {10}^{44}$	$-6.00\times {10}^{-5}$	$7.93\times {10}^{43}$	$8.90\times {10}^{41}$
$i\left(\Gamma ={\Gamma }_0\right)$	$33.1$	$9.60\times {10}^{-5}$	$2.28\times {10}^{-3}$	$2.28\times {10}^{-3}$	$9.18\times {10}^{44}$	$-2.21\times {10}^{-3}$	$1.95\times {10}^{43}$	$3.95\times {10}^{42}$

Equations

$$p=p(n_{\mathrm{b}},Y_p),$$

$${\Gamma }_0\equiv {\left(\frac{\partial \ln p}{\partial \ln n_{\mathrm{b}}}\right)}_{\beta -\mathrm{stable}}=\frac{ϵ+p}{p}\frac{\mathrm{d}p}{\mathrm{d}ϵ},$$

$${\tau }_{\mathrm{M}}\sim \frac{2\text{ month }}{T_9^6},\text{and}{\tau }_{\mathrm{D}}\sim \frac{20\mathrm{s}}{T_9^4},$$

$${\Gamma }_1$$

$$+{\frac{\partial ({\mu }_e-{\mu }_{\mu })}{\partial n_{\mathrm{b}}}|}_{Y_i}\frac{\mathrm{d}{Y}_{\mu }}{\mathrm{d}{n}_{\mathrm{b}}}].$$

$$\mu =0.1194{\left(\frac{4\pi }{3}\right)}^{1/3}(Ze{)}^2{\left(\frac{1-X_{\mathrm{n}}}{A}\right)}^{4/3}{n}_{\mathrm{b}}^{4/3},$$

$$v_{\mathrm{s}}={\left(\frac{\mu }{ϵ+p}\right)}^{1/2}.$$

$$\mathrm{d}{s}^2=-e^{\nu }\mathrm{d}{t}^2+e^{\lambda }\mathrm{d}{r}^2+r^2\mathrm{d}{\theta }^2+r^2{\sin }^2\theta \mathrm{d}{\varphi }^2,$$

$$T_{ab}^{\mathrm{fluid}}=\left(ϵ+p\right)u_a{u}_b+p{g}_{ab},$$

$$G_{ab}=8\pi T_{ab}^{\mathrm{fluid}},$$

$$\Delta Q=\delta Q+{\mathcal{L}}_{\mathbf{\xi }}Q,$$

$$\delta g_{ab}=-r^{\ell}\,\left(\begin{array}[]{cccc}H_{0}e^{\nu}&i\omega rH_{1}&0&0\\ i\omega rH_{1}&H_{2}e^{\lambda}&0&0\\ 0&0&r^{2}K&0\\ 0&0&0&r^{2}\sin^{2}\theta K\end{array}\right)Y_{\ell}^{{m}}e^{i\omega t},$$

$$\xi^{a}=r^{\ell}\,\left(\begin{array}[]{c}0\\ Wr^{-1}e^{-\lambda/2}\\ -Vr^{-2}\partial_{\theta}\\ -\displaystyle\frac{V}{r^{2}\sin^{2}\theta}\partial_{\phi}\end{array}\right)Y_{\ell}^{{m}}e^{i\omega t},$$

$$\delta G_{ab}=8\pi \delta T_{ab}^{\mathrm{fluid}}.$$

$$\delta T_a{}^{b\text{fluid}}$$

$$+(ϵ+p)\left(u_a\delta u^b+\delta u_a{u}^b\right).$$

$$\Delta g_{ab}=\delta g_{ab}+{\mathcal{L}}_{\mathbf{\xi }}{g}_{ab}=\delta g_{ab}+{\nabla }_a{\xi }_b+{\nabla }_b{\xi }_a,$$

$$\Delta u^a=\frac{1}{2}{u}^a{u}^b{u}^c\Delta g_{bc}.$$

$$\delta u^a=q_b{}^a{\mathcal{L}}_{\mathbf{u}}{\xi }^b+\frac{1}{2}{u}^a{u}^b{u}^c\delta g_{bc},$$

$$q^{ab}=g^{ab}+u^a{u}^b,$$

$$\frac{\Delta \rho }{\rho }=$$

$$-(\ell +1)e^{-\lambda /2}W-r{e}^{-\lambda /2}{W}^{\prime }]Y_{\ell }^m{e}^{\mathrm{i}\omega t}.$$

$$\Delta ϵ=(ϵ+p)\frac{\Delta \rho }{\rho },\delta ϵ=\Delta ϵ-{\rho }^{\prime }{\xi }^r,$$

$$\frac{\Delta p}{p}={\Gamma }_0\frac{\Delta \epsilon }{\epsilon +p},\delta p=\Delta p-p^{\prime }{\xi }^r.$$

$$\Delta s_a{}^b=\frac{1}{2}\left(q_a{}^c{q}^{db}-\frac{1}{3}{q}_a{}^b{q}^{cd}\right)\Delta g_{cd}.$$

$$\delta {\pi }_a{}^b=-2\mu \delta s_a{}^b,$$

$$\delta T_a^b=\delta T_a^{b\text{fluid}}+\delta {\pi }_a^b.$$

$$T_1{Y}_{\ell }^m{e}^{i\omega t}\equiv \delta {\pi }_r^r,\text{and}\frac{1}{r}{T}_2{\partial }_{\theta }{Y}_{\ell }^m{e}^{i\omega t}\equiv \delta {\pi }_r^{\theta },$$

$${\omega }^2=k^2{c}_l^2=k^2\left(c_{\mathrm{s}}^2+\frac{4}{3}{v}_{\mathrm{s}}^2\right)=k^2\left(\frac{\Gamma p}{ϵ+p}+\frac{4}{3}{v}_{\mathrm{s}}^2\right).$$

$${\omega }^2=k^2{v}_{\mathrm{s}}^2,$$

$$f\simeq 10^3\left(\frac{R-R_{\mathrm{cc}}}{\lambda }\right)\left(\frac{1\mathrm{km}}{R-R_{\mathrm{cc}}}\right)\left(\frac{v_{\mathrm{s}}}{10^{8}cm/s}\right)\mathrm{Hz},$$

$$N^2=e^{\nu -\lambda }\frac{p^{\prime }}{ϵ+p}A,A=\frac{p^{\prime }}{p}\left(\frac{1}{{\Gamma }_0}-\frac{1}{{\Gamma }_1}\right),$$

Full text

I Introduction

As the densest stable form of matter known in the Universe, neutron stars (NSs) are expected to have complex interiors dominated by strong interactions inaccessible through terrestrial experiments. Despite extensive experimental and theoretical efforts in probing the nuclear equation of state (EOS), our current knowledge about the EOS at the densities most relevant for determining the internal structure of NSs remains poor [102, 122, 19].

Compact binary mergers involving NSs provide multiple pathways to probe the dense matter, from subtle finite-size effects in the late inspiral [65, 43, 3, 63, 71, 135, 150, 101, 83, 64, 155] to the dramatic dynamics of merger remnants [149, 162, 18, 23, 38, 70], tidal disruption [99, 100, 98], and mass ejection [69, 81, 148]. The landmark detection of gravitational waves (GWs) from the binary NS inspiral GW170817 by the LIGO-Virgo-KAGRA Collaboration [3], along with its postmerger electromagnetic (EM) counterparts [4], marked the beginning of GW multimessenger astronomy, and has helped shed a new light on the EOS problem. However, the constraints on the EOS set by GW170817 are limited by the insufficient sensitivity of current GW detectors at the frequency band of the late-inspiral and merger, where matter effects on the GW waveform become most salient.

An upgrade to the current ground-based detectors is pending for the coming years, and the next observing run of the International Gravitational-Wave Observatory Network, tentatively planned for 2028, will provide an opportunity to detect GW signals emanated during the merger from GW170817-like events. A central question to be addressed is then: To what extent can observations of binary mergers truly reveal the composition and state of matter inside NSs? Because a NS’s characteristic oscillation frequencies are determined by its internal structure, identifying quasiperiodic oscillations (QPOs) in some observations offers a way to probe these extreme environments. This article focuses on a subset of the stellar spectrum that may imprint observable signatures in gravitational waveforms or EM flares through tidal interactions.

Much work has focused on the tidal deformation of NSs during binary interaction in the static limit [65, 43, 63]. The secular imprint on the GW signal associated to it is characterized by the so-called tidal deformability. The GW170817 event placed a constraint on the binary’s mass-weighted average tidal deformability, disfavoring extremely stiff EOSs [15, 3, 33, 185, 49]. Combined with lower bounds on the maximum mass from radio pulsar timing [34, 16, 44, 32] and joint mass–radius constraints from NICER X-ray timing [113, 112, 141], these observations form the current picture of EOS constraints from the astrophysical side.

The tidal response characterized by the tidal deformability is only the leading order effect obtained under the assumption that the tidal interaction is adiabatic. As the binary approaches the merger, this assumption gradually breaks down, and dynamical, frequency-dependent contributions in the tidal response function must be taken into account [64, 155, 130, 62]. In particular, NSs respond dynamically to the tidal field through the excitation of various oscillation modes as the orbital frequency sweeps through the detector band [150, 101, 83, 170, 28, 94]. Therefore, dynamical tidal interactions provide a realization of NS asteroseismology [13, 14, 7]. That is, the study of the interior composition the star through its oscillation spectrum.

Dynamical tides transfer energy and angular momentum from the orbit into oscillation modes of the star, resulting in a GW phase shift. Many studies (see, for instance, [155, 64, 9, 133]) have shown that the fundamental ($f$-) mode — acoustic oscillations driven by pressure gradients — can be appreciably excited even without reaching exact resonance during the inspiral. Although the resulting impact on the GW signal is relatively modest, incorporating dynamical tides has been demonstrated to improve waveform accuracy [64]. In rotating NSs, the $f$-mode spectrum splits into prograde and retrograde ones [87]. Of particular interest is the retrograde $f$-mode, which appears at lower frequencies in the inertial frame of the observer and can enter into resonance with the orbital frequency when the spin of the NS is antialigned with respect to the orbital angular momentum. This resonance is expected to produce a much larger GW phase shift [67, 108, 156, 91], and was recently captured in high-precision numerical simulations for the first time [90].

The interior of NSs is more complicated than a prescribed density-pressure relation. Matter densities span more than 14 orders of magnitude from the surface to the core, accompanied by changes in composition, from ordinary baryons to possibly exotic particles such as hyperons or quarks, as well as by distinct phase structures, such as solid and superfluid. Many low-frequency modes (compared to $f$-mode) restored by different physical mechanism may become tidally excited. Examples include gravity ($g$-) modes that arise from compositional gradients [150, 140, 180, 104, 31, 183, 182, 181], interface ($i$-) modes associated with discontinuities in density [111, 30, 128, 181, 103] or shear modulus [170, 123, 186, 103], and shear ($s$-) modes restored by shear stresses of the solid crust [86].

The tidal coupling of the foregoing modes is typically much smaller than that of the $f$-mode, but their lower frequencies make resonance more accessible earlier in the inspiral, and thus potentially detectable with current GW detectors. Several works showed that the GW phase shift induced by compositional $g$-modes is generally smaller than $10^{-3}$ radians; see, e.g., Refs. [180, 31, 93]. However, Refs. [97, 96] suggested that nonlinear effects could lock the $g$-mode resonance resulting in several radians of dephasing. Moreover, Refs. [170, 123, 186, 111, 30, 128, 186] suggested that the $i$-mode associated with the crust-core interface or a first-order phase transition may produce approximately $1$ to $10$ radians of GW dephasing. This offers a promising probe of the existence of phase transitions to quark matter, as well as into nuclear parameters around the saturation density.

Mode excitation may still yield observable signatures even if the associated GW dephasing is negligible, provided that the mode’s energy is efficiently converted into coherent EM emission. In particular, if a mode acquires a large amplitude during resonance, the resulting stresses can break the star’s crust, excite high-frequency crustal oscillations, and trigger bursts of EM waves radiation through coupling with the magnetosphere. This mechanism, first proposed for $i$-modes [170, 119], has been invoked to explain the precursor emission observed before some short gamma-ray bursts (sGRB) [169, 85, 169, 114, 184, 29, 173, 105, 174, 177, 168]. It has since been extended to other modes, such as $g$- [93, 94, 95, 92] and $f$-modes; see Ref. [161] for a review. Although the overstressing of the crust is triggered by resonant tidal excitation rather than slow magnetic-field evolution, the subsequent magnetic coupling and EM energy release are analogous to the mechanisms proposed for magnetar giant flares [163]. The recent detection of a 22 Hz QPO [177] in a sGRB may have originated from shear modes [160] or Alfven modes excited following a crust breaking, echoing similar interpretations proposed for the QPOs in the giant-flare tails of magnetars [73, 158].

In this paper, we will investigate the dynamical tides of nonrotating NSs and their potential observational consequences. Our work improves upon the previous literature in three aspects. First, we use a unified EOS that includes both a solid crust and compositional stratification, while earlier works often treated them separately or introduced stratification phenomenologically. Second, we employ the framework of relativistic linear perturbation theory, without approximations such as the Cowling approximation or the use of Newtonian perturbations on relativistic backgrounds. Third, we analyze how resonant and nonresonant modes contribute to the breaking of the crust. Together, these improvements give a more realistic description of dynamical tides in compact binaries.

The paper is organized as follows. In Sec. II, we describe the EOSs and the formalism to calculate the full set of quasi-normal modes (QNMs) in relativistic linear perturbation theory. The properties of both high frequency (fundamental and shear) and low frequency (gravity and interface) modes are discussed in Secs. III and IV, respectively. We then investigate the tidal excitation of different modes in Sec. V. Therein, we include a discussion of energy transfer and GW dephasing in Sec. V.1, and a detailed study of crust breaking in Sec. V.2. We summarize our finding in Sec. VI. In Appendix A, we present the master equations and describe the numerical methods we used to compute the mode frequencies. Throughout this paper, we use geometrical units $c=G=1$, where $G$ is the gravitational constant and $c$ is the speed of light.

II Neutron star model and its oscillation spectrum

II.1 The equation of state

We employ a unified EOS predicted by the relativistic density functional TW99, that adopts density-dependent couplings [120]. The EOS describes cold $npe\mu$ matter (i.e., consisting of neutrons, protons, electrons, and muons) which fulfills the (global) charge neutrality condition. The adopted relativistic-mean-field theory parameters and the corresponding nuclear saturation properties are given in Tables I and II of Refs. [120, 175, 176]. The EOS is formulated in a two-parameter form,

[TABLE]

where $p$ is the pressure, $n_{\rm b}$ is the baryon number density, and $Y_{p}=n_{p}/n_{\mathrm{b}}$ is the proton fraction with $n_{p}$ the proton number density. Similarly, we also express the energy density $\epsilon=\epsilon(n_{\mathrm{b}},Y_{p})$ and chemical potential $\mu_{i}=\mu_{i}(n_{\mathrm{b}},Y_{p})$, for a fermion $i$, in a two-parameter form. Charge neutrality requires that $Y_{p}=Y_{e}+Y_{\mu}$, while the chemical equilibrium of leptons requires that $\mu_{e}=\mu_{\mu}$. For $\beta$-stable NS matter, that is, matter fulfilling $\mu_{n}-\mu_{p}=\mu_{e}=\mu_{\mu}$, the thermodynamic quantities such as the pressure and energy density become functions of density alone.

II.1.1 Adiabatic index and stratification

The adiabatic index $\Gamma_{0}$ of NS matter is defined as

[TABLE]

where the derivative is taken along the $\beta$-stability. This adiabatic index can be deemed as a measure of an EOS’s stiffness as a function of $n_{\rm b}$ (larger $\Gamma_{0}$ is stiffer), and is thus closely related to the bulk properties of an NS in equilibrium, such as its mass $M$ and radius $R$. In Fig. 1, we show $\Gamma_{0}$ as a function of the rest-mass density $\rho=m_{u}n_{\mathrm{b}}$, where $m_{u}=931.494$ MeV is the atomic mass unit, for the TW99 EOS. We see that at the outer crust, where densities are smaller than the neutron drip density $\rho_{\rm nd}\simeq 3.0\times 10^{11}\,\mathrm{g/cm^{3}}$, the adiabatic index is roughly constant $\Gamma_{0}\simeq 4/3$ as the pressure is dominated by the relativistic electron Fermi gas.111The density discontinuity arising from the phase transition between different species of nuclei [20] has been smoothed out in this region since the nucleon numbers in the nuclei vary smoothly with density and take noninteger values in our EOS. As the density increases toward the neutron drip point, the EOS softens because of the emergence of free neutron gas, leading to a sharp drop in $\Gamma_{0}$. Beyond this point, the matter stiffens as the neutron gas contributes to increase the Fermi pressure with additional contributions from nuclear force and Coulomb interactions. At the crust-core interface where $\rho_{\rm cc}\simeq 1.03\times 10^{14}\,\rm g/cm^{3}$, the adiabatic index $\Gamma_{0}$ first drops abruptly and then rises sharply as nuclei dissolve. The adiabatic index increases to $\Gamma_{0}\simeq 3$ at densities $\rho\gtrsim 5\times 10^{14}\,\mathrm{g/cm^{3}}$ as repulsive nucleon interactions become progressively stronger, before gradually decreasing at higher densities wherein the EOS softens due to the density dependence of nuclear interactions.

If the nuclear reactions occur rapidly enough to maintain chemical equilibrium as a fluid element is perturbed, $\Gamma_{0}$ can also be used to describe the perturbations. However, current studies suggest that the equilibration times for nonsuperfluid matter are of the order of

[TABLE]

for the modified and direct Urca reactions, respectively [139, 58, 59]. Here, $T_{9}$ is the temperature in units of $10^{9}$ K. Inspiralling NSs are old and cold, typically with $T_{9}\lesssim 0.01$, which makes both $\tau_{\mathrm{M}}$ and $\tau_{\mathrm{D}}$ much longer, and suggests that the adiabatic index for perturbations is different from that for the background, unperturbed matter. Quantitatively, the typical timescale $\tau_{\rm osc}$ of the oscillations considered in this work ranges from $10^{-4}$ to $10^{-1}$ s. Comparing these timescales with Eq. (3) we conclude that the weak reactions are not fast enough to maintain $\beta$-equilibrium between the perturbed fluid element and its surroundings. For our purposes, it is therefore reasonable to assume that the composition of the fluid element remains effectively frozen during oscillations. The adiabatic index for the perturbed matter $\Gamma_{1}$ is then:

[TABLE]

The derivatives are taken at fixed particle fractions with $Y_{i}=n_{i}/n_{\mathrm{b}}$, and the second equality assumes $npe\mu$ matter [74, 159].

As shown in Fig. 1, $\Gamma_{1}$ generally follows the behavior of $\Gamma_{0}$ across most densities. However, it notably deviates from the latter near the neutron drip point and outer core close to the crust–core interface, when free neutrons and muons emerge, respectively. Around the neutron drip, the difference $\Gamma_{1}-\Gamma_{0}$ can be of order unity, while near the crust–core boundary it remains approximately $0.2$.

The difference between the adiabatic index of the background matter and that of perturbed fluid elements results in a displaced fluid element with a density different from its surroundings. The magnitude of this difference sets the strength of the buoyancy and determines the spectrum of compositional $g$-modes; see Sec. IV.2.

II.1.2 Elasticity and shear modulus

The crust of NSs behaves as a solid due to the presence of a Coulomb lattice. Assuming a body-centered cubic lattice, the shear modulus of a zero-temperature crust can be estimated as

[TABLE]

where $Z$ is the mean charge number of the nuclei, $A$ is the mean nucleus number, and $X_{\rm n}$ denotes the fraction of “free” neutrons [157, 154].

The shear speed $v_{\rm s}$ characterizes the velocity of elastic waves in the crust; its relativistic value can be computed as [27]

[TABLE]

As illustrated in Fig. 2, the shear speed is of the order of $10^{-3}$–$10^{-2}$ of the speed of light for densities $\sim 10^{8}$–$10^{14}\,\rm g/cm^{3}$. Below the neutron drip density, the shear speed gradually increases with $\rho$ because of the increasing mean charge number as well as the density. However, the fraction of free neutrons grows rapidly when $\rho$ exceeds $\rho_{\rm nd}$, which progressively suppresses the shear speed.

The solid crust is important to NS phenomenology for several reasons. Its ability to sustain large elastic stresses allows the formation of nonaxisymmetric “mountains” [75, 54, 56], which can both modulate pulse profiles through free precession [76, 51, 50, 35] and act as sources of continuous GWs [51, 77]. Because the crust supports shear motion, it hosts oscillation families that cannot exist in a perfect fluid, most notably shear ($s$-) and torsional ($t$-) modes [109, 146, 158, 144, 152, 151]. Moreover, the sharp transition between the solid crust and the underlying fluid core gives rise to additional $i$-modes [109, 170, 129].

II.2 Polar perturbation and oscillation spectrum

Having described the EOS we will use and presented some of its properties, we now explain how we construct our equilibrium NS models and how we compute their oscillation spectra.

The spacetime line element of a static spherically symmetric NS in hydrostatic equilibrium can be expressed as

[TABLE]

where $\nu$ and $\lambda$ are the metric functions that depend only on the radial coordinate $r$. We assume that the matter content of star can be described as a perfect fluid, with energy-momentum tensor

[TABLE]

where $u^{a}$ is the four-velocity and $g_{ab}$ the spacetime metric.

Through Einstein’s equations,

[TABLE]

where $G_{ab}$ is the Einstein tensor with $u^{a}=u^{t}\delta^{a}_{~t}$, these assumptions lead to the Tolman–Oppenheimer–Volkoff (TOV) equations that describe an NS in equilibrium [166, 121]. These equations can be readily integrated once an EOS has been specified. We note that strictly speaking, the perfect fluid assumption is inconsistent with the presence of a solid component. Nevertheless, it still makes sense to assume the background is in a relaxed state, which may be reasonable for NSs close to merger.

Nonrotating NSs modeled with the TW99 EOS have a maximum mass of $M_{\rm TOV}=2.08\,M_{\odot}$, consistent with the most massive pulsars observed to date [44, 16, 34, 32]. Besides, for an NS with mass $1.35\,M_{\odot}$, the corresponding radius is $R_{1.35}=12.29\,\mathrm{km}$, and the dimensionless tidal deformability is $\Lambda_{1.4}=510.0$, in agreement with observations X-ray timing of NSs from NICER [112, 141, 172, 143] and from the binary NS inspiral GW170817 [3].

The oscillations of nonrotating NSs are typically studied by linearizing Einstein’s equations on a background equilibrium solution. We employ the relativistic Lagrangian perturbation formalism developed in Refs. [45, 11, 10, 46]. The key dynamical quantities in this framework include the metric perturbation, $\delta g_{ab}$, and the fluid displacement vector, $\xi^{a}$. For any spacetime quantity $Q(t,x)$, the relationship between the Lagrangian perturbation, $\Delta Q(t,x)$, and the Eulerian perturbation, $\delta Q(t,x)$, is given by

[TABLE]

where $\mathcal{L}_{\bm{\xi}}$ denotes the Lie derivative along the displacement vector $\xi^{a}$. Because of the spherical symmetry of the background, we decompose the perturbations into spherical harmonic and assume they have an harmonic time dependence, that is, $Y_{\ell}^{{m}}e^{i\omega t}$, where $\omega$ is the mode’s oscillation angular frequency and $Y_{\ell}^{{m}}$ is the spherical harmonics with quantum numbers $\ell$ and $m$. By adopting the Regge–Wheeler gauge [138], we can express the metric perturbation of polar parity as

[TABLE]

where $H_{0}$, $H_{1}$, $H_{2}$, and $K$ are functions of $r$ only, and that dictate the perturbations in the metric. The associated displacement vector has the form

[TABLE]

where $W$ and $V$ are functions of $r$.

For a perfect fluid, the master equations for $\{H_{0}$, $H_{1}$, $H_{2}$, $K\}$ and $\{W$, $V\}$ can be derived by solving the perturbed Einstein’s equations,

[TABLE]

The perturbed Einstein tensor $\delta G_{ab}$ can be obtained from the background metric in Eq. 7 and the perturbed metric Eq. 11. The perturbed energy-momentum tensor $\delta T_{ab}^{\text{fluid}}$ can be written in terms of the background four-velocity, pressure, energy density, and their corresponding perturbations as

[TABLE]

The Lagrangian perturbation of the metric, given by,

[TABLE]

is related to the Lagrangian perturbation of the four-velocity as [45]

[TABLE]

Using Eq. 10, we obtain

[TABLE]

where

[TABLE]

is the projection tensor orthogonal to the four-velocity $u^{a}$.

To obtain the expressions for $\delta\epsilon$ and $\delta p$, a useful relation is the Lagrangian perturbation of the rest-mass density [164],

[TABLE]

where a prime indicates differentiation with respect to $r$. With this expression, the Lagrangian and Eulerian perturbations of the energy density are

[TABLE]

respectively. The pressure perturbation is related to the the energy density perturbation by the adiabatic index. However, which adiabatic index, $\Gamma_{0}$ or $\Gamma_{1}$ should we use? If we assume that matter is in $\beta$-equilibrium during the perturbation, the pressure perturbation is then given by

[TABLE]

On the other hand, for frozen composition, the pressure perturbation can be obtained by replacing $\Gamma_{0}$ with $\Gamma_{1}$ in the equation above. In reality, the adiabatic index should be somewhere in between because of the finite weak reaction timescale. However, as we argued in Sec. II.1.1, the oscillation periods of the modes under consideration are likely to be much shorter than the typical weak reaction times in inspiralling NSs for our EOS. We will therefore focus on the QNM spectrum computed with $\Gamma_{1}$ for the perturbations, and use the spectrum with $\Gamma_{0}$ as a reference when investigating the impact of stratification.

What about the perturbation of solid crust? The foundations of elasticity and its perturbation in GR were formulated by Carter and Quintana in the 1970s [25, 27]; see also Refs. [26, 79, 21, 22] for subsequent developments and Ref. [11] for a review. For a perturbed spherical NS, the Lagrangian perturbation of the strain tensor is given by

[TABLE]

Assuming a Hookian elastic solid, the stress tensor is given by

[TABLE]

where $\delta s_{a}{}^{b}$ is the Eulerian perturbation of the strain tensor. The Eulerian perturbation of the total energy-momentum tensor is then:

[TABLE]

We follow Refs. [42, 86, 109] to further define the radial traction $T_{1}$ and tangential traction $T_{2}$ as

[TABLE]

respectively. Both of them vanish in the fluid and are useful for the determination of the junction conditions at the solid-fluid interface. The remaining components of $\delta\pi_{a}{}^{b}$ can be expressed in terms of $\{T_{1}$, $T_{2}$, $V\}$, the shear modulus $\mu$, and the background metric functions. In summary, we have two additional variables $\{T_{1}$, $T_{2}\}$ describing the perturbations of the solid crust besides the six variables $\{H_{0}$, $H_{1}$, $H_{2}$, $K$, $W$, $V\}$ describing the perturbation of the fluid core. To avoid diverting further from the main topic, we present the master equations for both fluid and solid perturbations, the junction conditions at the solid–fluid interface, and the numerical methods we used to compute the oscillation frequencies $\omega$ in Appendix A. These frequencies are complex valued, known as QNMs. The real part of the frequencies are related to the oscillations of the system, and the imaginary part of the frequencies are related with the energy dissipation through GW emission. In our problem, this damping timescale is much longer than the star’s oscillation and binary-inspiral timescales. For this reason, we only consider the real part of the QNMs in this work.

In Fig. 3 we plot the real parts of the mode spectra for NSs with $M=1.35\,M_{\odot}$ (top) and $M=1.8\,M_{\odot}$ (bottom), both with a solid crust. The $1.35\,M_{\odot}$ case will be referred as the canonical model, and the background parameters for both stars are listed in Table 1. We adopt the adiabatic index $\Gamma_{1}$ for the perturbations of stratified NSs. The vertical axis shows the amplitude $A_{\rm in}$ of the ingoing GWs at the stellar surface. As we explain in Appendix A, a frequency $f$ for which $|A_{\rm in}|=0$ represents a QNM of the star. These are easily visible in the figure as narrow kinks. We separately display the low-frequency band ($\lesssim 300\,\rm Hz$), which consists of $g$-modes, and the high-frequency band ($\gtrsim 800\,\rm Hz$), which includes acoustic waves ($f$- and $p$-modes) and $s$-modes. To examine the influence of stratification on the spectrum, we overlay our results for a nonstratified (i.e., $\Gamma=\Gamma_{0}$) NS of the same mass. We observe that stratification has little impact on the high-frequency band as the restoring forces for these modes are much stronger than buoyancy. Unsurprisingly, in the low frequency band, $g$-modes, that have buoyancy as restoring their force, cease to exist. However, an $i$-mode is revealed in the nonstratified scenario. We explore in detail the interplay between $i$- and $g$-modes in Section IV.3.

For comparison, we also show the results for a more massive NS with $M=1.8\,M_{\odot}$ in the lower panel of Fig. 3. The overall picture still applies, though the high-frequency part of the spectrum shifts to higher frequencies, while the low-frequency part shifts to lower frequencies. We also observe that the $s$-mode spectrum is noticeably sparser in the heavier NS than in the lighter one as we discuss in Section III.

In the next two sections, we will discuss in detail the high and low-frequency modes separately.

III Characteristics of high-frequency modes

The high-frequency band in the right panel of Fig. 3 accommodates acoustic modes (i.e., the $f$-mode and $p_{1}$-mode) that are restored by pressure gradient as well as the shear modes that are restored by the shear force of the solid crust. In the short-wavelength approximation, the dispersion relation of acoustic waves in Newtonian gravity is given by [109]

[TABLE]

Here, $k$ is the wave number, $c_{l}$ is the wave propagation speed, and $c_{\rm s}$ is the adiabatic sound speed where $\Gamma$ refers either to $\Gamma_{0}$ or $\Gamma_{1}$ for nonstratified and stratified cases, respectively. On average, the shear velocity $v_{\rm s}$ is about 5% of the adiabatic sound speed $c_{\rm s}$, so the solid crust has little effect on the acoustic modes. Indeed, we verified that the frequencies of the $f$- and $p_{1}$-modes are altered by less than 0.01% for the canonical model when we artificially set the shear modulus $\mu$ to zero. Stratification also has a negligible effect on the $f$-mode, but changes the $p_{1}$-mode frequency by about 1%. In the upper panel of Fig. 4, we show the $f$-mode frequency as a function of the NS mass, which increases monotonically as heavier NSs are more compact and have a higher average sound speed. In the lower panel of Fig. 4, we show the real parts of the radial, $\xi^{r}$, and tangential, $\xi^{\theta}$, eigenfunctions for the canonical model and the model with $M=1.8\,M_{\odot}$. The radial eigenfunction grows from the center of the star located at the right end of the plot and slightly decreases in the outer crust. On the other hand, the tangential eigenfunction is nearly constant in the fluid core, and its absolute value increases sharply in the outer crust. This behavior is important to study the breaking the outer crust because the shear strain is proportional to the gradient of the eigenfunction with respect to the radial coordinate, as we will discuss in Section V.2.

As we saw in Fig. 3, the $s$-modes are approximately evenly-spaced in frequency. Their dispersion relation in Newtonian gravity is given by

[TABLE]

where $v_{\rm s}$ is the shear wave speed [109]. The corresponding frequency can be estimated as

[TABLE]

where $R-R_{\rm cc}$ is the the thickness of the crust, that is, the difference between the star’s radius $R$ and the location of the crust-core interface $R_{\rm cc}$. For low overtones of $s$-modes, the wavelength $\lambda$ scales with the crust thickness, resulting in a characteristic frequency of approximately 1 kHz. Because of this scaling, the frequency of a given $s$-mode is higher in the heavier NS, as its crust is thinner. Moreover, the s-mode spectrum is sparser for the more massive star than in the lighter one. This is expected because the spacing between neighboring $s$-modes scales approximately as $1/(R-R_{\rm cc})$.

In the upper panel of Fig. 5, we show the frequencies of the $s_{1}$ and $s_{2}$ modes as functions of the mass. The mode frequencies increase with mass due to the progressive thinning of the crust, which reduces the wavelength and thus increases the frequency. As shown in Fig. 3, stratification has a negligible effect on the $s$-mode spectrum, as the shear stresses in the solid crust are typically an order of magnitude larger than the local buoyancy forces in this region. In the lower panel of Fig. 5, we present the real parts of the radial and tangential eigenfunctions of the $s_{1}$ and $s_{2}$ modes for the canonical model. The oscillations are confined almost completely within the crust because only the solid can sustain shear stresses. The tangential displacement is one order of magnitude larger than the radial displacement. For this reason we divide it by a factor of 10 for for the purpose of presentation.

IV Characteristics of low-frequency modes

We now turn our attention to the low-frequency band in the left panel of Fig. 3 that includes the $i$-mode and a subset of $g$-modes (viz., $g_{1}$ to $g_{11}$). We examine the spectral properties and eigenfunctions of these modes in Secs. IV.1 and IV.2, respectively. As we will see, an interesting feature of the low-frequency modes is that the $i$-mode appears only in nonstratified models, while it is eliminated by stratification. However, its eigenfunction characteristics are imprinted on the $g$-modes, particularly at higher overtones. This phenomenon is further discussed in Sec. IV.3.

IV.1 $i$-mode of nonstratified NSs

The $i$-modes were first studied in geophysics where they are known as Scholte waves at the interfaces between elastic and fluid media [145], and have since been investigated in the context of NSs (e.g., Refs. [109, 170, 171, 178, 88, 153, 124, 186]). The basic picture of the crust-core interfacial mode is a surface wave propagating along the core-crust interface. As a pressure perturbation travels across the crust-core boundary, it compresses the crust, leading to a pronounced radial displacement at the interface. This interficial mode of oscillation is eliminated by buoyancy when NS stratification is taken into account. Nevertheless, its characteristics manifest in the overtones of $g$-modes as will be further demonstrated later. Therefore, all results we present about the $i$-modes are obtained by setting the buoyancy to zero, in other words, using $\Gamma=\Gamma_{0}$.

We depict $i$-mode’s eigenfunctions in the bottom panel of Fig. 6 for the canonical model. The radial displacement $\xi^{r}$ decays into both the core and the crust, as has a characteristic kink-like structure at the core-crust interface. The tangential displacement $\xi^{\theta}$ is significantly excited in the crust but is highly suppressed in the core due to the abrupt change in shear modulus, resulting in a discontinuity at the interface. Since the stress is related to the spatial derivative of displacements, the large displacement gradient at the the base of the crust can induce large shear stresses there. Some literature (e.g., Refs. [170, 123]) thus suggests that the resonantly excited $i$-mode can potentially yield the crust breaking if the amplitude is large enough.

We also present the $i$-mode frequency as a function of mass in the upper panel of Fig. 6. The frequency decreases with increasing mass, primarily due to crust thickening. For NSs with $M\gtrsim 1.0\,M_{\odot}$, the $i$-mode frequency remains below approximately 40 Hz. These frequencies are of the same order of magnitude as those reported in Refs. [186, 88, 153], but are in tension with Refs. [116, 117], which report $i$-mode frequencies in the range of 100–200$\,\rm Hz$ for a broad set of nuclear parameters. The inconsistencies may arise from the EOS prescriptions or the Newtonian calculations used in these latter works. Understanding this inconsistency is left for future. A difference of an order of magnitude in these mode frequencies would significantly impact the resonance onset during binary inspirals, and would thus render very different observational implications. Should an observation be made, the difference could also result in biased Bayesian inferences about the underlying nuclear parameters, which determine the onset density of the crust-core transition that the $i$-mode frequency is sensitive to. Therefore, this discrepancy warrants further investigation. A systematic study of the $i$-mode frequency band across a broad range of EOSs will be carried out elsewhere.

IV.2 $g$-mode of stratified NSs

Once stratification is included, a perturbed fluid element experiences buoyancy due to the composition gradient, giving rise to a spectrum of $g$-modes. The characteristic angular frequency of the oscillating fluid parcel, known as the Brunt–Väisälä frequency $N$, is given by [72]

[TABLE]

where $A$ is the relativistic Schwarzschild criterion [36, 93, 48]. If $\Gamma_{1}>\Gamma_{0}$ (i.e., $A<0$ and $N^{2}>0$), the fluid is stable against convection and supports the propagation of $g$-modes.

As shown in Fig. 1, the TW99 EOS satisfies the condition $\Gamma_{1}-\Gamma_{0}\geqslant 0$, with the difference becoming particularly pronounced near the neutron drip point and in the outer core around $\rho\sim 2\times 10^{14}\,\mathrm{g/cm^{3}}$. Neglecting the solid phase, we expect two distinct classes of $g$-modes propagating in the crust and core, respectively, as previously discussed by Counsell et al. [31] using the BSK family of EOSs [57, 126, 132, 40, 147]. We perform a similar analysis for the TW99 EOS by setting the shear modulus to zero. The resulting mode spectrum is shown in Fig. 7, where the core $g$-modes and crust $g$-modes are labeled as $g^{\,\rm co}_{n}$ and $g^{\,\rm cr}_{n}$, respectively. To distinguish between the two families, we artificially suppress the buoyancy in a selected region by setting $\Gamma_{1}=\Gamma_{0}$. For instance, when the buoyancy around the neutron-drip layer is removed, the crustal $g$-modes disappear from the spectrum while the core $g$-modes persist. Conversely, eliminating the buoyancy in the outer core leaves only the crustal $g$-modes. The subscript $n$ labels the overtone number of the mode. The frequencies of the $g$-modes decrease monotonically with respect to the mass, except near the maximal mass, but the dependence on the mass is weak. This weak dependence can be understood by noting that $g$-mode frequencies scale approximately with the Brunt–Väisälä frequency $N$ defined in Eq. 29, which itself scales as $N^{2}\sim gA$, where $g$ is the local gravitational acceleration and $A$ characterizes the stratification. Since $N$ is primarily determined by the composition gradient in the outer core, and $g$ does not vary significantly across models with different masses, the $g$-mode frequencies remain relatively insensitive to the stellar mass.

The eigenfunctions of the two lowest-order modes for the same background model as in Fig. 3 are displayed in Fig. 8. For the crust $g$-modes, the amplitudes of the radial and tangential eigenfunctions are primarily concentrated in the outer crust, especially for the $g^{\rm cr}_{1}$-mode. The morphology of the eigenfunctions near the neutron drip point resembles that expected for $i$-modes, with a kink structure in $\xi^{r}$ and a sharp jump in $\xi^{\theta}$ at the neutron drip point. For the core $g$-modes, the amplitude of the radial eigenfunctions peaks in the core and becomes small in the crust, while the amplitude of the tangential eigenfunctions resides in both the core and crust, indicating that the buoyancy force in the core is sufficiently strong to penetrate into the crust region if we neglect the elasticity.

Once elasticity is included in the crust, the crust $g$-modes disappear from the spectrum, and the core $g$-modes dominate the entire $g$-mode spectrum. To distinguish this case from the perfect-fluid model, we refer to the $g$-modes in NSs with a solid crust as “full $g$-modes” and denote them as $g_{n}$. As shown in Fig. 7, the frequencies of the full $g$-modes closely follow those of the core $g$-modes, with increasing deviations for overtones. We further illustrate the eigenfunctions of the $g_{1}$–$g_{6}$ modes in Fig. 9. Compared to the $g^{\rm co}_{1}$ and $g^{\rm co}_{2}$ modes shown in Fig. 8 for the fluid case, the amplitudes of both $\xi^{r}$ and $\xi^{\theta}$ are heavily suppressed in the crust region. The amplitude within the crust is further quenched for higher overtones (corresponding to lower mode frequencies), and the eigenfunctions near the crust-core interface begin to exhibit features of $i$-modes.

To understand the behavior of the $g$-mode eigenfunctions in the solid crust, we compare the mode’s transverse momentum with the contribution from crustal shear stresses. For simplicity, we adopt the short-wavelength approximation, in which the eigenfunction varies as $\exp(ikr)$. According to Eqs. (22) and (23), the dominant tangential contribution from the solid can be approximated as

[TABLE]

where $\mu$ is the shear modulus. This term strongly modifies the eigenfunctions when it becomes comparable to the transverse acceleration, which scales as $(\epsilon+p)\omega^{2}\xi^{\theta}$. This yields a critical frequency,

[TABLE]

where $\lambda$ is the mode’s wavelength and $v_{\rm s}$ is the shear speed. If the mode frequency $f$ is comparable to $f_{\rm crit}$, the mode can still propagate into the crust but is strongly affected by elasticity. On the other hand, for $f\ll f_{\rm crit}$, the mode is effectively excluded from the crust, or its frequency band is substantially altered.

The shear velocity in the crust is typically of the order of $10^{8}\,\mathrm{cm\,s^{-1}}$. For the lowest-order core $g$-modes, the frequency is around $200\,\mathrm{Hz}$, with a characteristic wavelength $\lambda\sim R$, leading to a critical frequency $f_{\rm crit}\sim 100\,\mathrm{Hz}$. In this case, the mode can still penetrate into the crust, but its amplitude is suppressed by a factor of $\sim 2$ due to the elasticity. For higher-order core $g$-modes, the frequency becomes increasingly small relative to $f_{\rm crit}$, resulting in stronger suppression of the eigenfunction in the crust. In contrast, the lowest-order crust $g$-modes have frequencies on the order of $300\,\mathrm{Hz}$, but much shorter wavelengths, $\lambda\sim 0.05R$, yielding a much higher critical frequency, $f_{\rm crit}\sim 2000\,\mathrm{Hz}$. As a result, these modes are unlikely to propagate in the crust once elasticity is included. Indeed, we find that crust $g$-modes vanish from the spectrum in the elastic case. However, due to numerical noise dominating the spectrum below $10\,\mathrm{Hz}$, we cannot rule out the possibility that they survive with very low, but finite, frequencies.

IV.3 Interplay between $g$- and $i$-mode

As discussed earlier, the $i$-mode seemingly “disappears” when stratification is included. At the same time, the eigenfunctions of high-order $g$-modes — whose frequencies are close to that of the $i$-mode — begin to resemble the characteristic behavior of $i$-modes near the crust-core interface. This interesting interplay between $i$- and $g$-modes arises from the competition between their respective restoring forces. A simple way to assess the relative strength of these forces is to compare the characteristic frequencies of the mode families they support. To this end, we introduce a phenomenological modification to the adiabatic index,

[TABLE]

and examine how the strength of the buoyancy force, and consequently the $g$-mode frequency, changes with varying $\delta$.

In Fig. 10, we present the relationship between the frequencies of the first three $g$-modes and the control parameter $\delta$. As $\delta$ decreases from $10^{-2}$ to $10^{-4}$, the buoyancy force weakens, leading to a corresponding decrease in the frequencies of the $g$-modes. The $g_{1}$-mode, which has no nodes in the interior of the NS, gradually transitions into an $i$-mode, with its frequency approaching the $i$-mode frequency of a nonstratified star, $f=32.9\,\rm Hz$. This transition is illustrated in Fig. 11 for three representative values of $\delta$.

When $\delta=10^{-2}$, the buoyancy force of the $g_{1}$-mode is much stronger than the restoring force of the $i$-mode, and the eigenfunction clearly exhibits the characteristics of a $g_{1}$-mode. For $\delta=9.4\times 10^{-4}$, the frequency of the $g_{1}$-mode is slightly higher than that of the $i$-mode, and the eigenfunction displays a mixed character, with a dominant contribution from the $i$-mode. At $\delta=10^{-4}$, the frequency of the $g_{1}$-mode is nearly equal to that of the $i$-mode, and the eigenfunction behaves as a pure $i$-mode; compare with Fig. 6. Since the frequency of $i$-modes is in the range of $30$–$40\,\rm Hz$, comparable to $f_{\rm crit}\sim 100\,\rm Hz$, the mode can penetrate into the crust. As the difference between $\Gamma_{1}$ and $\Gamma_{0}$ in the outer core for the TW99 EOS is much larger than $10^{-2}$, the $i$-mode naturally disappears and only leaves imprints near the crust-core interface in the higher-order $g$-modes.

The frequencies of the $g_{2}$ and $g_{3}$ modes drop well below the $i$-mode frequency as $\delta\sim 10^{-4}$, and these modes gradually transition into the core $g_{1}$ and core $g_{2}$ modes, respectively. As shown in Fig. 11, the $g_{2}$ mode becomes nearly confined to the core and exhibits the characteristic structure of a core $g_{1}$ mode. This is consistent with expectations, as the mode frequency of the $g_{2}$-mode, $f=12.7\,\mathrm{Hz}$, is much smaller than the critical frequency $f_{\rm crit}=100\,\mathrm{Hz}$ (for $\lambda\sim R$), and thus the mode is effectively excluded from the crust by the restoring shear forces.

It is worth emphasizing here that compositional $g$-modes are different from thermal $g$-modes that are supported by temperature gradients. For an old NS with a typical temperature of $\sim 10^{7}\,\mathrm{K}$, the frequencies of low-order thermal $g$-modes are of the order of $0.001$–$0.01\,\mathrm{Hz}$, well below $f_{\rm crit}$, and are therefore completely excluded from the solid crust. These modes are confined to the fluid core or the ocean layer, and are referred to as core $g$-modes and ocean $g$-modes in McDermott et al. [109].

V Tidal resonance kinetics and orbital energy transfer

The spectrum of an NS almost completely dictates its internal fluid motion to the linear order (though see Refs. [130, 131, 90, 137, 96, 97] for progress in nonlinearity). In particular, the perturbative flows $\xi^{i}$ inside an inspiraling NS driven by the tidal field of its companion can be approximately modeled as a sum of harmonic oscillators,

[TABLE]

where $q^{\alpha}$ is the amplitude of a specific mode with $\alpha$ denoting the quantum number $\{n,\,\ell,\,m\}$. Hereafter we use the indices $\alpha$ and $\{n,\,\ell,\,m\}$ interchangeably. How these oscillators are forced and excited as the binary orbiting around at a chirping frequency is governed by the Hamiltonian

[TABLE]

Here $H_{\rm orb+gw}$ captures the inspiral motion for point particles, which includes the conservative orbital dynamics and the dissipative part caused by back reaction of GWs. They can be well described by post-Newtonian dynamics via the effective-one-body formalism [24]. The Hamiltonian for the harmonic oscillators is denote by $H_{\rm osc}$, and $H_{\rm tid}$ denotes the tidal work done on the NS through the excitation of its internal modes of oscillation.

In Newtonian gravity, $\xi_{\alpha}^{i}$ satisfies a self-adjoint eigenvalue equation (see, e.g., Refs. [101, 134]),

[TABLE]

where $\mathcal{C}$ is appropriate potential operator. Different modes are orthogonal to each other

[TABLE]

where $\mathcal{A}_{\alpha}^{2}$ is the normalization factor, and overhead bar denotes complex conjugation. The Hamiltonian $H_{\rm osc}$ can be written as [6, 83]

[TABLE]

where $p^{\alpha}$ is the canonical momenta associated with $q_{\alpha}$.

In general relativity, the eigenvalue problem for stellar oscillations is generally non-self-adjoint due to the presence of GW emission, which leads to complex mode frequencies. However, for certain modes such as $g$-modes and interface $i$-modes, the GW damping timescales are extremely long [86, 88]. For these modes, the imaginary part of the frequency can be safely neglected with their frequencies and eigenfunctions essentially unchanged. Assuming a real spectrum leads to treating the eigenvalue problem as Hermitian system when computing mode energies and couplings. By adopting the variational principle used in Refs. [36, 146, 42], and retaining only the conservative contributions, we arrive at

[TABLE]

with $g$ being the determinant of the metric $g_{\mu\nu}$. Here, on the left-hand side, the first term represents the kinetic energy of the perturbed fluid while the terms involving metric perturbations have no Newtonian counterpart and may be broadly interpreted as the kinetic energy of the gravitational field. The right-hand side consists of potential terms with lengthy expressions that are not needed in this paper; we therefore refer the reader to Refs. [36, 42] for more details. By setting the normalization factor

[TABLE]

the corresponding Hamiltonian retains the canonical oscillator forms as given in Eq. 37. We omit the subscript “$\alpha$” on the metric functions for simplicity. For the literature adopting the Cowling approximation (e.g., Ref. [31]), the metric perturbations are omitted completely throughout (not only restricted to the real frequencies). In this limit, the terms associated with the metric perturbation are absent from the normalization factor. Within the relativistic framework, an alternative choice of the normalization factor used in Ref. [94] normalizes eigenfunction by its components in fluid motion while the contributions of metric functions are suppressed since their contributions are considerably less important than the former. In this work, we take the normalization factor in Section V to be

[TABLE]

The Hamiltonian $H_{\rm tid}$ characterizes the gravitational energy absorbed by the NS due to its coupling to the companion’s tidal field. A full treatment requires solving for the dynamical tidal response in the buffer zone surrounding the star [165], which calls for further development along the lines of Refs. [62, 80, 155], and is beyond the scope of this work. Instead, we adopt a simplified approximation following Kuan et al. [94], in which the interaction Hamiltonian takes the form

[TABLE]

where $\delta\bar{\epsilon}$ is the conjugation of the Eulerian perturbation of the energy density, $\delta\epsilon$, and the external tidal potential is given by the Newtonian one [134]

[TABLE]

Here, $M^{\prime}$ is the mass of the companion, $a$ is the orbital separation, $\Phi(t)$ is the orbital phase, and $W_{lm}$ is zero for odd $\ell+m$, and otherwise

[TABLE]

For quadrupolar perturbations, we have

[TABLE]

The Hamiltonian $H_{\rm tid}$ can then be recasted into

[TABLE]

with

[TABLE]

the relativistic analogue of the tidal overlap integral.

In Fig. 12, we present $Q_{n\ell}$ as a function of the NS mass for several representative modes: the $g_{1}$-, $g_{2}$-modes, the $i$-mode in the nonstratified case, the $s_{1}$-, $s_{2}$-modes, and the $f$-mode. The results corresponding to the canonical model is shown in Table 2. Among these, the $f$-mode exhibits the strongest coupling, with $Q_{n\ell}\sim 0.2$ – $0.5$ [150, 140, 101, 83, 94, 31]. In contrast, the $g$-modes possess significantly weaker couplings (3 to 5 orders of magnitude smaller than that of the $f$-mode) that is smaller for higher overtones. The $s_{1}$- and $s_{2}$-modes yield overlap integrals of the order of $10^{-4}$. Unlike the $g$-modes, their values do not show a clear trend of decreasing with increasing radial order. The $s_{1}$-mode has larger tidal overlap than $s_{2}$-mode for small masses and their dominance switches when $M\approx 1.5\,M_{\odot}$. The $i$-mode in the nonstratified case has a tidal overlap comparable to that of the $g_{1}$-mode. Neglecting the metric perturbation terms in the normalization [Section V] and in the tidal overlap integral [Section V] affects the computed value of $Q_{n\ell}$ by no more than $10\%$ across most NS masses. While our simplified treatment of $H_{\rm tid}$ may introduce some quantitative deviations, it is not expected to qualitatively alter our results.

Given the mode frequency, tidal overlap integral, and binary parameters, the evolution of the orbital phase and mode amplitude can be determined by numerically integrating the Hamiltonian equations of motion.

V.1 Energy transfer and estimates of GWs

To have a semi-quantitative estimation on the evolution of modes in the late inspiral phase, we treat the orbit as Newtonian and adopt the quadrupole formula of GW dissipation, which provides a reasonable approximation for our purposes. The corresponding Hamiltonian can be found in Ref. [83]. The orbital separation, $a$, of GW-driven quasi-circular inspirals decreases at the rate

[TABLE]

where the dot means the derivative with respect to time $t$. The orbital angular frequency of the binary,

[TABLE]

increases with time. A resonance between a given mode with $l=|m|=2$ and tidal force occurs when $\omega_{\alpha}\simeq 2\Omega=2\pi f_{\rm gw}$, where $f_{\rm gw}$ is the GW frequency. After evolving the Hamiltonian equations of motion, the energy transferred from the orbit to the modes can be directly computed from the oscillator Hamiltonian $H_{\rm osc}$, which includes both kinetic and potential energy contributions. From the Hamiltonian, one can readily derive the amplitude evolution equation for nonrotating NSs,

[TABLE]

This equation describes a forced harmonic oscillator, and is formally identical to the Newtonian result; see, e.g., Refs. [101, 134, 180].

In Fig. 13, we show the evolution of the mode energy $E_{\rm mode}$ computed from Eq. 37 during the inspiral for a binary NS system with equal mass $M^{\prime}=M=1.35\,M_{\odot}$ in the unit of $E_{\rm b}=GM^{2}/R$. The modes are gradually excited and the amplitude grows as the separation of the binary decreases. Before reaching the resonance, i.e., when $\omega_{\alpha}^{2}\gg(m\Omega)^{2}$, the solution of $q^{\alpha}$ scales as

[TABLE]

and consequently, the mode energy scales as

[TABLE]

The $f_{\rm gw}^{4}$ scaling of the $f$-mode energy is clearly visible. We also see how tidal resonances excite different oscillation modes of the NS at different values of $f_{\rm gw}$. After the orbital frequency goes off resonance with a given oscillation mode, the mode’s energy nearly remains constant with minor modulations. This is because the tidal force continues to act upon the oscillation, although the net energy transfer is negligible. For the stratified NS, the $f$-mode always dominates over the $g$-modes, even though it does not enter into the resonance regime during the inspiral in this model. This result is consistent with the Newtonian study in Ref. [124]. We also show the energy evolution of the $i$-mode for the nonstratified model. Since its frequency is approximately one order of magnitude smaller than the frequency of the $g_{1}$-mode, the $i$-mode has more time to be in resonance with binary’s orbital frequency and hence to accumulate energy. Around the resonance point of the $i$-mode, its mode amplitude is larger than that of $f$-mode. The “post-resonance” energy of the $i$-mode lies between that of the $g_{1}$ and $g_{2}$ modes.

To give a better physical interpretation, we employ approximate methods to quantify the post-resonance amplitude and energy of the modes in the following. The orbital decay timescale can be estimated as

[TABLE]

where $q=M^{\prime}/M$ is the mass ratio of the binary. Since this timescale is much longer than the resonance timescale, we can determine $q^{\alpha}$ with the stationary-phase approximation. In particular, the Green-function solution of Eq. 49 gives the following estimate for the peak amplitude of the excited modes:

[TABLE]

Note that both $m=\pm 2$ modes contribute equally to the tidal response. In Table 2, we present the maximum amplitudes computed using Eq. (53), alongside the exact results obtained from direct numerical evolution of Hamilton’s equations. The discrepancy between the two approaches is less than 5%.

The kinetic energy and the potential energy are equal after the resonance, and the total energy of the mode can be estimated as

[TABLE]

where $\Omega_{0}=(GM/R^{3})^{1/2}$, $R_{12}=R/12\,\rm km$, $M_{1.35}=M/1.35\,M_{\odot}$, and $E_{\rm b}=GM^{2}/R$. The dependence on binary’s component masses in the expression for the mode’s maximal energy enters primarily through the factor $q\left(2/(1+q)\right)^{-5/3}$, which is relatively weak compared to the dependence on the frequency and the tidal overlap integral. We will discuss this below in this section.

Because a fraction of the orbital energy is absorbed by stellar oscillations, the binary’s orbital phase will be shifted. The characteristic energy scale of the binary is

[TABLE]

An estimate of the shift in the orbital phase $\Delta\Phi$ due to the energy transfer $E_{\rm mode}^{\rm max}$ can be estimated as [101, 66]

[TABLE]

where $t_{\rm orb}=2\pi/\Omega$ is the orbital period. In Fig. 14, we plot the relation between mode frequency and the tidal overlap integral for a range of NS masses from $0.4\,M_{\odot}$ up to the maximum mass of $2.08\,M_{\odot}$. As visual representations of the analytic expressions given in Section V.1 and Section V.1, we overlay contour lines of the mode energy $E_{\rm mode}^{\rm max}$ and the GW phase shift $\Delta\Phi$ for an NS with $M=1.35\,M_{\odot}$ and $R=12\,\rm km$ in an equal-mass binary. Note that the $s$-modes do not reach resonance in all models, and the contour lines represent an upper limit. We also find that the mode energy remains below $\sim 10^{47}\,\rm erg$ for the $s$-modes, $g$-modes, and $i$-modes in NSs with masses $M\gtrsim 1.1\,M_{\odot}$. Correspondingly, the GW phase shift induced by resonant excitation of these modes stays below $10^{-3}\,\rm rad$. This extent of dephasing is quite challenging to be observable even with the next generation ground-based GW observatories [136].

So far we have focused on equal-mass binaries. Here, we briefly examine how the maximal mode energy depends on the mass ratio. The mass ratio $q$ enters into $E_{\rm mode}^{\rm max}$ in Section V.1 through the factor $q[2/(1+q)]^{-5/3}$, which equals to unity for equal-mass binaries, reaches a maximum at $q=1.5$, and decreases to zero for large values of $q$. For asymmetric NS–NS binaries, the typical mass ratio is smaller than 1.5, and thus the mode energy deposited through a resonance is comparable to that in equal-mass binaries. However, for NS–BH binaries, the mass ratio is typically much greater than $1.5$ [5]. We therefore expect a significant suppression of the amount of energy transferred into mode energy after resonance. We confirm this expectation in Fig. 15, for a representative NS–BH system with masses $M=1.35\,M_{\odot}$ for the NS and $M^{\prime}=10\,M_{\odot}$ for the BH. Comparing to the equal-mass NS–NS scenario, the maximum mode energy of the resonantly excited mode is reduced by approximately $21\%$. Another key effect is that NS–BH binaries merge at lower merger frequency compared to NS–NS systems. As shown in Fig. 15, the $s$-modes can no longer reach resonance, and the energy of the nonresonantly excited $f$-mode is suppressed by three orders of magnitude relative to the NS–NS case shown in Fig. 13.

V.2 Crust breaking and the von Mises criterion

The resonance-induced oscillations generate stresses in the crust, which may cause it to break. For a relaxed background star, the local maximum strain is given by

[TABLE]

Remember that $q^{\alpha}$ is the amplitude of an oscillation mode and $\Delta s_{ab}$ is the Lagrangian perturbation of the strain tensor; cf. Eq. (22). The prefactor $\sqrt{2}$ accounts for contributions from both $m=2$ and $m=-2$ components. For perturbations with $\ell=|m|=2$, this expression takes the explicit form

[TABLE]

Remember again that $T_{1}$ and $T_{2}$ are the radial and tangential tractions, defined in Eq. (25); $V$ arises in the angular components of the displacement vector, defined in Eq. (12); $\lambda$ is a background metric function; and $\mu$ the shear modulus. We note that our expression for the strain differs from that in Ref.[124] (Eq.(58)), although the source of this discrepancy is not yet understood. For reference, we have made a Mathematica notebook with our derivation available online [2].

In Fig. 16, we illustrate the strain amplitude for different oscillation modes of the canonical model. Since the $f$-mode does not undergo resonance, and the $s_{1}$- and $s_{2}$-modes reach the resonance point very close to merger, we present the results for these high-frequency modes at the amplitude $|q^{\alpha}|$ corresponding to the orbital phase at which the $g_{1}$-mode is in resonance. In contrast, the results for the low-frequency modes ($g_{1}$-, $g_{2}$-, and $i$-modes) are shown at their respective resonance points, with amplitudes evaluated at resonance. Except the $i$-mode, a general feature is that the strain is largest along the direction to the companion ($\theta=90^{\circ}$) compared to other $\theta$ values. Owing to the linear nature of the problem, the strain amplitude at any other orbital phase can be obtained by rescaling with the value of $|q^{\alpha}|$ provided in Fig. 17. The time before merger is related to the frequency of GW as

[TABLE]

Here $t_{\rm merger}$ denotes the merger time, and $M_{c}$ is the chirp mass of the binary, which equals to $1.175\,M_{\odot}$ for an equal-mass binary NS with $M^{\prime}=M=1.35\,M_{\odot}$.

We also observe in Fig. 16 that the $f$-mode’s strain amplitude exhibits several peaks and valleys in the inner crust, and increases nearly monotonically in the outer crust toward the star’s surface. This behavior is driven by the characteristics of the horizontal displacement, which oscillates in the inner crust and grows steadily toward the surface in the outer crust, as shown in Fig. 4.

The crust is broken when the strain amplitude exceeds a critical threshold, i.e., $\sigma>\sigma_{\rm break}$, which is the so-called von Mises criterion. The breaking strain is commonly taken to be $\sigma_{\rm break}=0.1$, based on molecular dynamics simulations of the NS crust [68]. As shown in the contour plot, the outer crust has already been broken at the position where $g_{1}$-mode reaches resonance. If we consider a more conservative estimate of $\sigma_{\rm break}=0.04$, as suggested by Ref. [17], a part of the inner crust will also be broken. To our knowledge, we have thus demonstrated for the first time that the nonresonant excitation of the $f$-mode can break the crust at $f_{\rm gw}\sim 100\text{--}200\,\rm Hz$. A Newtonian study in [124] also indicated that the $f$-mode can break the crust, but it only occurs very close to the merger. This inconsistency may be caused by the simplified EOS and the different strain expression [Eq. (89) in [124]] adopted therein.

To obtain a clearer understanding of crustal breaking induced by the nonresonant $f$-mode, we analyze the contributions of the three terms appearing in Section V.2. We find that the term involving the horizontal traction, $T_{2}$, is almost always dominant—typically 1–2 orders of magnitude larger than the other two terms for the $f$-mode. Consequently, the radial gradient of the horizontal displacement is the primary driver of crust breaking, and the leading-order contribution to the von Mises strain is given by

[TABLE]

Here, the amplitude of the $f$-mode is set to its value at the $g_{1}$-mode resonance. In Fig. 18, we plot the radial profile of $rV^{\prime}$ within the crust for a $1.35\,M_{\odot}$ NS. Combined with the estimate in Section V.2, this analysis indicates that the outer crust is already broken in the vicinity of the $g_{1}$ resonance. Meanwhile, the peak-valley structure in the inner crust region observed in Fig. 16 arises from the radial variation of $V^{\prime}$, as illustrated in Fig. 18. It is worth noting that static tides are generally ineffective at fracturing the crust, as demonstrated in [127, 55]. This is primarily because, under static tidal forces, the solid crust deforms much like a fluid. As a result, the star can experience significant global deformation while the internal shear strain remains small, making static tides inefficient at breaking the crust.

For the $s_{1}$-mode, the strain amplitude peaks in the inner crust; however, the maximum amplitude at the $g_{1}$ resonance is insufficient to break the crust. Based on the mode amplitude evolution shown in Fig. 17, the inner crust is expected to break at GW frequencies of approximately $f_{\rm gw}\sim 400\,\mathrm{Hz}$ and $\sim 600\,\mathrm{Hz}$ for breaking strains of $\sigma_{\rm break}=0.04$ [17] and $\sigma_{\rm break}=0.1$ [68], respectively. In contrast to the $s_{1}$-mode, the $s_{2}$-mode exhibits its maximum strain in the outer crust, with an overall pattern resembling that of the $f$-mode. As shown in the contour plot, at the orbital phase corresponding to the $g_{1}$-mode resonance, only a small region of the outer crust approaches the breaking strain $\sigma_{\rm break}=0.04$. The GW frequency at which the $s_{2}$-mode breaks the inner crust is approximately similar to that of the $s_{1}$-mode.

For the $g_{1}$- and $g_{2}$-modes, the strain distribution is similar, as their eigenfunctions share a comparable structure throughout the star. The strain is clearly larger along the direction of the companion relative to other directions, with the largest values occurring near the crust-core boundary. This is consistent with the fact that the buoyancy force peaks at these locations. The maximum strain of the $g_{1}$-mode reaches approximately $0.025$, which is close to the typical breaking strain. However, we cannot conclude that $g$-modes are incapable of breaking the crust, as their frequencies and tidal overlap integrals are highly sensitive to the nuclear physics near the crust-core interface. A more detailed investigation into the possibility of crust breaking by $g$-modes across different EOSs will be presented in a forthcoming paper.

Finally, for the $i$-mode in the nonstratified case, the strain is largest at the base of the crust and concentrated near the equator. Although the mode energy at resonance is smaller than that of the $g_{1}$-mode, the $i$-mode is more efficient at inducing strain in the crust. It can completely break the crust if we adopt a breaking strain of $\sigma_{\rm break}=0.04$, with the maximum strain approaching $\sigma_{\rm break}=0.1$. However, as previously discussed, this mode does not exist for stratified NSs. The characteristic motion of $i$-mode is encoded in the high overtones of $g$-modes, while these modes are only feebly susceptible to tidal pushing force and can be barely excited to leave observables.

V.3 Implication in EM precursors

Although the energy deposited into oscillation modes generally has a negligible effect on the GW phase evolution except for the $f$-mode, a given mode could produce observable precursor flares before the merger if it manages to break the crust and the energy thus liberated can be efficiently converted into EM radiation. Once a mode grows to the point where the induced stress exceeds the crust’s elastic limit, the crust may break. This shattering process redistributes the mode energy into high-frequency oscillations, which can couple more effectively to the star’s magnetic field and generate radiation. This mechanism was first proposed by Tsang et al. [170] and has since been explored in various studies [170, 94, 95, 160, 161, 115]. In this context, even modes that are energetically subdominant in GW emission may play an important role in generating EM counterparts through crustal breaking and magnetic coupling.

A complete description of the nonlinear dynamics following crust breaking and the subsequent radiative processes in the magnetosphere is highly complex. Since our study focuses on a single EOS and is primarily concerned with establishing the formalism, we restrict ourselves to presenting the energy budget available $E_{\rm crust}$ (the mode energy confined within the crust) and the elastic potential energy $E_{\rm ela}$. These contributions can be expressed as

[TABLE]

where the integrals are evaluated over the crust region only.

In Fig. 19, we present $E_{\rm crust}$ and $E_{\rm ela}$ for low-frequency modes at their respective resonance as functions of NS masses. Both quantities decrease with increasing stellar mass, primarily due to the reduction in the tidal overlap integral. For the $g_{1}$- and $g_{2}$-modes, the elastic potential energy constitutes only about $1\%$ of the mode energy confined in the crust, whereas for the $i$-mode (in the nonstratified case) this fraction increases to approximately $10\%$. For high-frequency modes, not all models reach resonance during inspiral; therefore, we only show results for the canonical model as a function of orbital evolution in Fig. 20. Among these, the mode energy confined in the crust for the $f$-mode remains dominant, with the elastic contribution again being only $\sim 1\%$ of $E_{\rm crust}$. In contrast, for the $s_{1}$- and $s_{2}$-modes, $E_{\rm crust}$ is several orders of magnitude smaller than that of the $f$-mode, and yet $E_{\rm ela}$ can still reach about $10\%$ of $E_{\rm crust}$, similar to the $i$-mode. This feature reflects the stronger shear character of the $s$- and $i$-modes compared to the globally coherent $f$-mode and composition-driven $g$-modes. The enhanced shear motion can be directly seen in the eigenfunctions or inferred from the von Mises stress criterion shown in Fig. 16.

Once the crust breaks, the energy stored in various oscillation modes can potentially power EM precursor emissions. Based on our canonical model, we propose a plausible scenario outlining this process. As the binary orbit decays, low-frequency $g$-modes are resonantly excited first, while high-frequency modes such as the $f$- and $s$-modes are nonresonantly driven, with the $f$-mode contributing to the dominant share of energy. At the resonance of the $g_{2}$-mode, the crust is not significantly overstrained by any of the excited modes. The oscillation energy remains confined within the star, with $E_{\rm crust}\sim 10^{44}\,\mathrm{erg}$. We do not consider the $i$-mode here, as it is expected to mix with $g$-modes in a stratified NS, which is a more realistic scenario during the late inspiral phase [8].

As the orbit further shrinks, the $g_{1}$-mode reaches resonance. By this stage, the outer crust has already broken because of the accumulated stresses from the nonresonantly excited $f$-mode and $s_{2}$-mode, as indicated in Fig. 16. This could initiate the onset of EM precursor emission. Shortly afterward, the system evolves to a GW frequency of $f_{\rm gw}\sim 300\,\mathrm{Hz}$ within a timescale of $\sim 0.1\,\mathrm{s}$. At this point, the crust is likely fully broken, according to Figs. 17 and 16. The total mode energy confined in the crust, dominated by contributions from the $f$- and $g_{1}$-modes, reaches $\sim 10^{46}\,\mathrm{erg}$ and becomes available to power EM emissions. If the energy release occurs on a timescale of $0.01$–$0.1\,\mathrm{s}$, the corresponding luminosity could reach up to $10^{45}$–$10^{46}\,\mathrm{erg/s}$ assuming the 1% energy conversion efficiency [115]. It may be worth mentioning that the typical luminosity of observed precursors is $10^{45\text{--}46}\,\mathrm{erg/s}$, compatible with our proposed scenario.

VI Summary and discussions

In this paper, we investigated nonradial oscillations of NSs with an elastic crust self-consistently modeled in the EOS, and a special attention has been paid of their tidal resonances during the late stage of binary inspiral. The QNMs of a spherically symmetric background were computed using a fully relativistic perturbation framework that incorporates both the elasticity of the solid crust and composition stratification. We modeled the dynamical tidal response of the various modes as driven harmonic oscillators in an external tidal field, and examined the associated energy transfer, resulting GWs dephasing, and potential crust breaking.

In the low-frequency band, the polar-parity oscillations include the $i$-mode associated with the solid–liquid transition at the crust–core boundary, and composition $g$-modes. For the specific EOS considered here, we find that the $i$-mode exists only in nonstratified NSs with frequencies of the order of several tens of Hz, regardless of the NS masses. In reality, composition stratification is expected to be present, which smears out the $i$-mode and withdraws the clear distinction between $i$- and the compositional $g$-modes. Instead, a number of modes with mixed gravity–interfacial character are revealed: the eigenfunctions of the leading-order $g$-modes are predominantly shaped by buoyancy from the composition gradient, while the feature of $i$-mode – a cusp in the radial motion and a discontinuous tangential motion at the crust-core boundary – gradually manifests in overtones (Section IV.2).

To better understand how the transition in the motion characteristics relates to the competition between the buoyancy and shear stress, we adopted a simple phenomenological stratification model, $\Gamma_{1}=\Gamma_{0}(1+\delta)$, whereby the buoyancy is smoothly controlled by varying the parameter $\delta$. (The buoyancy is completely switched off when $\delta=0$.) As $\delta\to 0$, the global $g_{1}$-mode gradually transitions into an $i$-mode, while higher-order $g$-modes become more and more confined into the core (Fig. 11). In the high-frequency regime, we analyzed the $f$-mode and $s$-modes. Both exhibit significant amplitude in the crust, with the latter being almost entirely confined within the crust.

The crust not only supports new families of oscillations, such as shear $s$-modes, but also exerts shear stresses that significantly influence mode propagation in this region. We established a penetration criterion by comparing the mode’s transverse acceleration with the restoring force from crustal shear stresses [see Eqs. (30) and (31)]. We found that the crustal $g$-mode identified in Ref. [31], associated with neutron drip at $\sim 3$–$4\times 10^{11}\,\rm g/cm^{3}$, is absent from the QNM spectrum due to strong shear forces. A similar disappearance of density-discontinuity $g$-modes [41], associated with phase transitions between nuclear species, was previously reported by Krüger et al. [86]. This can also be accommodated by our criterion. In contrast, the $i$-mode in the nonstratified case, core composition-gradient $g$-modes, and the $f$-mode can all penetrate the crust due to their large transverse momentum. In particular, the noticeable permeation of the composition-gradient $g$-modes into the crust essentially removes the concern raised by Neill et al. [118] (and its series of works) that these $g$-modes of oscillations are unable to break the crust. Understanding such mode penetration is essential, especially when evaluating the shear strain in the crust and the possibility of crust breaking triggered by mode excitation.

For the dynamical tidal excitation of various modes, we found that the energy of the $f$-mode is always dominant although the exact resonance is not reached. The transferred energy into $s$-modes and $g$-modes is smaller than $\sim 1\%$ of $f$-mode for typical NSs since their tidal overlapping is by 3–4 orders of magnitude smaller than that of $f$-mode. Correspondingly, the GW phase shift induced by resonant excitation of these modes stays well below $10^{-3}\,\rm rad$ for typical NSs.

We also analyzed the strain induced in the crust by various oscillation modes and found that nonresonant $f$- and $s$-modes can trigger crust breaking already around the first gravity-mode resonance at $f_{\rm gw}\sim 200\,\rm Hz$ for an equal-mass binary with component masses $1.35\,M_{\odot}$. Interestingly, the $f$-mode and $s_{2}$-mode produce the largest strain in the outer crust with similar spatial patterns, while the $s_{1}$-mode peaks in the inner crust. Our results show slight tension with fully Newtonian studies [124, 125], which also report crust breaking from nonresonant $f$- and $s$-modes excitation, but only very close to merger. This discrepancy may stem from different strain expressions: Eq. (68) in [124] versus our Section V.2. Moreover, our calculations of the modes are relativistic, which is more accurate. The $g_{1}$-mode resonance concentrates strain at the crust–core boundary but yields only marginal stress for breaking. These results suggest that resonant and nonresonant mode excitations can break the crust, and potentially contribute to pre-merger energy transfer into the magnetosphere via crust shattering, offering a potential mechanism for EM precursor signals.

For asymmetric NS–NS binaries, our conclusions remain largely unaffected, as the typical mass ratio lies within the range of $\sim 1\text{--}2$, and the mode energy given in Section V.1 should be on the same order. In contrast, for NS–BH binaries with the BH mass of about $10\,M_{\odot}$, the mass ratio is much higher. This leads to the suppression of the mode energy after resonance by approximately one order of magnitude. Moreover, due to the substantially lower merger frequency compared to NS–NS binaries, the $s$-modes may no longer be resonantly excited, and the amplitude of the nonresonantly excited $f$-modes is reduced by several orders of magnitude. As a result, crust breaking is less likely to occur in NS–BH binaries, or if it does occur, it happens very close to the merger.

Our work represents a step toward more realistic modeling of dynamical tides in binary NSs. In this study, we mainly establish the formalism using a single EOS, and there remains a lot to be done. For example, as is well known, the frequencies and tidal overlap integrals of $g$-, $i$-, and $s$-modes depend sensitively on nuclear parameters near saturation density. A systematic investigation of these dependencies is crucial for clarifying the disappearance of the $i$-mode, verifying the tensions in the predicted frequencies and tidal overlaps of $i$-modes compared to Refs. [170, 123, 186], assessing the role of $g$-mode resonances, and evaluating crust breaking triggered by both resonant and nonresonant excitations. It is worthwhile to note that, once the crustal strain exceeds the elastic limit, plastic flow may occur, converting oscillation energy into thermal one [20, 82, 123, 186]. It will be interesting to examine how different oscillation modes contribute to such heating.

It is also important to note that further advances are needed in both the microphysical description of matter and the relativistic formulation of the problem. More specifically, for matter fields, superfluidity in the outer core is expected to increase the frequencies of the gravity modes [78, 125, 180], potentially modifying the mode spectrum and the crustal breaking behavior predicted in our work. Going beyond $npe\mu$ matter studied here, phase transitions to exotic degrees of freedom — such as hyperons or deconfined quarks — may lead to new phase structures, accompanied by additional $g$- and $i$-modes [167, 182, 111, 30, 179]. Recent studies suggest that $i$-mode resonances associated with first-order phase transitions can produce significant phase shifts in GWs, provided that such phase transitions do exist. In addition, bulk viscosity may be important for exotic phases of matter and could potentially result in appreciable tidal heating during resonances [52, 53, 142]. Regarding the modeling of dynamical tides, our treatment of the tidal coupling is based on a relativistic extension of the Newtonian formalism. It would be interesting to repeat our calculations in a full self-consistent relativistic approach (e.g., Refs. [155, 62, 61, 130, 110, 60] for recent progress). Finally, nonlinear effects induced by tidal resonances may also lead to prolonged mode excitations, amplified GW phase shifts [97, 96, 131], and induce differential fluid flows within the NS interior [137, 90]. We will investigate these aspects in future work.

Acknowledgements

We thank the anonymous referee for comments and suggestions. Y. G. thanks Zexin Hu, Jan Steinhoff, and Yanbei Chen for useful discussions. H. O. S. acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG) - Project No.: 386119226. Numerical computations were performed on the Sakura cluster at the Max Planck Computing and Data Facility. This work was in part supported by Grant-in-Aid for Scientific Research (grant No. 23H04900) of Japanese MEXT/JSPS.

Appendix A Polar perturbation for multi-layer NSs

In this Appendix, we present the polar-sector perturbation equations for multi-layer NSs and numerical procedures to solve the eigenvalue problem following Refs. [42, 86, 37, 107]. A Mathematica notebook that derives the perturbation equations is available at [2].

A.1 Master equations

A.1.1 Perturbation equations for perfect fluid

For the perturbations of perfect fluid, we use the form of the perturbation equations derived by Detweiler and Lindblom [37, 107],

[TABLE]

where a prime indicates a derivative with respect to $r$, and $n_{\ell}=(\ell+2)(\ell-1)/2$ for convenience. Here $X$ is related to the Lagrangian variation of the pressure $\Delta p$ as $X=-e^{\nu/2}\Delta p/r^{\ell}$, and the fluid perturbation variable $V$ can be determined by

[TABLE]

The variable $H_{0}$ can be obtained from the following algebraic equation deriving from the $rr$ and $r\theta$ components of the perturbed Einstein equations:

[TABLE]

The metric perturbation function $H_{2}$ satisfies

[TABLE]

These equations have been extensively employed in the analysis of high-frequency oscillation modes, such as the $f$, $p$, and $w$ modes (see [84] for a review).

However, the formalism is problematic for low-frequency modes (e.g., higher-order $g$-modes and $i$-modes). This is purely numerical and mainly stems from the Eq. 67 to calculate $V$, where numerical cancellation is inevitable and hence leads to inaccurate solutions (see Krüger et al. [86] for details). Krüger et al. [86] pointed out that using $V$ instead of $X$ as a “fundamental variable” can solve this problem; the perturbation equations are essentially the same, except for changing Eq. 66 into

[TABLE]

As a summary, in the fluid core, we use Eqs. 63, 64, 65, 66, 67, A.1.1 and 69 if the oscillation frequency $\omega\geq 0.01/M$. Otherwise, we use Eq. 70 to replace Eq. 66. For the fluid ocean (if added), the systematic equations Eqs. 63, 64, 65, 66, 67, A.1.1 and 69 are applied.

A.1.2 Perturbation equations for elastic solid

The full set of perturbation equations for elastic solid can be expressed as [86]222The perturbation equations presented in Ref. [86] miss a factor of 1/2 in the shear modulus $\mu$, which has been also corrected in Ref. [88].

[TABLE]

In addition to the foregoing six ordinary differential equations, we have three algebraic relations

[TABLE]

The first algebraic equation can be used to determine $H_{2}$, while the last two algebraic equations return $T_{1}$ and $X$.

A.2 Joint and boundary conditions

To integrate the ODEs and obtain the interior solutions for a given oscillation angular frequency $\omega$ and EOS, we require matching conditions at the solid-liquid interface, as well as the boundary conditions at both the center and the surface of the NS. The matching conditions at the solid-liquid interface arise from the requirement that the true degrees of freedom of spacetime remain continuous. We use $[Q_{r}]$ to denote the jump of a quantity $Q$ across the interface,

[TABLE]

The matching conditions at the solid-liquid interface are

[TABLE]

Note that the horizontal displacement $V$ is not necessarily continuous across the fluid-solid interface. The equations of motion allow us to derive more matching conditions, though not independent. Eq. 77 indicates that $H_{2}$ is not continuous across the interface, and the jump of $H_{2}$ is given by

[TABLE]

The solutions must remain regular at the center of the star, necessitating a Taylor expansion of the perturbation variables near the origin to approximate the ODE solutions. It has been shown in Ref. [37] that not all the quantities are independent. Once $K(r=0)$ and $W(r=0)$ are specified, all other coefficients in the Taylor expansion can be uniquely determined (see details in Ref. [89, pp. 27–30, 59–60]). The surface of the star is a matter-vacuum interface. If a two-layer NS with a solid surface is considered, the boundary conditions at the surface are

[TABLE]

where the latter one is derived from the algebraic relation Eq. 78 by imposing $[H_{0}]_{R}=[H_{1}]_{R}=[K]_{R}$. While for a three-layer star with a fluid ocean, the boundary condition at the surface is simply

[TABLE]

Since the systematic equations are linear, one can still choose one overall normalization for the solutions. The boundary and joint conditions of the perturbation equations for the multi-layer NS, as well as the overall normalization, are summarized in Fig. 21. It should be noted that the QNM calculations are performed with the normalization $W(R)=1$, whereas the results reported in the main text are rescaled using $\mathcal{A}_{\alpha}^{2}=MR^{\ell}$ [cf. Section V].

A.3 Numerical procedures

A.3.1 The interior solution

Our numerical procedures to solve the perturbation equations for a multi-layer NS is similar with Lin et al. [106] and Krüger et al. [86]. Considering an NS with $n$ layers, where the interfaces including the center and the surface of the NS, are located at the radii $0=R_{0}<R_{1}<R_{2}<\cdots<R_{n-1}<R_{n}=R$. The dynamical equations that govern the perturbations at each layer depends on the nature of the matter field at this layer. For $i$-th layer, the perturbation equations can be written as

[TABLE]

Here ${}^{(i)}\mathbf{Y}=\left(y_{1},\ldots,y_{k_{i}}\right)$ is the vector of perturbation variables and $k_{i}$ is the number of variables in the $i$-th layer to be integrated, ${}^{(i)}\mathbf{M}={}^{(i)}\mathbf{M}(\texttt{star},r,\ell,\omega)$ is a $k_{i}\times k_{i}$ matrix of coefficients, with star represents the background fields.

To obtain the general solution in layer $i$, we need to integrate the perturbation equations $x_{i}$ times starting from the boundary, where $x_{i}$ is the number of linearly independent solutions at each layer. The general solution in layer $i$ is the linear combination of these solutions

[TABLE]

where the coefficients $c_{i,j}$ is determined by the boundary and joint conditions, with $i$ denoting the layer and $j$ counting the different solutions in the layer.

For a two-layer NS, as illustrated in Fig. 21, the fluid interior admits two independent solutions, depending on the choices of $(K(0),W(0))$. In contrast, within the solid crust, the tangential traction at the crust-core interface vanishes $(T_{2}=0)$, reducing the number of independent solutions from six to five. Consequently, there are seven independent solutions in total, demanding seven conditions to determine the coefficients $c_{i,j}$. These conditions comprise four continuity conditions at the crust-core interface: $\left[H_{0}\right]_{R_{cc}}=\left[H_{1}\right]_{R_{cc}}=[K]_{R_{cc}}=[W]_{R_{cc}}=0$, along with three boundary conditions at the NS surface: $T_{2}(R)=0,\ W(R)=1,\ \text{and}\ R^{2}e^{-\nu(R)/2}X(R)+T_{1}(R)=0$.

For a three-layer model, the approach to integrating the perturbation equations in the fluid core and solid crust remains unchanged. However, the presence of a fluid ocean introduces an additional set of perturbation equations. To handle this, we integrate the system of ODEs from the surface to the crust-ocean interface, imposing the condition $X(R)=0$, which yields three independent solutions in the fluid ocean. Consequently, we obtain a total of ten independent solutions. Then we have ten independent solutions in total. Except for the four joint conditions at the crust-core interface, we have four joint conditions at the crust-ocean interface: $\left[H_{0}\right]_{R_{\mathrm{co}}}=\left[H_{1}\right]_{R_{\mathrm{co}}}=[K]_{R_{\mathrm{co}}}=[\mathrm{W}]_{R_{\mathrm{co}}}$ along with one constraint $T_{2}(R_{\rm co})=0$. Furthermore, at the surface, we impose additional one normalization $W(R)=1$. With these conditions in place, the coefficients $c_{i,j}$ can be fully determined.

A.3.2 The exterior solution

Outside of the star, the matter field vanish and the system of perturbation equations reduces to a second-order system of equations for metric perturbation $H_{1}$ and $K$. These equations can be transformed into a standard Zerilli equation [39, 107], which allows two independent solutions. Far away from the star, one of the solutions can be identified as an outgoing wave, whereas the other solution describes an ingoing wave. For a given oscillation frequency $\omega$, the solution for the interior of the NS leads to a mixture of ingoing and outgoing waves. Physically, it means that a general frequency $\omega$ will not corresponds to a resonance of the stellar model, the star is forced to oscillate at this frequency, driving by the ingoing gravitational waves. The QNMs are those special frequencies at which the ingoing wave vanishes, and the star oscillates “freely”.

We follow the phase-amplitude method [47, 12] to solve the exterior problem and identify the QNMs. Instead of solving the Zerilli function $Z$, they introduce a new independent variable $\Psi$

[TABLE]

to transform the linear wave equation to a non-linear differential equation with a slowly varying variable $q$ respect to $r$. The physically acceptable solution for the interior corresponds to a wave at spatial infinity,

[TABLE]

where $A_{\rm in}$ and $A_{\rm out}$ are the amplitude of the ingoing and outgoing waves, respectively. A QNM corresponds to a solution possessing no incoming component (i.e., with $A_{\rm in}=0$).