2. Correlation and spectral properties

We consider a two-component NS with angular velocities $\boldsymbol{\Omega}_t = \left(\Omega_p(t), \Omega_n(t)\right)$ , where one component (labelled p) represents the EM observable ‘normal’ component (protons, electrons, and all constituents strongly coupled to the magnetosphere), and the other (labelled n) corresponds to the neutron superfluid, as is standard in glitch models (see Haskell & Melatos Reference Haskell and Melatos2015; Antonelli et al. Reference Antonelli, Montoli and Pizzochero2022, for a review).

Following Meyers et al. (Reference Meyers, Melatos and O’Neill2021), we model the spin-down of the NS in terms of the stochastic equation

(1) \begin{equation}\dot{\boldsymbol{\Omega}}_t = B \, \boldsymbol{\Omega}_t +\mathbf{A} + M \dot{\mathbf{W}}_t \, ,\end{equation}

where $\dot{\mathbf{W}}_t$ is a vector consisting of two independent standard white noise processes,

(2) \begin{equation}\langle \dot{\mathbf{W}}_t \rangle = 0 \, ,\qquad \quad\langle \dot{W}^i_t \dot{W}^j_s\rangle = \delta_{ij} \delta(t-s) \, ,\end{equation}

and

(3) \begin{equation}B = \begin{pmatrix} -\dfrac{x_n}{\tau} & \dfrac{x_n}{\tau} \\[1em] \dfrac{x_p}{\tau} & -\dfrac{x_p}{\tau} \end{pmatrix} ,\quad M = \begin{pmatrix} \dfrac{\sigma_\infty}{x_p} & -\dfrac{\sigma_\mathcal{T}}{x_p} \\[1em] 0 & \dfrac{\sigma_\mathcal{T}}{x_n} \end{pmatrix} ,\end{equation}

with

(4) \begin{equation}\mathbf{A}= \begin{pmatrix}-\dfrac{\dot{\Omega}}{x_p}\\[1em] 0 \end{pmatrix} .\end{equation}

In the above expressions $x_n$ is the moment of inertia fraction of the superfluid responsible for the internal fluctuating torque, $x_p=1-x_n$ is the fractional moment of inertia of the rest of the star, $\tau$ is the relaxation timescale due to internal friction, while $\sigma_\mathcal{T}$ sets the intensity of the fluctuations in the internal mutual friction torque, and $\sigma_\mathcal{\infty}$ sets the strength of fluctuations in the external braking torque (Antonelli et al. Reference Antonelli, Basu and Haskell2023). Finally, $\dot{\Omega}\gt0$ is the absolute value of the average spin-down rate, measured over a long period.

To study long-term fluctuations due to TN, we are interested in the timing residuals, that is, the deviations from a deterministic spin-down model. We therefore define the angular velocity residuals $\delta\boldsymbol{\Omega}_t=\left(\delta{\Omega}_p(t) , \delta{\Omega}_n(t)\right)$ as:

(5) \begin{equation} \delta\boldsymbol{\Omega}_t = \boldsymbol{\Omega}_t - \langle \boldsymbol{\Omega}_t \rangle \, ,\end{equation}

where the expected value $\langle \boldsymbol{\Omega}_t \rangle$ coincides with the steady-state deterministic drift of $\boldsymbol{\Omega}_t$ due to $\mathbf{A}_t$ , namely

(6) \begin{equation} \langle \boldsymbol{\Omega}_t \rangle =\begin{pmatrix} \Omega-\dot{\Omega} \, t\\ \Omega+\tau\dot{\Omega}/x_p-\dot{\Omega}\,t \end{pmatrix} .\end{equation}

The evolution of the angular velocity residuals is then given by the stochastic equation:

(7) \begin{equation}\delta\dot{\boldsymbol{\Omega}}_t = B \, \delta\boldsymbol{\Omega}_t +M\dot{\mathbf{W}}_t \, .\end{equation}

The above equation provides a simple framework to derive the expected features of TN, like its PSD matrix:

