Glitch Rises as a Test for Rapid Superfluid Coupling in Neutron Stars

Neutron stars provide the unique opportunity to study matter under extreme conditions. Learning about their unknown nuclear equation of state (EoS) relies on understanding the connection between the macroscopic observables and microphysics. One possibility is to probe the interior physics with glitches. These sudden spin-ups interrupt the regular pulsar spin-down (Espinoza et al. 2011) and are typically associated with the transfer of angular momentum from a crustal superfluid, decoupled from the lattice (and everything tightly coupled to it) due to vortex pinning (Anderson & Itoh 1975). Upon reaching a critical lag, the glitch is triggered and a large number of vortices simultaneously unpin. The frictional forces acting on free vortices on small scales and mechanisms causing their gradual repinning subsequently govern the macroscopic post-glitch response (Pines et al. 1980). The latter are typically associated with an exponential recovery and modeled within vortex-creep theory (Alpar et al. 1984b, 1993; Akbal et al. 2017), whereas the former dominate the behavior at early times. In this paper, we focus on the glitch rise.

Observations of the Vela pulsar suggest that crust coupling is very efficient: initial constraints for the spin-up timescale (Dodson et al. 2002, 2007) have been recently improved showing that the crust accelerates within $\sim 5\,{\rm{s}}$ (Palfreyman et al. 2018) after the glitch is initiated. Within hydrodynamical models, this rapid recoupling is captured via a dimensionless mutual friction coefficient ${ \mathcal B }$ (directly related to the vortex dynamics) as the timescale to recouple the bulk superfluid is ∝1/2Ω_sf${ \mathcal B }$ (Alpar & Sauls 1988; Andersson et al. 2006). Provided that neutron stars are continuously monitored, spin-ups can be "caught in the act," allowing access to the early transient dynamics and the corresponding mutual friction coefficients, which in turn are controlled by the underlying small-scale processes. The most promising candidate to study this connection between macro- and microphysics is the Vela pulsar and dedicated observation campaigns have been performed, for example, at the Mount Pleasant Radio Observatory, Tasmania and the Hartebeesthoek Radio Astronomy Observatory, South Africa. Using the former, Palfreyman et al. (2018) have recently reported the first single-pulse observations of a sudden spin-up, providing the most detailed information of the glitch rise to date.

As a first step toward making realistic predictions for such an observation and constraining neutron star microphysics, we review two existing calculations (Epstein & Baym 1992; Jones 1992), as well as highlighting inconsistencies between them, analyzing the mechanism that is held responsible for rapidly recoupling the crustal superfluid: the excitation of Kelvin waves along superfluid vortices. Instead of following previous work using constant mutual friction coefficients (e.g., Haskell et al. 2012), we subsequently determine the Kelvin wave coupling strength for a realistic crust model, discussing uncertainties in the microscopic parameters. These new density-dependent couplings are further combined with a simplified treatment of the core coupling and implemented in a three-component neutron star model to make a prediction of the initial glitch response of a Vela-like pulsar. This is followed by a comparison between our predicted glitch rise and single-pulse observations of the 2016 December Vela pulsar glitch (Palfreyman et al. 2018).

Following the large-scale unpinning initiating a pulsar glitch, vortices move with a local velocity Δv relative to the crustal lattice. Provided that Δv is sufficiently large (see below for more details), the excitation of circularly polarized Kelvin waves dominates the dissipation. On small scales, these dynamics are fully characterized by a dimensionless drag parameter ${ \mathcal R }$, because an individual vortex feels a resistive force per unit length,

$$\begin{eqnarray}&&{f}_{\mathrm{res}}={\rho }_{{\rm{s}}}\kappa { \mathcal R }{\rm{\Delta }}v,\end{eqnarray} \tag{ 1 }$$

where ρ_s ≡ m_un_s is the mass density of the free crustal superfluid, m_u the atomic mass unit, n_s the superfluid number density, and κ ≈ 2.0 × 10⁻³ cm² s⁻¹ the quantum of circulation. Assuming that a large number of vortices moves freely and experiences f_res, the microscopic drag is related to the large-scale hydrodynamic mutual friction coefficient by (Glampedakis et al. 2011)

$$\begin{eqnarray}&&{ \mathcal B }\equiv \displaystyle \frac{{ \mathcal R }}{1+{{ \mathcal R }}^{2}}.\end{eqnarray} \tag{ 2 }$$

To obtain ${ \mathcal B }$, the drag coefficient ${ \mathcal R }$ has to be known. Kelvin wave dynamics have been addressed by Epstein & Baym (1992) and Jones (1992), albeit arriving at different results for the corresponding dissipation. To provide context for these papers and to discuss the origin of the discrepancy, we use a simplified version of Epstein & Baym's argument to derive the expected scalings for ${ \mathcal R }$. The equation of motion for forced vortex oscillations reads

$$\begin{eqnarray}&&{\rho }_{{\rm{s}}}\kappa \hat{z}\times \displaystyle \frac{\partial {\boldsymbol{\epsilon }}}{\partial t}+T\,\displaystyle \frac{{\partial }^{2}{\boldsymbol{\epsilon }}}{\partial {z}^{2}}={\boldsymbol{f}},\end{eqnarray} \tag{ 3 }$$

where ${\boldsymbol{\epsilon }}$ is the displacement of a vortex aligned with the z-direction, T the vortex tension, and ${\boldsymbol{f}}$ the driving force per unit length. In the absence of forces, a plane wave ansatz shows that the vortex supports Kelvin waves with characteristic frequency (Thomson 1880; Jones 1990)

