High priority targets for transient gravitational waves from glitching pulsars

Glitching pulsars are expected to be important sources of gravitational waves. In this paper, we explore six different models that propose the emission of transient continuous waves, lasting days to months, coincident with glitches. The maximal gravitational wave energy is calculated for each model, which is then used to determine whether associated gravitational waves could be detectable with LIGO-Virgo-KAGRA's O4 detectors. We provide an analytical approximation to calculate the signal-to-noise ratio which includes information about the source's sky position, improving on previous estimates that assume isotropic or sky and orientation averaged sensitivities. By analysing the entire glitching population, we find that certain models predict detectable signals in O4, whereas others do not. We also rank glitching pulsars by signal-to-noise ratio, based on archival data, and we find that for all models, the Vela pulsar (PSR J0835$-$4510) would provide the strongest signal. Moreover, PSR J0537$-$6910 is not expected to yield a detectable signal in O4, but will start becoming relevant for next generation detectors. Our analysis also extends to the entire pulsar population, regardless of whether they have glitched, and we provide a list of pulsars that would present a significant signal, if they were to glitch. Finally, we apply our analysis to the April 2024 Vela glitch and find that a signal should be detectable under certain models. The non-detection of a supposedly detectable signal would provide an efficiency factor that quantifies a model's contribution to gravitational wave emission, eventually leading to a differentiation of models and independent constraints on physical parameters.


INTRODUCTION
With LIGO-Virgo-KAGRA's (LVK's) Observing Run 4 (O4) currently underway, each day lies a new opportunity to detect a gravitational wave (GW) signal from the Universe's most energetic events.So far, we have detected at least 90 GW signals using ground-based detectors which have all come from the collision of binary compact objects such as black holes and neutron stars (NSs; Abbott et al. 2023).Moreover, several pulsar timing array collaborations have recently reported evidence of a stochastic GW background (Agazie et al. 2023;Antoniadis et al. 2023;Reardon et al. 2023;Xu et al. 2023).It is evident that the era of GW astronomy has well and truly begun.
Although we have successfully made several detections, there are still some types of signals that we have not yet detected, namely, continuous GWs (CWs) and burst signals.Conventionally, CWs are quasi-monochromatic and have durations lasting much longer than a given observational period, such that they can be treated as quasiinfinite.These types of signals are therefore modelled as sinusoids so are relatively straightforward to understand.Burst signals on the other hand do not necessarily have such a simple form.They can be ★ E-mail: g.yim@pku.edu.cn(GY) modelled, e.g. a rapidly decaying sinusoid, or they can be unmodelled, which is what is typically done in practice.
In recent years, more interest has developed in transient CWs, which are sometimes called "long transient signals" (Prix et al. 2011;Abbott et al. 2022a,b).Like conventional CWs, transient CWs are modelled, having the same quasi-monochromatic behaviour but having a finite duration which is longer than a burst signal.Typically, this means that they have a duration of minutes to months.Having a finite duration means that there must also be a start time.Often, transient CWs are treated like conventional CWs but with a window function that encompasses the signal, causing a time evolution to the amplitude that otherwise would not have been there.
Hypothetically, any transient event could trigger transient CWs (or a burst signal), but only certain astrophysical cases have been looked at in detail.One case comes from magnetar bursts and giant flares, which are radiative events categorised by having peak luminosities of ≲ 10 43 erg s −1 and > 10 44 erg s −1 , respectively (Kashiyama & Ioka 2011;Kaspi & Beloborodov 2017).
The other well-known case, which will be the case that we will focus on here, is from NS glitches.A glitch refers to the sudden increase in the spin of a NS, thought to be due to the unpinning of superfluid vortices (Anderson & Itoh 1975).This unpinning causes a rapid transfer of angular momentum from the fast superfluid interior to the crust, which the magnetosphere and electromagnetic radiation are tied to.There is also another popular model, though disfavoured for NSs with large and frequent glitches, which is called the starquake model (Ruderman 1969;Baym & Pines 1971).This model attributes glitches to a rapid decrease in the moment of inertia of the NS when the crust's elastic limit is exceeded during the spin-down of a NS.A review of these glitch models and others can be found in Haskell & Melatos (2015).For some glitches, there is also a glitch recovery where the spin recovers back to, but not necessarily reaching, the preglitch spin frequency.This process normally takes days to months (Antonopoulou et al. 2022;Zhou et al. 2022), similar to the timescales of interest that define transient CWs.
Several groups have extended these glitch models by considering the GWs that are emitted either during or immediately after a glitch.A recent review can be found in Haskell & Jones (2024).Several of these models predict bursts of GWs (e.g.Keer & Jones 2015;Ho et al. 2020a;Yim & Jones 2023;Lopez et al. 2022;Pradhan et al. 2023;Wilson & Ho 2024), but we will only focus on models that predict transient CWs from now on.
Transient CWs from glitches can come from: the liberation of pinned superfluid energy (Prix et al. 2011), Ekman pumping (van Eysden & Melatos 2008;Bennett et al. 2010;Singh 2017), transient mountains (Yim & Jones 2020;Moragues et al. 2023) or trapped ejecta (Yim et al. 2024).In this work, we will explore these models and what will be most important to us will be the amount of emitted GW energy,  GW .We will evaluate model-independent estimates of  GW as well as model-specific values.
The ultimate goal of this work is to provide a list of high priority targets for GW searches.This would allow for efficient use of limited computing power as well as provide electromagnetic observers more reason to obtain more accurate and frequent timing data.Furthermore, having a list of high priority targets will help guide GW searches, especially those using open data where archival data can (and should) be searched over.
In fact, there have already been two GW searches conducted for the signals that we are focusing on here.The first transient CW search from glitching pulsars was done by Keitel et al. (2019) and they looked specifically at the Crab and Vela pulsars which glitched once each during Observing Run 2 (O2).No GW signal was detected so upper limits were set on the GW strain as a function of signal duration.A similar search was done by LVK and other astronomers (Abbott et al. 2022b;Modafferi et al. 2021) for Observing Run 3 (O3), where there were 9 glitches from 6 pulsars, but no GW signals were detected in that effort either.
Alongside trying to detect these signals, we can also make statements about the prospects of detecting them with current and future GW detectors.This is what Yim & Jones (2020) did for the Crab and Vela pulsars using the transient mountain model.Their analysis was based on the signal-to-noise ratio (SNR), a proxy for the coherent F -statistic (Jaranowski et al. 1998;Cutler & Schutz 2005;Dreissigacker et al. 2018), with a fixed threshold for detection.Moragues et al. (2023) extended their work by considering the entire known glitching pulsar population using a model-agnostic approach, as well as specialising to the transient mountain model.In their analysis, they used the semi-coherent F -statistic (Prix et al. 2011) as their detection statistic and the detection threshold was carefully calculated according to Tenorio et al. (2022).
In this work, we return to the relatively simple method of using the SNR to quantify detectability, but we will make some improvements that allow us to include previously neglected features, such as the source's sky position and orientation as well as that of the GW detector.Doing this allows the detectability to be approximated more accurately than before, c.f. Yim & Jones (2020), but also allows for a relatively easy analysis of theoretical models without having to go into the complexities related to the GW search itself, like what was done in Moragues et al. (2023).We will follow Moragues et al. (2023) by generalising our analysis to the entire glitching pulsar population as well as choosing to be agnostic to the specific transient CW model, in addition to specialising to specific models.We go one step further and also look at pulsars that have not glitched, assuming that one day they might do.For the first time, we assess the detectability of transient CWs from Ekman pumping.We also apply our analysis to the recent April 2024 glitch observed in the Vela pulsar (Zubieta et al. 2024;Campbell-Wilson et al. 2024;Grover et al. 2024;Palfreyman 2024;Wang et al. 2024).
In Section 2, we provide an overview of the different models used in this analysis.Then, in Section 3, we present an analytic approximation for the SNR of transient CWs, which accounts for the position and orientation of the source and GW detector.In Section 4, we use this approximation to calculate the detectability of all known glitching pulsars, but we also provide an analysis that is independent of whether the pulsar has been observed to glitch.Moreover, in Section 5, we use our analytic approximation to assess the detectability of the recent April 2024 glitch observed in the Vela pulsar.Finally, Section 6 contains some points of discussion and Section 7 concludes this paper.