(8) \begin{equation} P_{ij}(\omega) = \left[ \left( i \omega\,\mathbb{I} - B \right)^{-1} M M^\top \left( -i \omega \,\mathbb{I}- B^\top \right)^{-1} \right]_{ij}\, ,\end{equation}

where $i,j=n,p$ (the PSD matrix is Hermitian since the residuals are real). The existence of internal superfluid layers leaves its mark on the PSD of the two components (Antonelli et al. Reference Antonelli, Basu and Haskell2023), obtained by setting $i=j=p$ and $i=j=n$ in the above equation:

(9) \begin{eqnarray}\begin{split} P_{p}(\omega) &= \dfrac{x_p^2 \sigma_\infty^2 + \omega^2 \left(\sigma_\mathcal{T}^2+\sigma_\infty^2\right) \tau^2}{ \omega^2 x_p^2 \left(1+ \omega^2 \tau^2\right)}\propto \frac{\omega^2+\mu_p^2}{\omega^2+\xi^2} \, \omega^{ -\beta} ,\\ P_{n}(\omega) &= \dfrac{x_n^2 \sigma_\infty^2+\omega^2 \sigma_\mathcal{T}^2 \tau^2}{ \omega^2 x_n^2 \left(1+ \omega^2 \tau^2\right)}\propto \frac{\omega^2+\mu_n^2}{\omega^2+\xi^2} \, \omega^{ -\beta} ,\\\xi&=\frac{1}{\tau}\, ,\qquad \quad\mu_n = \frac{x_n \sigma_{\infty} }{\tau \, \sigma_\mathcal{T} } \, ,\\\mu_p &= \frac{x_p \sigma_{\infty} }{\tau \left(\sigma_\mathcal{T}^2+\sigma_{\infty}^2\right)^{1/2} } \, .\end{split}\end{eqnarray}

In the present specific case where white noise is injected, we have that $\beta=2$ . It is possible, however, to directly inject red noise and obtain a higher spectral index $\beta\gt2$ for the asymptotic decay of $P_n$ and $P_p$ . For example, if instead of injecting $\dot{\mathbf{W}}_t$ (white noise) we were to inject $\mathbf{W}_t$ (Wiener noise), we would obtain exactly the two PSDs in (9) but with $\beta=4$ , recovering a more red behaviour that is similar to the one obtained with the alternative stochastic model for TN proposed by Vargas & Melatos (Reference Vargas and Melatos2023). An explicit example is carried out in Appendix A.

Whatever the injected noise model, the main message is that the presence of an internal superfluid component leaves its imprint in the PSD $P_{p}$ of the residuals $\delta\Omega_p$ as a whiterFootnote ^a region $\mu_p\lt\omega\lt\xi$ defined by the two corner frequencies $\mu_p$ and $\xi$ . Note that the corner frequencies are always ordered as $0 \lt \mu_p \lt\xi$ , see equation (9).

Similarly, if $\delta\Omega_n(t)$ were observable in the GW channel, we could cross-validate the presence of the corner frequency $\xi$ in $P_n$ and, possibly, of $\mu_n$ . In this case, however, there is no unique possible ordering of the corner frequencies: if $0\lt\mu_n\lt\xi$ then the intermediate frequency region tends to be whiter ( $P_n\sim \omega^{2-\beta}$ ), if $0\lt\xi\lt\mu_n$ then the intermediate region is redder ( $P_n\sim \omega^{-2-\beta}$ ).

Possible extraction of $\delta\Omega_n(t)$ from CW data would still be valuable even without transitioning to the frequency domain, which would require sufficiently high resolution to distinguish the corner frequencies $\xi$ and $\mu_n$ . An alternative to studying $P_n$ could be to extract information from the correlation between the CW and EM signals, as measured by the parameter $\alpha$ (Jones Reference Jones2004). By imposing the initial condition $\delta\boldsymbol{\Omega}_{t=0}=0$ at a certain arbitrary time $t=0$ , we solve the stochastic system (7) and compute the equal-time correlation matrix

(10) \begin{equation}c^{ij}_{t}=\langle \delta{\Omega}_t^i \delta{\Omega}_t^j \rangle \, ,\end{equation}

where the brackets indicate the average over noise realisations. Explicit calculation gives

