Pulsar glitches from quantum vortex networks

Neutron stars or pulsars are very rapidly rotating compact stars with extremely high density. One of the unsolved long-standing problems of these enigmatic celestial bodies is the origin of pulsars’ glitches, i.e., the sudden rapid deceleration in the rotation speed of neutron stars. Although many glitch events have been reported, there is no consensus on the microscopic mechanism responsible for them. One of the important characterizations of the glitches is the scaling law \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P(E) \sim E^{-\alpha }$$\end{document}P(E)∼E-α of the probability distribution for a glitch with energy E. Here, we reanalyse the accumulated up-to-date observation data to obtain the exponent \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha \approx 0.88$$\end{document}α≈0.88 for the scaling law, and propose a simple microscopic model that naturally deduces this scaling law without any free parameters. Our model explains the appearance of these glitches in terms of the presence of quantum vortex networks arising at the interface of two different kinds of superfluids in the core of neutron stars; a p-wave neutron superfluid in the inner core which interfaces with the s-wave neutron superfluid in the outer core, where each integer vortex in the s-wave superfluid connects to two half-quantized vortices in the p-wave superfluid through structures called “boojums”.

Neutron stars or pulsars are very rapidly rotating compact stars with extremely high density.One of the unsolved long-standing problems of these enigmatic celestial bodies is the origin of pulsars' glitches, i.e., the sudden rapid deceleration in the rotation speed of neutron stars.Although many glitch events have been reported, there is no consensus on the microscopic mechanism responsible for them.One of the important characterizations of the glitches is the scaling law P (E) ∼ E −α of the probability distribution for a glitch with energy E. Here, we reanalyse the accumulated up-to-date observation data to obtain the exponent α ≈ 0.88 for the scaling law, and propose a simple microscopic model that naturally deduces this scaling law without any free parameters.Our model explains the appearance of these glitches in terms of the presence of quantum vortex networks arising at the interface of two different kinds of superfluids in the core of neutron stars; a p-wave neutron superfluid in the inner core which interfaces with the s-wave neutron superfluid in the outer core, where each integer vortex in the s-wave superfluid connects to two half-quantized vortices in the p-wave superfluid through structures called "boojums." Neutron stars (NSs) and in particular pulsars are compact stars with the highest known density in our universe (about one solar mass within 10 3 km 3 ) [1], thereby providing an astrophysical laboratory to study phases of matter under extraordinary conditions: not only at very high density but also under rapid rotation and extremely strong magnetic fields (see refs.[2,3] for recent reviews).The study of NSs attracts great interest from researchers in diverse fields as recently there has been observations of highly massive NSs [4,5] and gravitational waves from a binary NS merger [6].One of the most important unsolved longstanding problems in NS physics is pulsar's glitches, i.e., the abrupt deceleration of the rotation speed of NSs [7].Although many glitch events have been reported, there is no consensus on the microscopic mechanism responsible for them.However, several ideas have been theoretically proposed: starquakes [8,9] and the catastrophic unpinning of quantum superfluid vortices [10,11].However, neither of these proposals has been rigorously proven or accepted widely in the community, as explained below.
One of the important characterizations of the glitches is the scaling law of the cumulative probability distribution of glitch sizes P (E) ∼ E −α with an exponent α, giving the probability of the occurrence of a glitch with energy E [12] (see also ref. [13]).The value of α ≈ 0.14 was obtained in ref. [12].Adopting the scaling law seems to be natural if one regards the glitches occurrence as consequences of starquakes in the crust region at the surface of NSs.In fact, this scaling law resembles the Gutenberg-Richter law expressing the probability distribution of the total number of earthquakes of certain magnitude [14].This is known as an instance of power-law distributions commonly found in diverse subjects in the natural, human and social sciences; for example the Pareto distribution drawn from economics [15] and scale-free networks such as the World Wide Web in network science [16].The "network" is one of the key ingredients in our study.
The other key ingredient of our study is quantum vortices in superfluids.Superfluids are fluids possessing zero viscosity which can support states with persistent flows.In laboratory experiments, 4 He atoms become a superfluid at low temperatures.In the interiors of NSs, neutrons form Cooper pairs and exhibit superfluidity [17] (see ref. [18] for a recent review) [see Figure 1 (a)], which is consistent with the thermal evolution of NSs, the long relaxation time after each glitch [8,9,19], and the neutrino emissivity [20].When superfluids undergo rotation, vortices are created.One of remarkable properties of superfluids is the fact that these vortices carry a quantized amount of angular momentum unlike classical fluids, hence the name quantum vortices.Since neutron stars are rapidly rotating, the superfluid component is pierced by quantum vortices lying along the rotation axis, which typically form an Abrikosov triangular lattice whose number can reach 10 19 .There is also an accompanying ordinary component (not forming Cooper pairs) which lowers its rotation speed by releasing gravitational waves and electromagnetic pulses, whose period is observed to be a continuous quantity.On the other hand, the superfluid component maintains a constant rotation speed, as a consequence of superfluidity.Therefore, a gap between the rotation speeds of the two components grows over time.The superfluid component can lower its rotation speed by releasing vortices.If the vortices can be released one by one, then the superfluid component's rotation speed can catch up with the ordinary component immediately once this difference reaches the amount corresponding to one quantized vortex, and the change of the rotation should be (almost) smooth at any given time, which does not explain the appearance of glitches.
In order to overcome this problem, the hypothesis of catastrophic unpinning of quantum vortices [10,11] assumes that all vortices are pinned to nuclei, which play the role of impurities in the crust [see the inset of Figure 1 (a)], analogously to metallic superconductors where Abrikosov vortices are energetically favored to be pinned to impurities.Then, in order to explain the occurrence of glitches, an avalanche of unpinning of a large number of vortices is assumed to occur spontaneously.Several models have been suggested, but one usually needs some phenomenological parameters in order to account for the momentum transfer from the core region to the crust, a subject around which there is quite some uncertainty, fueling a long debate [21][22][23][24][25][26].
A more serious drawback comes from a recent work based on microscopic calculations [27] showing that, unlike in metallic superconductors, pinning of vortices to impurities is energetically disfavored in the case of nuclear superfluids.
Here, we first reanalyse the accumulated up-to-date observation data [28,29] to obtain the scaling law P (E) ∼ E −α with the exponent α ≈ 0.88.We then propose a simple microscopic model that naturally deduces the scaling law with this exponent without any additional free parameters.Our model explains the origin of the glitches in terms of a quantum vortex network that arises at the interfaces of the two different kinds of superfluids in the cores of NSs.In contrast to the catastrophic vortex unpinning hypothesis, we need no (un)pinning for the accumulation of quantum vortices.
The two different kinds of superfluids explained above were theoretically predicted as two different types of neutron Cooper pairs: s-wave and p-wave parings.While an s-wave paring [17] is dominant in the low-density regime relevant for the neutron star outer core, a p-wave paring is dominant in the high-density regime relevant for the inner core [30][31][32][33][34][35][36][37][38][39]. 1 Therefore, we will assume that the interior of neutron stars consists of a layer structure with a p-wave inner core surrounded by an s-wave outer core forming a spherical shell [Figure 1 (b)].
We show that the glitch mechanism can be explained by the crucial property of p-wave superfluids: the existence of halfquantized vortices (HQVs) carrying half-quantized circulations [40], as opposed to integer-quantized vortices (IQVs) [35,37,38,41,42].HQVs are energetically favored so that one IQV is split into two HQVs [denoted as blue and red vortices in Figure 1  (b)] with additional topological charges cancelling each other (see Supplementary Information).Quantum vortices in the s-wave and p-wave superfluids are connected through junctions called "boojums" [43] at the interface; one IQV in the s-wave superfluid is connected to two HQVs in the p-wave superfluid.We thus have a picture of a large number of vortices penetrating the p-wave inner core surrounded by the s-wave outer core as in Figure 1 (b).Then, it is possible that the same pair of two HQVs in the p-wave core are connected to a single IQV in the upper and lower s-wave regions [see (i) in Figure 1 (b)].If this is the case, a neutron star superfluid can change its rotation speed by releasing vortices one by one from the core as is the case without vortex pinning at the crust, which would not explain the appearance of glitches.However, more generally, as illustrated in (ii) in Figure 1 (b), we can have two HQVs connected to a single IQV in the upper (lower) s-wave region which are connected to a different IQV in the lower (upper) s-wave region, thereby leading to the formation of a vortex network composed of a cluster of connected vortices as in Figure 2. As we show below, in this case, each cluster can contain a large number of vortices, which exhibits the power-law distribution of glitches.Observed distribution of energy.-Let us analyze the glitch dataset reported by ref. [28]; we additionally use the pulsar catalogue [29] to retrieve the pulsar periods.Let P o (E) denote the cumulative probability of observed glitch energies as in [12] (see Appendix for details).Figure 3 displays the log-log plot of P o (E) obtained from the analysis of 533 glitches, where the line, determined by the least squares method using all 533 data points, shows the cumulative distribution of the power law One can see that in the central region (away from extremely small or extremely large glitches), the cumulative distribution is well approximated by the power law.
Cluster size and energy distribution.-The radius of the whole core is typically 10km.There is an uncertainty for the size of the p-wave inner core but it is thought to be a few km.The mean intervortex distance is of order 10 −6 m, and the number of vortices is of order 10 19 .We assume that the s-p interfaces are mostly flat and parallel for simplicity although the outer and inner cores would be almost spherical.However, our results do not depend on the precise details of the underlying shapes and sizes at all.We further assume a triangular vortex lattice which is rigidly rotating, as is typically the case for superfluids.We consider a rotating frame in which this lattice is static.Then, at the s-p interface, two HQVs with different topological charges in the p-wave region will pair and connect to an integer vortex in the s-wave region via a boojum.
In order to study the cluster size distribution we employ the following model.We assume that the HQVs in the p-wave region form an Abrikosov-like triangular lattice with periodic boundary conditions; however we do not assume any restriction on the blue/red vortex pattern except that for a red vortex it is always possible to find a blue nearest neighbor and vice versa.We then consider a triangular vortex lattice of rhombic shape and simulate two random blue-red pairings, that would correspond to the boojums at the top and bottom ends of the vortices.The superposition of the two pairings generates loops of various sizes, which corresponds to clusters in which all vortices are connected via boojums.This was schematically depicted in Figure 2, and an example of configurations taken from our simulation is presented in Figure 4 (note that a minimal loop made of two vortices appears in the form of a segment).Statistical analysis of the cluster size distribution shows a power-law behaviour (for details of the analysis see Appendix).The best estimate of the probability distribution of cluster size is where s is the cluster size and the subscript t stands for "theoretical"; the simulated data points and the best fit to Eq. ( 2) is displayed in Figure 5.Then, we can define the cumulative probability as P t (s) = smax s p t (u)du.A cluster of size s defines a region inside of which the number of vortices is of order s 2 .Since there is no reconnection between HQVs with different topological charges [44,45], when a cluster is expelled from the neutron star core, it necessarily drags all the other vortices enclosed by that cluster.It is therefore safe to assume that the energy associated with the emission of a vortex cluster satisfies the relation E = cs 2 for some constant c.By using this relation, we can translate the size distribution in Eq. ( 2) to the energy probability distribution, p t (E) ∼ E −1.8±0.2 and the corresponding cumulative distribution (see Appendix) (3) Thus, our model gives a description of the scaling law in Eq. ( 1) for the set of observed glitches without any free parameters.
Conclusion.-We have obtained the power-law scaling law of Eq. ( 1) for glitches from recent observational data, and proposed a simple model to explain the scaling law of Eq. ( 3) based on vortices penetrating the p-wave superfluid core surrounded by the s-wave superfluid.Boojums connecting two HQVs in the p-wave superfluid and one integer vortex in the s-wave superfluid give rise to clusters of vortices, whose distributions realise the scaling law.The appearance of the vortex network naturally explains a similarity with power-law distributions in other systems widely discussed in network science.The key strength of our approach is that our model contains no free parameters to explain the observed data.
Figure 3 displays the log-log plot of P o (E) obtained from the analysis of 533 glitches, where the line, determined by the least squares method using all 533 data points, shows the cumulative distribution of the power law in Eq. ( 1).Note that a change of the overall energy scale in Eq. ( 4) would not affect this behaviour.
This scaling law is significantly different from the previously found scaling behaviour P o (E) ∼ E −0.14 [12].The reason can be considered as follows: ref. [12] used 27 data points of glitches in the energy range 10 39 -10 42 erg, for which the upper limit came from the observational limitations at that time, while our analysis with 533 data sets covers a broader range of energies than theirs per Figure 3, which yields the observed differences of the exponents.Another important point is that several datum corresponding to extremely giant glitches with energy higher than 10 43.6 erg are rare events, with corresponding lower probability, giving a smaller contribution to our fitting.Cluster configurations are generated by simulating two uncorrelated random pairings in a triangular lattice with periodic boundary conditions.This is achieved by a version of the worm algorithm.First, we choose an arbitrary pairing, e.g. the simple pairing in Figure 6, we then perform a series of decorrelation steps to obtain random pairings.Each step consists of the following operations.First we generate a path defined by a sequence of sites as follows: we randomly choose a site, that will be the first in the sequence; then we add the site to which the first one is paired; next, we randomly choose one of the five available nearest neighbours of the latter site (if it is a paired site already, which is already in the sequence, then it cannot be chosen) and add it to the sequence, followed by the corresponding paired site, and so on, until we encounter a site that is already in the sequence.This defines a path that contains a loop.If the loop has an odd number of segments, we discard it; otherwise we flip the pairing along the loop.This completes one decorrelation step; we perform 10 3 steps for each site to obtain the first configuration (interpreted as the boojums at the top s-p-wave interface) and the same number of steps to obtain the second configuration (interpreted as the boojums at the bottom s-p interface).Superposing the two configurations gives a schematic top view of the HQV system, in which loops represent vortex clusters that are kept together by boojums.The size of the clusters is calculated by a standard analysis of the adjacency matrix, which is the accepted procedure for complex networks.See Supplementary Information (Sec.c) for possible realisations in condensed matter setups of s-p-s interfaces.In this section we statistically analyze 100 independent simulations of size with 48 2 = 2304 sites.We use a Bayesian technique based on Dirichlet distributions [46], which is appropriate for discrete random variables, such as our cluster size.The likelihood of observing a data set is given by the multinomial distribution

