Persistent gravitational radiation from glitching pulsars. II. Updated scaling with vortex number

Superfluid vortices pinned to nuclear lattice sites or magnetic flux tubes in a neutron star evolve abruptly through a sequence of metastable spatial configurations, punctuated by unpinning avalanches associated with rotational glitches, as the stellar crust spins down electromagnetically. The metastable configurations are approximately but not exactly axisymmetric, causing the emission of persistent, quasimonochromatic, current quadrupole gravitational radiation. The characteristic gravitational wave strain $h_0$ as a function of the spin frequency $f$ and distance $D$ from the Earth is bounded above by $h_0 = 1.2\substack{+1.3 \\ -0.9} \times 10^{-32} (f/30\;{\rm Hz})^{2.5} (D/1\;{\rm kpc})^{-1}$, corresponding to a Poissonian spatial configuration (equal probability per unit area, i.e. zero inter-vortex repulsion), and bounded below by $h_0 = 1.8\substack{+2.0 \\ -1.5} \times 10^{-50} (f/30\;{\rm Hz})^{1.5} (D/1\;{\rm kpc})^{-1}$, corresponding to a regular array (periodic separation, i.e.\ maximum inter-vortex repulsion). N-body point vortex simulations predict an intermediate scaling, $h_0 = 7.3\substack{+7.9 \\ -5.4} \times 10^{-42} (f/30\;{\rm Hz})^{1.9} (D/1\;{\rm kpc})^{-1}$, which reflects a balance between the randomizing but spatially correlated action of superfluid vortex avalanches and the regularizing action of inter-vortex repulsion. The scaling is calibrated by conducting simulations with ${N_{\rm v}} \leq 5\times10^3$ vortices and extrapolated to the astrophysical regime ${N_{\rm v}} \sim 10^{17} (f/30\;{\rm Hz})$. The scaling is provisional, pending future computational advances to raise ${N_{\rm v}}$ and include three-dimensional effects such as vortex tension and turbulence.


INTRODUCTION
Pulsars that undergo spin-up glitches are predicted to emit gravitational waves at amplitudes that may be detectable by long-baseline interferometers like the Laser Interferometer Gravitational-wave Observatory (LIGO).The predicted signals fall into two categories: (i) bursts lasting ≲ 1 s, which occur during spin-up events due to the excitation of stellar oscillation modes (Andersson & Comer 2001;Sidery et al. 2010;Ho et al. 2020;Andersson 2021;Yim & Jones 2022) or the reorganization of the superfluid velocity field inside the star (Warszawski & Melatos 2012;De Lillo et al. 2022); or (ii) quasimonochromatic signals lasting ≲ weeks, which occur after the spin-up event due to relaxation processes such as Ekman circulation (van Eysden & Melatos 2008;Bennett et al. 2010;Singh 2017) or subsidence of a transient mountain (Yim & Jones 2020;Moragues et al. 2023), for which there is evidence in radio timing data (van Eysden & Melatos 2010).Both categories have been the target of LIGO searches (Abadie et al. 2011;Abbott et al. 2021;Prestegard 2016) and the development of search algorithms (Thrane et al. 2011;Prix et al. 2011;Thrane et al. 2015;Keitel & Ashton 2018;Miller et al. 2018;Oliver et al. 2019;Macquet et al. 2021;Lopez et al. 2022;Moragues et al. 2023).
★ E-mail: tcheunchitra@student.unimelb.eduSuperfluid vortex avalanches are a standard theoretical explanation for neutron star glitches (Anderson & Itoh 1975;Warszawski & Melatos 2011) although there are other viable mechanisms as well, such as starquakes (Middleditch et al. 2006;Chugunov & Horowitz 2010;Kerin & Melatos 2022;Morales & Horowitz 2023) and hydrodynamical instabilities (Andersson et al. 2003;Peralta et al. 2006;Glampedakis & Andersson 2009); see Haskell & Melatos (2015), Antonopoulou et al. (2022), Zhou et al. (2022), and Antonelli et al. (2022) for recent reviews.In a superfluid vortex avalanche, a large number of vortices collectively and spontaneously unpin and migrate outward, transferring angular momentum to the crust before they repin or annihilate at the crust.The vortices transition impulsively from one metastable configuration before the avalanche to another after the avalanche, i.e. the time-scale on which the vortices adjust their configuration is much shorter than the time-scale between avalanches.The configuration between avalanches does not change due to pinning (Link & Epstein 1991;Link et al. 1993;Donati & Pizzochero 2006;Avogadro et al. 2008;Seveso et al. 2016).None of the metastable configurations is the stable, ground-state configuration, i.e. the regular, triangular, Abrikosov array; they are all slightly nonaxisymmetric perturbations to the ground-state 1 .Interestingly, this means that the superfluid velocity field is nonaxisymmetric at all times, not just during or immediately before or after a glitch, which in turn means that the star always has a nonzero, time-dependent current quadrupole moment (Melatos et al. 2015;Haskell et al. 2022).The nonaxisymmetry is small, but it may be large enough to produce gravitational wave signals detectable using long integrations of data from a detector such as LIGO in the medium-term future.Calculating this current quadrupole signal is the goal of this paper.
Previous efforts to calculate the current quadrupole moment have used Gross-Pitaevskii simulations to evolve the pinned vortex array during and between unpinning events (Warszawski & Melatos 2012;Melatos et al. 2015).These calculations are limited by computational cost to ∼10 2 vortices.This is much less than the  v ∼ 10 15 (Ω/1 rad s −1 ) vortices in a typical neutron star.A natural next step is to calculate the current quadrupole moment analytically.However, this faces a fundamental challenge: there is no general analytic theory of far-from-equilibrium systems like vortex avalanches, which exhibit self-organized criticality and long-range correlations between widely separated vortices, not just nearest neighbors (Jensen 1998;Aschwanden 2011).This fundamental problem lies outside the scope of this paper.Instead, we estimate how the current quadrupole moment scales with  v for vortex distributions that approximate plausibly the realistic situation in a neutron star and can be analysed analytically.We test the results against an N-body code (Howitt et al. 2020;Howitt & Melatos 2022), which simulates avalanches containing up to ∼10 4 vortices and therefore extends the baseline over which the  v scaling can be verified beyond the Gross-Pitaevskii regime -although, not as far as  v ∼ 10 15 (Ω/1 rad s −1 ) as in a realistic neutron star.The results are compared with the hydrodynamic limit calculated by Haskell et al. (2022).
The paper is structured as follows.In Section 2, we summarize the physics of a vortex avalanche.In Section 3, we relate the gravitational wave strain to the current quadrupole moment for an arbitrary array of rectilinear vortices.In Section 4, we use the N-body simulations to calculate the current quadrupole moment and hence the gravitational wave strain as functions of  v in the regime 5 × 10 2 ≤  v ≤ 5 × 10 3 .In Section 5, we consider two idealized vortex configurations as upper and lower limits, namely a uniform Poisson point process and a regular periodic array respectively, motivated by the output of the Nbody simulations.The current quadrupole moment and gravitational wave strain are calculated analytically for both configurations as functions of  v , so that the N-body results for 5 × 10 2 ≤  v ≤ 5 × 10 3 lie between them.In Section 6, we extrapolate the N-body and analytic scalings obtained in Sections 4 and 5 to realistic values of  v for neutron stars.We comment briefly on the astrophysical implications and uncertainties, including the possibility of superfluid turbulence, in Section 7.