(11) \begin{equation} c^{ij}_t = c_0^{ij} + c_1^{ij} e^{-t/\tau} + c_2^{ij} e^{-2t/\tau} \approx c_0^{ij} \quad(\text{for}\,\,\,t\gg\tau)\end{equation}

where $c_1^{ij}$ and $c_2^{ij}$ are constant coefficients and

(12) \begin{equation}\begin{split}& c_0^{pp}= t \, \sigma_\infty^2 -\frac{\tau}{2 x_p^2} \left(\sigma_\mathcal{T}^2 + \sigma_\infty^2 (1+2 x_p +3 x_p^2)\right) \, ,\\& c_0^{nn}= \frac{\tau}{2 x_n^2} \left(\sigma_\mathcal{T}^2 -3 \sigma_\infty^2 x_n^2 \right) \, ,\\& c_0^{np}= t \, \sigma_\infty^2 -\frac{\tau}{ x_n x_p} \left(\sigma_\mathcal{T}^2 + \sigma_\infty^2 x_n(3 x_p-1)\right) \, .\end{split}\end{equation}

This allows us to find the long-term average slope $\alpha$ of the trajectory $\delta{\boldsymbol{\Omega}}_t$ in the plane of the possible residual values $(\delta{\Omega}^n_t, \delta{\Omega}^p_t)$ as

(13) \begin{equation} \alpha \, = \, c^{np} \, / c^{pp} \approx \, c_0^{np}/ c_0^{pp} \quad(\text{for}\,\,\,t\gg\tau) \end{equation}

where the approximation is valid for long times. For $t \gg \tau$ (so that the system loses memory of the arbitrary initial condition $\delta\boldsymbol{\Omega}_{t=0}=0$ ), we find two interesting limits: $\alpha=1$ if there are no fluctuations in the internal torque (i.e. $\sigma_\mathcal{T}=0$ ), and $\alpha=-x_p/x_n$ when there are no fluctuations in the external torque (i.e. $\sigma_\mathcal{\infty}=0$ ).

3. Synthetic timing noise

In order to validate the above analysis and get a sense of the numbers involved, we simulate the stochastic dynamics in (7). In the following, we will adopt the common parametrisation of the relaxation timescale $\tau$ in terms of the dimensionless friction parameter $\mathcal{B}$ ,

(14) \begin{equation} \tau \, = \, x_p/(2\, \Omega \, \mathcal{B}) \, ,\end{equation}

see Montoli et al. (Reference Montoli, Antonelli, Magistrelli and Pizzochero2020) for of how to interpret this body-averaged parameter starting from models that account for local stratification and gradients. Moreover, rather than using the parameters $\sigma_\mathcal{T}$ and $\sigma_{\infty}$ to set the intensity of the internal and external torque fluctuations, we will use two dimensionless parameters $\alpha_\mathcal{T}$ and $\alpha_\infty$ . For a given pulsar of angular velocity $\Omega$ and absolute value of the spin-down rate $\dot{\Omega}$ , they are defined via

(15) \begin{equation}\sigma_\mathcal{T}^2 \,=\, \dfrac{ \alpha_\mathcal{T}^2 \, x_1^2 \,\dot{\Omega}^2 }{2 \, \mathcal{B} \, \Omega} \, ,\qquad \qquad \quad\sigma_\infty^2 \, = \, \alpha_\infty^2 \, \Omega \, \dot{\Omega} \, .\end{equation}

In this way, the dimensionless parameters $0\leq \alpha_\mathcal{T}\lt1$ and $0\leq \alpha_\infty\lt1$ set the relative strength of fluctuations with respect to the deterministic part of the corresponding torque (Antonelli et al. Reference Antonelli, Basu and Haskell2023).

