Transient gravitational waves from pulsar post-glitch recoveries

Traditionally, searches for GWs have been directed towards detecting compact binary coalescences (CBCs), continuous GWs (CWs), burst GWs and stochastic GWs. Of most relevance here are CWs which are GWs with long durations compared to the time spent observing them and are quasi-monochromatic in nature. CBC GWs and CWs have modelled waveforms meaning the method of matched-filtering can be used to detect them. Both burst and stochastic GWs are generally unmodelled making them particularly difficult to detect. Out of these four GW groups, only CBC GWs have been detected so far (Abbott et al., 2019a).

In this work, we will be focusing on transient CWs from individual NSs (e.g. Prix et al., 2011). Transient CWs, a type of long-duration transient GW, are CWs which have durations between the time-scales of CBCs and conventional CWs. Roughly speaking, this means they have durations on the order of hours to weeks, well within the duration of any Advanced Virgo or Advanced LIGO (aLIGO) observation run which is on the order of several months (Acernese et al., 2015; Aasi et al., 2015).

So far, two types of complementary long-duration transient GW searches have been developed. One from Thrane et al. (2011), who took a burst search (looking for unmodelled GWs on the duration of $\mathcal{O}$(seconds)) and extended it to include GWs in our time-scales of interest. The other type came from the creators of the term “transient CWs”, Prix et al. (2011), who took a conventional (modelled) CW search and shortened the CW duration instead. There are other examples of the latter in Abbott et al. (2019d).

There have been a couple of attempts to use these searches to detect long-duration transient GWs. Most recently, Keitel et al. (2019) tried to find transient CWs in the post-glitch recoveries of the Crab and Vela pulsars but unfortunately, no transient CWs were detected. Abbott et al. (2019d) did a similar search but for unmodelled transient GWs after the binary NS merger GW170817, for durations of up to 8.5 days. Again, no transient GWs were detected.

There have been many other searches conducted by aLIGO, but these were for transient GWs with durations shorter than the time-scales we are interested in here, which includes burst searches. These include transient GWs after magnetar bursts (Abbott et al., 2019c), during the ringdown (up to 500 s after) of GW170817 (Abbott et al., 2017), after long gamma-ray bursts (Abbott et al., 2019e) and those due to NS oscillations after the Vela pulsar glitched in 2006 (Abadie et al., 2011a). There was also an all-sky search within O2 data for transient GWs which were 2 – 500 s in duration but besides from recovering GW170817, no other transients were detected (Abbott et al., 2019b).

On the other end of the GW spectrum, the CW spin-down limit for the Crab and Vela pulsars were beaten in 2008 and 2011 respectively (Abbott et al., 2008; Abadie et al., 2011b) with recently updated GW strain upper limits of $h_{0}^{\text{Crab}}\lesssim 2\times 10^{-26}$ and $h_{0}^{\text{Vela}}\lesssim 2\times 10^{-25}$ (Abbott et al., 2020). These upper limits were calculated assuming a coherent CW signal for the duration of the observing time. It is perfectly acceptable for transient CWs to have values of $h_{0}$ larger than these upper limits due to their shorter duration.

It is evident that there are searches being conducted for transient CWs but thus far, there has been little theoretical modelling of how these transient CWs come about. Pulsar glitches, spontaneous step-like increases in a NS’s spin frequency, have often been explored as possible triggers. The idea of Ekman pumping after pulsar glitches, which is the bulk movement of fluids due to a tangential force (in the case of a glitching two-component NS, it is the viscous shear force between the crust and core), has been suggested (van Eysden & Melatos, 2008; Bennett et al., 2010; Singh, 2017) as well as polar Alfvén waves from magnetar flares (Kashiyama & Ioka, 2011). Additionally, Prix et al. (2011) created a toy model proposing the excess superfluid energy from a pulsar glitch could seed a non-axisymmetric deformation leading to transient CW emission. One obstacle for this model is it depends on knowing the lag between the superfluid core and the crust or the superfluid angular velocity, both of which we cannot currently observe.

Here, we propose another transient CW model which is simple and makes precise and falsifiable predictions. Like other models, it uses pulsar glitches as the trigger and with our model, we are able to naturally explain the post-glitch recovery found immediately after a pulsar glitch. We concern ourselves with only the post-glitch recovery and not the cause of the glitch.