Statistical analysis
where n i is the number of observed clusters with size s i = 2i, N = i n i and q i are the corresponding probabilities.Since our simulations are independent, we simply consider the data from all simulations at the same time.We can infer the the probabilities q = {q i } given the observed counts n = n i using Bayes' theorem, p(q|n) = p(n|q) p(q) dq p(n|q) p(q) , where p(n|q) is the likelihood defined in Eq. ( 6) and p(q) is the prior distribution of the probabilities q.The appropriate choice of prior is the Dirichlet distribution with a i = 1, ∀i; this is chosen since we do not assume any specific knowledge about the probabilities before doing the simulations.For the properties of the Dirichlet and multinomial distribution the posterior distribution will also be of Dirichlet type, namely The expectation value of the probabilities will then be The error bars for the q i 's can be estimated by taking the square root of the variance Var[q i ] = E[(q i − qi ) 2 ] over the distribution Eq. ( 9); they turn out to be quite small, essentially invisible in Figure 7.We try to fit the data {s i , q i } with two different models, namely exponential and power law, with a least-squares method, in order to understand the general behaviour of the cluster size distribution.By looking at the Akaike information criterion (see Supplementary Information (Sec.d) for details), it is clear that that the power-law model is preferred and the best fit gives where ±0.3 denotes the statistical error.In Figure 5 in the main text and in Figure 7 the data {s i , q i } and the inferred law of Eq. ( 11) are represented both in linear and logarithmic scales, respectively.Figure 7 shows the log-log plot, in which the line is determined by the least square method using all 100 simulations.The large size configurations on the right part deviate from the line, but these deviations are quite tiny (≈ 10 −4 -10 −3 ) giving a smaller contribution to our fitting.In fact, one can find that the fitting is very good in the linear plot in Figure 5.
We have also generated simulations with smaller N and get similar results, thereby implying the N independence of our results for large N .
Translation from the cluster size distribution to the glitch energy distribution By using the relation E = cs 2 between the glitch energy E and the vortex cluster size s, we can translate the size distribution in Eq. ( 11) to the (cumurative) energy distribution of glitchs as with some constants c and c , which defines the energy probability distribution, p t (E) = c E −1.8±0.2 and the corresponding cumulative distribution with a spatial angle θ (0 ≤ θ < 2π) around the vortex core [35,37,38,41], where A 0 is a constant matrix representing the ground state configuration.On the other hand, a HQV takes the form of with the D 4 biaxial nematic ground state A 0 ∼ diag.(1,−1, 0), where the sign ± (∓) corresponds to HQVs with different topological charges cancelling each other [40], denoted by red and blue in the main text.The exponential factors e iθ and e i θ 2 for IQV and HQV are the origin of integer and half quantizations, respectively.