Figure 1. Simulated timing residuals $\delta a_p(t)$ and $\delta a_n(t)$ for the canonical pulsar J2043+2740 (solid and dotted lines, respectively), shown for four noise realisations. Left panel: Internal noise scenario with parameters $x_p=0.8$ , $\mathcal{B}=10^{-9}$ , $\alpha_\mathcal{T}=10^{-1}$ , $\alpha_\infty=10^{-6}$ , yielding a relaxation time $\tau\approx 70\,$ days. Despite the large $\alpha_\mathcal{T}$ , the simulated peak-to-peak amplitude underestimates the observed $\sim2\,$ s variation over $\sim10$ yr reported in Shaw et al. (Reference Shaw2022). Right panel: External noise scenario with $\alpha_\mathcal{T}=10^{-13}$ and $\alpha_\infty=10^{-5}$ , keeping all other parameters the same. This configuration reproduces the observed TN peak-to-peak amplitude.

In view of the analysis outlined in Section 4, the synthetic TN signal extracted from simulations corresponds to the rotational parameters $\Omega$ and $\dot{\Omega}$ of the canonical pulsar J2043+2740 and the millisecond pulsar J1024−0719. In particular, we aim to reproduce the peak-to-peak TN amplitudeFootnote ^b observed in J2043+2740 and J1024 $-$ 0719 by integrating the stochastic system in equation (1) for different choices of physical parameters. This procedure yields the angular velocity residuals $\delta\boldsymbol{\Omega}_t$ , or equivalently, the frequency residuals. The corresponding phase residuals, $\left(\delta \Phi_p(t), \delta \Phi_n(t)\right)$ , and time-of-arrival (TOA) residuals – hereafter referred to simply as timing residuals – $\left(\delta a_p(t), \delta a_n(t)\right)$ , are then computed using equations (76) and (85) of Antonelli et al. (Reference Antonelli, Basu and Haskell2023). Within this scheme, the timing and phase residuals are proportional to each other and carry the same information. Nonetheless, we use both representations, as EM observations of TN often employ some measure of the typical amplitude of the timing residuals $\delta a_p(t)$ to quantify the TN strength, while searches for CW with interferometers generally rely on the phase wandering $\delta\Phi_n(t)$ of the signal (Jones Reference Jones2004; Piccinni Reference Piccinni2022).

3.1 PSR J2043+2740

PSR J2043+2740 is the shortest period object in the sample of 17 canonical pulsars of Shaw et al. (Reference Shaw2022), completing one rotation every $96\,$ ms (i.e. $\Omega\approx 65.4\,$ rad/s, corresponding to $\nu_{GW}\approx20.8\,$ Hz assuming GW emission from a non-precessing triaxial star). This pulsar exhibits a strongly correlated variation in both emission and the spin-down rate $\dot \Omega$ (with a change in $\dot \Omega\approx 8.4\times 10^{-13}$ rad/s $^2$ of about 6%), occurring approximately once every 10 yr (Shaw et al. Reference Shaw2022). This $\sim6\%$ change in $\dot\Omega$ does not visibly affect the timing residuals we compare with and is included in the behaviour reported by Shaw et al. (Reference Shaw2022), which we use as a reference. This change is interpreted as evidence of magnetospheric variability, suggesting an external origin for at least part of the TN, which, within our framework, is expected to be more consistent with simulations driven by a fluctuating external torque.

To test whether our model can reproduce the peak-to-peak TN amplitude observed in this pulsar, we simulate (7) in two limiting cases: one where the source of noise is almost entirely internal, and one where it is almost entirely external.

According to Shaw et al. (Reference Shaw2022), the peak-to-peak TN amplitude reaches approximately $2\,200\,$ ms over $\sim10$ yr. As shown in Figure 1, the internal noise scenario (left panel), with $\alpha_\mathcal{T} = 0.1$ and $\alpha_\infty = 10^{-6}$ , fails to reproduce the observed TN amplitude despite the relatively high level of internal torque fluctuations. In contrast, the external noise scenario (right panel), with $\alpha_\mathcal{T} = 10^{-13}$ and $\alpha_\infty = 10^{-5}$ , successfully matches the observed variation in timing residuals. Hence, in our toy model, even small stochastic variations in the external torque can lead to large-amplitude TN features when the rotational parameters of PSR J2043+2740 are used, in line with the interpretation that the TN is primarily magnetospheric in origin for this pulsar.