GRAVITATIONAL WAVE ENERGY BUDGET
Broadly speaking, the energy budget for GW emission can come from two distinct parts of a glitch: the glitch rise and/or the glitch recovery.One common aspect of glitch rise models is that they are able to explain the short duration spin-up part of the glitch by enforcing the conservation of angular momentum.Post-glitch models are typically agnostic to what causes the spin-up; they are only concerned about the spin-down part which is seen as the recovery.The post-glitch models outlined here also consider angular momentum conservation, but only for the post-glitch regime.There has also been discussion about potential pre-glitch GW emission (Yim & Jones 2020) but we will not consider that possibility here.A recent review of GWs from glitching pulsars can be found in Haskell & Jones (2024).
In addition to these two types of models, there are also "agnostic" models which only consider energy conservation but not angular momentum conservation.These often give energy budgets that are larger than ones found by self-consistent calculations, but they are still useful in that they give firm upper limits on how much energy can be emitted as GWs.
In this section, we outline six models, the first two being "glitch rise" models, the next two being "post-glitch" models, and the last two are "agnostic" models: • Excess superfluid energy model (two components).
Of most interest to us will be the amount of energy that can be released in these models for GW emission, so for each model, we will approximate the (maximal) energy available for GW emission.The energy for each model will be presented in a table in Section 2.4.The energies will then be used to calculate the SNR later in Section 3.

Starquake model (one component)
The starquake model (Ruderman 1969;Baym & Pines 1971) begins with a newly-born rotating NS that is assumed to be initially entirely fluid.The centrifugal force causes an equatorial bulge and after some time, when the NS cools, the body starts to solidify eventually forming a crust that is oblate in shape.The shape at which this crust solidifies is known as the "reference oblateness" and it is defined as when there is no strain in the crust.
As the NS spins-down (primarily through magnetic dipole radiation), the centrifugal force becomes weaker, leading to the NS wanting to become more spherical.However, due to the rigidity of the crust, this is prevented from happening at the expense of elastic strain being built up.
When the strain becomes too large and the elastic limit is reached (at the breaking strain), there is a catastrophic event, called a starquake, which causes the crust to break and the NS rapidly changes shape into one that is less oblate.The moment of inertia decreases when this happens and because angular momentum needs to be conserved, we observe that the NS must spin-up, i.e. the NS has a glitch.
In the simplest starquake model, the NS can be treated as just being one component behaving as an elastic body.The moment of inertia about the rotation axis is  and the angular frequency is given by Ω = 2, where  is the spin frequency of the NS.If we assume that all the excess rotational kinetic energy after the glitch goes into GW production, we find where we have ensured that angular momentum has been conserved where subscript '0' represents quantities before the glitch and 'Δ' quantities show the perturbed value immediately after the glitch, which are assumed to be much smaller than their respective preglitch values, i.e.Δ/ 0 ≪ 1 for a generic variable .To first order in small quantities, we find that the observed glitch size is given by a reduction in the moment of inertia, shown as Note that eq. ( 2) is an exact result and does not require one to consider only first order terms.It is also worth comparing this energy budget with the "naïve" estimate where only energy is conserved, which will be covered later.One finds that  GW from the starquake model is half the value expected from just energy conservation alone, c.f. eq. ( 21).
The results calculated here agree with earlier calculations by Sidery et al. (2010) and Abadie et al. (2011).

Vortex unpinning model (two components)
We will now move onto a more realistic two component set-up for describing NSs.There is great reason to consider NSs with (at least) two components with (at least) one of them being superfluid; without superfluidity, there are several observed phenomena that are difficult to explain such as: large and frequent glitches (Baym & Pines 1971;Pines & Alpar 1985), long post-glitch recovery timescales (Baym et al. 1969a,b;Alpar et al. 1984;Pines & Alpar 1985), and long thermal cooling timescales (Shibazaki & Lamb 1989;Page & Applegate 1992;Larson & Link 1999).More recently, observations of glitches with a "delayed spin-up" (Shaw et al. 2018;Ashton et al. 2019;Shaw et al. 2021) give observational evidence for there being three components (Graber et al. 2018),1 but we will not cover this here and leave it for future work (Yim & Jones, in prep.).For our purposes, two components are enough as this is all that is needed to explain glitch rises and post-glitch recoveries (Anderson & Itoh 1975;Alpar et al. 1984).
The way in which superfluidity describes glitches is as follows (Anderson & Itoh 1975).In the two component model, one component is a "pinned" superfluid inside the NS, and the other component represents everything else, including any "unpinned" superfluid in the crust or core, as well as the normal baryonic crust (Baym et al. 1969b).
A rotating NS means that the superfluid component also rotates, but it can only carry angular momentum by forming an array of vortices where the vortices are assumed to be parallel to the rotation axis.The areal density of these vortices then determines the local rotational velocity of the superfluid, which is valid for regions encasing a sufficiently large number of vortices (Andersson & Comer 2021).
As the NS spins-down on secular timescales, the areal density of vortices must also decrease, so it does this by having vortices migrate radially outwards (on average).However, during this migration, the vortices can become pinned in the core or in the inner crust, preventing further spin-down of the superfluid.This leads to the other (unpinned) component being the sole component that is spinning down.As a result, a lag develops between the two components.
As the NS spins-down more, the lag grows and it becomes harder for the vortices to remain pinned.At some critical lag, the vortices collectively unpin and rapidly transfer their angular momentum to the unpinned component which is what we observe as a glitch.Consequently, the pinned component spins-down and the unpinned component spins-up, resulting in the lag decreasing.Usually, the lag is assumed to be zero after the glitch so that both components co-rotate temporarily.However, this assumption is still unclear and modeldependent.For simplicity, we will assume co-rotation immediately after the glitch.
With this information, we can now calculate how much energy remains after accounting for the changes in rotational kinetic energy of the two components.To do this, we first write down the total angular momentum of the NS system where the subscripts 'p' and 'u' represent the pinned and unpinned components, respectively.The magnetosphere which is responsible for the observed rotational frequency is assumed to be frozen to the unpinned component so we have Ω u = Ω.
We then perturb the above equation and set Δ = 0 to find where we have assumed that the moment of inertias remain constant throughout the glitch.This equation verifies our previous qualitative description, where it was said that the pinned component spins-down as the unpinned component spins-up.Next, we look at the change in rotation kinetic energy of each component, labelled by  = p, u.If both components co-rotate after the glitch at angular frequency Ω co , then the change in rotational kinetic energy of each component is and the total change in energy is the sum of these, giving which is a result also obtained by Prix et al. (2011).Note that the change in energy is negative showing that the final state has less energy than the initial state.We assume this released energy goes directly into GW production meaning where we have used  =  p +  u .Typically, the fractional moment of inertia associated with the pinned superfluid  p /, is comparable to the observed fractional change in spin-down rate Δ Ω/ Ω 0 (Link et al. 1999;Andersson et al. 2012), and from statistical analyses of glitch data, one finds  p / ∼ 0.01 − 0.1 (Andersson et al. 2012;Ho et al. 2015;Fuentes et al. 2017;Ho et al. 2022).It is also possible to calculate  p / from the perspective of microphysics, where an equation of state and superfluid pairing gap model is needed (Andersson et al. 2005(Andersson et al. , 2012)).Depending on model assumptions, it is possible to get reasonable agreement between microphysical theory and glitch observations.However, there is ongoing debate on whether the superfluid in the inner crust is enough to provide the angular momentum reservoir required to power glitch activity when the effect of entrainment is included (Chamel 2012;Andersson et al. 2012;Piekarewicz et al. 2014).For our purposes here, we neglect entrainment.
One notable point about eq. ( 9) is that the energy obtainable in the vortex unpinning model is quadratic in the glitch size, ΔΩ, whereas it is linear for the starquake model.This means calculated energies will be a few orders of magnitude smaller.Aside from Prix et al. (2011), this quadratic relation was also found by Sidery et al. (2010) and Abadie et al. (2011), however, their coefficients were slightly different from what was derived here.We attribute this to a difference in the interpretation of the two components in the case for Sidery et al. (2010), or simply from rounding for Abadie et al. (2011).
Finally, one can also write eq. ( 8) in terms of the lag between the two components,  = Ω p − Ω u , giving where we have assumed co-rotation after the glitch, so that  0 ≡ Ω p,0 − Ω u,0 = ΔΩ u − ΔΩ p , where  0 is the lag just before the glitch.We choose to use eq.( 9) for our expression for  GW as there is only one uncertain parameter,  p , rather than two as in the case for eq.( 10) where  u and  0 are uncertain.