SUPERFLUID VORTEX AVALANCHES
The neutron superfluid inside a neutron star attempts to rotate uniformly and match the angular velocity Ω of the crust by forming vortices, each of which carries a quantum of circulation  = ℎ/2  , where   is the mass of the neutron, and ℎ is Planck's constant.The number of vortices,  v , is therefore determined by the total circulation, 2 2 * Ω, where  * is the radius of the neutron condensate challenges (Greenstein 1970;Peralta et al. 2005Peralta et al. , 2006;; (Anderson & Itoh 1975).As the crust decelerates, a fictitious Magnus force pulls vortices radially outward.However, vortices are embedded in a nuclear lattice and an array of magnetic flux tubes, to which they pin (Srinivasan et al. 1990;Pethick & Smith 2001;Seveso et al. 2016;Drummond & Melatos 2017, 2018;Thong et al. 2023).When the Magnus force exceeds a threshold, one or more vortices unpin at random, through single-or multi-site breakaway (Link & Epstein 1991;Link et al. 1993;Link & Levin 2022), and knock-on their neighbours through proximity or acoustic mechanisms (Warszawski et al. 2012).Knock-on triggers a runaway process, in which anywhere from one to  v vortices unpin -that is, an avalanche.The associated transfer of angular momentum to the crust produces a spin-up glitch.
Gross-Pitaevskii simulations of vortex avalanches are used to study pulsar glitches.They demonstrate key aspects of the physics, such as the operation of proximity and acoustic knock-on, and yield estimates of the probability density functions of key stochastic variables, e.g.avalanche waiting times (typically an exponential), avalanche sizes (typically a power law), radial vortex displacements (typically a Gaussian), and so on (Warszawski & Melatos 2011;Melatos et al. 2015).However, computational cost limits such simulations to  v ≲ 10 2 vortices.
Recently, Howitt et al. (2020) developed N-body simulations which show glitch-like behaviour for larger systems with  v ≤ 5 × 10 3 .The N-body simulations produce avalanches whose internal structure (e.g.wedge-like shape) and externally observable stochastic properties (e.g.waiting times and sizes) resemble qualitatively those produced by Gross-Pitaevskii simulations (Howitt et al. 2020;Howitt & Melatos 2022).Furthermore, the N-body simulations predict scalings between observables and control parameters (e.g.spin-down torque, pinning site separation, pinning potential), which agree with those produced by Gross-Pitaevskii simulations; compare, for example, Table 2 in Howitt et al. (2020) with Table 9 in Warszawski & Melatos (2011).Encouragingly, this occurs even though the N-body simulations neglect certain physics that is present in Gross-Pitaevskii simulations (e.g.acoustic knock-on) and include certain physics that is absent from Gross-Pitaevskii simulations (e.g.mutual friction).
Figure 1 depicts a snapshot of a typical vortex avalanche generated by the N-body code.The snapshot is taken after the avalanche ends and is displayed in the frame corotating with the pinning sites.Vortices that move during the avalanche are marked in red, with blue tracers tracking their discretized paths during the avalanche.Vortices that do not move before, during, or after the avalanche are drawn in grey.The avalanche is triggered by a single unpinning event at the top left of the group of red dots, slightly displaced from the others, which itself triggers the rest of the avalanche by proximity knockon.The unpinned vortices occupy a narrow "wedge", which curves clockwise to the bottom left in response to the condensate's local departure from corotation and the drag due to mutual friction, which corresponds to rotation through a dissipation angle ; see Section 2 in Howitt et al. (2020).The simulation is initialized with  v = 2000 vortices.A total of 21 vortices (out of the  v = 1838 vortices remaining before the avalanche) move during the avalanche, which is the 34th in a sequence starting from a uniform distribution.The local inhomogeneities in the distribution of grey vortices are caused by successive avalanches; many of the grey vortices participate in one of the 33 avalanches leading up to the one depicted in Fig. 1.Time-lapse movies of avalanches in action can be viewed by the interested reader in fig. 5 2020).The snapshot is displayed in the frame corotating with the pinning sites.The vortex distribution is shown after the avalanche ends.Vortices that move in the avalanche are drawn as red dots, with a blue tracer showing their discretized path throughout the avalanche.Vortices that do not move are drawn as grey dots.The simulation is initialized by drawing positions of 2×10 3 vortices from a spatially uniform distribution within container radius  * = 10.0.Pinning sites are arranged in a rectangular lattice (lattice size  pin = 0.1, site width  = 0.02) (all lengths measured in simulation units; see Appendix A).The snapshot is taken at simulation time-step  = 826 (midway between the 34 th and 35 th glitches) and the track shows movement starting from  = 798 (midway between the 33 th and 34 th glitches; time in simulation units).Control parameters are recorded in Appendix A; see also Sections 2 and 3 in Howitt et al. (2020).Vortices often annihilate at the boundary during an avalanche, although none does so in this figure.

GRAVITATIONAL RADIATION
Figure 1 makes it clear that pinned vortices are not arranged axisymmetrically between (let alone during) avalanches.The metastable, pinned configuration is clustered and features filaments interspersed with voids, even though there is no preferred axis overall.The sizes and locations of the nonaxisymmetric features are a stochastic outcome of the avalanche history.Indeed, the system exhibits many of the features of self-organized critical systems, such as long-range spatial correlations and hysteresis (Jensen 1998;Aschwanden et al. 2016Aschwanden et al. , 2018)).In this section, we calculate the gravitational wave signal expected from the nonaxisymmetric superfluid velocity field generated by the nonaxisymmetric vortex array.The calculation divides into two parts: writing the metric perturbation in terms of the current quadrupole moment (Section 3.1), and calculating the current quadrupole moment as a function of the vortex positions from the simulation output (Section 3.2).

