Revisiting the post-glitch relaxation of the 2000 Vela glitch with the neutron star equation of states in the Brueckner and relativistic Brueckner theories

We revisit the short-term post-glitch relaxation of the Vela 2000 glitch in the simple two-component model of pulsar glitch by making use of the latest realistic equations of states from the microscopic Brueckner and the relativistic Brueckner theories for neutron stars, which can reconcile with the available astrophysical constraints. We show that to fit both the glitch size and the post-glitch jumps in frequency derivatives approximately one minute after the glitch, the mass of the Vela pulsar is necessarily small, and there may be demands for a stiff equation of state (which results in a typical stellar radius larger than $\sim\negthickspace12.5\,\rm km$) and a strong suppression of the pairing gap in the nuclear medium. We discuss the implications of this result on the understanding of pulsar glitches.


INTRODUCTION
Glitches, i.e., the sudden spin-up of pulsars with subsequent recovery, were discovered during pulsar timing studies in the Vela pulsar (PSR J0835-4510 or PSR B0833-45) more than fifty years ago (Radhakrishnan & Manchester 1969). They are generally agreed to be the result of the complicated interplay between the superfluid component and the rest of the star (Baym et al. 1969;Anderson & Itoh 1975), i.e., the two-component framework. Glitch observations are important for understanding the properties of neutron stars and provide a unique way to probe neutron star internal structure and dynamics. Many unknown questions concerning both the trigger and the relaxation of a glitch are closely related to the microphysical state of the matter and hence to the equation of state (EOS) and the structure of the neutron stars.
With a total of 22 glitches in about 50 years, the Vela pulsar is one of the most frequently glitching pulsars of over 2000 known pulsars. Among the many efforts made to try to unravel the glitches in Vela-like pulsars, it has been shown that the sizes and subsequent recoveries can be understood to some extent from the hydrodynamic Corresponding author: Ang Li liang@xmu.edu.cn evolution of the vortices (Sidery & Alpar 2009;Sidery et al. 2010;van Eysden & Melatos 2010;Pizzochero 2011;Haskell et al. 2012;Seveso et al. 2012;Howitt et al. 2016;Haskell et al. 2018;Gügercinoglu & Alpar 2020); furthermore, a "snowplow" model has been proposed (Pizzochero 2011), which incorporates microscopically based pinning forces throughout the stellar crust. Seveso et al. (2012) extended the snowplow model from a polytropic EOS to realistic EOSs obtained in the popular relativistic mean-field (RMF) model. They found that the model could reproduce the main features of giant glitches in the Vela pulsar, including the size and post-glitch jump infrequency derivatives of the January 2000 glitch (Dodson et al. 2002). Note that a fourth exponential decay term (approximately 1 minute) was found, which represents the fastest decay term to date. Haskell et al. (2013) further extended the model to account for possible pinned vorticity in the core due to the interaction between vortices and flux tubes.
In those previous studies on short-term transient behaviours around glitch epoch, the proton fraction was adopted from a microscopic EOS in the Brueckner theory (Zuo et al. 2004), which is different from the RMF EOSs used for the pulsar. Due to the critical dependence of the star composition (including the proton fraction) on the EOS, in the present study, we intend to revisit the 2000 Vela glitch based on a consistent determination of both the EOS and the proton fraction from the microscopic (relativistic) Brueckner theory. Formerly, consistent EOS and composition on the basis of the Skyrme nuclear effective force have been adopted by Antonelli et al. (2018) for the study of the general relativistic corrections to pulsar glitch amplitudes.
The Brueckner theory has been proven to be one of the most advanced microscopic theories (Baldo 1999), as it allows an ab initio calculation of neutron star EOS from relativistic nucleon-nucleon (NN) potentials. In particular, the exact solution of the Brueckner-Bethe-Goldstone equation has been made available recently, which is free from the total momentum approximation that was commonly adopted in previous works (see more discussions in Shang et al. 2021). Furthermore, high-precision NN potentials (i.e., the pvCD-Bonn potentials) that are applicable in the relativistic Brueckner theory (Machleidt 1989) have been generated recently and can precisely describe the charge dependence in NN scattering data (see more discussions in Wang et al. 2019). The present calculations are done from both the Brueckner and relativistic Brueckner theories.
The aim of this paper is to better understand the microscopy and macroscopy of a neutron star with glitch observations. In particular, we aim to examine to what extent a simple two-component model of glitch can be applicable to observational constraints from the Vela pulsar. Section 2 is dedicated to formulating the twocomponent glitch model. Section 3 demonstrates the inner structure of a pulsar based on the microscopic EOSs within the Brueckner and relativistic Brueckner theories. Section 4 presents an application of the EOS and the glitch model to the Vela pulsar. Finally, we present a discussion and the conclusions in Section 5.