Transient mountain model
We now outline the first of two post-glitch models.For glitches that show a glitch recovery, the rate in change of angular frequency Ω < 0 (also called the spin-down rate | Ω|), usually shows a decrease immediately after the glitch, i.e.Δ Ω < 0. Or, in other words, the spin-down rate increases after the glitch.
In the transient mountain model (Yim & Jones 2020), the increase in spin-down rate is ascribed to the formation of a NS mountain at the moment of the glitch, which radiates away GWs and causes a braking torque on the system.Then, as the glitch exponentially recovers, the mountain dissipates away, reducing the braking torque and causes the spin-down rate to recover back towards the pre-glitch spin-down rate.It is exactly analogous to the case when calculating the "upper spin-down limit" for secularly spinning-down pulsars, but here, it is specific for glitches.
During this recovery, transient CWs are emitted at twice the spin frequency of the NS, i.e.  = 2.If the recovery of the observed spin frequency is exponential with a recovery timescale of  EM , then the amplitude of the transient CWs will also recover exponentially, but with a timescale of  GW = 2 EM .Explicitly, the GW amplitude evolves as where  is the distance to the source and Δ =  −  g is the time that has elapsed after the glitch.The above equation only holds when there is little or no permanent change in the spin-down rate caused by the glitch.In the case that there is, then Δ Ω should be replaced by the "transient" part of Δ Ω that fully recovers on long timescales, i.e. when Δ ≫  EM .
The timing model of a glitch is often modelled as an instantaneous increase in the spin frequency at the time of the glitch  g followed by the sum of an exponential recovery, a linear recovery as well as a permanent offset (Edwards et al. 2006;Yu et al. 2013).Mathematically, this is shown as where Δ p is the permanent offset in spin due to the glitch, Δ  p is the permanent offset in the spin-down rate due to the glitch, and Δ t is the transient change in spin that will fully recover long after the glitch.
One can see that the total spin-up at  =  g is Δ( g ) = Δ p + Δ t and it is standard to define the healing parameter as which tells us the fraction of the total spin-up that recovers on timescales much greater than  EM .For  = 0, there is no recovery so the glitch is simply seen as a step function, whereas  = 1 corresponds to a full recovery with no permanent offset.With this definition, we can write down the expected GW energy from the transient mountain model as (Yim & Jones 2020) If one uses  = 1/2, we get the same energy as the starquake model2 and for  = 1, we get the same as the "naïve" model that will be covered later in Section 2.3.1,where only energy conservation is considered.The value of  varies for different pulsars and even for different glitches from the same pulsar, but to get an idea of their values, the average  for the Crab and Vela pulsars are around ⟨⟩ ∼ 0.8 and ⟨⟩ ∼ 0.2, respectively (Crawford & Demiański 2003;Yim & Jones 2020).

Ekman pumping model
Ekman pumping is a concept taken from fluid mechanics that concerns what happens to a viscous fluid when acted upon by a tangential force at its boundary with a container.For a review, see Benton & Clark (1974).The analogy for NSs is the fluid interior reacting to the rapid spin-up of the solid crust during a glitch.At the fluid-container boundary, the fluid has a no-slip condition, meaning that if the fluid were to have a relative velocity with respect to the container, then there must be a tangential fluid flow gradient that develops.The region in which this gradient exists is called the Ekman layer and the gradients eventually cause meridional flows in the fluid, "pumping" fluid in cells around the container.
There are three phases of Ekman pumping that occur when the container spins-up.Firstly, there is the formation of the Ekman layer which occurs on the spin-up timescale, then there is the spin-up of the fluid in the Ekman layer on the order of the Ekman timescale  E , and lastly, there is the viscous decay of the Ekman layer which occurs on some viscous timescale.For post-glitch recoveries, we are concerned about the second phase where the fluid is increasing its velocity on the Ekman timescale.Greenspan & Howard (1963) were pioneers in this field for the modern era, developing and writing down the governing equations for an incompressible and uniform fluid contained in a cylinder.Walin (1969) improved on this by exploring stratified fluids, and Abney & Epstein (1996) further extended to compressible fluids.The main conclusions were that stratification and compressibility decreased the size of the Ekman layer, effectively decreasing the amount of fluid that needs to spin-up, resulting in smaller Ekman timescales.
Over a decade later, van Eysden & Melatos (2008) realised that it was possible to get GWs from Ekman pumping if non-axisymmetric meriodional flows were excited.They calculated the GW strain from the mass quadrupole and found that it peaked resonantly around frequencies of  = 2 and  =  (the  radiation is only observable for a non-polar observer), with characteristic dimensionless GW strain where  0 is the equilibrium mass density of the fluid ( 0 ∼ 10 17 kg m −3 ),  is some typical length scale of the system ( ∼ 10 4 m),  is the gravitational acceleration at the boundary ( ∼ 10 12 m s −2 ) and  is the distance to the source.The typical timescale for the GW emission is given by the rescaled Ekman timescale where  =  vis /( 2 Ω) is the dimensionless Ekman number with  vis being the viscosity of the fluid interior, and   is some dimensionless parameter that determines how quickly the Ekman flow decays away which depends on stratification, compressibility and how fast the system is spinning (van Eysden & Melatos 2008;Singh 2017). and  label what mode is excited and the leading quadrupole term is represented by  = 2. Soon after, Bennett et al. (2010) found that the current quadrupole contributes more to GW radiation than the mass quadrupole does.The characteristic dimensionless GW strain of the current quadrupole is (Singh 2017) and the GW timescale is the same as the rescaled Ekman timescale provided above in eq. ( 16).Like the mass quadrupole case, emission is at  = 2 and  = , the latter being true only for non-polar observers.The only way to distinguish whether the GW radiation comes from a mass or current quadrupole is by analysing the difference in GW polarisation for a non-polar observer (Singh 2017).Following these calculations, Singh (2017) was able to relax some assumptions and explored a wider parameter space.Crucially, he provided an estimate of the amount of GW energy expected to be radiated from Ekman pumping.First, he estimated the maximum energy that is available from the glitch, which is on the order of how much rotational kinetic energy the crust gains where Γ is the fraction of the total mass contained within the crust which is around Γ ∼ 0.01 (Lorenz et al. 1993;Ravenhall & Pethick 1994).For simplicity, the system evaluated by Singh (2017) and others was assumed to be cylindrical in shape (of radius  and height 2) so  = 2 0  3 .As a result, the glitch energy can now be written as Finally, after simulating the system, Singh (2017) found that only a small fraction  Ek ∼ 10 −7 -10 −5 (for  = 1 -10 kpc) of the glitch energy can be radiated away as GWs making the chance of detection slim.Instead, most of the glitch energy goes into the kinetic and potential energy of the fluid.Nevertheless, we can use eq.( 20) later in our analysis, and we will show quantitatively that it will indeed be difficult to detect GWs from such a mechanism.