Characteristic wave strain
In the transverse-traceless gauge, the far-field metric perturbation generated by a superposition of time-varying current multipole mo-ments   is given by (1) (Thorne 1980;Warszawski & Melatos 2012;De Lillo et al. 2022), where  is the retarded time,  is the distance from the source to the observer,  is the Newtonian constant of gravitation, and  is the speed of light in vacuum. B2,   is the (, ) angular beam pattern, a function of the observer's orientation relative to the source, defined in equation (2.30f) in Thorne (1980).In the Newtonian limit, the current multipole moments simplify to (Thorne 1980) where  *  is the complex conjugate of the scalar spherical harmonic of order (, ),  is the mass density, x is the position vector, and v is the velocity field.The integral is performed over the entire source volume.
For rectilinear vortices, the leading order contribution to equations (1) and ( 2) are  = 2,  = ±1.In cylindrical polar coordinates (, , ), the corresponding spherical harmonic is where only the -component of vorticity [∇ × (v)]  contributes to the current quadrupole due to the rectilinearity of vortex lines.Each vortex generates a vorticity field where x  = (  ,   ) is the position of each vortex, expressed in polar coordinates in the midplane  = 0, and e  is the unit vector along the rotation axis.The integral in equation ( 4) reduces therefore to a sum over all the vortices in the star, which yields In equation ( 6),  B2,2,±1

𝑗 𝑘
are Cartesian components of the tensor spherical harmonic where  denotes the inclination angle between the rotation axis of the star and the line of sight,  denotes the azimuth of the line of sight relative to the  = 0 axis, and ι and ζ denote unit vectors in the -and -directions, respectively. 2Equation ( 6) corresponds to equation ( 6) in Melatos et al. (2015).
The positions of each vortex are known from the N-body simulation, so we can evaluate equation ( 6) directly.The vortex positions in equation ( 6) refer to the inertial frame, not the corotating frame, i.e. one has   =  ,0 + Ω, where  ,0 is the azimuth of the vortex in the frame corotating with the pinning sites [e.g. in Fig. (1)].We can therefore write ℎ TT   =  B2,2,±1 ℎ 0 cos(Ω + Φ) where Φ is an astrophysically irrelevant phase (which depends on the vortex positions in the corotating frame), with and The characteristic wave strain ℎ 0 depends on the vortex configuration purely through the dimensionless quantity  in equations ( 8) and ( 9).The value of  is constant between avalanches, because the vortices are pinned, but it varies stochastically from one inter-avalanche interval to the next3 .On time-scales that are short compared to the star's spin-down time-scale, the spatial vortex configuration exhibits approximately stationary statistics, and hence so does .Moreover, Fig. 1 demonstrates that the system does not have a preferred axis, i.e.
it is isotropic and homogeneous globally albeit not locally.Consequently, the statistics of , such as its first and second moments ⟨⟩ and ⟨ 2 ⟩ respectively, depend only on  v .The prefactor ∝ Ω 2 ∝  2 v in equation ( 8) also depends on  v , of course, but the latter scaling does not depend on the spatial configuration of the vortices.The pinned, metastable vortex configuration between avalanches exemplified by Fig. 1 is a snapshot of a far-from-equilibrium, selforganized critical process.No general analytic theory exists for the statistics of such processes, due to the profound complications introduced by long-range spatiotemporal correlations, "memory", and hysteresis (Jensen 1998).In this section, we calculate the scaling of  versus  v empirically, by evaluating equation ( 9) from the output of N-body simulations such as the one that produced Fig. 1 (Howitt et al. 2020).The simulations are restricted to  v ≤ 5 × 10 3 due to computational cost.

N-body simulations
The N-body solver developed by Howitt et al. (2020) integrates the equations of motion for  v point vortices in two dimensions.The velocity dx  ()/d of a vortex at x  () equals the bulk velocity of the inviscid superfluid condensate at x  () induced by the other vortices (the N-body term) plus motion arising from interactions with a grid of pinning potentials, the container wall (via image vortices), and a viscous superfluid component [drag and mutual friction, implemented by rotating through a dissipation angle (Schwarz 1985)].The simulations are conducted in the frame that corotates with the stellar crust, which corotates with the pinning grid and viscous component.The angular frequency Ω of the crust is adjusted at every time step, in response to an external spin-down torque (astrophysically the magnetic dipole braking torque, denoted by  ext ) and the decrement in the total angular momentum of the superfluid condensate, if one or more vortices change position during the time-step (as happens during an avalanche).The N-body equations of motion and simulation parameters are summarized briefly in Appendix A for the reader's convenience.
For the -versus- v studies in this section, we run eight independent simulations, each starting with  v = 5 × 10 3 .As avalanches occur, and vortices progressively exit the system, we calculate  across a range of  v values, e.g. for 5 × 10 2 ≤  v ≤ 5 × 10 3 , as  v decreases with time4 .The simulations are initialized by drawing {x  } from a spatially uniform probability distribution.They are evolved with  ext = 0 until all vortices are pinned, as described in Section 5.1 in Howitt et al. (2020).Subsequently, the simulations are evolved with  ext ≠ 0 in order to trigger avalanches.Avalanches are detected automatically using the glitch-finding algorithm described in Section 3.3 in Howitt et al. (2020), after Ω() is smoothed over a few adjacent time-steps with a top-hat function (details about the smoothing are recorded in Appendix A).A value of  representative of each interglitch interval is calculated by evaluating equation ( 9) using the vortex positions at the midpoint of the interval.