Fig. 1 shows the usual picture of a pulsar glitch. Tracing from the left, we see the NS’s spin frequency, $\nu$, decreases on a secular time-scale until the point of the glitch which occurs at time $t_{\text{g}}$. The glitch is defined as an instantaneous increase in a NS’s spin frequency ($\Delta\nu(t_{\text{g}})>0$) whilst the time derivative of the spin frequency simultaneously decreases ($\Delta\dot{\nu}(t_{\text{g}})<0$). Then, immediately after the glitch is the post-glitch recovery where the spin frequency and its time derivative are observed to be well-modelled by an exponential recovery (e.g. Lyne et al., 2000; Espinoza et al., 2011). Both these parameters generally do not fully recover to pre-glitch values (e.g. Haskell & Melatos, 2015).

In our model, we propose the decrease in $\dot{\nu}$ (i.e. an increase in the spin-down rate, $\Delta|\dot{\nu}|>0$) during the post-glitch recovery is due to an external braking torque caused by a NS “mountain”. This corresponds to a change to the shape of the NS into a tri-axial ellipsoid (a $l=2$, $m=2$ deformation in the language of spherical harmonics). Along with rotation, there is enough within the model to give the time-varying mass quadrupole moment required for GW emission. The rest of the post-glitch recovery is explained by the mountain slowly dissipating away, on a similar time-scale as the glitch recovery time-scale, $\tau$, leading to a reduction in the braking torque until the spin-down rate recovers. The dissipation could be caused by plastic flow (e.g. Baiko & Chugunov, 2018) of the NS crust or magnetic diffusion (e.g. Pons & Viganò, 2019) but we do not focus on any particular mechanism here. Similarly, we do not focus on the formation mechanism of transient mountains. We only explore the consequences of such mountains. These transient mountains would give rise to transient CWs.

It should be noted that our model is not the usual explanation for the post-glitch recovery which is normally attributed to the model of vortex creep, pioneered by Alpar et al. (1984). The vortex creep model builds on the vortex unpinning glitch model of Anderson & Itoh (1975) which relies on there being two components to a NS, a normal component and a superfluid component (Baym et al., 1969). It is generally accepted that large glitches require some superfluid contribution (Alpar et al., 1993) but in our simple model, we choose to ignore any internal effects.

In this paper, we take our transient mountain model and calculate whether the post-glitch transient CWs given off would be detectable with current and future GW detectors. Our model is largely motivated by the fact we have observations of pulsar glitches at radio wavelengths, and so we have a set of observables to shape our model. This is particularly topical with the rapid developments seen recently within multi-messenger astronomy. We also aim to interpret and provide physical context to why searches such as Keitel et al. (2019) should continue in the future.

In Section 2, we introduce and develop relevant equations from radio pulsar astronomy that describe glitches and the subsequent recovery. In Section 3, we look at the dynamics of the system which allows us to calculate the ellipticity required to cause the post-glitch recovery. This, in turn, tells us the GW strain achievable. In Section 4, we find an expression for the total GW energy radiated away by a transient mountain. In Section 5, we relate this GW energy to the signal-to-noise ratio (SNR) achievable in a GW detector. In Section 6, we apply our model to radio data and give a table of results including the SNRs for Crab and Vela glitches for both aLIGO and the Einstein Telescope (ET). In Section 7, we discuss the predictions and limitations of our model. Finally, Section 8 summarises all the results of this work and concludes with some final remarks on future prospects of this research.

Timing a pulsar requires accurate modelling of the rotational dynamics of a NS. One well-known property of NSs is that they spin-down on long (secular) time-scales. This is thought to be primarily due to magnetic dipole radiation (e.g. Lyne & Graham-Smith, 2012) though searches are being performed to see if CWs could also play a role (e.g. Abbott et al., 2020). This secular spin-down is represented by the lower dashed line in Fig. 1. Since we are only interested in glitches (which have much shorter time-scales), we can subtract the secular spin-down to leave the change in the spin frequency due to the glitch, $\Delta\nu(t)$, represented by