$$\begin{eqnarray}&&{\omega }_{k}=\displaystyle \frac{{{Tk}}^{2}}{{\rho }_{{\rm{s}}}\kappa }=\displaystyle \frac{{\hslash }{k}^{2}}{2\mu (k)}.\end{eqnarray} \tag{ 4 }$$

Here, k is the wave number along a vortex, ℏ is the reduced Planck constant, and μ(k) is an effective mass that varies slowly with k. This dispersion relation provides the tension associated with a specific mode, T = ρ_sκℏ/2μ.

Consider a point interaction with a lattice nucleus in which a force f ∼ δ(z)E_p/ℓ is exerted on the vortex over a time, τ ∼ ℓ/Δv. E_p and ℓ are the pinning energy and typical interaction scale. This will excite Kelvin waves of characteristic frequencies ω ≲ τ⁻¹ and wave numbers $k\lesssim {k}_{* }\equiv {(2\mu /{\hslash }\tau )}^{1/2}$, related to the Fourier-transformed amplitude $\tilde{\epsilon }(k\lesssim {k}_{* })\sim {E}_{{\rm{p}}}\tau /{\rho }_{{\rm{s}}}\kappa {\ell }$ (see also Link 2003). The energy associated with the perturbations is

$$\begin{eqnarray}&&{\rm{\Delta }}E\sim \int {{Tk}}^{2}\tilde{\epsilon }{(k)}^{2}\,{dk}\sim \displaystyle \frac{{\hslash }{k}_{* }^{3}}{{\rho }_{{\rm{s}}}\kappa \mu }{\left(\displaystyle \frac{{E}_{{\rm{p}}}}{{\rm{\Delta }}v}\right)}^{2}.\end{eqnarray} \tag{ 5 }$$

As we are concerned with the scalings, numerical prefactors have been dropped. They are reintroduced below.

According to Epstein & Baym (1992), the power transferred into Kelvin waves per unit length is

$$\begin{eqnarray}&&p\sim {n}_{{\rm{l}}}{\rm{\Delta }}v\int {\rm{\Delta }}E(b){db}\sim \displaystyle \frac{{\hslash }{k}_{* }^{3}}{{\rho }_{{\rm{s}}}\kappa \mu }{\left(\displaystyle \frac{{E}_{{\rm{p}}}}{{\rm{\Delta }}v}\right)}^{2}{n}_{{\rm{l}}}{\ell }{\rm{\Delta }}v,\end{eqnarray} \tag{ 6 }$$

where we have n_l nuclei per unit volume and the integral over impact parameters b is cut off at the scale ℓ. A vortex, hence, experiences the resistive force f_res = p/Δv per unit length. With Equation (1) and k_*, we obtain

$$\begin{eqnarray}&&{ \mathcal R }\sim {\left(\displaystyle \frac{\mu }{{\hslash }}\right)}^{1/2}{\left(\displaystyle \frac{{E}_{{\rm{p}}}}{{\rho }_{{\rm{s}}}\kappa }\right)}^{2}\displaystyle \frac{1}{{\rm{\Delta }}{v}^{3/2}}\displaystyle \frac{{n}_{{\rm{l}}}}{{{\ell }}^{1/2}}.\end{eqnarray} \tag{ 7 }$$

Epstein & Baym (1992) consider a vortex-nucleus interaction potential,

$$\begin{eqnarray}&&E(s)=\displaystyle \frac{{E}_{{\rm{s}}}}{{(1+{s}^{2}/{R}_{{\rm{N}}}^{2})}^{4}}+\displaystyle \frac{{E}_{{\rm{l}}}}{1+{s}^{2}/{R}_{{\rm{N}}}^{2}},\end{eqnarray} \tag{ 8 }$$

where s is the vortex-nucleus separation, and E_s (E_l) the short-range (long-range) contribution (Epstein & Baym 1988). The potential falls off on the scale of the nuclear radius R_N, corresponding to ℓ ≃ R_N and

$$\begin{eqnarray}&&{{ \mathcal R }}_{\mathrm{EB}}\simeq 1.4{\left(\displaystyle \frac{\mu }{{\hslash }}\right)}^{1/2}{\left(\displaystyle \frac{{E}_{{\rm{p}}}}{{\rho }_{{\rm{s}}}\kappa }\right)}^{2}\displaystyle \frac{1}{{\rm{\Delta }}{v}^{3/2}}\displaystyle \frac{{n}_{l}}{{R}_{{\rm{N}}}^{1/2}},\end{eqnarray} \tag{ 9 }$$

which is in agreement with the Epstein & Baym (1992) scalings, and we include an appropriate numerical prefactor. Performing a more detailed analysis of the Kelvin wave excitation process and employing an interaction potential of the form (8), E_p is found to be a mixture of E_s and E_l, with coefficients that depend on the scalings of each term with s. We obtain

$$\begin{eqnarray}&&{E}_{{\rm{p}}}^{2}\simeq {E}_{{\rm{s}}}^{2}+{E}_{{\rm{l}}}{E}_{{\rm{s}}}+0.5{E}_{{\rm{l}}}^{2}.\end{eqnarray} \tag{ 10 }$$

Note that these coefficients, as well as the numerical prefactor in Equation (9), disagree with the results of Epstein & Baym (1992). Repeating the calculation outlined in their Appendix B and Section 3, we determine drag coefficients that are about one order of magnitude smaller than those corresponding to Equation (3.18) in Epstein & Baym (1992). We trace the disagreement and different coefficients in Equation (10) back to an erroneous integration in the energy associated with the Kelvin wave excitations and/or power dissipated.³