Probability density function of 𝑄
Figure 2 displays snapshots of  versus  v as a scatter plot.The snapshots are taken at instants within the range 5×10 2 ≤  v ≤ 5×10 3 during the eight independent simulations discussed in Section 4.1.The instants correspond to the immediate aftermath of every 10-th glitch, starting from the 10-th.That is, the time series () and  v () are decimated to ensure that the plotted data points are reasonably independent statistically; 10 avalanches occur between one data point and the next, during which the vortex configuration "scrambles" to a degree (see footnote 3).Across the eight runs,  varies by a factor of ∼20 (smallest to largest) at any fixed  v value in the plotted range.One may speculate that the initial vortex configuration affects the magnitude of  for a given run, so that some runs have consistently higher , and others consistently lower .We inspect visually the  of all eight runs, and find that such behaviour is not present, suggesting that  loses memory of initial conditions.For clarity, we join the dots for only two out of eight runs in Fig. 2 as the red and blue lines.The spread of  decreases with  v , as quantified below.There is no strong trend between the central tendency of  and  v , but the broad spread in  makes it hard to identify a trend visually.
What is the probability density function (PDF) of  as a function of  v ?We develop an approximate answer to this question by noting from equation ( 9) that  is the length of a random vector, which can be calculated for each vortex configuration.Clearly, each component of S is distributed symmetrically around zero.The exact PDF of each component is challenging to predict analytically, for the reasons discussed in the second paragraph of Section 4.However, it is unimodal with a steep tail, so we approximate each component of S as a zero-mean Gaussian with standard deviation  =  v  , a power-law parameterization which is deliberately scale-invariant and motivated by the analytically tractable configuration discussed in Section 5.1.Consequently, at fixed  v ,  follows a Rayleigh distribution We sample the posterior distribution of  and  assuming equation ( 11) and uninformative priors using the Markov Chain Monte Carlo (MCMC) sampler in the emcee Python package (Foreman-Mackey et al. 2013).The resulting posterior probability distribution is presented as a corner plot in the middle three panels of Fig. 2 5 .We estimate  = 35 +20 −13 and  = −0.15+0.06 −0.06 where central values correspond to the median, and the error bars delineate the 90 per cent credible intervals of the marginalized posterior distribution.
The dark grey line and the grey band in the top panel of Fig. 2 indicate the median and the 90 per cent credible interval of the posterior predictive distribution at each  v , respectively.This posterior predictive distribution includes both the scatter due to avalanche stochasticity and uncertainty in parameter estimates.
A natural question is to ask whether the model in equation ( 11) 5 The corner plots are produced by the corner.pyPython package (Foreman-Mackey 2016).describes the data well.Ideally, we would perform a large number (≳10 3 , say) of simulations at each  v , to compare an empirical distribution of  at a given  v with equation ( 11) directly.However, this is computationally prohibitive.Instead, in the bottom panel of Fig. 2, we perform the following quantitative consistency check, using data at all values of  v simultaneously.We sample 10 3 values of (, ) from the posterior probability distribution.For each sample, and each datum in the top panel of Fig. 2, we compute /, with  =   v .We compare a histogram of all such values with a Rayleigh distribution described by (/) = (/) exp − 2 /(2 2 ) in the bottom panel of Fig. 2 (grey bars and solid black curve respectively).If the distribution of  at each  v deviates significantly from equation ( 11), the empirical histogram of / would fail to follow the Rayleigh curve.By visual inspection, the histogram and theoretical distribution are broadly consistent, giving us some confidence that equation ( 11) describes the data well.
In an ensemble of vortex configurations described by equation ( 11) with median estimates of  and , a typical  is given by where the central value and error bars correspond to the median and the central 90 percentile of equation ( 11), respectively.It is important to note that a posterior distribution is not a point estimate, but we derive equation ( 12) using the median of the posterior distribution to quantify the scatter due to stochasticity of the avalanches.This scatter is distinct from the credible interval for the estimates of  and  (vertical dashed lines in the middle three panels of Fig. 2).At a given  v , there are many possible configurations of vortices depending on the random initial configuration and avalanche history.Individual configurations typically fall within the error bars of equation ( 12), and there is no way to know where the actual configuration for a specific astronomical object lies within the range.We emphasize that equation ( 12) is not the posterior predictive distribution displayed as the grey band in the top panel of Figure 2. Rather, it can be thought of as a best-fit distribution of equation ( 11) to the data, which is quoted in order to distinguish the uncertainty in the estimates of  and  from the inherent stochasticity and hence dispersion of the avalanches which generate the vortex configurations.

𝑄 VERSUS 𝑁
The N-body simulations in Section 4 are restricted to  v ≤ 5 × 10 3 by computational cost, whereas a neutron star contains ∼10 15 (Ω/1 rad s −1 ) vortices.The -versus- v data in Fig. 2 are approximated well by equations ( 11) and ( 12) over one decade in  v .This empirical fact offers some encouragement.Needless to say, extrapolating over 12 decades of  v into the neutron star regime is dangerous even for scale-invariant, power-law expressions like equation ( 12).One would prefer to have a theoretical justification for the empirical scaling ( 12), but such a justification is hard to derive from first principles, due to the well-documented and unsolved theoretical challenges posed by the statistical mechanics of far-from-equilibrium systems, as discussed in Section 2 (Jensen 1998).Instead, in this section, we calculate  as a function of  v for two analytically tractable point processes which "bracket" the self-organized critical dynamics observed in the simulated vortex avalanches in Section 4, i.e. which bound the simulated  above and below.The uniform Poisson point process (Section 5.1), where vortex positions are drawn independently from a uniform probability distribution, approximates the regime where vortices do not repel each other, and vortex positions are uncorrelated.The uniform Poisson configuration is less symmetric than the avalanche configuration in Section 4 and produces median  greater than the central value in equation ( 12).The regular array configuration (Section 5.2) realizes the other extreme, where vortex repulsion dominates pinning, and vortex positions are highly correlated; in fact, the vortices are arranged periodically in an Abrikosov-like array.The regular array configuration is more symmetric than the avalanche configuration in Section 4 and produces median  less than the central value in equation ( 12).The analytic -versus- v scalings in Sections 5.1 and 5.2 bracket the avalanche-driven N-body scaling and can therefore be employed to place bounds on ℎ 0 in astrophysical applications, where conservatism is preferred, and the extrapolation of equation ( 12) is deemed unreliable 6 .

Uniform Poisson configuration: upper bound on 𝑄
A uniform Poisson point process generates a configuration of points in a two-dimensional region by placing each point at an independent location, with all locations having equal probability per unit area.In this paper, where vortices are assumed to be rectilinear, the region of interest is an equatorial disk of radius  * .Of course, this is not a realistic generative model for vortex configurations such as in Fig. 1.Vortex positions resulting from avalanches are correlated due to vortex-vortex repulsion and knock-on unpinning, whereas the uniform Poisson point process involves zero vortex-vortex repulsion and zero correlation between vortex positions.For these reasons, a random realization of a uniform Poisson point process is less symmetric typically, with higher , than a realistic avalanche snapshot like Fig. 1.One can calculate the PDF and moments of  analytically for arbitrary  v and derive an upper bound on  for astrophysical applications.
Consider the random vector S defined in equation ( 10).We show in Appendix B using the central limit theorem that the uniform Poisson point process leads to each component of S being distributed according to a normal distribution with standard deviation ( v /40) 1/2 , viz.