Naïve model (one component)
We will now move onto the energy budgets that are agnostic to the detailed mechanisms at play.They only consider energy conservation, as opposed to the aforementioned models, and they result in optimistic energy budgets.Nevertheless, they are interesting to calculate as they provide an upper limit on how much energy will be available for GW radiation during and immediately after a glitch.We begin by looking at the simplest case of a one component system with constant moment of inertia .If it were to spin up, then the amount of rotational kinetic energy that it has gained is simply to leading order in ΔΩ.One might believe that this is the maximal energy that could then go into GWs if all the gained rotational kinetic energy is converted into GWs.In such a case, the increase in rotational kinetic energy must come from somewhere but we remain agnostic to where this energy comes from.It is conceivable that this rotational kinetic energy is gained from mechanisms that occur in the NS's interior, but could equally be due to external factors too.The above equation is identical to the glitch energy used by Ho et al. (2020a) who were concerned about the excitation of  -modes by glitches.The transient CW sources considered here have durations that last much longer than the burst sources considered there but they are similar in the fact that both utilise glitches as their energy source.

Excess superfluid energy model (two components)
Next, we look at the two component model where there is a relatively faster rotating pinned superfluid component contained within an unpinned component.The fact that the pinned superfluid is not co-rotating with the unpinned component means that there is potentially a reservoir of rotational kinetic energy that can be tapped into (or "excess" energy), given by (Prix et al. 2011) to leading order in the lag  = Ω p − Ω, and where  p and Ω p (> Ω) are the moment of inertia and angular frequency of the pinned component.Note that the above energy estimate does not require a glitch to occur.One could imagine that this reservoir of energy is depleted slowly by some dissipative mechanism where energy is lost from the system but angular momentum is not.Clearly, if GW radiation was the underlying mechanism for depletion, then one must account for the angular momentum carried away by the GWs for a self-consistent calculation, but the point of the above calculation is that it provides a firm upper limit for  GW whilst disregarding angular momentum conservation.In this case, the only variable that changes is the lag (angular momentum conservation would require Ω ≠ 0), and so the amount of energy obtainable in this agnostic model is where '' represents the change in a variable over timescales much longer than the glitch rise timescale (which is indicated by 'Δ').We also know  < 0 and the minus signs have been added so that  GW is positive.
It is interesting to note that if one uses the conservation of angular momentum at this point and suppose that the change in the lag is due to a glitch,  = Δ, then by conserving angular momentum in eq. ( 5), we find which leads to which is the same as what was calculated using the naïve model in eq. ( 21), making the two component model no different from the one component model.This is precisely what was done in eq. ( 13) of Prix et al. (2011) and eq.( 6) of Moragues et al. (2023).However, if one wanted to self-consistently account for the conservation of angular momentum, then it should have been done from the beginning, in which case, the GW energy should also include the effect of Ω leading to and when we assume  = Δ and Ω = ΔΩ, and using eq.( 24), we get If the pinned superfluid fraction  p / and/or the normalised preglitch lag  0 /Ω 0 is much smaller than one, then we recover the earlier estimate given by eq. ( 25).
In reality, any lag over a sufficiently long time will tend towards the steady-state lag  ∞ which is non-zero so long as there is a spindown torque acting on the NS.It is usually assumed that the lag is in steady-state just before the glitch, where some instability or event triggers the glitch (Alpar et al. 1994;Andersson et al. 2003;Warszawski et al. 2012).Using the two component model (Baym et al. 1969b), one finds that the normalised steady-state lag is given by where  coup is a coupling timescale between the pinned and unpinned component which is on the order of the glitch recovery timescale (days to months), and  age = −Ω 0 /2 Ω 0 is the characteristic age of the NS, see also Yim & Jones (in prep.).Inputting a range of possible parameters,  p / = [0.01,0.1] and  coup / age = [10 −3 , 10 −4 ], we find that the second term in the parenthesis in eq. ( 27) is smaller by a factor of 10 −8 − 10 −5 compared to the first term, and so can be safely ignored.This results in the GW energy being given by eq. ( 25), regardless of whether angular momentum conservation is self-consistently accounted for or not.

Table of energy budgets
We have now looked at all the energy budget models in detail and we summarise the results in Table 1.In the table, we also provide expressions for  for the different models, which are the coefficients when compared to the agnostic models, defined as where  = 1 for the agnostic models of Section 2.3.One general pattern is that the glitch rise and post-glitch models all produce energies lower than the agnostic models.In this sense, we call the associated agnostic energies the "upper energy limits" that can be obtained from glitches.The "upper spin-down limit" for glitches is given by the transient mountain model.
We can also rank the models by how much GW energy can be obtained.In descending order, we have the agnostic models (both one and two components) being the most optimistic, followed by the transient mountain and starquake models, followed by vortex unpinning and finally the most pessimistic model comes from Ekman pumping, which would typically give energies 10 −9 − 10 −7 smaller than the agnostic case.
In the following sections, we follow the energies obtained from three models: the agnostic case, the transient mountain case and vortex unpinning case.The starquake and Ekman pumping  GW can be simply calculated by multiplying the agnostic case by a factor of 1/2 and 10 −9 − 10 −7 , respectively, whereas the other cases have dependencies on observables which we have data for.

SIGNAL-TO-NOISE RATIO
To connect the GW energy to the SNR, we follow a calculation first done by Prix et al. (2011), as well as Yim & Jones (2020), but we modify it to include geometric effects from the source's sky position and orientation as well as the GW detector's position and orientation on Earth.

Exact description
As usual, the optimal SNR  is defined as (Jaranowski et al. 1998) Table 1.A summary of the GW energies calculated using different models.A comparison to agnostic estimates is given by , where agnostic models have  = 1 and is defined in eq. ( 29).A detailed description of the model parameters can be found under the respective sections in the main text.