b. A story of the boojum
Originally the boojum is a particular variety of the fictional animal species called snarks created by Lewis Carroll in his nonsense poem "The Hunting of the Snark."The similar structures found in helium superfluids were named boojums by Mermin [43].See [70] for the story how he made "boojum" an internationally accepted scientific term.Now boojums have been predicted to occur in 3 He superfluids [71] in particular at the A-B phase boundary [72,73], liquid crystals [74], Bose-Einstein condensates [75], quantum field theory [76][77][78] and high density quark matter relevant for neutron star cores [79,80].

c. Proposals for laboratory experiments
We further point out that our mechanism for glitches by vortex networks through boojums can be simulated by laboratory experiments.One is with ultracold atomic gases, a mixture of a scalar Bose-Einstein condensate (BEC) and a spin-2 nematic BEC [81][82][83] (see ref. [84] for a review of spinor BEC), where the Gross-Pitaevskii equation for the latter is the same as the Ginzburg-Landau equation for a 3 P 2 superfluid: a scalar BEC sandwiched between two spin-2 BECs.The other is 3 He superfluids with the A-B phase boundary on which one vortex in the A-phase is connected to two vortices in the B-phase through a boojum [72], therefore an A-B-A phase configuration is relevant (a B-A-B phase configuration was experimentally realized [73]).