𝑝(𝑆
and similarly for (  ), where  and  denote Cartesian coordinates.Therefore, the PDF for  = |S| is the Rayleigh distribution which matches equation ( 11) with  = (1/40) 1/2 and  = 1/2.The median and central 90 percentile are given by Let us verify the above scaling numerically.The top panel of Fig. 3 shows the -versus- v scaling for a randomly realized set of vortex configurations generated by the uniform Poisson point process.Configurations are generated for 10 3 samples of  v drawn from a discrete uniform distribution on the interval [5 × 10 2 , 5 × 10 3 ], with  calculated via equation ( 9).The dark grey line and grey band 6 The least symmetric vortex configuration with the highest  corresponds to placing every vortex at one off-axis point at radius .This artificial configuration has  = (/ * ) (1 −  2 / 2 * ) 3/2  v from equation ( 9) but is not relevant to a neutron star.9) for each configuration.The dark grey curve correspond to the median of the theoretical PDF in equation ( 14) at each  v .The grey band shows the central 90 percentile of the same PDF.
Bottom panel: Samples of  for 10 3 realizations of a regular periodic array with randomized centre for 5 × 10 2 ≤  v ≤ 5 × 10 3 .For each configuration, the intervortex spacing  is drawn from a continuous uniform distribution on the interval [2 × 10 −2 , 9 × 10 −2 ],   and   are drawn from a uniform distribution on [0, ), and  is calculated by evaluating equation ( 9) for each configuration.
in the top panel of Fig. 3 displays the median and the central 90 percentile of equation ( 14), respectively.Approximately 90 per cent of the data lie within the grey band.As a quantitative consistency check, we use an MCMC sampler in the emcee Python package to sample the posterior distribution of  and  assuming equation ( 11) and uninformative priors.We find the 90 per cent credible intervals of the marginalized posterior distribution of  and  to be [0.13,0.25] and [−0.53, −0.44] respectively.These credible intervals contain the theoretical values  ≈ 0.15 and  = 0.5, validating equation ( 14).

Regular array: lower bound on 𝑄
The uniform Poisson point process produces zero correlation between vortex positions.It cannot be responsible for spatially correlated configurations containing filaments and voids, such as the one in Fig. 1.Long-range spatial correlations are a standard feature in self-organized critical systems such as sand piles (Jensen 1998).In the application studied here, they arise from vortex-vortex repulsion and the randomizing action of avalanches, which lead to capacitive depletion zones as discussed by previous authors (Cheng et al. 1988;Alpar et al. 1996;Alpar & Baykal 2006;Melatos & Warszawski 2009).The uniform Poisson process in Section 5.1 leads to a higher value of  than for spatially correlated configurations, such as the one in Fig. 1.One can place a lower bound on  for spatially correlated vortex positions by considering regular configurations dominated by vortex-vortex repulsion.In the latter regime, the vortices form a periodic array akin to an Abrikosov vortex array.The periodicity causes near-cancellation between vortices reflected about the origin.The cancellation becomes more perfect as  v increases.
In Appendix C, we derive an exact, closed-form expression for  as a function of  v for a regular square array with intervortex spacing  and centre (  ,   ) relative to the origin.The square array is not necessarily the same as the true equilibrium configuration in the complicated environment of the neutron star, but it is a reasonable and analytically tractable approximation.We note that the vortex array need not match the symmetry of pinning sites in general, as demonstrated by the vortex configurations in Section 4 and Howitt et al. (2020).In Appendix C, we find where   () is the spherical Bessel function of order .Unlike the uniform Poisson configuration of Section 5.1, a regular array configuration is not stochastic; given (  ,   ),  can be computed deterministically from equation ( 16).
It is clear from equation ( 16) that  vanishes for special choices of the centre, such as (  ,   ) = (0, 0).This follows intuitively from equation ( 9).If the regular array is centred at the origin, every vortex at (  ,  ,0 ) has a counterpart at (  ,   + ).The terms corresponding to these two vortices cancel perfectly in the sums in equation ( 9).Such perfect cancellation relies on the exact periodicity of the configuration, and would be unlikely to occur in a neutron star.In general, one has 0 <   < , where  ∈ {, }, the cancellation is imperfect, and one finds  ≠ 0. We show in Appendix C that one obtains  ∝  v −1/2 for  v ≫ 1. Notably, the scaling applies to a regular array of any shape, not just a square array.
Given  * and ,  ranges from zero to a maximum value depending on (  ,   ).In this paper, we explore the range of  by drawing   and   from a uniform distribution on the interval [0, ).The bottom panel of Fig. 3 shows the -versus- v scaling for a set of regular array configurations, with randomized (  ,   ).Vortex configurations are generated for a set of 10 3 samples of  drawn from a continuous uniform distribution on the interval [2 × 10 −2 , 9 × 10 −2 ].For each vortex configuration,  v is counted, and  is computed using equation (9).To match the range of  v in Sections 4 and 5.1, only vortex configurations with 5 × 10 2 ≤  v ≤ 5 × 10 3 are included in the plot.
Since the scatter of  in the bottom panel of Fig. 3 is a result of randomizing (  ,   ), not the stochasticity of the vortex configurations, the distribution of  at each  v does not follow equation ( 11).Deriving the analytic PDF that characterizes this scatter is difficult, due to the complicated form of equation ( 16).Instead, we calculate numerically a set of heuristic scalings comparable to the typical  of equations ( 12) and ( 15).Motivated by the analysis in Appendix C, we consider a family of curves  =  v −1/2 , and find numerically the values of  that lead to curves which lie above 5 per cent, 50 per cent, and 95 per cent of the data in the bottom panel of Fig. 3.We find empirically that the typical  for a regular array configuration is where, by construction, the central value bisects the data, and the error bars bracket the central 90 percentile.We plot the central value and the uncertainty in equation ( 17) in Fig. 3 as the dark grey line, and the grey band respectively.

