Dynamical tides in neutron stars: The impact of the crust
MMon. Not. R. Astron. Soc. , 000–000 (0000) Printed 18 December 2020 (MN L A TEX style file v2.2)
Dynamical tides in neutron stars: The impact of the crust
A. Passamonti , N. Andersson and P. Pnigouras , , Via Greve 10, 00146, Roma, Italy School of Mathematics and STAG Research Centre, University of Southampton, Southampton SO17 1BJ, UK Dipartimento di Fisica, “Sapienza” Universit`a di Roma & Sezione INFN Roma1, Piazzale Aldo Moro 2, 00185 Roma, Italy Department of Physics, Aristotle University of Thessaloniki, 54124 Thessaloniki, Greece
18 December 2020
ABSTRACT
We consider the dynamical tidal response of a neutron star in an inspiralling binary,focussing on the impact of the star’s elastic crust. Within the context of Newtoniangravity, we add the elastic aspects to the theoretical formulation of the problem andquantify the dynamical excitation of different classes of oscillation modes. The resultsdemonstrate the expectation that the fundamental mode dominates the tidal responseand show how the usual tidal deformability (and the Love number) emerge in thestatic limit. In addition, we consider to what extent the different modes may be ex-cited to a level where the breaking strain of the crust would be exceeded (locally).The results show that the fundamental mode may fracture the crust during the latestages of inspiral. This is also the case for the first gravity mode, which reaches thebreaking threshold in strongly stratified stars. In our models with a fluid ocean, inter-face modes associated with the crust-ocean transition may also induce crust fracture.If this happens it does so earlier in the inspiral, at a lower orbital frequency.
Key words: stars: neutron, neutron star mergers, gravitational waves
The breakthrough detections of gravitational waves from binary neutron star inspirals (Abbott et al 2017a, 2018, 2019) haveled to renewed focus on the elusive neutron star equation of state. The problem has a number of complicating aspects—bothrelating to the observational data and the theoretical underpinning—but the essential question is quite simple: To what extentcan we use observations to constrain the state and composition of matter under the extreme conditions that neutron starsrepresent?Much of the recent focus has been on the neutron star tidal deformability, essentially the extent to which the tidalinteraction with a binary companion deforms the neutron star fluid. This is a useful measure as it can be extracted from(or at least, constrained by) the gravitational-wave signal (Flanagan and Hinderer 2008; Hinderer et al 2010). Notably, thecelebrated GW170817 event has led to a constraint on a suitable weighted average tidal deformability, corresponding (roughly)to a neutron star radius in the range 10-13 km (Abbott et al 2018) (the result is somewhat model dependent). Morever, asthis radius range agrees well with the constraints obtained from the NICER observations of PSR J0030+0451 (Miller et al2019; Riley et al 2019) a consistent picture is beginning to emerge.The deformability (often expressed in terms of the dimensionless Love number, k l ) represents the static contribution to theneutron star’s tidal response. In addition, there is a dynamical tide. This is traditionally represented by the excitation of thedifferent oscillation modes of the star. The resonance problem was first considered some time ago (Lai 1994; Reisenegger andGoldreich 1994; Kokkotas and Schaefer 1995; Andersson and Ho 2018), but the issue is back in focus following the suggestionthat the (fundamental) f-mode of the star may be excited to a relevant level, even though it may not reach resonance duringthe inspiral (Hinderer et al 2016; Steinhoff et al 2016; Andersson and Pnigouras 2020b). The associated effect on the inspiralsignal is weak, but its inclusion has been demonstrated to improve waveform models.When we turn to dynamical features of the tide, we need to be mindful of the fact that a neutron star interior is a littlebit more complicated than a prescribed pressure-density relation. Nuclei in the lower density region are expected to freezeto form the so-called crust, neutrons and protons likely form superfluid/superconducting condensates at high densities, there © a r X i v : . [ a s t r o - ph . H E ] D ec Passamonti, Andersson and Pnigouras may be phase transitions to states with net strangeness (involving hyperons or deconfined quarks) and so on. These issuesare important because the tidal response can be expressed as a sum over the (presumably complete set of) stellar oscillationmodes (Lai 1994; Reisenegger and Goldreich 1994; Kokkotas and Schaefer 1995) and the additional features of the interiorphysics may bring new modes into play and shift existing ones. On the one hand, this makes the problem more complicated.On the other hand, it raises the question of whether the additional features may be observable. Even if this is not expected,one should take the care to demonstrate it (as there may be surprises!). A useful recent step in this direction (Anderssonand Pnigouras 2020a) demonstrates that variations in the composition of the neutron star core affects the dynamical tideat the few percent level. This discussion also provides the (conceptually important) link between the static tide and thestar’s oscillation modes. Another important step was taken by Yu and Weinberg (2017) who considered the impact of moderesonances in a superfluid neutron star core, while Poisson (2020) recently revisited the problem of (slowly) rotating stars,demonstrating that this issue also requires further thought.How do we do better? In principle, we know that the problem requires a fully relativistic analysis (one of the mostdetailed seismology models to date was developed by Kr¨uger, Andersson and Ho (2015)). However, the framework required toquantify the role of tidal resonances in a relativistic star has not yet been developed. There are technical issues to consider,including the fact that the modes of the system—which is now dissipative as gravitational waves are emitted—can not be(formally) complete. Quantitatively, this may not be important. Conceptually, it is a hurdle. While we ponder this issue, wemay progress using phenomenological Newtonian models. In this direction, perhaps the most “complete” effort is provided bythe results of Passamonti and Andersson (2012) (see also Andersson, Haskell and Samuelsson (2011)), which represent a starwith an elastic crust penetrated by superfluid neutrons as well as a superfluid core (taking into account the entrainment effectin both regions) and an ocean. A reasonable aim then would be to consider the tidal response of such a star. This involvestwo steps. First, we need to develop the mode-coupling framework for a star with an elastic crust and ocean, and ensure thatthe mode calculation is sufficiently precise that we can quantify the (presumably) weak tidal excitation of (say) the crustalshear modes, the ocean g-modes and the interface modes (McDermott et al 1985; McDermott, van Horn and Hansen 1988).The second step involves adding the superfluid components. With this paper, we take the first of these steps. We also considerthe issue of tidally induced crust fracturing (Tsang et al 2012; Pan et al 2020)—to what extent resonant oscillation modesmay reach the amplitude required to exceed the breaking strain at some point in the crust.
In order to set the stage for the analysis, let us remind ourselves of the context. The tide raised by a binary companion(treated as a point particle, which should be a good enough approximation for our purposes) induces a linear response in theprimary. If we want to quantify the fluid response to this external agent—in essence, model the dynamical tide—we need tosolve the linearised fluid equations in Newtonian gravity. As we are dealing with an (at least partly) elastic body, it is naturalto do this in the framework of Lagrangian perturbation theory (Friedman and Schutz 1978). In fact, many of the formalitieshave already been dealt with by Andersson, Haskell and Samuelsson (2011) and Passamonti and Andersson (2012).
Assuming that the star is non-rotating, which makes sense on astrophysical grounds as the two partners in a binary wouldhave had plenty of time to spin down before the system enters the sensitivity band of a ground-based gravitational-wavedetector (although, see the recent analysis of Poisson (2020) for effects associated with spinning stars), we first of all have theperturbed continuity equation ∂ t (cid:16) ∆ ρ + ρ ∇ i ξ i (cid:17) = 0 , (1)where ξ i is the displacement vector associated with the Lagrangian perturbation (noting that, following Friedman and Schutz(1978) we express vector components in a coordinate basis)∆ = δ + L ξ (2)(with δ the corresponding Eulerian variation and L ξ the Lie derivative along ξ i ), such that the perturbed velocity is given by∆ v i = δv i = ∂ t ξ i . (3)Provided the background is in hydrostatic equilibrium, the perturbed Euler equation is ∂ t ξ i + 1 ρ ∇ i δp − ρ δρ ∇ i p + ∇ i δ Φ = −∇ i χ (4)and we also have the Poisson equation for the gravitational potential ∇ δ Φ = 4 πGδρ , (5) © , 000–000 mpact of the crust while the tidal potential, χ , due to the presence of the binary partner (which generates the fluid perturbation), is a solutionto ∇ χ = 0.As we progress, it is useful to note that (1) leads to∆ ρ + ρ ∇ i ξ i = 0 = ⇒ δρ = −∇ i (cid:16) ρξ i (cid:17) , (6)which allows us to write the Poisson equation as ∇ i (cid:16) g ij ∇ j δ Φ + 4 πGρξ i (cid:17) = 0 . (7)This will be useful later, as we can integrate to get ∇ i δ Φ = − πGρξ i + S i (8)for some vector S i such that ∇ i S i = 0. It is easy to see that, if ρ → S i = ∂ i δ Φ , at r = R , (9)and we know that this should not vanish. Moreover, we know that the solution has to match to the vacuum exterior, whichmeans that (again, as long as the density vanishes as r → R ) ddr δ Φ l + l + 1 r δ Φ l = 0 at r = R , (10)where we have expanded δ Φ in spherical harmonics and focussed on a single l multipole. That is, we must haveˆ n i S i = ∂ r δ Φ = ( ∂ r δ Φ l ) Y ml = − l + 1 r δ Φ at r = R (11)where ˆ n i is the normal to the star’s spherical surface (i.e. parallel to the radial basis vector). We will have reason to recallthis result later. We now want to express the driven tidal response of the star—the solution to (4)—in terms of a set of normal modes (Lai 1994;Reisenegger and Goldreich 1994; Kokkotas and Schaefer 1995; Andersson and Pnigouras 2020a), corresponding to solutions ξ in (with n labelling the modes). Letting the (real) mode frequency be ω n , each individual mode then satisfies the homogeneousversion of (4), with χ = 0.In order to provide a complete picture, it is useful to first consider a barotropic model for which we have (see Andersson,Comer and Grosart (2004) for a similar analysis) δp = nδµ = ρδ ˜ µ , (12)where µ = m B ˜ µ is the chemical potential and ρ = m B n , with m B the mass of each baryon and n the baryon number density.This means that. for a background in hydrostatic equilibrium, we have1 ρ ∇ i δp − ρ δρ ∇ i p = ∇ i δ ˜ µ . (13)As a result, we have ∂ t ξ i + ∇ i ( δ Φ + δ ˜ µ ) = 0 , (14)or, if we take the baryon number density to be the primary matter variable; ∂ t ξ i + ∇ i (cid:20) δ Φ + (cid:18) ∂ ˜ µ∂n (cid:19) δn (cid:21) = ∂ t ξ i + ∇ i (cid:20) δ Φ − (cid:18) ∂ ˜ µ∂n (cid:19) ∇ j (cid:16) nξ j (cid:17)(cid:21) = 0 . (15)Following Friedman and Schutz (1978) we now write the Euler equation as − ω n Aξ in + Cξ in = 0 , (16)where A = ρ (17)and Cξ in = ρ ∇ i (cid:20) δ Φ − (cid:18) ∂ ˜ µ∂n (cid:19) ∇ j (cid:16) nξ j (cid:17)(cid:21) (18)and introduce the inner product (cid:104) η, ρξ (cid:105) = (cid:90) ρη ∗ i ξ i dV , (19) Later, when we consider the contribution from the elastic crust, we assume that the crust is covered by a shallow fluid ocean. Thismakes sense on physical grounds, which is convenient as the surface boundary conditions then remain those of a perturbed fluid star. ©000
Assuming that the star is non-rotating, which makes sense on astrophysical grounds as the two partners in a binary wouldhave had plenty of time to spin down before the system enters the sensitivity band of a ground-based gravitational-wavedetector (although, see the recent analysis of Poisson (2020) for effects associated with spinning stars), we first of all have theperturbed continuity equation ∂ t (cid:16) ∆ ρ + ρ ∇ i ξ i (cid:17) = 0 , (1)where ξ i is the displacement vector associated with the Lagrangian perturbation (noting that, following Friedman and Schutz(1978) we express vector components in a coordinate basis)∆ = δ + L ξ (2)(with δ the corresponding Eulerian variation and L ξ the Lie derivative along ξ i ), such that the perturbed velocity is given by∆ v i = δv i = ∂ t ξ i . (3)Provided the background is in hydrostatic equilibrium, the perturbed Euler equation is ∂ t ξ i + 1 ρ ∇ i δp − ρ δρ ∇ i p + ∇ i δ Φ = −∇ i χ (4)and we also have the Poisson equation for the gravitational potential ∇ δ Φ = 4 πGδρ , (5) © , 000–000 mpact of the crust while the tidal potential, χ , due to the presence of the binary partner (which generates the fluid perturbation), is a solutionto ∇ χ = 0.As we progress, it is useful to note that (1) leads to∆ ρ + ρ ∇ i ξ i = 0 = ⇒ δρ = −∇ i (cid:16) ρξ i (cid:17) , (6)which allows us to write the Poisson equation as ∇ i (cid:16) g ij ∇ j δ Φ + 4 πGρξ i (cid:17) = 0 . (7)This will be useful later, as we can integrate to get ∇ i δ Φ = − πGρξ i + S i (8)for some vector S i such that ∇ i S i = 0. It is easy to see that, if ρ → S i = ∂ i δ Φ , at r = R , (9)and we know that this should not vanish. Moreover, we know that the solution has to match to the vacuum exterior, whichmeans that (again, as long as the density vanishes as r → R ) ddr δ Φ l + l + 1 r δ Φ l = 0 at r = R , (10)where we have expanded δ Φ in spherical harmonics and focussed on a single l multipole. That is, we must haveˆ n i S i = ∂ r δ Φ = ( ∂ r δ Φ l ) Y ml = − l + 1 r δ Φ at r = R (11)where ˆ n i is the normal to the star’s spherical surface (i.e. parallel to the radial basis vector). We will have reason to recallthis result later. We now want to express the driven tidal response of the star—the solution to (4)—in terms of a set of normal modes (Lai 1994;Reisenegger and Goldreich 1994; Kokkotas and Schaefer 1995; Andersson and Pnigouras 2020a), corresponding to solutions ξ in (with n labelling the modes). Letting the (real) mode frequency be ω n , each individual mode then satisfies the homogeneousversion of (4), with χ = 0.In order to provide a complete picture, it is useful to first consider a barotropic model for which we have (see Andersson,Comer and Grosart (2004) for a similar analysis) δp = nδµ = ρδ ˜ µ , (12)where µ = m B ˜ µ is the chemical potential and ρ = m B n , with m B the mass of each baryon and n the baryon number density.This means that. for a background in hydrostatic equilibrium, we have1 ρ ∇ i δp − ρ δρ ∇ i p = ∇ i δ ˜ µ . (13)As a result, we have ∂ t ξ i + ∇ i ( δ Φ + δ ˜ µ ) = 0 , (14)or, if we take the baryon number density to be the primary matter variable; ∂ t ξ i + ∇ i (cid:20) δ Φ + (cid:18) ∂ ˜ µ∂n (cid:19) δn (cid:21) = ∂ t ξ i + ∇ i (cid:20) δ Φ − (cid:18) ∂ ˜ µ∂n (cid:19) ∇ j (cid:16) nξ j (cid:17)(cid:21) = 0 . (15)Following Friedman and Schutz (1978) we now write the Euler equation as − ω n Aξ in + Cξ in = 0 , (16)where A = ρ (17)and Cξ in = ρ ∇ i (cid:20) δ Φ − (cid:18) ∂ ˜ µ∂n (cid:19) ∇ j (cid:16) nξ j (cid:17)(cid:21) (18)and introduce the inner product (cid:104) η, ρξ (cid:105) = (cid:90) ρη ∗ i ξ i dV , (19) Later, when we consider the contribution from the elastic crust, we assume that the crust is covered by a shallow fluid ocean. Thismakes sense on physical grounds, which is convenient as the surface boundary conditions then remain those of a perturbed fluid star. ©000 , 000–000 Passamonti, Andersson and Pnigouras where the asterisk indicates the complex conjugate and η i is another solution to the perturbation equations. It is then fairlyeasy to show that mode solutions are orthogonal if the mode frequencies are real. As we will need to extend the proof of thisto the elastic case, let us have a look at the argument—taking the opportunity to keep careful track of “surface terms” thatcome into play if we allow internal phase transitions.First of all, it is obvious from the definition of the inner product that (cid:104) η, Aξ (cid:105) = (cid:104) ξ, Aη (cid:105) ∗ . (20)The argument for the C operator is a bit more involved. We need to use integration by parts to show that (cid:90) nη i ∗ ∇ i (cid:20)(cid:18) ∂ ˜ µ∂n (cid:19) ∇ j ( nξ j ) (cid:21) dV = (cid:90) n (cid:18) ∂ ˜ µ∂n (cid:19) (cid:104) η i ∗ ∇ j ( nξ j ) − ξ i ∇ j ( nη j ∗ ) (cid:105) ˆ n i dS + (cid:90) nξ j ∇ j (cid:20)(cid:18) ∂ ˜ µ∂n (cid:19) ∇ i ( nη i ∗ ) (cid:21) dV = (cid:90) nξ i ∇ i (cid:20)(cid:18) ∂ ˜ µ∂n (cid:19) ∇ j ( nη j ∗ ) (cid:21) dV , (21)where ˆ n i is the (outwards pointing) normal to the surface of the volume we are integrating over. The last equality holds ifcan ignore the contribution from the surface term. This would certainly be the case if we integrate over the entire star, aslong as n → r = ¯ R , say, we need the surface term n (cid:18) ∂ ˜ µ∂n (cid:19) (cid:104) η i ∗ ∇ j ( nξ j ) − ξ i ∇ j ( nη j ∗ ) (cid:105) ˆ n i (22)to be continuous. To see that this should be the case, use (again for a barotrope) dp = ndµ to get1 m B (cid:18) ∂p∂n (cid:19) (cid:104) η i ∗ ∇ j ( nξ j ) − ξ i ∇ j ( nη j ∗ ) (cid:105) ˆ n i = − m B (cid:18) ∂p∂n (cid:19) (cid:104) η i ∗ δ ξ n − ξ i δ η n ∗ (cid:105) ˆ n i = − m B (cid:104) η i ∗ δ ξ p − ξ i δ η p ∗ (cid:105) ˆ n i , (23)where δ ξ and δ η distinguishes the Eulerian perturbations associated with ξ i and η i , respectively. Continuity across a surfacerequires (i) the radial component of the displacement, and (ii) the Lagrangian pressure variation, to be continuous. The lastterm of equation (23) can be rewritten as (cid:104) η i ∗ δ ξ p − ξ i δ η p ∗ (cid:105) ˆ n i = (cid:104) η i ∗ ∆ ξ p − ξ i ∆ η p ∗ − (cid:16) η i ∗ ξ j ∇ j p − ξ i η ∗ j ∇ j p (cid:17)(cid:105) ˆ n i . (24)For a spherical star, this expression is clearly continuous as the background pressure depends only on the radial coordinate.Turning to the gravitational potential, we have (cid:90) ρη i ∗ ∇ i δ ξ Φ dV = (cid:90) (cid:104) ρη i ∗ δ ξ Φ (cid:105) ˆ n i dS − (cid:90) δ ξ Φ ∇ i ( ρη i ∗ ) dV = (cid:90) (cid:20) ρ (cid:16) η i ∗ δ ξ Φ − ξ i δ η Φ ∗ (cid:17) + 14 πG g ij ( δ ξ Φ ∇ j δ η Φ ∗ − δ η Φ ∗ ∇ j δ ξ Φ) (cid:21) ˆ n i dS + (cid:90) ρξ i ∇ i δ η Φ ∗ dV = (cid:90) ρξ i ∇ i δ η Φ ∗ dV , (25)where the last identity (again) holds as long as we may ignore the surface terms.In particular, at the surface of the star we need to ensure that ρ (cid:16) η i ∗ δ ξ Φ − ξ i δ η Φ ∗ (cid:17) + 14 πG g ij ( δ ξ Φ ∇ j δ η Φ ∗ − δ η Φ ∗ ∇ j δ ξ Φ) = 0 . (26)The first term (obviously) vanishes as long as ρ → δ ξ Φ ∂ r δ η Φ ∗ − δ η Φ ∗ ∂ r δ ξ Φ = − l + 1 R ( δ η Φ ∗ δ ξ Φ − δ ξ Φ δ η Φ ∗ ) r = R = 0 . (27)The conditions that apply when the density does not vanish at the surface, can be inferred from the results for an internaldensity jump. Again, assume that there is a density jump at an internal point, r = ¯ R , such that (in a small neighbourhoodof the phase transition) ρ ( r ) = (cid:40) ρ − , r < ¯ Rρ + , r > ¯ R (28)and integrate (7) over a small volume (with height = 2 (cid:15) ) across the (spherical) surface to get[ ∂ r δ Φ + 4 πGρξ r ] ¯ R + (cid:15) ¯ R − (cid:15) = 2 (cid:15)S (29) © , 000–000 mpact of the crust for some constant S . Then let (cid:15) → ∂ r δ Φ + + 4 πGρ + ξ r = ∂ r δ Φ − + 4 πGρ − ξ r , (30)where we have used the fact that the radial component of the displacement must be continuous. What does this mean for thesurface terms in (25)? Well, we integrate up to ¯ R and then continue on the other side. The surface terms associated with theinterface will vanish as long asˆ n i (cid:20) ρ (cid:16) η i ∗ δ ξ Φ − ξ i δ η Φ ∗ (cid:17) + 14 πG g ij ( δ ξ Φ ∇ j δ η Φ ∗ − δ η Φ ∗ ∇ j δ ξ Φ) (cid:21) ¯ R + (cid:15) ¯ R − (cid:15) → (cid:15) → . (31)We can rewrite this as [ δ ξ Φ ( ∂ r δ η Φ + 4 πGρη r ) ∗ − δ η Φ ∗ ( ∂ r δ ξ δ Φ + 4 πGρξ r )] ¯ R + (cid:15) ¯ R − (cid:15) → (cid:15) → . (32)Combining the fact that the perturbed gravitational potential is continuous with (30), we see that these contributions mustvanish. It is also easy to adjust the argument to arrive at the standard result for the gravitational potential matching at afinite-density surface. In effect, we have provided a (somewhat lengthy) demonstration that we do not have to worry aboutdensity discontinuities in the following. As long as the relevant junction conditions are respected, the usual orthogonalityargument remains unchanged.At the end of the day, we have the expected result (Friedman and Schutz 1978): (cid:104) η, Cξ (cid:105) = (cid:104) ξ, Cη (cid:105) ∗ . (33)This means that, if we consider two mode solutions (now letting ξ n be associated with frequency ω n and ξ m be the solutionassociated with ω m ) we have0 = (cid:104) ξ m , − ω n Aξ n + Cξ n (cid:105) = − ω n (cid:104) ξ n , Aξ m (cid:105) ∗ + (cid:104) ξ n , Cξ m (cid:105) ∗ = − ω n (cid:104) ξ n , Aξ m (cid:105) ∗ + ( ω ∗ m ) (cid:104) ξ n , Aξ m (cid:105) ∗ = (cid:2) ( ω ∗ m ) − ω n (cid:3) (cid:104) ξ n , Aξ m (cid:105) ∗ . (34)We see that, if the frequencies are distinct and real, the mode solutions must be orthogonal. That is, we may use (cid:104) ξ n , Aξ m (cid:105) = A n δ mn , (35)leaving the normalisation A n unspecified for the moment.So far, the arguments are standard. The only aspect we have added relates to internal phase transitions. Before we moveon let us, without detailed calculation, comment on the case of non-barotropic perturbations, relevant whenever we want toaccount for internal matter stratification. It is easy to argue that this case applies to neutron star tides as the relevant nuclearreactions are too slow to keep the perturbed matter in beta equilibrium (Andersson and Pnigouras 2020a). In particular, itleads to the presence of gravity g-modes in the oscillation spectrum.In the limit of slow reactions we may describe the equation of state in terms of two parameters. The most natural choiceis to complement the density ρ with either parameter representing the deviation from chemical equilibrium or the protonfraction x p = ρ p /ρ , the Lagrangian perturbation of which will vanish for slow reactions. However, the mode orthogonalitycan be established without assuming an explicit equation of state (keeping both the pressure and the density perturbations).The argument follows from Friedman and Schutz (1978), and we arrive at (cid:90) η i ∗ (cid:20) ∇ i δ ξ p − ρ δ ξ ρ ∇ i p (cid:21) dV = (cid:90) (cid:104) η i ∗ δ ξ p − ξ i δ η p ∗ (cid:105) ˆ n i dS + (cid:90) ξ i (cid:20) ∇ i δ η p ∗ − ρ δ η ρ ∗ ∇ i p (cid:21) dV . (36)This is the result we need to prove the mode orthogonality for stratified stars. As the argument for the gravitational po-tential remains unchanged, the stated symmetry relation (33) will hold as long as the radial displacement and the pressureperturbation are both continuous (for a spherical surface). We now take the “natural” next step towards a more realistic neutron star model, adding the elastic crust that is expectedto reach from (close to) the surface up to about 60-70% of the nuclear saturation density (roughly, the outer kilometer ofthe star). Assuming that the crust of the equilibrium star is “relaxed”—not associated with strains, the elasticity impactsonly on the perturbations. The argument is similar to that used in previous work on the tidal deformability and neutron star“mountains”, see for example Gittins, Andersson and Jones (2020). The Euler equation then changes to − ω n Aξ i + Cξ i + Eξ i = 0 , (37)where Eξ i = ∇ j σ ij (38)with the elastic stress tensor given by σ ij = ˇ µ ( ∇ i ξ j + ∇ j ξ i ) −
23 ˇ µg ij (cid:16) ∇ k ξ k (cid:17) . (39) ©000
23 ˇ µg ij (cid:16) ∇ k ξ k (cid:17) . (39) ©000 , 000–000 Passamonti, Andersson and Pnigouras
Note that the shear modulus ˇ µ must not to be confused with the chemical potential.Now, we need (cid:104) η i ∗ , Eξ i (cid:105) = (cid:90) η i ∗ ∇ j σ ij dV = (cid:90) (cid:104) η i ∗ σ ij (cid:105) ˆ n j dS − (cid:90) σ ij ∇ j η i ∗ dV . (40)The second term has two pieces: (cid:90) ˇ µ ( ∇ i ξ j + ∇ j ξ i ) ∇ j η i ∗ dV = (cid:90) ˇ µ (cid:0) ∇ i η ∗ j + ∇ j η ∗ i (cid:1) ∇ j ξ i dV (41)and (cid:90) ˇ µg ij (cid:16) ∇ k ξ k (cid:17) ∇ j η i ∗ dV = (cid:90) ˇ µ (cid:16) ∇ k ξ k (cid:17) (cid:16) ∇ j η ∗ j (cid:17) dV = (cid:90) ˇ µg ij (cid:16) ∇ k η ∗ k (cid:17) (cid:16) ∇ j ξ i (cid:17) dV . (42)Introducing ¯ σ ij = ˇ µ (cid:0) ∇ i η ∗ j + ∇ j η ∗ i (cid:1) −
23 ˇ µg ij (cid:16) ∇ k η ∗ k (cid:17) , (43)we see that (cid:104) η i ∗ , Eξ i (cid:105) = (cid:90) (cid:104) η i ∗ σ ij (cid:105) ˆ n j dS − (cid:90) ¯ σ ij ∇ j ξ i dV = (cid:90) (cid:104) η i ∗ σ ij − ξ i ¯ σ ij (cid:105) ˆ n j dS + (cid:90) ξ i ∇ j ¯ σ ij dV = (cid:90) (cid:104) η i ∗ σ ij − ξ i ¯ σ ij (cid:105) ˆ n j dS + (cid:104) ξ i , Eη ∗ i (cid:105) . (44)In order to deal with the surface terms we need to consider the traction conditions (Andersson, Haskell and Samuelsson2011). Hence, we impose the continuity of the gravitational potential and the perpendicular and radial components of thetraction vector t i = − (cid:16) g ij ∆ p − σ ij (cid:17) ˆ n j (45)noting that the traction reduces to the pressure perturbation in the fluid regions. In order to establish that the surface termsvanish at the crust-core transition, we combine the above result with the relevant term from (21). This leads to the requirementthat − (cid:104) η i ∗ δ ξ p − ξ i δ η p ∗ (cid:105) ˆ n i + (cid:104) η i ∗ σ ij − ξ i ¯ σ ij (cid:105) ˆ n j = − (cid:104) η i ∗ ( g ij ∆ ξ p − σ ij ) − ξ i ( g ij ∆ η p ∗ − ¯ σ ij ) − g ij ( η i ∗ ξ k − ξ i η ∗ k ) ∇ k p (cid:105) ˆ n j (46)is continuous at the interface. We see that, given the traction conditions, this should always be the case for a spherical star.In essence, the addition of an elastic region does not complicate the formal analysis. Finally, we are set to return to the tidal problem. This part of the analysis is straightforward given the properties we haveestablished. Perhaps not surprisingly, the results we need are identical to those from the fluid case, see, for example, Anderssonand Pnigouras (2020a). In particular, it is clear that the oscillation modes of the fluid+elastic problem are orthogonal andthat (35) still holds. As in the fluid problem, we can use this fact to rewrite the Euler equation ρ∂ t ξ i + Cξ i + Eξ i = − ρ ∇ i χ (47)as an equation for individual mode amplitudes. Introducing the mode expansion ξ i = (cid:88) n a n ( t ) ξ in , (48)where the eigenfunctions ξ in are time-independent (as they are obtained in the frequency domain), we easily arrive at¨ a n + ω n a n = − A n (cid:104) ξ n , ρ ∇ χ (cid:105) = − A n (cid:90) χδρ ∗ n dV (49)making use of the continuity equation and integrating by parts in the last step.Expanding the tidal potential in spherical harmonics (Andersson and Pnigouras 2020a) we have χ = (cid:88) l,m v l r l Y lm , (50)where the v l coefficients will be given later. Assuming an adiabatic inspiral and working in the frequency domain (with timedependence e iωt ), the tidal driving only has support at frequency ω = m Ω. It is useful to keep this in mind. Anyway, we arriveat a set of driven modes with amplitude a n = 1 ω n − ω Q n A n v l , (51) © , 000–000 mpact of the crust where we have introduced the overlap integral Q n = − (cid:90) δρ ∗ n r l +2 dr . (52)As discussed by Andersson and Pnigouras (2020a) the overlap integral can also be expressed in terms of the matching ofthe gravitational potential at the star’s surface4 πG (cid:90) R r l +2 δρdr = − (2 l + 1) R l +1 δ Φ( R ) , (53)so we have Q n = 2 l + 14 πG R l +1 δ Φ n ( R ) = I n , (54)where I n represents the contribution each of the star’s oscillation modes makes to the mass multipole moment.Finally, we connect the mode expansion to the tidal deformability and the effective Love number by introducing arepresentation for the perturbed gravitational potential in terms of the mode eigenfunctions. Expressing the displacementvector for each multipole as ξ i = (cid:18) W ( r ) ∇ i rr (cid:19) Y lm + V ( r ) ∇ i Y lm , (55)it follows from the θ -component of (4) that (since ∆ p = 0 at the surface) − ω V ( R ) − p (cid:48) ρ W ( R ) R + δ Φ( R ) = − χ ( R ) . (56)That is, we have δ Φ( R ) = − χ ( R ) + ω V ( R ) − g W ( R ) R (57)and we arrive at an expression for the dynamical tidal response k eff l = 12 δ Φ( R ) χ ( R ) = −
12 + 12 v l R l (cid:20) ω V ( R ) − g W ( R ) R (cid:21) . (58)Moreover, making contact with (51) we see that ξ i = (cid:88) n ω n − ω Q n v l A n ξ in , (59)which means that we have (simply adding a label, n , to indicate the individual mode eigenfunctions) k eff l = −
12 + 12 R l (cid:88) n Q n A n ω n − ω (cid:20) ω V n ( R ) − GM (cid:63) R W n ( R ) (cid:21) (60)with M (cid:63) being the primary’s mass. Keeping the normalisation of the modes unspecified—which may be useful when weconsider the numerical aspects of the problem—we have A n = (cid:90) R ρ (cid:2) W n + l ( l + 1) V n (cid:3) dr . (61)Introducing the dimensionless frequency ˜ ω = ω R GM (cid:63) (62)as well as the dimensionless overlap integral ˜ Q n = Q n M (cid:63) R l , (63)we get k eff l = − − (cid:88) n ˜ Q n ˜ ω n − ˜ ω M (cid:63) W n ( R ) A n (cid:20) − ˜ ω (cid:18) V n W n (cid:19) R (cid:21) . (64)This is the final result, but we can massage it a bit by recalling (54) and using δ Φ n ( R ) = ω n V n ( R ) − g W n ( R ) R . (65)We have (Andersson and Pnigouras 2020a) W n ( R ) = − π l + 1 ˜ Q n R (cid:20) − ˜ ω n (cid:18) V n W n (cid:19) R (cid:21) − (66)which, when combined with (64), leads to the final expression k eff l = −
12 + 2 π l + 1 (cid:88) n ˜ Q n ˜ ω n − ˜ ω (cid:18) M (cid:63) R A n (cid:19) (cid:20) − ˜ ω (cid:18) V n W n (cid:19) R (cid:21) (cid:20) − ˜ ω n (cid:18) V n W n (cid:19) R (cid:21) − . (67) ©000
12 + 2 π l + 1 (cid:88) n ˜ Q n ˜ ω n − ˜ ω (cid:18) M (cid:63) R A n (cid:19) (cid:20) − ˜ ω (cid:18) V n W n (cid:19) R (cid:21) (cid:20) − ˜ ω n (cid:18) V n W n (cid:19) R (cid:21) − . (67) ©000 , 000–000 Passamonti, Andersson and Pnigouras -5 -4 -3 -2 -1 0 1 log ( 1 - r/R )0.00.20.40.60.81.0 f modeVW G = 7/3crust coreocean Figure 1.
Lagrangian displacements for the l = 2 f mode. The black and red curves show, respectively, the functions W = W/ ( rR ) and V = V/ ( rR ) for a model with Γ = 7 / In the low-frequency limit, we obtain a mode sum for the tidal Love number k l ≈ −
12 + 2 π l + 1 (cid:88) n ˜ Q n ˜ ω n (cid:18) M (cid:63) R A n (cid:19) (cid:20) − ˜ ω n (cid:18) V n W n (cid:19) R (cid:21) − = −
12 + (cid:88) n k nl . (68)The careful reader will note that the expressions for the effective Love number remain unchanged from Andersson andPnigouras (2020a). This is to be expected, as long as we consider a neutron star with a fluid (ocean) surface. The resultswould change if we were to impose a direct transition from the elastic region to the exterior vacuum as we would then haveto consider elastic terms in, for example, (56). A fluid ocean is expected on physical grounds (even though it may be veryshallow for mature neutron stars). Moreover, as the ocean may sustain its own (more or less distinct) set of oscillation modes(g-modes) it is relevant to include this aspect in the model. Having discussed the theoretical framework, let us move on to consider the numerical results for a sample of neutron starmodels. In order to facilitate a direct comparison, we have chosen to focus on the stratified polytropic models already consideredby Andersson and Pnigouras (2020a). The stellar models we consider then involve three distinct regions: the single fluid core,the elastic crust and a (shallow) fluid ocean. In effect, there will be new features, like shear modes associated with the elasticityand interface modes linked to the core-crust and crust-ocean transitions (McDermott et al 1985; McDermott, van Horn andHansen 1988). We do not consider the impact of superfluidity in the core or the crust at this point, but plan to return to thisproblem later.In the first instance, we are interested in the static and dynamical aspects of the tide, as represented (via the differentoscillation modes) by the tidal deformability and the Love number k eff l . As we anticipate the impact of the elasticity on alreadyexisting fluid modes to be slight and the tidal excitation of, for example, the shear oscillations in the crust to be weak, wehave to make sure our numerical approach is robust. We need to do better than the proof-of-principle analysis of Anderssonand Pnigouras (2020a). Indeed, in order to reach the desired precision in the relevant eigenfunctions and the overlap integrals,we have to improve on the approach from Passamonti and Andersson (2012). The results we present in the following wereobtained using a “classic” shoot-and-relax approach, where a shooting method was used to identify each oscillation modeand the precision was subsequently improved by a relaxation step. That this implementation allows us to reliably extract thequantities we are interested in should be evident from the results we provide.Specifically, we consider a simple polytropic equation of state, p = kρ Γ , with Γ = 2. To account for the effects ofstratification on the oscillation modes and the tide, we introduce an adiabatic index, Γ , for the perturbations:∆ ρρ = 1Γ ∆ pp . (69)Following Andersson and Pnigouras (2020a) we consider three models, a barotropic star with Γ = 2 and weakly and stronglystratified stars with, respectively, Γ = 2 .
05 and Γ = 7 /
3. For simplicity, we assume that this equation of state describes © , 000–000 mpact of the crust -5 -4 -3 -2 -1 0 log ( 1 - r/R )-100-50050100150200250300 W G = 2.05 crust coreocean g modeg mode - fluid model -5 -4 -3 -2 -1 0 log ( 1 - r/R )-100-50050100150200250300 V G = 2.05 crust coreocean g modeg mode - fluid model Figure 2.
Lagrangian displacements for the g mode. The two panels show, respectively, W = W/ ( rR ) (left panel) and V = V/ ( rR ) (rightpanel). The solid lines indicate the solutions for the Γ = 2 .
05 star with a crust while the dashed blue line represents the correspondingfluid case. the entire star, from the core to the ocean. This may not be particularly realistic, but as we are working in the frameworkof Newtonian gravity we have to settle for a somewhat phenomenological set-up. From equation (69) it follows that we canwrite the relation between the Eulerian perturbations of pressure and mass density δρρ = 1Γ δpp − A i ξ i , (70)where A i = ∇ i ln ρ − ∇ i ln p Γ (71)is the Schwarzschild discriminant, which quantifies the presence of buoyancy and determines the properties of gravity modes.All numerical results presented in the following were obtained with the crust-core transition taken to be at R cc = 0 . R and the crust-ocean interface at R co = 0 . R . The elastic properties of the crust are described by the shear modulus ˇ µ .Following Douchin and Haensel (2001) we assume that the specific shear modulus is nearly constant across the crust. Thiswould mean taking ˇ µρ (cid:39) cm / s . (72)but in practice we use the (further) simplified version from Passamonti and Andersson (2012) and letˇ µ = ˜ µ ρ (cid:18) GM (cid:63) R (cid:19) , (73)where ˜ µ is a dimensionless parameter, which must not be confused with the chemical potential defined in equation (12). Fora star with M (cid:63) = 1 . M (cid:12) and R = 10 km we then find GM (cid:63) /R = 1 . × cm / s . Linking to (72) we assume the value˜ µ = 10 − for most of our examples, exploring other cases only for the interface modes. Even though it is fairly simple, the stellar model we consider sustains many (more or less) distinguishable classes of oscillationmodes: the fundamental mode (f mode), the pressure modes (p modes), the gravity modes (g modes), the shear modes (smodes) and the interface modes (i modes) (see McDermott et al (1985) and McDermott, van Horn and Hansen (1988) fordetailed discussions). The f and p modes are present also in a fluid star, while the g modes appear in stratified models (in ourcase when Γ (cid:54) = Γ), see Andersson and Pnigouras (2020a). The presence of the crust impacts on these classes of modes andintroduces, due to the elasticity, the shear modes. Moreover, any transition between different regions in the star (core-crustand crust-ocean) tends to be associated with a set of interface modes.Since we are interested in the tidal excitation of the different modes by a binary companion and the possible impact onthe gravitational-wave signal, we focus our attention on the quadrupole ( l = 2) modes. We then have the sets of oscillationfrequencies provided (in dimensionless units) in Tables 1, 2 and 3. The tabulated data include other relevant quantities requiredfor the tidal problem (e.g. the overlap integrals, the ratio at the surface between the tangential and radial components of ©000
05 star with a crust while the dashed blue line represents the correspondingfluid case. the entire star, from the core to the ocean. This may not be particularly realistic, but as we are working in the frameworkof Newtonian gravity we have to settle for a somewhat phenomenological set-up. From equation (69) it follows that we canwrite the relation between the Eulerian perturbations of pressure and mass density δρρ = 1Γ δpp − A i ξ i , (70)where A i = ∇ i ln ρ − ∇ i ln p Γ (71)is the Schwarzschild discriminant, which quantifies the presence of buoyancy and determines the properties of gravity modes.All numerical results presented in the following were obtained with the crust-core transition taken to be at R cc = 0 . R and the crust-ocean interface at R co = 0 . R . The elastic properties of the crust are described by the shear modulus ˇ µ .Following Douchin and Haensel (2001) we assume that the specific shear modulus is nearly constant across the crust. Thiswould mean taking ˇ µρ (cid:39) cm / s . (72)but in practice we use the (further) simplified version from Passamonti and Andersson (2012) and letˇ µ = ˜ µ ρ (cid:18) GM (cid:63) R (cid:19) , (73)where ˜ µ is a dimensionless parameter, which must not be confused with the chemical potential defined in equation (12). Fora star with M (cid:63) = 1 . M (cid:12) and R = 10 km we then find GM (cid:63) /R = 1 . × cm / s . Linking to (72) we assume the value˜ µ = 10 − for most of our examples, exploring other cases only for the interface modes. Even though it is fairly simple, the stellar model we consider sustains many (more or less) distinguishable classes of oscillationmodes: the fundamental mode (f mode), the pressure modes (p modes), the gravity modes (g modes), the shear modes (smodes) and the interface modes (i modes) (see McDermott et al (1985) and McDermott, van Horn and Hansen (1988) fordetailed discussions). The f and p modes are present also in a fluid star, while the g modes appear in stratified models (in ourcase when Γ (cid:54) = Γ), see Andersson and Pnigouras (2020a). The presence of the crust impacts on these classes of modes andintroduces, due to the elasticity, the shear modes. Moreover, any transition between different regions in the star (core-crustand crust-ocean) tends to be associated with a set of interface modes.Since we are interested in the tidal excitation of the different modes by a binary companion and the possible impact onthe gravitational-wave signal, we focus our attention on the quadrupole ( l = 2) modes. We then have the sets of oscillationfrequencies provided (in dimensionless units) in Tables 1, 2 and 3. The tabulated data include other relevant quantities requiredfor the tidal problem (e.g. the overlap integrals, the ratio at the surface between the tangential and radial components of ©000 , 000–000 Passamonti, Andersson and Pnigouras -5 -4 -3 -2 -1 0 log ( 1 - r/R )-20-16-12-8-4048 i modeV x 10 -3 VW G = 2crust coreocean -5 -4 -3 -2 -1 0 log ( 1 - r/R )-10-8-6-4-202468 i modeV x 10 -3 VW G = 2.05crust coreocean -5 -4 -3 -2 -1 0 log ( 1 - r/R )-20-15-10-505 i modeV x 10 -3 VW G = 7/3crust coreocean Figure 3.
Eigenfunctions for the i mode, which is associated with the crust-ocean transition (note the distinctive cusp in the radialdisplacement). The black curves represent the radial component W = W/ ( rR ) while the red curves are for the tangential component V = V/ ( rR ). The left-hand panel shows the results for a barotropic model (Γ = 2). The middle panel represents the Γ = 2 .
05 case,while the right-hand panel is for a stratified model with Γ = 7 / -6 -5 -4 -3 -2 -1 0 log ( 1 - r/R )-12-8-404812 g s1 modeV x 10 -3 VW G = 2.05crust coreocean -6 -5 -4 -3 -2 -1 0 log ( 1 - r/R )-50510152025 g s2 modeV x 10 -3 VW G = 2.05crust coreocean Figure 4.
Eigenfunctions for the first two surface g modes for the stratified model with Γ = 7 /
3. The left-hand panel shows the W and V eigenfunctions for the g s1 mode, while right-hand panel displays the same quantities for the g s2 mode. the Lagrangian displacement and the mode contribution to the Love number, all discussed in Section 2.4). The p n modes(with n labelling each mode) cover the high frequency range of the spectrum—above the frequency of the f mode—andtheir frequencies increase for higher order ( n ) modes. In contrast, the (gravity) g modes cover the low frequency range ofthe oscillation spectrum—below the f mode—and their frequencies decrease for higher order modes. As a result, the g-modespectrum becomes progressively dense towards zero frequency. This makes the identification of the individual high-order gmodes challenging. A stratified ocean also sustains a family of g modes, here referred to as surface g modes (g s ). Thesealso have low frequencies and their eigenfunctions are mainly confined to the star’s ocean. The shear mode frequency scalesapproximately as ω n ∼ (cid:112) ˇ µ/ρ and increases for higher order modes. Finally, the interface modes are characterised by theirlow frequency and a distinctive cusp in the radial displacement at the relevant interface.The main contribution to the dynamical tide is provided by the f mode (Hinderer et al 2016; Steinhoff et al 2016;Andersson and Pnigouras 2020a,b). This mode is known to be weakly affected by composition stratification (see the results ofAndersson and Pnigouras 2020a) and crust elasticity (McDermott, van Horn and Hansen 1988). This is demonstrated by theresults in Figure 1, where we show the l = 2 f-mode eigenfunction for a strongly stratified model with Γ = 7 /
3. The radialcomponent ( W ) of the Lagrangian displacement is practically equal to the barotropic case, while the tangential component( V ) exhibits oscillations in the crust region (compared to the result for a fluid star). Note that, in stars with a crust thetangential displacement may be discontinuous at the fluid-elastic transitions. This is also evident from Figure 1.Qualitatively similar conclusions apply to the g and p modes. As an example, we consider the g mode for the Γ = 2 . © , 000–000 mpact of the crust g g g g g f p p p p -5 -4 -3 -2 -1 | Q n | ~ G = 2 G = 2.05 G = 7/3 g g g g g f p p p p -5 -4 -3 -2 -1 | k n l | G = 2 G = 2.05 G = 7/3 Figure 5.
Overlap integrals (left panel) and mode contributions to the Love number k nl (right panel) for the fundamental, pressure andgravity modes for the three stellar models we consider. Different markers indicate the different values of Γ (see the legend). we show the Lagrangian displacements for this mode. As before, the crust mainly impacts on the tangential component ( V )which appears to be almost constant in the crust.Associated with each fluid-elastic interface we find a single interface mode. For these modes, the radial displacement ispeaked at the interface. This feature is evident from Figure 3, where we show the displacements for the crust-ocean interfacemode (i ) for the barotropic model as well as for the weakly and strongly stratified stars. As anticipated, the W eigenfunctionhas a notable cusp at the transition between the crust and the ocean. For the barotropic case, there is also another cusp(albeit with smaller amplitude) at the crust-core interface. As in the case of other modes, the tangential eigenfunction V isdiscontinuous at the interfaces, and reaches a very large amplitude at the star’s surface. Note that the function V has beenreduced, in Figure 3, by a factor 10 − in the ocean (only). Similar large amplitudes in the ocean were noted by McDermott,van Horn and Hansen (1988).The stratification associated with the matter composition affects the i-mode eigenfunctions, especially in the core, seeFigure 3. In this region, the character of the interface modes is very similar to that of the (core) g modes. This behaviourusually occurs when two modes lie in the same frequency range and the resulting oscillation mode exhibits a mixed character.In this case, we find g-mode features in the core and the characteristic interface mode cusp at the transition density.The other interface mode (i ), associated with the crust-core transition, is readily determined for the barotropic case(Γ = 2), for which it has the expected properties. However, we find it difficult to identify this mode for the stratified models.The problem is most likely due to the fact that, in the core, the i and (core) g modes have similar eigenfunctions and it isdifficult to numerically separate them. This behaviour is shown in Figure 2 for the g mode, where it is clear that the W eigenfunction is not differentiable at the crust-core transition. For higher-order g modes this behaviour is even more apparent.The mixing of modes belonging to different classes is also notable for the surface g modes. This is illustrated in Figure 4,where we show, for the stratified model with Γ = 7 /
3, the W and V eigenfunctions for the first two g s modes. We see that,in the ocean we have the characteristic eigenfunction of a surface g mode, while the eigenfunctions W and V show featuresassociated with a high-order g mode in the core. The crust region appears to “isolate” the two regions which makes it difficultto numerically establish which set of modes a given solution belongs to. It is not even clear that this is a meaningful questionfor the high overtone modes. The tidal response of a neutron star is closely related to the nature of the different oscillation modes. This is natural sincethe modes form a complete set and, hence, can be used as a basis to express the behaviour of the stellar fluid. In the staticlimit the sum over the star’s oscillation modes leads to the Love number, as explicitly demonstrated by Andersson andPnigouras (2020a). The mode-sum also provides a handle on the dynamical tide (through the effective, frequency dependent,Love number k eff l from Section 2.4). In particular, during a binary inspiral some of the oscillation modes may pass throughresonance with the tidal driving and as a result reach a significant amplitude.Let us first consider the static limit. In this case oscillation modes are (clearly) not in resonance with the orbital motion.Nevertheless, we can use the mode-sum to represent the tidal response. This may seem a somewhat odd way to go about ©000
3, the W and V eigenfunctions for the first two g s modes. We see that,in the ocean we have the characteristic eigenfunction of a surface g mode, while the eigenfunctions W and V show featuresassociated with a high-order g mode in the core. The crust region appears to “isolate” the two regions which makes it difficultto numerically establish which set of modes a given solution belongs to. It is not even clear that this is a meaningful questionfor the high overtone modes. The tidal response of a neutron star is closely related to the nature of the different oscillation modes. This is natural sincethe modes form a complete set and, hence, can be used as a basis to express the behaviour of the stellar fluid. In the staticlimit the sum over the star’s oscillation modes leads to the Love number, as explicitly demonstrated by Andersson andPnigouras (2020a). The mode-sum also provides a handle on the dynamical tide (through the effective, frequency dependent,Love number k eff l from Section 2.4). In particular, during a binary inspiral some of the oscillation modes may pass throughresonance with the tidal driving and as a result reach a significant amplitude.Let us first consider the static limit. In this case oscillation modes are (clearly) not in resonance with the orbital motion.Nevertheless, we can use the mode-sum to represent the tidal response. This may seem a somewhat odd way to go about ©000 , 000–000 Passamonti, Andersson and Pnigouras g g g g g -5 -4 -3 -2 -1 | Q n | ~ G = 2.05 G = 7/3 G = 2.05 - fluid model G = 7/3 - fluid model g g g g g -5 -4 -3 -2 -1 | k n l | G = 2.05 G = 7/3 G = 2.05 - fluid model G = 7/3 - fluid model Figure 6.
Comparing models with and without a crust. We show the (dimensionless) overlap integral Q n (left panel) and the modecontribution to the Love number k nl (right panel) for the first five gravity modes. The results for purely fluid stars are shown as emptymarkers, while the results for models with crust and ocean are shown as filled markers (see the legend). As a general trend, each quantitydecreases for the higher order modes, but there are clearly exceptions to this. i i s s s s s s -6 -5 -4 -3 -2 | Q n | ~ G = 2 G = 2.05 G = 7/3 i i s s s s s s -6 -5 -4 -3 -2 | k n l | G = 2 G = 2.05 G = 7/3 Figure 7.
Overlap integrals (left panel) and mode contributions to the Love number k nl (right panel) for the shear and interface modesfor the three stellar models we consider. Different markers indicate the different values of Γ (see the legend). it, given that the result we want is contained in the usual Love number k l which can be calculated in a much simpler way(Hinderer et al 2010). However, the mode representation provides valuable additional insight. In particular, it brings out theexpectation that the main contribution to the tidal response is provided by the f mode, which has the best overlap withthe tidal driving force. This also leads to the question of the level at which other modes, which may depend on the mattercomposition etcetera, contribute. For example, we know that the elastic crust sustains shear modes. These are expected tohave small overlap integrals and therefore have little impact. However, even if these expectations are true, it is useful toquantify what we mean by “small” and to what extent we can safely neglect the contribution of these modes. We need to bemindful of the fact that, even if the contribution from each mode is too small to impact on the gravitational-wave signal, thepresence of an excited mode may have other repercussions. For example, it is possible that a mode passing through resonancereaches an amplitude where it leads to fracturing of the crust (Tsang et al 2012; Pan et al 2020). We will return to thispossibility later. Finally, we already know from the eigenfunctions we have provided that the presence of the crust affects theproperties of all modes (see Section 3.1). This means that there will be a (probably weak, but nevertheless) impact on therespective overlap integral and the contribution to the tidal response.As explained in Section 2.4, in the low-frequency limit (low in the sense that ω (cid:28) ω n ), the tidal Love number may be © , 000–000 mpact of the crust -3 -2 -1 w n -6 -5 -4 -3 -2 -1 | k n l | f ~ p g s i i i i s g g g p g g g g G = 2 G = 2.05 G = 7/3 Figure 8.
Summary results for all the modes we have considered. We show the Love number contribution, k nl , against the modefrequencies (in dimensionless units). The illustrated modes are the fundamental (f) mode, the pressure (p n ) and gravity (g n ) modes, theshear (s n ) modes and finally the interface (i) modes which arise from the core-crust and crust-ocean interfaces. written as a sum over individual mode contributions k nl (see Equation (68)). Hence, we report in Tables 1, 2 and 3 the valueof k nl for the different oscillation modes (along with other relevant quantities). The results confirm that the f mode dominatesthe tidal response. It has the largest value for the overlap integral and makes the dominant contribution to the Love number.The contributions from the pressure and gravity modes become less important for the higher overtones (increasing n , seeFigures 5 and 8), in accordance with the results of Andersson and Pnigouras (2020a). For strongly stratified models, withΓ = 7 /
3, the first g mode has a value for k nl similar to that of the first p mode.Comparing the results for stellar models with and without a crust, we find that the pressure modes are essentially notaffected at all, while the gravity modes have a less regular behaviour. As shown in Figure 6, the contribution of the g modesto the Love number, k nl , tends to be smaller for the model with a crust (although it is practically unchanged for the g mode).In strongly stratified models, the g modes are less influenced by the crust, and the corresponding values of Q n and k nl aresimilar to the fluid case, most likely because the buoyancy dominates the elastic restoring force. We also note an interestingresult for the g mode. In the Γ = 2 .
05 model, the overlap integral, Q n , is about an order of magnitude larger than in thefluid case, while the value for k nl does not change correspondingly (see Figures 5 and 6). This demonstrates that the overlapintegral is not the only factor that determines the contribution to the tidal response.The shear modes in the crust are, as expected, more or less irrelevant for the tidal problem, as confirmed by the smallvalues for Q n and k nl . In Figure 7 we show the two quantities for the first six s modes. As in the case of the g modes, thebehaviour of the overlap integral is less regular with respect to the radial number n than what we find for k nl . Meanwhile, theinterface modes linked to the core-crust and crust-ocean transitions may contribute to the tidal response at a level similar tothat of the first or second g modes (see Figures 7 and 8). The importance of these modes increases for models with strongerstratification, although their overlap integrals remain very small (see the left panel of Figure 7).An overall summary view of the results for the Love number is provided in Figure 8, where we show the quantity k nl forall modes considered in this work. From this figure we infer which are the most relevant modes in the static limit and in whichorder the modes will be excited during a binary inspiral (as the driving frequency increases). In this figure we also show k nl for the first three surface g modes, which mainly reside in the ocean.Finally, from a technical perspective, it is worth noting that some of the oscillation modes have very small overlapintegrals, the calculation of which may be subject to numerical errors. This is particularly the case for high-order modeswhich tend to have many nodes in their eigenfunctions, leading to cancellations in the calculation of the overlap integral. Tomonitor the numerical errors we determine Q n from equations (52), (54) as well as the (equivalent) expression (Lai 1994) Q n = l (cid:90) ρ [ W + ( l + 1) V ] r l dr . (74)As we have already mentioned, in order to increase the precision of the calculation, we solve the perturbation equations firstwith a multiple shooting method and then with a relaxation step. As expected, the solutions obtained from the relaxation aremore accurate, which allows us to extract the high-order modes. A few additional comments on the technical aspects of thecalculation are provided in Appendix A. ©000
05 model, the overlap integral, Q n , is about an order of magnitude larger than in thefluid case, while the value for k nl does not change correspondingly (see Figures 5 and 6). This demonstrates that the overlapintegral is not the only factor that determines the contribution to the tidal response.The shear modes in the crust are, as expected, more or less irrelevant for the tidal problem, as confirmed by the smallvalues for Q n and k nl . In Figure 7 we show the two quantities for the first six s modes. As in the case of the g modes, thebehaviour of the overlap integral is less regular with respect to the radial number n than what we find for k nl . Meanwhile, theinterface modes linked to the core-crust and crust-ocean transitions may contribute to the tidal response at a level similar tothat of the first or second g modes (see Figures 7 and 8). The importance of these modes increases for models with strongerstratification, although their overlap integrals remain very small (see the left panel of Figure 7).An overall summary view of the results for the Love number is provided in Figure 8, where we show the quantity k nl forall modes considered in this work. From this figure we infer which are the most relevant modes in the static limit and in whichorder the modes will be excited during a binary inspiral (as the driving frequency increases). In this figure we also show k nl for the first three surface g modes, which mainly reside in the ocean.Finally, from a technical perspective, it is worth noting that some of the oscillation modes have very small overlapintegrals, the calculation of which may be subject to numerical errors. This is particularly the case for high-order modeswhich tend to have many nodes in their eigenfunctions, leading to cancellations in the calculation of the overlap integral. Tomonitor the numerical errors we determine Q n from equations (52), (54) as well as the (equivalent) expression (Lai 1994) Q n = l (cid:90) ρ [ W + ( l + 1) V ] r l dr . (74)As we have already mentioned, in order to increase the precision of the calculation, we solve the perturbation equations firstwith a multiple shooting method and then with a relaxation step. As expected, the solutions obtained from the relaxation aremore accurate, which allows us to extract the high-order modes. A few additional comments on the technical aspects of thecalculation are provided in Appendix A. ©000 , 000–000 Passamonti, Andersson and Pnigouras
Table 1.
Mode results for the barotropic model with Γ = 2. The first four columns, respectively, report the oscillation mode, thedimensionless mode frequency, the absolute value of the dimensionless overlap integral and the ratio between the radial and tangentialcomponents of the displacement vector at the surface. The last column provides the contribution each mode makes to the static Lovenumber, k nl . We also provide, in the final row, the summed result which should be compared to the anticipated value k l = 0 . ω n | ˜ Q n | ( V/W ) R k nl p . × − . × − . × − p . × − . × − − . × − p . × − . × − . × − p . × − . × − − . × − f 1.2269 5 . × − . × − . × − s . × − . . × − s . × − . − . × − s . × − . . × − s . × − . − . × − s . × − . − . × − s . × − . . × − i . × − . × . × − i . × − . × − . × − k l Table 2.
Same as Table 1, for the weakly stratified model with Γ = 2 . ω n | ˜ Q n | ( V/W ) R k nl p . × − . × − . × − p . × − . × − − . × − p . × − . × − . × − p . × − . × − − . × − f 1.2274 5 . × − . × − . × − g . × − . × . × − g . × − . × . × − g . × − . × . × − g . × − . × − . × − g . × − . × − . × − s . × − . − . × − s . × − . − . × − s . × − . . × − s . × − . − . × − s . × − . − . × − s . × − . . × − i . × − . × . × − g s1 . × − . × − . × − . × − g s2 . × − . × − . × . × − g s3 . × − . × − . × . × − k l As different oscillation modes pass through resonance during a binary inspiral, their amplitude may become large enoughthat the motion in the crust induces (local) fracturing of the nuclear lattice. It has been suggested that the interface modesare particularly relevant in this respect (Tsang et al 2012; Pan et al 2020). It also known that the static tide is unlikely tobreak the crust before the binary merger (Penner et al 2012; Gittins, Andersson and Pereira 2020). In order to consider thisproblem we need to complement our mode analysis with the energy deposited in each mode during inspiral. This analysisclosely follows, for example, Lai (1994) so we will only outline the steps here. Some further details are provided in AppendixB. © , 000–000 mpact of the crust Table 3.
Same as Table 1, for the strongly stratified model with Γ = 7 / ω n | ˜ Q n | ( V/W ) R k nl p . × − . × − . × − p . × − . × − − . × − p . × − . × − . × − p . × − . × − − . × − f 1.2294 5 . × − . × − . × − g . × − . . × − g . × − . × . × − g . × − . × . × − g . × − . × . × − g . × − . × . × − s . × − . − . × − s . × − . . × − s . × − . − . × − s . × − . . × − s . × − . . × − s . × − . − . × − i . × − . × . × − g s1 . × − . × − . × − g s2 . × − . × − . × . × − g s3 . × − . × − . × − . × − k l The binary separation, D , shrinks at a rate which can be described, to leading order, by˙ D = − G c M (cid:48) M (cid:63) ( M (cid:63) + M (cid:48) ) D , (75)where M (cid:48) is the mass of the companion. This leading-order expression should be sufficient as long as we are not trying toresolve the fine details of the problem. We already know that the effects of the tide on the orbital evolution enter at (much)higher post-Newtonian order (Lai 1994; Kokkotas and Schaefer 1995; Flanagan and Hinderer 2008). The orbital frequency Ωfollows from Kepler’s law, so we have Ω = (cid:20) G ( M (cid:63) + M (cid:48) ) D (cid:21) / . (76)We know from the perturbation analysis that the mode amplitude can be calculated from Equation (49), where we needthe tidal potential χ = − GM (cid:48) (cid:88) l ≥ m = l (cid:88) m = − l W lm r l D ( t ) l +1 Y lm e − im Φ( t ) (77)with Φ = (cid:82) Ω( t )d t and explicitly defining the v l coefficient we used earlier. For l = 2 the W lm coefficients have the followingvalues: W = − (cid:114) π , W ± = 0 , W ± = (cid:114) π . (78)For a given mode ( l, m ), Equation (49) leads to¨ a n + ω n a n = GM (cid:48) R ˜ Q n ˜ A n (cid:18) RD ( t ) (cid:19) l +1 W lm e − im Φ( t ) (79)where ˜ Q n is the dimensionless overlap integral defined in equation (63) and ˜ A n = A n / ( M (cid:63) R ).Equation (79) is a forced harmonic oscillator for the mode amplitude a n , where the forcing term is provided by the tidalpotential. An oscillation mode is in resonance with the orbit when ω n (cid:39) m Ω. For the most relevant modes, the resonanceoccurs at the late stages of binary inspiral, when the separation changes rapidly and hence the energy transfer to the modeis limited. Using the resonance condition in equation (76) we have the (dimensionless) resonance distance (Lai 1994) DR = (cid:20) m (1 + q )˜ ω n (cid:21) / (80)where q = M (cid:48) /M (cid:63) and ˜ ω n is the dimensionless mode frequency. ©000
Same as Table 1, for the strongly stratified model with Γ = 7 / ω n | ˜ Q n | ( V/W ) R k nl p . × − . × − . × − p . × − . × − − . × − p . × − . × − . × − p . × − . × − − . × − f 1.2294 5 . × − . × − . × − g . × − . . × − g . × − . × . × − g . × − . × . × − g . × − . × . × − g . × − . × . × − s . × − . − . × − s . × − . . × − s . × − . − . × − s . × − . . × − s . × − . . × − s . × − . − . × − i . × − . × . × − g s1 . × − . × − . × − g s2 . × − . × − . × . × − g s3 . × − . × − . × − . × − k l The binary separation, D , shrinks at a rate which can be described, to leading order, by˙ D = − G c M (cid:48) M (cid:63) ( M (cid:63) + M (cid:48) ) D , (75)where M (cid:48) is the mass of the companion. This leading-order expression should be sufficient as long as we are not trying toresolve the fine details of the problem. We already know that the effects of the tide on the orbital evolution enter at (much)higher post-Newtonian order (Lai 1994; Kokkotas and Schaefer 1995; Flanagan and Hinderer 2008). The orbital frequency Ωfollows from Kepler’s law, so we have Ω = (cid:20) G ( M (cid:63) + M (cid:48) ) D (cid:21) / . (76)We know from the perturbation analysis that the mode amplitude can be calculated from Equation (49), where we needthe tidal potential χ = − GM (cid:48) (cid:88) l ≥ m = l (cid:88) m = − l W lm r l D ( t ) l +1 Y lm e − im Φ( t ) (77)with Φ = (cid:82) Ω( t )d t and explicitly defining the v l coefficient we used earlier. For l = 2 the W lm coefficients have the followingvalues: W = − (cid:114) π , W ± = 0 , W ± = (cid:114) π . (78)For a given mode ( l, m ), Equation (49) leads to¨ a n + ω n a n = GM (cid:48) R ˜ Q n ˜ A n (cid:18) RD ( t ) (cid:19) l +1 W lm e − im Φ( t ) (79)where ˜ Q n is the dimensionless overlap integral defined in equation (63) and ˜ A n = A n / ( M (cid:63) R ).Equation (79) is a forced harmonic oscillator for the mode amplitude a n , where the forcing term is provided by the tidalpotential. An oscillation mode is in resonance with the orbit when ω n (cid:39) m Ω. For the most relevant modes, the resonanceoccurs at the late stages of binary inspiral, when the separation changes rapidly and hence the energy transfer to the modeis limited. Using the resonance condition in equation (76) we have the (dimensionless) resonance distance (Lai 1994) DR = (cid:20) m (1 + q )˜ ω n (cid:21) / (80)where q = M (cid:48) /M (cid:63) and ˜ ω n is the dimensionless mode frequency. ©000 , 000–000 Passamonti, Andersson and Pnigouras
D / (10 km) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 E / E f g g g g i s s s model G = 2.05 s b = 0.1 s b = 0.04 D / ( 10 km ) -14 -12 -10 -8 -6 -4 -2 E / E g g g g f i model G = 7/3g s s s b = 0.1 s b = 0.04 Figure 9.
The evolution of the mode energies for varying binary separation. The energy is given in units of E = GM (cid:63) /R . The modelwith Γ = 2 .
05 is shown in the left-hand panel, while the Γ = 7 / E b of the crust is represented by a circle or a square(see legend), which, respectively, denote the breaking energy for two choices of the breaking strain, ¯ σ b = 0 . σ b = 0 . For each oscillation mode, we can determine the kinetic and potential energy from E k ( t ) = 12 (cid:90) ρ ∂ξ ∗ i ∂t ∂ξ i ∂t dV = 12 (cid:88) n A n | ˙ a n ( t ) | , (81) E p ( t ) = 12 (cid:90) ( ξ ∗ i Cξ i ) dV = 12 (cid:88) n A n ω n | a n ( t ) | , (82)where we have used (cid:104) ξ n , Cξ n (cid:105) = ω n (cid:104) ξ n , ρξ n (cid:105) = A n ω n . (83)The total tidal energy is therefore given by E = E k + E p = 12 (cid:88) n A n (cid:0) | ˙ a n ( t ) | + ω n | a n ( t ) | (cid:1) (84)and the rate of energy transfer from the tide to the oscillation modes is obtained from (Lai 1994)˙ E = − (cid:90) ρ ∂ξ i ∂t ∇ χ ∗ i dV = (cid:88) n (cid:18) GM (cid:48) W lm Q n D l +1 e im Φ( t ) (cid:19) ˙ a n ( t ) . (85)Differently from Lai (1994), we consider only the energy of a single mode and not the couple of m = ± l = 2).The maximum mode energy just after the resonance can be estimated as˜ E max (cid:39) π Q n ˜ A n ˜ ω / n (cid:18) Rc GM (cid:63) (cid:19) / q (cid:18)
21 + q (cid:19) / . (86)In practice, we find that equation (86) provides values which are about 25% smaller than the maximum energy determinedfrom numerical solutions (see Yu and Weinberg 2017, for a similar result).We have quantified the mode resonances for an equal-mass binary system with M (cid:63) = 1 . M (cid:12) , and R = 10 km. In Figure9 we show the evolution of the mode energy during the inspiral for the two stratified models with Γ = 2 .
05 and Γ = 7 / D = 2 R . The first modes to be excited are low-frequency modes, like the surface g modes (notshown in the figure), the interface modes and high-overtone g modes. As the inspiral proceeds, the orbital angular velocityincreases and low-order g modes and shear modes are progressively satisfying the resonant condition. However, it is clear fromthe results that the f mode always dominates the dynamical tide, even though it does not become resonant until after merger(in this example). The energy of other resonant modes does not at any point reach above about 1% of the f-mode energy.This accords well with the discussion of Andersson and Pnigouras (2020a), and supports the assumptions made by Anderssonand Pnigouras (2020b). The maximum g-mode energy grows larger as the stratification becomes stronger, but the resonanceoccurs at a later time, closer to the merger. It is worth noting that, after the resonance, the modes keep oscillating due to theabsence of dissipative processes in our model. This effect would be (at least to some extent) suppressed if we were to includeviscous damping (see Lai 1994). © , 000–000 mpact of the crust D / (10 km) -12 -11 -10 -9 -8 -7 -6 E / E i i i model G = 2 ~ m = 10 -4 ~ m = 5 x10 -4 m = 10 -3 ~ s b = 0.1 s b = 0.04 Figure 10.
Illustrating the dependence of the core-crust interface (i ) mode resonance excitation on the shear modulus of the crust.The results are for the barotropic model (Γ = 2). The three curves show the results for three different values of the shear modulus:˜ µ = 10 − , 5 × − and 10 − (see legend). The shear modulus is given in units of GM (cid:63) /R and the energy in units of E = GM (cid:63) /R .The breaking energy limit E b of the crust is represented by a circle or a square (see legend), which, respectively, denote the breakingenergy for two choices of the breaking strain, ¯ σ b = 0 . σ b = 0 . For the i interface mode the maximum energy is very small, although it increases with stratification. It is about an orderof magnitude larger for the Γ = 7 / = 2 .
05 case. Having a lower frequency compared to the othermodes (see Figure 9) the interface mode enters resonance earlier, when the orbital separation is about 320 km for the modelwith Γ = 2 .
05 and 230 km for the model with stronger stratification. In this earlier phase, the binary evolution is slower andas a result the interface modes have more time to accumulate energy. This will be important when we consider the issue ofcrust failure. Before we consider this question, it is useful to assess the impact of the shear modulus on the interface modes.For the barotropic model (Γ = 2), we determine the i mode—which originates at the crust-core interface—for various valuesof the shear modulus parameter ˜ µ , respectively, ˜ µ = 10 − , 5 × − and 10 − (see Equation (73)). The results are shown inFigure 10, from which it is evident that that the i mode depends strongly on the crust rigidity. Let us now quantify the resonant mode amplitudes relative to the level required to exceed the breaking strain of the crust.The mode energy is not enough to answer this question, we also need to evaluate the elastic strain associated with the modeand this depends on the detailed eigenfunctions.The crust problem is complex (including aspects that are difficult to pin down, like the impact of possible pasta regionsclose to the core-crust transition), but we can make progress by combining the mode eigenfunctions, an estimate for thebreaking strain and the standard von Mises criterion. We first of all need the elastic strain tensor. Hence, we define the tensorfield ¯ σ ij = σ ij / ˇ µ , such that ¯ σ ij = ( ∇ i ξ j + ∇ j ξ i ) − g ij (cid:16) ∇ k ξ k (cid:17) . (87)Through the von Mises criterion, the crust breaking is then established by comparing¯ σ ≡ (cid:114)
12 ¯ σ ∗ ij ¯ σ ij . (88)to the breaking strain. Based on the molecural dynamics simulations of Horowitz and Kadau (2009) the breaking strain iscommonly taken to be ¯ σ b = 0 .
1, so we naturally focus on this case. At the same time it is useful to ask how the results dependon the this assumption. Hence, we also consider the results from Baiko and Chugunov (2018) which suggest the slightly lowervalue of ¯ σ b = 0 . l = m = 2 mode:¯ σ = 5 π sin θr (cid:34) (cid:18) dWdr (cid:19) − (cid:18) dVdr (cid:19) + (cid:18) V − W (cid:19) r dWdr + (cid:18) V − W (cid:19) r dVdr + 14 (cid:18) Wr (cid:19) (cid:35) (89)From this expression we can then determine, for each mode, the oscillation energy required to satisfy the ¯ σ = ¯ σ b condition at ©000
1, so we naturally focus on this case. At the same time it is useful to ask how the results dependon the this assumption. Hence, we also consider the results from Baiko and Chugunov (2018) which suggest the slightly lowervalue of ¯ σ b = 0 . l = m = 2 mode:¯ σ = 5 π sin θr (cid:34) (cid:18) dWdr (cid:19) − (cid:18) dVdr (cid:19) + (cid:18) V − W (cid:19) r dWdr + (cid:18) V − W (cid:19) r dVdr + 14 (cid:18) Wr (cid:19) (cid:35) (89)From this expression we can then determine, for each mode, the oscillation energy required to satisfy the ¯ σ = ¯ σ b condition at ©000 , 000–000 Passamonti, Andersson and Pnigouras f mode cc co r/R cc co r / R g mode cc co r/R cc co r / R i mode cc co r/R cc co r / R Figure 11.
The strain field ¯ σ for the resonant modes which reach the breaking limit, ¯ σ b , during inspiral. For the stratified model withΓ = 7 /
3, we show, from the left to the right panel, meridional 2D cross sections for the f, g and i modes. The f mode can reach thebreaking limit ¯ σ b = 0 .
1, while the g and i modes can fracture the crust only if we consider the lower breaking limit ¯ σ b = 0 .
04 (see thebar legend). Lighter colours indicate larger strain. some point in the crust. For each mode we have considered, and the three stellar models, we report in Tables 4, 5 and 6 thebreaking energy in the crust, E b , and the maximum energy reached during an inspiral. To determine the energy for the case¯ σ b = 0 .
04, one can easily rescale the results given for ¯ σ b = 0 . E b | ¯ σ =0 . = 0 . E b | ¯ σ =0 . .The results for the crust fracturing show, first of all, that the f mode (which dominates the tidal response, see Figure9) may break the crust when the system is close to merger, roughly at a separation of D (cid:39) −
70 km (depending on thechosen value for ¯ σ b ). For the weakly stratified model, we find that the interface i mode reaches an energy slightly lowerthan the breaking limit. In contrast, for the strongly stratified model, the g and i modes both reach an energy above E b and hence may impact on the crust. The former mode reaches the breaking amplitude during the late stages of the inspiral,even later than the f mode. Perhaps more interesting, in this respect, is the i mode which may break the crust at a muchlarger separation, D (cid:39)
230 km. As discussed in Section 3.1, the g modes and the core-crust interface mode have similareigenfunctions in the stratified models. To study the impact of the i mode on the crust, we therefore consider the barotropicstar and explore the effect of the shear modulus on the mode dynamics. For the three values of the shear modulus consideredin Section 4.1, the i mode overcomes the breaking energy E b only in models with ˜ µ ≥ × − .Crust failure due to the resonance of the i mode has been previously studied by Tsang et al (2012) and Pan et al (2020).In general, we find that the maximum i mode energy is close to the crust breaking limit, but the crust only fractures forsome of the models. In our models the overlap integral for the i mode is about an order of magnitude smaller than thosereported in Tsang et al (2012) and Pan et al (2020). This difference is likely due to the different equation of state and shearmodulus prescription. Tsang et al (2012) and Pan et al (2020) use the Newtonian perturbation equations, as in our work,but the stellar model is described by tabulated equations of state. In particular, in Tsang et al (2012) the equilibrium staris determined by using the relativistic structure equations. As the i mode eigenfunctions depend on the properties of thestellar models, as we have demonstrated, it is not surprising that we arrive at slightly different results. The main implicationis clear. If we want to draw firm conclusions on the likelihood of crust failure due to the interface mode, we need to use amore realistic model and a complete relativistic formulation of the problem.Turning to the question of the location at which the crust first fails, we note that the largest stress is reached at theequator for all the modes we consider. The exact location where a mode reaches the breaking limit is indicated in Figure11. The f mode first reaches the threshold for crust failure at r (cid:39) . R , in the low density region close to the surface. Incontrast, the g mode stresses the crust predominantly at the crust-core and crust-ocean interfaces, while the i mode strainreaches its maximum value at the crust-ocean transition, as anticipated given the distinct cusp at the transition density. Forthe barotropic model (Γ = 2), the strain tensor ¯ σ ij for the i mode has a larger magnitude at the top of the crust but reachesa significant level throughout the equatorial plane. This property is more pronounced in the ˜ µ = 10 − case, and may besignificant as it could indicate a global, rather than local, crust failure. We have studied the tidal response of a binary neutron star during the inspiral phase, considering spherically symmetric modelswith crust and ocean. First, we revisited the theoretical formulation of the problem to understand whether the presence ofa crust or density discontinuities require changes to the formalism so far developed. Our analysis shows that these newingredients do not change substantially the Newtonian formalism used for fluid models as long as the fundamental boundaryand junction conditions, at the different interfaces, are satisfied.Considering the various oscillation modes sustained by our stellar models, we have extended the previous analysis of © , 000–000 mpact of the crust Table 4.
Mode excitation and breaking energy for the stratified model with Γ = 2 .
05. We provide the maximum resonant energy E max (second column) and the breaking energy E b (third column; normalised to E = GM (cid:63) /R ) for each mode (first column) we haveconsidered. The breaking energy is determined from the von Mises criterion for ¯ σ b = 0 . E max /E E b /E ¯ σ b = 0 .
1f 4 .
28 1 . × − g . × − . × − g . × − . × − g . × − . × − g . × − . × − g s . × − . × − g s . × − . × − g s . × − . × − i . × − . × − s . × − . × − s . × − . × − s . × − . × − Table 5.
The same as Table 4, but for the model with Γ = 7 / E max /E E b /E ¯ σ b = 0 .
1f 4 .
04 1 . × − g . × − . × − g . × − . × − g . × − . × − g . × − . × − g . × − . × − g s . × − . × − g s . × − . × − g s . × − . × − i . × − . × − s . × − . × − s . × − . × − s . × − . × − Andersson and Pnigouras (2020a) focusing on the effect of the crust. In particular we have studied the Love number in thestatic limit, considering the contribution from the most relevant oscillation modes. We have shown that the presence of thecrust does not significantly affect the fundamental, pressure and gravity modes. As expected, the contribution to the Lovenumber from the shear modes is negligible, while the interface and surface gravity modes have an impact similar to the firstcore gravity modes. The influence of these modes, albeit small compared to the fundamental mode, increases for stronglystratified models.Oscillations may be amplified by tidal resonances during the binary inspiral. This amplification is not only importantfor the gravitational wave signal, but also for the impact that mode resonances can have on the crust. We have studiedthe dynamical tidal evolution and determined the mode energy during the orbital shrinking. Our results confirm that thefundamental mode dominates the dynamical tides even when it is far from resonance. In our models, the f mode would enter
Table 6.
Mode excitation and breaking energy for the barotropic model (Γ = 2) and the i interface mode. We provide the maximumresonant energy E max (second column) and the breaking energy E b (third column; normalised to E = GM (cid:63) /R ) for three different valuesof the shear modulus parameter ˇ µ (first column). The breaking energy is determined from the von Mises criterion for ¯ σ b = 0 . µ E max /E E b /E ¯ σ b = 0 . − . × − . × − × − . × − . × − − . × − . × − ©000
Mode excitation and breaking energy for the barotropic model (Γ = 2) and the i interface mode. We provide the maximumresonant energy E max (second column) and the breaking energy E b (third column; normalised to E = GM (cid:63) /R ) for three different valuesof the shear modulus parameter ˇ µ (first column). The breaking energy is determined from the von Mises criterion for ¯ σ b = 0 . µ E max /E E b /E ¯ σ b = 0 . − . × − . × − × − . × − . × − − . × − . × − ©000 , 000–000 Passamonti, Andersson and Pnigouras resonance for a binary separation
D/R (cid:39) .
74, i.e. after merger. Among the other modes, the first gravity mode reachesthe largest oscillation energy during the late stages of the inspiral, roughly when
D/R (cid:39) −
10 (depending on the stellarmodel). The interface modes are resonantly excited at an earlier stage and may accumulate enough energy to fracture thecrust. This is mainly due to their peaked radial displacement at the crust-core or crust-ocean surface transition. We have usedthe von Mises criterion to determine the minimum energy required to fracture the crust and compare the result to the energygained by a given mode during the binary evolution. The interface modes do not break the crust for all our models. Stronglystratified cases are favoured in this respect. This is not surprising, because it is well known that the interface mode propertiesdepend on the equation of state, shear modulus, density discontinuities, etcetera. A variation in these quantities can lead todifferent conclusions. We also found that the interface modes mainly break the crust at the equator and predominantly at thecrust-ocean transition.The fundamental mode reaches the crust breaking limit in all our models, but not until the final part of the inspiral. Thef-mode eigenfunctions have a more regular behaviour at the crust boundaries than the interface modes, therefore it needs toreach a larger energy in order to fracture the crust. We find that the strain tensor for the f mode reaches its largest valuearound the middle of the crust. Finally, we have shown that the first gravity mode can fracture the crust only for stronglystratified stars and (again) in the very final phase of inspiral.A natural extension of this work would be the inclusion of superfluid and superconducting constituents in the core andthe inner crust. For such models we expect that shear and gravity modes will be shifted towards higher frequencies (see forinstance Passamonti and Andersson 2012; Yu and Weinberg 2017), but the problem is complicated as superfluid entrainmentcomes into play (and may have a particularly large effect in the crust). The corresponding mode resonances would then beexpected very close to the merger and might have negligible impact on the tidal problem. However, in order to quantifythe effect, we need to add the superfluid degree of freedom to the analysis. We have already worked through the formalaspects of this and expect to complete the analysis with a sample of numerical results before too long. Another importantstep would be the development of a relativistic formulation of the problem. This would allow us to study the crucial influenceof realistic/tabulated equations of state and relativistic effects on the problem.
ACKNOWLEDGEMENTS
N.A. and P.P. acknowledge support from the Science and Technology Facilities Council (STFC) via grant ST/R00045X/1.P.P. acknowledges support from the MIUR PRIN 2017 programme (CUP: B88D19001440001) and from the Amaldi ResearchCenter funded by the MIUR programme “Dipartimento di Eccellenza” (CUP: B81I18001170001).
DATA AVAILABILITY
Additional data related to this article will be shared on reasonable request to the corresponding author.
REFERENCES
Abbott, B.P. et al, 2017, Phys. Rev. Lett. 119, 161101Abbott, B.P. et al, 2018, Phys. Rev. Lett., 121, 161101Abbott, B.P. et al, 2019, Phys. Rev. X, 9, 011001Andersson, N., and K.D. Kokkotas, 1998, MNRAS, 299, 1059Andersson, N., G.L. Comer and K. Grosart, 2004, MNRAS, 355, 918Andersson, N., B. Haskell and L. Samuelsson, 2011, MNRAS 416, 118Andersson, N., and W.C.G. Ho, 2018, Phys. Rev. D, 97, 023016Andersson, N., and P. Pnigouras, 2020, Phys. Rev. D, 101, 083001Andersson, N., and P. Pnigouras, 2020,
The phenomenology of dynamical neutron star tides submitted to MNRASBaiko, D.A., and A.I. Chugunov, 2018, MNRAS, 480, 5511Douchin, F., and P. Haensel, 2001, Astron. Astrop., 380, 151Flanagan, E.E., and T. Hinderer, 2008, Phys. Rev. D, 77, 021502Friedman, J.L. and B.F. Schutz, 1978, Ap. J. 221, 937Gittins, F., N. Andersson and J.P. Pereira, 2020, Phys. Re. D, 101, 103025Gittins, F., N. Andersson and D.I. Jones,
Modelling neutron star mountains , to appear in MNRASHinderer, T., B.D. Lackey, R.N. Lang and J.S. Read, 2010, Phys. Rev. D, 81, 123016Hinderer, T., et al, 2016, Phys. Rev. Lett. 116, 181101Horowitz, C.J., and K. Kadau, 2009, Phys. Rev. Lett. 102, 191102 © , 000–000 mpact of the crust Kokkotas, K.D., and G. Schaefer, 1995, MNRAS, 275, 301Kr¨uger, C.J., W.C.G. Ho and N. Andersson, 2015, Phys. Rev. D, 92, 063009Lai, D., 1994, MNRAS 270, 611McDermott, P.N., C.J. Hansen, H.M. van Horn and R. Buland, 1985, Ap. J. Lett., 297, L37McDermott, P.N., H.M. van Horn and C.J. Hansen, 1988, Ap. J., 325, 725Passamonti, A., and N. Andersson, MNRAS 419, 638Poisson, E., 2020, Phys. Rev. D, 102, 064059Miller, M.C. et al, 2019, Ap. J. Lett., 887, 24Pan, Z. et al,
Probing Crust Meltdown in Inspiraling Binary Neutron Stars , preprint arXiv:2003.03330Penner, A.J. et al, 2012, Ap. J. Lett., 749, L36Reisenegger, A., and P. Goldreich, 1994, Ap. J., 426, 688Riley, T.E. et al, 2019, Ap. J. Lett. 887, 21Steinhoff, J., T. Hinderer, A. Buonanno and A. Taracchini, 2016, Phys. Rev. D, 94, 104028Tsang, D. et al, 2012, Phys. Rev. Lett. 108 011102Yu, H., and N. Weinberg, 2017, MNRAS, 464, 2622 ©000
Probing Crust Meltdown in Inspiraling Binary Neutron Stars , preprint arXiv:2003.03330Penner, A.J. et al, 2012, Ap. J. Lett., 749, L36Reisenegger, A., and P. Goldreich, 1994, Ap. J., 426, 688Riley, T.E. et al, 2019, Ap. J. Lett. 887, 21Steinhoff, J., T. Hinderer, A. Buonanno and A. Taracchini, 2016, Phys. Rev. D, 94, 104028Tsang, D. et al, 2012, Phys. Rev. Lett. 108 011102Yu, H., and N. Weinberg, 2017, MNRAS, 464, 2622 ©000 , 000–000 Passamonti, Andersson and Pnigouras p p p p f g g g g g -1.0-0.50.00.51.0 D Q n l / Q n l x ( |Q | - |Q | ) / |Q | x 100( |Q | - |Q | ) / |Q | x 100( |Q | - |Q | ) / |Q | x 100 i s s s s s s -8-6-4-202468 D Q n l / Q n l x ( |Q | - |Q | ) / |Q | x 100( |Q | - |Q | ) / |Q | x 100( |Q | - |Q | ) / |Q | x 100 Figure A1.
Percentage errors for the calculation of the overlap integral with the three expressions (52), (54) and (74), which are,respectively, denoted as Q , Q and Q in this figure. The results represent the stratified star with Γ = 7 /
3. The various oscillationmodes are indicated on the horizontal axis. The overlap integral calculations agree to better than 1% for most of the oscillation modes.The largest difference, (cid:46)
APPENDIX A: NUMERICAL CODE
We determine the oscillation mode properties by solving the perturbation equations, obtained from the single-fluid limitof the equations given by Passamonti and Andersson (2012) for superfluid stars. The stellar model we consider has threeregions: core, crust and ocean. Therefore, we must impose junction conditions at the origin and the star’s surface, as well asboundary/junction conditions at the crust-core and crust-ocean transitions.We solve the linearised equations as an eigenvalue problem by using both multiple shooting methods and a relaxationapproach. The latter was necessary to increase the accuracy of the calculated overlap integral. Basically, some oscillationmodes have very small overlap integrals, the calculation of which may be subject to numerical errors. This is the case, forinstance, for higher order modes which have many nodes in their eigenfunctions, and for shear and interface modes, whichare mainly present in the crust region. To monitor the numerical errors we determine Q n from Equations (52), (54) and (74).The solutions obtained after the relaxation step are much more accurate which allow us to study high-order oscillation modes.In Figure A1 we show the relative difference between the overlap integrals calculated from these three equations. The resultsagree to better than 1% for all modes, with the exception of the interface mode which can have an error at most of order afew %. To reach accurate results we have used a very high grid resolution with 1 . × points. All other relevant quantities,as for instance mode frequencies and the eigenfunction ratio at the surface ( V /W ) R , agree with the results from the literaturealso for lower resolutions. APPENDIX B: MODE RESONANCE DYNAMICS
In order to quantify the mode excitation during binary inspiral, we need to solve equation (79). We do this by introducing anew variable (this differs slightly from Lai (1994)) a = b e − im Φ( t ) . (B1)In terms of this new variable the amplitude equation takes the form¨ b n − im Ω˙ b n + (cid:16) ω n − m Ω − im ˙Ω (cid:17) b n = GM (cid:48) R ˜ Q n ˜ A n W lm (cid:18) RD ( t ) (cid:19) l +1 . (B2)Introducing the decomposition b = b R + ib I we obtain equations for the real and imaginary parts, b R and b I respectively, ofthe scalar function b . These are¨ b Rn + 2 m Ω˙ b In + (cid:0) ω n − m Ω (cid:1) b R + m ˙Ω b In = GM (cid:48) R ˜ Q n ˜ A n W lm (cid:18) RD ( t ) (cid:19) l +1 , (B3)¨ b In − m Ω˙ b Rn + (cid:0) ω n − m Ω (cid:1) b I − m ˙Ω b Rn = 0 . (B4) © , 000–000 mpact of the crust Following Lai (1994), we determine the initial conditions for b R and b I from the static limit of equation (B3), as the modesare then not resonant. By neglecting the time derivative in equation (B3) we obtain for the real part: b R = GM (cid:48) R ˜ Q l ˜ A n W nl (cid:18) RD ( t ) (cid:19) l +1 ω n − m Ω , (B5)where quantities with the subscript “0” are calculated at t . Inserting equation (B5) into (B4), we have˙ b R (cid:39) (cid:20) − ( l + 1) ˙ DD + 2 m Ω ˙Ω ω n − m Ω (cid:21) b R (cid:39) − (cid:20) ( l + 1) + 3 m Ω ω n − m Ω (cid:21) ˙ DD b R , (B6)where we have neglected ¨ b In in equation (B4) and used ˙ΩΩ = −
32 ˙
DD . (B7)For the imaginary part we can choose b I (cid:39) mω n − m Ω (cid:16) b R + ˙Ω b R (cid:17) , (B8)˙ b I (cid:39) . (B9)This then allows us to solve the problem, leading to the results presented in Section 4. ©000