To assess the recoverability of the correlation slope $\alpha \approx \delta \Phi_n / \delta \Phi_p$ from synthetic data, we examine the trajectory of the residuals $(\delta\Omega_p(t), \delta\Omega_n(t))$ , $(\delta\nu_p, \delta\nu_n)$ , and $(\delta\Phi_p, \delta\Phi_n)$ for both internal and external noise scenarios. These are shown in Figure 2, which corresponds to the same four noise realisations used in Figure 1. In the internal noise case (left panel), the residual trajectories tend to cluster around a line with slope $\alpha = -x_p/x_n = -4$ , as expected from the conservation of angular momentum in the absence of significant external torque noise. However, even in this limiting case, even small external fluctuations or the simple presence of a long relaxation timescale $\tau$ cause noticeable deviations from the ideal correlation line.

Figure 2. Four examples of typical trajectories in the $(\delta \nu_p,\delta \nu_n)$ and $(\delta \Phi_p,\delta \Phi_n)$ planes for the non-recycled pulsar J2043+2740, corresponding to the same four noise realisations shown in Figure 1. Left panel: Internal noise scenario. Although the noise originates mainly from internal torque, the presence of small external torque fluctuations causes the trajectories to deviate from the black line of slope $\alpha = -x_p/x_n = -4$ (for comparison, the red line has slope $\alpha = 1$ ). Right panel: External noise scenario. In this case, a small internal torque noise component is still present, causing deviations from the red line with slope $\alpha = 1$ (the black line again marks $\alpha = -x_p/x_n = -4$ ).

In the external noise case (right panel), the trajectories follow a line close to $\alpha = 1$ , as expected from (13). Yet, the presence of a small internal noise contribution and the finite relaxation timescale $\tau$ still produces deviations from the ideal linear behaviour. These deviations underscore an important point: even when the dominant noise source is well understood, accurately recovering the correlation slope $\alpha$ from data requires sufficiently long observational baselines to average out the stochastic spread. The relatively wide trajectory envelopes in both scenarios highlight the challenge of inferring the internal dynamics from phase residuals in any single realisation.

Taken together, Figures 1 and 2 demonstrate that although idealised models provide good guidance, the interplay of internal and external noise – as well as finite coupling times – produces non-negligible scatter in realistic datasets. Careful statistical treatment will be required to extract meaningful information from future GW-resolved TN.

3.2 PSR J1024-0719

Regarding MSPs, Perera et al. (Reference Perera2019) analysed the TN of PSR J0437-4715, PSR J1024-0719, and PSR J1939+2134. The one with the strongest TN is J1939+2134 (B1937+21). However, its timing residuals draw an almost perfect sinusoid of period 30 yr, likely due to a planetary companion or precession (Vivekanand Reference Vivekanand2020). For this reason, we focus instead on PSR J1024 $-$ 0719 (Kaplan et al. Reference Kaplan2016), which has the second-largest peak-to-peak amplitude in the TOA residuals among the MSPs studied by Perera et al. (Reference Perera2019).

This pulsar has a rotational period of $P \approx 5.2$ ms (i.e. $\Omega \approx 1\,200$ rad/s), with long-term observations spanning over 18 yr. The timing residuals exhibit a peak-to-peak variation of $\sim10^{-4}$ s, with an RMS of the TOA of about $7.3\,\unicode{x03BC}$ s (Perera et al. Reference Perera2019).

As for the canonical pulsar J2043+2740, we explore whether our two-component model can account for the observed TN amplitude in two contrasting regimes: one where the fluctuations are predominantly internal and one where they are mostly external. Figure 3 shows the simulated TOA residuals $\delta a_p(t)$ and $\delta a_n(t)$ for both scenarios. In the left panel, the internal noise case, we adopt $\alpha_\mathcal{T} = 10^{-1}$ and $\alpha_\infty = 10^{-8}$ , yielding a long relaxation time of approximately $3\,800$ days (these parameter values are chosen as simple powers of ten for illustrative purposes and are not fine-tuned to reproduce the observed amplitude in detail). This long timescale is necessary to maximise the effect of the fluctuations in the internal torque, as it can be seen in the scaling of $P_p$ in (9). This configuration leads to TN amplitudes that systematically overshoot the observed variation over a 10-yr span (a more realistic peak-to-peak amplitude can be achieved by tuning the $\mathcal{B}$ parameter). We deliberately retain this choice of a very long $\tau$ to test what happens when the relaxation timescale is comparable to the observational baseline (recall that (13) is valid in the opposite limit).