Pinning potential
The -versus- v bounds in Sections 5.1 and 5.2 are calculated via geometric arguments and therefore do not depend on the pinning potential  0 (see equation (A4) and Table A1), as long as  0 occupies the astrophysically relevant regime, where vortex avalanches (and hence rotational glitches) occur.If  0 is too low, then vortices are weakly pinned.As the crust decelerates, the Magnus force pulls vortices outward one-by-one, and few vortices in the configuration are on the verge of unpinning at any instant.Hence, vortices unpinned by knock-on move slowly and repin before they knock-on other vortices, and avalanches are curtailed.Angular momentum transfers gradually to the crust, the superfluid closely matches the angular velocity of the crust at all times, and there is no glitch.If  0 is too high, then vortices cannot unpin until the angular velocity differential becomes large, and one large avalanche occurs7 .
It is interesting to ask how the value of  0 in normalized code units compares with the values of  0 in unnormalized, physical units in a realistic neutron star.It is clear from the outset that one should not expect agreement, because the length-and time-scales in the simulation are far from the neutron star regime for computational tractability.Specifically, we have (Link & Epstein 1991, see also equation ( 43) in Antonelli & Haskell 2020) where  p is the interaction energy per nucleus, and  pin is the pinning lattice spacing.The simulations in Section 4 have  0 = 2 in normalized code units, which corresponds to  0 = / in physical units and hence This is two orders of magnitude larger than the typical predictions ( p ∼ few MeV) of quantum simulations (Link & Epstein 1991;Donati & Pizzochero 2006;Wlazłowski et al. 2016), and the values inferred from glitch activity (Gügercinoğlu et al. 2022;Melatos & Millhouse 2023).This confirms that the simulations are outside the neutron star regime, as expected.Note that  0 = 2 is close to the minimum pinning potential where avalanches still occur in the Nbody solver (Howitt et al. 2020).Seveso et al. (2016) argued that the effective pinning strength for long vortex lines (∼10 3  pin ) is less than typical estimates due to the varying orientation of the pinning lattice along the vortex line.They proposed the alternative scaling i.e. equation ( 44) in Antonelli & Haskell (2020), where 10 −3 ≲  ≲ 10 −1 is a dimensionless factor quantifying the reduction in pinning force.Upon replacing equation ( 18) with equation ( 20), we find that  p in equation ( 19) is increased by one to three orders-ofmagnitude.Therefore, if the arguments of Seveso et al. (2016) apply astrophysically, then the simulations move further from the neutron star regime.

CHARACTERISTIC WAVE STRAIN
There is sufficient uncertainty about vortex microphysics in neutron stars that it is worthwhile to evaluate the gravitational wave amplitude for all three -versus- v scalings in Sections 4 and 5, that is, equation ( 12) for extrapolated numerical simulations of vortex avalanches, equation ( 15) for a uniform Poisson configuration, and equation ( 17) for a regular array.It is likely that the avalanche ℎ 0 lies between the uniform Poisson and regular array results.The wide range is a fair reflection of the uncertainties surrounding the subtle far-fromequilibrium physics of vortex avalanches (Jensen 1998).
Let us begin by considering the extrapolated simulation scaling.Upon substituting  from equation ( 12) into equation ( 8) and applying the Feynman condition The central value in equation ( 21) refers to the median ℎ 0 over random realizations of the simulated vortex configuration, and the error bars correspond to the central 90 percentile.In an astronomical context, one observes a single realization (namely the actual one), so the median is not observable directly.Instead, the characteristic wave strain for a given interglitch interval is sampled from a Rayleigh distribution with median given by the central value of equation ( 21).Equation ( 21) is ∼10 orders of magnitude lower than the central-limittheorem estimate calculated in Section 5 of Melatos et al. (2015), which assumes ⟨ℎ 0 ⟩ ∝ Ω 2  −1/2 v .The upper bound corresponding to a uniform Poisson configuration (Section 5.1) is calculated by substituting  from equation (15) into equation ( 8) to obtain The lower bound corresponding to a regular array configuration (Section 5.2) is calculated by substituting  from equation (17) into equation ( 8), to obtain

