Dynamical tides in neutron stars: The impact of the crust

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 Newtonian gravity, we add the elastic aspects to the theoretical formulation of the problem and quantify the dynamical excitation of different classes of oscillation modes. The results demonstrate the expectation that the fundamental mode dominates the tidal response and show how the usual tidal deformability (and the Love number) emerge in the static limit. In addition, we consider to what extent the different modes may be excited 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 late stages of inspiral. This is also the case for the first gravity mode, which reaches the breaking threshold in strongly stratified stars. In our models with a fluid ocean, interface 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.


INTRODUCTION AND SCOPE
The breakthrough detections of gravitational waves from binary neutron star inspirals (Abbott et al 2017a(Abbott et al , 2018(Abbott et al , 2019 have led to renewed focus on the elusive neutron star equation of state. The problem has a number of complicating aspects-both relating to the observational data and the theoretical underpinning-but the essential question is quite simple: To what extent can we use observations to constrain the state and composition of matter under the extreme conditions that neutron stars represent? Much of the recent focus has been on the neutron star tidal deformability, essentially the extent to which the tidal interaction 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, the celebrated 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, as this radius range agrees well with the constraints obtained from the NICER observations of PSR J0030+0451 (Miller et al 2019;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 the neutron star's tidal response. In addition, there is a dynamical tide. This is traditionally represented by the excitation of the different oscillation modes of the star. The resonance problem was first considered some time ago (Lai 1994;Reisenegger and Goldreich 1994;Kokkotas and Schaefer 1995;Andersson and Ho 2018), but the issue is back in focus following the suggestion that the (fundamental) f-mode of the star may be excited to a relevant level, even though it may not reach resonance during the inspiral Steinhoff et al 2016;Andersson and Pnigouras 2020b). The associated effect on the inspiral signal 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 little bit more complicated than a prescribed pressure-density relation. Nuclei in the lower density region are expected to freeze to form the so-called crust, neutrons and protons likely form superfluid/superconducting condensates at high densities, there may be phase transitions to states with net strangeness (involving hyperons or deconfined quarks) and so on. These issues are important because the tidal response can be expressed as a sum over the (presumably complete set of) stellar oscillation modes (Lai 1994;Reisenegger and Goldreich 1994;Kokkotas and Schaefer 1995) and the additional features of the interior physics 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 (Andersson and Pnigouras 2020a) demonstrates that variations in the composition of the neutron star core affects the dynamical tide at the few percent level. This discussion also provides the (conceptually important) link between the static tide and the star's oscillation modes. Another important step was taken by Yu and Weinberg (2017) who considered the impact of mode resonances 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 most detailed seismology models to date was developed by Krüger, Andersson and Ho (2015)). However, the framework required to quantify 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, we may progress using phenomenological Newtonian models. In this direction, perhaps the most "complete" effort is provided by the results of Passamonti and Andersson (2012) (see also Andersson, Haskell and Samuelsson (2011)), which represent a star with an elastic crust penetrated by superfluid neutrons as well as a superfluid core (taking into account the entrainment effect in both regions) and an ocean. A reasonable aim then would be to consider the tidal response of such a star. This involves two steps. First, we need to develop the mode-coupling framework for a star with an elastic crust and ocean, and ensure that the mode calculation is sufficiently precise that we can quantify the (presumably) weak tidal excitation of (say) the crustal shear 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 consider the issue of tidally induced crust fracturing (Tsang et al 2012;Pan et al 2020)-to what extent resonant oscillation modes may reach the amplitude required to exceed the breaking strain at some point in the crust.

FORMULATING THE PROBLEM
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 the primary. If we want to quantify the fluid response to this external agent-in essence, model the dynamical tide-we need to solve the linearised fluid equations in Newtonian gravity. As we are dealing with an (at least partly) elastic body, it is natural to do this in the framework of Lagrangian perturbation theory (Friedman and Schutz 1978). In fact, many of the formalities have already been dealt with by Andersson, Haskell and Samuelsson (2011) and Passamonti and Andersson (2012).

The perturbation equations
Assuming that the star is non-rotating, which makes sense on astrophysical grounds as the two partners in a binary would have had plenty of time to spin down before the system enters the sensitivity band of a ground-based gravitational-wave detector (although, see the recent analysis of Poisson (2020) for effects associated with spinning stars), we first of all have the perturbed continuity equation 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) (with δ the corresponding Eulerian variation and L ξ the Lie derivative along ξ i ), such that the perturbed velocity is given by Provided the background is in hydrostatic equilibrium, the perturbed Euler equation is and we also have the Poisson equation for the gravitational potential while the tidal potential, χ, due to the presence of the binary partner (which generates the fluid perturbation), is a solution to ∇ 2 χ = 0. As we progress, it is useful to note that (1) leads to which allows us to write the Poisson equation as This will be useful later, as we can integrate to get for some vector S i such that ∇iS i = 0. It is easy to see that, if ρ → 0 at the star's surface, then we must have and we know that this should not vanish. Moreover, we know that the solution has to match to the vacuum exterior, which means that (again, as long as 1 the density vanishes as r → R) d dr where we have expanded δΦ in spherical harmonics and focussed on a single l multipole. That is, we must havê whereni is the normal to the star's spherical surface (i.e. parallel to the radial basis vector). We will have reason to recall this result later.

Mode orthogonality
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 ξ i n (with n labelling the modes). Letting the (real) mode frequency be ωn, each individual mode then satisfies the homogeneous version 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) where µ = mBμ is the chemical potential and ρ = mBn, with mB the mass of each baryon and n the baryon number density. This means that. for a background in hydrostatic equilibrium, we have 1 ρ ∇iδp − 1 ρ 2 δρ∇ip = ∇iδμ .
As a result, we have or, if we take the baryon number density to be the primary matter variable; Following Friedman and Schutz (1978) we now write the Euler equation as where and and introduce the inner product where the asterisk indicates the complex conjugate and η i is another solution to the perturbation equations. It is then fairly easy to show that mode solutions are orthogonal if the mode frequencies are real. As we will need to extend the proof of this to the elastic case, let us have a look at the argument-taking the opportunity to keep careful track of "surface terms" that come into play if we allow internal phase transitions. First of all, it is obvious from the definition of the inner product that The argument for the C operator is a bit more involved. We need to use integration by parts to show that whereni is the (outwards pointing) normal to the surface of the volume we are integrating over. The last equality holds if can ignore the contribution from the surface term. This would certainly be the case if we integrate over the entire star, as long as n → 0 at the surface, but there may be situations where one would need to be more careful.
For example, as we may want to consider models with internal phase transitions-e.g. at the crust-core interface-it is important to establish what happens if such a transition is associated with a discontinuity in the density. It is easy to see that (21) implies that, if there is an internal density jump, at r =R, say, we need the surface term to be continuous. To see that this should be the case, use (again for a barotrope) dp = ndµ to get 1 mB where δ ξ and δη distinguishes the Eulerian perturbations associated with ξ i and η i , respectively. Continuity across a surface requires (i) the radial component of the displacement, and (ii) the Lagrangian pressure variation, to be continuous. The last term of equation (23) can be rewritten as 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 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 The first term (obviously) vanishes as long as ρ → 0 at the surface. Moreover, noting that the matching to the exterior potential requires (10) to hold, we see that the second term vanishes, as well. Specifically, we have (for each multipole) The conditions that apply when the density does not vanish at the surface, can be inferred from the results for an internal density jump. Again, assume that there is a density jump at an internal point, r =R, such that (in a small neighbourhood of the phase transition) and integrate (7) over a small volume (with height = 2 ) across the (spherical) surface to get for some constant S. Then let → 0 to get the junction condition where we have used the fact that the radial component of the displacement must be continuous. What does this mean for the surface terms in (25)? Well, we integrate up toR and then continue on the other side. The surface terms associated with the interface will vanish as long aŝ We can rewrite this as Combining the fact that the perturbed gravitational potential is continuous with (30), we see that these contributions must vanish. It is also easy to adjust the argument to arrive at the standard result for the gravitational potential matching at a finite-density surface. In effect, we have provided a (somewhat lengthy) demonstration that we do not have to worry about density discontinuities in the following. As long as the relevant junction conditions are respected, the usual orthogonality argument remains unchanged. At the end of the day, we have the expected result (Friedman and Schutz 1978): This means that, if we consider two mode solutions (now letting ξn be associated with frequency ωn and ξm be the solution associated with ωm) we have We see that, if the frequencies are distinct and real, the mode solutions must be orthogonal. That is, we may use leaving the normalisation A 2 n unspecified for the moment. So far, the arguments are standard. The only aspect we have added relates to internal phase transitions. Before we move on let us, without detailed calculation, comment on the case of non-barotropic perturbations, relevant whenever we want to account for internal matter stratification. It is easy to argue that this case applies to neutron star tides as the relevant nuclear reactions are too slow to keep the perturbed matter in beta equilibrium (Andersson and Pnigouras 2020a). In particular, it leads 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 choice is to complement the density ρ with either parameter representing the deviation from chemical equilibrium or the proton fraction xp = ρp/ρ, the Lagrangian perturbation of which will vanish for slow reactions. However, the mode orthogonality can 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 This is the result we need to prove the mode orthogonality for stratified stars. As the argument for the gravitational potential remains unchanged, the stated symmetry relation (33) will hold as long as the radial displacement and the pressure perturbation are both continuous (for a spherical surface).

Adding the crust
We now take the "natural" next step towards a more realistic neutron star model, adding the elastic crust that is expected to reach from (close to) the surface up to about 60-70% of the nuclear saturation density (roughly, the outer kilometer of the star). Assuming that the crust of the equilibrium star is "relaxed"-not associated with strains, the elasticity impacts only 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 where with the elastic stress tensor given by Note that the shear modulusμ must not to be confused with the chemical potential. Now, we need The second term has two pieces: and Introducingσ we see that In order to deal with the surface terms we need to consider the traction conditions (Andersson, Haskell and Samuelsson 2011). Hence, we impose the continuity of the gravitational potential and the perpendicular and radial components of the traction vector noting that the traction reduces to the pressure perturbation in the fluid regions. In order to establish that the surface terms vanish at the crust-core transition, we combine the above result with the relevant term from (21). This leads to the requirement that 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.

Expanding the tidal response
Finally, we are set to return to the tidal problem. This part of the analysis is straightforward given the properties we have established. Perhaps not surprisingly, the results we need are identical to those from the fluid case, see, for example, Andersson and Pnigouras (2020a). In particular, it is clear that the oscillation modes of the fluid+elastic problem are orthogonal and that (35) still holds. As in the fluid problem, we can use this fact to rewrite the Euler equation as an equation for individual mode amplitudes. Introducing the mode expansion where the eigenfunctions ξ i n are time-independent (as they are obtained in the frequency domain), we easily arrive aẗ 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 where the v l coefficients will be given later. Assuming an adiabatic inspiral and working in the frequency domain (with time dependence e iωt ), the tidal driving only has support at frequency ω = mΩ. It is useful to keep this in mind. Anyway, we arrive at a set of driven modes with amplitude where we have introduced the overlap integral As discussed by Andersson and Pnigouras (2020a) the overlap integral can also be expressed in terms of the matching of the gravitational potential at the star's surface so we have where In 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 a representation for the perturbed gravitational potential in terms of the mode eigenfunctions. Expressing the displacement vector for each multipole as it follows from the θ-component of (4) that (since ∆p = 0 at the surface) That is, we have and we arrive at an expression for the dynamical tidal response Moreover, making contact with (51) we see that which means that we have (simply adding a label, n, to indicate the individual mode eigenfunctions) with M being the primary's mass. Keeping the normalisation of the modes unspecified-which may be useful when we consider the numerical aspects of the problem-we have Introducing the dimensionless frequencyω as well as the dimensionless overlap integralQ we get This is the final result, but we can massage it a bit by recalling (54) and using We have (Andersson and Pnigouras 2020a) which, when combined with (64), leads to the final expression  In the low-frequency limit, we obtain a mode sum for the tidal Love number The careful reader will note that the expressions for the effective Love number remain unchanged from Andersson and Pnigouras (2020a). This is to be expected, as long as we consider a neutron star with a fluid (ocean) surface. The results would change if we were to impose a direct transition from the elastic region to the exterior vacuum as we would then have to consider elastic terms in, for example, (56). A fluid ocean is expected on physical grounds (even though it may be very shallow 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.

THE DYNAMICAL TIDE
Having discussed the theoretical framework, let us move on to consider the numerical results for a sample of neutron star models. In order to facilitate a direct comparison, we have chosen to focus on the stratified polytropic models already considered by 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 elasticity and interface modes linked to the core-crust and crust-ocean transitions (McDermott et al 1985;McDermott, van Horn and Hansen 1988). We do not consider the impact of superfluidity in the core or the crust at this point, but plan to return to this problem later.
In the first instance, we are interested in the static and dynamical aspects of the tide, as represented (via the different oscillation modes) by the tidal deformability and the Love number k eff l . As we anticipate the impact of the elasticity on already existing fluid modes to be slight and the tidal excitation of, for example, the shear oscillations in the crust to be weak, we have to make sure our numerical approach is robust. We need to do better than the proof-of-principle analysis of Andersson and 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 were obtained using a "classic" shoot-and-relax approach, where a shooting method was used to identify each oscillation mode and the precision was subsequently improved by a relaxation step. That this implementation allows us to reliably extract the quantities 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 of stratification on the oscillation modes and the tide, we introduce an adiabatic index, Γ1, for the perturbations: Following Andersson and Pnigouras (2020a) we consider three models, a barotropic star with Γ1 = 2 and weakly and strongly stratified stars with, respectively, Γ1 = 2.05 and Γ1 = 7/3. For simplicity, we assume that this equation of state describes  the entire star, from the core to the ocean. This may not be particularly realistic, but as we are working in the framework of Newtonian gravity we have to settle for a somewhat phenomenological set-up. From equation (69) it follows that we can write the relation between the Eulerian perturbations of pressure and mass density δρ where 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 Rcc = 0.9R and the crust-ocean interface at Rco = 0.999R. 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. This would mean takingμ but in practice we use the (further) simplified version from Passamonti and Andersson (2012) and leť whereμ is a dimensionless parameter, which must not be confused with the chemical potential defined in equation (12). For a star with M = 1.4M and R = 10 km we then find GM /R = 1.86 × 10 20 cm 2 /s 2 . Linking to (72) we assume the valuẽ µ = 10 −4 for most of our examples, exploring other cases only for the interface modes.

The oscillation modes
Even though it is fairly simple, the stellar model we consider sustains many (more or less) distinguishable classes of oscillation modes: the fundamental mode (f mode), the pressure modes (p modes), the gravity modes (g modes), the shear modes (s modes) and the interface modes (i modes) (see McDermott et al (1985) and McDermott, van Horn and Hansen (1988) for detailed discussions). The f and p modes are present also in a fluid star, while the g modes appear in stratified models (in our case when Γ1 = Γ), see Andersson and Pnigouras (2020a). The presence of the crust impacts on these classes of modes and introduces, due to the elasticity, the shear modes. Moreover, any transition between different regions in the star (core-crust and 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 on the gravitational-wave signal, we focus our attention on the quadrupole (l = 2) modes. We then have the sets of oscillation frequencies provided (in dimensionless units) in Tables 1, 2 and 3. The tabulated data include other relevant quantities required for the tidal problem (e.g. the overlap integrals, the ratio at the surface between the tangential and radial components of  Figure 3. Eigenfunctions for the i 1 mode, which is associated with the crust-ocean transition (note the distinctive cusp in the radial displacement). The black curves represent the radial component W = W/(rR) while the red curves are for the tangential component The left-hand panel shows the results for a barotropic model (Γ 1 = 2). The middle panel represents the Γ 1 = 2.05 case, while the right-hand panel is for a stratified model with Γ 1 = 7/3.  the Lagrangian displacement and the mode contribution to the Love number, all discussed in Section 2.4). The pn modes (with n labelling each mode) cover the high frequency range of the spectrum-above the frequency of the f mode-and their frequencies increase for higher order (n) modes. In contrast, the (gravity) g modes cover the low frequency range of the oscillation spectrum-below the f mode-and their frequencies decrease for higher order modes. As a result, the g-mode spectrum becomes progressively dense towards zero frequency. This makes the identification of the individual high-order g modes challenging. A stratified ocean also sustains a family of g modes, here referred to as surface g modes (g s ). These also have low frequencies and their eigenfunctions are mainly confined to the star's ocean. The shear mode frequency scales approximately as ωn ∼ μ/ρ and increases for higher order modes. Finally, the interface modes are characterised by their low 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 Steinhoff et al 2016;Andersson and Pnigouras 2020a,b). This mode is known to be weakly affected by composition stratification (see the results of Andersson and Pnigouras 2020a) and crust elasticity (McDermott, van Horn and Hansen 1988). This is demonstrated by the results in Figure 1, where we show the l = 2 f-mode eigenfunction for a strongly stratified model with Γ1 = 7/3. The radial component (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 the tangential 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 g4 mode for the Γ1 = 2.05 model, which, as we will see later, has a significant variation in the overlap integral compared to the fluid star. In Figure 2  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 is peaked at the interface. This feature is evident from Figure 3, where we show the displacements for the crust-ocean interface mode (i1) for the barotropic model as well as for the weakly and strongly stratified stars. As anticipated, the W eigenfunction has 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 is discontinuous at the interfaces, and reaches a very large amplitude at the star's surface. Note that the function V has been reduced, in Figure 3, by a factor 10 −3 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, see Figure 3. In this region, the character of the interface modes is very similar to that of the (core) g modes. This behaviour usually 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 (i2), associated with the crust-core transition, is readily determined for the barotropic case (Γ1 = 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 is difficult to numerically separate them. This behaviour is shown in Figure 2 for the g4 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 Γ1 = 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 features associated with a high-order g mode in the core. The crust region appears to "isolate" the two regions which makes it difficult to numerically establish which set of modes a given solution belongs to. It is not even clear that this is a meaningful question for the high overtone modes.

The Love number
The tidal response of a neutron star is closely related to the nature of the different oscillation modes. This is natural since the modes form a complete set and, hence, can be used as a basis to express the behaviour of the stellar fluid. In the static limit the sum over the star's oscillation modes leads to the Love number, as explicitly demonstrated by Andersson and Pnigouras (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 through resonance 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  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 the expectation that the main contribution to the tidal response is provided by the f mode, which has the best overlap with the tidal driving force. This also leads to the question of the level at which other modes, which may depend on the matter composition etcetera, contribute. For example, we know that the elastic crust sustains shear modes. These are expected to have small overlap integrals and therefore have little impact. However, even if these expectations are true, it is useful to quantify what we mean by "small" and to what extent we can safely neglect the contribution of these modes. We need to be mindful of the fact that, even if the contribution from each mode is too small to impact on the gravitational-wave signal, the presence of an excited mode may have other repercussions. For example, it is possible that a mode passing through resonance reaches an amplitude where it leads to fracturing of the crust (Tsang et al 2012;Pan et al 2020). We will return to this possibility later. Finally, we already know from the eigenfunctions we have provided that the presence of the crust affects the properties of all modes (see Section 3.1). This means that there will be a (probably weak, but nevertheless) impact on the respective 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 ω ωn), the tidal Love number may be . Summary results for all the modes we have considered. We show the Love number contribution, k n l , against the mode frequencies (in dimensionless units). The illustrated modes are the fundamental (f) mode, the pressure (pn) and gravity (gn) modes, the shear (sn) 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 n l (see Equation (68)). Hence, we report in Tables 1, 2 and 3 the value of k n l for the different oscillation modes (along with other relevant quantities). The results confirm that the f mode dominates the 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, see Figures 5 and 8), in accordance with the results of Andersson and Pnigouras (2020a). For strongly stratified models, with Γ1 = 7/3, the first g mode has a value for k n l 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 not affected at all, while the gravity modes have a less regular behaviour. As shown in Figure 6, the contribution of the g modes to the Love number, k n l , tends to be smaller for the model with a crust (although it is practically unchanged for the g1 mode). In strongly stratified models, the g modes are less influenced by the crust, and the corresponding values of Qn and k n l are similar to the fluid case, most likely because the buoyancy dominates the elastic restoring force. We also note an interesting result for the g1 mode. In the Γ1 = 2.05 model, the overlap integral, Qn, is about an order of magnitude larger than in the fluid case, while the value for k n l does not change correspondingly (see Figures 5 and 6). This demonstrates that the overlap integral 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 small values for Qn and k n l . In Figure 7 we show the two quantities for the first six s modes. As in the case of the g modes, the behaviour of the overlap integral is less regular with respect to the radial number n than what we find for k n l . Meanwhile, the interface modes linked to the core-crust and crust-ocean transitions may contribute to the tidal response at a level similar to that of the first or second g modes (see Figures 7 and 8). The importance of these modes increases for models with stronger stratification, 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 n l for all modes considered in this work. From this figure we infer which are the most relevant modes in the static limit and in which order the modes will be excited during a binary inspiral (as the driving frequency increases). In this figure we also show k n l 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 overlap integrals, the calculation of which may be subject to numerical errors. This is particularly the case for high-order modes which tend to have many nodes in their eigenfunctions, leading to cancellations in the calculation of the overlap integral. To monitor the numerical errors we determine Qn from equations (52), (54) as well as the (equivalent) expression (Lai 1994) Qn = l ρ [W + (l + 1) V ] r l dr .
As we have already mentioned, in order to increase the precision of the calculation, we solve the perturbation equations first with a multiple shooting method and then with a relaxation step. As expected, the solutions obtained from the relaxation are more accurate, which allows us to extract the high-order modes. A few additional comments on the technical aspects of the calculation are provided in Appendix A. Mode results for the barotropic model with Γ 1 = 2. The first four columns, respectively, report the oscillation mode, the dimensionless mode frequency, the absolute value of the dimensionless overlap integral and the ratio between the radial and tangential components of the displacement vector at the surface. The last column provides the contribution each mode makes to the static Love number, k n l . We also provide, in the final row, the summed result which should be compared to the anticipated value k l = 0.259909.

CRUST FRACTURING
As different oscillation modes pass through resonance during a binary inspiral, their amplitude may become large enough that the motion in the crust induces (local) fracturing of the nuclear lattice. It has been suggested that the interface modes are particularly relevant in this respect (Tsang et al 2012;Pan et al 2020). It also known that the static tide is unlikely to break the crust before the binary merger (Penner et al 2012;Gittins, Andersson and Pereira 2020). In order to consider this problem we need to complement our mode analysis with the energy deposited in each mode during inspiral. This analysis closely follows, for example, Lai (1994) so we will only outline the steps here. Some further details are provided in Appendix B.

The mode excitation
The binary separation, D, shrinks at a rate which can be described, to leading order, bẏ where M is the mass of the companion. This leading-order expression should be sufficient as long as we are not trying to resolve 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 We know from the perturbation analysis that the mode amplitude can be calculated from Equation (49), where we need the tidal potential with Φ = Ω(t)dt and explicitly defining the v l coefficient we used earlier. For l = 2 the W lm coefficients have the following values: For a given mode (l, m), Equation (49) leads tö whereQn is the dimensionless overlap integral defined in equation (63) andÃ 2 n = A 2 n /(M R 2 ). Equation (79) is a forced harmonic oscillator for the mode amplitude an, where the forcing term is provided by the tidal potential. An oscillation mode is in resonance with the orbit when ωn mΩ. For the most relevant modes, the resonance occurs at the late stages of binary inspiral, when the separation changes rapidly and hence the energy transfer to the mode is limited. Using the resonance condition in equation (76) we have the (dimensionless) resonance distance (Lai 1994) where q = M /M andωn is the dimensionless mode frequency. For each oscillation mode, we can determine the kinetic and potential energy from where we have used ξn, Cξn = ω 2 n ξn, ρξn = A 2 n ω 2 n .
The total tidal energy is therefore given by and the rate of energy transfer from the tide to the oscillation modes is obtained from (Lai 1994) Differently from Lai (1994), we consider only the energy of a single mode and not the couple of m = ±2 modes (for l = 2). The maximum mode energy just after the resonance can be estimated as In practice, we find that equation (86) provides values which are about 25% smaller than the maximum energy determined from 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 = 1.4M , and R = 10 km. In Figure  9 we show the evolution of the mode energy during the inspiral for the two stratified models with Γ1 = 2.05 and Γ1 = 7/3. The results show how the various oscillation modes are resonantly excited as the binary separation decreases towards merger, roughly corresponding to D = 2R. The first modes to be excited are low-frequency modes, like the surface g modes (not shown in the figure), the interface modes and high-overtone g modes. As the inspiral proceeds, the orbital angular velocity increases and low-order g modes and shear modes are progressively satisfying the resonant condition. However, it is clear from the 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 Andersson and Pnigouras (2020b). The maximum g-mode energy grows larger as the stratification becomes stronger, but the resonance occurs at a later time, closer to the merger. It is worth noting that, after the resonance, the modes keep oscillating due to the absence of dissipative processes in our model. This effect would be (at least to some extent) suppressed if we were to include viscous damping (see Lai 1994). For the i1 interface mode the maximum energy is very small, although it increases with stratification. It is about an order of magnitude larger for the Γ1 = 7/3 model compared to the Γ1 = 2.05 case. Having a lower frequency compared to the other modes (see Figure 9) the interface mode enters resonance earlier, when the orbital separation is about 320 km for the model with Γ1 = 2.05 and 230 km for the model with stronger stratification. In this earlier phase, the binary evolution is slower and as a result the interface modes have more time to accumulate energy. This will be important when we consider the issue of crust 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 i2 mode-which originates at the crust-core interface-for various values of the shear modulus parameterμ, respectively,μ = 10 −4 , 5 × 10 −4 and 10 −3 (see Equation (73)). The results are shown in Figure 10, from which it is evident that that the i2 mode depends strongly on the crust rigidity.

Breaking the crust
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 mode and 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 regions close to the core-crust transition), but we can make progress by combining the mode eigenfunctions, an estimate for the breaking strain and the standard von Mises criterion. We first of all need the elastic strain tensor. Hence, we define the tensor fieldσij = σij/μ, such thatσ Through the von Mises criterion, the crust breaking is then established by comparinḡ to the breaking strain. Based on the molecural dynamics simulations of Horowitz and Kadau (2009) the breaking strain is commonly 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 depend on the this assumption. Hence, we also consider the results from Baiko and Chugunov (2018) which suggest the slightly lower value ofσ b = 0.04. Introducing the vector expansion of the Lagrangian displacement (55) into equation (88) we obtain the following expression for an l = m = 2 mode:  Figure 11. The strain fieldσ for the resonant modes which reach the breaking limit,σ b , during inspiral. For the stratified model with Γ 1 = 7/3, we show, from the left to the right panel, meridional 2D cross sections for the f, g 1 and i 1 modes. The f mode can reach the breaking limitσ b = 0.1, while the g 1 and i 1 modes can fracture the crust only if we consider the lower breaking limitσ b = 0.04 (see the bar 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 the breaking energy in the crust, E b , and the maximum energy reached during an inspiral. To determine the energy for the casē σ b = 0.04, one can easily rescale the results given forσ b = 0.1 by using E b |σ=0.04 = 0.16 E b |σ=0.1.
The results for the crust fracturing show, first of all, that the f mode (which dominates the tidal response, see Figure  9) may break the crust when the system is close to merger, roughly at a separation of D 50 − 70 km (depending on the chosen value forσ b ). For the weakly stratified model, we find that the interface i1 mode reaches an energy slightly lower than the breaking limit. In contrast, for the strongly stratified model, the g1 and i1 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 i1 mode which may break the crust at a much larger separation, D 230 km. As discussed in Section 3.1, the g modes and the core-crust interface mode have similar eigenfunctions in the stratified models. To study the impact of the i2 mode on the crust, we therefore consider the barotropic star and explore the effect of the shear modulus on the mode dynamics. For the three values of the shear modulus considered in Section 4.1, the i2 mode overcomes the breaking energy E b only in models withμ ≥ 5 × 10 −4 .
Crust failure due to the resonance of the i2 mode has been previously studied by Tsang et al (2012) and Pan et al (2020). In general, we find that the maximum i2 mode energy is close to the crust breaking limit, but the crust only fractures for some of the models. In our models the overlap integral for the i2 mode is about an order of magnitude smaller than those reported in Tsang et al (2012) and Pan et al (2020). This difference is likely due to the different equation of state and shear modulus 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 star is determined by using the relativistic structure equations. As the i2 mode eigenfunctions depend on the properties of the stellar models, as we have demonstrated, it is not surprising that we arrive at slightly different results. The main implication is clear. If we want to draw firm conclusions on the likelihood of crust failure due to the interface mode, we need to use a more 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 the equator for all the modes we consider. The exact location where a mode reaches the breaking limit is indicated in Figure  11. The f mode first reaches the threshold for crust failure at r 0.96R, in the low density region close to the surface. In contrast, the g1 mode stresses the crust predominantly at the crust-core and crust-ocean interfaces, while the i1 mode strain reaches its maximum value at the crust-ocean transition, as anticipated given the distinct cusp at the transition density. For the barotropic model (Γ1 = 2), the strain tensorσij for the i2 mode has a larger magnitude at the top of the crust but reaches a significant level throughout the equatorial plane. This property is more pronounced in theμ = 10 −3 case, and may be significant as it could indicate a global, rather than local, crust failure.

CONCLUDING REMARKS
We have studied the tidal response of a binary neutron star during the inspiral phase, considering spherically symmetric models with crust and ocean. First, we revisited the theoretical formulation of the problem to understand whether the presence of a crust or density discontinuities require changes to the formalism so far developed. Our analysis shows that these new ingredients do not change substantially the Newtonian formalism used for fluid models as long as the fundamental boundary and junction conditions, at the different interfaces, are satisfied.

Mode
Emax 1.63 × 10 −5 3.59 × 10 −5 g 2 2.72 × 10 −7 1.96 × 10 −6 g 3 1.02 × 10 −8 3.71 × 10 −6 g 4 4.85 × 10 −10 3.19 × 10 −6 g 5 2.71 × 10 −11 3.24 × 10 −6 g s 1 4.10 × 10 −14 8.39 × 10 −11 g s 2 3.78 × 10 −15 2.05 × 10 −10 g s 3 6.41 × 10 −23 6.22 × 10 −6 i 1 6.06 × 10 −12 1.24 × 10 −11 s 1 1.63 × 10 −9 1.65 × 10 −8 s 2 9.06 × 10 −10 8.16 × 10 −9 s 3 1.68 × 10 −8 5.64 × 10 −9 Andersson and Pnigouras (2020a) focusing on the effect of the crust. In particular we have studied the Love number in the static limit, considering the contribution from the most relevant oscillation modes. We have shown that the presence of the crust does not significantly affect the fundamental, pressure and gravity modes. As expected, the contribution to the Love number from the shear modes is negligible, while the interface and surface gravity modes have an impact similar to the first core gravity modes. The influence of these modes, albeit small compared to the fundamental mode, increases for strongly stratified models. Oscillations may be amplified by tidal resonances during the binary inspiral. This amplification is not only important for the gravitational wave signal, but also for the impact that mode resonances can have on the crust. We have studied the dynamical tidal evolution and determined the mode energy during the orbital shrinking. Our results confirm that the fundamental 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 (Γ 1 = 2) and the i 2 interface mode. We provide the maximum resonant energy Emax (second column) and the breaking energy E b (third column; normalised to E 0 = GM 2 /R) for three different values of the shear modulus parameterμ (first column). The breaking energy is determined from the von Mises criterion forσ b = 0.1.
10 −4 8.36 × 10 −11 3.99 × 10 −9 5 × 10 −4 5.47 × 10 −8 4.35 × 10 −8 10 −3 8.00 × 10 −7 1.72 × 10 −7 resonance for a binary separation D/R 1.74, i.e. after merger. Among the other modes, the first gravity mode reaches the largest oscillation energy during the late stages of the inspiral, roughly when D/R 5 − 10 (depending on the stellar model). The interface modes are resonantly excited at an earlier stage and may accumulate enough energy to fracture the crust. This is mainly due to their peaked radial displacement at the crust-core or crust-ocean surface transition. We have used the von Mises criterion to determine the minimum energy required to fracture the crust and compare the result to the energy gained by a given mode during the binary evolution. The interface modes do not break the crust for all our models. Strongly stratified cases are favoured in this respect. This is not surprising, because it is well known that the interface mode properties depend on the equation of state, shear modulus, density discontinuities, etcetera. A variation in these quantities can lead to different conclusions. We also found that the interface modes mainly break the crust at the equator and predominantly at the crust-ocean transition.
The fundamental mode reaches the crust breaking limit in all our models, but not until the final part of the inspiral. The f-mode eigenfunctions have a more regular behaviour at the crust boundaries than the interface modes, therefore it needs to reach a larger energy in order to fracture the crust. We find that the strain tensor for the f mode reaches its largest value around the middle of the crust. Finally, we have shown that the first gravity mode can fracture the crust only for strongly stratified 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 and the inner crust. For such models we expect that shear and gravity modes will be shifted towards higher frequencies (see for instance Passamonti and Andersson 2012;Yu and Weinberg 2017), but the problem is complicated as superfluid entrainment comes into play (and may have a particularly large effect in the crust). The corresponding mode resonances would then be expected very close to the merger and might have negligible impact on the tidal problem. However, in order to quantify the effect, we need to add the superfluid degree of freedom to the analysis. We have already worked through the formal aspects of this and expect to complete the analysis with a sample of numerical results before too long. Another important step would be the development of a relativistic formulation of the problem. This would allow us to study the crucial influence of realistic/tabulated equations of state and relativistic effects on the problem.  (52), (54) and (74), which are, respectively, denoted as Q 1 , Q 3 and Q 2 in this figure. The results represent the stratified star with Γ 1 = 7/3. The various oscillation modes are indicated on the horizontal axis. The overlap integral calculations agree to better than 1% for most of the oscillation modes. The largest difference, 8%, is found for the interface modes.

APPENDIX A: NUMERICAL CODE
We determine the oscillation mode properties by solving the perturbation equations, obtained from the single-fluid limit of the equations given by Passamonti and Andersson (2012) for superfluid stars. The stellar model we consider has three regions: core, crust and ocean. Therefore, we must impose junction conditions at the origin and the star's surface, as well as boundary/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 relaxation approach. The latter was necessary to increase the accuracy of the calculated overlap integral. Basically, some oscillation modes have very small overlap integrals, the calculation of which may be subject to numerical errors. This is the case, for instance, for higher order modes which have many nodes in their eigenfunctions, and for shear and interface modes, which are mainly present in the crust region. To monitor the numerical errors we determine Qn 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 results agree to better than 1% for all modes, with the exception of the interface mode which can have an error at most of order a few %. To reach accurate results we have used a very high grid resolution with 1.92 × 10 6 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 literature also 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 a new 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 bn − 2imΩḃn + ω 2 n − m 2 Ω 2 − imΩ bn = . (B2) Introducing the decomposition b = b R + ib I we obtain equations for the real and imaginary parts, b R and b I respectively, of the scalar function b. These arë b R n + 2mΩḃ I n + ω 2 n − m 2 Ω 2 b R + mΩb I n = Following Lai (1994), we determine the initial conditions for b R and b I from the static limit of equation (B3), as the modes are then not resonant. By neglecting the time derivative in equation (B3) we obtain for the real part: where quantities with the subscript "0" are calculated at t0. Inserting equation (B5) into (B4), we havė where we have neglectedb I n in equation (B4) and usedΩ For the imaginary part we can choose This then allows us to solve the problem, leading to the results presented in Section 4.