Glitch rise Post-glitch Agnostic Starquake
Vortex unpinning Transient mountain Ekman pumping One component Two components where the (•|•) notation represents an inner product defined as where   (  ) is the one-sided power spectral density at GW frequency  , the tilde (˜) represents a Fourier transform, and the asterisk ( * ) represents the complex conjugate.
The GW strain ℎ() (in general relativity) can always be decomposed into two independent polarisations ℎ + and ℎ × as where  + and  × are the detector's antenna pattern functions for the plus and cross polarisations, respectively.ℎ + and ℎ × can be further decomposed into a GW amplitude ℎ 0 , a GW phase Ψ, and some functions that depend on the inclination angle  (the angle between the source's angular momentum vector and our line of sight) and also the wobble angle  (the angle between the source's angular momentum vector and its deformation axis), as illustrated on the right hand side of Fig. 1.The relevant equations can be found in eqs.( 21) and ( 22) of Jaranowski et al. (1998) but we will not explicitly need them here.Additionally, the GW amplitude can be written in terms of the ellipticity where the ellipticity is defined as where   represents the source's moment of inertia about its 'th principal axis.
If there is a time dependence to the GW amplitude, then substituting eq. ( 32) into eq.( 30) for a signal of duration  GW and using Parseval's theorem, we get where we have specialised to GW radiation at twice the source's spin frequency, i.e.  = 2.If one wanted to calculate the SNR for  identical detectors with detector arm opening angles of , then one should multiply the right hand side of this equation by  sin 2 .

Analytic approximations
If the timescale associated with the evolution of the GW amplitude is long compared to the periodic sidereal variation of the antenna pattern functions, i.e.  GW ≫ 1 d, then we can treat ℎ 0 () as constant and pull it out from the integrals to get but we will keep in mind that there is still some time dependence in ℎ 0 in reality.This will be important soon.Jaranowski et al. (1998) wrote down analytical expressions for the antenna pattern functions and after substituting them into the above equation, they found where  2 and  2 are functions of geometric parameters like the right ascension of the source , the declination of the source , the polarisation angle , the source's inclination angle , the detector's latitude  and the detector's orientation , measured as the angle between the bisector of the detector arms anticlockwise from East.These geometric parameters are shown visually in Fig. 1 and the full expressions for  2 and  2 can be found in Appendix B of Jaranowski et al. (1998) but can also be found here in Appendix A for convenience.
There are several assumptions we can make to simplify eq. ( 37) and they generally lead to the SNR being of the form but now we reinstate the time dependence of ℎ 0 by using ℎ 2 0  GW → ∫  GW 0 ℎ 2 0 () , which gives where the coefficient  is determined by what assumptions are being made.For instance,  = 1 corresponds to the assumption that the GW detector is equally sensitive to all parts of the sky such that ) with the addition of  = 0°.This can be seen directly from eq. ( 35).In addition, there is the case where  = 4/25 which is the well-known coefficient derived by Jaranowski et al. (1998) when , ,  and , i.e. the source's sky position and orientation, are averaged over.This latter option for  was exactly what was used by Prix et al. (2011) and Yim & Jones (2020).
Here, we derive yet another  but specific for transient CW signals.It begins by realising that  2 in eq. ( 37) is an oscillating term, with a period of two sidereal days.Therefore, it will always be bound, so for a sufficiently long  GW , the  2  GW term will eventually dominate, leading to  ≈  2 (, , , , ) .
(40) The angle between the angular momentum vector and our line of sight is the inclination angle  and the angle between  and  def is the inclination .
The explicit form of  2 is provided in Appendix A. Eq. ( 40) is what we propose as a suitable choice for calculating the SNR of transient CWs, given  GW is sufficiently large.Moreover, this choice of  allows for the inclusion of additional information, like the declination of the source, polarisation angle, inclination, detector latitude and detector orientation.Note that when using this approximation, the right ascension is not important.This makes sense as for a sufficiently long signal, the Earth would rotate about its axis at least once, meaning each value of the right ascension would carry equal weighting.Table 2 summarises the different options for , their respective assumptions as well as where each choice is applicable.

Accuracy of approximation
Next, one should ask how good of an approximation eq. ( 40) is.To do this, we explicitly evaluate the difference between the approximated SNR, where only the  2  GW term is considered in eq. ( 37), and the actual SNR, given precisely by eq. ( 37).The fractional difference is therefore given by where we have suppressed some variable dependencies for clarity.In particular, recall that  2 and  2 depend on the detector's latitude and orientation, so we expect there to be different fractional differences for different detectors, even when the GW signal has the same properties.The values of  and  for LIGO Hanford, LIGO Livingston and Virgo are provided in Table 3.
Using this, we plot the fractional difference in SNR as a function of GW duration  GW for different detectors.We find that for sufficiently long signals, the error of the approximation tends towards and oscillates around zero.See Fig. 2 for an example.
One can go further and explore, for a given  and , all (, ) parameter space and find the minimum  GW required to have less than some percentage error, say 10%.We will call this the threshold GW duration,  thres .For instance, in Fig. 2,  thres ≈ 1 d for all three detectors.For a given detector, each (, ) maps onto a single value of  thres and in Fig. 3, we show the maximum  thres as a function of , where the maximisation has been performed over all possible  and for all three detectors.In other words, for a given , no  thres would be larger than max( thres ) regardless of what value of  or detector one chooses.To create these figures, we selected  = 0°resulting in there being no dependence on .It also means that we only consider circularly polarised GW radiation, i.e. ℎ + = ℎ × = ℎ 0 .
Working backwards, one can now use Fig. 3 to find the minimum time the GW signal must last for the approximation to be accurate to more than 90%.For example, a source with a declination of  = 70°r equires a GW signal to last for at least  GW ≈ 0.75 d for the approximation error to be less than 10%.One can also conclude that as long as the GW signal lasts for at least  GW ≈ 1.74 d, then the error would be less than 10% regardless of sky position or detector used.Since glitch recoveries typically have timescales on the order of several days to months, then one can safely use the approximation proposed here to estimate the SNR.
For ease of use, we plot  2 (≈ ) as a function of  for different  and for different detectors in Fig. 4, for  = 0°.One can directly read off the appropriate value of  2 and use it in eq. ( 39) to estimate the SNR.

Connecting to E GW
Taking eq. ( 39) as our expression for the SNR, we can rewrite it in terms of the GW energy,  GW .To do this, one recalls that the GW luminosity can be written as and if we impose ℎ 0 = ℎ 0 (), i.e. the GW amplitude has a time dependence, then  = () too, by eq. ( 33).One can then integrate the above equation with respect to time to give and the integral on the right hand side can be eliminated by using the integral of (the square of) eq. ( 33), resulting in recovering eq. ( 3) of Prix et al. (2011).Rearranging the above for the integral and substituting into our expression for the SNR, eq. ( 39), gives the final result where  = 1 if assuming isotropic sensitivity,  = 4/25 if averaging over all sky positions and orientations, or  =  2 for a sufficiently long transient CW signal.Putting in eq. ( 29) and setting  = 2, we find which describes the SNR for models based on glitches.There are a few interesting features to note.Firstly, there is no other dependency on the frequency Ω besides as a normalisation factor in the glitch size and in the power spectral density of the detector.Therefore, at fixed ΔΩ/Ω, the detectability would be determined by the sensitivity of the detector, given all other parameters are fixed.Secondly, there is no explicit dependence on the spin derivative Ω.It is generally accepted that younger NSs, typically with larger spin-down rates, glitch more often and this should be accounted for if talking about the average chance of detecting a signal from a glitch, however, if we are only concerned about an individual glitch, then the detectability should not depend on Ω.
We will use eq.( 46) in the subsequent section to calculate the approximate SNRs of all the models covered in Section 2 for all appropriate observed pulsars, using the assumption of  ≈  2 .

POPULATION STUDY
In this section, we combine all the previous information about different GW energy models and how to approximate the SNR and apply it to the observed pulsar population.The first subsection considers pulsars that have previously glitched, as recorded by the JBCA Glitch Catalogue3 (Espinoza et al. 2011;Basu et al. 2022) and the ATNF Glitch Table 4 (Manchester et al. 2005), and the second subsection drops this criterion, focusing on all 3000+ pulsars found in the ATNF Pulsar Catalogue5 (Manchester et al. 2005).
In this analysis, we use the latest version of the ATNF Pulsar Catalogue, version 2.1.1,which contains 3534 pulsars.Since it is imperative to have the spin frequency and distance to calculate the SNR, we remove all pulsars that do not have these recorded, leaving 3402 pulsars remaining in our dataset.3. The inclination of the source was set to  = 0°here.