where $\Delta t=t-t_{\text{g}}$ is the time elapsed after the glitch with $t_{\text{g}}$ being the time of the glitch. $\Delta\nu_{\text{p}}$ refers to the change in spin frequency which is permanent due to the glitch and $\Delta\nu_{\text{t}}$ is the same but refers to the transient part, i.e. the change in the spin frequency which fully recovers on time-scales much larger than the recovery time-scale of the glitch, $\tau$. The product $\Delta\dot{\nu}_{\text{p}}\cdot\Delta t$ represents the contribution to $\Delta\nu(t)$ due to a permanent change in the spin-down rate caused by the glitch, $\Delta\dot{\nu}_{\text{p}}$. From this formulation, it is seen at $t=t_{\text{g}}$, we get the relation

We can then differentiate Equation (1) to get the time derivative of the spin frequency, $\Delta\dot{\nu}(t)$, and in the regime of $t\geq t_{\text{g}}$, we have

From this, we associate the coefficient of the time-dependent term to be equal to the transient change in the spin-down rate, $\Delta\dot{\nu}_{\text{t}}$, which is given by

such that $\Delta\dot{\nu}(t)=\Delta\dot{\nu}_{\text{p}}+\Delta\dot{\nu}_{\text{t}}e^{-\frac{\Delta t}{\tau}}=\Delta\dot{\nu}_{\text{p}}+\Delta\dot{\nu}_{\text{t}}(t)$. There is a slight subtlety in the notation here since we have explicitly used parentheses to represent a time-dependence, e.g. $\Delta\dot{\nu}_{\text{t}}(t)$ depends on time but $\Delta\dot{\nu}_{\text{t}}$ is a constant. The phenomenological glitch model above (or slight variants of it) have been established for a long time, having been used to model the first few glitches of the Crab pulsar (Boynton et al., 1972) and the Vela pulsar (Downs, 1981). Even now, it is still being used in pulsar timing softwares such as TEMPO2 (Edwards et al., 2006).

Our analysis uses publicly available glitch data from the JBCA Glitch Catalogue (Espinoza et al., 2011). Within this catalogue are values for $\frac{\Delta\nu(t_{\text{g}})}{\nu_{0}}$ and $\frac{\Delta\dot{\nu}(t_{\text{g}})}{\dot{\nu}_{0}}$ for each glitch of a given pulsar. Respectively, $\nu_{0}$ and $\dot{\nu}_{0}$ are the measured values of the spin frequency and the time derivative of the spin frequency immediately before the glitch. However, there is an issue using this data since we are not told exactly how much of $\Delta\nu(t_{\text{g}})$ and $\Delta\dot{\nu}(t_{\text{g}})$ is due to a transient part and how much of it is due to a permanent part. We therefore require more data and for us, it will be in the form of the healing parameter, $Q$.

By definition, $Q$ is defined as

where the numerator represents the change in spin frequency which will “heal” at $t\gg t_{\text{g}}$ and the denominator represents the total change in spin frequency at $t=t_{\text{g}}$ which is the sum of permanent and transient parts, as seen in Equation (2). Therefore, a glitch which recovers completely has $Q=1$ and a glitch showing no recovery whatsoever has $Q=0$.

We can substitute Equation (5) into Equation (4) to relate our unknowns ($\Delta\dot{\nu}_{\text{t}}$, $\Delta\nu_{\text{t}}$) to our observables ($Q$, $\Delta\nu(t_{\text{g}})$) to get the relations

$\tau$ makes an appearance in Equation (6) and it can either be treated as known or unknown. It is known when radio observations are frequent enough to acquire $\tau$ directly from the data. However, this is not the case for most glitches as the exponential recovery is sometimes missed, or if too few observations are made during the post-glitch recovery to allow a reliable fit. If it is unknown, then we can use Equation (6) to approximate $\tau$, see Section 7. Throughout this analysis, we will treat $\tau$ as unknown but if it is known in reality, then our model provides an independent value of $\tau$ which can be checked for consistency and guide future transient CW searches.

In this section, we focus on calculating the ellipticity and GW strain obtainable from transient mountains in terms of our observables. It will help at this point if we remind ourselves of the steps of the model: 1) a NS glitches, 2) a transient mountain immediately forms, 3) transient CWs are emitted whilst the transient mountain dissipates at a similar rate as the post-glitch recovery. We do not specify the cause of the glitch (e.g. superfluid unpinning (Anderson & Itoh, 1975), starquakes (Ruderman, 1969)) and we do not attempt to explain the mechanism behind how the mountain forms or dissipates, but for a recent review on possible mechanisms for the formation of mountains, see Sieniawska & Bejger (2019).

