{'fraction_non_alphanumeric': 0.03221462817504846, 'fraction_numerical': 0.009364480261144548, 'mean_word_length': 4.617235846894339, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 4, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Properties of a quantum vortex in neutron matter

We have studied systematically microscopic properties of a quantum vortex in neutron matter at finite temperatures and densities corresponding to different layers of the inner crust of a neutron star. To this end and in preparation of future simulations of the vortex dynamics, we have carried out fully self-consistent 3D Hartree-Fock-Bogoliubov calculations, using one of the latest nuclear energy-density functionals from the Brussels-Montreal family, which has been developed specifically for applications to neutron superfluidity in neutron-star crusts. By analyzing the flow around the vortex, we have determined the effective radius relevant for the vortex filament model. We have also calculated the specific heat in the presence of the quantum vortex and have shown that it is substantially larger than for a uniform system at low temperatures. The low temperature limit of the specific heat has been identified as being determined by Andreev states inside the vortex core. We have shown that the specific heat in this limit does not scale linearly with temperature. The typical energy scale associated with Andreev states is defined by the minigap, which we have extracted for various neutron-matter densities. Our results suggest that vortices may be spin-polarized in the crust of magnetars. Finally, we have obtained a lower bound for the specific heat of a collection of vortices with given surface density, taking into account both the contributions from the vortex core states and from the hydrodynamic flow.


I. INTRODUCTION
The existence of superfluids in neutron stars was conjectured long before the discovery of those compact objects [1] (see, e.g., Ref. [2] for a recent overview). The first evidence came from accurate radio timing measurements of the rotational frequency of pulsars, revealing sudden spin-ups (whose duration is still not resolved nowadays) followed by relaxations lasting days to months. Such phenomena occurring on a timescale considerably longer than that usually encountered in nuclear-physics experiments can be naturally explained by neutron superfluidity [3]. The frequency 'glitches' themselves are generally interpreted as macroscopic manifestations of the unpinning of neutron quantized vortices in the crust of neutron stars [4,5].
While there is a scientific consensus that superfluidity plays a key role in the glitch phenomenon, many questions still remain open (see, e.g., Ref. [6] for a recent review). This stems from the fact that unlike terrestrial superfluids such as liquid helium and ultracold atomic gases, whose properties can be measured and even experimentally controlled, the cold dense neutron liquid present inside neutron stars cannot be produced on Earth. Although astrophysical observations can indirectly reveal the properties of nuclear superfluids, their interpretation poses a great challenge for theorists. Indeed, it is expected that the global dynamics of neutron stars is de- * daniel.pecak@pw.edu.pl † nicolas.chamel@ulb.be ‡ piotrm@uw.edu § gabriel.wlazlowski@pw.edu.pl termined by the dynamics of about 10 18 vortices (whose typical core size extends over a length of no more than about few tens of femtometers), averaged over a scale of the order of one kilometer for which general-relativistic effects come into play [7,8].
The very large differences in scales suggest a bottomup approach: the motion of large collections of vortices at mesoscopic scales (large compared to intervortex spacing but small compared to the stellar radius) could thus be followed using effective Vortex Filament Models (VFM), which have proved their usefulness for the description of vortices in superfluid 4 He [9][10][11]. However, one has to remember that there are two important differences between vortices in 4 He and in neutron stars. Contrary to bosonic vortices, fermionic ones have their core filled with matter in a normal state. Moreover, fermionic vortices allow for yet another degree of freedom, which plays a crucial role for vortex structure. Namely, the spin imbalance, which in the case of ultracold gases is routinely investigated in laboratories [12], will affect both the internal structure of the core and its size [13,14]. In neutron stars, spin polarization can be potentially induced by strong magnetic fields, but it remains an open question whether it is strong enough to affect the structure of the vortex core. Although the finite size of vortices can be ignored at mesoscopic scales, their quantum structure (as effectively embedded in model parameters) may still have important implications for the stellar dynamics [15]. A better understanding of the microscopic physics underlying effective VFM is therefore highly desirable.
At the smallest scales of interest, which are related to the vortex core size and intervortex spacing, the motion of individual vortices is nonrelativistic (their velocity is small compared to the speed of light [16] and space-time curvature is negligible [17]) so that the methods developed in condensed-matter physics to study laboratory superfluids can be directly adapted to dense stellar environments. The microscopic structure of a single neutron vortex was previously studied at zero temperature solving the Bogoliubov-de Gennes equations [18] (i.e. the single-particle Hamiltonian is that of noninteracting particles) with an effective density-dependent contact pairing interaction fitted to many-body calculations in homogeneous neutron matter using bare nucleon-nucleon potentials. More realistic types of calculations involved the nuclear-energy density functional (EDF) theory within the Hartree-Fock-Bogoliubov (HFB) method [19] using the FaNDF 0 functional [20,21] supplemented with a pairing functional adjusted to 1 S 0 pairing gaps in homogeneous neutron-matter. However, these calculations focused on very dilute neutron matter at densities corresponding to the shallow layers of the inner crust of neutron stars. The pinning of a vortex by a nuclear cluster was later studied in Ref. [22] by solving the HFB equations in a cylindrical cell using the Skyrme SII functional [23] for the normal part, and the phenomenological parametrization of Ref. [24] for the pairing part. More recently, studies of the dynamics of a vortex have been undertaken by solving the fully three-dimensional and symmetry-unrestricted time-dependent HFB equations [25].
As a first step towards realistic simulations of the vortex dynamics in neutron-star crusts, we present in this paper finite-temperature HFB calculations of a single vortex in an otherwise homogeneous neutron superfluid using the Brussels-Montreal EDF BSk31 [26], which was specifically constructed for astrophysical applications. The contributions from time-odd mean fields, which may play an important role in the superfluid dynamics [27,28], are taken into account for the first time.
The paper is organized as follows: in Sec. II we describe the basic ingredients of the approach and the techniques used together with our numerical setup. Then we describe our main results: in Sec. III we describe the length scales associated with the vortex at different temperatures, subsequently the discussion of velocity fields and superfluid fraction is included in Sec. IV. Finally, in Sec. V the specific heat of neutron matter for densities corresponding to the inner crust in the presence of the vortex is discussed. We summarize this article in Sec. VI.

II. HARTREE-FOCK-BOGOLIUBOV EQUATIONS AND NUMERICAL SETUP
A. Nuclear energy-density functional theory The EDF theory has been successfully applied for the description of both the static and dynamic properties of finite nuclei, as well as for nuclear reactions [29][30][31][32].
The building blocks of the semi-local energy-density functional E we consider here, comprise of the fol-lowing local densities and currents (for each nucleon species): particle number density ρ(r r r), kinetic density τ (r r r), anomalous density ν(r r r), and momentum density j j j(r r r), all of which are defined through the Bogoliubov quasiparticle amplitudes u n,σ (r r r), v n,σ (r r r), where n represents the relevant set of quantum numbers while σ = {↑, ↓ } denotes the spin components. The spin-orbit coupling, which plays a minor role in neutron-star crusts (see, e.g., the discussion in Ref. [33]), will be omitted hereafter. Consequently, the quasiparticle amplitudes corresponding to spin-up and spin-down components are related to each other and one needs to calculate u n,↑ (r r r), v n,↓ (r r r) amplitudes only. The densities and currents read (see Refs. [32,34] for details) where the summations are performed over positive quasiparticle energies with a suitable regularization to avoid divergences, as will be discussed in Section II C. In order to allow for the description at finite temperatures, the following thermal occupation factors have been introduced: (setting Boltzmann's constant k B = 1). The quasiparticle amplitudes and quasiparticle energies E n fulfill the HFB equations, which in the coordinate-space representation take the following form for all temperatures: where µ is the chemical potential. The single-particle fields are defined via the variational principle: ∆(r r r) = −2 δE δν * , where {., .} denotes the anticommutator, and δE/δj j j means a vector constructed by variation over three components of the current j j j. The contribution of the mean potential vector field induced by currents to the singleparticle Hamiltonian h(r r r) was ignored in previous calculations. The EDF is the central element of the theory and encodes information about nuclear interactions. It is expressed through the densities and currents (ρ, τ, ν, j j j) and has the generic form: The first term corresponds to the kinetic energy density for a nucleon with mass M , while consecutive terms represent the energy density due to nuclear interactions. The energy densities E ρ and E ∆ρ are related to the interaction of the nucleons with the background density and its fluctuations respectively. The next term E τ is associated with a density-dependent effective mass and also gives rise to current-current couplings (so called entrainment effects) [28]. The last term describes the pairing energy density E π . The contribution related to the spinorbit coupling has been omitted in the presented calculations. This term has been shown to play a minor role in the context of stellar environment where spatial density fluctuations are much smaller than in the case of finite nuclei [33,35,36].

B. Brussels-Montreal nuclear-energy density functionals
The Brussels-Montreal functionals, based on generalized Skyrme effective interactions, were not only precision-fitted to experimental nuclear data (e.g. binding energies, radii), but were specifically constructed for applications to extreme astrophysical environments. In particular, properties of uniform neutron matter, as determined from many-body calculations using bare nucleon-nucleon potentials, were included in their fit. Some of these functionals have been already employed to calculate in a unified and thermodynamically consistent way the internal constitution of neutron stars and their equation of state, from the surface to the core [33,37]. Except for BSk19, most recent Brussels-Montreal functionals have been shown to be compatible with existing astrophysical observations, including the latest constraints inferred from analyses of the gravitational-wave signal GW170817 [38]. For our present purpose, we have adopted the functional BSk31 [26]. This functional was not only accurately fitted to the 2353 measured atomic masses of nuclei with proton and neutron numbers ≥ 8 (taken from the 2012 Atomic Mass Evaluation [39]) with a root-mean-square deviation of about 0.6 MeV, but it was also adjusted to microscopic calculations of the equation of state, effective masses and most importantly 1 S 0 pairing gaps of pure neutron matter. Therefore, this functional appears to be particularly well-suited for the study of neutron superfluidity in neutron-star crusts.
The pairing gaps of Ref. [40], to which the BSk31 functional was fitted, are depicted in Fig. 1 and were obtained from diagrammatic calculations taking into account medium-polarization and self-energy effects. The figure also indicates the densities that will be considered in this paper. These representative densities span different layers of the crust of neutron stars. Focusing on neutron matter, the form of the respective contributions to the total energy in Eq. (9) for the BSk31 functional (see Ref. [26] for details) is the following: where the term κ|∇ρ| 2 (see Eq. (9) in Ref [41]) in the pairing energy density will be neglected hereafter. The coupling coefficients C ρ , C τ , C ∆ρ depend on density ρ, and are defined in Appendix A, see Eqs. (A8)-(A10). The pairing strength is defined according to Ref. [41]: where the effective mass M ⊕ is defined via Eq. (A6). We use the following approximate analytic formula of Ref. [42] for the denominator: where ε Λ is a cutoff that must be introduced to avoid ultraviolet divergences, as will be discussed in the next section. The function Λ(x) reads: The advantage of the present formulation as compared to previous calculations using phenomenological parametrizations such as that of Garrido et al. [24] fitted to a specific value of the cutoff is that the pairing strength adopted here is automatically renormalised for different choices of the cutoff.

C. Convergence and cutoff
The regularization scheme associated with the local pairing functional, defined via Eqs (14)- (16), was originally constructed [43] by cutting off single-particle energies lying above some limiting value ε Λ (as measured with respect to the chemical potential), which was set to 6.5 MeV for the functional BSk31 [26]. Although such a low value allows for fast static HFB calculations, it is not suitable for reliable time-dependent HFB simulations due to the violation of energy conservation [32]. In anticipation of future studies, we shall adopt here values of order ε Λ ∼ 100 MeV. In this case, we thus have ε Λ |∆| therefore to a good approximation we can impose the same cutoff on the quasiparticle energies directly, i.e. by summing over quasiparticle states such that E n < E c ≈ ε Λ .
Let us also recall that the analytical pairing strength, Eqs (14)- (16), was obtained under the assumptions that |∆| µ and ε Λ |∆|. To assess the reliability of these approximations and the validity of our computer code, we have recalculated the pairing gaps ∆ in uniform neutron matter for different cutoff energies ε Λ (approximated by E c ) and we have compared our results with the reference pairing gaps of Ref. [40] shown in Fig 1. We have solved the HFB equations in a cubic box with periodic boundary conditions on a Cartesian grid of 128 3 points with a spacing δx = δy = δz = 1 fm. Results are shown in Fig. 2 for different densities. With increasing cutoff starting from E c = 16.5 MeV, one can see that the error in the pairing gap varies rapidly within 2 − 15% depending on the density but remains fairly independent of the cutoff above E c ≈ 40 MeV. The sudden deterioration of the precision around E c = 180 MeV arises from the discretization of space. Indeed, a finite grid spacing δx prohibits wave numbers higher than k c = π/δx, which translates into an

D. Numerical setup
Our numerical setup to simulate a vortex is similar to that previously adopted in Ref. [25]. Namely, we introduce an axially symmetric external potential to confine the neutron superfluid in a tube. The potential is parametrized as follows (distances are given in fm, and As expected from the weak-coupling approximation |∆| µ underlying the adopted pairing functional, the error is the highest for the lowest density. In the vicinity of the maximum pairing gap, the overall precision of the selfconsistent calculations lies within 5%. energies in MeV): where by r we denote the distance from the z axis. The smooth switching function s in the transition zone is defined as follows: In Ref. [25], the volume of the box was chosen to be V = 60 × 60 × 75 fm 3 , and the parameters of the external potential were set to r 1 = 30 fm, r 2 = 35 fm. To reduce as much as possible the influence of this potential on the superfluid properties, we consider here a cubic box with a larger volume V = (90 fm) 3 . For our setup, r 1 and r 2 are now 40 and 44 fm respectively. The spatial resolution is set to 1.5 fm in each direction, as in Ref. [25]. This corresponds to a maximum accessible momentum k c ≈ 400 MeV/c, or equivalently, a cutoff energy E c ≈ 90 MeV. A vortex is generated at the center of the tube by imposing constraints on the phase of the pairing field ∆(r r r) = |∆(r)|e iφ , where φ = arctan(y/x) is the azimuthal angle. Although we consider here a single vortex (with a winding number equal to 1), our computer code is flexible enough to allow for multiple vortices with higher winding numbers at no additional cost. Due to the symmetry along the z direction, the system is effectively twodimensional. In Fig. 3 we show red translucent contours corresponding to ≈ 0.19MeV of the pairing gap. Due to the presence of the vortex and boundaries, the pairing field is nonuniform and vanishes in the vortex core and at the boundaries of the tube. Therefore, the thin red line in Fig. 3 can be associated with the vortex line, while the red tube roughly corresponds to the external potential. The amplitude of the pairing field |∆(r r r)| is shown at the bottom of the tube. The line at the center shows the vortex core, while the outer tube is connected to boundary effects and roughly delimits the range of the external potential. Value of |∆(r r r)| are plotted at the bottom of the cube (the system is translationally invariant along the z axis). Arrows at mid-height show the vector current j j j(r r r).
In the calculations of the specific heat we changed the geometry described above. Namely, we have considered the volume of V = 180 × 180 × 36 fm 3 , for which we set the external potential parameters accordingly higher to r 1 = 80 fm, r 2 = 87.5 fm. These modifications were introduced in order to minimize the influence of boundary effects and to make sure that the density oscillations induced by the vortex are not influenced by the boundary conditions.

III. LENGTH SCALES
Fermionic systems are characterized by two different characteristic length scales. One is connected to the inverse of the Fermi momentum k −1 F = (3π 2 ρ) −1/3 . The other is specific to superfluid systems: the coherence length ξ = 2 k F /(M π|∆|). In the context of neutron-star crusts, these two length scales are well separated. Even for the largest pairing gap predicted in dilute neutron matter, the coherence length is a few times larger than the average interparticle distance measured by k −1 F . The only system known in Nature, where these two scales become comparable is the unitary Fermi gas realized in ultracold atomic systems [44]. Nevertheless, dilute neutron matter in neutron-star crusts is the nuclear system that shares many similarities with the unitary Fermi gas. Separation of length scales has two important consequences for the vortex structure. First, the spatial length scales over which the pairing field changes will be much larger than those related to the normal density fluctuations. Second, the density of Andreev states inside the vortex core will be significant [14,45,46]. This has important implications for the calculation of the specific heat [47], as will be shown explicitly in Section V.
The top and bottom panels of Fig. 4 display the density profile of a vortex and the associated pairing field, respectively. Characteristic length scales are indicated by vertical lines. In bosonic superfluids, the order parameter is directly related to the density, and as a consequence both of them vanish in the vortex core. This is not the case for fermionic superfluids: although the pairing field, which is related to the order parameter (see, e.g., Eq.(29) of Ref. [28]) vanishes, only a partial depletion of the density is observed in the vortex core, as previously shown for a neutron vortex in Ref. [19]. The existence of a finite density can be traced back to the occupation of Andreev states. Namely, all Andreev states have a nonzero component of angular momentum along the zaxis (vortex line), and consequently the density distribution associated with these states exhibits a minimum in the core. The situation changes when one allows for spin polarization, as could be induced by a strong magnetic field in magnetars [48]. In this case, the majority of the spin particles start to occupy states with reversed angular momenta including the state with 0 . Consequently, the density minimum in the vortex core immediately disappears [14].
As can be seen in Fig. 4, the density variations within the vortex core are relatively small, thus justifying the neglect of the spin-orbit coupling (as in previous studies). The vanishing of the pairing field at the vortex center suggests that the vortex core is a "normal" fluid, as will be explicitly shown in the next section by computing the superfluid fraction. It may be noted that the behavior of the pairing field in the core seems to change within two length scales. Deep inside the core the variation is more rapid, and occurs within a length scale set by k −1 F , whereas at larger distances it varies at distance of the order of the coherence length. The presence of these two length scales in the vicinity of the core has been previously discussed in Ref. [45] in the context of atomic gases. Far enough from the center of the vortex core, the local density and pairing field saturate at values characteristic for a uniform system. This can be clearly seen in Fig. 4: within few coherence lengths from the vortex core, the system tends to be homogeneous. The oscillations seen near the boundaries of the system are due to the imposed external potential. They produce fluctuations on a scale of the order of k −1 F and behave similarly to Friedel oscillations.
While Fig. 4 presents zero-temperature results, we have also performed calculations at finite temperatures. To better characterize the structure of a vortex, we have computed the core radius R core defined as the distance r at which the pairing field increases to 90% of the bulk value far from the vortex. The systematic study of R core as a function of temperature T and density ρ is summarized in Fig. 5(a). Our zero-temperature results are of the same order as those obtained in previous studies [18,49]. In general, thermal effects lead to a reduction of the pairing field thus resulting in a larger vortex core, whose size eventually diverges at the critical temperature T crit = |∆|/1.76 (the vortex disappearing in that case). In the finite system we are considering, this means that R core increases beyond the size of the box.
We have also carried out systematic calculations of the vortex tension per unit length. The tension is defined as the difference of energies between the vortex and uniform matter for a given density and temperature. The results are plotted in Fig. 5(b). The tension of a vortex line is found to increase with density and temperature. The tension grows with density and vanishes at the critical point.

IV. VELOCITY FIELD AND SUPERFLUID FRACTION
At a distance of a few coherence lengths from the vortex, both the density and absolute value of the pairing field becomes essentially uniform. Still, the presence of superflow makes the critical difference between the system with and without the vortex even at large distances from the core. The presence of the superflow can be characterized by the velocity field which decays like 1/r and is a consequence of a particular phase pattern of the pairing field. The presence of the superflow is responsible for long range interaction between vortices and, as can be seen in the next section, produces also a correction to the specific heat of the system. The superfluid velocity field v sf (r r r) is related to the gradient of the phase φ(r r r) of the condensate wave function, which is related to the pairing field ∆(r r r) = |∆(r r r)|e iφ(r r r) . The velocity thus acquires the following form: Note that we have used the mass 2M of the Cooper pair which is twice the mass of the neutron. Another way to describe the flow is to define the velocity through its relation to the mass current [28]: In case of absence of the normal component, it is identical to the superfluid velocity Eq. (19). Both definitions agree at large distances from the center for r ξ, r k −1 F , where the impact of the vortex core is negligible, and v(r r r) ≈ v sf (r r r), but they differ if one approaches the core, see Fig. 6(c). Namely, at a certain distance the velocity of the superflow becomes comparable with Landau's critical velocity and the modulus of the pairing field starts to decrease as the center of the core is approached. This is precisely the transition point where the influence of the core is no longer negligible and affects the motion of neutrons. The identification of the transition point is crucial if one intends to properly define the VFM. The VFM ignores the complexity of the core structure, representing a vortex as a line of vanishing thickness. It assumes that the velocity field decays with the distance from the vortex line as ∼ r −1 . The finite thickness of the vortex line, still sits in the VFM as a parameter, which cannot be set to zero, as this would lead to divergences (the divergence is logarithmic in the so-called local induction approximation [50]). Therefore, we address below the question of determining the suitable vortex core size for the VFM.
In Figs. 6(a) and 6(b), we show the norm of the velocity v v v, and of the current j j j, respectively. Due to the geometry imposed, both quantities lay within a two-dimensional plane. To compare the values of the different length scales, we plot circles centered in the vortex core, having radii of: k −1 F , ξ, R VFM (this is the radius for which the two definitions (19)(20) coincide). Within the radius of ξ the velocity doubles to roughly 2% of the speed of light c, and then within the radius k −1 F the velocity vanishes in the vortex core. Fig. 6(c) presents the cross section of the velocity field through the vortex core. At large distances from the center, the expected behavior v sf (r) ∼ r −1 is reproduced. However, deviations are seen at short distances. In particular, the velocity vanishes in the vortex core, instead of diverging like ∼ r −1 . This reflects the appearance of a normal component, even in the zero temperature limit. We compare both velocities in Fig. 6(c), and mark the area where the two velocities differ by less than 10%. The discrepancy far from the core at x ≈ ±35 fm comes purely from the boundary effects and the finite size of our simulation. However, we define the radius R VFM ≈ 10.5 fm While the former has a finite value in the core, the latter diverges. The green area shows the region where the difference between the two velocities is smaller than 10%. The radius RVFM is defined as the distance from the vortex core to the point where the curves start to match. from the vortex core, at which both definitions of velocity agree. For other densities, see Table I. Our calculations indicate that R VFM is the relevant radius for the VFM, and is typically of order of a few coherence lengths ξ.
From the two definitions of the velocity field, one may introduce the superfluid fraction as the amount of matter that locally takes part in the superflow. Namely, the superfluid fraction η = ρ s /ρ can be estimated by the ratio η = v(r r r)/v sf (r r r) We depict the value of η across the cross-section of the vortex in Fig. 7(a). Combining this with the density distribution, Fig. 7(b), we clearly demonstrate that the vortices in neutron matter (and general in fermionic superfluids) carry some normal component in their core, even in the zero temperature limit. In the next section, we will demonstrate that this component induces modifications to the heat transport in the  7. (a) Superfluid fraction η = v(r r r)/v sf (r r r) for different densities ρ in the zero temperature limit. Far from the vortex core it saturates to 1, while it drops to zero in the core. (b) Density of the system relative to the bulk density ρ∞ far away from the vortex core. Note that the vertical axis does not start with 0, ie. the core is filled with matter. crust.

V. HEAT CAPACITY
The thermal electromagnetic emission from young neutron stars and from transiently accreting neutron stars (see, e.g., Refs. [51,52] for recent reviews) does not directly reveal the properties of their dense core. This stems from the fact that the core remains thermally insulated by the hot crust, whose outermost layers form a heat-blanketing envelope. The (observable) time it takes to reach thermal equilibrium -typically several decadesthus depends on the thermal properties of the crust, and in particular the ratio between the thermal conductivity (mainly due to electrons) and the specific heat. The latter is a measure of the number of degrees of freedom, which in turn is determined by the composition and the structure of the crust. Although the internal constitution of the outer crust is fairly well known, that of the inner region is more uncertain and has been the subject of considerable theoretical efforts [53]. The major part of the inner crust is expected to be made of a Coulomb lattice of spherical neutron-proton clusters immersed in a neutron superfluid and in a charge compensating relativistic elec-tron gas. Due to the delicate interplay between nuclear and Coulomb interactions, the densest layers may consist of more exotic phases referred to as nuclear "pasta".
Since the specific heat coming from single quasiparticle excitations in a uniform superfluid is exponentially suppressed at temperatures much lower than the critical temperature (see, e.g., Ref. [54] for HFB calculations over the whole range of temperatures), the specific heat of superfluid neutron-star crusts is usually thought to be mainly determined by electrons and lattice vibrations (phonons) [51] . The (volumetric) specific heat is thus approximately given by where ε eF is the electron Fermi energy, v is the transverse phonon speed, and the temperature T is assumed to be much lower than the Debye temperature. Longitudinal phonons have much higher velocities and therefore can be neglected. However, pairing is a nonlocal phenomenon so that nuclear clusters may influence the quasiparticle excitations of the neutron superfluid, hence also the specific heat and cooling time [55]. Moreover, the inhomogeneity of neutron superfluid admits the presence of localized in-gap states at the Fermi surface originating from the quasiparticle scattering on the pairing field [56]. The effect of nonlocality has been addressed in the self-consistent HFB calculations in spherical Wigner-Seitz cells [57][58][59] as well as fully three-dimensional bandstructure calculations [60,61]. More importantly, the lattice vibrations are influenced by the neutron superfluid e.g. originating from in-medium renormalization of the nuclear cluster masses [62,63]. The effective low-energy theory describing couplings between superfluid and solid matter has been formulated in Ref. [64]. In particular, the transverse phonon speed is reduced and the corresponding specific heat, the second term in Eq. (21), is increased [65]. Besides, longitudinal lattice vibrations are mixed with the low-energy collective excitations of the superfluid so that their combined contribution to the crustal specific heat may become comparable or even larger than that of electrons [65,66]. The superfluid excitations may be also coupled to nuclear-shape vibrations [67,68]. The latter become energetically favorable in the inner crust due to significant decrease of the nuclear surface tension which is expected to drop by order of magnitude as compared to normal nuclei. At the bottom layers of the inner crust the situation become even more complex. Increased susceptibility towards deformation of nuclear clusters results eventually in the appearance of exotic structures involving very deformed nuclei (pasta phases) [43,69,70]. Moreover due to shell effects [70][71][72] it is expected that the long range order is lost. Consequently it is still not a priori known what are the degrees of freedom of nuclear matter that give the dominant contribution to specific heat. In this paper, we focus on the influence of superfluid vortices on the neutron specific heat. The role of vor-tices is twofold. First, the neutron matter in the vortex core may be excited relatively easily as compared to the surrounding superfluid medium. Second, the superflow around the vortex modifies the quasiparticle excitations and thus introduces corrections to the specific heat of a uniform superfluid. Clearly, whether the vortex contribution is significant or not depends on the density of topological excitations. Our aim, however, is to determine the modification of the specific heat associated with a single vortex in a way that will allow to make predictions for an arbitrary vortex surface density, provided it is much smaller than 1/ξ 2 , where ξ is the coherence length. Before embarking on numerical calculations let us first discuss the source of the increased specific heat associated with a vortex.
Fermionic vortices admit the existence of so called Caroli-de Gennes-Matricon in-gap states [46,47], which give rise to a finite density in the core. In Fig. 8 we show a typical relation between the angular momentum per particle along the vortex axis m (m is the magnetic quantum number) of a given quasiparticle state, and its energy E n . In a superconductor the pairing gap ∆ sets the lowest energy scale. Here, due to the presence of the normal phase in the vortex core, the spectrum is modified and subgap states of lower energy |E n | < |∆| are allowed. The energy of states associated with small enough angular momenta ≤ k F ξ can be quite accurately reproduced within the semiclassical (Andreev) approximation, and follows with a good accuracy the formula [14] E m = 4 3 where ε F = 2 k 2 F /(2M ) is the Fermi energy. For larger angular momenta this band becomes flatter as a function of angular momentum and their energies tend to |∆| from below. Therefore, the contribution to specific heat coming from large angular momentum states becomes smaller and is exponentially suppressed for T |∆|. As a consequence the largest contribution to the specific heat comes from excitations localized inside the vortex core. Note that due to the fact that these states have energies lying within the gap, they are formed by mixtures of particles and holes in almost equal proportions. The number of these states increases with the coherence length as k F ξ reaching large values at the BCS side and vanishing in the vicinity of the unitary point [14,45].
In a uniform superfluid, the existence of a pairing gap ∆ suppresses exponentially the specific heat C V ∝ e −|∆|/T . The presence of quantum vortices introduces states inside the gap. The lowest energy of subgap states so called minigap, which according to Eq. (22) is given by thus becomes an important energy scale, see Fig. 8. It defines the characteristic temperature below which the specific heat will be exponentially suppressed. Moreover, the minigap also determines the minimum strength of the magnetic field that will effectively polarize the vortex core: for this to happen, a quasineutron excitation in the core of energy 2E mg is needed, leading to a critical magnetic field of the order of B crit ≈ 2E mg /(|µ n |), where µ n −1.913µ N is the neutron magnetic moment (µ N being the nuclear magneton). Values of E mg and B crit extracted numerically for densities corresponding to neutron-star crusts are listed in Table I. These values are comparable to those expected to be found in magnetars [73]. Vortex-core polarization may thus occur in these strongly-magnetized neutron stars. Within the temperature range E mg ≤ T ≤ |∆|, the specific heat is expected to have a different form. In the far BCS limit where the interlevel spacing of core states is negligibly small it should follow a linear dependence on T , the same as for a Fermi gas [47]. More precise estimates can be performed based on the density of Andreev states. Namely, let us consider the regime in which the temperature is much smaller than the critical temperature (T T crit ) and at the same time is large compared to the minigap, T E mg , so that excitations of subgap states are allowed. It follows that the heat capacity (per unit of length L) is linear in T and is given by the relation: In practice, however, as will be shown below, the finite interlevel spacing and the increase of the size of the vortex core with temperature lead to departure from this linear behavior.
In order to evaluate the specific heat, we have applied the numerical setup described in Sec. II D. The energy of the vortex configuration as a function of temperature was computed, and next by taking derivative (using finite difference method) the corresponding specific heat C V was extracted. In the case of the tube with the vortex, the contribution to the energy per particle can be expressed as the sum of three terms: The first term corresponds to the uniform fluid with a velocity field generated by the vortex, i.e., ∝ 1/r. The second term is related to the core structure, filled with Caroli-de Gennes-Matricon states. Finally, the last term is due to the boundary imposed in the calculations. It leads to the fluctuation of neutron density close to the boundary. The wavelength of these spatial oscillations is related to the inverse of Fermi momentum 1/k F , like in Friedel oscillations. The geometry of our setup is adjusted to ensure that oscillations due to boundary effects and those due to the presence of the vortex core are spatially separated. Therefore, the first two terms have physical origin, whereas the last one is associated with the numerical setup. The same kind of oscillations arise in the absence of a vortex: the density ρ(r r r) is almost uniform except near the boundary of the cylinder. To get rid of the spurious contribution to the specific heat coming from boundary effects, we will thus consider the following difference: where E uniform is the energy of the same system in the absence of a vortex. The specific heat per unit length of a vortex is shown in Fig. 9(a) for a uniform system (lines), and for a vortex solutions (points) for different densities. The specific heat of the system with a vortex is found to be systematically higher in comparison to the uniform system over the whole range of temperatures. However, the deviations are the most pronounced at very low temperatures T T crit , as can be clearly seen in Fig. 9(b). This expected result is due mostly to the presence of the vortex core inside of which superfluidity disappears, |∆| → 0, as shown in Fig. 4. Clearly at temperatures T approaching the critical temperature T crit , the vortex ceases to exist and therefore the specific heats with and without the vortex practically coincide. For the lowest density we considered, this corresponds to T 0.2 MeV. For the sake of clarity, we have not displayed results for temperatures such that the pairing gap has dropped by 25 % or more compared to its value at T = 0. In the logarithmic scale adopted in Fig. 9(b), the uniform solution varies as the specific heat of a vortex solution is always considerably larger, and exhibits a different type of behaviors dictated by subgap states. To better assess the relative contribution of the vortex on the specific heat, we have plotted in Fig. 10 the ratio between ∆C V /L and the specific heat C uniform V of the corresponding uniform system. The results show clearly that the presence of a vortex may increase the specific heat of the neutron superfluid by several orders of magnitude at all densities provided the temperature is sufficiently low.
The result presented above are obtained for a cylinder of a certain radius, which in practice is much smaller than the typical intervortex distance expected in neutron stars. However, in order to determine the specific heat of a system with vortices, one needs to take also into account the correction coming from large distances from the vortex. This contribution is related to the presence of the superflow and arises from the modification of the quasiparticle spectrum. Here, we present the procedure that allows to infer a lower bound for the specific heat for a given vortex density. Namely, let us consider a large volume S × L, where S is the section area, containing a macroscopic piece of neutron-star crust corresponding to a neutron density n. Let us suppose that the area is threaded by N vor vortices so that their surface density is σ vor = N vor /S. The (volumetric) specific heat difference between the configuration with and without vortices is proportional to σ vor : The first term in the bracket arises from the quantum structure of each vortex core (as discussed above). It is calculated by considering a small cylinder with a typical radius of several tens of femtometers surrounding a single vortex. The second term represents the contribution coming from the vortex flow at large distances r ξ. In the absence of any vortex and ignoring mean-field effects, the quasiparticle energies for a uniformly moving superfluid with velocity v s are given by In order to estimate the modification of the specific heat due to the vortex flow v s = /(2M r), we apply the Thomas-Fermi approximation, namely, we calculate the modification of the specific heat locally, at the distance r from the vortex and integrate the results over the cylindrical shell defined by the radii R in and R out (see Appendix B): where C uniform V is the specific heat of the uniform superfluid characterized by the pairing gap ∆. In practice, R in and R out are of the order of the coherence length and the intervortex spacing respectively, therefore R in R out ∼ S/N vor /π. Note that Eq. (27) provides only a lower bound since the interactions between vortices, which we have neglected by treating each vortex independently, bring additional contributions to the specific heat.
From the equation (29), it is clear that the correction of the specific heat due to the vortex flow is positive. Assuming T |∆|, this contribution is significant whenever or equivalently, whenever the temperature is of the order of T v as given by where we have approximated µ by the Fermi energy ε F . However, it should be remarked that the contribution to the specific heat induced by the vortex flow is proportional to C uniform V , therefore it is exponentially suppressed in comparison to the specific heat contribution associated with the vortex core.
The procedure described above allows to estimate a lower bound for the specific heat associated with a given vortex density in neutron-star crusts, using the microscopic results presented in this paper and the formula (27). The relative importance of the two terms in Eq. (27) depends on the actual vortex density. Given a typical surface vortex density of a neutron star σ vor ≈ 10 −21 fm −2 , and the temperature of the crust T /∆ = 0.1, Eq. (29) yields that the flow related heat capacity ∆C flow V is 16 orders of magnitude smaller than C uniform V . In view of the results shown in Fig. 10, ∆C flow V is therefore negligible compared to ∆C V . This stems from the very low surface density of vortices. However, the clustering of vortices due to the pinning might increase the contribution of vortex flow to the specific heat. Let us remark that this contribution may be substantial in physical systems studied in terrestrial laboratories, such as superfluid helium and ultracold fermionic condensates and for which the vortex core size is comparable with the intervortex spacing.

VI. CONCLUSIONS
In summary, we have carried out systematic fully selfconsistent 3D HFB calculations of a quantum vortex in neutron matter for different temperatures and for various densities corresponding to different layers of the inner crust of a neutron star. We have adopted the BSk31 EDF from the latest series of the Brussels-Montreal family, which was specifically designed for such applications. We have determined the effective radius relevant for the VFM. Moreover, we have extracted the vortex core size as a function of temperature and have shown that it diverges at the critical temperature. We have shown that the superfluid fraction drops to zero inside the vortex core, thus indicating the presence of a normal phase. On the contrary, the density remains finite due to the presence of so-called Caroli-de Gennes-Matricon in-gap states. We have determined their spectrum numerically and we have found that the magnetic field prevailing in magnetars may spin polarize the vortex core.
We have also shown that the neutrons located in the vortex core are responsible for an additional contribution to the heat capacity. This effect considerably modifies the heat capacity at low temperatures, i.e., temperatures comparable to the energy of minigap E mg . At temperatures above the minigap but below the critical temperature the specific heat does not increase linearly with T , as expected [47]. We attribute this effect to the dependence of the core size on temperature and to the relatively low density of subgap states as compared to the BCS limit. Besides the contribution coming from vortex core state, we have obtained a lower bound for the specific heat associated with the hydrodynamic flow induced by vortices. This contribution turns out to be negligibly small in neutron stars unless the surface density of vortices is sufficiently large. On the other hand, the associated contribution may be substantial in terrestrial superfluids, where the density of vortices is large, even comparable with 1/ξ 2 . Moreover, an additional increase of this term would stem from the fact that the flow farther than R out from the vortex cannot be neglected.
With the relatively large quasiparticle energy cutoff we have adopted here (thus ensuring energy conservation on long time scales), the wavefunctions we have computed could be used as initial data to study the time evolution of neutron vortices in neutron stars, in particular vortex collisions and dynamical instabilities. I. Properties of neutron matter for different bulk densities ρ∞. The first three rows give characteristic length scales: the inverse Fermi momentum k −1 F , the coherence length ξ, and the effective radius for the VFM RVFM. The following four rows give the energy scales: the bulk pairing gap ∆∞ and corresponding critical temperature Tcrit, the Fermi energy εF, the chemical potential µ, and the minigap Emg. The last row give the critical magnetic field Bcrit needed to spin polarize the vortex core.
ρ∞ the other terms and therefore will be neglected. Clearly C V (v s ) = C V (v s = 0) + ∆C V , where the first term corresponds to the specific heat of the static superfluid, whereas ∆C V is the correction due to the superflow, and since T /|∆| 1 it is positive.
To apply this formula to the flow induced by a vortex, we apply the local density approximation. We consider a single vortex inside a cylinder of radius R out . In order to evaluate the correction to the specific heat coming from the region between radii R in and R out we use the fact that the magnitude of superfluid velocity behaves like v s (r) = /(2M r), where r is the distance from the vortex axis. The radius R in denotes the distance from the core, where the fluctuations of density and pairing field become negligible. In this paper we take R in as the radius of the cylinder where the microscopic calculations have been performed. Subsequently, we divide the volume between R in and R out into infinitesimal concentric cylindrical shells (of length L) in which the superfluid velocity is well defined. Integrating the contribution dC V (v s (r)) of each shell from r = R in to r = R out yields Equivalently, one may express this correction through the specific heat C uniform V of the uniform static superfluid, characterized by the density of states N (0) and the pairing gap ∆: With the above expression, we can estimate the minimal contribution to the specific heat coming from the hydrodynamic flow outside the vortex core for any given surface density of vortices σ vor = N vor /S. Assuming that the vortices are independent and the flow contributes only up to the distance R out , we get a lower bound for the heat capacity by multiplying Eq. (B9) by the number N vor of vortices and substituting R out = 1/ √ πσ vor :