# Rotating self-gravitating Bose-Einstein condensates with a crust: A model for pulsar glitches

Akhilesh Kumar Verma<sup>1,*,†</sup> Rahul Pandit,<sup>1,†</sup> and Marc E. Brachet<sup>2,‡</sup>

¹Centre for Condensed Matter Theory, Department of Physics, Indian Institute of Science, Bangalore 560012, India

²Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université and Université de Paris, 24 Rue Lhomond, 75005 Paris, France

<figure><img src="image_1.png" /></figure>

(Received 24 May 2020; accepted 14 December 2021; published 11 January 2022)

We develop a minimal self-gravitating model for pulsar glitches by introducing a solid-crust potential in the three-dimensional Gross-Pitaevskii-Poisson equation, which we have used earlier to study gravitationally bound Bose-Einstein condensates, i.e., bosonic stars. In the absence of the crust potential, we show that, if we rotate such a bosonic star, it is threaded by vortices. We then show, via extensive direct numerical simulations, that the interaction of these vortices with the crust potential yields (a) stick-slip dynamics and (b) dynamical glitches. We demonstrate that, if enough momentum is transferred to the crust from the bosonic star, then the vortices are expelled from the star, and the crust's angular momentum $J_c$ exhibits features that may be interpreted as glitches. From the time series of $J_c$, we compute the cumulative probability distribution functions (CPDFs) of event sizes, event durations, and waiting times, which are consistent with the previous work. We show that these CPDFs have signatures of self-organized criticality, which are similar to those seen in observations of pulsar glitches and are consistent with previous work.

DOI: 10.1103/PhysRevResearch.4.013026

## I. INTRODUCTION

Rotating magnetized neutron stars [1,2], or pulsars, display glitches, which are sudden increases in their rotational frequencies. These observations have a long history [3–5], and they indicate that glitches are associated with the transfer of angular momentum, which is carried by quantum vortices in the superfluid interior, to the solid crust, in the outer layers of the pulsar. This transfer occurs because of vortex-crust interactions, as suggested in Refs. [6,7]. The quantitative modeling of pulsar glitches is complex [8], so different models have been suggested [1,9]: Some involve avalanches of superfluid vortices [10]; other mainstream models are based on neutron superfluidity in the crust [11]. General relativity effects have been considered in Ref. [12].

Neutron Cooper pairs [13], which comprise a major component of the nuclear matter in a pulsar, form a superfluid. Therefore Refs. [14,15] have proposed simply to model this superfluid by using the two-dimensional (2D) Gross-Pitaevskii equation (GPE); in addition, they have included an externally imposed potential or container and a pinning potential for the crust. Similar models have been recently studied in three dimensions (3D) [16,17]. However, gravitational effects

are important on stellar scales. It is important, therefore, to account for these crucial features in a model for pulsars and the glitches they exhibit. In this paper, we construct a minimal self-gravitating model for pulsar glitches by (a) accounting for gravitational effects via the Gross-Pitaevskii-Poisson equation (GPPE) for a self-gravitating superfluid (see, e.g., Refs. [18–20] for a nonrotating bosonic star) and (b) including rotation and an interacting solid crust.

We carry out extensive pseudospectral direct numerical simulations (DNSs) to show that our minimal self-gravitating model yields pulsar glitches with properties that are akin to those seen in observations. In particular, the time series of the angular momentum $J_c$ of the crust shows the hallmarks of self-organized criticality (SOC) [21–28]. To obtain these results, we develop a sophisticated algorithm to find the ground state of the GPPE, with rotation: This algorithm uses an ancillary advective real Ginzburg-Landau-Poisson equation (ARGLPE), an imaginary-time version of the GPPE. The resulting ground states contain vortices and yield uniformly rotating solutions of the GPPE. When we include the crust and its dynamics, we find a transfer of the angular momentum from the star to the crust, where it is dissipated by friction; if this transfer is large enough, vortices move outwards, and glitches are observed.

## II. THE MODEL

Self-gravitating GPPE superfluids are described by a complex wave function $\psi(x, t)$, governed by the following partial differential equation (PDE):

$$i\hbar\partial_t\psi = -\frac{\hbar^2}{2m}\nabla^2\psi + [V_\theta + G\Phi + g|\psi|^2]\psi, \quad (1)$$

$$\nabla^2\Phi = |\psi|^2 - \langle |\psi|^2 \rangle. \quad (2)$$

<footnote>*akhilesh@iisc.ac.in</footnote>

<footnote>†rahul@iisc.ac.in; also at Jawaharlal Nehru Centre For Advanced Scientific Research, Jakkur, Bangalore, India.</footnote>