It is well-known that a NS mountain will emit CWs at twice the NS’s spin frequency, $f=2\nu$ (e.g. Shapiro & Teukolsky, 1983). A mountain also creates an extra braking torque on the system which we assume will explain the post-glitch recovery. Basic mechanics tells us that the torque, $\mathcal{N}$, is related to power (or in our case GW luminosity) by the equation $L=-\mathcal{N}\Omega$ (for $L>0$) where $\Omega$ is the angular velocity of the NS. It can also be shown (e.g. Shapiro & Teukolsky, 1983) that the GW luminosity due to a mountain is

where $I$ represents the moment of inertia about the rotation axis, $f$ is the GW frequency of the emitted CWs and $\varepsilon(t)$ is the dimensionless equatorial ellipticity of the NS (see Equation (62) for the definition in terms of moment of inertias about the principal axes).

By taking out a factor of $\Omega=2\pi\nu=\pi f$ from Equation (7), multiplying throughout by $-1$ and using $f\approx 2\nu_{0}$, we get the torque due to a NS mountain

Also from mechanics, we know $\mathcal{N}=\dot{I}\Omega+I\dot{\Omega}$ in general. However, as we show in Appendix A, during the post-glitch recovery we can associate the change in torque solely to a change in the spin-down rate caused by a transient mountain which means

where we have ignored the effect of $\Delta\dot{I}$, $\Delta\Omega$ and $\Delta I$ and specialised to transient mountains only. It is worth noting that the internal superfluid is important during glitches and can cause changes to the NS’s moment of inertia (e.g. Link et al., 1999; Andersson et al., 2012). However, our model only concerns the change of the shape of the NS during the post-glitch recovery and omits details about the interior, as the change in shape is what we ascribe to the post-glitch recovery.

Since our model associates the change in torque purely to a transient mountain, we equate the left hand sides of Equations (8) and (9). This allows us to find an expression for the ellipticity for $t\geq t_{\text{g}}$ which is

Note that this means our model will only apply to pulsars rotating fast enough, as slow rotators will give a value of $\varepsilon(t)$ which would be problematically large in terms of the physics. Since $\Delta\dot{\nu}_{\text{t}}(t)$ is in the expression for $\varepsilon(t)$ and $\Delta\dot{\nu}_{\text{t}}(t)$ is unknown, we will set $\Delta\dot{\nu}_{\text{t}}(t)=\Delta\dot{\nu}(t)$ to give us an approximation for the ellipticity from a transient mountain. This means we assume the entire change in the spin-down rate is purely transient, i.e. no permanent mountains are formed which would be associated with $\Delta\dot{\nu}_{\text{p}}$ within the framework of our model. As mentioned later in a footnote of Section 7, this is true for most glitches besides the largest glitches from the Crab where there appears to be a linear relationship between $\Delta\dot{\nu}_{\text{p}}$ and $\Delta\nu$ (Lyne et al., 2015). We can then rearrange in terms of how the JBCA Glitch Catalogue reports glitch values, as well as setting $t=t_{\text{g}}$, to get the ellipticity at the moment of the glitch which is

Now that we have the ellipticity, we can calculate the GW strain, $h_{0}(t)$. In general, a given $\varepsilon(t)$ sources a GW strain of

where $d$ is the distance to the source (e.g. Jaranowski et al., 1998). We can then substitute the ellipticity, Equation (10), into Equation (12) to find the corresponding GW strain at $t\geq t_{\text{g}}$ which is

where we have used $\Delta\dot{\nu}_{\text{t}}(t)=\Delta\dot{\nu}_{\text{t}}e^{-\frac{\Delta t}{\tau}}$. Note that $h_{0}(t)\propto e^{-\frac{\Delta t}{2\tau}}$ and not $h_{0}(t)\propto e^{-\frac{\Delta t}{\tau}}$ as one might originally think. This implies $\tau_{\text{GW}}=2\tau_{\text{radio}}$, see Section 7. This is due to $\Delta\dot{\nu}_{\text{t}}(t)$ being square-rooted in the expression for $\varepsilon(t)$, as per Equation (10). Finally, we can get an approximation for the GW strain at the moment of the glitch by using $\Delta\dot{\nu}_{\text{t}}=\Delta\dot{\nu}(t_{\text{g}})$ which gives