Figure 3. Simulated timing residuals $\delta a_p(t)$ and $\delta a_n(t)$ for the millisecond pulsar J1024 $-$ 0719 (solid and dotted lines, respectively), shown for four noise realisations. Left panel: Mostly internal noise scenario with parameters $x_p = 0.8$ , $\mathcal{B} = 10^{-12}$ , $\alpha_\mathcal{T} = 10^{-1}$ , and $\alpha_\infty = 10^{-8}$ , resulting in a relaxation time $\tau \approx 3\,800$ days. This configuration tends to overestimate the observed peak-to-peak timing residual amplitude of $\sim 10^{-4}$ s over $\sim 10$ yr, as reported in Perera et al. (Reference Perera2019). Right panel: Mostly external noise scenario with $\mathcal{B} = 10^{-11}$ , $\alpha_\mathcal{T} = 10^{-8}$ , and $\alpha_\infty = 10^{-5}$ (all other parameters unchanged), which reproduces the observed TN amplitude. This supports an external (e.g. magnetospheric) origin for the residuals in this MSP.

On the other hand, the external noise scenario (right panel), with $\alpha_\mathcal{T} = 10^{-8}$ and $\alpha_\infty = 10^{-5}$ , matches more closely the observed peak-to-peak amplitude of $\sim 10^{-4}$ s: as in the case of the canonical pulsar, small but persistent stochastic variations in the external torque can drive long-term TN consistent with observations. Given the generally low levels of internal superfluid activity expected in MSPs (they tend not to exhibit glitches, and the few detected are of very small amplitude), this result seems to favour a magnetospheric (external) origin of the noise in PSR J1024 $-$ 0719 as well.

To examine whether the correlation slope $\alpha \approx \delta \Phi_n / \delta \Phi_p$ can be extracted from phase-resolved data, we analyse the correlation trajectories in both frequency and phase space. Figure 4 shows the trajectories in the $(\delta \nu_p, \delta \nu_n)$ and $(\delta \Phi_p, \delta \Phi_n)$ planes for the same noise realisations used in Figure 3. In the mostly internal noise scenario (left panel), the trajectories exhibit noticeable scatter about the red reference line of slope $\alpha = 1$ , due to the relatively long relaxation time with respect to the observational baseline.

Figure 4. Trajectories in the $(\delta \nu_p,\delta \nu_n)$ and $(\delta \Phi_p,\delta \Phi_n)$ planes for the millisecond pulsar J1024 $-$ 0719, corresponding to the same typical noise realisations shown in Figure 3. Left panel: Mostly internal noise scenario. Although the internal torque dominates the fluctuations, a small amount of external torque noise and the finite relaxation time $\tau$ lead to deviations from the red line of slope $\alpha = 1$ (the black line marks $\alpha = -x_p/x_n = -4$ for comparison). Right panel: Mostly external noise scenario. Despite the dominance of external fluctuations, the residual internal noise and finite $\tau$ still cause slight deviations from the red slope line $\alpha = 1$ . In both panels, the black line serves as a reference for $\alpha = -x_p/x_n = -4$ .

In the mostly external noise case (right panel), the trajectories are more tightly clustered around the red line ( $\alpha = 1$ ), reflecting a nearly coherent response of both components to the common external forcing. However, slight deviations remain due to the residual presence of internal noise and a finite relaxation timescale $\tau$ . As in the canonical case, the black reference line with slope $\alpha = -x_p/x_n = -4$ is included for comparison, though it plays a minor role in the external noise regime.

Overall, the behaviour observed in Figures 3 and 4 suggests that even in MSPs with high coupling and low glitch activity, measurable TN can arise from small stochastic fluctuations in the external torque. Moreover, while correlation slopes can, in principle, be recovered from residual trajectories, the accuracy of such measurements will depend critically on two ratios: how the relaxation timescale compares to the observation baseline and the amplitude ratio of internal and external torque fluctuations.

