Collective excitations in the neutron star inner crust
aa r X i v : . [ nu c l - t h ] O c t Collective excitations in the neutron star inner crust
Luc Di Gallo, ∗ Micaela Oertel, † and Michael Urban ‡ LUTH, Observatoire de Paris, CNRS, Universit´e Paris Diderot, 5 place Jules Janssen, 92195 Meudon, France Institut de Physique Nucl´eaire, CNRS/IN2P3 and Universit´e Paris Sud 11, 91406 Orsay, France
We study the spectrum of collective excitations in the inhomogeneous phases in the neutron starinner crust within a superfluid hydrodynamics approach. Our aim is to describe the whole rangeof wavelengths, from the long-wavelength limit which can be described by macroscopic approachesand which is crucial for the low-energy part of the spectrum, to wavelengths of the order of thedimensions of the Wigner-Seitz cells, corresponding to the modes usually described in microscopiccalculations. As an application, we will discuss the contribution of these collective modes to thespecific heat of the “lasagna” phase in comparison with other known contributions.
PACS numbers: 26.60.Gj
I. INTRODUCTION
Neutron stars are fascinating objects, containing mat-ter under extreme conditions of temperature, density andmagnetic field. In order to study these celestial bodies,theoretical modeling has to be confronted with observa-tions. A prominent observable is the thermal evolution ofisolated neutron stars. Properties of the crust thereby in-fluence the cooling process mainly during the first 50-100years, when the crust stays hotter than the core whichcools down very efficiently via neutrino emission (see e.g.[1]). Heat transport in the crust is the key ingredientto explain the afterburst relaxation in X-ray transients,too [2, 3]. Concerning the models for the thermal relax-ation of the crust, the most important microscopic ingre-dients are thermal conductivity and heat capacity, andto less extent neutrino emissivities. Here, as a first ap-plication, we will concentrate on the heat capacity. Moredetails about the evaluation of the heat capacity and adiscussion of the usually considered contributions can befound in [4]. In what follows, we will concentrate onthe particularly interesting case of the neutron star innercrust.The core of neutron stars is composed most probablyof homogeneous neutron rich matter, whereas the crustcontains different inhomogeneous structures. The innercrust is thereby characterized by the transition from a lat-tice of atomic nuclei in the outer crust to homogeneousmatter in the core. Ravenhall et al. [5] and Hashimoto etal. [6] predicted that this transition passes via more andmore deformed nuclei. Starting from an almost spheri-cal shape, they could form rods or slabs immersed in aneutron gas at the different densities. These “spaghetti”and “lasagna” phases are commonly called the nuclear“pasta”. At higher densities, even closer to the core,other phases such as neutron-gas bubbles inside the densematter (“swiss cheese” phase) etc. are expected. The ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] formation of the different structures strongly dependson the relative strength of the nuclear surface energy,Coulomb energy and bulk energy, such that it dependson the nuclear interaction. This prediction has been con-firmed within different models for the nuclear interaction,see, e.g., [7–10]. These evaluations have been performedat zero temperature. It is clear that at some criticaltemperature the pasta structures will disappear due tothermal excitations. However, this melting temperatureis of the order of several MeV (see e.g. [11, 12]).Another point is that in neutron stars older than sev-eral minutes, matter becomes superfluid. A first evidencefor superfluidity in neutron stars has been discussed al-ready in 1969 [13], shortly after the discovery of the firstpulsars, in connection with the observation of “glitches”.Since then much effort has been devoted to the ques-tion of superfluidity and superconductivity in neutronstar matter, for the inner crust as well as for the homo-geneous core, see for example [8]. There is no consen-sus on the exact value of the energy gaps ∆ in the innercrust [14, 15], but the common agreement is that they areof the order of 1 MeV [16]. A pairing gap much largerthan the temperature strongly suppresses the contribu-tion of individual neutrons to the specific heat which isthus very much dependent on the pairing strength [4, 17].For moderate and strong pairing, the main contributionsto the heat capacity considered so far in the crust arethus electrons and lattice vibrations as well as collectiveexcitations of nuclei. However, the superfluid characterof neutron star matter induces collective excitations ofthe neutron gas, not considered before, which can givean important contribution to the heat capacity in cer-tain regions, see [18–20].The aim of the present paper is to study these col-lective excitations in the inner crust employing a super-fluid hydrodynamics approach. Naturally, there is a vastliterature on hydrodynamics for neutron stars in differ-ent contexts, including the effects of superfluidity [21–27]. Most of these models are dedicated to the study ofmacroscopic neutron star properties, whereas our mainaim is to study the excitation spectrum of the crust onmuch smaller length scales. In spirit this is similar toRefs. [19, 20, 28], where hydrodynamic equations are de-veloped to study the superfluid Goldstone boson and (lat-tice) phonons in the long wavelength limit. However, weare interested in shorter wavelengths at which effects ofthe inhomogeneous structure will manifest themselves,too. In this sense, our approach is situated in betweenthe long wavelength limit and the completely microscopiccalculations (see e.g. [17, 29, 30]) employing the Wigner-Seitz approximation [31]. In the latter case, the wave-lengths are limited to the size of the Wigner-Seitz cell,because of the imposed boundary conditions which donot include the coupling between neighboring cells.The paper is organized as follows. In Sec. II, we de-scribe the superfluid hydrodynamics approach. We sum-marize the hydrodynamic equations, discuss the bound-ary conditions and the microscopic input. In Sec. III, weshow our first results. For simplicity, we restrict ourselvesto one-dimensional inhomogeneities (lasagna phase) inthis exploratory study. We discuss the spectrum of thecollective modes and their contributions to the specificheat. A summary and perspectives of our work are ex-posed in Section IV.Throughout the article, c , ¯ h , and k B denote the speedof light, the reduced Planck constant, and the Boltzmannconstant, respectively. II. MODELA. Superfluid hydrodynamics approach
In this paper, we are interested in temperatures be-low ∼ K, which is very small compared to the gapenergy ∆. Therefore we can use the zero temperatureapproximation, thus assuming that there are no normalfluids but only superfluids. In this limit, the dynamics ofa superfluid system with slow temporal and spatial vari-ations is completely determined by the dynamics of thephase of the order parameter: If the superfluid order pa-rameter is written as ∆( r ) = | ∆( r ) | e iφ ( r ) , the superfluidvelocity is given by v s = (¯ h/ m ) ∇ φ [32], m being thenucleon mass. Actually, as pointed out in Ref. [33], thephase of the order parameter determines the momentumper particle p = m v s and not the fluid velocity v . Thisdistinction is important in the context of “entrainment”in a system containing protons and neutrons, see below.An important length scale is the superfluid coherencelength, ξ = ¯ hv F /π ∆ [34], where v F denotes the Fermivelocity. It varies from several fm up to tens of fm for typ-ical values of the densities, neutron fractions, and gaps inthe inner crust. As can be shown by deriving the equa-tions of superfluid hydrodynamics from the microscopictime-dependent Bogoliubov-de Gennes (or Hartree-Fock-Bogoliubov) equations [35, 36], the hydrodynamic ap-proach is valid if the length scale of spatial variationsis larger than ξ , and frequencies are small comparedwith ∆ / ¯ h . As will be discussed later on, for the concreteexamples we consider, we are at the limits of validity ofthe approach. However, considering the tremendous dif- ficulties to perform completely microscopic calculationsbeyond the Wigner-Seitz approximation, we leave suchinvestigations for the future and consider our approachsufficient for the moment.In addition, we will neglect the Coulomb interaction ofthe protons. This represents an enormous simplification,but at the same time it implies that we cannot correctlyreproduce the phonons of the Coulomb lattice. In homo-geneous matter, too, the Coulomb interaction plays animportant role for the collective modes, in particular thecoupling of the proton plasmon mode with the electronsas discussed in [37, 38], but it is beyond the scope of thepresent paper. Our main focus lies therefore on the dy-namics of the neutron gas, which is however coupled tothe proton dynamics due to the nuclear interaction.In principle, the equations of superfluid hydrodynam-ics can be derived from the underlying microscopic the-ory, as it was done for the case of ultracold trappedfermionic atoms in [35, 36]. Here, we follow the simplerway to derive them from local conservation laws. Sincethe fluid velocities and the densities are low enough, wewill use a non-relativistic formulation.The first conservation law is neutron and proton num-ber conservation [51]. This results in two continuityequations, one for neutrons ( a = n ) and one for protons( a = p ), ∂ t n a + ∇ · ( n a v a ) = 0 , (1)where n a denotes the particle number density of species a . The second conservation law, the conservation of mo-mentum, results in the Euler equations, which can bewritten as n a ( ∂ t p a + ∇ ˜ µ a ) = 0 , (2)where ˜ µ a is the rest-frame chemical potential defined asthe conjugate momentum with respect to the particledensity n a in a variational approach [39]. The explicitexpression is ˜ µ a = µ a + v a · p a − m a v a , (3)where µ a is the local chemical potential of species a . Dueto the interaction between neutrons and protons, µ a de-pends on the densities of both species.In pure neutron matter, the momentum p n is simplygiven by p n = m n v n . However, in a system containingneutrons and protons, the two species drag each otherdue to their interaction. In the theory of superfluids, thiseffect is known as entrainment [40]. As a consequence,fluid momenta are misaligned with particle velocities.The relationship between the velocity and the momen-tum can be expressed via the entrainment matrix (alsocalled Andreev-Bashkin or mass-density matrix) [33]: m a n a v a = X b = n,p ρ ab p b m b . (4)In practice, at densities which are relevant in the innercrust, the non-diagonal elements of ρ are small [33], i.e., ρ ab ≈ m a n a δ ab . (5) B. Microscopic input
As microscopic input, we need the equation of state,i.e., the relation between the densities n a and the chem-ical potentials µ a , and the entrainment matrix ρ . In ourconcrete numerical examples, we will use the results ofthe work by Avancini et al. [10] for the equilibrium config-urations. They evaluate the structure of the pasta phasesfor charge neutral matter in β equilibrium using a den-sity dependent relativistic mean-field model, the DDH δ model (originally called DDH ρδ ) [10, 41, 42], for the nu-clear interaction. In order to be consistent, we shall cal-culate the chemical potentials µ a and the entrainmentmatrix ρ with the same interaction. For the entrainmentmatrix, we closely follow Gusakov et al. [43], who gen-eralized the determination of the entrainment matrix forneutron-proton mixtures based on Landau-Fermi liquidtheory [44] to relativistic models. The only modifica-tion of the expressions in Ref. [43] we have to performis due to the presence of the isovector-scalar δ mesonin the DDH δ model, which modifies the Dirac effectivenucleon mass. In particular, the latter is no longer thesame for neutrons and protons. Since our hydrodynamicequations are formulated non-relativistically, we consideronly the non-relativistic limit of the entrainment matrix( ρ ab = m a m b c Y ab in the notation of [43]). C. Linearization around stationary equilibrium
In order to proceed we will linearize Eqs. (1) and (2)around stationary equilibrium. Let us write the differentquantities as a sum of their equilibrium value and a per-turbation, X = X eq + δX (in the case of the velocitiesand momenta we will write the perturbation simply as v a and p a , respectively, since the equilibrium values of thesequantities are zero). The equations can be simplified alot, since all temporal and spatial derivatives of equi-librium quantities vanish (except at phase boundaries,which will be treated in the next subsection). Eqs. (1)and (4) then reduce to ∂ t δn a = − X b = n,p ρ ab, eq m a m b ∇ · p b , (6)and Eqs. (2) and (3) become ∂ t p a = − ∇ δµ a . (7)We will now express the variation of the densities inEq. (6) in terms of the variation of the chemical poten-tials, δn a = X b = n,p J ab δµ b , (8) where J ab = (cid:16) ∂n a ∂µ b (cid:17) eq . (9)Inserting the resulting equation into the divergence ofEq. (7) one obtains the following system of two coupledwave equations for δµ n and δµ p : X b = n,p ( KJ ) ab ∂ t δµ b = ∇ δµ a , (10)where K is the inverse of the matrix( K − ) ab = ρ ab, eq m a m b . (11)The coupling arises from the non-diagonal elements ofthe matrices J and K due to the neutron-proton interac-tion. Let us now make the ansatz that the perturbationshave the form of a plane wave, δµ a ( r , t ) = U a e − iωt + i k · r .Eq. (10) can then be written as a 2 × X b = n,p ( KJ ) ab U b = 1 u U a , (12)with u = ω/k denoting the sound velocity. The twoeigenvalues give two sound velocities which we will label u ± . Note that the corresponding eigenvectors, U ± a , donot describe pure proton or neutron waves, but combina-tions of both. We denote by + and − the modes whereneutrons and protons oscillate in phase and out of phase,respectively.In the special case of pure neutron matter, there isonly one mode, which can be obtained from the aboveequations by setting n p = 0. Its sound velocity is givenby u = (cid:16) n n m n ∂µ n ∂n n (cid:17) eq . (13) D. Boundary conditions
In our model, we consider the inhomogeneous phasesin the inner crust as mixed phases where a neutron gas(phase 1) coexists with a dense phase (phase 2) contain-ing neutrons and protons. However, in order not to haveto write everything separately for phase 1 and phase 2,we will write all equations, unless otherwise stated, forthe general case that neutrons and protons are present inboth phases. The equations relevant for phase 1 can eas-ily be obtained by considering the special case n p = 0.The fact that both phases coexist implies that in equi-librium the chemical potentials and pressures are equalin both phases: µ a = µ a and P = P . The descrip-tion of the interface between the two phases requires amicroscopic formalism and is beyond the scope of thiswork.In our model, we assume that the hydrodynamic equa-tions are valid in both the gas and the dense phase, butsince they do not say anything about the behavior atthe interface, they have to be supplemented by appropri-ate boundary conditions. The first boundary conditionarises from the obvious requirement that contact has tobe maintained at all times at the interface [45]. There-fore, the displacement normal to the interface has to becontinuous and equal for all components ( a = n, p ) at alltimes. Hence, the velocities normal to the interface mustsatisfy: v ⊥ n ( r ) = v ⊥ p ( r ) = v ⊥ n ( r ) = v ⊥ p ( r ) . (14)The second boundary condition arises from the require-ment that the pressure P on both sides of the interfacemust be equal [46]: P ( r ) = P ( r ) . (15)If we linearize this condition, it can be written as X a = n,p n a ( r ) δµ a ( r ) = X a = n,p n a ( r ) δµ a ( r ) , (16)where the index eq after n a and n a has been droppedfor brevity.Before applying our model to the neutron-star innercrust, let us see whether these boundary conditions givereasonable results for collective modes in isolated nuclei.For simplicity, we will consider a nucleus with equal num-bers of neutrons and protons ( N = Z = A/ r = R . The proton andneutron densities inside the nucleus are n n = n p = n / n = 0 .
153 fm − is the saturation density withinthe DDH δ model. As a first example, we consider theisoscalar monopole mode, where neutrons and protonsoscillate together in radial direction. The solution ofthe wave equation inside the nucleus is δµ n = δµ p ∝ j ( ωr/u + ), where j l is a spherical Bessel function and u + = 0 . c is the sound velocity for the in-phase oscil-lation of neutrons and protons. Since protons and neu-trons move together, the first boundary condition (14) isautomatically satisfied, while the second one, Eq. (16),requires δµ ( r = R ) = 0. Consequently, the energy of themonopole mode is ¯ hω = π ¯ hu + /R ≈
90 MeV /A / .Another interesting simple case is the isovector giant-dipole resonance (GDR), where neutrons and protons os-cillate against each other in z direction. In this case,our approach is identical to the Steinwedel-Jensen modelof the GDR [47]. Again, the solution of the waveequation is straight-forward and gives δµ n = − δµ p ∝ j ( ωr/u − ) cos θ , where θ is the angle between r andthe z axis, and u − = 0 . c is the sound velocity forthe out-of-phase oscillation of neutrons and protons. Inthis case, the second condition (16) is automatically ful-filled, but now the first one, Eq. (14), becomes rele-vant. Using the Euler equation (7), one can show thatthe radial component of the velocity field is proportionalto v rn = − v rp ∝ ∂δµ/∂r ∝ j ′ ( ωr/u − ) cos θ , so thatEq. (14) gives ¯ hω = 2 . hu − /R ≈
82 MeV /A / . FIG. 1: Diagram representing the slab structure.
These results for the isoscalar monopole and the isovec-tor GDR are quite reasonable, at least for heavy nuclei,although their energies are much higher than the pair-ing gap ∆, such that superfluid hydrodynamics shouldstrictly speaking not be applicable. The reason is thatfor these particular resonances (contrary to, e.g., thequadrupole mode [48]) the Fermi-surface distortion doesnot play any role, so that hydrodynamics works even inthe normal phase without pairing. We conclude that, atleast in some cases, the limits of validity of the hydrody-namic approach may be interpreted very generously.
E. Collective modes in a periodic slab structure
Because of their electric charge, the droplets (or rods,or slabs) of the dense phase arrange in a regular peri-odic lattice in order to minimize the Coulomb energy.Charge neutrality on a macroscopic scale is guaranteedby the presence of an almost uniform, strongly degener-ate electron gas. The size and form of the structures isdetermined by the interplay of Coulomb energy (favor-ing small structures) and surface energy (favoring largestructures). Since both the Coulomb and the surface en-ergy are neglected in our approach, the determination ofthe size and form of the structures in equilibrium is be-yond the scope of our work. Instead, we will considerthe equilibrium geometry as input and calculate the col-lective oscillations in this geometry. For simplicity, wewill restrict ourselves to the simplest geometry which isa structure of periodically alternating slabs with differ-ent proton and neutron densities as illustrated in Fig. 1(lasagna phase). To be specific, we will consider the slabsto be perpendicular to the z axis.Our aim is to describe the collective modes of thisstructure. The equilibrium properties of the structureitself, i.e., the densities n n , n p , n n and n p , and theslab thicknesses L and L , are input parameters whichwe take from Ref. [10]. The excitations are then obtainedby solving in each slab the wave equation (10) togetherwith the boundary conditions (14) and (16) at the inter-faces between neighboring slabs.At each phase boundary, the waves will be partially (ortotally) reflected. It is therefore not sufficient to make aplane-wave ansatz in each slab, but one has to considerthe reflected wave, too. Thus, we can make the followingansatz in each slab: δµ a ( r , t ) = X σ = ± e − iωt + i k k · r k (cid:16) α σ e ik σz z + β σ e − ik σz z (cid:17) U σa , (17)where U ± a denote the normalized eigenvectors of Eq. (12), k k = ( k x , k y ,
0) and r k = ( x, y,
0) are the components of k and r parallel to the slab, and the k ± z have to satisfy k ± z = ω u ± − k k . (18)Note that k ± z can be real or imaginary. The velocities v a can be expressed in terms of the coefficients α ± and β ± ,too. By using the Euler equation (7), one finds that for aplane wave with wave vector k the velocity field is givenby v a = k n a ω X b ( K − ) ab δµ b . (19)If we define V ± a = 1 n a X b ( K − ) ab U ± b , (20)the superposition of plane waves according to Eq. (17)gives v za = X σ = ± k σz ω e − iωt + i k k · r k (cid:16) α σ e ik σz z − β σ e − ik σz z (cid:17) V σa (21)and a similar relation for v k .The next step is to determine the coefficients α ± and β ± by matching the solutions in neighboring slabs ac-cording to the boundary conditions. If we use indices1, 2, 3 in order to indicate the quantities in three con-secutive slabs, we have perturbations δµ a valid for 0 Let us now investigate the resulting excitation spec-trum for a specific example. As mentioned before, thevalues for the equilibrium quantities will be taken fromthe work by Avancini et al. [10], who have studied thestructure of pasta phases in a relativistic mean fieldmodel. Our geometry corresponds to the lasagna phase,appearing close to the transition to uniform matter inthe core, which has been found in Ref. [10] in the case ofzero temperature and β -equilibrium for baryon numberdensities 0 . 077 fm − < ∼ n B < ∼ . 084 fm − , in good agree-ment with the results by Oyamatsu [7]. For our examplewe have chosen an intermediate density, n B = 0 . 08 fm − .The corresponding properties of the two phases 1 and 2are listed in Table I.With the actual numbers for the densities and the di-mensions of the structure, the coherence length for agap of 1 MeV is of the same order of magnitude asthe size of the layers, i.e. the scale for spatial varia-tions. That means that our superfluid hydrodynamicsapproach touches its limit of validity for this example.Strictly speaking, we should also limit ourselves to ener-gies which are small compared to ∆. However, there aremany cases where the hydrodynamic approach works rea-sonably well although its initial assumptions are not ful-filled. Examples are the dipole and monopole resonancesin ordinary nuclei mentioned in the preceding section,or the “supergiant resonances” in spherical Wigner-Seitzcells used to model the neutron-star inner crust, whose excitation energies agree well with an estimate obtainedfrom the sound velocity of the hydrodynamic Bogoliubov-Anderson mode [30].After this remark of caution, let us discuss the solu-tions for the energies ω shown in Fig. 2 as functions of q ≡ | q | for three different angles θ between q and the z axis (i.e., q z = q cos θ and q k = q sin θ ). The left panelshows the dispersion relation for waves propagating in z -direction, i.e. perpendicular to the interfaces betweenthe different slabs. One observes an acoustic branch withan approximately linear dispersion law ω = u s q (27)at low energies, and several optical branches with a finiteenergy for q = 0, analogously to phonons branches in acrystal.Note that within the Wigner-Seitz approximation,which is usually employed in microscopic calculations[17, 29, 30], we would only obtain a discrete spectrumcorresponding to our spectrum in the case q = 0. Thereason is that in this approximation the coupling betweencells is neglected, and thus each cell has the same excita-tion spectrum. The degeneracy of the modes in each cellis lifted by the coupling between cells, which gives riseto a momentum dependent spectrum as obtained in ourapproach.The slope of the acoustic branch, i.e., the speed ofsound, coincides (see dashed line in Fig. 2) with the usualthermodynamic expression for the sound velocity u s = 1 m ∂P∂n B (cid:12)(cid:12)(cid:12) Y p , (28)where n B is the average baryon density of the inhomo-geneous phase. To evaluate this derivative, we squeezeor expand our unit cell of length L by a small amount δL = δL + δL . From the requirement δP = δP = δP we can determine δL and δL and thus δP . The finalresult can be written in a compact form as Ln B u s = L n B u s + L n B u s , (29)where we have defined for each phase i = 1 , u si = 1 m ∂P i ∂n Bi (cid:12)(cid:12)(cid:12) Y pi . (30)Note that u s is identical to the sound velocity u [cf.Eq. (13)], whereas u s is different from the two soundvelocities u ± .This linear branch corresponds roughly to the longwavelength limit discussed in Ref. [19, 20] although, ofcourse, the numerical value of the sound speed is not thesame because we neglect elastic effects of the proton lat-tice due to Coulomb interaction. At higher wave vectors q , there are deviations from the linear behavior relatedto the inhomogeneous structure. At these energies thelong wavelength limit is no longer valid. - h ω [ M e V ] |q| [fm -1 ] θ = 0 0 1 2 3 4 5 6 0 0.05 0.1 0.15 0.2 0.25|q| [fm -1 ] θ = π /4 0 1 2 3 4 5 6 0 0.05 0.1 0.15 0.2 0.25|q| [fm -1 ] θ = π /2 FIG. 2: Dispersion relation of the modes propagating along the z -axis ( θ = 0, left), at an angle of 45 ◦ ( θ = π/ 4, center), and inthe plane parallel to the slabs ( θ = π/ 2, right). The dashed line in the left panel corresponds to the approximation Eq. (27). The central and right panels of Fig. 2 show the excita-tion spectrum for different values of the angle, θ = π/ π/ 2, respectively. One observes that the slope of theacoustic branch discussed before changes: in the presentexample, u s increases from 0 . c in the case θ = 0 to0 . c in the case θ = π/ 2. The reason is that the wave,which is perfectly longitudinal ( v k q ) in the case θ = 0,becomes more complicated in the case θ = 0 and thenucleons oscillate now in both longitudinal and trans-verse directions. But the most important consequence ofnon-zero angle θ is the appearance of a second acousticbranch, whose slope is strongly angle dependent. In fact,if one writes the energy of this new branch as ω = u ′ s q k = u ′ s q sin θ , (31)the “two-dimensional sound velocity” u ′ s defined by thisequation depends only weakly on q z and q k : in thepresent example, u ′ s varies between 0 . c for q z ≪ q k and 0 . c for q k ≪ q z . A detailed analysis of the solu-tions for the coefficients α and β corresponding to thisbranch shows that in this mode, the protons and neu-trons oscillate practically only in the direction parallel tothe slabs (i.e., v z ≈ B. Application to specific heat We are interested here in the contribution of the abovediscussed excitation modes to the specific heat. The spe-cific heat, the heat capacity for constant volume per unit volume, is defined as c v ( T ) = ∂ǫ∂T (cid:12)(cid:12)(cid:12) n , (32)where ǫ denotes the energy density. The contribution ofthe collective modes can be calculated as follows: ǫ ( T ) = Z π/L − π/L dq z π Z d q k (2 π ) X i ¯ hω i ( q ) 1 e ¯ hω i ( q ) /k B T − . (33)Note that we suppose here that the energies ω i ( q ) de-pend only very weakly on temperature such that it isjustified to neglect their temperature dependence. Thisshould be a good approximation as long as the tempera-ture stays well below the value of the energy gap and wetherefore have no significant contribution from a normalfluid. Another type of temperature dependence couldarise from a change in the structure of the pasta phases.At the temperatures considered here, however, we donot expect a significant effect either since the structurestarts to be modified considerably only at higher temper-atures [11, 12].In Fig. 3 we show the different contributions to thespecific heat in the density range where the model byAvancini et al. [10] predicts the lasagna phase, for a typ-ical temperature of 10 K. Besides the contribution ofthe collective modes (solid line), we display for compari-son the contribution of the electrons (dashed line), whichare considered as a practically uniform ultra-relativistic( µ e ≫ m e c ) ideal Fermi gas with number density n e = n p . At low temperature, the electron gas is strongly de- c v [ e r g / ( c m K )] n B [fm -3 ]Collective excitationsElectronsweakly paired neutrons FIG. 3: Different contributions to the specific heat in thedensity range where one expects to find the lasagna phase,for T = 10 K. Solid: collective modes, dashed: electrons,dotted: neutron quasiparticles (from Ref. [17]). Concerningthe conversion between astrophysical and nuclear units, notethat 10 K = 86 . k − B keV and 10 erg K − cm − = 7 . × − k B fm − . generate and its contribution to the specific heat reads c el .v = k B µ e T (¯ hc ) . (34)The importance of the collective modes becomes clearif one considers the contribution of the gapped neutronquasiparticles (dotted curve), taken from Ref. [17]: Inthe absence of collective modes, an excitation of the neu-tron gas requires the breaking of Cooper pairs, which issuppressed by a factor of the order of e − ∆ /k B T . Evenin the case of weak pairing, at the present temperature,this contribution is suppressed by approximately one or-der of magnitude with respect to the contribution of thecollective modes.In Fig. 4, we show the temperature dependence of thespecific heat corresponding to the intermediate-densitycase discussed in Sec. III A (solid line). For comparison,we again display the specific heat due to the electrons(dashed line). Due to its linear temperature dependence,Eq. (34), the electron contribution is always dominantat low temperature, but at higher temperature, the con-tribution of the collective modes is comparable or evenlarger than the electron contribution.At the low temperatures considered here, which arewell below the energy of the first optical branch, the con-tribution of the collective modes to the specific heat iscompletely dominated by the two linear branches dis-cussed in the preceding subsection. As is well known[50], the specific heat due to an acoustic branch with alinear dispersion relation, Eq. (27), reads c v = 2 π k B T h u s ≡ bT . (35)In the present case, however, we have seen that there isin addition a “two-dimensional” branch which propagates c v [ e r g / ( c m K )] T [10 K]Collective excitationsElectronsApproximation a T + b T FIG. 4: Temperature dependence of the contribution of thecollective modes (solid line) and of the electrons (dashes) tothe specific heat for the the example studied in the precedingsubsection (see Table I). The approximate formula , aT + bT , see Eqs. (35) and (36), is shown as a dotted line. only parallel to the slabs and whose dispersion relationis approximately given by Eq. (31). The contribution ofsuch a mode to the specific heat is readily shown to be c v = 3 ζ (3) k B T π ¯ h u ′ s L ≡ aT , (36)where ζ is the Riemann zeta function [ ζ (3) = 1 . . . . ].Due to its quadratic temperature dependence, this is thenext dominant contribution at low temperatures after theelectrons. The result of the simple formula aT + bT ,where a and b have been calculated with the averagevalues u s = 0 . c and u ′ s = 0 . c , is shown as thedotted line in Fig. 4. Up to the temperatures consideredhere, it fits reasonably well the full calculation. IV. SUMMARY In this paper, we have presented a formalism of super-fluid hydrodynamics to treat density-wave propagation ininhomogeneous pasta-like nuclear structures which ap-pear in the inner crust of neutron stars. To accountfor the periodicity of the structure, we incorporate theFloquet-Bloch boundary conditions. The idea is some-where in between the approaches of Refs. [19, 20], con-sidering only the long-wavelength limit, averaging overthe microscopic details of the structure, and microscopiccalculations of the crust within the Wigner-Seitz approx-imation [29, 30], valid for wavelengths smaller than theradius of the Wigner-Seitz cell. Concerning the micro-scopic input for the nuclear equation of state and the ge-ometry of the structure, we followed the work by Avanciniet al. [10].Within this approach, we have calculated the excita-tion spectrum of a periodic structure of parallel slabs,the lasagna phase. We have shown that the structurecan indeed induce non-negligible effects on the excita-tion spectrum. In particular, we found that the soundvelocity of the usual acoustic mode depends on the di-rection of the propagation, and, more surprisingly, thatthere is a second acoustic mode whose dispersion relationis almost independent of q z . In addition, we found dif-ferent optical branches, similar to the phonon spectrumof ordinary crystals.We have calculated the specific heat corresponding tothis excitation spectrum and found that its contributionis much more important than that of individual neutrons,which is strongly suppressed due to the superfluid gap∆. At temperatures relevant for neutron stars, the maincontributions to the specific heat come from the electronsand from the acoustic collective modes. The latter cannotbe obtained within the Wigner-Seitz approximation. Dueto the curious sound mode whose energy is independentof q z , the specific heat due to the collective modes goeslike T instead of T . [With the same arguments, onewould predict that in a rod structure (spaghetti phase),the specific heat should be linear in T .] However, it isnot clear whether this feature survives when the Coulombinteraction, which has been neglected here, will be takeninto account.Of course, in order to treat the complex geom-etry, we were obliged to make a couple of ap-proximations. Contrary to the microscopic ap-proaches based on the Quasiparticle-Random-Phase-Approximation (QRPA) [29, 30], we rely on the assump-tion that the modes can be described hydrodynamically,which implies in particular that the local neutron andproton Fermi surfaces stay spherical at all times. Thisassumption is justified if all spatial variations are slowcompared to the superfluid coherence length and the tem-poral variations are slow compared to the superfluid gap. Both assumptions are not very well fulfilled. However,we have cited examples where hydrodynamics gives rea-sonable answers even beyond these very restrictive lim-its, and we believe that the results are at least qualita-tively correct. The most serious limitation of the presentwork is probably that the Coulomb interaction has beenneglected. The Coulomb interaction between the pro-tons results in an additional coupling between neighbor-ing cells, which can have important consequences for theexcitation spectrum. In the approaches of Refs. [19, 20],it was accounted for by including the elasticity of theCoulomb lattice. In our more microscopic approach, theCoulomb potential would have to be included from thebeginning into the proton chemical potential µ p ( r , t ) inthe Euler equation (2). This is a difficult task which willbe left for future studies.It has to be stressed that the contribution of the col-lective modes studied here is potentially more importantthan other contributions, notably the contribution fromindividual neutrons. Therefore it is interesting to pursuetheir investigation and to include the additional contribu-tion to the specific heat in studies of neutron star thermalevolution. Acknowledgments We are indebted to C. Da Providencia for providingus with the data for densities and geometries of the dif-ferent pasta phases in the model of Ref. [10] as well asto M. Fortin for providing us with the data of Ref. [17].We thank S. Chiacchiera, D. Pe˜na, and M. Fortin fordiscussions. This work was supported by ANR (projectNEXEN), and by CompStar, a research networking pro-gramme of the European Science Foundation. [1] D. G. Yakovlev, O. Y. Gnedin, A. D. Kaminker, andA. Y. Potekhin, AIP Conf. Proc. , 379 (2008).[2] P. S. Shternin, D. G. Yakovlev, P. Haensel, and A. Y.Potekhin, Mon. Not. R. Astron. Soc. , L43 (2007).[3] E. F. Brown and A. Cumming, Astrophys. J. , 1020(2009).[4] O. Y. Gnedin, D. G. Yakovlev, and A. Y. Potekhin,Mon. Not. R. Astron. Soc. , 725 (2001), arXiv:astro-ph/0012306.[5] D. G. Ravenhall, C. J. Pethick, and J. R. Wilson, Phys.Rev. Lett. , 2066 (1983).[6] M. Hashimoto, H. Seki, and M. Yamada, Prog. Theor.Phys. , 320 (1984).[7] K. Oyamatsu, Nucl. Phys. A , 431 (1993).[8] C. J. Pethick and D. G. Ravenhall, Ann. Rev. Nucl. Part.Sci. , 429 (1995).[9] G. Watanabe, K. Sato, K. Yasuoka, and T. Ebisuzaki,Phys. Rev. C , 035806 (2003).[10] S. S. Avancini, L. Brito, J. R. Marinelli, D. P. Menezes,M. M. W. de Moraes, C. Providˆencia, and A. M. Santos,Phys. Rev. C , 035804 (2009). [11] G. Watanabe, K. Iida, and K. Sato, Nucl. Phys. A ,455 (2000), arXiv:astro-ph/0001273.[12] S. S. Avancini, S. Chiacchiera, D. P. Menezes, andC. Providencia, Phys. Rev. C , 055807 (2010).[13] G. Baym, C. Pethick, and D. Pines, Nature , 673(1969).[14] U. Lombardo and H. J. Schulze, Lect. Notes Phys. ,30 (2001), arXiv:astro-ph/0012209.[15] A. Gezerlis and J. Carlson, Phys. Rev. C , 025803(2010), arXiv:0911.3907.[16] N. Chamel and P. Haensel, Living Reviews in Relativity , 10 (2008), arXiv:0812.3955.[17] M. Fortin, F. Grill, J. Margueron, D. Page, and N. San-dulescu, Phys. Rev. C , 065804 (2010).[18] D. N. Aguilera, V. Cirigliano, J. A. Pons, S. Reddy, andR. Sharma, Phys. Rev. Lett. , 091101 (2009).[19] C. J. Pethick, N. Chamel, and S. Reddy, Prog. Theor.Phys. Suppl. , 9 (2010), arXiv:1009.2303.[20] V. Cirigliano, S. Reddy, and R. Sharma (2011),arXiv:1102.5379.[21] R. I. Epstein, Astrophys. J. , 880 (1988). [22] B. Carter and D. Langlois, Nucl. Phys. B , 478(1998), arXiv:gr-qc/9806024.[23] N. Andersson and G. L. Comer, Mon. Not. R. Astron.Soc. , 1129 (2001), arXiv:astro-ph/0101193.[24] B. Carter, N. Chamel, and P. Haensel, Nucl. Phys. A , 441 (2005), arXiv:astro-ph/0406228.[25] M. E. Gusakov and P. Haensel, Nucl. Phys. A , 333(2005), arXiv:astro-ph/0508104.[26] B. Carter and E. Chachoua, Int. J. Mod. Phys. D15 ,1329 (2006), arXiv:astro-ph/0601658.[27] B. Carter and L. Samuelsson, Class. Quant. Grav. ,5367 (2006), arXiv:gr-qc/0605024.[28] A. D. Sedrakian, Astrophysics and Space Science ,267 (1996).[29] N. Sandulescu, Phys. Rev. C , 025801 (2004).[30] E. Khan, N. Sandulescu, and N. V. Giai, Phys. Rev. C , 042801 (2005).[31] J. W. Negele and D. Vautherin, Nucl. Phys. A207 , 298(1973).[32] E. M. Lifshitz and L. P. Pitaewskii,