We now proceed on to calculating how much energy is available for the emission of GWs due to the loss of kinetic energy during the post-glitch recovery. We do not include in our calculation any energy which might be liberated during the glitch or required for the formation of transient mountains since these values are uncertain. We begin with the simple expression for the rotational kinetic energy, $E_{\text{rot}}=\frac{1}{2}I\Omega^{2}$. Any instantaneous loss (hence minus sign) of rotational kinetic energy we can write as a luminosity $L(>0)$ by differentiating with respect to time leading to

In particular, we are interested in the GW luminosity achievable due to a mountain which we attribute to having a negative $\Delta\dot{\nu}$. We can also get a change in the luminosity if there is a change in $I$ or $\nu$. However, as demonstrated in Appendices A and B, we are allowed to ignore the contributions due to $\Delta I$ and $\Delta\nu$ as these will be negligible.

The change in luminosity due to the glitch, $\Delta L$, can therefore be written as

where we have put time-dependence into the equation and specialised to transient mountains only. This means we are only tracking the GWs due to transient mountains and not from any permanent mountains which arise due to a permanent change in the spin-down rate (these would give conventional CWs).

We can now integrate Equation (16) across all time to get the GW energy due to a transient mountain and using $\Delta\dot{\nu}_{\text{t}}(t)=\Delta\dot{\nu}_{\text{t}}e^{-\frac{\Delta t}{\tau}}$, we get

and using Equation (6) we can relate the unknowns to our observables giving the final result of

or when put in terms of how the JBCA Glitch Catalogue reports measurements, we get

Note that Equation (18) can be seen as the answer you would get if you were to naïvely find the change in the rotational kinetic energy, but with an additional factor of $Q$ to account for the partial recovery of the glitch.

Next, we relate the GW energy to the optimal¹¹1The SNR is optimal when the matched-filter exactly describes the data. Whenever there is a mismatch, the SNR becomes sub-optimal. See Fig. 1 of Prix et al. (2011). SNR achievable. We follow the steps of Prix et al. (2011) but apply it to our transient mountain model instead of their superfluid excess energy model. Additionally, we account for the effects of having multiple detectors and non-perpendicular interferometer arms which were not considered previously.

For a CW with a time-varying GW amplitude, $h_{0}(t)$, the geometrical-average of the optimal SNR squared, $\langle\rho_{0}^{2}\rangle$, is commonly written as

where the angled brackets represents an average over geometrical parameters which are the right ascension, declination, polarisation angle and inclination of the source. This is what is responsible for the factor of $\frac{4}{25}$ which otherwise would not be there if the SNR was not averaged (Jaranowski et al., 1998; Prix et al., 2011). $S_{\text{n}}(f)$ is the GW detector’s noise power spectral density which is a measure of how much noise the detector picks up at a given GW frequency, $f$. $T_{\text{obs}}$ is the duration for which we observe the source. We eventually set $T_{\text{obs}}\rightarrow\infty$ to ensure we capture the entire transient CW.

Equation (20) is used frequently within the literature though it does not fully capture the effects which apply to us here in this analysis. The ET is a 3rd-generation GW detector which will consist of three identical interferometers, each with a pair of arms with an opening angle of 60°, arranged in a triangle such that two arms share the same side of the triangle (Freise et al., 2009). The fact that we have multiple detectors increases the SNR and an opening angle less than 90° decreases the SNR. Combining both these effects changes Equation (20) into

where $N$ is the number of independent interferometers and $\zeta$ is the opening angle between the arms of the interferometer. The $\sin^{2}\zeta$ factor can also be found in Jaranowski et al. (1998). It is important we account for this as the $S_{\text{n}}(f)$ data for the ET is for a single 90° interferometer²²2The ET-D sensitivity curve was taken from http://www.et-gw.eu/index.php/etsensitivities., even though in reality the ET is made up of three 60° interferometers (Hild et al., 2011).

We now need to evaluate the integral in Equation (21). To do this, we directly substitute for $h_{0}(t)$ using Equation (13) followed by Equation (6) to get

which is the SNR (squared) in terms of the observables within our model.