The second study of Kelvin wave dynamics adopts a different prescription for the vortex-nucleus interaction: according to Jones (1992), this process dissipates the power p ∼ ΔE/τa per unit length, where a denotes the bcc lattice constant. The drag coefficient now reads

$$\begin{eqnarray}&&{ \mathcal R }\sim {\left(\displaystyle \frac{\mu }{{\hslash }}\right)}^{1/2}{\left(\displaystyle \frac{{E}_{{\rm{p}}}}{{\rho }_{{\rm{s}}}\kappa }\right)}^{2}\displaystyle \frac{1}{{\rm{\Delta }}{v}^{3/2}}\displaystyle \frac{1}{a{{\ell }}^{5/2}}.\end{eqnarray} \tag{ 11 }$$

Further, Jones (1992) does not account for a long-range contribution and uses a short-range potential

$$\begin{eqnarray}&&E(s)={E}_{{\rm{p}}}\exp \left(-\displaystyle \frac{{s}^{2}}{2{\xi }^{2}}\right)\end{eqnarray} \tag{ 12 }$$

that falls off on a much larger scale, the coherence length ξ. The appropriate choice is now ℓ ≃ ξ, and we find

$$\begin{eqnarray}&&{{ \mathcal R }}_{{\rm{J}}}\simeq \displaystyle \frac{1}{2\sqrt{\pi }}{\left(\displaystyle \frac{\mu }{{\hslash }}\right)}^{1/2}{\left(\displaystyle \frac{{E}_{{\rm{p}}}}{{\rho }_{{\rm{s}}}\kappa }\right)}^{2}\displaystyle \frac{1}{{\rm{\Delta }}{v}^{3/2}}\displaystyle \frac{1}{a{\xi }^{5/2}},\end{eqnarray} \tag{ 13 }$$

reproducing the scalings of Jones (1992) and we added his numerical prefactor.

Equations (9) and (13) would suggest that the dissipation associated with Kelvin wave excitations depends sensitively on E_p and Δv. These quantities are, however, not independent. On small scales, vortex unpinning is initiated once the Magnus force, generated by the background superfluid, exceeds the pinning force. As these dynamics are governed by the forces acting on a vortex per unit length, microscopic pinning interactions have to be modified to account for the finite length of vortices. In the inner crust, these structures remain straight over a distance $L\sim {10}^{3}\,{R}_{\mathrm{WS}}$ (Seveso et al. 2016) and interact with many randomly orientated nuclei. By geometrically averaging over this mesoscopic scale, Seveso et al. (2016) determine a decrease in the pinning force per unit length by about two orders of magnitude in agreement with Jones (1990, 1992). We include this by accounting for a constant reduction factor, δ ≈ 10⁻², and introduce effective pinning energies E_p → E_pδ. By balancing the Magnus force and pinning force per unit length, the critical velocity lag Δv_cr between a vortex and the background superfluid flow at which unpinning takes place can thus be related to the microscopic parameters characterizing the pinning interaction. The local relative velocity Δv between a free vortex and the nuclear lattice is typically of the same order as Δv_cr, and we estimate

$$\begin{eqnarray}&&{\rm{\Delta }}v\simeq {\rm{\Delta }}{v}_{\mathrm{cr}}=\displaystyle \frac{f}{{\rho }_{{\rm{s}}}\kappa }\sim \displaystyle \frac{{E}_{{\rm{p}}}\delta }{{la}{\rho }_{{\rm{s}}}\kappa }.\end{eqnarray} \tag{ 14 }$$

Substituting this into Equations (9) and (13) gives

$$\begin{eqnarray}&&{{ \mathcal R }}_{\mathrm{EB}}\simeq 2.8{\left(\displaystyle \frac{\mu }{{\hslash }}\right)}^{1/2}{\left(\displaystyle \frac{{E}_{{\rm{p}}}\delta }{{\rho }_{{\rm{s}}}\kappa }\right)}^{1/2}\displaystyle \frac{{R}_{{\rm{N}}}}{{a}^{3/2}},\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{{ \mathcal R }}_{{\rm{J}}}\simeq \displaystyle \frac{1}{2\sqrt{\pi }}{\left(\displaystyle \frac{\mu }{{\hslash }}\right)}^{1/2}{\left(\displaystyle \frac{{E}_{{\rm{p}}}\delta }{{\rho }_{{\rm{s}}}\kappa }\right)}^{1/2}\displaystyle \frac{{a}^{1/2}}{\xi }.\end{eqnarray} \tag{ 16 }$$

The two expressions differ by

$$\begin{eqnarray}&&\displaystyle \frac{{{ \mathcal R }}_{{\rm{J}}}}{{{ \mathcal R }}_{\mathrm{EB}}}\simeq 0.1\,\displaystyle \frac{{a}^{2}}{{R}_{{\rm{N}}}\xi }.\end{eqnarray} \tag{ 17 }$$

In the next section, we calculate these coefficients for a realistic crust model and show that the different choices for the vortex-nucleus interaction affect the strength of the crustal mutual friction.

