{'fraction_non_alphanumeric': 0.034145841267981134, 'fraction_numerical': 0.01498549790525298, 'mean_word_length': 4.508068420203324, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Towards real neutron star seismology: Accounting for elasticity and superfluidity

We study the effects of an elastic crust on the oscillation spectrum of superfluid neutron stars. Within the two fluid formalism, we consider Newtonian stellar models that include the relevant constituents of a mature neutron stars. The core is formed by a mixture of superfluid neutrons and a conglomerate of charged particles, while the inner crust is described by a lattice of nuclei permeated by superfluid neutrons. We linearise the Poisson and the conservation equations of nonrotating superfluid stars and study the effects of elasticity, entrainment and composition stratification on the shear and acoustic modes. In both the core and the crust, the entrainment is derived from recent results for the nucleon effective mass. Solving the perturbation equations as an eigenvalue problem, we find that the presence of superfluid neutrons in the crust and their large effective mass may have significant impact on the star's oscillation spectrum.


INTRODUCTION
Neutron stars represent very complex systems, the understanding of which relies on uncertain theory combined with largely indirect observational evidence. The nature of the deep neutron star core, and the state of matter within, remains poorly constrained. However, it is clear that the quality of the gathered data, and our interpretation of it, is improving. As evidence of recent progress, with direct impact on the work discussed in this paper, we may consider the quasiperiodic oscillations observed in the tails of giant magnetar flares (Israel et al. 2005;. The interpretation of the data in terms of torsional oscillations of the neutron star's crust provides us (at least in principle) with a handle on the associated physics. These observations have led to a renewed interest in neutron star seismology, and work addressing issues ranging from the crust equation of state (Steiner & Watts 2009), to the role of the magnetic field (and the coupling to the fluid core) (Glampedakis et al. 2006;van Hoven & Levin 2011;Colaiuda & Kokkotas 2011;Gabler et al. 2011) and even the relevance of the superfluid neutrons that permeate the nuclear lattice at densities beyond neutron drip Andersson et al. 2009). A similar breakthrough is hinted at by the observed cooling rate of the compact remnant in Cassiopeia A. The data provides clear evidence for the presence of superfluid neutrons in the star's core (associated with additional cooling due to the Cooper pair breaking mechanism), leading to a constraint on the superfluid transition temperature (Page et al. 2011;Shternin et al. 2011). This constrains the theoretical models, which depend crucially on many-body interactions at extreme densities. Combined, the information gleaned about the crust region and the superfluid core provides us with a clear motivation to improve our seismology models accordingly. This is not a trivial task, but considerable progress has been made on the key aspects, realistic models would now seem to be within reach.
The work reported in this paper takes a serious step towards the modelling of real neutron star seismology. Our effort should be considered in the context of previous work, such as McDermott et al. (1988) and Strohmayer (1991) who accounted for the crust elasticity and finite temperature effects, the modelling of superfluid dynamics by (for example) Lindblom & Mendell (1994) and most recently Passamonti et al. (2009) and Passamonti & Andersson (2011) in Newtonian gravity and Lin et al. (2008),  and Kantor & Gusakov (2011) in the framework of General Relativity, and finally the effort aimed at including the magnetic field (obviously a key aspect for magnetars!) (Colaiuda & Kokkotas 2011;Gabler et al. 2011). These studies provide key insights into the individual pieces of physics, but do not combine them. It is obviously desirable to do so, and the present work takes important steps in this direction. We consider the oscillations of a neutron star model with an elastic crust permeated by superfluid neutrons and a fluid core with two dynamically distinct components. We develop, and implement, the relevant conditions at the interface between these two regions. We do not, however, account for the magnetic field. Neither do we (fully) account for thermal effects or consider the problem in General Relativity. These are obvious shortcomings of our final model, especially since we cannot meaningfully make use of a truly realistic equation of state unless we work in relativistic gravity. However, our model still represents the state-of-the-art for the coupled multi-fluid-elastic aspects and the analysis provides valuable insights that will assist future developments in this problem area.

THE TWO "FLUID" MODEL
We want to build a neutron star model that accounts for the degrees of freedom associated with the crust elasticity and the superfluid components, both in the crust region and in the core. It is well-known that this involves modelling the core as a two-fluid system, with the superfluid neutron being treated as distinct from a charge-neutral conglomerate of protons and electrons. If the latter components are strictly co-moving, and we ignore issues associated with superconductivity and the presence of magnetic fluxtubes (Glampedakis, Andersson & Samuelsson 2011), then we may ignore magnetic field effects altogether. This is obviously not a true representation of a neutron star core, which should be magnetized, but the model is nevertheless useful as it allows us to explore the relevance of the additional degree of freedom implied by superfluidity.

The dynamics
In the outer core of a neutron star, we expect superfluid neutrons and superconducting protons to be present. Meanwhile, in the inner crust, at densities above the neutron drip threshold, a fraction of neutrons are superfluid and coexist with the lattice of heavy nuclei. As we approach the neutron drip density, these "free" neutrons disappear and all the constituents of the crust are bound in nuclei. A two "fluid" formalism is sufficient to describe these key aspect of neutron star dynamics, but we have to carefully consider the nature of each component in the crust and in the core. In the core it is natural to distinguish the superfluid neutrons (denoted by n in the following) from a conglomerate of charged particles (protons/electrons, denoted by p). The distinction in the inner crust is not quite so obvious, see the discussion by . In this case, we consider the dynamics of baryons confined in the nuclei (denoted by c) permeated by a gas of free superfluid neutrons (represented by f).
With these definitions, the component fractions of the crust and the core are in general discontinuous at the crust/core interface, even if the total number of neutrons and protons is continuous. Therefore, the choice of junction conditions at the crust/core transition is very important for the oscillation dynamics of the star.
The dynamics of a two-component star are governed by mass and momentum conservation equations for each constituent, together with the Poisson equation for the gravitational potential (Prix 2004). The conservation equations take the form: where v i x and w i xy = v i x − v i y are, respectively, the constituent and the relative velocities, whileμx = µx/mx represents the chemical potential (we will assume that the neutron and proton masses are equal, i.e., we take m = mp = mn). The mass densities are defined by ρx = mnx, where nx is the constituent number density. The gravitational potential Φ obeys the Poisson equation, where ρ = x ρx is the total mass density. In the interest of economy we will, whenever we write down an equation that describes both the dynamics in the core and the crust, imply that the two crust components have x either "c" or "f". In the core, these constituent indices should be replaced by "p" and "n" and the shear modulus should obviously be set to zero.
To complete the dynamical equations, we need to provide the quantity εx, that accounts for the entrainment, and the "external" force density f x i that acts on each fluid component. As we will see later, the force can be used to account for the crust elasticity (and, in general, also the magnetic field). Finally, the system of equations is closed by an equation of state (EoS). It can be described by an energy functional that ensures Galilean invariance; The chemical potentialμx and the entrainment parameter εx are then determined from 1

Equilibrium model
Let us first consider the equations that need to be solved in order to determine a non-rotating superfluid neutron star in dynamical and chemical equilibrium. We consider a static, nonrotating, background star in which matter is in β-equilibrium, which means thatμ ≡μp =μn. Moreover, we assume that the crust is unstrained, which means that the elasticity affects only the linear perturbations of the system. Under these assumptions, the star is described by the following equations: Equation (7) is obtained from the stationary Euler equations, the β-equilibrium condition and the standard relation between the pressure and chemical potential ∇P = ρ∇μ. The system of equations (7) and (8) can be written as a single ordinary differential equation The background configuration is therefore completely determined once we choose an EoS that provides a second relation between the mass density and the chemical potential. More details on the EoSs used in this work are given below, in Sec. 3.

Lagrangian perturbations
In order to model the dynamics of the crust, where the elastic restoring force depends on the deviation from the unstrained configuration, it is natural to turn to Lagrangian perturbation theory. In the problem considered in this paper, the relevant perturbations are trivially related to the Eulerian variations of the various quantities, but in the general (rotating and possibly pre-strained) case the Lagrangian description becomes more involved. However, there has been recent progress in developing the required theory , and we draw on these results in the following. In fact, our analysis provides the first application of the new formalism. In the two-fluid situation, it is natural to introduce a Lagrangian displacement vector ξ i x for each fluid component (Andersson, Comer & Gr 2004). In order to describe the motion associated with these vectors, we define the Lagrangian perturbation of a general (scalar or vector) quantity Q along the vector field ξ i where δQ is the Eulerian perturbation of Q and L ξx denotes the Lie derivative along ξ i x . The dynamics of neutron star oscillations can be studied by linearising the system of equations (1)-(3). The linearised mass conservation and Poisson equations are then given by; ∇ 2 δΦ = 4πG (δρx + δρy) .
Meanwhile, for non-rotating stars, the linearised momentum conservation equations take the form ) ( where the elastic stress tensor is defined as The shear modulus is denoted byμ, and must not be confused with the chemical potential (5). For nonrotating configurations in β-equilibrium, equations (13) and (14) assume the simpler form; (1 − εc) ( As noted earlier, Equations (16) and (17) also describe the dynamics of the neutron star core provided we replace, respectively, the indices c and f with p and n, and setμ to zero.

Chemical gauge and crust/core junction conditions
In the inner neutron star crust, the protons and a sizeable fraction of the neutrons are bound by the strong interaction. They form a lattice of nuclei that become heavier towards the crust/core transition, and which supports elastic stresses. The remaining fraction of (unbound) neutrons form a superfluid that may flow through the confined nucleons. This superfluid component first appears at the neutron drip density, ρND = 4.3 × 10 11 g cm −3 , which marks the transition between the inner and outer crust, where there are no free neutrons. In the inner crust there are, at least, two dynamical degrees of freedom; a confined component of protons and neutrons and the free superfluid neutrons. However, the distinction between superfluid and confined neutrons is not obvious. It depends on the local density and the dynamical time-scale of the process under consideration. Adopting the operational definition due to Chamel & Carter (2006), a neutron may be considered free when (on the time-scale considered) its energy is sufficient to overcome the potential barriers separating the nuclei, either classically or by quantum tunnelling. At the bottom of the crust, where the nuclei become denser and form exotic structure (the so-called pasta phase), the potential wells become closer and marginally bound states may exist. These states may penetrate the potential barriers on time-scales that are macroscopically long (but cosmologically short), see Chamel & Carter (2006) for discussion.
Within the two-fluid formalism, this uncertainty may be expressed in terms of the chemical gauge parameter ac introduced by (Carter & Samuelsson 2006;Chamel & Carter 2006;): where n f and nc are, respectively, the number density of the free neutrons and the confined component, while nn and np are the total number densities of neutrons and protons. One can show that the neutron conjugate momentum and the chemical potential are independent of the chemical gauge choice (Carter & Samuelsson 2006;Chamel & Carter 2006; as long as ac is either held fixed or depends only on the nuclear charge number Z. This gauge independence is important for the derivation of the crust/core junction conditions. In order to study the oscillations of a given neutron star model, we must prescribe boundary conditions at the star's origin and surface. In a model with a crust, we must also specify conditions at the crust/core interface and at the outer/inner crust transition. In this work, we implement the static limit of the crust/core conditions derived by . From the mass conservation equations, we obtain the conditions for the radial component of the Lagrangian displacements at r = Rcc: As in the single fluid case the other junction conditions may be derived from the total momentum conservation equation, which leads to the following vertical and horizontal projections : where gij is the (flat) metric, N i is a unit vector orthogonal to the crust/core interface and . . . denotes the change of a physical quantity across the interface. Equations (21) and (22) reduce to the standard traction conditions in the single fluid case. Finally, from the chemical gauge independence we can determine the condition for the chemical potential of the superfluid neutron components; From the junction condition (19), it is clear that for nonrotating models in β-equilibrium equation (23) is equivalent to δμ f = δμn. These conditions complete the description of the problem that needs to be solved in order to determine that star's oscillation spectrum (see Sec. 4.2 for the inner/outer crust boundary conditions).

MODEL EQUATION OF STATE
In order to make maximal use of the theoretical framework, we ought to build our models using a realistic equation of state for supranuclear matter. We will, however, resort to using a simple model EoS. There are several reasons for this. Most importantly, we do not yet have a realistic EoS that provides all the different parameters required in our analysis. Most tabulated EoS simply provide the pressure vs density relation, whereas we need (at the very least) detailed information about the composition, the superfluid pairing gaps and the entrainment between neutrons and protons. Progress in this direction is being made, c.f., Chamel (2008), but the models are not yet complete. Still, the developments have reached the point where it would be meaningful to consider a fully relativistic analysis (as is needed to make the use of a realistic EoS meaningful). This is, however, beyond the scope of the present work. Our main interest here is to probe the key issues and consider the associated phenomenology. For this exercise it makes sense to focus on a simple analytic model with freely adjustable parameters.
We model the neutron star matter in terms of two polytropic EoSs. Despite the simplicity of the construction, combinations of polytropes enable us to study the effects of composition stratification and symmetry energy on the oscillation spectrum. Moreover, we have already used these EoSs in previous work (Passamonti & Andersson 2011;Sidery et al. 2010), which means that we can build on, and compare to, the previous results. As a new application of the polytropic EoS, we construct stratified two-fluid models with different parameters in the crust and the core in order to approximate the qualitatively different compositions in these two regions.
In most astrophysical systems the relative velocity between the two fluids is small. Therefore, equation (4) can be expanded in a series: This approximation has the advantage that the bulk EoS, E0, and the entrainment parameter, α0, can be independently specified at w i fc = 0. From equation (6) it follows that the entrainment parameter εx is related to the function α0 by ρxεx = 2α0 .
In a co-moving background, the density perturbations can be expressed in terms of the chemical potential perturbations as: where Sxy ≡ ∂ρx ∂μy .

The bulk equation of state
Given the above relations, we need to first of all provide the energy E0. The first EoS we consider has also been used by ; Yoshida & Eriguchi (2004); Passamonti et al. (2009): where K is a polytropic constant, xc is the confined nucleon fraction and σsym is a parameter that can be related to the symmetry energy . In this EoS, both xc and σsym are taken to be constant (Passamonti et al. 2009). Using equation (28) and (29), we readily obtain S cf = S fc = σsym Scc .
As discussed by Passamonti et al. (2009), this EoS leads to non-stratified stellar models. The global properties of these models are similar to those of the N = 1 polytrope in the single fluid case. For instance, the pressure is related to the total mass density by the standard equation P = Kρ 2 . The set of models constructed from this EoS will be referred to as models A, c.f. Passamonti et al. (2009). At this point, it is useful to pause and make some comments on the numerical implementation and the resultant stellar models. In our numerical code, we use dimensionless units based on the gravitational constant G, the central mass density ρ0 and the stellar radius R. Unless directly specified, all models A considered in this paper are such that the crust/core transition is at Rcc = 0.9R. Therefore, from the expected crust/core transition density ρcc = 1.2845 × 10 14 g cm −3 (Douchin & Haensel 2001), we can determine the central mass density of the star. At r = Rcc the dimensionless mass density is ρcc = 0.109 ρ0, hence ρ0 = 1.175 × 10 15 g cm −3 . Meanwhile, the dimensionless mass of models A is M/(ρ0R 3 ) = 1.273. Therefore, the physical The right panel shows (on a semi-logarithmic scale) the fraction of the crust/core constituents (vertical axis) and the total mass density (horizontal axis). In the right panel the vertical dotted lines denotes the density of the crust/core transition ρcc and the neutron drip density ρ ND . Note that the total number of neutrons and protons is continuous at r = Rcc, as in the crust a fraction of neutrons is confined in the nuclei of the lattice. Table 1. This table provides the main parameters for the stellar model sequence D. The core is described by the EoS (33) with Nn = 1, Np = 1.4 and kn = 0.766, kp = 2.256. The crust/core transition has been set to Rcc = 0.88R. In the crust, we choose the same polytropic indices for the free neutrons and confined nucleons, respectively, N f = 2.9 and Nc = 1.2, but vary xc at r = Rcc. For different values of xc at the crust/core interface (first column), we show in the second and third column the dimensionless EoS parameters kc and k f , respectively. The dimensionless mass is given in the fourth column. Note that the coefficients kx are given in units of mass and the radius can also be determined. For instance, in the case of the canonical neutron star mass, M = 1.4 M⊙, we obtain the radius R = 12.29 km. It is obviously not the case that the component fractions are constant in a real neutron star. We can account for composition gradients by considering another simple combination of polytropes; where kx and γx are constants. This EoS does not have the symmetry energy term, but it can leads to models with varying composition when γ f = γc. From equations (5) and (33) it follows that the chemical potential and the corresponding mass density are related by where the polytropic index is given by Nx = (γx − 1) −1 . From this result we can determine the proton fraction for a given stellar model by imposing β-equilibrium. After some calculations, we obtain: Forμ → 0, it is clear from equation (35) that xc vanishes for N f > Nc and tends to unity when N f < Nc.
In a "realistic" neutron star, the proton fraction decreases from the centre to the crust/core interface, while the confined component increases from the bottom of the crust and reaches unity at the neutron drip. This behaviour can be approximated using the EoS (33)  In the left panel, we show the neutron effective mass for both the core (neutrons) and the crust (free neutrons). In the crust, the filled circles denote the data set determined by Chamel and collaborators (Carter et al. 2005;Chamel 2005Chamel , 2006. These data have been fitted with the curve (dashed-line) shown in the left panel. For model D with xc(Rcc) = 0.3, we display in the right panel the effective mass for the charged component (core) and the confined particles (crust).
(N f , Nc). For simplicity, we also assume that the equilibrium star is such that the total mass density, the chemical potential and its first derivative are continuous at the crust/core interface. We have found a good compromise that satisfies all these conditions by choosing Nn = 1, Np = 1.4 and kn = 0.766, kp = 2.256, in the core and N f = 2.9, Nc = 1.2 and k f = 0.281, kc = 1.431 in the crust. The crust/core transition has been fixed at Rcc = 0.88R. Note that the coefficients kx are given in units of GR 2 ρ 2−γx 0 . With these parameters, the proton fraction in the core ranges from xp = 0.2 at r = 0 to xp = 0.1 at the crust/core boundary. These values are on the high side, but this is not too concerning here. In the crust, the charged component varies from xc = 0.3 at r = Rcc to xc = 1 at the inner/outer crust interface, r = RND. In this model, the total mass density at the interface is ρcc = 0.115 ρ0 in dimensionless units. As in the case of model A, we can therefore determine the central mass density, ρ0 = 1.117 × 10 15 g cm −3 . For a star with mass M = 1.4 M⊙, we can determine the stellar radius from the dimensionless mass of this model; from M/(ρ0R 3 ) = 1.155 we find that the radius is R = 12.92 km. Fig. 1 illustrates the properties of this crust-core model. The mass density and the chemical potential are shown in the left panel and the constituent fractions are provided in the right panel. The crust/core and neutron drip transitions are also indicated. We have assumed a neutron drip density ρND = 4.3 × 10 11 g cm −3 (Shapiro & Teukolsky 1983;Douchin & Haensel 2001) which for this stellar model corresponds to RND = 0.9988R. For the same core model, we have also studied solutions with different values of xc at the crust/core interface. From equation (35) and the continuity conditions for the total mass density and chemical potential, it is evident that the coefficients kc and k f depend on the parameter xc at r = Rcc. Their values and the total mass of the stellar model are listed in Table 1. We refer to this sequence of stratified stars as models D in order to distinguish them from models B and C considered in previous work (Passamonti et al. 2009;Sidery et al. 2010;Passamonti & Andersson 2011).

Entrainment
In a real neutron star, superfluid neutrons may flow through the lattice of nuclei in the crust and also coexist with superconducting protons in the core. The different constituents are coupled through entrainment, mutual friction, crust/core boundary conditions, EoS, etc. The entrainment is a non-dissipative process that arises from strong interaction. Owing to this effect each superfluid constituent can carry along a fraction of the other component. From a dynamical point of view, it is natural to account for the entrainment by introducing an effective mass associated with each component of the system and which is in general different from the bare particle mass. Meanwhile, in the two-fluid formalism, the entraiment is usually described by the parameter εx defined by equation (6). It is related to the effective mass m ⋆ x as Entrainment affects the constituent conjugate momentum in such a way that it is not longer aligned with the transport velocity of the fluid constituent.
In recent work, Chamel (2008) determined the nucleon effective masses in a neutron star core. The results are based on the two-fluid formalism and describe the strong interaction between nucleons using a two-body force of the Skyrme type with different parametrisations. In a non-relativistic two-fluid model the effective masses can be described by the following relation (Chamel 2008): where ρ b is the baryon mass density, andβ3 is defined bŷ where is the reduced Planck's constant. The quantity B3 is a coefficient that depends on the parametrisation of the Skyrme force (Chamel 2008).
In this work, the superfluid neutron star core is modeled by two fluid polytropes and the effective neutron mass is determined by equation (37). For the proton effective mass, we use the EoS constraint (25), i.e.
together with the definition (36). This model is obviously not consistent, as the effective mass relation (37) has been determined for "realistic" EoS and not for polytropes. However, this is the first time that the complete oscillation spectrum of superfluid neutron star with an elastic crust is considered. Hence, it is natural to focus on the phenomenology sacrificing (to some extent) realism. Our stellar model has the key features expected of a real neutron star, but the "proportions" may not be quite realistic. This should, however, be easy to fix once the problem is considered in General Relativity. The relevant developments are, in fact, already under way. The effective mass of the free superfluid neutrons in the inner crust has been determined within the band theory approximation, but only for a few specific densities (Carter, Chamel & Haensel 2005;Chamel 2005Chamel , 2006. The results indicate that, in the deep crust, the effective mass can be as large as m ⋆ f = 15.4 m f . As already discussed by Andersson et al. (2009), who studied the problem within a local (plane-wave) approximation, these large entrainment values may have considerable effect on the properties of the oscillation spectrum. The currently available entrainment data is shown in the left panel of Fig. 2. We interpolate these data points with a curve, and determine the proton effective mass with equation (39) by applying the method that was used for the core. The results are shown in Fig. 2, where the right panel displays the proton effective mass for a model D with xc = 0.3 at the crust/core interface.

THE EIGENVALUE PROBLEM
For non-rotating stellar models, the study of the linearised equations (11)-(12) and (16)-(17) is considerably simplified if we expand the perturbation variables in vector harmonics. For each set of harmonic indices (l, m) (associated with the familiar spherical harmonics Y m l ), the relevant equations then depend only on the radial coordinate r and time t. In order to set up an eigenvalue problem, we assume that the perturbation variables have a harmonic time dependence e iσt (this essentially amounts to taking a Fourier transform of the equations). The mode eigenfrequency σ is in general a complex number, where the real part describes the oscillation frequency and the imaginary part the growth/damping time of the mode under consideration. However, in this work the eigenfrequency σ is real since we are not accounting for dissipative processes or mechanisms that may trigger instabilities.
The perturbations can be decomposed into two classes of different parity, namely the polar (spheroidal) and the axial (torsional) perturbations. The spectrum of non-rotating stars is also degenerate with respect to the azimuthal harmonic m. This degeneracy is removed when we account for additional physics like rotation or the star's magnetic field. Oscillation modes of a non-rotating star are commonly classified by the dominant restoring force that acts on the displaced fluid elements, the harmonic index l and the radial number n (Cowling 1941) that labels the "overtones" of a given set of oscillations. In simple models the latter parameter can be associated with the number of radial nodes in the relevant eigenfunctions, but this identification no longer holds for more complex systems.
The axial part of the oscillation spectrum is comparatively simple. For non-rotating stars without magnetic fields and crust, the axial perturbations form a degenerate set of zero frequency modes. The crust elasticity breaks this degeneracy as it sustains a set of torsional shear modes. We denote the torsional modes (t-modes) by the symbol l tn, where l is the harmonic index and n labels the overtones. For any l, the frequency of the torsional modes assumes its lowest value for the fundamental torsional mode and increases with n. In compact stars, the typical frequency of the 2 t0-mode is a few tens of Hz (McDermott et al. 1988).
The polar oscillation spectrum is more complex. In addition to the shear waves, the presence of pressure and composition/thermal gradients may generate acoustic and gravity waves. The acoustic modes are mainly restored by pressure variations (p-modes) and cover the high frequency band of the spectrum, above 1 kHz. At lower frequencies, typically below 100 Hz, composition and thermal gradients can sustain the class of gravity modes (g-modes), which are restored by buoyancy. For a given l, there is an infinite series of pressure (gravity) modes with n 1, whose frequency is increasing (decreasing) with n. The ξ r eigenfunction of the pressure, l p n , and gravity modes, l g n , have n radial nodes. The fundamental mode (f-mode), l f, separates the class of p-and g-modes. Its frequency scales with the averaged stellar density and no radial nodes are present in its ξ r eigenfunction.
The crust elasticity affects the polar oscillation spectrum in two ways. Firstly, it generates a class of spheroidal shear modes (s-modes) and secondly it slightly modifies the acoustic mode frequencies (McDermott et al. 1988). Spheroidal shear modes are represented by l sn, with n 1. Their frequencies increase with n, and for small n they lie below the f-mode frequency. In the very low frequency band of the spectrum, there are also interface modes,which are denoted by l i. These modes arise from the presence of internal interfaces and they strongly depend on the physical conditions near these interfaces. For any l, there is an l i mode associated with each internal interface.
The mode classification presented above allows us to describe the complete oscillation spectrum of non-rotating, nonmagnetised single-fluid stars. In superfluid neutron stars, the multi-component dynamics of the system lead to additional degrees of freedom and the oscillation spectrum is therefore richer. In fact, in a two-fluid model we can (usually) decompose the dynamics into a co-moving (counter-moving) motion, where the fluid elements oscillate in (and out of) phase. As a result, we can distinguish two classes of acoustic modes that are commonly called "ordinary" and "superfluid" modes. The ordinary modes have spectral properties similar to the acoustic modes of single fluid stars, as their dynamics is mainly dominated by the co-moving degrees of freedom (and they are therefore restored by the total pressure gradients etcetera). Meanwhile, the superfluid modes exist only in multi-component stars, and the counter-phase motion of the perturbed fluid elements dominates their dynamics. Apart from in particular cases, a neat decomposition of co-and counter-moving motion into "ordinary" and "superfluid" modes is not possible. A typical acoustic mode tends to exhibit both degrees of freedom. Nevertheless, in order to distinguish between these two families of acoustic modes, we denote the ordinary (superfluid) f-and p-modes with the symbols l f o and l p o n ( l f s and l p s n ), respectively. The spectrum of superfluid neutron stars exhibits another important difference with respect to the standard single fluid models. Superfluid neutron stars do not support an independent class of g-modes, even if the model has composition gradients (Lee 1995;Andersson & Comer 2001;Prix & Rieutord 2002). In a perturbed two-fluid system, deviations from chemical equilibrium tend to excite the counter-moving degree of freedom. The buoyancy only introduces small corrections to the superfluid acoustic mode frequencies (see Andersson & Comer 2001, for more details).

The polar perturbation problem
The expansion of the polar perturbations in vector harmonics is given by the following expressions: All scalar perturbations are expanded in terms of the spherical harmonics Y lm , c.f., the perturbed gravitational potential. The problem becomes more tractable if we work with dimensionless quantities. In the core, we choose to work with the following set of perturbation variables: where g is the gravitational acceleration. For spherical stars, this quantity can be related to the chemical potential by the equilibrium equation g = Φ ′ = −μ ′ . The core is described by a system of six perturbation equations that we can write (formally) as where the variables yi are defined in (41), and the function b(r) denotes a generic background field. The complete system of linearised equations is given in Appendix A1.
In the crust, the perturbation equations are evidently more complex due to the presence of the elastic strain tensor. As in single fluid models, it is important to choose a set of perturbation variables that makes the solution of the problem, including the implementation of the crust/core junction conditions, easier. In the inner crust, where a mixture of superfluid free neutrons coexist with a lattice of heavy nuclei, we define the following set of eight dimensionless perturbation variables: where qc has the same units as the shear modulus, and has been introduced to make z2 and z4 dimensionless. In this work, we choose qc = p0, where p0 is the pressure at the origin r = 0. In the function z2, the divergence can be expressed in terms of the harmonic expansion variables since The variable z2 and z4 are directly related to the radial and angular traction components (see Sec 2.4 ). This simplifies the implementation of the junction conditions. Another difference with respect to the functions used in the core (y k ) is that z2 and z6 depend on the Lagrangian perturbation of the chemical potential. This choice provides a simpler boundary-value problem and less complex equations in the core. In the single fluid limit, i.e. when xc = 1, the perturbation variables (43) reduce to the quantities used by McDermott et al. (1988) provided we set qc = p. 2 Naturally, the variables z5 and z6 are obsolete in this limit. The linearised equations of the crust form a system of eight ordinary differential equations that we can write (again, formally) as: We provide the complete form of these equations in Appendix A2.
In the outer crust, all the neutrons are confined to nuclei, i.e. we have x f = 0. Therefore, the perturbation problem reduces to the single-fluid case. The oscillation dynamics is then described by a system of six ordinary differential equations, where the variables z5 and z6 are no longer used.

Boundary Conditions
To solve equations (42) and (45) as an eigenvalue problem, we must impose regularity conditions at the centre, boundary conditions at the surface of the star, and junction conditions at the two internal interfaces. The first internal interface marks the transition between the outer and the inner crust, and the second represents the point where the inner crust gives way to the fluid core.
At the origin (r = 0), the regularity conditions lead to the following relations: where g0 and y 0 k are, respectively, the gravitational acceleration and the perturbation variables at the centre. In Sec. 2.4, we determined the junction conditions at the crust/core surface r = Rcc. In terms of the dimensionless perturbation variables these conditions take the form The boundary conditions that must be satisfied at the inner/outer crust interface are the continuity of the radial displacement z1, the traction components z2 and z4, the gravitational potential perturbation z7 and its radial derivative z8. The component of free superfluid neutrons vanishes on this surface. Given this, we impose a zero boundary condition for the Lagrangian perturbation of the free neutron chemical potential (∆μ f ), i.e. take z6 = 0.
Beyond the transition to the outer crust, the problem is equivalent to the single fluid approximation. Therefore, at the stellar surface we must impose the continuity of the traction perturbation and a condition for the gravitational perturbation. The r → R limit of the linearised Poisson equation leads to the following expression (Unno et al. 1989;Prix & Rieutord 2002): For a model such that the density vanishes at the surface, this means that we should have At the star's surface, the perturbation variables must then satisfy the following relations: where we have imposed a normalization condition for the confined nucleon Lagrangian displacement, z1 = 1. The outer crust of the stellar models studied in this work is very thin. Neutrons start to drip out from the nuclei at r = RND ∼ = 0.9988R. The mode frequency calculation is not really affected if we impose the boundary condition ∆μ f = 0 (z6 = 0) relevant to the inner/outer crust interface at the surface of the star. This is not surprising. The fact that the low density region has little effect on the global oscillations of the star also means that we can ignore the thin ocean that should exist below the density at which the crust is expected to melt (at finite temperatures).

The axial perturbation problem
Axial perturbations of spherical stars describe torsional oscillations which are sustained by the crust's shear. The corresponding vector harmonic expansion of the Lagrangian displacement is given by while all Eulerian perturbations of the scalar variables are zero. From equation (62), it is clear that the vector ξx has zero divergence, ∇iξ i x = 0. Therefore, we have ∆ρx = ∆μx = 0. Substitution of the expansion (62) into the linearised equations (16)-(17) leads to a system of two equations: where we have defined As before, we have assumed a harmonic time dependence with frequency σ, as before. The system of equations (63)-(64) describes the linear dynamics of torsional oscillations in the inner crust. In the outer crust, there are no free superfluid neutrons, which means that we have x f = 0. Therefore, equation (64) becomes obsolete, while equation (63) reduces to the single fluid approximation provided xc = 1 and the entrainment is set to zero, i.e. εx = 0. It is worth noting that, in both the outer and inner crust, equation (63) depends only on the Uc variable, which therefore completely determines the torsional oscillation modes. That this is the case also in General Relativity has already been demonstrated by . Once the eigenvalue problem is solved for the confined nucleons, we can determine the eigenfunctions for the free neutrons in the inner crust from equation (64).
As in the polar case, we introduce a new system of dimensionless variables in order to simplify the boundary value problem. These variables are defined as: where s2 is proportional to the transverse component of the traction. Equation (63) now becomes a system of two first order ODE: The relevant boundary conditions are; the continuity of the traction at the crust/core interface, inner/outer crust interface and at the stellar surface. We therefore impose s2 = 0 both at the crust/core boundary and at the star's surface. In addition, we normalize the radial displacement by setting s1 = 1 at the stellar surface.

RESULTS
Let us now consider the detailed oscillations of the multi-component neutron star model developed in the previous sections. We focus on non-rotating, polytropic two-fluid models, and study the effects of varying component fractions, crust elasticity and entrainment on both the shear modes and acoustic waves. Building on the discussion in Sec. 4, we decompose the problem into polar and axial perturbations. The spectrum of each family of perturbations can be studied as a boundary value problem by solving the system of equations (42) and (45) for the polar modes and (67)-(68) for the axial modes. In both cases we use a shooting method to a fixed point. For the polar sector, we set the fixed point at the crust/core interface, where the core and crust solutions are matched using the junctions conditions discussed in Sec. 2.4 and 4.2.

Shear modes
In simple plane-wave models, the frequency of the crust shear waves is proportional to the shear velocity; However, in the case of a two-fluid model the shear velocity depends on the mass density of the confined nucleons and not on the total mass density, as in the single fluid case ). Still, considering typical parameters, it is easy to see that the shear waves have low frequencies. One would also expect the shear model to be predominantly transverse. They should be present in both the polar and the axial sector, and we will refer to the two classes as spheroidal and torsional shear modes, respectively. The crust elasticity will also affect the acoustic mode frequencies by introducing very small corrections (roughly of the order of the ratio between the shear wave speed and the speed of sound) to the acoustic modes. Furthermore, the free superfluid neutrons in the inner crust may interact with the lattice of nuclei through entrainment and modify the spectrum. The main impact of the entrainment is expected in the deep crust, where the effective mass of the free neutrons may be very large (Chamel 2005), c.f., Figure 2. From a plane-wave analysis of the two-fluid problem, we can derive the frequency dependence of the shear modes on the mass density, the entrainment and the shear modulus ): where k is the magnitude of the wave-vector. We are interested in testing the accuracy of this rough approximation against full numerical results for stratified and non-stratified superfluid models. The shear modulus of the crust can be determined from : . The numerical and estimated frequencies are represented using the same notation as in Fig. 3. The results show that the shear mode eigenfrequencies are well approximated by equation (70).
This result was inferred from Monte Carlo simulations of the Coulomb interactions in the neutron star crust. The quantity ni is the ion density, Z the atomic number of the nuclei and a = (3/4πni) 1/3 is the average inter-ion spacing. The ratio between the Coulomb and thermal energy is parameterised by Γ = (Ze) 2 /akBT , where T is the temperature and kB the Boltzmann constant. Farouki & Hamaguchi (1993) have discussed plausible values for the parameter Γ0.
To be truly consistent with (71), a stellar model should obviously be obtained from a realistic EoS. As we have already discussed, we adopt a different strategy in this work. We focus on simple polytropic two-fluid models and determine the shear modulus from a phenomenological relation. Detailed calculations for realistic EoSs (see Douchin & Haensel (2001)) suggest that the specific shear modulus is almost constant across the crust, i.e. µ ρ ≃ 10 16 cm 2 s −2 .
Hence, it is reasonable to (as a first approximation) consider simple models such thať where µ⋆ is a constant parameter. We consider only this simple model, and in the numerical code, we use its dimensionless formμ * , which is defined as follows: The constant parameter (74) may be estimated at the crust/core interface by using equation (72), the transition mass density ρcc = 1.2845×10 14 g cm −3 at r = Rcc, and assuming typical neutron star values for the central mass density ρ0 = 10 15 g cm −3 and radius R = 10 km. This way we obtainμ * ≃ 10 −4 at r = Rcc.

Spheroidal shear and interface modes
We consider first the spheroidal shear modes of the set of A models; stellar configurations with constant component fractions.
In order to explore the dependence on the shear modulus, we focus on a model with xp = 0.1, xc = 0.3 and zero entrainment. The resulting frequencies of the l = 2 shear modes and the interface mode are shown in the left panel of Fig. 3 for different values of the shear modulus (μ * ). We compare our numerical results to the plane-wave relation (70). To do this we first "infer" the effective value of the wave vector k in equation (70) for a single numerical data point (σ,μ) and then use this value in equation (70) to extend the relation to other parameter values. From the results in Fig. 3, we see that the numerical results scale according to the phenomenological relation. The agreement is, in fact, very good for the shear modes. Meanwhile, the interface mode frequencies deviate from the plane-wave relation. This is expected since these modes have a different origin; they depend sensitively on the physical conditions near the crust/core interface. Interface modes typically have small frequencies and amplitudes concentrated around the internal interfaces. A characteristic property of these modes is a cusp in the radial displacement eigenfunction at the internal boundary layer. In this work, we do not consider neutron stars with  (76). We show also the shear mode frequencies (denoted by stars) for models without free superfluid neutrons in the crust (xc = 1).
an ocean, which generates another interface mode associated with the ocean/crust transition (see McDermott et al. 1988;Piro & Bildsten 2005, for more details).
Let us now apply the same strategy to the dependence on the fraction of confined nucleons xc. In this case, we consider a model A with xp = 0.1, zero entrainment andμ * = 10 −4 . The results shown in the right panel of Fig. 3 confirm the good agreement between the numerical shear-mode frequencies and the analytical relation (70) also in this case.
We can also study the dependence on the entrainment parameter ε⋆. To do this we construct Model A to be a star with fractions xp = 0.1, xc = 0.3, and shear modulusμ * = 10 −4 . We consider two distinct entrainment configurations. In the first case, we set the parameterε constant and determine εx from the relation εx =ε/xy. We have already considered this configuration in a different context (Passamonti et al. 2009;Passamonti & Andersson 2011). In the second case, we investigate the entrainment model described in Sec. 3.2, which is based on the work of Chamel (2005Chamel ( , 2006Chamel ( , 2008. Whenε is taken to be constant, equation (70) describes the frequencies of the shear modes very well, see figure 4. However, the plane-wave approximation is not accurate for the model with a realistic entrainment profile. This result is (obviously) expected, as the plane-wave approximation follows from a local analysis (based on constant parameters). Nevertheless, we find that an "average" of equation (70) over the crust reproduces the shear frequencies of this model rather well. This means that we can use the following shear mode relation: where av denotes an average over the crust. In the left panel of Fig. 5, we compare the shear mode frequencies of model A without entrainment, i.e. withε = 0, to a model based on "realistic" entrainment. The results show that the relation between these two sets of data is accurately described by equation (75), i.e., where σ0 is the frequency for the zero entrainment model. Next we consider the effect of composition gradients on the shear modes, making use of the D models introduced in Sec. 3. With two distinct values for the polytropic indices of the core and the crust, these models mimic the composition of realistic EoS, i.e., the proton fraction decreases from a central value, while the fraction of confined nucleons in the crust increases away from the interface. We study a sequence of models with the same core configuration, such that xp = 0.2 at r = 0, and different constituent fractions in the crust. More details on these models D are given in Sec. 3. Before quantifying the effect that a composition gradient has on the shear modes, we consider the realistic entrainment profile for these stratified models. To this end, we select a model D with xc = 0.3 at the crust/core interface and letμ⋆ = 10 −4 . The obtained shear mode eigenfrequencies are shown in the right panel of Fig. 5. Approximate frequencies determined from the formula (76), are (again) in good agreement with the numerical values.
In order to emphasize the importance of entrainment on the spectrum, we also show in Fig. 5 the shear mode frequencies for stellar models without free superfluid neutrons in the crust (xc = 1). With a large effective mass, the relative motion between the free superfluid neutrons and the crust is restricted. Therefore, the shear mode frequencies tend to the values expected for models without a superfluid component. This effect is clear when we consider the plane-wave approximation (70) in the limit m ⋆ n ≫ m. In this case, we have ε⋆ ≈ xc and equation (70) becomes ) which is the formula expected for pure elastic crust without superfluid neutrons (McDermott et al. 1988). For the stellar models and effective masses considered in this work, there is roughly a 10% difference between the shear mode spectrum determined with and without a superfluid component in the crust. Finally, we study the effects of composition stratification on the spectrum. We consider a sequence of models D, where we setμ⋆ = 10 −4 and take the entrainment to vanish. The original plane-wave relation (70) is (obviously) not adapted to describe the shear modes in stratified stars. However, we obtain a useful description of the spectrum from equation (75), according to which we have; where σ1 is the frequency for the non-stratified, xc av = 1, case. Figure 6 shows the frequencies of the first six shear modes for the set of D models. In the left panel, the numerical eigenfrequencies are compared to the values estimated from equation (78). The accuracy of the averaged plane-wave relation decreases when xc becomes smaller at the crust/core interface. This behaviour is better illustrated in the right panel of Fig. 6, where we show the relative error of the approximation. In particular, we note that equation (75) provides more accurate results for the realistic entrainment configuration effects than for the composition gradients. Therefore, when both the component fraction and the entrainment vary, the shear modes can be described by equation (75) with accuracy limited by the composition gradient effects. Tables 2 and 3 provide the eigenfrequencies of the l = 2 interface mode and the first three l = 2 spheroidal shear modes for some selected models A and D (see the table captions for details). In Table 2, the mode frequencies are given for different confined nucleon fractions xc. In Table 3, the stellar models include the realistic entrainment profile. In order to test these results, we can convert the dimensionless eigenfrequencies of Table 2 into physical units and compare them to results from literature. As far as we know, there are no published results for shear and interface modes in single fluid stars with polytropic EoS and Newtonian gravity. However, there have been a number of analyses of the shear and interface mode spectrum for various other neutron star models. We can asses to what extent our results are consistent with these studies. To this end, we select from the stellar models A and D described in Sec. 3 the model with xc = 1, as the results in the literature rely on the single fluid "approximation". From Table 2 and the central mass densities calculated in Sec. 3, the oscillation frequencies of the 2 s1 and 2 s2 modes are, respectively, ν = 534.68 Hz and ν = 986.49 Hz for model A, and ν = 444.69 Hz and ν = 812.43 Hz Table 2. Frequencies of the l = 2 interface mode and the first three l = 2 spheroidal and toroidal shear modes for some of the stellar models studied in this work. Frequencies are given in units of σ/ √ Gρ 0 , where G is the gravitational constant and ρ 0 represents the central mass density. The first column labels the stellar model. The non-stratified models A have constant proton fraction, xp = 0.1 and the confined nucleon fraction xc shown in the second column. The stratified models D have the same central proton fraction, xp = 0.2 at r = 0, while the values of xc at r = Rcc are given in the second column. All models have zero entrainment andμ⋆ = 10 −4 .
Model xc 2 i 2 s 1 2 s 2 2 s 3 2 t 0 2 t 1 2 t 2  Table 3. This table provides the frequencies of the l = 2 interface mode and the first three spheroidal and toroidal shear modes for models A and D with the "realistic" entrainment determined by Carter et al. (2005);Chamel (2005Chamel ( , 2006. See Sec. 3.2 for more details.

Torsional modes
Motivated by magnetar seismology, the spectrum of toroidal shear modes has already been studied in relativistic two-fluid neutron stars ). Our Newtonian analysis may therefore seem somewhat superfluous. However, if we want to improve our current stellar models by considering the effects of rotation and magnetic field, the perturbation problem becomes more complex (as the dynamics couples the polar and the axial perturbations). As the problem is more tractable in the context of Newtonian gravity, it makes sense to provide the corresponding non-rotating results here. Following the strategy adopted for the spheroidal shear modes, we study the dependence of the l = 2 toroidal modes on the shear modulus, composition fraction and entrainment. As previously, the plane-wave relation (70) describes the spectrum of the non-stratified models A accurately, while the averaged version (75) is more appropriate for both models D and models with realistic entrainment. It is not surprising that the main conclusions drawn for the s-modes remain valid also for the t-modes, as both sets of oscillations are restored by the crust elasticity. Given the similarities to the s-modes, we do not illustrate the torsional mode properties further, but the frequencies of the first three l = 2 t-modes are given in Tables 2  and 3.
The frequency of the fundamental torsional mode, 2 t0, is typically between 10 and 50 Hz (McDermott et al. 1988;Steiner & Watts 2009), depending on the chosen model for the crust. Using the stellar parameters in Table 2, we find that ν = 30.31 Hz for model A and 30.09 Hz for model D. These values depend on the shear modulus asμ 1/2 ⋆ , and hence the frequencies assume lower values forμ⋆ < 10 −4 . The results are certainly consistent with the frequencies expected for the fundamental torsional modes.

Acoustic modes
Superfluid two-constituent stars can sustain two classes of fundamental and pressure modes, the ordinary and superfluid modes, respectively. The ordinary modes, for which the two fluid components tend to move in phase, have spectral properties similar to the f-and p-modes of single fluid stars. Given previous results for that problem, we expect their frequencies to be weakly affected by the presence of the elastic crust (McDermott et al. 1988). The effects of the crust on the superfluid modes, for which the fluid elements oscillate in counter-phase, have not been studied previously. In this section, we provide the first analysis of this problem. We begin by considering the acoustic spectrum of non-stratified stellar models and then extend the analysis to stars with composition gradients.  For a non-stratified model A with standard parameters, i.e. zero entrainment, xp = 0.1, Rcc = 0.9R andμ⋆ = 10 −4 , we determine the l = 2 f-and p-modes for different confined nucleon fractions in the crust, xc. Results for the ordinary and superfluid modes are shown in the left panel of Fig. 7 and in Table 4. The dependence of the 2 f o and 2 p o n modes on xc is very weak, as expected. For instance, the variation of the p-mode frequency for the two models, respectively, with xc = 0.1 and xc = 1.0 is less than 0.01%. Focusing our attention on the ordinary modes, we also study the effects of the crust thickness on the eigenfrequencies. We consider a model A with xp = 0.1 and xc = 0.2, and change the interface location from Rcc = 0.1R to Rcc = 0.99R. The deviation away from the ordinary mode frequency is very small even in the extreme cases, less than 0.05%.
The behaviour of the superfluid modes is different. Figure 7 and Table 4 show that the 2 f s and 2 p s n modes depend significantly on xc, although in different ways. For instance, the frequencies of the 2 f s and 2 p s 3 2 p s 2 modes are increasing (decreasing) functions of xc, while the frequency of the 2 p s 1 mode remains almost constant. The behaviour becomes even more complex when we vary the crust thickness. In the right panel of Fig. 7, we show the superfluid l = 2 f-and p-modes for three models A with different crust thickness, i.e. Rcc/R = 0.85, 0.90 and 0.95, respectively. The crust region clearly affects the spectrum of superfluid modes and it is not (yet) clear to us whether a simple relation can be found to describe the observed behaviour. For the 2 f s mode, we see that thicker crusts and larger fractions of confined nucei both tend to increase the mode frequency.
In order to understand the 2 f s mode better, we study its dependence on the change in composition (due to the chosen chemical gauge) between the core and crust, expressed in terms of ∆x ≡ xc − xp. For a given Rcc, we construct two sets of models A. In the first set, we choose xp = 0.1 and vary xc, while in the second model we choose xc = 0.9 and vary xp. We Model A 2 f s mode x p * = 0.5 Figure 8. This figure displays the dependence of the superfluid l = 2 f-mode frequency on the crust/interface position Rcc and on the composition step ∆x = xc − xp. The dimensionless eigenfrequencies are shown on the vertical axis. The background star is a model A with zero entrainment andμ⋆ = 10 −4 . In the left panel, the 2 f s mode frequencies are determined for two different sequences of models A and for four different crust thickness, Rcc (see legend). The first sequence of stellar models has constant xp ≡ x * p = 0.1 and different xc, while the second sequence has constant xc = 1 − x * p = 0.9 and different xp. The quantity x * p is a parameter that identifies the two sequence of models A. The mode frequencies for the first sequence are denoted with dashed-lines, while those for the second sequence with symbols. The right panel shows the effect of ∆x and Rcc on the 2 f s mode frequencies for two sets of models A with x * p = 0.5. This figure shows that for any two families of models A, denoted by x * p , the eigenfrequencies are described by the same function of ∆, i.e. σ =σx p (∆x).
find that for both these classes of models, the 2 f s mode frequency is described by the same function of ∆x. Results for this function, which we denote σ =σ (∆x), are shown in the left panel of Fig. 8.
We find a similar behaviour if we extend this analysis to other stellar models, e.g., a sequence of models A with constant xp ≡ x * p and a second sequence with constant xc = 1 − x * p . The parameter x * p identifies the two sets of models. For a given x * p , the 2 f s mode frequency has the same functional dependence on ∆x. Note that the frequencies are different for distinct classes of stellar models, i.e. depend also on the parameter x * p , σ =σx p (∆x) .
In the right panel of Fig. 8, we show the 2 f s mode spectrum for two families of models A with x * p = 0.5. In this case, we vary the composition fraction and the crust thickness extending the parameter space towards unrealistic configurations (Rcc = 0.1R) in order to highlight the qualitative behaviour.
Let us now consider the acoustic mode spectrum of stratified stellar models. For models D with zero entrainment, we construct a sequence of stellar configurations with different composition xc at the crust/core interface. The frequencies of the l = 2 f-mode and the first three p-modes for these models are shown in Fig. 9 and Table 4. The main effect of composition stratification is the coupling between the co-and counter-moving degrees of freedom of the oscillations. The interaction between ordinary and superfluid modes is evident in the spectrum (see Fig. 9), where avoided crossings appear as the composition is varied. The chemical coupling also modifies the frequencies of the ordinary modes, as seen from the data in Table 4. For the 2 f o mode, there is a frequency variation of 2.8% between models D with xc = 0.1 and xc = 1.0. However, this change may be due to the properties of the D models being somewhat different. In order to impose the continuity conditions at the crust/core boundary, the total mass density of the crust is a function of xc (see Table 1). This mass variation of the crust can modify the 2 f o mode, as its frequency is proportional to the mean density of the stellar model, σ 2 ∼ M/R 3 . We have not investigated this further, since the determined changes in the mode frequencies is not very significant. Now we turn to the effects that the large effective mass expected in the inner crust has on the acoustic spectrum. We consider a set of models D with different composition xc at the crust/core interface and determine the frequencies of the ordinary and superfluid l = 2 fundamental modes. In the right panel of Fig. 9 we compare the mode frequencies for the vanishing and realistic entrainment (Chamel 2005) cases. The results show that the effective mass has a noticeable impact on the frequency of the fundamental mode. For the 2 f s mode the frequency changes significantly for the entire sequence of models D, while for the 2 f o mode we find a variation of 12.4% for a model D with xc = 0.7, but only 0.01% for xc = 0.3. This is an interesting result that deserves further investigation. However, in order to study more realistic neutron star models we need to work in the framework of General Relativity. Developments in this direction are in progress, building on Lin et al. (2008).  Figure 9. Frequencies of the ordinary and superfluid l = 2 acoustic modes for stellar models with composition stratification. The left panel shows the mode frequency dependence on the crust component xc determined at the crust/core interface, r = Rcc. On the vertical axis the frequencies are in dimensionless units. The background model D has zero entrainment,μ⋆ = 10 −4 and Rcc = 0.88R. The p-mode spectrum exhibits avoided crossings between the superfluid and the ordinary modes. These arise from the interaction between the coand counter-moving degrees of freedom that are coupled in stellar models by composition gradients. The right panel shows the effect of a "realistic" entrainment on the f-mode frequencies. The solid lines represent the results for zero entrainment models (the same as in the left panel), while the dashed lines denote the frequencies for models D with the Chamel entrainment configuration. The large entrainment of the deep regions of the crust affect both the superfluid and the ordinary l = 2 f-mode. We find that the variation of the 2 f o mode frequency is about 12.4% for a model D with xc = 0.7, but decreases to 0.01% for xc = 0.3 Finally, we test the results against the acoustic spectrum obtained from time evolutions of superfluid neutron stars without crust (Passamonti et al. 2009;Passamonti & Andersson 2011). For the non-rotating model A, the frequencies of the ordinary and superfluid acoustic modes agree to better than 0.5%. This comparison is limited by the accuracy of the mode-extraction from the time-evolutions, so the result is acceptable.

CONCLUSIONS
In this paper we have studied the effects of an elastic crust on the oscillation spectrum of superfluid neutron stars. Building on recent results for the perturbations of two-fluid systems with both elasticity and superfluidity , we considered Newtonian models of mature neutron stars. The stellar core was modelled as a mixture of superfluid neutrons and a conglomerate of charged particles, while the inner crust was described by a lattice of nuclei permeated by a gas of free superfluid neutrons. Focussing on nonrotating superfluid stars, we studied for the first time the effects of elasticity, entrainment and composition stratification on both shear modes and acoustic modes. Our results show that the superfluid neutrons in the crust may have considerable impact on the star's spectrum. The superfluid component clearly cannot be neglected if we want to determine an accurate oscillation spectrum of mature neutron stars. If we compare stellar models with and without superfluid neutrons in the crust, the difference between the shear mode frequencies is very large, but becomes smaller (about 10%) when we consider a "realistic" effective neutron mass. A superfluid component in the crust affects also the frequency of the superfluid acoustic modes, which are mainly dominated by a counter-phase motion between the two matter constituents. These superfluid modes are very sensitive to the composition difference between the core and the crust, and to the crust thickness. Ordinary acoustic modes, in which the stellar constituents oscillate mainly in phase, are less affected by the presence of the crust. However, the spectrum of ordinary modes becomes more complex when we consider stellar models with composition stratification and more "realistic" entrainment. The spectrum exhibits avoided crossings between superfluid and ordinary pressure modes due to the chemical interaction between the co-and counter-moving degrees of freedom. Furthermore, for some specific stellar models we have found that a large effective mass may affect also the frequency of the ordinary f-mode at the 10% level.
The results presented in this paper are important as they represent the first analysis of the polar modes of the coupled elastic superfluid system, complementing previous work on the corresponding axial mode problem. Our understanding of the latter problem remains more detailed, as it has also been considered in General Relativity. This is important since one must work in the context of relativity in order to be able to consider realistic neutron star equations of state in a meaningful way. The results we have presented were obtained for "relatively" simple, more phenomenological, stellar models. These models have the advantage that one can explore the dependence on the key parameters, e.g., the composition, but they are obviously not "realistic". Nevertheless, the insights gained from the analysis provide useful guidelines for developing and testing fully relativistic models in the future. The state-of-play as far as such developments are concerned is that the required linear perturbation framework is complete and ready to be used in various applications of astrophysical interest ). This step is, obviously, non-trivial but we expect to make progress in the near future. In order to make these models truly realistic, there also needs to be progress on various equation of state issues. We need more detailed models for the superfluid entrainment and the related effective masses, especially in the crust region. Further work in this direction should certainly be encouraged. After all, the entrainment is one of the key parameters that we may be able to infer from future neutron star seismology.