Alternatively, we can utilise the GW energy we calculated in Section 4. Equation (7) gives the GW luminosity as a function of $\varepsilon(t)$, but $\varepsilon(t)$ is related to $h_{0}(t)$ through Equation (12). We can therefore write the GW luminosity as a function of $h_{0}(t)$ and integrate to get an expression for the GW energy. This gives

We rearrange for the integral and substitute the integral into Equation (21) to get

Finally, we substitute in the GW energy calculated in Equation (18) and we get the same answer as Equation (22). This alternative method is perhaps more powerful than the first as Equations (23) and (24) are completely agnostic with regard to the time-dependence of $h_{0}(t)$ (since we substituted the entire integral). For the same reason, Equation (24) is valid even if we integrate over a short period of time. In other words, it means that we accumulate more SNR the longer the observation is. Equation (24) agrees with Equation (4) in Prix et al. (2011) besides the factor of $N\sin^{2}\zeta$ which we added in here.

To bring familiarity to all this, we can also express our results in terms of a measure of the GW strain used in burst searches, since after all, transient CWs are essentially just a long burst. The measure used is called the root-sum-squared of the GW amplitude, $h_{0,\text{rss}}$, which is defined as

Similarly, we can define the geometrically-averaged $h_{0,\text{rss}}^{2}$ as

such that

where the last equalities in Equations (26) and (27) come from using our transient mountain model and setting $T_{\text{obs}}\rightarrow\infty$. The first equality in Equation (27) can be seen from combining Equation (21) and the definition in Equation (26).

We can check for consistency by substituting our expression for $h_{0}(t_{\text{g}})$ from Equation (12) and then using Equation (6) to get the SNR in terms of observables. After doing that, we find that the right hand side of Equation (27) gives the same answer as Equation (22). In terms of our observables, $\langle h_{0,\text{rss}}^{2}\rangle$ can explicitly be written as

The benefit of calculating $\langle h_{0,\text{rss}}^{2}\rangle$ is that it has the same units as $S_{\text{n}}(f)$ which are units of $\text{Hz}^{-1}$. This then allows us to plot both quantities on the same set of axes allowing a visual comparison between the two terms. In fact, if we plot $\sqrt{\langle h_{0,\text{rss}}^{2}\rangle}$ and $\sqrt{\frac{S_{\text{n}}(f)}{N\sin^{2}\zeta}}$ on the same axes, then the ratio of the two gives exactly the SNR.

In reality, there is some SNR threshold, $\rho_{\text{thres}}$, which we must exceed before confidently accepting a detection. It varies depending on what search method is used, whether the data is stacked coherently and how large the search parameter space is (Walsh et al., 2016; Dreissigacker et al., 2018). For a targeted coherent CW search, this threshold³³3This particular SNR threshold gives a single trial false alarm rate of 1% and a false dismissal rate of 10%. is $\rho_{\text{thres}}\approx 11.4$ (Abbott et al., 2004). Therefore, if we plot $\sqrt{\langle h_{0,\text{rss}}^{2}\rangle}$ and $\sqrt{\frac{S_{\text{n}}(f)\rho_{\text{thres}}^{2}}{N\sin^{2}\zeta}}$ on the same axes, any signal which lies above the modified sensitivity curve will be confidently detectable, with a SNR greater than $11.4$.

We will now take our model and apply it to observed radio data. First we need to select which pulsars to use in our model. Obviously the pulsar must glitch so that constrains us to 190 pulsars⁴⁴4Value taken from http://www.jb.man.ac.uk/pulsar/glitches/gTable.html on 24th August 2020.. However, the most limiting factor requires us to resolve the recovery of the glitch which requires a high cadence of observations, i.e. pulsars which are observed frequently enough to see changes due to the glitch recovery. There are two outstanding candidates which satisfy these constraints and they are the Crab pulsar (B0531+21) and the Vela pulsar (B0833-45). These two pulsars are observed daily at the Jodrell Bank Observatory (Lyne et al., 2015) and at the Mount Pleasant Radio Observatory (Dodson et al., 2007) respectively. Therefore, our analysis will focus on the Crab and Vela pulsars, though the model is applicable to any glitching pulsar where $Q$ can be obtained and is rotating fast enough to prevent unphysically large mountains, see Equation (10).