Microscopic parameters for five inner crustal regions are summarized in Table 1. These are based on the EoS of Negele & Vautherin (1973), calculated in the Wigner–Seitz approximation. Adopting this ground-state composition, several authors have studied the vortex-nucleus interaction by analyzing the energy gain/loss of superimposing a vortex and a single lattice site. Table 1 shows estimates for R_N, E_s and E_l determined by Epstein & Baym (1992) using a Ginzburg–Landau approximation. Whereas E_s changes sign and can be repulsive or attractive, E_l is of hydrodynamical origin (related to the change in the superfluid's kinetic energy due to the nucleus as a result of the Bernoulli effect) and always repulsive (Shaham 1980; Epstein & Baym 1988), leading to different pinning geometries. More recently, Donati & Pizzochero (2006) have studied the pinning problem within a semi-classical model. Their estimates for Δ, ξ and E_p are also shown in Table 1. Due to the competition of the superfluid's internal, kinetic and condensation energy (Donati & Pizzochero 2004), they also obtain distinct configurations: interstitial pinning at lower densities (E_p > 0 implies vortex-nucleus repulsion) and nuclear pinning at high densities (vortex-nucleus attraction due to E_p < 0). Additionally, Donati & Pizzochero (2006) argue that the pinning strength decreases significantly toward the crust-core boundary due to collective pinning: each vortex contains several nuclei because the Wigner–Seitz radius R_WS is smaller than ξ, effectively weakening the interaction.

Table 1.  Composition for Five Crustal Domains and Corresponding Vortex-nucleus Interaction Parameters

	nb	Z	N	$\tilde{x}$	ns	ρ	A	RWS	nl	a	RN	Es	El	Δ	ξ	Ep
	(10−4×				(10−4×	(1012×			(10−6×
	fm−3)				fm−3)	g cm−3)		(fm)	fm−3)	(fm)	(fm)	(MeV)	(MeV)	(MeV)	(fm)	(MeV)
I	8.8	40	280	0.53	4.8	1.5	115	44.3	2.7	90.0	5.9	0.42	0.16	0.21	15.6	0.21
II	57.7	50	1050	0.45	47.0	9.6	161	35.7	5.2	72.5	6.7	−0.13	0.94	0.68	10.1	0.29
III	204.0	50	1750	0.35	184.0	33.9	193	27.6	11.3	56.1	7.2	−1.64	1.40	0.91	12.0	−2.74
IV	475.0	40	1460	0.28	436.0	78.9	183	19.6	31.7	39.8	7.3	−1.00	1.00	0.56	26.1	−0.72
V	789.0	32	950	0.16	737.0	131.0	232	14.4	80.3	29.2	7.2	−0.78	0.49	0.19	90.8	−0.02

Notes. Baryon density n_b, proton number Z, total neutron number N within a Wigner–Seitz sphere, proton-to-neutron ratio $\tilde{x}$, inside a nucleus and free neutron density n_s, are taken from Negele & Vautherin (1973). We calculate the total mass density ρ ≃ m_un_b, number of baryons inside a nucleus $A\simeq Z(1+1/\tilde{x})$, Wigner–Seitz radius ${R}_{\mathrm{WS}}\simeq {[3(N+Z)/(4\pi {n}_{{\rm{b}}})]}^{1/3}$, lattice nucleus density ${n}_{{\rm{l}}}\simeq 3/(4\pi {R}_{\mathrm{WS}}^{3})$ and lattice constant $a\simeq {(2/{n}_{{\rm{l}}})}^{1/3}$. Estimates for the nuclear radius R_N and short-range (long-range) contribution E_s (E_l) to the pinning interaction are from Epstein & Baym (1992).^{a, b} neutron gap Δ and coherence length ξ, related as $\xi ={{\hslash }}^{2}{k}_{\mathrm{Fs}}/(\pi {m}_{{\rm{u}}}{\rm{\Delta }})$ (k_Fs is the free neutron Fermi wave number), and microscopic pinning energies E_p (corresponding to β = 3) are taken from Donati & Pizzochero (2006)^c.

^aDonati & Pizzochero (2006) use different values for R_N in their calculation of E_p. We choose the parameters of Epstein & Baym (1992) because they are in better agreement with the results of Negele & Vautherin (1973). ^bWe follow Epstein & Baym (1992) and choose a short-range contribution that is reduced by a factor 10. ^cNote that Donati & Pizzochero (2006) use a different definition of ξ, resulting in a factor $1/\sqrt{6}$ instead of 1/π.

Download table as:  ASCIITypeset image

Combining these microscopic parameters with Equation (14), the relative vortex-nucleus velocities and pinning forces can be evaluated for the five crustal layers. To assess the impact of different assumptions about the pinning interaction, and allow a comparison between the formalisms of Epstein & Baym (1992) and Jones (1992), we calculate Δv and f for three distinct cases: (A) the full expression (10) including E_s,l together with l ≃ R_N, (B) E_p with l ≃ R_N, and (C) E_p with l ≃ ξ. The resulting estimates are given in Table 2. We observe that Δv changes significantly between the microscopic models and can differ by up to two orders of magnitude from the macroscopically averaged velocity lag Δv_av (see also Gügercinoglu & Alpar 2016).⁴ Note that according to Jones (1992), Kelvin wave dissipation is only effective if Δv ≳ 10² cm s⁻¹. For lower relative velocities, the energy loss proceeds via the excitations of lattice phonons, which is much weaker (Jones 1990, 1992). One velocity estimate in Table 2 (Δv ≈ 50 cm s⁻¹ for case (C) in zone V) drops slightly below this limit, suggesting that strong Kelvin wave drag might be replaced by weaker phonon drag close to the crust-core interface. As all but one Δv-estimates indicate a Kelvin wave-dominated regime, we focus on the friction coefficients presented in Section 2 and postpone an analysis of the impact of phonon-related dissipation to future work. Table 2 further provides the pinning forces per unit length given by Seveso et al. (2016). Comparison with our estimates shows that accounting for the geometric average over a mesoscopic vortex segment via a reduction factor δ gives reasonable agreement and discrepancies can be attributed to different microscopic treatments of the pinning energy and the corresponding interaction length scales.

To evaluate the drag coefficients of Section 2, we require the effective mass. Because the vortex tension can be approximated with $T\simeq -{\rho }_{{\rm{s}}}{\kappa }^{2}\mathrm{ln}(k\xi )/4\pi$ (Sonin 1987), we obtain μ(k) ≃ −2m_u/ln kξ.⁵ To simplify the calculation, we neglect the density-dependence of μ and evaluate the logarithmic factor at a characteristic wave number ${k}_{* }\approx 1.5\times {10}^{-3}\,{\mathrm{fm}}^{-1}$, and a typical coherence length, ξ ≈ 12 fm, which gives μ ≃ μ(k_*) ≈ m_u/2. The error introduced by keeping μ constant is small, due to the weak dependence on k and ξ. Using fiducial values for domain III, Equations (15) and (16) yield

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal R }}_{\mathrm{EB}} & \approx & 3.5\times {10}^{-3}{\left(\displaystyle \frac{| {E}_{{\rm{p}}}| }{3\mathrm{MeV}}\right)}^{1/2}{\left(\displaystyle \frac{\delta }{{10}^{-2}}\right)}^{1/2}\\ & & \ \times {\left(\displaystyle \frac{2\times {10}^{-2}{\mathrm{fm}}^{-3}}{{n}_{{\rm{s}}}}\right)}^{1/2}{\left(\displaystyle \frac{56\mathrm{fm}}{a}\right)}^{3/2}\left(\displaystyle \frac{{R}_{{\rm{N}}}}{7\,\mathrm{fm}}\right)\end{array}\end{eqnarray} \tag{ 18 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal R }}_{{\rm{J}}} & \approx & 1.3\times {10}^{-2}{\left(\displaystyle \frac{| {E}_{{\rm{p}}}| }{3\mathrm{MeV}}\right)}^{2}{\left(\displaystyle \frac{\delta }{{10}^{-2}}\right)}^{2}\\ & & \ \times {\left(\displaystyle \frac{2\times {10}^{-2}{\mathrm{fm}}^{-3}}{{n}_{{\rm{s}}}}\right)}^{2}{\left(\displaystyle \frac{a}{56\mathrm{fm}}\right)}^{1/2}\left(\displaystyle \frac{12\,\mathrm{fm}}{\xi }\right).\end{array}\end{eqnarray} \tag{ 19 }$$