THE TWO-COMPONENT MODEL
In the simple two-component framework, one aspect is the charged component (labelled as index "p"), and the other is the superfluid component (labelled as index "n"). The visible p component includes crustal lattices, core protons, and neutrons (core neutrons can be partially coupled, as discussed later). The n component is usually assumed to be crustal neutrons, which may indeed be the case according to recent studies on the type of proton superconducting in the core (Haskell et al. , 2013 and on the interaction between neutron rotational vortices and proton flux tubes (Babaev 2009).
A glitch is triggered by the sudden transfer of angular momentum from n to p at some critical spin lag between the two. When the superfluid vortices are pinned to the crustal lattice, the charged component evolves according to the following (Glampedakis & Andersson 2011): where the term on the right side represents the standard torque due to a magnetic dipole (a is a coefficient). From balancing the pinning force to the Magnus force, one determines the critical spin lag at which the vortices will unpin (as shown in more details later in Section 4). Normally, the lag is small; we then have the following: where ∆t g is the interglitch time (∼2 yr for Vela), and τ c = −Ω p /2Ω p is the characteristic age of the pulsar (∼11.3 kyr for Vela). The fractional moment of inertia I n /I p is related to the glitch activity A g = (1/T )( ∆Ω p )/Ω p as follows: where T is the total data span for glitch monitoring, and (∆Ω g ) is the sum of all observed glitch frequency jumps. The glitch activity A g can be estimated for systems that have exhibited at least two glitches of similar magnitude, such as the Vela pulsar. Using the observed value of Vela's activity parameter A g = 1.91 × 10 −9 per day , the glitch observations have put a constraint on the fractional moment of inertia, which is I n /I p 1.6% (Link et al. 1999;Andersson et al. 2012;Montoli et al. 2021). It has been argued that the "entrainment" of the neutron superfluid by the crystalline structure of the crust greatly reduces its mobility; namely, Equation (3) is modified to the following (Andersson et al. 2012;Chamel 2013): where m * n is the average effective neutron mass. m * n /m n is in the range of ∼1.5−6 (Chamel 2005(Chamel , 2012(Chamel , 2017Watanabe & Pethick 2017).
The glitch rise time is usually very short. Presently, the best constraint comes from the 2016 Vela glitch, which is less than 12.6 seconds (Ashton et al. 2019). Correspondingly, it has been suspected that only a fraction of the core neutrons ( 30%, as seen later in Section 4) might be coupled to the charged component of the crust on glitch-rise time scales . Previously, Li (2015); Li et al. (2016) argued the presence of the glitch crisis in the framework of the microscopic Brueckner theory, namely, if when accounting for the entrainment effect, there is not enough angular momentum transferred to trigger the big Vela-like glitches from only the crustal superfluid neutrons. That is, Equation (4) cannot be guaranteed. Such calculations were performed under the assumption that the whole core part is coupled to the charged visible part during a glitch. Here, we relax this assumption and consider a more realistic scenario of a partially coupled core, which equivalently decreases the denominator I p on the right side of Equation (4). Therefore, the glitch crisis is currently not a problem. That the glitch crisis can be solved with the core superfluid had been already shown in Gügercinoglu & Alpar (2014); Haskell & Melatos (2015); Hooker et al. (2015).
We now present the inner structure of a pulsar based on the modern EOSs within the (relativistic) Brueckner theory and confront the theoretical results in the snowplow model with the only one available observation of the short-term recovery component of the Vela glitch (Dodson et al. 2002).

Brueckner-Hartree-Fock (BHF) approach
The Brueckner theory is based on the re-summation of the perturbation expansion of the ground-state energy of nuclear matter (Baldo 1999). For the present study, the EOS is calculated at the Brueckner-Hartree-Fock (BHF) level and adopts the realistic Argonne V 18 with a microscopic three-body force (Zuo et al. 2002). The details are described elsewhere (Zuo et al. 1999). Here, we outline it briefly.
In the BHF approach, the infinite series of certain diagrams in perturbation expansion for the total energy of the system is embodied in the so-called G-matrix operator, which is adopted instead of the bare nucleon interaction (with care taken for double counting of the same contributions). The main advantage of the G-matrix is that the elements of G do not diverge, as could happen when using a bare nucleon potential. Accordingly, one can start the calculation from a realistic nucleon-nucleon potential with no artificial parameters in the approach.
The G-matrix satisfies the Bethe-Goldstone equation: with the particle spectrum in BHF approach: where v is the bare nucleon potential, and Q is the Pauli operator, which represents the nuclear medium effect.
Here we remark the bare nucleon mass as m N and neglect the bare mass difference between neutrons and protons. U (p) is the auxiliary self-consistent potential and is determined by the following: where |12 A means that the wave function is properly antisymmetrized and p F 2 denotes the Fermi momentum. After several self-consistent iterations of Equations (5)−(7), the effective interaction matrix G is obtained. Using this G-matrix, the total energy per nucleon can be expressed as follows: Following this, the EOS of uniform nuclear matter, pressure P as a function of the density ρ, or energy density ε, can be calculated. The saturation properties and the symmetry energy parameters are found to be consistent with the empirical constraining bands. More results and discussions can be found in Shang et al. (2020Shang et al. ( , 2021.

Relativistic BHF approach
Different from the BHF approach, one needs a Dirac spinor instead of the free nucleon wave function (plane wave function) for the description of single particle motion in a nuclear medium within the RBHF approach: Here, M * N = m N + U S and E * 2 p = M * 2 N + p 2 . U S represents the scalar potential, and χ s is the Pauli spinor.
In the RBHF calculation, the G matrix, which serves as an effective interaction, can be obtained by solving the in-medium relativistic scattering equation. One of the most widely used covariant scattering equations is the Thompson equation (Thompson 1970), which is a relativistic three-dimensional reduction of the Bethe-Salpeter equation (Salpeter & Bethe 1951). In the rest frame of nuclear matter, the Thompson equation takes the following form: where E p is the eigenvalue of the Diarc equation in the medium. P = (k 1 + k 2 )/2 and k = (k 1 − k 2 )/2 are the total and relative momentum of the two interacting nucleons with momenta k 1 and k 2 , respectively. q, k, and q are the initial, intermediate, and final relative momenta of the two nucleons in the medium, respectively. The Pauli operator Q is the same as that used in the BHF approach. Using this G matrix, the self-energy for the positive energy solution can be calculated as follows: In addition, the expectation value of the single-particle potential for nucleons with momentum p can be expressed by the following: In general, the scalar and vector potentials U S and U V are momentum dependent; however, this dependence is very weak and can be neglected. The constants U S and U V can be adjusted to the positive energy solution in Equation (11) at the Fermi momentum. Accordingly, the eigenvalues E p can be obtained by solving the relativistic Hatree-Fock equation: where α = γ 0 γ and β = γ 0 are the Dirac matrices. Equations (10)−(13) should be solved self-consistently by iterating until their convergence. Afterwards, the binding energy per nucleon is evaluated by: (14) Then, the EOS of nuclear matter can be calculated in the same manner as that used in the BHF approach. As mentioned in the introduction, in the present RBHF calculations, we adopt charge-dependent Bonn potentials with pseudovector coupling between the pion and nucleon, pvCD-Bonn A (Wang et al. 2019), to evaluate the EOS. The saturation properties of nuclear matter for the two employed EOSs in the Brueckner and relativistic Brueckner theories are collected in Table 1, together with the constraints from terrestrial experiments (Shlomo et al. 2006;Piekarewicz 2010;Oertel et al. 2017;Li et al. 2019;Drischler et al. 2020;Reed et al. 2021).

Stellar structure
To study the structure of neutron stars, we must calculate the composition and EOS of charge-neutral neutron star matter, which consists of neutrons, protons, and leptons (e − , µ − ) in beta equilibrium . Once the EOS P (ε) is specified, the stable configurations of neutron stars can be obtained from solving the Tolman-Oppenheimer-Volkoff (TOV) equation for pressure P and the enclosed mass M : where G is the gravitational constant. The moment of inertia can be calculated under rigid-body rotation in Pizzochero (2011): Thus, the total moment of inertia is I tot = I(0, R), with R being the radius of the star. For later use, R ic is defined as the radius of the inner crust. The BPS (Baym et al. 1971) and NV (Negele & Vautherin 1973) EOS has been employed for nonuniform matter present in the outer and inner crust parts of the stars (below 0.08 fm −3 or ∼0.5ρ 0 ), respectively. That is, the pressures from the BHF/RBHF calculations are replaced by that of Negele & Vautherin (1973) from the nucleonic density 0.08 fm −3 . We mention here that the NV EOS is based on quantal Hartree-Fock calculations for Wigner-Seitz cells for the inner crust.
In Figure 1, we present the binding energy of pure neutron matter, as well as the EOS and the proton fraction of neutron star matter for the BHF and RBHF models described above, together with the mass-radius relations of stable neutron stars from solving the TOV equations. It is important to note here that the neutron star sequences such described, without any adjusted parameter, can, on the one hand, agree with the current constraints on the saturation properties of nuclear matter (as reported earlier in Table 1), as well as the results obtained for pure neutron matter within the chiral interactions at next-to-next-to-next-to-leading order (N3LO) level (Drischler et al. 2020), and on the other hand, meet the latest astrophysical mass and radius constraints from PSR J0030+0451 (Miller et al. 2019;Riley et al. 2019) and MSP J0740+6620 (Cromartie et al. 2020;Fonseca et al. 2021;Miller et al. 2021;Riley et al. 2021), as well as the binary tidal deformability constraint from GW170817 (Abbott et al. 2017(Abbott et al. , 2018. Based on these EOSs and consistently determined stellar composition from the BHF and RBHF models, in the following we shall proceed to discuss the glitch observations of the Vela pulsar in the snowplow model of Pizzochero (2011). Table 1.

Pinning force
In the snowplow scenario (Pizzochero 2011), the superfluid vorticity is assumed to accumulate in a thin sheet in the inner crust and expected to be released simultaneously at some cylindrical radius x max , when the critical lag between the two components reaches the maximum value ∆Ω cr max . The critical lag ∆Ω cr (x) is de- termined by the competition of the pinning force: and the Magnus force: i.e., ∆Ω cr ( Here, the vortex line should be parallel to the rotational axis of the star with a distance x. l(x) = 2 R 2 ic − x 2 represents the length of the vortex tube with the radius of inner crust R ic . The constant κ = π /m N is the quantum of circulation of a neutron fluid. The function ρ s (r) is the density of the superfluid component at radius r. We assume that the neutrons are superfluid throughout the star: ρ s = (1 − x p )ρ with x p being the proton fraction calculated within BHF/RBHF (see Figure 1).
To determine the critical lag, the pinning force per unit length f pin (ρ), which acts on the vortex line stemming from its interaction with the lattice in the inner crust, should be evaluated from the pinning energy. Numerical simulations to calculate this quantity have been performed by considering the orientation of the bcc lattice (Grill & Pizzochero 2012). The results have shown that the order of magnitude of the maximum pinning force is approximately 10 15 dyn cm −1 and that the density where the pinning force reaches its maximum is in the inner crust density realm, i.e., from the neutron drip density to half of the nuclear saturation density ρ 0 . Seveso et al. (2012) later simplified the function f pin (ρ); i.e., the pinning force linearly increases up to the maximum magnitude and then linearly decreases. They also assumed that the pinning force would vanish at the boundary of the inner crust because the lattice exists only in the crust and that in the outer crust, there are no free neutrons to produce vortices. The location of the maximum pining force crucially depends on the pairing gap. Yet, owing to the complexity of the nuclear medium polarization effects, referring to the self-energy (which affects the density of states for pairing) and the correction of the interaction vertex (which affects the pairing strength) (Schulze et al. 1996;Shen et al. 2003), the exact information of the pairing gap remains an open problem in nuclear physics. However, different calculations indicate the polarization effects correspond to a reduction of the pairing gap calculated with a realistic interaction by a factor between two and three (Lombardo et al. 2001). Therefore a reduction factor β, i.e., ∆ = ∆ 0 /β, is introduced to account for the polarization effect in Ref. (Donati et al. 2006). Where the gap ∆ 0 is calculated by using the free single particle energy with microscopic nucleon-nucleon potential (Argonne). The maximum polarization effect corresponds to β = 3 with the maximal pairing gap of about 1 MeV, while β = 1 corresponds to nonpolarized effect with the maximal pairing gap of about 3 MeV. This definition is adopted in Grill & Pizzochero (2012); Seveso et al. (2012) to evaluate the pining force. In the present calculation, we take the same f pin (ρ) in two extreme configurations, i.e., the nonpolarized effect (indicated by β = 1) and the maximum polarization effect (indicated by β = 3), as in Seveso et al. (2012). Also, the shape of f pin (ρ) is adopted from Seveso et al. (2012), whereas the magnitude of the pinning force in the both configuration is chosen to reproduce the max critical lag value ∆Ω cr max of the Vela-like pulsars.
Here, we follow Seveso et al. (2012) to evaluate the critical lag. In Figures 2 and 3, we report the pinning force F pin (x), the Magnus force F * m (x), and the criti-cal lag ∆Ω cr (x) by adopting the BHF/RBHF input in Section 3. Figure 2 demonstrates the differences in the BHF and RBHF results in the case without polarization, while in Figure 3, the results with or without the inclusion of polarization are compared in the RBHF case. We see that there is no apparent difference in the shapes of the pinning force, the Magnus force, and the critical lag for different cases, except for their maximum values and the related locations. The pinning force and the critical lag ∆Ω cr (x) for M = 1.4 M are sharper than those for M = 1.2 M since a larger stellar mass leads to a more compressed stellar structure and a smaller (x max , R ic ) region. Similarly, the softer BHF EOS results in a more compressed stellar structure and a smaller (x max , R ic ) region than the RBHF EOS. When the polarization effect is considered, the density for the maximum pining force shifts to lower densities (Grill & Pizzochero 2012;Seveso et al. 2012). Accordingly, the region (x max , R ic ) is compressed. Since the angular momentum reservoir is assumed to accumulate in the thin sheet from x max to R ic in the snowplow model, a smaller (x max , R ic ) region should correspond to less available transferred angular momentum, as will be expressed in Equation (23) in Section 4.2. In addition, because the Magnus force F * m (x) is free of polarization, F * m (x) in the cases with or without polarization is identical.

Frequency jump and fast-decaying component of the glitch
Immediately before a glitch, a lag of ∆Ω cr max will have been created at x = x max . For an axial symmetry system, the angular velocity Ω n (x) of the superfluid component of the star is in fact proportional to the number N (x) of vortices enclosed in a cylindrical region with radius x, and can be expressed as follows: One can use this expression to evaluate the number of vortices stored at the peak just before a glitch, i.e.: Moreover, this max critical lag value ∆Ω cr max can also be evaluated by the time interval between two glitches ∆t g and the spin-down rate: Due to the particular shape of ∆Ω cr (x) in Figures 2  and 3, one can assume that just before a glitch, the excess vortices in the region x > x max have been entirely removed by the Magnus force, and consequently only the N υ vortices respond to the angular momentum transfer. Following the procedure in the "snowplow" scenario (Seveso et al. 2012), we employ dL g = Ω n (x)dI n to evaluate this angular momentum transfer. One should note that only the N υ vortices at the peak of the pinning potential ought to be considered in the calculation. After performing the integration on the cylindrical region x max < x < R ic , the requested angular momentum transfer can be obtained as follows: As mentioned above, during the glitch, only a fraction of the core neutrons might be coupled to the charged component. The parameter Y gl is introduced to describe the fraction that contributes to the charged component of the star at the time of glitch (Pizzochero 2011;Seveso et al. 2012;Haskell et al. 2013;Hooker et al. 2015). Y gl is related to the trigger mechanism of glitches and cannot be determined from the current simple "bulk" model. Presently, we do not have a good estimation from the hydrodynamic simulation, and we treat it as a free parameter to be determined by reproducing the observed glitch size. Y gl = 1 means that the core superfluidity is fully coupled to the charged component during a glitch.
Defining Q as the superfluid fraction of the total moment of inertia: the observed jump in the angular velocity is thus: with the total moment of inertia I tot = I n + I p . If the recoupling of the rest of the core superfluid is responsible for the relative angular deceleration of the charged component immediately after the glitch, from angular momentum conservation ∆L g = ∆(I p Ω p ) = 0, the relative deceleration of the crust immediately after the glitch follows: whereΩ p is the steady state pre-glitch spin-down rate.
In addition to a partially coupled core, we also consider that the protons in neutron star cores possibly form a type II superconductor in which the magnetic field is carried by flux tubes; thus, the strong interaction between the flux tubes and the neutron rotational vortices could impede the vortex motion and effectively pin some of the neutron vorticity in the core (Ruderman et al. 1998;Sidery & Alpar 2009). Haskell et al. (2013) proposed a parameter ξ, which is defined as the fraction of pinned vorticity in the core by the flux tubes, to account for this effect. Namely, this ξ vorticity does not contribute to the angular momentum transfer, i.e., the total number of vortices stored at x max just before a glitch is instead: One should note that the effects of the flux tube are not included completely, especially, the effect of the flux tube on the pinning force, which can influence all stages of the glitch dynamics (Ruderman et al. 1998;Gügercinoglu & Alpar 2014Graber et al. 2018;Sourie et al. 2020), is omitted. Substituting the above expression of N υ for Equation (21), one can obtain the angular momentum exchanged during the glitch. Accordingly, the coupled fraction of superfluid Y gl can be derived by fitting the size of a glitch, ∆Ω p : The instantaneous step in the frequency derivative then follows: For the 2000 January glitch from the Vela pulsar, ∆Ω p = 2.2 × 10 −4 rad s −1 . The 2000 Vela glitch also yielded the first measurement of the fractional change in the spin-down rate on a short time scale of approximately one minute, ∆Ω p /Ω p = 18 ± 6 (1σ error) (Dodson et al. 2002).
In Figure 4, the obtained Y gl in the present study is shown as a function of the stellar mass under various choices of ξ for the BHF and RBHF models. It is seen that Y gl 15% for BHF and Y gl 30% for RBHF in the case of an unpolarized gap. A maximally polarized gap would result in a much loosely coupled core, with Y gl 4% for RBHF.
In Figure 5, we confront a consistent and microscopic neutron star description with the observed short time scale relaxation immediately after the 2000 Vela glitch. We see that consistent with previous works (Pizzochero 2011;Seveso et al. 2012;Haskell et al. 2013), the pinned core vertices are disfavoured by the observed ∆Ω p /Ω p ; hence, most of the vorticity in the core should be free. One main point here is that when the pairing gap is considered to be not polarized, one cannot fit the step in the frequency derivative at the 1σ level even in the wholly unpinned case ξ = 0. This is the case for both BHF and RBHF models. It is only barely achievable at the 1σ level if less than ∼30% of the vorticity in the core is pinned with for a maximally polarized gap in the stiffer RBHF case because the RBHF EOS allows a large stellar radius (∼12.5 km compared to ∼11.1 km with BHF; see Figure 1) and the correspondingly thicker accumulating sheet. One should realize, however, that the above results are obtained without considering the entrainment effect previously mentioned. The inclusion of the entrainment effect would further worsen the situation. For example, reducing the available angular momentum by a factor of m * n /m n ∼1.35 would make it impossible to reconcile with the observations even for RBHF, especially that the actual neutron superfluid should be more strongly entrained by the crust, with m * n /m n greater than ∼1.5 (Chamel 2005(Chamel , 2012Watanabe & Pethick 2017;Chamel 2017). The current framework also provides stringent constraints on the mass of the Vela pulsar, i.e., less than 1.16 M , as seen in Figure 5 (b), if both the observed glitch size and the post-glitch jump are to be reproduced. The inferred mass of the Vela pulsar is not only much smaller than previous values, either from the largest observed glitch in a pulsar (Pizzochero et al. 2017) or from combining X-ray with radio data (Ho et al. 2015), but also approaches that of the possible lightest measured neutron star, the Pulsar J0453+1559 companion of mass 1.174 ± 0.004 M (Martinez et al. 2015). How such light neutron stars are formed and their implications on progenitors and binary stellar evolution are presently open problems. Finally, our calculations on ∆Ω p /Ω p are based on a strong assumption that the post-glitch deceleration is caused entirely by the increase in the effective moment of inertia of the star; however, the crust-core mutual friction is suggested to play an important role . Future studies are planned by incorporating both effects to reproduce the glitch magnitudes and the post-glitch rotational evolution.

SUMMARY AND DISCUSSIONS
The glitch observation of Vela-like pulsars (especially the details on short-term morphology) may provide a test of the standard superfluid glitch theory. To do so, it is crucial to make use of consistent and microscopic descriptions of the matter state and the structure of neutron stars. In the present study, our input star properties are derived from the microscopic BHF/RBHF theory, including a consistent determination of the star EOS and composition. The two representative EOSs (without any adjusted parameter) have been tested from laboratory nuclear experiments and the latest astrophysical observations and thus are regarded as reliable at ∼0.5−3 ρ 0 times the nuclear saturation density, which is crucial to determine the structure of the crust and therefore the pinning vortices.
Within the particular EOS models we use, we found that, although the snowplow scenario of the simple twocomponent model has no problem matching the amplitudes of the glitches from the Vela pulsar, it could hardly explain the one available observation of short-term relaxation from the January 2000 glitch. Additionally, the Vela mass is unexpectedly small. These results thus challenge to some extent the present understanding of the simple two-component framework of superfluid glitch theory, and strongly support the conclusions that the mutual friction may closely relate to post-glitch recoveries. Extended studies on the post-glitch rotational evolution with a large sample of EOSs will be reported in a separate work. Further efforts should also be made to construct a unified theory able to describe on a microscopic level the complete structure of neutron stars from the outer crust to the core. Upcoming higher time resolution measurements with post-glitch relaxation should provide more insight into these issues, especially because some real-time automated glitch detection pipelines have been presented (Singha et al. 2021), which can measure the postglitch recovery phase at a high cadence.
On the other hand, there is increasing evidence of radiative changes accompanying glitch events in magnetars (Dib et al. 2008;Dib & Kaspi 2014;Kaspi & Beloborodov 2017) and radio pulsars (Weltevrede et al. 2011;Keith et al. 2013;Kou et al. 2018;Palfreyman et al. 2018). For example, there are changes in the spin-down state, emission mode changes and narrowing of the pulse associated with glitch activity in PSR B2035+36 (Kou et al. 2018). Palfreyman et al. (2018) detected sudden changes in the pulse shape coincident with the glitch event of the Vela pulsar. In anomalous Xray pulsar 1E 2259+586, a glitch was observed preceding an X-ray burst (Woods et al. 2004). One of the brightest radioquiet gamma-ray pulsars, PSR J2021+4026, showed a sudden decrease in gamma-ray emission at the glitch (Zhao et al. 2017). Recently, Feng et al. (2020) reported an association between the glitch and polarization change of the Crab pulsar. They proposed that the sudden spin-up of the crust may cause a change in the configuration of magnetic fields threaded in it and consequently in the co-rotating magnetosphere.
Therefore, although no variation could be seen either before or after the 2000 Vela glitch in profile shape, magnitude, or polarization in the radio observations on any timescale, for future studies aiming for a thorough understanding of various observed glitch behaviours, a coherent picture involving stellar interiors and their magnetosphere may be necessary. Further information will soon be obtained from high-accuracy facilities, e.g., FAST or SKA, with the help of multiband collaborative observations. We expect this information to give additional constraints to the viability of the superfluid model or to identify the glitch (together with preceding slowdown) that arises by other means, such as magne-tospheric activities, hydrodynamic instabilities and turbulence, or starquakes.