Including glitch data
As for the glitch catalogues, we first analyse the JBCA Glitch Catalogue.Unprocessed, it contains 672 glitches from 207 pulsars. 6We go through a process of cleaning which includes: removing glitches that do not provide Δ/, removing antiglitches, removing pulsars with no distance measurement, and editing some J-names to be consistent with the other databases.The specific details of the cleaning procedure can be found in Appendix B. After cleaning, we are left with 662 glitches from 202 pulsars.
We do a similar cleaning process for the ATNF Glitch Table, which originally contains 626 glitches from 211 pulsars, 7 but this is reduced to 616 glitches from 205 pulsars after cleaning.Again, the details of the cleaning procedure can be found in Appendix B. 6 When arranged by J-name, there are 208 pulsars, but this is due to a typo where the J-name of the MJD 55367 glitch of PSR J1913+0446 has been accidentally entered as "J1914+0446".After correcting this, we count 207 unique J-names, not 225 as stated at the top of the catalogue. 7We have not counted any extra recovery components as a glitch, as they are already attributed to an existing glitch in the dataset.Explicitly, the dataset contains 626 glitches, but 15 of these also have an extra recovery component (so are modelled by 2 exponential recoveries) and there is also a glitch that has 3 recovery components, and another with 4 recovery components.In total, the dataset has 626 glitches with 18 recovery components (making a total of 644 entries) from 211 pulsars.
Clearly, there are some differences between the two glitch catalogues but one of the most obvious is the fact that the ATNF Glitch Table provides values of the healing parameter  alongside the glitches where  has been measured.This is particularly important for the transient mountain model (Section 2.2.1).As such, we create a separate dataset, which is a subset of the ATNF Glitch Table , that contains all the glitches that have .These glitches can sometimes have more than one recovery component, in which case, each component carries its own  and recovery timescale  EM .In the end, there are 114 independent glitches, with some of these showing multiple glitch recovery components (18 extra components in total), making a total of 132 entries from 57 pulsars.
Of course, there will be a great deal of overlap between these glitch catalogues.We therefore need to choose which will be our "base" catalogue, where we obtain Δ/ from by default, and which we will add data from.In this work, we choose the JBCA Glitch Catalogue as our base and add onto it any missing glitches that come from the ATNF Glitch Table .In Appendix C, we provide the list of pulsars that are missing from the JBCA Glitch Catalogue but are present in the ATNF Glitch Table (and vice versa).
Dataset A contains Δ/ from 686 glitches from 219 pulsars and Dataset B contains both Δ/ and  from 114 glitches (132 entries) from 57 pulsars.To both datasets, we append pulsar information such as , ,  and  from the ATNF Pulsar Catalogue.
As mentioned at the end of Section 2, we will explicitly consider only three models: the agnostic case ( = 1), the vortex unpinning case ( = (1/2) (Δ/) (/ p − 1)) and the transient mountain case ( = ).The other models can be obtained by simple rescaling of the agnostic case.Dataset A will be used to calculate the first two cases and we will use  p / = 0.01 for the vortex unpinning model (Andersson et al. 2012).Dataset B will be used for the transient mountain model.
We present the results of the SNR calculation as histograms in Fig. 8.When using  =  2 , each GW detector has a slightly different SNR due to the dependence of  2 on (, ).Therefore, we choose to provide histograms for the largest SNR obtained from the three detectors.The sensitivities that we use are O4 sensitivities and we obtain the sensitivity curves from LIGO Document T2200043-v3. 8 The orange-filled histograms represent the case where sky positions have been accounted for, by using  =  2 in the most optimal orientation of  = 0°.In this case, one sees that the agnostic models (top panel) and the transient mountain model (bottom panel) predict a GW signal would have been detectable ( ≳ 10) if the glitch were to happen during O4.From eq. ( 46), one sees that  ∝ √ , and so to get to the starquake and Ekman pumping m odels, one rescales the agnostic SNR by 1/ √ 2 and 10 −4.5 − 10 −3.5 , respectively.It is therefore apparent that the SNRs are still detectable for the starquake model, but Ekman pumping would give SNRs no larger than ∼ 0.1, rendering the associated GWs completely undetectable with O4 detectors.Even with next generation detectors, where the sensitivity 8 https://dcc.ligo.org/LIGO-T2200043/public is an order of magnitude better, Ekman pumping will not give detectable GWs.Similarly, the vortex unpinning model (middle panel) will not give detectable GWs with O4 detectors, but with next generation detectors, there could potentially be a marginal detection.Not only will the improvement of the sensitivity matter, but the location of the detectors will also affect the detectability.
Also plotted (as solid black lines) are the histograms that would be obtained if one averaged over sky position and orientation, corresponding to  = 4/25.One finds that using this approximation shifts the SNR distribution towards lower values when compared to the  =  2 ( = 0°) case.This shifting to lower SNRs can also be achieved by having a larger value for .Setting  = 60°(and  = 0°) leads to a SNR distribution that is similar to the sky and orientation averaged case (dashed blue line).This is consistent with the  = 60°c urves in Fig. 4 being close to the value of 4/25.Therefore, one can conclude that unless  ≳ 60°, the  =  2 approximation leads to larger SNRs than the sky and orientation averaged case.However, note that the SNR here only captures the  = 2 radiation but when  > 0°, the  =  radiation will contribute too.We do not worry about this here but it will have the effect of changing the value of  that gives an equivalence to the  = 4/25 case.
We will now look at the most favourable pulsars that we expect to yield a detection.In Tables 4 -6, we present the top 15 pulsars that give that largest GW energies and SNRs assuming different energy models using  =  2 ( = 0°).For each pulsar, we also present its right ascension, declination, spin frequency, spin frequency derivative, distance, number of observed glitches, glitch size corresponding to the largest SNR, fractional change in the spin-down rate corresponding to the largest SNR, and for the transient mountain model, the healing parameter corresponding to the largest SNR.
One can immediately see results that are consistent with the conclusions drawn from Fig. 8, namely, the agnostic and transient mountain models give detectable signals with O4 but the vortex unpinning model cannot.There are 15 pulsars that have max() ≥ 10 using the agnostic model, none for the vortex unpinning model and two pulsars for the transient mountain model, PSR J0835−4510 (Vela) and PSR J0205+6449.
It is also worth looking at individual pulsars too, particularly the ones that have already been targeted by LVK.These are the Crab pulsar (PSR J0534+2200), the Vela pulsar (PSR J0835−4510) and the "big glitcher" (PSR J0537−6910), of which the latter has been suggested to also emit GWs via -modes (e.g.Fesik & Papa 2020;Abbott et al. 2021b).Firstly, we track the Crab pulsar.It ranks 8th, 15th and 3rd in the agnostic, vortex unpinning and transient mountain models, respectively.The drop in ranking for the vortex unpinning model is due to the relatively small glitch sizes that the Crab exhibits, as the GW energy in this model is quadratic in the glitch size rather than linear.The Crab ranks relatively highly for the transient mountain model due to its large .Moving onto Vela, it remains the highest priority target for all models primarily due to its large glitches and close distance, holding also true for PSR J0940−5428.
Finally, we found that PSR J0537−6910 is not included in the top 15 for any of our models due primarily to its large distance from us ( = 49.7 kpc).This is despite it having the largest  GW ∼ 1 × 10 44 erg in our dataset, which comes from it having the 3rd largest spin frequency, only to be beaten by millisecond pulsars PSR J1824−2452A ( ≈ 327.4 Hz; Lyne et al. 1987;Reardon et al. 2021) and PSR J0613−0200 ( ≈ 326.6 Hz; Lorimer et al. 1995;Reardon et al. 2021).However, these millisecond pulsars exhibited the smallest glitches of our dataset, of sizes Δ/ = 8 × 10 −12 (Cognard & Backer 2004;Espinoza et al. 2011) and Δ/ = 2.5 × 10 −12 (McKee et al. 2016), respectively, meaning their GW energies are only  GW ∼ 10 40 erg.One major reason that PSR J0537−6910 is interesting is that it glitches on average every ∼ 100 d but with max() = 2.7 (using the agnostic models), any GWs are unlikely to be detected, even if we coherently stack multiple signals together.A detectable signal would require stacking more than 15 events which would take over 4 years.It is clear that glitches from PSR J0537−6910 will only start being probed by GWs with next generation detectors, a conclusion in agreement with Yim (2022).