As before, we compare the assumptions about the microphysics in the formalisms of Epstein & Baym (1992) and Jones (1992) (long- plus short-range interactions versus short-range interaction only) by determining the friction coefficient in three different ways: (A) ${ \mathcal R }$_EB calculated with the full expression (10) including E_s,l, (B) ${ \mathcal R }$_EB calculated with E_p only, and (C) ${ \mathcal R }$_J calculated with E_p. We refer to these cases with the labels (A), (B), and (C).

The resulting profiles of ${ \mathcal B }$ are illustrated in Figure 1.⁶ We employ a spline function to interpolate results for the five domains from neutron drip at ${\rho }_{{\rm{D}}}\approx 4.0\times {10}^{11}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ to the crust-core interface at ${\rho }_{\mathrm{cc}}\approx 1.3\times {10}^{14}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. The fit gives unphysical results when extrapolating below 1.5 × 10¹² g cm⁻³; instead, we take ${ \mathcal B }$ constant for simplicity. We observe that mutual friction varies strongly with density and differs significantly close to the crust-core interface. This region carries the majority of the superfluid's mass as shown in the bottom panel of Figure 1, where ${ \mathcal B }$ is given as a function of relative mass fraction ΔM/M. Here, ΔM denotes the overlying mass and M the total mass taken to be M ≈ 1.41 M_⊙ (M_⊙ is the solar mass) in our model. The range of ${ \mathcal B }$ suggests different post-glitch behavior.

Zoom In Zoom Out Reset image size

To analyze the effects of density-dependent friction on the post-glitch response, we use a simple time-dependent, three-component model to determine the shape of the glitch rise. The star is decomposed into a crust neutron superfluid, core neutron superfluid and a non-superfluid "crust" component, representing the nuclear lattice and tightly coupled charged conglomerate in the core (Easson 1979). Angular velocities and moments of inertia are denoted by Ω_x and I_x with x ∈ {sf,core,crust}, respectively. We also assume the crust and core components to be rigidly rotating, and incorporate crust-core coupling by assigning a constant mutual friction coefficient ${ \mathcal B }$_core. To account for uncertainties in our understanding of the underlying mechanism, we determine the glitch rise for two fiducial values. If the dynamics are dominated by the scattering of electrons off magnetized vortices (Alpar et al. 1984a), ${{ \mathcal B }}_{\mathrm{core}}\approx 5\times {10}^{-5}$ (Andersson et al. 2006) is a suitable choice. Stronger friction could be present if the crust-core coupling is mediated by the interactions between vortices and superconducting fluxtubes (Link 2003; Sidery & Alpar 2009), provided that the core protons form a type-II state (Baym et al. 1969). To study this possibility, we follow Link (2003) and Haskell et al. (2014) and assume that Kelvin waves are excited along the vortices as they cut through fluxtubes. The formalism discussed in Section 2 can be directly translated provided that the microphysics are adjusted. The vortex-fluxtube interaction is predominantly magnetic, resulting in pinning energies on the order of E_p ≈ 5 MeV (Link 2003) active over a length scale l ≃ λ_* ≈ 100 fm, i.e., the London penetration depth. See e.g., Graber et al. (2017) for typical length scale estimates. For a fiducial magnetic field strength of B ≈ 10¹² G, the inter-fluxtube distance is ${d}_{\mathrm{ft}}\approx {10}^{3}\,\mathrm{fm}$, leading to pinning forces per unit length of f ≈ 8 × 10¹⁵ dyn cm⁻¹. As for the crust, a local force balance provides the means to estimate the relative vortex-fluxtube velocity to Δv ≈ 4 × 10⁴ cm s⁻¹ for ${\rho }_{{\rm{s}}}\approx {10}^{14}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. For an effective mass of μ ≈ m_u/2, this subsequently gives ${{ \mathcal B }}_{\mathrm{core}}\simeq {{ \mathcal R }}_{\mathrm{core}}\approx {10}^{-2}$.