4. GW detection prospects

To determine whether joint EM and GW observations can be used to study the nature of the torques acting on a NS and the coupling between its components, it is crucial to assess whether the phase fluctuations discussed above are detectable, at least with the next generation of GW detectors.

To do this, we will consider the sensitivity of the planned next-generation ground-based GW detector, the ET, considering in particular the configuration with two parallel 15 km long L-shaped detectors, as described in Branchesi et al. (Reference Branchesi2023). The detectability of the signal depends on the exact GW emission mechanisms, the distance, and the frequency of the source. However, it has been shown (Abac et al. Reference Abac2025) that for the currently observed galactic NSs, several hundred would be observable after one year by ET, assuming emission due to a quadrupolar mountain, with a characteristic amplitude $h_0$ of the form:

(16) \begin{equation} h_0 = 10^{-25}\left(\frac{10\mbox{ kpc}}{d}\right)\bigg(\frac{\epsilon}{10^{-6}}\bigg)\bigg(\frac{\nu_s}{500\,\mbox{Hz}}\bigg) I_{45} \, ,\end{equation}

where d is the distance of the source, $\epsilon$ the ellipticity, $I_{45}$ the moment of inertia of the star in units of $10^{45}$ g cm $^2$ , and $\nu_s$ its spin frequency (we recall that the GW frequency $\nu_{GW}=2\nu_s$ in this quadrupolar mountain scenario).

If the source is detectable by ET, the error in measuring the phase $\delta\Phi_{GW}$ will be (Jones Reference Jones2004):

(17) \begin{equation}\delta\Phi_{GW} = \frac{N_{sim}}{\rho_d} \, ,\end{equation}

where $N_{sim}$ is a factor of order unity that depends on the number of parameters involved in the search. For our estimates, we use the conservative value of $N_{sim}=3$ (see Jones Reference Jones2004 and references therein for a discussion of this point) and is $\rho_d$ the signal-to-noise ratio with which the signal is detected, calculated as (Watts et al. Reference Watts, Krishnan, Bildsten and Schutz2008)

(18) \begin{equation}\rho_d=h_0\sqrt{\frac{N_d T_\textrm{obs}}{S_n(\nu)}} \, ,\end{equation}

where $N_d$ is the number of detectors $T_{ \rm obs}$ is the observation length, and $S_n(\nu)$ the noise spectral density of the detector at frequency $\nu$ .

To measure a phase difference greater than the minimum in (17) it is, of course, necessary for the source to be loud enough for detection, with louder sources allowing for more precise measurements of the TN. We set the limit for detection at $\rho_d \gt 11.4$ (Watts et al. Reference Watts, Krishnan, Bildsten and Schutz2008) and consider two cases: first a 1 yr coherent integration, and a $T_\textrm{obs}=$ 10 yr total observation, summed incoherently over $T_\textrm{coh}=$ 1 yr stretches of data, for which (Branchesi et al. Reference Branchesi2023):

(19) \begin{equation}\rho_d\approx h_0\sqrt{\frac{N_d T_\textrm{coh}}{S_n}} \,\left( \frac{T_{ \rm obs}}{T_\textrm{coh}} \right)^{1/4} .\end{equation}

In the following, we will consider the 2L configuration for ET (Branchesi et al. Reference Branchesi2023) and therefore take $N_d=2$ , and one has $\rho_d^{10\,\text{yr}}\approx1.78\rho_d^{1\,\text{yr}} $ . With the above definitions, we can calculate the minimum phase wandering $\delta\Phi_{GW}$ that could be detected in a one-year observation, which is presented in Figure 5 as a function of frequency. We plot our limits for a fixed reference value of $\epsilon=10^{-6}$ for all pulsars, for a range of distances between $d=0.1$ kpc and $d=5$ kpc, that would thus roughly encompass most observed pulsars. We see that for most galactic NSs, with a GW frequency around the KHz range, we would be able to measure, at best, phase wandering of the order of $\delta\Phi_{GW} \gtrsim 10^{-4}$ rad.