We used radio data from three main sources: JBCA Glitch Catalogue for $\frac{\Delta\nu(t_{\text{g}})}{\nu_{0}}$ and $\frac{\Delta\dot{\nu}(t_{\text{g}})}{\dot{\nu}_{0}}$ (Espinoza et al., 2011), Crawford & Demiański (2003) for $Q$ and the ATNF Pulsar Database for $\nu_{0}$, $\dot{\nu}_{0}$ and $d$ (Manchester et al., 2005). The value for $\nu_{0}$ needed to be modified to account for the secular spin-down of the NS⁵⁵5We found that if one used the value of $\nu_{0}$ for the Crab from the ATNF Pulsar Catalogue without accounting for secular spin-down, then the GW frequency fell within the 60 Hz mains power line spike in the sensitivity curves of the O2 detectors., but $\dot{\nu}_{0}$ and $d$ were taken as what was reported in the ATNF Pulsar Catalogue. The secular spin-down model we used included the first time derivative of the frequency only, namely

where the subscript “ATNF” represents the value taken from the ATNF Pulsar Catalogue. $t_{\text{ATNF}}$ is the epoch corresponding to when the ATNF values were calculated. This was mid-1991 for the Crab and early-2000 for Vela (Manchester et al., 2005). The left hand side of Equation (29) is what was used for $\nu_{0}$ in all calculations.

$t_{\text{g}}$ was set as the date of the glitch when calculating $E_{\text{GW}}$, $\varepsilon_{\text{approx}}(t_{\text{g}})$ and $h_{0,\text{approx}}(t_{\text{g}})$. However, when calculating SNRs for the different detectors, $t_{\text{g}}=\text{MJD}~{}57856$ for O2 detectors, $t_{\text{g}}=\text{MJD}~{}59761$ for aLIGO at design sensitivity and $t_{\text{g}}=\text{MJD}~{}64693$ for the ET. These correspond to the middle of the O2 run, the middle of 2022 and the middle of the 2030s respectively, which are the approximate times for when these detectors may be operational. This step is required since we have $S_{\text{n}}(2\nu_{0})$ in the denominator of the SNR and $\nu_{0}$ depends on the epoch you observe it.

The values for $\nu_{0}$, $\dot{\nu}_{0}$ and $d$ from the ATNF Pulsar Catalogue are found in Table 1 and the remaining data from the other two data sources can be found in Tables 2 and 3. Any unknown $Q$’s were set to the average $Q$ from existing known values. For the Crab, this was $Q\approx 0.84$ and for Vela, it was $Q\approx 0.17$.

Crawford & Demiański (2003) provides a comprehensive list of $Q$’s compiled from a large collection of literature. Each glitch has a $Q$ value and for some glitches, there are several $Q$’s, each from a different research group. Whenever there was more than one $Q$ for a single glitch, we used the average which is the value we report in Tables 2 and 3.

As for the GW side of the analysis, we will look at the SNRs achievable in the Hanford and Livingston detectors (during the O2 observation run), aLIGO at its design sensitivity and the ET in its ‘D’ configuration⁶⁶6The sensitivity curves of Hanford, Livingston and aLIGO at design sensitivity were taken from https://dcc.ligo.org/LIGO-G1701570/public, https://dcc.ligo.org/LIGO-G1701571/public and https://dcc.ligo.org/LIGO-T1800044/public respectively.. Like the sensitivity curve of the ET, the aLIGO design sensitivity curve is for a single interferometer (Aasi et al., 2015) and since two 90° interferometers make up aLIGO, we set $N_{\text{aLIGO}}=2$ and $\zeta=90\degree$. As established in Section 5, $N_{\text{ET}}=3$ and $\zeta=60\degree$ for the ET. We will treat each of the Hanford and Livingston detectors in O2 as individual interferometers since they have slightly different sensitivity curves, hence, $N=1$ and $\zeta=90\degree$ for both. One can simply multiply the SNR of either Hanford or Livingston by $\sqrt{2}$ to get an approximate combined SNR (or add in quadrature for a more accurate value), but as we will see, this makes little difference in detecting transient CWs outlined in our model for these two detectors.

Together, this data was used to calculate $E_{\text{GW}}$, $\sqrt{\langle\rho_{0}^{2}\rangle}$, $\varepsilon_{\text{approx}}(t_{\text{g}})$, $h_{0,\text{approx}}(t_{\text{g}})$ and $\sqrt{\langle h_{0,\text{rss}}^{2}\rangle}$. The results are shown in Tables 2 and 3.