Generalizing the results of Haskell & Melatos (2015) to three components as well as neglecting entrainment, the equations of motion read

$$\begin{eqnarray}&&{\dot{{\rm{\Omega }}}}_{\mathrm{sf}}={ \mathcal B }\left[2{{\rm{\Omega }}}_{\mathrm{sf}}+\tilde{r}\,\displaystyle \frac{\partial {{\rm{\Omega }}}_{\mathrm{sf}}}{\partial \tilde{r}}\right]({{\rm{\Omega }}}_{\mathrm{crust}}-{{\rm{\Omega }}}_{\mathrm{sf}}),\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{\dot{{\rm{\Omega }}}}_{\mathrm{core}}=2{{ \mathcal B }}_{\mathrm{core}}{{\rm{\Omega }}}_{\mathrm{core}}\,({{\rm{\Omega }}}_{\mathrm{crust}}-{{\rm{\Omega }}}_{\mathrm{core}}),\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{\dot{{\rm{\Omega }}}}_{\mathrm{crust}}=-\displaystyle \frac{{N}_{\mathrm{ext}}}{{I}_{\mathrm{crust}}}-\displaystyle \frac{{I}_{\mathrm{core}}}{{I}_{\mathrm{crust}}}\,{\dot{{\rm{\Omega }}}}_{\mathrm{core}}-\displaystyle \frac{\int \rho {\tilde{r}}^{2}{\dot{{\rm{\Omega }}}}_{\mathrm{sf}}\,{dV}}{{I}_{\mathrm{crust}}},\end{eqnarray} \tag{ 22 }$$

where N_ext is the external spin-down torque, $\tilde{r}$ the cylindrical radius, and the integral is performed over the inner crust. For simplicity, a cylindrical geometry is used: we solve the problem in the equatorial plane and rescale the results so that the total crustal moment of inertia in cylindrical coordinates matches that in spherical ones. We assume a total moment of inertia of ${I}_{\mathrm{tot}}\approx 0.35\,{{MR}}^{2}$ (Lattimer & Prakash 2001), with the core neutrons and charged particles constituting 95% and 5% of the core's moment of inertia, respectively. To relate ρ and $\tilde{r}$ in the crust, we integrate the TOV equations assuming a core radius and mass of 10 km and 1.4 M_⊙. For consistency, we consider the Negele & Vautherin (1973) EoS for the inner crust. As this EoS does not apply in the outer crust, we take the pressure below neutron drip to be dominated by relativistic electrons with Y_e ≈ 0.4.

Using initial conditions which are typical for the Vela pulsar, i.e., ${{\rm{\Omega }}}_{\mathrm{crust},\mathrm{core}}(0)\,\approx \,70.34\,\mathrm{rad}\,{{\rm{s}}}^{-1}$ and ${\rm{\Delta }}{{\rm{\Omega }}}_{\mathrm{crit}}\equiv {{\rm{\Omega }}}_{\mathrm{sf}}(0)\,-{{\rm{\Omega }}}_{\mathrm{crust}}(0)\approx 6.3\times {10}^{-3}\,\mathrm{rad}\,{{\rm{s}}}^{-1}$ to allow a comparison with the glitch observations of Palfreyman et al. (2018), we evolve the equations of motion (20)–(22) for $120\,{\rm{s}}$ to encompass observational constraints on the spin-up timescale (Dodson et al. 2002, 2007; Palfreyman et al. 2018). Note that our crustal mutual friction profiles can vary in space, but remain constant over time. We further ignore the external spin-down torque when integrating the equations of motion, as it has negligible effect on the short-term post-glitch response within $120\,{\rm{s}}$.

As illustrated by the characteristic shape of ${{\rm{\Omega }}}_{\mathrm{sf}}(t,\tilde{r})$ in Figure 2, the superfluid's differential rotation is mainly driven by the ${ \mathcal B }(\tilde{r})$-dependence.⁷ We show results for the drag profiles (A) and (C), strongest and weakest at high densities, respectively. For case (A), we observe that the superfluid starts to couple within ∼100 ms, eventually transferring all its excess angular momentum and spinning down to a new steady state, where all three components are corotating. The superfluid's evolution looks qualitatively similar for the other two drag profiles, albeit exhibiting stronger differential rotation because they cover a larger range of mutual friction strengths than the profile (A). Moreover, for cases (B) and (C), the bottom of the crust recouples on longer timescales as ${ \mathcal B }$ is weaker in this region. Note that the angular momentum transfer from the superfluid to the crust component is not completed within $120\,{\rm{s}}$ for model (C). Figure 2 also indicates that for model (A), a stronger ${ \mathcal B }$_core causes the superfluid to reach the equilibrium faster. The impact of different core couplings on the superfluid's differential rotation is less pronounced for cases (B) and (C).

Zoom In Zoom Out Reset image size