d. Akaike Information Criterion (AIC)
The Akaike Information Criterion (AIC) is defined as where L is the maximum of the likelihood function L (often simply called the likelihood), which is formed from the joint probability distribution of a sample of data given a set of model parameter values; it is viewed and used as a function of the parameters given by the data sample.Generally, a lower AIC indicates a better fit.

FIG. 1 .
FIG. 1. Inner structures of NSs.(a) Pinned vortices in the s-wave superfluid in the core surrounded by the crust (dark blue) in the outermost region.(b) Vortex network in the p-wave inner core (pink) surrounded by the s-wave outer core of the spherical shell (grey).(a) Vortices form a lattice and are pinned at nuclei in the crust (see the inset).(b) Vortices are not pinned at nuclei in the crust.A single integer quantum vortex (IQV) in the s-wave outer core is split into two half-quantized vortices (HQVs) in the p-wave inner core.(i) The same pair of two HQVs in the p-wave inner core are connected to a single IQV in the upper and lower s-wave regions, forming a cluster of a certain minimum size.(ii) Two HQVs connected to a single vortex in the upper (lower) s-wave region are connected to different IQVs in the lower (upper) s-wave region, forming a cluster with a larger size.For (a) and (b), the orange arrows denote that vortices are released when a NS decreases its rotation speed, in which case all vortices in the same cluster have to be released at the same time for the example presented in (b).

FIG. 2 .
FIG. 2. Schematic views of pairings of vortices at s-p-s interfaces resulting in vortex clusters.(a) top view of a vortex network, (b) top view and (c) view from diagonally above of a 3D configuration.In (a), black and grey dots denote integer vortices in the s-wave region at the top and bottom, respectively, while red and blue dots and lines denote HQVs in the p-wave region forming an Abrikosov lattice.Clusters with sizes one, three and four are shown.

FIG. 3 .
FIG.3.log P (E) vs. log E plot of observed glitches of energy E.

FIG. 4 .FIG. 5 .
FIG. 4. Snapshot of pairings at the top (a) and bottom (b) interface, and the resulting vortex clusters (c), which appears as loops in the top view.

FIG. 6 .
FIG. 6.Schematic representation of the worm algorithm: (a) starting configuration; (b) random generation of the worm, which terminates when it crosses itself; (c) only the loop section of the worm is retained; (d) pairs are flipped along the worm, giving rise to the new configuration.In (a) and (d), only pairs before and after flipping are coloured by red and blue.

FIG. 7 .
FIG.7.Plot in logarithmic scale of the simulated cluster-size probability data points {si, qi = p(si)} and the inferred power law pt(s) ∼ s −2.6 .