Figure 5. The minimum $\delta\Phi_{GW}$ detectable by ET with pulsar observations. The shaded region between the blue lines is for a 1-yr coherent integration of a signal for sources with $\epsilon=10^{-6}$ , with a distance between 0.1 and 5 kpc. The crosses and diamonds are for the observed pulsars described in the text, considering both a 1-yr coherent integration (crosses) and 10-yr in-coherent integration (diamonds). The horizontal line represents the limit $\delta\Phi_{GW}=N/11.4 \approx 0.18$ , which gives us an estimate of detectability. Pulsars above this line would not be detectable. Note, however, that pulsars below the line may very well have an intrinsic timing noise that is larger than the minimum detectable plotted here.

We also show our limit for detection, which for our cutoff at signal-to-noise ratio $\rho_d=11.4$ , gives us $\delta\Phi_{GW} \lt N/11.4\approx 0.18$ rad. This means that pulsars above this line would not be detectable with our assumptions on $\epsilon$ . It does, however, indicate that for all detectable pulsars we would be able to measure a minimum $\delta\Phi\lt1\,$ rad, therefore below the value that would cause our templates to drift out of phase with the real signal and hinder detection. We caution the reader, however, that we are simply estimating the smallest phase shift that could be measured. The shift in phase due to a specific pulsar’s TN will depend on the strength of said noise, and as we have seen, it can be significantly larger.

In Figure 5 we also show the minimum $\delta\Phi_{GW}$ that would be measurable by ET for a subset of observed pulsars, which are known to exhibit strong TN (Lyne Reference Lyne and van Leeuwen2013): 17 non-recycled radio pulsars analysed by Shaw et al. (Reference Shaw2022), Lyne et al. (Reference Lyne, Hobbs, Kramer, Stairs and Stappers2010), and 3 ms pulsars that exhibit stronger TN (Perera et al. Reference Perera2019), namely PSR J0437-4715, PSR J1024-0719, and PSR J1939+2134: (using equation 17 and 18)

(20) \begin{equation} \delta\Phi_{GW} = \frac{2}{h_0 (d,\epsilon,\nu_s)} \sqrt{\frac{S_n(2 \nu_s)}{2 \, T_{obs}}} \, ,\end{equation}

where $I=10^{45}\,$ g cm $^2$ is assumed for all objects.

Note that for the 17 standard radio pulsars, we assume a physically motivated ellipticity of $\epsilon = 10^{-6}$ for our analysis and use their measured distance. However, for the 3 ms pulsars, this value exceeds their spin-down upper limit, $\epsilon_{sd}$ , which represents the ellipticity required for the observed spin-down to be entirely due to GW emission

(21) \begin{equation} \epsilon_{sd}= 1.9\times 10^{-5} \! \left(\frac{100\,\mbox{Hz}}{\nu_s}\right)^{\frac{5}{2}} \!\! \left(\frac{\dot{\nu}_s}{10^{-10}\,\mbox{Hz/s}}\right)^{\frac{1}{2}} \!\! I_{45}^{-\frac{1}{2}} .\end{equation}

For these 3 systems, we therefore assume GW emission occurs at their spin-down limit. Note that this smaller value of $\epsilon$ is the reason that the estimated minimum value of $\delta\Phi_{GW}$ that could be measured is larger than those in the shaded band. We see that, despite the smaller ellipticities, the MSPs in our sample still represent the best targets for detection with ET, and would allow us to measure variations in phase due to TN of the order of $\delta\Phi_{GW}\approx 10^{-2}$ rad. On the other hand, for our assumed mountain size of $\epsilon=10^{-6}$ , only the fastest ‘regular’ pulsars would be loud enough for detection (and therefore allow to attempt a TN measurement).

In fact, not only may the TN be detectable, but it may also be strong enough, in both of the pulsars we consider, to hinder detection, causing the signal to drift out of phase with our model by several tens of radians in the case of the MSP PSR J1024-0719, for both internal and external noise models. This points towards the fact that TN must be accounted for in searches for these signals with ET, and mitigation strategies must be considered (Ashton, Jones, & Prix Reference Ashton, Jones and Prix2015).