CONCLUSIONS
Neutron star glitches may be caused by the sudden unpinning and collective movement of vortices in the superfluid condensate inside the star, also known as vortex avalanches.The metastable vortex configuration between avalanches is determined by the far-from-equilibrium avalanche dynamics and is nonaxisymmetric in general, producing a small but nonzero current quadrupole moment which generates gravitational waves as it rotates (Melatos et al. 2015).In this paper, we use the N-body solver developed by Howitt et al. (2020) to simulate vortex avalanches with  v ≲ 5 × 10 3 vortices.We find that the current quadrupole moment scales ∝  v −0.15 , implying ℎ 0 = 7.3 +7.9 −5.4 × 10 −42 (/1.4⊙ ) (  /30 Hz) 1.9 ( * /10 km) 0.7 (/1 kpc) −1 upon extrapolating to the neutron star regime, with  v ∼ 10 17 (  /30 Hz).Cautious about extrapolating over ≳10 decades in  v , we also develop for safety two analytic scalings for the uniform Poisson and the regular array configurations, which represent upper and lower bounds on ℎ 0 given by equations ( 22) and ( 23) respectively.The bounds correspond to vortex pinning dominating mutual repulsion (uniform Poisson) and vice versa (regular array) and bracket the empirical scaling extrapolated from the simulations.The wide range is a fair reflection of the uncertainties surrounding the subtle question of long-range correlations between vortex positions in self-organized critical systems like vortex avalanches, which remains a fundamental unsolved problem in statistical mechanics.Due to the stochasticity of vortex avalanches, the current quadrupole moment at a fixed  v scatters around an average according to a distribution that is wellapproximated by the Rayleigh distribution [equation ( 11)].
The mechanism in this paper predicts the wave strain to increase strongly with pulsar frequency.Consequently, millisecond pulsars, which have  ≳ 100 Hz, are favored as search candidates.For a millisecond pulsar with  = 500 Hz, equation ( 21) predicts median wave strain ℎ 0 = 1.3 × 10 −39 .If vortex dynamics are dominated by pinning, then the appropriate scaling is equation ( 22), which predicts median wave strain ℎ 0 = 1.4 × 10 −29 .Recent narrow-band searches using data from the LIGO-Virgo third observing run place an upper limit for continuous gravitational waves from known pulsars at ℎ 0 ≲ 2 × 10 −26 (Abbott et al. 2022b).However, it is unclear if vortex avalanches occur in most millisecond pulsars.Only two glitches have ever been observed in millisecond pulsars, one in PSR J1824-2452 (Cognard & Backer 2004) and the other in PSR J0613-0200 (McKee et al. 2016).If the paucity of glitches in millisecond pulsars is because some unknown physics intervenes to prevent vortex avalanches (e.g.due to a temperature or age threshold excluded from the N-body simulations in Section 4.1), then the prospects for detection diminish accordingly.On the other hand, it may be that most millisecond pulsars do glitch, but only a few glitches have been detected because the average waiting time between events is longer than in ordinary pulsars, e.g.once every 10 3 years instead of once per year (Shemar & Lyne 1996;Lyne et al. 2000;Espinoza et al. 2011;Millhouse et al. 2022).In the latter scenario, the avalanche process is ongoing (albeit slowly), and millisecond pulsars may have nonaxisymmetric vortex distributions like in Fig. 1 today, while they are being observed.This scenario is explored in Section 4.2 of Melatos et al. (2015).If it is viable, millisecond pulsars are plausible targets for continuous gravitational wave searches (Abbott et al. 2010;Vargas & Melatos 2022;Abbott et al. 2022a,c).
The two-dimensional simulations and analytic calculations in this paper assume that vortices are rectilinear.In reality, the vortices in a neutron star are likely to be curved, whereupon vortex tension plays a vital role (Link & Epstein 1991;Hirasawa & Shibazaki 2001).Furthermore, the vortices in a neutron star are likely to be tangled, both because differential rotation drives macroscopic, Kolmogorovlike turbulence (Greenstein 1970;Peralta et al. 2005Peralta et al. , 2006;;Peralta & Melatos 2009;Melatos & Peralta 2010;Khomenko et al. 2019), and because pinning and differential rotation combine to generate a microscopic vortex tangle via Kelvin-wave and other instabilities (Glaberson et al. 1974;Donnelly 1991;Andersson et al. 2007;Mongiovì et al. 2017;Drummond & Melatos 2017, 2018;Haskell et al. 2020;Thong et al. 2023;Levin & Link 2023).It is unclear whether turbulence increases or decreases ℎ 0 .On the one hand, it may be argued that the randomizing and hence homogenizing action of eddy-like motions overwhelms mutual vortex repulsion and pushes the system towards a uniform Poisson configuration, increasing ℎ 0 relative to equation ( 21).On the other hand, it may be argued that eddy-like motions imprint macroscopic, eddy-related lengthscales on the vortex configuration, which are unrelated to and much larger than the microscopic vortex separation.The eddies are axisymmetric when averaged over time but not instantaneously, so they contribute to the current quadrupole moment and indeed introduce another time-scale, the eddy turnover time-scale, which competes with  −1 .Previous studies of gravitational waves from neutron star turbulence are purely hydrodynamic (Melatos & Peralta 2010;Lasky et al. 2013); they analyze the eddies only, not the quantized vortices that the eddies contain.The combined problem is subtle, requires more sophisticated and expensive simulations to be studied with confidence, and lies outside the scope of this paper.
tex  and the image of vortex , and   ,im = |x  ,im |8 .The third terms take into account the rotating reference frame, rotating with angular velocity Ω, and comoving with the container.The fourth term incorporates pinning, with  (x  − x  ) denoting the pinning potential at x  due to a pinning site at x  . pin denotes the total number of pinning sites.The pinning sites corotate with the container, so x  is fixed.Pinning sites are arranged in a square lattice centred at the origin, and the pinning potential is an isotropic Gaussian, where  is the distance from the pinning site to the vortex, and  is the characteristic width of the pinning site (Howitt et al. 2020).
The rotation matrix R  in equation (A1) embeds the effect of the dissipation for viscous superfluids by rotating the vortex velocity by an angle , as in Schwarz (1985); see also Section 2.6 of Howitt et al. (2020).
The superfluid communicates with the container through the angular momentum equation where  rel =  s / c is the ratio between the superfluid moment of inertia  s and container moment of inertia  c , and  ext is the external spin-down torque divided by  c .Ω s is the angular velocity of the superfluid, computed from the superfluid angular momentum  s =  s Ω s , which in turn is related to the vortex configuration through The equations of motion (A1) -(A5) are solved numerically using the Runge-Kutta Cash-Karp (RKCK) scheme (Press et al. 1992).
The units of the simulations are set such that  = 1,  * = 10.0, and  rel = 1.0.The spacing between pinning sites is  pin = 0.1, the characteristic width of the pinning sites is  = 0.02, and the pinning strength is  0 = 2.To accelerate the code, we only consider the closest pinning site to a given vortex in equations (A2) and (A3).The dissipation angle is set to be  = 0.1.
The simulations are initialized by drawing positions of 5 × 10 3 vortices from a spatially uniform probability distribution and evolving without spindown, until all the vortices are pinned.During this evolution, ≲ 60 vortices leave the system.Then, we turn on the external spin-down torque  ext = −2.5 × 10 −2 in simulation units.The time-step is  = 2 × 10 −3 in simulation units.All simulations run for 2 × 10 6 time-steps or until the angular frequency of the container reaches zero, whichever comes first.For each run, x 1 , . . ., x  v and Ω are recorded every 50 time-steps to reduce the data volume.Table A1 includes a summary of control parameters for the simulation.For more details, see Sections 2 and 3 and Table 1 in Howitt et al. (2020).
The angular frequency time-series is smoothed with a top-hat function of width 500 time-steps to reduce fluctuations from vortex "jiggling" (Warszawski & Melatos 2011).Glitches are found with an algorithm which scans the smoothed Ω time-series for time-steps, when the incremental change in Ω switches sign.with x = (, ) symbolizing cylindrical coordinates as in Section 3.2.The integral is performed over the disc of radius  * centred at the origin.The real and imaginary parts of S correspond to the and -components of S, respectively.One has  = | S| for a regular array defined by equation (C1).Equation (C1) can be rewritten in terms of its reciprocal lattice (or equivalently, a Fourier series), viz.
where k  1  2 =  1 b 1 +  2 b 2 is the reciprocal lattice position vector, b 1 and b 2 are the reciprocal primitive translation vectors, and  1 and  2 are integers.The Fourier amplitudes   1  2 are given by where   = |d  | for  ∈ {1, 2}.Substituting equation (C1) into equation (C4), we obtain Substituting equation (C5) into equation (C2) gives The exponential term can be expanded using the Jacobi-Anger identity, where we write k  1  2 =   1  2 (cos   1  2 , sin   1  2 ), and   () denotes the -th order Bessel function of the first kind.Substituting equation (C8) into equation (C6) and integrating with respect to , we find that every term but  = 1 vanishes.The result is and hence where   () is the spherical Bessel function of order .Evaluating equation (C10) for a square array gives equation ( 16).
The sum in equation (C10) vanishes for c = 0, because the summand is an odd function of  1 and  2 .This corresponds to no emission of gravitational waves.Physically, the wave strain vanishes because c = 0 corresponds to a configuration in which a vortex at ( 1 ,  1 ) has a counterpart opposite the centre of the star at ( 1 ,  +  1 ).These vortex pairs cancel perfectly.Needless to say, perfect cancellation (e.g.c = 0) is unlikely in a real neutron star, where one expects 0 < |c| < ( 2 1 +  2 2 ) 1/2 in general.We now calculate the scaling of -versus- v .Let us take  1 ∼  2 ∼  and hence If the number of vortices is large, with  v ≈  2 * / 2 ≫ 1, and hence   1  2  * ≫ 1 for all ( 1 ,  2 ), the spherical Bessel function is approximated by Hankel's expansion, whose leading term is Upon substituting equation (C12) into (C10), we find The summand of equation (C13) depends on  v only through the factor  * / ∼  v 1/2 .The dependence lies in the argument of the cosine, which does not contribute to the magnitude of S. The summand decreases with ( 2 1 +  2 2 ) −3/2 , so the sum evaluates to O (1) at most.Both the real and complex parts of S are then of order O (/ * ) = O ( v −1/2 ), yielding  ∝  v −1/2 .The -versus- v scaling is insensitive to array geometry, as long as the array is regular, and one has  v ≫ 1.
This paper has been typeset from a T E X/L A T E X file prepared by the author.
Figure1depicts a snapshot of a typical vortex avalanche generated by the N-body code.The snapshot is taken after the avalanche ends and is displayed in the frame corotating with the pinning sites.Vortices that move during the avalanche are marked in red, with blue tracers tracking their discretized paths during the avalanche.Vortices that do not move before, during, or after the avalanche are drawn in grey.The avalanche is triggered by a single unpinning event at the top left of the group of red dots, slightly displaced from the others, which itself triggers the rest of the avalanche by proximity knockon.The unpinned vortices occupy a narrow "wedge", which curves clockwise to the bottom left in response to the condensate's local departure from corotation and the drag due to mutual friction, which corresponds to rotation through a dissipation angle ; see Section 2 inHowitt et al. (2020).The simulation is initialized with  v = 2000 vortices.A total of 21 vortices (out of the  v = 1838 vortices remaining before the avalanche) move during the avalanche, which is the 34th in a sequence starting from a uniform distribution.The local inhomogeneities in the distribution of grey vortices are caused by successive avalanches; many of the grey vortices participate in one of the 33 avalanches leading up to the one depicted in Fig.1.Time-lapse movies of avalanches in action can be viewed by the interested reader in fig.5(strong mutual friction) and fig.9 (weak mutual friction) ofHowitt et al. (2020).