To illustrate the effects of density-dependent crustal profiles and a variable crust-core coupling strength on observables, we compute the change in crustal frequency with time. For comparison, Δν is also determined for four constant coefficients, ${ \mathcal B }$ ≈ 10⁻¹, 10⁻², 10⁻³, 10⁻⁴. Results for Δν(t) are shown in the left panels of Figure 3, highlighting the effect of the relative strength between the crust and the core couplings: for ${{ \mathcal B }}_{\mathrm{core}}\approx 5\times {10}^{-5}$, angular momentum transfer from the superfluid to the crust is very effective and acting on timescales shorter than the crust-core coupling timescale, ${\tau }_{{cc}}\approx 7.5\,{\rm{s}}$. This causes the crust's rotation frequency to increase above the asymptotic value, generating a characteristic "overshoot." As soon as crust and core recouple, Δν decreases and eventually approaches the new steady state. For the stronger scenario, ${{ \mathcal B }}_{\mathrm{core}}\approx {10}^{-2}$, crust-core coupling proceeds on τ_cc ≈ 37.4 ms, which is faster than the density-dependent crust couplings. Thus, the superfluid transfers angular momentum to the combined crust-core system, resulting in a slower, monotonic rise of Δν. Note that as illustrated in the top-right panel of Figure 3, the onset of crust-core coupling in the case ${{ \mathcal B }}_{\mathrm{core}}\approx 5\times {10}^{-5}$ is clearly visible as a break in the phase shift ϕ, accumulating after the glitch. For ${{ \mathcal B }}_{\mathrm{core}}\approx {10}^{-2}$, however, the break in phase shift ϕ moves to the left and becomes basically invisible. As we will illustrate in the next section, observing such a feature in ϕ could provide important information about the frictional processes in the core.

Zoom In Zoom Out Reset image size

Assuming that the superfluid reservoir is completely depleted, angular momentum conservation dictates for the equilibrium lag

$$\begin{eqnarray}&&{\rm{\Delta }}{\nu }_{\mathrm{equi}}=\displaystyle \frac{{I}_{\mathrm{sf}}}{{I}_{\mathrm{total}}}\,\displaystyle \frac{{\rm{\Delta }}{{\rm{\Omega }}}_{\mathrm{crit}}}{2\pi }\approx 16.0\,\mu \mathrm{Hz},\end{eqnarray} \tag{ 23 }$$

which is in agreement with the glitch step size determined by Palfreyman et al. (2018) for the 2016 Vela pulsar glitch. Details of the approach to steady state depend crucially on the strength of ${ \mathcal B }$ close to the crust-core interface. Note that for profile (C), weakest at high densities, the crust does not reach the new equilibrium within the $120\,{\rm{s}}$ integration window and additional time would be required for the superfluid to transfer all its angular momentum (also illustrated in the right panels of Figure 2). Such a slow recovery could in principle be misinterpreted as the new spin-equilibrium, leading to incorrect initial conditions or moment of inertia estimates.

The first pulse-to-pulse glitch observation was recently published by Palfreyman et al. (2018) for a glitch that occurred on 2016 December 12 in the Vela pulsar. To test the potential of our model as a tool to constrain microphysics, a preliminary comparison between the predictions and the new data is presented. Because various processes introduce noise into single-pulse observations (e.g., Shannon et al. 2014), we average the timing residuals (the difference between the observed pulse arrival times and those expected from a timing model) into $2\,{\rm{s}}$ bins and center the data around the glitch epoch t_g given by Palfreyman et al. (2018). We subsequently determine the timing residuals corresponding to our predicted spin-up. Provided that residuals are small, they are proportional to −ϕ. As demonstrated by Figure 3 (right panels), the residuals thus start at zero and become increasingly negative with time, reproducing the characteristic glitch signature. The observation however reveals positive (approximately constant) timing residuals after t_g. In concordance with the magnetospheric changes accompanying the glitch (Palfreyman et al. 2018), we do not interpret this as a spin-down of the pulsar but instead as a phase shift. We include this by applying a constant shift Δt to the residuals, so that theoretical predictions and observation agree at t = 0. A comparison between the resulting timing residuals as well as the cumulative residuals is shown in Figure 4 for the first $120\,{\rm{s}}$ after the glitch based on ${{ \mathcal B }}_{\mathrm{core}}\approx 5\times {10}^{-5}$.

Zoom In Zoom Out Reset image size

Figure 4 highlights that the initial post-glitch response is rather insensitive to the crustal mutual friction profile provided that ${ \mathcal B }$ ≳ 10⁻³. We observe that models (A), (B) and the constant couplings ${ \mathcal B }$ ≈ 10⁻¹, 10⁻², 10⁻³ fit the data similarly well. Only profile (C) (the weakest close to the crust-core interface) and ${ \mathcal B }$ ≈ 10⁻⁴ decrease slower than what is observed, suggesting that strong mutual friction prevails in a sizable fraction of the inner crust. Together with the fiducial choice ${{ \mathcal B }}_{\mathrm{core}}\approx 5\times {10}^{-5}$ this ensures that a large portion of the superfluid's angular momentum can be transferred to the crust component before the core recoupling takes place. To illustrate the interplay of the crust and core coupling strengths in more detail, we focus on the strongest coupling profile (A) and determine the glitch response for a range of crust-core coupling coefficients. As depicted in Figure 5, the neutron star's rotational evolution is very sensitive to ${ \mathcal B }$_core. Disagreement between the model predictions and data is amplified as soon as ${ \mathcal B }$_core diverges from the fiducial value: For stronger (weaker) mutual friction, the phase shifts become much smaller (larger), which results in smaller (larger) timing residuals. Our comparison thus suggests that the dominant core mutual friction mechanism covers a rather narrow range $3\,\times {10}^{-5}\lesssim {{ \mathcal B }}_{\mathrm{core}}\lesssim {10}^{-4}$, as typical for electron scattering off magnetized vortices (Alpar et al. 1984a).