Ranking all pulsars
In this section, we look at all pulsars from the ATNF Pulsar Catalogue (with known  and ) and assume that at some point in the future, it will glitch, even if it has not glitched before.As a result, we will be able to rank which pulsars would give the most significant signal if a glitch were to occur.
To do this, we fix  GW (and ), and we calculate the parameter  for each pulsar, defined as such that the SNR can be calculated straightforwardly as which comes from eq. ( 45).For  = /2,  is simply the constant of proportionality between  and √  GW .Note that  is a property that is unique for each pulsar and GW detector combination, as it depends on the source's sky location and orientation, as well as the detector's location, orientation and sensitivity.
For definiteness, we set  = 2 and  = 0°, and we calculate  for every pulsar in the ATNF catalogue.A given pulsar will have a different  for each GW detector but we only report the largest value.The ranked results of the top 100 pulsars are presented in Table 7.
From here, we see that the Vela pulsar is the pulsar with the 8th highest , with a value of  ≈ 3.99×10 −20 erg − 1 2 .We also report the values for the Crab pulsar and PSR J0537−6910 which are  ≈ 5.61× 10 −21 erg − 1 2 and  ≈ 2.60×10 −22 erg − 1 2 , respectively.We can verify our earlier calculations, for example, for Vela, by substituting  GW = 1.53 × 10 43 erg, 2.35 × 10 39 erg and 1.50 × 10 42 erg into eq.( 48) for the agnostic, vortex unpinning and transient mountain models and we recover the maximum SNRs in Tables 4 -6 to within some rounding error.This makes  a powerful tool to quickly estimate the SNR if  GW is provided.Additionally, the calculation no longer needs to be restricted to glitches but could also be used for other transient phenomena such as magnetar bursts or fast radio bursts, if one assumes  = 2 radiation is emitted.
We can also do the same by fixing ΔΩ/Ω and in a similar way, we define such that the SNR can be calculated as In Table 8, we present the top 100 pulsars with the highest Υ, assuming  = 2 and  = 0°.Note that none of these top 100 pulsars are known glitchers and in fact, the first glitcher in the ranked list comes in 165th place, which is the Vela pulsar with Υ ≈ 8.86 × 10 4 .The Crab pulsar has Υ ≈ 3.34 × 10 4 and PSR J0537−6910 has Υ ≈ 3.21 × 10 3 .Better electromagnetic timing of these top 100 pulsars would unveil better constraints on the glitch mechanism if a glitch were to occur and it was coincident with observing GW detectors.
Table 4.A table showing the top 15 pulsars that give the largest SNRs when using agnostic models.Starting with the left-most column, we have: pulsar J-name, right ascension, declination, spin frequency, spin frequency derivative, distance, number of glitches observed, glitch size that gives the largest SNR, fractional change in spin-down rate corresponding to the largest SNR, largest GW energy, largest SNR.Columns 2 -6 were taken from the ATNF Pulsar Catalogue (Manchester et al. 2005), columns 7 -9 were taken from the JBCA Glitch Catalogue (Espinoza et al. 2011;Basu et al. 2022) and columns 10 -11 were calculated using eqs.( 21) and ( 45 (Manchester et al. 2005), columns 7 -9 were taken from the JBCA Glitch Catalogue (Espinoza et al. 2011;Basu et al. 2022) and columns 10 -11 were calculated using eqs.( 9) and ( 45

APRIL 2024 GLITCH IN VELA
In this section, we apply our analysis to the recent glitch that was observed from the Vela pulsar (Zubieta et al. 2024;Campbell-Wilson et al. 2024;Grover et al. 2024;Palfreyman 2024;Wang et al. 2024).
Although the glitch was first announced by Zubieta et al. (2024), we use the glitch parameters reported by Palfreyman (2024) as his data and analysis was able to narrow the glitch time to within a ∼ 7 second window.The glitch occurred on 29th April 2024 between 20:52:11.4and 20:52:18.1 and the glitch size was ΔΩ/Ω ≈ 2.4 × 10 −6 .
Fortunately, LIGO Hanford, LIGO Livingston and Virgo were all observing when the glitch happened, as part of O4.At that time, the sensitivities of the LIGO Hanford, LIGO Livingston and Virgo detectors, at  = 2 ≈ 22 Hz, were √︁   (  ) ≈ 8 × 10 −23 Hz − 1 2 , √︁   (  ) ≈ 4 × 10 −23 Hz − 1 2 and √︁   (  ) ≈ 3 × 10 −22 Hz − 1 2 , respectively. 9sing the approximation of  =  2 , with  = 0°, we calculate the maximum SNR (maximised over the three detectors), max(), for this event and we find • Agnostic model: max() = 137.8,Table 6.A table showing the top 15 pulsars that give the largest SNRs when using the transient mountain model.Starting with the left-most column, we have: pulsar J-name, right ascension, declination, spin frequency, spin frequency derivative, distance, number of glitches observed, glitch size that gives the largest SNR, fractional change in spin-down rate corresponding to the largest SNR, healing parameter that gives the largest SNR, largest GW energy, largest SNR.Columns 2 -6 were taken from the ATNF Pulsar Catalogue (Manchester et al. 2005), columns 7 -9 were taken from the JBCA Glitch Catalogue (Espinoza et al. 2011;Basu et al. 2022), column 10 was taken from the ATNF Glitch Table (Manchester et al. 2005) and columns 11 -12 were calculated using eqs.( 14 • Vortex unpinning model: max() = 1.5, • Transient mountain model: max() = 61.6,where we used  p / = 0.01 for the vortex unpinning model and  = 0.2 for the transient mountain model, which is the historical average for the Vela pulsar.Excitingly, some of these transient CW models can now be probed with the latest Vela glitch.In the event of a non-detection, upper limits can be placed on how much each of these models can contribute towards GW radiation which in turn puts constraints on the physical parameters that govern the models.However, it should be noted that the SNRs calculated here are likely to be overestimates since there are other ways for the glitch energy to be used up.For instance, it is believed that  -modes of the NS could be excited during a glitch (Ho et al. 2020a;Yim & Jones 2023;Wilson & Ho 2024).Interestingly, this could also radiate GWs.Following Yim & Jones (2023), we calculate the associated SNR coming from the excitation and decay of  -modes and find that one could get a burst signal, with frequency  ∼ 2 kHz and duration  ∼ 0.1 s, with a SNR of  = 25, 50, 7 for LIGO Hanford, LIGO Livingston and Virgo, respectively.Note that this  -mode calculation used  = 1 as the signal is not long enough for the  =  2 approximation to be valid.This can only give a rough estimate of the SNR so the  -mode calculation should only be taken as indicative.

DISCUSSION
In this work, we have often taken the most optimistic scenario for detecting transient CWs, such as the inclination angle being zero, i.e.  = 0°.Since we observe pulsars pulse, we know that we cannot have  = 0°, meaning that we will never truly have just pure  = 2 radiation but rather a mix containing  =  radiation as well.One should be able to replicate the analysis presented here but track the  =  radiation too, but we will leave this for a future study.
If  is non-zero, then one might ask what a suitable value would be.We can obtain this from electromagnetic observations, specifically from radio polarisation data, which are often fit with the "rotating vector model" (Radhakrishnan & Cooke 1969;Everett & Weisberg 2001;Jones 2007).Using this model, one obtains a value of the magnetic inclination angle ( RVM , the angle between the pulsar's rotation axis and magnetic axis) and the viewing angle ( RVM , the angle between the magnetic axis and the observer's line of sight).The inclination  is therefore just the sum of these two angles,  =  RVM +  RVM (or  = 180°−  RVM −  RVM , see Jones (2015) and Abbott et al. (2017)).Future work should incorporate this inclination dependence into the analysis to get a more faithful SNR calculation.
In any case, one can still take the assumption of  = 0°and define upper limits.The agnostic models ( = 1) define an "upper energy limit" and the transient mountain model defines an "upper spin-down limit".In the case where one of these models predicts a detectable signal but is, in reality, not detected, then one can define an efficiency factor  that tells us the largest fraction of, say, the upper energy limit that goes into GW emission.For instance, if the value of  GW corresponds to a detectable signal ( >  thres ), but a signal was not detected in reality, then we know the energy emitted must be smaller than  GW , where  < 1.Specifically, the efficiency factor can be calculated as where the subscript  shows variables that are model-dependent and  thres is the threshold SNR used to claim a detection.For instance, for the recent Vela glitch covered in Section 5, the agnostic model predicted a maximum SNR of 137.8 (where the maximisation refers to maximising across all three detectors) and the transient mountain model predicted a maximum SNR of 61.6.Let us take  thres = 10.If after searching, we do not find a signal, then we are able to place limits on the efficiencies of each model:  en ≤ 0.00526 and  sd ≤ 0.02635 for the agnostic (energy upper limit) and transient mountain (spin-down limit) models, respectively.Note that one can easily convert from upper energy efficiencies to than this, we argue that the approximation is suitable, especially if only getting a sense of detectability.
Using this novel approximation, we applied the detectability calculation to the known population of glitching pulsars.We found that some models (agnostic, transient mountain, starquakes) predict a detectable signal in O4, whereas the vortex unpinning and Ekman pumping models do not predict detectable signals in O4.However, with next generation detectors, it is possible that the vortex unpinning model could provide a marginal detection, but the Ekman pumping model would not.
We compared the SNRs obtained from the novel  =  2 approximation with a commonly-used assumption of sky and orientation averaging ( = 4/25) and we found that the  =  2 approximation leads to more optimistic SNRs, so long as the source's inclination is smaller than around ∼60°.
From these SNR calculations, we were also able to provide tables of the highest priority targets, dependent on the energy model assumed.In all models, the Vela pulsar was the top priority as it gave the largest detectability estimates.The Crab pulsar was also interesting, especially with the number of glitches that has been observed coming from it, but we also suggest PSR J0940−5428, PSR J1952+3252, PSR J0205+6449 and PSR J1813−1246 could potentially be even more interesting.Unfortunately, glitch-induced transient CWs from PSR J0537−6910 do not appear to be detectable with current generation detectors, with it only becoming constraining for models during the next generation of detectors.
Next, we no longer restricted ourselves to the population of glitching pulsars and applied the analysis to the entire pulsar population, assuming that one day, some of these pulsars may glitch.We provided two lists of high priority sources (Tables 7 and 8), and in the latter, the top 164 pulsars have not been observed to glitch before.We advocate for better timing programmes for these potentially interesting pulsars.
Also, we formulated two new parameters,  and Υ, which can be used to quickly calculate the SNR for a given pulsar-detector combination.The SNR can be calculated by multiplying  by the square root of the GW energy, or similarly, from multiplying Υ by the square root of the dimensionless glitch size.
Finally, we applied our analysis to the recent April 2024 glitch from the Vela pulsar, which was coincident with O4.We found that the agnostic and transient mountain models give a detectable signal, whereas the vortex unpinning model does not.It is likely that LVK will conduct a search for this GW signal and in the event of a nondetection, we will be able to place upper limits on how much energy is radiated away as GWs, placing constraints on the models.In this case, this event would correspond to the first time the upper energy (and spin-down) limit for glitches would have been beaten.
From this analysis and Vela's most recent glitch, it is clear that the study of GWs from pulsar glitches is promising and timely.It will offer a new opportunity to independently probe NSs and with current and upcoming GW detectors, we should be able to place constraints on physical parameters that allow us to test our understanding of superfluidity, elasticity/plasticity, viscosity, magnetic fields, temperature gradients, and so on.More resources need to be dedicated towards this effort so that hopefully, when a signal is detected, we will be ready to maximise the science that can come out from it.

Figure 1 .
Figure 1.A figure showing the different geometric parameters that concern the GW source and detector.Starting from the left, the diagram shows an L-shaped ( = /2) GW detector as bold solid lines.The angle between the bisector of the detector measured anticlockwise from East is .The central dot represents the centre of the Earth and  is the detector's latitute on Earth.Moving onto the central diagram,  and  represent the source's right ascension and declination respectively, with the right ascension measured from the meridian that passes through the vernal equinox. is the polarisation angle and tells us how rotated the source is about our line of sight.Finally, on the right is a diagram of a deformed NS, which has a deformation axis  def and angular momentum vector .The angle between the angular momentum vector and our line of sight is the inclination angle  and the angle between  and  def is the inclination .

Figure 2 .
Figure 2. A figure showing an example of how the fractional difference between the approximated SNR and the exact SNR vary with GW signal duration.It can be clearly seen that the error induced by the approximation decreases for increasing  GW .After around  GW ≈ 1 d, the error drops to less than 10% for all three detectors.The parameters used to create this figure were  = 40°,  = 30°and  = 0°.

Figure 3 .
Figure 3.A figure showing the threshold GW duration a signal must have in order for the SNR approximation to have less than 10% error, as a function of declination .The maximisation on  thres was done across all values of right ascension  and for all three detectors found in Table3.The inclination of the source was set to  = 0°here.

Figure 4 .
Figure 4.A figure showing  2 as a function of declination , for different inclinations  and different GW detectors.Linestyles of solid, dotted and dash-dots represent values of  = 0°, 30°, 60°, respectively, and colours of blue, red and grey represent LIGO Hanford, LIGO Livingston and Virgo, respectively.The horizontal black line at  2 = 4/25 represents the equivalent value one gets if one averages over all sky positions and orientations (Jaranowski et al. 1998).

Figure 5 .
Figure 5.A figure showing the fractional difference between the approximate and exact SNR as a function of GW signal duration  GW for the Crab pulsar.Values of  = 0°,  = 83.6°and = 22.0°were used, where the latter two parameters were obtained from the ATNF Pulsar Catalogue(Manchester et al. 2005).One finds  thres = 0.36 d (for < 10% error) for the Crab pulsar.

Figure 6 .
Figure 6.A figure showing the fractional difference between the approximate and exact SNR as a function of GW signal duration  GW for the Vela pulsar.Values of  = 0°,  = 128.8°and = −45.2°wereused, where the latter two parameters were obtained from the ATNF Pulsar Catalogue(Manchester et al. 2005).One finds  thres = 0.77 d (for < 10% error) for the Vela pulsar.

Figure 7 .
Figure 7.A figure showing the fractional difference between the approximate and exact SNR as a function of GW signal duration  GW for PSR J0537−6910.Values of  = 0°,  = 84.4°and = −69.2°wereused, where the latter two parameters were obtained from the ATNF Pulsar Catalogue (Manchester et 2005).One finds  thres = 0.30 d (for < 10% error) for PSR J0537−6910.

Figure 8 .
Figure 8. Histograms of the maximum SNR, referring to the largest SNR from the three GW detectors.Top panel: agnostic models, middle panel: vortex unpinning model, bottom panel: transient mountain model.The orange-filled histograms represent the SNRs that are obtained when  =  2 for  = 0°.The solid black line represents the case when  = 4/25, i.e. when the source's sky position and orientation are averaged over.The dashed blue line represents the case when  = 60°and  = 0°, which recovers a similar SNR distribution as the  = 4/25 case.

Table 2 .
A table containing the different options for , their respective assumptions and the scenarios where their use may be appropriate.

Table 3 .
A table containing the latitude  and orientation  parameters of LIGO Hanford, LIGO Livingston and Virgo, taken from Jaranowski et al.

Table 5 .
), respectively.A table showing the top 15 pulsars that give the largest SNRs when using the vortex unpinning model.Starting with the left-most column, we have: pulsar J-name, right ascension, declination, spin frequency, spin frequency derivative, distance, number of glitches observed, glitch size that gives the largest SNR, fractional change in spin-down rate corresponding to the largest SNR, largest GW energy, largest SNR.Columns 2 -6 were taken from the ATNF Pulsar Catalogue ), respectively.