Figure 1 .
Figure 1.Snapshot of the vortex distribution extracted from an N-body simulation of a pinned, decelerating superfluid with the code developed by Howitt et al. (2020).The snapshot is displayed in the frame corotating with the pinning sites.The vortex distribution is shown after the avalanche ends.Vortices that move in the avalanche are drawn as red dots, with a blue tracer showing their discretized path throughout the avalanche.Vortices that do not move are drawn as grey dots.The simulation is initialized by drawing positions of 2×10 3 vortices from a spatially uniform distribution within container radius  * = 10.0.Pinning sites are arranged in a rectangular lattice (lattice size  pin = 0.1, site width  = 0.02) (all lengths measured in simulation units; see Appendix A).The snapshot is taken at simulation time-step  = 826 (midway between the 34 th and 35 th glitches) and the track shows movement starting from  = 798 (midway between the 33 th and 34 th glitches; time in simulation units).Control parameters are recorded in Appendix A; see also Sections 2 and 3 in Howitt et al. (2020).Vortices often annihilate at the boundary during an avalanche, although none does so in this figure.

Figure 2 .
Figure2.Statistics of the dimensionless vortex configuration factor  evaluated in eight independent N-body simulations starting from  v = 5×10 3 with the parameters in Appendix A. Top panel: Scatter plot of snapshots (  v , ) taken after every 10-th glitch in the range 5 × 10 2 ≤  v ≤ 5 × 10 3 for all eight runs. is calculated by evaluating equation (9) at midpoints between avalanches.The red and blue lines each track snapshots from arbitrarily chosen independent runs.Snapshots at more than 90 per cent or less than 10 per cent of the initial container frequency are discarded both in the plot and the parameter estimation.The dark grey line and the grey band displays the mean and 90 per cent credible interval of the posterior predictive distribution at each  v , respectively.Middle three panels: Corner plot displaying the twodimensional and marginalized one-dimensional posterior distributions of the estimated PDF parameters  and  in equation (9).Parameter estimates in the corner plots are reported with their 5-th, 50-th, and 95-th percentiles (vertical dashed lines).Bottom panel: Normalized histogram of / with  =  v  for every (  v , ) point in the top panel for 10 3 samples of (, ) drawn from the posterior.The Rayleigh distribution  (/) = (/) exp − 2 /2 2 is overlaid as a solid black curve.

Figure 3 .
Figure3.Scaling of the dimensionless vortex configuration factor  as a function of the number of vortices  v .Top panel: Samples of  for 10 3 random configurations generated by a uniform Poisson point process for 5 × 10 2 ≤  v ≤ 5 × 10 3 .For each configuration,  v is drawn from a discrete uniform distribution on the interval [5 × 10 2 , 5 × 10 3 ], and  is calculated by evaluating equation (9) for each configuration.The dark grey curve correspond to the median of the theoretical PDF in equation (14) at each  v .The grey band shows the central 90 percentile of the same PDF.Bottom panel: Samples of  for 10 3 realizations of a regular periodic array with randomized centre for 5 × 10 2 ≤  v ≤ 5 × 10 3 .For each configuration, the intervortex spacing  is drawn from a continuous uniform distribution on the interval [2 × 10 −2 , 9 × 10 −2 ],   and   are drawn from a uniform distribution on [0, ), and  is calculated by evaluating equation (9) for each configuration.

Table A1 .
Summary of control parameters for the vortex avalanche simulations in Section 4.