Zoom In Zoom Out Reset image size

For the first time, we calculate mutual friction profiles resulting from Kelvin wave excitation for a realistic crust model and combine those with a simplified treatment of the crust-core coupling to develop a predictive model of the glitch rise. We find that density-dependent coupling affects the amount of angular momentum that can be exchanged on specific timescales and hence influences the glitch response of the crust. This illustrates that uncertainties in deriving the underlying ${ \mathcal B }$ and microscopic parameters have a crucial influence on observables.

We demonstrate that the ${ \mathcal B }$ profiles depend most sensitively on the assumed vortex-nucleus interaction. Model (A) accounts for the contributions E_s,l included by Epstein & Baym (1992), which remain almost constant at high densities. This causes stronger drag and thus faster recoupling. For the profiles (B) and (C), we instead considered E_p, which decreases significantly with density due to collective pinning, and results in longer coupling timescales. Nonetheless, differences remain between the glitch rise predictions based upon the formalisms of Epstein & Baym (1992) and Jones (1992) due to their respective assumptions on the interaction potentials and dissipation length scales (see Figure 1).

Other microphysical parameters of the crust also play an important role. Whereas the composition itself does not vary significantly between different EoSs, our results are sensitive to superfluid parameters such as the energy gap and in principle entrainment, which we have neglected to keep our introductory analysis tractable. Strong entrainment would reduce the size of the crustal angular momentum reservoir, causing difficulties for the "crust-only" glitch framework (Andersson et al. 2012; Chamel 2013) (see however Watanabe & Pethick 2017). Future work will be needed to address how entrainment impacts on the initial glitch response. Our results are further strongly affected by the pinning strength. Calculations of these parameters rely on many assumptions and are very uncertain: Whereas Epstein & Baym (1988) and Donati & Pizzochero (2006) employed a Ginzburg–Landau approach and semi-classical model, respectively, Avogadro et al. (2008) have examined the vortex-nucleus interaction using a quantum mean-field framework arriving at pinning energies of opposite signs. Future work is essential to reconcile these results. A correct description of vortex transport should also account for interactions with a nuclear pasta phase expected to be present close to the crust-core interface (Ravenhall et al. 1983). This high-density region carries the majority of the crustal mass and should strongly affect the post-glitch behavior. Real-time studies of the vortex-nucleus interaction (Bulgac et al. 2013; Wlazłowski et al. 2016) could help to address this issue, but it remains unclear how this microscopic picture relates to the dynamics of a mesoscopic vortex communicating with many nuclei.

Finally, note that we based our model on the assumption that Kelvin wave excitations dominate the dissipation. Other processes, such as vortex coupling to lattice defects or impurities (Harding et al. 1978) and dissipation due to lattice phonon excitations (Jones 1990, 1992), could similarly alter the glitch response and their effects studied as outlined above once the mutual friction profile is known.

In addition to crustal microphysics, the shape of the glitch rise is crucially influenced by the relative strength between crust coupling and core mutual friction. The amount of angular momentum that the superfluid transfers to the crust before the core is recoupled controls the size of the phase shifts, providing the means to constrain the core physics. This plays an important role in comparing our predictive model with the first resolved glitch rise observation of the 2016 December Vela glitch (Palfreyman et al. 2018). Although a more detailed analysis will be needed to systematically study the impact of the underlying microscopic parameters on the glitch rise, our comparison points toward strong crustal mutual friction satisfying ${ \mathcal B }$ ≳ 10⁻³ in combination with weaker core coupling in the range $3\times {10}^{-5}\lesssim {{ \mathcal B }}_{\mathrm{core}}\lesssim {10}^{-4}$. Such strengths are typical for electron scattering off the magnetized vortices (Alpar et al. 1984a; Andersson et al. 2006), but much weaker than the drag associated with excitations of vortex Kelvin waves in the neutron core (Link 2003; Haskell et al. 2014). The absence of strong dissipation (characteristic for the regime where vortex-fluxtube interactions dominate the dynamics) could be explained if the protons do not form a type-II superconductor. Coupling dynamics in a type-I state are however rather uncertain (Sedrakian 2005; Jones 2006). Furthermore, our predictive model only accounts for constant ${ \mathcal B }$_core values, and additional work incorporating density-dependent crust-core coupling would be needed to verify an absence of strong core friction. Our conclusions were further based on the assumption that the Vela pulsar undergoes a shift in phase at the time of the glitch. Future observations will be required to confirm if this is justified and the phase shift is indeed a real feature of pulsar glitches. Upcoming facilities like the Square Kilometer Array will play an important role in this endeavor as they may allow the glitch rises of other sources to be observed (Watts et al. 2014; Kramer & Stappers 2015).

We thank the referee for suggesting an improved treatment of the microscopic vortex velocity and Palfreyman et al. (2018) for making their data publicly available. This work also benefited from discussions with Robert Archibald, Evan Keane and Toby Wood. V.G. is supported by a McGill Space Institute postdoctoral fellowship and the Trottier Chair in Astrophysics and Cosmology. A.C. is supported by an NSERC Discovery grant, is a member of the Centre de Recherche en Astrophysique du Québec (CRAQ), and thanks Newcastle University for hospitality. N.A. acknowledges funding from STFC in the UK through grant number ST/M000931/1.

Software: IPython (Perez & Granger 2007), Matplotlib (Hunter 2007), NumPy (Oliphant 2006), Pandas (McKinney 2010), SciPy (Jones et al. 2001).