<footnote>‡brachet@phys.ens.fr</footnote>

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.</footnote>
Here, *m* is the mass of the bosons, *n* = |ψ|² is their number density, *G* = 4π*G*<sub>N</sub>*m*² (*G*<sub>N</sub> denotes Newton's gravitational constant), and *g* = 4π*aħ²/*m*, with *a* being the s-wave scattering length [29]. We describe the dynamics of the pulsar's solid crust by a single polar angle *θ*, which evolves as follows:

$$I_c \frac{d^2 \theta}{dt^2} = \frac{1}{N} \int d^3 x \partial_{\theta} V_{\theta} |\psi|^2 - \alpha \frac{d\theta}{dt}, \quad (3)$$

$$V_{\theta}(\mathbf{r}_p) = V_0 \exp \left( - \frac{|\mathbf{r}_p| - r_{\text{crust}}}{\Delta r_{\text{crust}}} \right)^2 \tilde{V}(x_{\theta}, y_{\theta}). \quad (4)$$

$I_c$ and $V_θ$ denote the moment of inertia of the crust and the crust potential, respectively; $α$ controls the frictional slowing down of the rotation of the crust, with $\sqrt{I_c/α}$ being the crust-friction decay time; for specificity, we choose $\tilde{V}(x_θ, y_θ) = 3 + \cos(n_{\text{crust}}x_θ) + \cos(n_{\text{crust}}y_θ)$, with $x_θ = \cos(θ)x_p + \sin(θ)y_p$ and $y_θ = -\sin(θ)x_p + \cos(θ)y_p$; here, $n_{\text{crust}}$ determines the number of pinning sites in the crust potential, $r_{\text{crust}}$ is the radius at which this potential is a maximum, and $\Delta r_{\text{crust}}$ determines the thickness of the crust. We have chosen to add a friction term, with coefficient $α$, in the equation for the angle to introduce the simplest model for the spin-down mechanism that removes angular momentum from the system; of course, the interplay of this friction with the rotation leads to stick-slip events. We use periodic boundary conditions (PBCs), so, to obtain a periodic version of the angular momentum $J_z$ of the GPPE condensate, we use $π$-centered, $2π$-periodic coordinates: $\mathbf{r}_p = (x_p, y_p, z_p)$, with $x_p = -\sum_{n=1}^{n=10} \exp(-\frac{16}{100}n^2)(-1)^n \frac{\sin[n(x-\pi)]}{n}$, $y_p = -\sum_{n=1}^{n=10} \exp(-\frac{16}{100}n^2)(-1)^n \frac{\sin[n(y-\pi)]}{n}$, and $z_p = π$ ($x_p$ and $y_p$ are nearly linear near the center because they are truncated Fourier expansions of sawtooth waves).

The crust acts on the GPPE superfluid and also reacts to it [30]. If $α = 0$, the GPPE-crust coupled system (3) and (5) (or the Lagrangian from which it is derived) obeys the following global conservation laws: (a) rotational invariance, in the whole space $\mathbb{R}^3$, which leads to the conservation of the total angular momentum $J = J_c + J_z$, with $J_c = I_c \frac{d\theta}{dt}$ being the crust angular momentum and $J_z = \int d^3x \psi(\hat{e}_z \times \mathbf{r}) \cdot (-i\hbar\nabla)\psi$; and (b) time-translation invariance, by virtue of which the total energy $E_{\text{tot}} = E + E_c$ is also conserved. Here, $E_c = \frac{1}{2}(\frac{d\theta}{dt})^2$ is the rotational energy of the crust, and $E$ is the energy of the GPPE system, which we rewrite as

$$i\hbar \frac{\partial \psi}{\partial t} = -\frac{\hbar^2}{2m} \nabla^2 \psi + [V_{\theta} + (G\nabla^{-2} + g)|\psi|^2]\psi. \quad (5)$$

Equation (5) conserves the number of particles $N = \int d^3x |\psi|^2$ and (for time-independent $θ$) the GPPE energy $E = E_{kq} + E_{\text{int}} + E_G + E_V$, where $E_{kq} = \frac{\hbar^2}{2m} \int d^3x |\nabla \psi|^2$, $E_{\text{int}} = \frac{g}{2} \int d^3x (|\psi|^2)^2$, $E_G = \frac{G}{2} \int d^3x (|\psi|^2) \nabla^{-2} (|\psi|^2)$, and $E_V = \int d^3x |\psi|^2 V_{\theta}$. For $V_{\theta} = 0$, the momentum $\mathbf{P} = \frac{i\hbar}{2} \int d^3x (\psi \nabla \psi - \bar{\psi} \nabla \psi)$ is also conserved.

We solve the GPPE (5) by using a 3D Fourier pseudospectral method [31–33], with $\psi(x) = \sum_{|k|<k_{\text{max}}} \hat{\psi}_k \exp(i\mathbf{k} \cdot \mathbf{x})$ and $k_{\text{max}} = [N/3]$, where $N$ is the resolution and $[·]$ denotes the integer part. In the absence of the crust potential, friction, and rotation, we obtain the conventional GPPE (see Ref. [20], where we discuss the gravitational Jeans instability

<figure><img src="image_2.png" /><figcaption>FIG. 1. Illustrative plot of the GPPE condensate angular momentum $J_z$ vs time $t$; $J_z$ is not conserved exactly (decay of ≈1% per turn) because rotational symmetry holds only approximately if we use PBCs. Parameters: $N = 64$, $g = 5$, $G = 50$, and $\Omega = 14$; the system has four vortices. Inset: Plot of the ARGLPE-converged value (see text) of $E - \Omega J_z$ vs the rotation speed $\Omega$. The differences between scans with increasing (red) and decreasing (green) $\Omega$ indicate hysteresis; the system has no vortices along the horizontal part of the red line and four vortices along the green one.</figcaption></figure>

and the pseudospectral method (see also the Supplemental Material [34]). In all three spatial directions, our DNS uses $2π$-periodic PBCs, which we also use for $J_z$ and $V_θ$; so, even if $α = 0$, the conservation of the total angular momentum holds only approximately, because the system does not have rotational invariance in $\mathbb{R}^3$.

We first obtain uniformly rotating states for the GPPE with rotational speed $\Omega$ by solving the imaginary-time equation

$$i\hbar \partial_t \psi = -\frac{\delta}{\delta \bar{\psi}} \left( E - \Omega J_z - \mu N - \lambda \left( \frac{\mathbf{P}}{N_0 m} \right)^2 \right). \quad (6)$$

Rotational ground states are minima of $E - \Omega J_z$; $μ$, the chemical potential, and $λ$ are Lagrange multipliers; at each time step, we tune $μ$ to keep the boson number fixed; and we choose a large value of $λ$ so that $\mathbf{P}$ is small. We now obtain the following advective real Ginzburg-Landau equation (AR-

<figure><img src="image_3.png" /><figcaption>FIG. 2. Isosurface plots illustrating a rotating, compact object with vortices obtained via the ARGLPE (see text): (a) Isosurfaces of the boson density (top view) (for the spatiotemporal evolution of these isosurfaces, see Video S1 in the Supplemental Material [34]) and (b) isosurfaces of $(\nabla \times (\rho v))^2$ (side view); here, $N = 256$, $G = 800$, $g = 80$, and $\Omega = 60$.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 3. Plots of crust-potential isosurfaces in blue, with $V_{\theta} = 450$, and of ten-level isosurfaces of $(\nabla \times (\rho v))^{2}$, from our DNS of the GPPE, for the representative parameter values $V_{0} = 180$, $n_{\text{crust}} = 12$, $I_{c} = 0.01$, $r_{\text{crust}} = 1.0$, $\Delta r_{\text{crust}} = 0.15$, $\Omega_{0} = 14$, and $\alpha = 0.007$ (as in Fig. 4 below) at times (a) $t = 0.06$, (b) $t = 6.48$, (c) $t = 7.38$, and (d) $t = 9.72$. For the spatiotemporal evolution of these isosurfaces, see Video S2 in the Supplemental Material [34].</figcaption></figure>

GLP):

$$ \begin{aligned} \frac{\partial \psi}{\partial t} = & \frac{\hbar^2}{2m} \nabla^2 \psi + \mu \psi - [V_{\theta} + (G \nabla^{-2} + g) |\psi|^2] \psi \\ & - i\hbar \left( \Omega \hat{\mathbf{e}}_z \times \mathbf{r}_p - \lambda \frac{\mathbf{P}}{N_0 m} \right) \cdot \nabla \psi, \end{aligned} \quad (7) $$

which we solve to obtain the rotational ground states (minima) mentioned above; to stabilize this minimization procedure, we reset the center of mass $\mathbf{r}_{cm} = \int d^3x \mathbf{r}_p |\psi|^2 / N$ to $(\pi, \pi, 0)$, after each time step.

### III. RESULTS

We obtain rotational ($\Omega \neq 0$) states by integrating Eq. (7) until we get convergence; given our initial data, the system contains $N$ bosons. For $\Omega = 0$ and $V_{\theta} = 0$, the solution of the ARGLPE (7) converges, at large times, to ground states $\psi_{0}$ that are spherically symmetric, compact objects of radius $R$. This radius can be estimated by using a variational ansatz [35], or it can be computed numerically at both zero temperature ($T = 0$) [18,19] and finite temperature ($T > 0$) [20]. We start from the $T = 0$ state of Ref. [20], and then we increase the value of $\Omega$ in steps of 1; at each such increase in the value of $\Omega$, we use the converged ARGLPE state, from the previous value of $\Omega$, as the initial data. In the inset of Fig. 1 we plot the ARGLPE-converged values of $E - \Omega J_{z}$ versus $\Omega$: from this plot it is apparent that the nonrotating state, with no vortices, loses its stability around $\Omega \simeq 14$, to a state with four vortices. We show scans in which $\Omega$ increases (red line) and decreases (green line). The differences between these scans indicate that this system exhibits hysteresis: It has no vortices along the horizontal part of the red line and four vortices along the green one, and it goes from the four-vortex branch back to the zero-vortex branch at $\Omega \simeq 9$. The hysteresis in the inset of Fig. 1 follows from the ARGLPE minimization of $E - \Omega J_{z}$. The nucleation of new vortices under ARGLPE dynamics not only increases $E$ but also increases $J_{z}$. So, when $\Omega$ is increased, new vortices are nucleated. They are destroyed when $\Omega$ is decreased. At large values of $\Omega$, an Abrikosov lattice of vortices is present; its lattice spacing decreases with increasing $\Omega$. In Fig. 1 we give an illustrative plot of $J_{z}$ versus time $t$ for the GPPE evolution of the four-vortex state at $\Omega = 14$; we see that $J_{z}$ is not conserved perfectly, but it decreases slightly ($\simeq 1\%$ per turn). $J_{z}$ is not conserved exactly because rota-

tional symmetry holds only approximately, given our PBCs; the closer we are to the center of our computational domain, the better is this rotational symmetry. We show, in Figs. 2(a) and 2(b), isosurface plots of the boson density (top view) and of $(\nabla \times (\rho v))^{2}$ (side view), respectively, for the converged solution of ARGLPE, with $N = 256$, $G = 800$, $g = 80$, and $\Omega = 60$. Clearly, our ARGLPE-based algorithm yields rotating, self-gravitating Bose-Einstein condensates (BECs), with vortices in the GPPE; this has not been possible hitherto [36]. Previous work finds vortex tangles and reconnections [16,17], but our study does not because of its limited resolution.

We turn now to the GPPE (5) coupled with the crust-rotation equation (3). The full dynamical system, comprising the crust angle and the GPPE condensate, starts from an unperturbed ARGLPE state with the same angular velocity $\Omega = \frac{d\theta}{dt}$ for both the condensate and the crust. It spins down because the friction term (with coefficient $\alpha$) in the crust-angle dynamics slows down the crust that is coupled to the GPPE dynamics. In Figs. 3(a)–3(d), we plot the crust-potential isosurface in blue, with $V_{\theta} = 450$, and ten-level isosurfaces of $(\nabla \times (\rho v))^{2}$ for a representative set of parameters and the times $t = 0.06$ [Fig. 3(a)], $t = 6.48$ [Fig. 3(b)], $t = 7.38$ [Fig. 3(c)], and $t = 9.72$ [Fig. 3(d)]. On average, the crust gains angular momentum from the superfluid; hence, at long times, these vortices move outwards after losing enough angular momentum to the crust [compare Figs. 3(a) and 3(d)]. However, the time series of the condensate and crust angular momenta, $J_{z}$ and $J_{c}$, respectively, are complicated, and, as we show below, they display the signatures of SOC [21–28]. We illustrate this temporal evolution in Fig. 4, for the representative parameter values $V_{0} = 180$, $n_{\text{crust}} = 12$, $I_{c} = 0.01$, $r_{\text{crust}} = 1.0$, $\Delta r_{\text{crust}} = 0.15$, $\Omega_{0} = 14$, and $\alpha = 0.007$. In Fig. 4(a) we plot, versus the scaled time $t\Omega_{0}$, $(J_{c}/J_{c0} + 14)$ (blue curve), $J_{z}/J_{c0}$ (red curve), and $(J_{c} + J_{z})/J_{c0}$ (green curve), where $J_{c0}$ is the initial angular momentum of the crust. This figure shows that, if we neglect the overall, gentle decay of the total angular momentum [37], fluctuations of $J_{c}$ compensate for those in $J_{z}$. In Figs. 4(b) and 4(c) we show expanded plots of $J_{c}/J_{c0}$ for $0 < t\Omega_{0} \lesssim 170$ and $80 \le t\Omega_{0} \le 100$, respectively, to illustrate the irregular nature of the time series of the angular momentum of the crust.

From the time series $J_{c}(t)$ [Figs. 4(a)–4(c)], we see that the crust can lose (gain) angular momentum to (from) the superfluid, because vortices stick to (slip from) the crust. To characterize the statistical properties of these stick-slip events,
<figure><img src="image_5.png" /><figcaption>FIG. 4. (a) Plots, vs the scaled time $t\Omega_0$, of ($J_c/J_{c0}$ + 14) (blue curve), $J_z/J_{c0}$ (red curve), and ($J_c + J_z$)/$J_{c0}$ (green curve), where $J_{c0}$ is the initial angular momentum of the crust, for the representative parameter values given in Fig. 3 above. Expanded plots of $J_c/J_{c0}$ for (b) $0 < t\Omega_0 \lesssim 170$ and (c) $80 \le t\Omega_0 \le 100$.</figcaption></figure>

we calculate the gain $\Delta J_c$ in the crust angular momentum, between successive minima and maxima of $J_c(t)$; we call $\Delta J_c$ the event size; we scale it by $J_{c0}$. In Fig. 5(a) we present a log-log (base 10) plot of the cumulative probability distribution function (CPDF) $Q(\Delta J_c/J_{c0})$; this yields the power-law behavior $Q(\Delta J_c/J_{c0}) \sim (\Delta J_c/J_{c0})^\beta$ for the part of the CPDF that lies in the region shaded gray. Thus the probability distribution function (PDF) $P(\Delta J_c/J_{c0}) \sim (\Delta J_c/J_{c0})^{\beta-1}$; by fitting the CPDF in the gray region, we find $\beta \simeq 0.7 \pm 0.1$.

The upper and lower cutoffs of the power-law regimes in the CPDFs are chosen by identifying flat regions in our local-slope plots. As we increase $\Omega$, the number of vortices increases, as does the range of the glitch sizes and, consequently, the extent of the power-law regimes in the CPDFs. Thus the lower and upper power-law cutoffs will depend on the number of vortices in the system.

Previous papers show that the power-law exponent of the glitch size has the following values: 0.02 with moat and 0.81 without moat in Ref. [38], between 1 and 2 in Ref. [39], and 1.31 in Ref. [40]. Previous papers such as Refs. [38–40] focus on the sizes and waiting times of the events, but not their durations, which are an element in this paper. Indeed, our work should be relevant to new investigations in which glitch durations are being resolved well [41–43].

Next, we calculate the event-duration time $t_{ed}$, i.e., the time difference between successive minima and maxima of $J_c(t)$, and thence the CPDF $Q(t_{ed}\Omega_0)$ [log-log, base 10, plot in Fig. 5(b)]. This plot shows, in the shaded gray region, that $Q(t_{ed}\Omega_0) \sim (t_{ed}\Omega_0)^\gamma$, with $\gamma \simeq 2.1 \pm 0.3$. Clearly, the PDF $P(t_{ed}\Omega_0) \sim (t_{ed}\Omega_0)^{\gamma-1}$, in this region.

Finally, we compute the CPDF $Q(t_{wt}\Omega_0)$ of the waiting time $t_{wt}$, i.e., the time difference between two successive maxima. From the shaded gray region in the semilog (base 10) plot in Fig. 5(c), we observe the exponential form $Q(t_{ed}\Omega_0) \sim \exp(-t_{ed}\Omega_0)$.

These power-law behaviors of $Q(\Delta J_c)$ and $Q(t_{ed})$ and the exponential tail of $Q(t_{wt})$ together show that the stick-slip motion, between superfluid vortices and the crust, yields the time series $J_c(t)$, which has all the signatures of SOC found in measurements of pulsar glitches [25–28]. Some classes of pulsars exhibit a glitch-size PDF of the type we have obtained in our model; in particular, they show a power-law behavior in this PDF, over a certain range of sizes, and the power-law exponents lie in the range $-0.13 \lesssim -(\beta - 1) \lesssim 2.4$. For the pulsar PSR J 1825-0935 [26], the exponent for the power-law glitch-size PDF is $\simeq 0.36$, nearly the same as our calculated exponent. In many pulsars, including PSR J 1825-0935, the waiting-time PDF has an exponential tail; this is in agreement with our result [26,44].

We have shown that uniformly rotating vortex-containing gravitationally bound solutions of the GPPE can be generated by starting the evolution from initial data obtained by integrating to convergence the (imaginary-time) ARGLPE (7). We have built on this GPPE and introduced a minimal self-gravitating model, with a single, angular, dynamical variable for a solid crust coupled with a rotating GPPE star. We have demonstrated that this model exhibits stick-slip dynamics, whose statistical properties we have characterized by computing the event-size and event-duration CPDFs $Q(\Delta J_c/J_{c0})$ and $Q(t_{ed}\Omega_0)$, which show power-law forms, and the waiting-time CPDF $Q(t_{wt}\Omega_0)$, which exhibits an exponential tail. These

<figure><img src="image_6.png" /><figcaption>FIG. 5. Log-log (base 10) plots of (a) the CPDF $Q(\Delta J_c/J_{c0})$ of the event size and (b) the CPDF $Q(t_{ed}\Omega_0)$ of the event duration. (c) Semilog (base 10) plot of the CPDF $Q(t_{wt}\Omega_0)$. $J_{c0}$ and $\Omega_0$ are the initial angular momentum and angular velocity of the crust, respectively. Our DNS data are shown in blue; the black lines show fits (power law or exponential) to these data in the shaded gray regions in the plots.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

FIG. 6. Isosurface plots illustrating a rotating, compact object with vortices obtained via the ARGLPE (see text): Isosurfaces of the boson density (top view) are shown in the top left panel, and isosurfaces of $(\nabla \times (\rho v))^2$ (side view) are shown in the top right panel; here, $N = 256$, $G = 900$, $g = 80$, and $\Omega = 60$. Time series of the angular momentum and the radius of gyration are shown in the bottom left and bottom right panels, respectively.

SOC-type desiderata are in consonance with measurements on a class of pulsars [26].

## IV. DISCUSSION

As self-gravity is one of our main contributions in this paper, we study briefly the effects of changes in the value of

$G$. In the bottom panels of Fig. 6, we compare the ARGLPE-minimization evolution of the angular momentum and the radius of gyration for two values of $G$, namely, 800 and 900, with $N = 256$, $g = 80$, and $\Omega = 60$. It is apparent from this figure that increasing $G$ decreases both the radius of gyration and the angular momentum. The final ($t = 0.3$), converged value that we obtain from this ARGLPE minimization, $G = 800$, leads to the isosurface plots shown in Fig. 2; the top panels of Fig. 6 depict their counterparts for $G = 900$. By comparing these isosurface plots we conclude that the decrease in the angular momentum (at $G = 900$) reduces the number of vortices. A full study of the effects of self-gravity is left for future work.

We plan to study pulsar-glitch models that are more realistic than our minimal self-gravitating model. Examples include models with (a) a solid crust with six degrees of freedom, three rotational and three translational, instead of only one angle of rotation, or (b) a superconducting component with magnetic flux tubes. Although fundamental uncertainties about the input physics will remain and despite the fact that the dynamic range of the astrophysical parameters will still be out of computational reach, we believe that such generalizations of our minimal model will help us to better understand the types of statistical properties that are displayed by pulsar glitches in different pulsars [25–28,45–48]. We have shown that it is computationally feasible to include self-gravitation in the GPPE model for pulsar glitches. In future work it will be interesting to examine the dependence of the size, duration, and waiting-time statistics for glitches in the presence of a harmonic trap and as a function of $G$.

We thank DST (India), CSIR (India), and the Indo-French Center for Applied Mathematics (IFCAM) for their support and SERC (IISc) for computational resources.

[1] B. Haskell and A. Melatos, Models of pulsar glitches, Int. J. Mod. Phys. D 24, 1530008 (2015).

[2] A. Basu, P. Char, R. Nandi, B. Chandra Joshi, and D. Bandyopadhyay, Glitch behavior of pulsars and contribution from neutron star crust, Astrophys. J. 866, 94 (2018).

[3] V. Radhakrishnan and R. N. Manchester, Detection of a change of state in the Pulsar PSR 0833-45, Nature (London) 222, 228 (1969).

[4] P. E. Boynton, I. Groth, E. J., R. B. Partridge, and D. T. Wilkinson, Precision measurement of the frequency decay of the Crab Nebula Pulsar, NP 0532, Astrophys. J. Lett. 157, L197 (1969).

[5] R. N. Manchester, Pulsar glitches, Proc. Int. Astron. Union 13, 197 (2017).

[6] M. Ruderman, Neutron starquakes and pulsar periods, Nature (London) 223, 597 (1969).

[7] G. Baym, C. J. Pethick, D. Pines, and M. Ruderman, Spin up in neutron stars: The future of the Vela Pulsar, Nature (London) 224, 872 (1969).

[8] V. Graber, N. Andersson, and M. Hogg, Neutron stars in the laboratory, Int. J. Mod. Phys. D 26, 1730015 (2017).

[9] T. Sidery, A. Passamonti, and N. Andersson, The dynamics of pulsar glitches: contrasting phenomenology with numerical evolutions, Mon. Not. R. Astron. Soc. 405, 1061 (2010).

[10] V. Khomenko and B. Haskell, Modelling pulsar glitches: The hydrodynamics of superfluid vortex avalanches in neutron stars, Publ. Astron. Soc. Aust. 35, e020 (2018).

[11] G. Watanabe and C. J. Pethick, Superfluid Density of Neutrons in the Inner Crust of Neutron Stars: New Life for Pulsar Glitch Models, Phys. Rev. Lett. 119, 062701 (2017).

[12] M. Antonelli, A. Montoli, and P. M. Pizzochero, Effects of general relativity on glitch amplitudes and pulsar mass upper bounds, Mon. Not. R. Astron. Soc. 475, 5403 (2018).

[13] A. Migdal, Superfluidity and the moments of inertia of nuclei, Nucl. Phys. 13, 655 (1959).

[14] L. Warszawski and A. Melatos, Gross-Pitaevskii model of pulsar glitches, Mon. Not. R. Astron. Soc. 415, 1611 (2011).

[15] L. Warszawski, A. Melatos, and N. G. Berloff, Unpinning triggers for superfluid vortex avalanches, Phys. Rev. B 85, 104503 (2012).

[16] L. V. Drummond and A. Melatos, Stability of interlinked neutron vortex and proton flux tube arrays in a neutron star:
equilibrium configurations, Mon. Not. R. Astron. Soc. 472, 4851 (2017).

[17] L. V. Drummond and A. Melatos, Stability of interlinked neutron vortex and proton flux-tube arrays in a neutron star – II. Far-from-equilibrium dynamics, Mon. Not. R. Astron. Soc. 475, 910 (2017).

[18] P.-H. Chavanis, Mass-radius relation of Newtonian self-gravitating Bose-Einstein condensates with short-range interactions. I. Analytical results, Phys. Rev. D 84, 043531 (2011).

[19] P.-H. Chavanis and L. Delfini, Mass-radius relation of Newtonian self-gravitating Bose-Einstein condensates with short-range interactions. II. Numerical results, Phys. Rev. D 84, 043532 (2011).

[20] A. K. Verma, R. Pandit, and M. E. Brachet, The formation of compact objects at finite temperatures in a dark-matter-candidate self-gravitating bosonic system, Phys. Rev. Research 3, L022016 (2021).

[21] P. Bak, C. Tang, and K. Wiesenfeld, Self-Organized Criticality: An Explanation of the 1/f Noise, Phys. Rev. Lett. 59, 381 (1987).

[22] J. M. Carlson and J. S. Langer, Properties of Earthquakes Generated by Fault Dynamics, Phys. Rev. Lett. 62, 2632 (1989).

[23] H. J. Jensen, Self-Organized Criticality: Emergent Complex Behavior in Physical and Biological Systems, Cambridge Lecture Notes in Physics Vol. 10 (Cambridge University Press, New York, 1998).

[24] D. L. Turcotte, Self-organized criticality, Rep. Prog. Phys. 62, 1377 (1999).

[25] P. D. Morley and R. García-Pelayo, Scaling law for pulsar glitches, Europhys. Lett. 23, 185 (1993).

[26] A. Melatos, C. Peralta, and J. S. B. Wyithe, Avalanche dynamics of radio pulsar glitches, Astrophys. J. 672, 1103 (2008).

[27] M. J. Aschwanden, Self-Organized Criticality Systems (Open Academic, Berlin, 2013).

[28] M. J. Aschwanden, 25 years of self-organized criticality: Solar and astrophysics, Space Sci. Rev. 198, 47 (2016).

[29] Note that the subtraction of the mean density on the right-hand side (RHS) of the Poisson equation (2) can be justified either by taking into account the cosmological expansion [49,50] or by defining a Newtonian cosmological constant [51].

[30] The first term on the RHS of Eq. (3) can be obtained from the GPPE Lagrangian augmented with the crust rotational energy $\frac{l_z}{2} (\frac{\partial \psi}{\partial t})^2$; the variation of this Lagrangian with respect to $\psi$ yields Eq. (1), and the variation with respect to $\theta$ yields the first term on the RHS of Eq. (3) (see Ref. [52] for a similar procedure involving active particles in the Gross-Pitaevskii equation with Newtonian particles).

[31] G. Krstulovic and M. Brachet, Energy cascade with small-scale thermalization, counterflow metastability, and anomalous velocity of vortex rings in Fourier-truncated Gross-Pitaevskii equation, Phys. Rev. E 83, 066311 (2011).

[32] V. Shukla, M. Brachet, and R. Pandit, Turbulence in the two-dimensional Fourier-truncated Gross-Pitaevskii equation, New J. Phys. 15, 113025 (2013).

[33] D. Gottlieb and S. A. Orszag, Numerical Analysis of Spectral Methods (SIAM, Philadelphia, 1977).

[34] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevResearch.4.013026 for details of the pseu-

dospectral method and videos of the spatiotemporal evolution of boson-density and crust-potential isosurfaces.

[35] J. Eby, M. Leembruggen, L. Street, P. Suranyi, and L. C. R. Wijewardhana, Approximation methods in the study of boson stars, Phys. Rev. D 98, 123013 (2018).

[36] For example, such states were not obtained in Ref. [53], probably because pure GPPE evolution was used, without the benefit of ARGLPE-based minimization.

[37] We have noted above that, because we use $2\pi$-periodic coordinates to define $J_z$, the conservation of angular momentum is only approximate, even if $\alpha = 0$; also, in neutron stars, the crust angular momentum is only a few percent of the total angular momentum [2].

[38] J. R. Lönnborn, A. Melatos, and B. Haskell, Collective, glitch-like vortex motion in a neutron star with an annular pinning barrier, Mon. Not. R. Astron. Soc. 487, 702 (2019).

[39] A. Melatos, J. A. Douglass, and T. P. Simula, Persistent gravitational radiation from glitching pulsars, Astrophys. J. 807, 132 (2015).

[40] L. Warszawski and A. Melatos, Knock-on processes in superfluid vortex avalanches and pulsar glitch statistics, Mon. Not. R. Astron. Soc. 428, 1911 (2012).

[41] G. Ashton, P. D. Lasky, V. Graber, and J. Palfreyman, Rotational evolution of the Vela pulsar during the 2016 glitch, Nat. Astron. 3, 1143 (2019).

[42] E. Güercinoğlu and M. A. Alpar, The 2016 Vela glitch: a key to neutron star internal structure and dynamics, Mon. Not. R. Astron. Soc. 496, 2506 (2020).

[43] A. Sourie and N. Chamel, Vortex pinning in the superfluid core of neutron stars and the rise of pulsar glitches, Mon. Not. R. Astron. Soc.: Lett. 493, L98 (2020).

[44] J. R. Fuentes, C. M. Espinoza, and A. Reisenegger, Glitch time series and size distributions in eight prolific pulsars, Astron. Astrophys. 630, A115 (2019).

[45] C. M. Espinoza, A. G. Lyne, B. W. Stappers, and M. Kramer, A study of 315 glitches in the rotation of 102 pulsars, Mon. Not. R. Astron. Soc. 414, 1679 (2011).

[46] J. B. Carlin and A. Melatos, Generating quasi-periodic pulsar glitches using a state-dependent Poisson process, Mon. Not. R. Astron. Soc. 483, 4742 (2018).

[47] I. O. Eya, J. O. Urama, and A. E. Chukwude, On the distributions of pulsar glitch sizes and the inter-glitch time intervals, Res. Astron. Astrophys. 19, 089 (2019).

[48] J. R. Fuentes, C. M. Espinoza, A. Reisenegger, B. Shaw, B. W. Stappers, and A. G. Lyne, The glitch activity of neutron stars, Astron. Astrophys. 608, A131 (2017); 618, C1(E) (2018).

[49] P. J. E. Peebles, The Large-Scale Structure of the Universe, 2011 ed. (Princeton University Press, Princeton, NJ, 1980).

[50] M. Falco, S. H. Hansen, R. Wojtak, and G. A. Mamon, Why does the Jeans Swindle work? Mon. Not. R. Astron. Soc.: Lett. 431, L6 (2013).

[51] M. K.-H. Kiessling, The “Jeans swindle”: A true story—mathematically speaking, Adv. Appl. Math. 31, 132 (2003).

[52] V. Shukla, M. Brachet, and R. Pandit, Sticking transition in a minimal model for the collisions of active particles in quantum fluids, Phys. Rev. A 94, 041602(R) (2016).

[53] E. J. M. Madarassy and V. T. Toth, Evolution and dynamical properties of Bose-Einstein condensate dark matter stars, Phys. Rev. D 91, 044041 (2015).