#  Stochastic processes for pulsar timing noise: fluctuations in the internal
and external torques

Marco Antonelli1, Avishek Basu2, Brynmor Haskell3
1CNRS/IN2P3, Laboratoire de Physique Corpusculaire, 14050 Caen, France
2Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy,
University of Manchester, M13 9PL, Manchester, United Kingdom
3Nicolaus Copernicus Astronomical Center, Polish Academy of Sciences, ul.
Bartycka 18, Warszawa, Poland
E-mail<EMAIL_ADDRESS>(M. A<EMAIL_ADDRESS>(A.B<EMAIL_ADDRESS>(B.H.)

###### Abstract

Young pulsars deviate from a perfectly regular spin-down by two non-
deterministic phenomena: impulsive glitches and timing noise. Both phenomena
are interesting per se, and may provide insights into the superfluid
properties of neutron stars, but they also act as a barrier to high-precision
pulsar timing and gravitational wave experiments. We study a minimal
stochastic model to describe the spin-down of a multicomponent neutron star,
with fluctuations in both the internal and external torques. The power
spectral density and timing noise strength of this kind of model can be
obtained analytically, and compared with known results from pulsar timing
observational campaigns. In particular, the presence of flat regions of the
power spectral density can be interpreted as a signature of the presence of
internal superfluid components. We also derive the expected scaling of the
timing noise strength with the pulsar’s rotational parameters (or
characteristic age). Therefore, the present framework offers a theoretical
guideline to interpret the observed features of timing noise in both single
pulsars and across the pulsar population.

##  1 Introduction

The continuous monitoring of rotation powered radio pulsars in a timing
campaign reveals both sudden glitches (impulsive spin-up events), as well as a
long-timescale smooth wandering of their rotation period known as timing noise
(D’Alessandro, 1996). Timing noise is a type of rotational irregularity
observed in several pulsars, where the measured pulse time of arrival (TOA)
wander stochastically about the expected TOA obtained from best-fit timing
model via fitted ephemeris (Boynton et al., 1972; Cordes & Downs, 1985). It is
often characterised as a random walk in the rotational phase, angular
velocity, or torque (Groth, 1975); its Fourier spectrum is always red,
implying a process autocorrelated on a time-scale of a few months to years
(Cordes & Helfand, 1980; van Haasteren & Levin, 2013; Liu et al., 2019),
making it a distinct phenomenon from the so-called “jitter” noise caused by
variations in pulse shape at the single-pulse level (Taylor et al., 1975;
Rickett, 1975; Lam et al., 2019).

Timing noise is likely to be a genuine feature of the rotational dynamics of a
pulsar (Hobbs et al., 2010). Hence, this phenomenon has been attributed to
various possible mechanisms at work in the neutron star interior or within its
surroundings like unresolved microglitches (Cordes & Downs, 1985; Janssen &
Stappers, 2006), recovery from unseen glitches (Alpar et al., 1986; Johnston &
Galloway, 1999), changes in the magnetosphere (Lyne et al., 2010; Shaw et al.,
2022), variable coupling between the crust and liquid interior (Alpar et al.,
1986; Jones, 1990), turbulence (Melatos & Link, 2014) and stochastic
variations in the star’s shape and fluctuations in the spin-down torque, see
(D’Alessandro, 1996) and references therein.

In this paper we set a up a framework to study how fluctuations in both the
_internal_ and _external_ torques in a rotating neutron star can affect the
spin wandering of the observed phase. To do this, we rely on the natural
extension of the seminal model proposed by Baym et al. (1969) for the
rotational dynamics of a two-component neutron star (see Montoli et al. 2020
for a discussion of its three-component extension): we complement it by
introducing additive noise in the braking torque and in the internal friction
between the components.

The minimal model of Baym and collaborators proved to be an invaluable tool to
interpret glitches, to the point that both linear (Pizzochero et al., 2020;
Sourie & Chamel, 2020; Montoli et al., 2020) and nonlinear (Gügercinoǧlu &
Alpar, 2017; Celora et al., 2020) extensions have been developed and widely
used to study the spin-up of a glitch and its subsequent recovery, including
the general relativistic corrections (Sourie et al., 2017; Antonelli et al.,
2018; Gavassino et al., 2020). In particular, the signature of an overshoot in
the phase residuals recorded during spin-up phase of a glitch in the Vela
pulsar (Ashton et al., 2019; Montoli et al., 2020) forces us to extend the
original two-component model of Baym et al. (1969) to account for more degrees
of freedom in the interior of a pulsar. This can be done in two ways, by
considering that the superfluid can develop continuous differential rotation
(Haskell et al., 2012; Antonelli & Pizzochero, 2017; Graber et al., 2018), or
by adding a third rigid component (Pizzochero et al., 2020). Given this
observational signature, we will consider a generic number of internal
components.

This kind of stochastic modelling can serve as a theoretical tool to study and
derive the expected features of timing noise, with particular attention to its
spectral properties (like the presence of corner frequencies that indicate a
spectral turnover), in the hypothesis that the origin of noise is intrinsic to
the pulsar. In fact, whatever the origin of noise, the presence of loosely
coupled superfluid layers may leave its imprint on the spectral properties of
the detected noise. Hence, we analytically derive the scaling of the timing
noise strength parameter (the root mean square of the time residuals, see e.g.
Shannon & Cordes 2010) with the rotational parameters of the pulsar and the
physical input of the stochastic model, like the moments of inertia and the
superfluid mutual friction parameters.

A similar stochastic approach has also been recently studied by Meyers et al.
(2021a); Meyers et al. (2021b) in the context of the possible detection of
continuous gravitational wave emission. The main difference with our work is
that they consider two external torques acting on two components, while our
theoretical setting is for a generic number $m+1$ of internal components
driven by one external torque and $m$ internal ones (further extension of the
framework is given in App. A). Clear distinction between the nature of the
torques makes the algebraic results more transparent and allows us to select
only the physical results for the timing noise spectrum.

The paper is organised as follows. In Sec. 2 we describe the general
construction of the minimal stochastic extension of the seminal rigid model of
Baym et al. (1969), for $m$ superfluid components (the special cases $m=1,2$
are solved explicitly). In Sec. 3 we discuss a prescription to set the
strength of the stochastic fluctuations in the internal and external torques.
Sec. 4 is devoted to deriving the spectral properties of the timing noise
expected for the models with $m=1,2$. We show how to extract the timing noise
strength $\sigma$ of Shannon & Cordes (2010) from simulations in Sec. 5: the
scaling of the noise timing noise strength with the rotational parameters of a
pulsar (see, e.g. Parthasarathy et al. 2019) is also discussed. Finally, in
Sec. 6, we perform numerical simulations with the double intent of exploring
the role of the various model parameters and validating the correctness of our
analytical results. The specific two-component model of Meyers et al. (2021a)
is discussed within our formalism in App. A. App. B is devoted to the
interesting limit in which our model reduces to an Ornstein–Uhlenbeck process.

_Notation:_ Throughout the paper, the instantaneous angular velocity will be
indicated as $\Omega_{p}(t)$. We adopt the convention that $\Omega>0$ and
$\dot{\Omega}>0$ are constant values for the _measured_ angular velocity and
the absolute value of the secular spin-down rate of a given pulsar:
$\Omega\approx\Omega_{p}(t)$ and $\dot{\Omega}\approx|\dot{\Omega}_{p}(t)|$.
The parameter $|\dot{\Omega}_{\infty}|>0$ always indicates the strength of the
spin-down torque normalised over the pulsar’s total angular momentum.
Practically speaking, $\dot{\Omega}$ is an estimate - obtained from long
pulsar timing observations - of the unknown torque strength
$|\dot{\Omega}_{\infty}|$. Hence, to keep notation light, we use
$\dot{\Omega}$ in place of $|\dot{\Omega}_{\infty}|$ in expressions where the
exact value of the spin-down torque is unimportant.

##  2 Multi-component systems with additive noise

Building on the seminal work of Baym et al. (1969), we consider a minimal
model for the evolution of $m+1$ rigid internal components with angular
velocities

$\mathbf{\Omega}_{t}=\left(\Omega_{p}(t),\Omega_{1}(t),...,\Omega_{m}(t)\right)\,.$
(1)

We identify its first component $\Omega_{p}(t)$ with the observable angular
velocity, which is typically identified with the rotation of normal matter in
the star and the magnetosphere (Heintzmann et al., 1973; Easson, 1979), where
the signal originates. The other $m$ variables refer to different internal
superfluid components (e.g., $m=1$ for the original model of Baym and
collaborators, $m=2$ for the Vela’s spin-up model described in Pizzochero et
al. 2020): in principle, they may be different superfluid species coexisting
in the same region, as well as the same kind of superfluid extending in $m$
spatially non-overlapping layers 111 The general setting presented here is
agnostic about the actual disposition of the superfluid components inside the
star, the only requirement being the constraints (2) and (5). However, in Sec.
2.3 we explicitly consider the particular model of Montoli et al. (2020),
where there are two distinct superfluid regions. .

Each of the $m+1$ components contributes to a certain fraction
$\mathbf{x}=(x_{p},x_{1},...,x_{m})$ of the total moment of inertia. We do not
consider the effect of possible starquakes, implying that both the total
moment of inertia and the moment of inertia ratios $\mathbf{x}$ are constant.
Therefore, only $m$ out of the $m+1$ values of $\mathbf{x}$ are independent,

$\sum_{i=1...m}x_{i}=1-x_{p}\,.$ (2)

The ratio $x_{p}$ includes all components, including some of the superfluid
components (if any), which are strongly coupled with the normal matter (non-
superfluid component) of the star.

Regardless of the nature of the internal torques, the total angular
momentum222 We will use the terms “angular momentum” and “torque” even if they
are divided by the total moment of inertia of the star.
$L(t)=\mathbf{x}^{\top}\cdot\mathbf{\Omega}$ is affected only by the
continuous action of the external torque,

$\dot{L}\,=\,\mathbf{x}^{\top}\cdot\dot{\mathbf{\Omega}}\,=\,\mathcal{T}_{\infty}+\eta_{t}^{\infty}\,,$
(3)

where $\mathcal{T}_{\infty}$ is the deterministic part of the total braking
torque (arising from electromagnetic and, possibly, gravitational wave
emission), while $\eta_{t}^{\infty}$ is its fluctuating part. Given the above
setting, we model the pulsar’s evolution in terms of a Itô equation of the
kind

$d\mathbf{\Omega}_{t}=(B\,\mathbf{\Omega}_{t}+\mathbf{A}_{t})dt+Md{\mathbf{W}_{t}}\,,$
(4)

that is the natural stochastic generalisation of the multicomponent models for
glitches with linear coupling between the components (Baym et al., 1969;
Montoli et al., 2020). Here, $B$ and $M$ are $(m+1)\times(m+1)$ real matrices
and $\mathbf{A}_{t}$ is a (possibly time dependent) real vector describing the
braking torque. The term $d{\mathbf{W}_{t}}$ is the increment of a standard
$(m+1)$-dimensional Wiener process, where all the components of
$\mathbf{W}_{t}$ are independent identically distributed (i.i.d.) normal
random variables of zero mean and variance $t$, which allows us to model the
fluctuating part of the external torque in (3), as well as possible
fluctuations of the other $m$ internal torques (for a total of $m+1$
independent noise sources). In fact, the mutual transfer of angular momentum
between each of the $m$ superfluid components and the observable one may
undergo stochastic fluctuations because of the intrinsic noisiness in the
vortex-mediated processes (Alpar et al., 1986) or internal turbulence (Melatos
& Link, 2014).

Equation (4) is very general, but there are physical constraints it should
satisfy: the first is that it must be consistent with (3), the second is that
there should be no internal friction between components that have the same
velocity. Therefore, the vector $\mathbf{\Omega}_{t}=(\Omega\,,...,\Omega)$
must be a right null eigenvector of $B$,

$B(\Omega\,,...,\Omega)=0\,,$ (5)

implying that $B$ is not invertible (at corotation, there can not be any
deterministic transfer of angular momentum between the components due to
friction).

The Itô process in (4) can also be written in the equivalent Langevin form,

$\begin{split}&\dot{\mathbf{\Omega}}_{t}=B\,\mathbf{\Omega}_{t}+\mathbf{A}_{t}+M\dot{\mathbf{W}}_{t}\,,\\\
&\langle\dot{\mathbf{W}}_{t}\rangle=0\\\
&\langle\dot{W}^{i}_{t}\dot{W}^{j}_{s}\rangle=\delta_{ij}\delta(t-s)\end{split}$
(6)

where the white noise $\dot{\mathbf{W}}_{t}$ is the formal time derivative of
the standard Wiener process, $d\mathbf{W}_{t}=\dot{\mathbf{W}}_{t}\,dt$, see
e.g. Gardiner (1994). Clearly, the system in (6) must be consistent with the
fundamental requirement (3), so that $\mathbf{A}_{t}$, $B$ and $M$ can not be
set independently, but they should satisfy the constraint

$\mathbf{x}^{\top}\cdot(B\,\mathbf{\Omega}_{t}+\mathbf{A}_{t}+M\dot{\mathbf{W}}_{t})=\mathcal{T}_{\infty}+\eta_{t}^{\infty}\,.$
(7)

The above equation defines a physical property of the coupling matrix $B$:
internal torques should satisfy the action-reaction law. In fact, to guarantee
that (7) is consistent with (3), $\mathbf{x}$ must be a null left eigenvector
of $B$,

$B^{\top}\mathbf{x}=0\qquad\quad\text{(action-reaction property)}\,,$ (8)

and, on the other hand, the total external torque must be given by

$\mathbf{x}^{\top}\cdot\mathbf{A}_{t}=\mathcal{T}_{\infty}\qquad\qquad\mathbf{x}^{\top}M\dot{\mathbf{W}}_{t}=\eta_{t}^{\infty}\,.$
(9)

The solution to (4) is not formally different from the usual solution for an
ordinary differential equation of the same form,

$\mathbf{\Omega}_{t}=e^{Bt}\mathbf{\Omega}_{0}+\int_{0}^{t}\\!\\!e^{B(t-z)}\mathbf{A}_{z}\,dz+\int_{0}^{t}\\!\\!e^{B(t-z)}M\,d\mathbf{W}_{z}\,.$
(10)

For each $t$, the stochastic integral above defines a random variable of zero
average, so that the expected value $\langle\mathbf{\Omega}_{t}\rangle$
coincides with the deterministic drift of $\mathbf{\Omega}_{t}$,

$\langle\mathbf{\Omega}_{t}\rangle=e^{Bt}\langle\mathbf{\Omega}_{0}\rangle+\int_{0}^{t}\\!\\!e^{B(t-z)}\mathbf{A}_{z}\,dz\,.$
(11)

Therefore, by imposing a deterministic initial condition
$\mathbf{\Omega}_{0}=\langle\mathbf{\Omega}_{0}\rangle$, we can define the
angular velocity residual as

$\delta\mathbf{\Omega}_{t}=\mathbf{\Omega}_{t}-\langle\mathbf{\Omega}_{t}\rangle=\int_{0}^{t}\\!\\!e^{B(t-z)}M\,d\mathbf{W}_{z}\,,$
(12)

that is the solution of the Langevin equation

$\delta\dot{\mathbf{\Omega}}_{t}=B\,\delta\mathbf{\Omega}_{t}+M\dot{\mathbf{W}}_{t}\,.$
(13)

We are interested in the statistical properties of the residuals
$\delta\mathbf{\Omega}_{t}$, in particular of its first observable component
$\delta\Omega^{p}_{t}$. Given the solution (12), the average and the
correlation matrix read333 Since we imposed a deterministic initial condition,
we have that $\delta\mathbf{\Omega}_{0}=0$. If this hypothesis is dropped,
then the results in (14) gain an extra term.

$\begin{split}&\langle\delta{\Omega}^{i}_{t}\rangle=0\qquad\forall\,i=p,1,...,m\\\
&\langle\delta\mathbf{\Omega}_{t}\delta\mathbf{\Omega}_{u}^{\top}\rangle=\int_{0}^{\min(t,u)}\\!\\!dz\,\,e^{B(t-z)}MM^{\top}e^{B^{\top}(u-z)}\,.\end{split}$
(14)

Since the correlation matrix is a function of both times $t$ and $u$, the
process is not wide-sense stationary: for $u<t$ and $t=u+\Delta t$, the matrix
$\langle\delta\mathbf{\Omega}_{u+\Delta
t}\delta\mathbf{\Omega}_{u}^{\top}\rangle$ depends on $u$. Therefore, to
calculate the Power Spectral Density (PSD) of the process we can _not_ rely on
the Wiener–Khinchin theorem. However, it is still possible to perform the
Fourier transform on the Langevin equation and use the frequency
representation of the white-noise (Rice, 1944). First, we introduce the
Fourier representation

$\delta\mathbf{\Omega}_{t}=\int_{-\infty}^{\infty}\\!e^{i\omega
t}\delta\mathbf{\Omega}_{\omega}\,\dfrac{d\omega}{2\pi}$ (15)

and define the matrix

$P(\omega,\omega^{\prime})=\langle\delta\mathbf{\Omega}_{\omega}\,\delta\mathbf{\Omega}^{\dagger}_{\omega^{\prime}}\rangle\,,$
(16)

where $\delta\mathbf{\Omega}^{\dagger}_{\omega}$ is the transpose complex
conjugate of the vector $\delta\mathbf{\Omega}_{\omega}$. In particular, we
are interested in the $i=j=p$ component for $\omega=\omega^{\prime}$, that is
the PSD of the signal associated to the observable component:

$P_{p}(\omega)\propto\langle\,|\delta\Omega_{p}(\omega)|^{2}\,\rangle\,.$ (17)

From (13) we have

$\delta\mathbf{\Omega}_{\omega}=(i\omega\,\mathbb{I}-B)^{-1}M\,\dot{\mathbf{W}}_{\omega}\,,$
(18)

that, together with the basic properties of the Wiener process,

$\begin{split}&\langle\dot{W}^{i}_{\omega}\rangle=\int_{-\infty}^{\infty}\\!dt\,e^{-i\omega
t}\langle\dot{W}^{i}_{t}\rangle=0\qquad\forall\,i=p,1,...,m\\\
&\langle\dot{\mathbf{W}}_{\omega}\dot{\mathbf{W}}_{\omega^{\prime}}^{\dagger}\rangle=2\pi\,\delta(\omega-\omega^{\prime})\,\mathbb{I}\,,\end{split}$
(19)

can be used to calculate the matrix in (16):

$P(\omega,\omega^{\prime})=2\pi\,(i\omega\,\mathbb{I}-B)^{-1}MM^{\top}(-i\omega^{\prime}\,\mathbb{I}-B^{\top})^{-1}\delta(\omega^{\prime}-\omega)\,.$
(20)

Therefore, $P_{p}(\omega)$ will contain a divergent Dirac delta factor that
is, however, just a consequence of introducing the (formal) Fourier
representation (15) for a non-normalizable signal (Priestley, 1965): following
Gardiner (1994), we define the PSD for the observable component as

$P_{p}(\omega)=\left[(i\omega\,\mathbb{I}-B)^{-1}MM^{\top}(-i\omega\,\mathbb{I}-B^{\top})^{-1}\right]_{pp}\,,$
(21)

where the subscript $pp$ means that only the $i=j=p$ component of the matrix
defined within the square brackets is taken. Note that $P_{p}(\omega)$ has the
physical dimensions expected for a PSD, namely that of a power (i.e., the
squared amplitude of the signal, $|\delta\Omega^{p}_{t}|^{2}$) per unit
angular frequency $\omega$. Since both $\delta\Omega^{p}_{t}$ and $\omega$
have the dimensions of rad/s, also the PSD in (21) is measured in rad/s (a
concrete example will be given in Fig. 6).

###  2.1 Evolution of the total angular momentum and rigid rotor $(m=0)$

Before moving to more complex cases, let us briefly discuss the evolution of
the total angular momentum residue $\delta L_{t}$ for $m>0$. Its evolution is
given by

$\delta
L_{t}=\mathbf{x}^{\top}\cdot\delta\mathbf{\Omega}_{t}\qquad\quad\delta\dot{L}_{t}=\eta_{\infty}\,,$
(22)

meaning that it is unaffected by the possible presence of internal torques, as
it should be. Its associated power spectral density
$P_{L}(\omega)=\langle|\delta L_{\omega}|^{2}\rangle$ can be found by
considering that

$\delta
L_{\omega}=\mathbf{x}^{\top}\cdot\delta\mathbf{\Omega}_{\omega}\qquad\qquad
i\omega\,\delta
L_{\omega}=\eta_{\omega}=\mathbf{x}^{\top}M\dot{\mathbf{W}}_{\omega}\,.$ (23)

Therefore, we recover the usual result for a pure Wiener process444 The PSD in
(24) can also be obtained from the matrix in (20) after double contraction
with $\mathbf{x}$: thanks to the action-reaction property (8) we have that
$(\omega\,\mathbb{I}-B^{\top})^{-1}\,\mathbf{x}=\mathbf{x}/\omega$. ,

$P_{L}(\omega)=\langle|\delta
L_{\omega}|^{2}\rangle=\mathbf{x}^{\top}(i\omega)^{-1}MM^{\top}(-i\omega)^{-1}\mathbf{x}=\frac{\sigma_{\infty}^{2}}{\omega^{2}}\,,$
(24)

where we have defined the positive parameter $\sigma_{\infty}^{2}$ as

$\sigma_{\infty}^{2}\,=\mathbf{x}^{\top}MM^{\top}\mathbf{x}\,.$ (25)

In the limit in which the star rotates rigidly (i.e, $m=0$), we have that
$L=\Omega_{p}$, $B=0$ and $M=\sigma_{\infty}$. Therefore, the PSD for the
observable component in (21) boils down to the one in (24):

$P_{p}(\omega)=\frac{\sigma_{\infty}^{2}}{\omega^{2}}\qquad\quad(\text{if}\quad
m=0).$ (26)

Therefore, in the limit of rigid rotator, we recover the usual red-noise PSD
of a pure Wiener process.

###  2.2 Two-component model $(m=1)$

The simplest $m=1$ model reads

$\begin{split}&x_{p}\dot{\Omega}_{p}=-\mathcal{T}+\mathcal{T}_{\infty}\\\
&x_{1}\dot{\Omega}_{1}=\mathcal{T}\,,\end{split}$ (27)

where $\Omega_{p}$ and $\Omega_{1}$ are the angular velocities of the observed
and superfluid components, $x_{p}$ and $x_{1}=1-x_{p}$ are the moment of
inertia ratios of each component, see (2), and $\mathcal{T}$ is the internal
torque. In the following, for simplicity we will set
$\mathcal{T}_{\infty}=-|\dot{\Omega}_{\infty}|$, a constant value that sets
the secular spin-down rate. Due to the presence of fluctuations, the intrinsic
value of the external torque $|\dot{\Omega}_{\infty}|$ may slightly differ
from the absolute value of the spin-down rate $\dot{\Omega}$ extracted from a
fit to the observed signal555 We recall that $\dot{\Omega}>0$ is a constant
benchmark value for the observed absolute value of the secular spin down rate,
$\dot{\Omega}\approx|\dot{\Omega}_{\infty}|$. .

Note that we are assuming that the physical mechanism responsible for the spin
down acts directly only on the observable component. This is reasonable as
long as this mechanism does not directly cause an outward motion of vortex
lines. The internal torque $\mathcal{T}$ in (27) is phenomenologically
modelled in terms of vortex-mediated friction between the components as

$\mathcal{T}=-\frac{x_{p}\,x_{1}}{\tau}(\Omega_{1}-\Omega_{p})=-x_{1}2\Omega\mathcal{B}(\Omega_{1}-\Omega_{p})$
(28)

where $\Omega\approx\Omega_{p}\approx\Omega_{1}$ is a constant benchmark value
for the observed angular velocity and $\mathcal{B}$ is a dimensionless
coupling constant that sets the intensity of the vortex-mediated mutual
friction; the connection of all these body-averaged parameters with the local
hydrodynamic parameters for nuclear matter is derived in App. C of Montoli et
al. (2020). The timescale

$\tau=x_{p}/(2\mathcal{B}\Omega)$ (29)

is the observed _relaxation_ timescale, namely the timescale with which the
$m=1$ system exponentially relaxes back to the steady state when perturbed
away from it.

The system admits a steady state in which the value of the lag
$\Omega_{1}-\Omega_{p}$ is constant over time: imposing the steady state
condition $\dot{\Omega}_{p}=\dot{\Omega}_{1}=-|\dot{\Omega}_{\infty}|$, we
have

$\begin{split}&\Omega_{1}-\Omega_{p}=\tau|\dot{\Omega}_{\infty}|/x_{p}=|\dot{\Omega}_{\infty}|/(2\Omega\mathcal{B})\\\
&|\mathcal{T}|=x_{1}\,|\dot{\Omega}_{\infty}|\end{split}\quad\qquad\text{(steady
state)}$ (30)

We can turn (27) into a stochastic equation by introducing an additive white
noise term to both the external torque and internal torques:

$-|\dot{\Omega}_{\infty}|\rightarrow-|\dot{\Omega}_{\infty}|+\eta^{\infty}_{t}\qquad\quad\mathcal{T}\rightarrow\mathcal{T}+\eta^{\mathcal{T}}_{t}$
(31)

Now, we can recast (27) in the form (6) by setting

$\mathbf{\Omega}=\begin{pmatrix}\dot{\Omega}_{p}\\\
\dot{\Omega}_{1}\end{pmatrix}\quad B=\begin{pmatrix}-x_{1}/\tau&x_{1}/\tau\\\
x_{p}/\tau&-x_{p}/\tau\end{pmatrix}\quad\mathbf{A}=\begin{pmatrix}-|\dot{\Omega}_{\infty}|/x_{p}\\\
0\end{pmatrix}$ (32)

The matrix $M$ can be found by writing the two independent fluctuations in
(31) in terms of the standard white noise
$\dot{\mathbf{W}}_{t}=(\dot{W}^{\infty}_{t},\dot{W}^{\mathcal{T}}_{t})$ as

$\eta^{\infty}_{t}=\sigma_{\infty}\,\dot{W}^{\infty}_{t}\qquad\quad\eta^{\mathcal{T}}_{t}=\sigma_{\mathcal{T}}\,\dot{W}^{\mathcal{T}}_{t}\,,$
(33)

implying that (only the relative signs of the components of $M$ are important)

$M=\begin{pmatrix}\sigma_{\infty}/x_{p}&-\sigma_{\mathcal{T}}/x_{p}\\\
0&\sigma_{\mathcal{T}}/x_{1}\end{pmatrix}$ (34)

Differently from Meyers et al. (2021a), see App. A.1, we do not include the
additive noise terms as independent fluctuations in the expressions for
$\dot{\Omega}_{p}$ and $\dot{\Omega}_{1}$. Rather, we ascribe them to
different independent physical mechanisms: the external (electromagnetic) and
internal (vortex-mediated) torques. By using (21) the PSD for the observable
component is

$P_{p}(\omega)=\dfrac{(\sigma_{\mathcal{T}}^{2}+\sigma_{\infty}^{2})\tau^{2}\omega^{2}+x_{p}^{2}\sigma_{\infty}^{2}}{x_{p}^{2}(\tau^{2}\omega^{4}+\omega^{2})}\,,$
(35)

that reduces to (24) in the double limit $\sigma_{\mathcal{T}}^{2}\rightarrow
0$, $x_{p}\rightarrow 1$. To compare this result with the one for the $m=2$
model described in the next subsection, we may express the PSD in terms of the
convenient coupling parameter

$b\,=\,2\,\Omega\,\mathcal{B}\,=\,x_{p}/\tau\,,$ (36)

obtaining ($\mu$ and $\xi$ are two parameters that define the PSD)

$\begin{split}&P_{p}(\omega)\propto\dfrac{1}{\omega^{2}}\cdot\dfrac{\mu^{2}+\omega^{2}}{\xi^{2}+\omega^{2}}\\\
&\mu^{2}=\dfrac{b^{2}\,\sigma_{\infty}^{2}}{\sigma_{\infty}^{2}+\sigma_{\mathcal{T}}^{2}}\qquad\qquad\xi^{2}=\dfrac{b^{2}}{x_{p}^{2}}=\dfrac{1}{\tau^{2}}\,.\end{split}$
(37)

The above expressions guarantee that the parameters $\xi$ and $\mu$ are
ordered as

$0<\mu<\xi\,,$ (38)

Moreover, $\mu$ and $\xi$ are the “corner frequencies” that define the two
spectral turnovers of the PSD, see Section 4.

###  2.3 Three-component model $(m=2)$

The straightforward $m=2$ extension of the model of Baym and collaborators
reads,

$\begin{split}&x_{p}\dot{\Omega}_{p}=-\mathcal{T}_{1}-\mathcal{T}_{2}-|\dot{\Omega}_{\infty}|\\\
&x_{1}\dot{\Omega}_{1}=\mathcal{T}_{1}\\\
&x_{2}\dot{\Omega}_{2}=\mathcal{T}_{2}\,,\end{split}$ (39)

where $x_{p}=1-x_{1}-x_{2}$ and

$\mathcal{T}_{i}=-x_{i}b_{i}(\Omega_{i}-\Omega_{p})\qquad i=1,2$ (40)

This is exactly the deterministic model in equation (1) of Montoli et al.
(2020) and its general solution and properties are discussed therein. The
coupling parameters $b_{i}=2\Omega\mathcal{B}_{i}$ set the strength of the
deterministic part of the mutual friction, while the meaning of the
phenomenological dimensionless friction parameter $\mathcal{B}_{i}$ in terms
of the microscopic input is the one presented in App. C of Montoli et al.
(2020). Note that for $m>1$ the relaxation timescales are complicated
functions of $x_{1}$, $x_{2}$, $b_{1}$ and $b_{2}$, see equation (A.19)
therein. Similarly to (30), the steady state lags are, for $i=1,2$,

$\begin{split}&\Omega_{i}-\Omega_{p}=|\dot{\Omega}_{\infty}|/b_{i}=|\dot{\Omega}_{\infty}|/(2\Omega\mathcal{B}_{i})\\\
&|\mathcal{T}_{i}|=x_{i}\,|\dot{\Omega}_{\infty}|\end{split}\quad\qquad\text{(steady
state)}$ (41)

We promote the deterministic model in (39) to a stochastic one by adding
fluctuating terms to the external torque and to the two internal torques.
Following the same procedure outlined in the previous subsection, we define
$\mathbf{\Omega}=(\dot{\Omega}_{p},\dot{\Omega}_{1},\dot{\Omega}_{2})^{\top}$
and

$B=\begin{pmatrix}\dfrac{-x_{1}b_{1}-x_{2}b_{2}}{x_{p}}&\dfrac{x_{1}b_{1}}{x_{p}}&\dfrac{x_{2}b_{2}}{x_{p}}\,\\\
b_{1}&-b_{1}&0\\\
b_{2}&0&-b_{2}\end{pmatrix}\qquad\mathbf{A}=\begin{pmatrix}-\dfrac{|\dot{\Omega}_{\infty}|}{x_{p}}\,\\\
0\\\ 0\end{pmatrix}$ (42)

The matrix $M$ is given by

$M=\begin{pmatrix}\sigma_{\infty}/x_{p}&-\sigma_{\mathcal{T}_{1}}/x_{p}&-\sigma_{\mathcal{T}_{2}}/x_{p}\\\
0&\sigma_{\mathcal{T}_{1}}/x_{1}&0\\\
0&0&\sigma_{\mathcal{T}_{2}}/x_{2}\end{pmatrix}$ (43)

Using (21), the PSD for the observable component is

$\begin{split}&P_{p}(\omega)=\dfrac{N_{0}+N_{2}\omega^{2}+N_{4}\omega^{4}}{D_{2}\omega^{2}+D_{4}\omega^{4}+D_{6}\omega^{6}}\propto\dfrac{n_{0}+n_{2}\omega^{2}+\omega^{4}}{d_{2}\omega^{2}+d_{4}\omega^{4}+\omega^{6}}\\\
&N_{0}=b_{1}^{2}b_{2}^{2}\sigma_{\infty}^{2}\\\
&N_{2}=b_{1}^{2}(\sigma_{\mathcal{T}_{2}}^{2}+\sigma_{\infty}^{2})+b_{2}^{2}(\sigma_{\mathcal{T}_{1}}^{2}+\sigma_{\infty}^{2})\\\
&N_{4}=\sigma_{\mathcal{T}_{1}}^{2}+\sigma_{\mathcal{T}_{2}}^{2}+\sigma_{\infty}^{2}\\\
&D_{2}=b_{1}^{2}b_{2}^{2}\\\
&D_{4}=b_{1}^{2}(1-x_{2})^{2}+b_{2}^{2}(1-x_{1})^{2}+2b_{1}b_{2}x_{1}x_{2}\\\
&D_{6}=x_{p}^{2}\\\ &n_{0}=\dfrac{N_{0}}{N_{4}}\qquad
n_{2}=\dfrac{N_{2}}{N_{4}}\qquad d_{2}=\dfrac{D_{2}}{D_{6}}\qquad
d_{4}=\dfrac{D_{4}}{D_{6}}\,.\end{split}$ (44)

In the double limit $(x_{2}\,,\,\sigma_{\mathcal{T}_{2}})\rightarrow 0$, the
above result reduces to the PSD of the $m=1$ model in (35) with
$b_{1}=x_{p}/\tau$. Again, we can introduce four positive corner frequencies,
$\mu_{\pm}$ and $\xi_{\pm}$, to obtain a form of the PSD that is analogous to
the one in (37):

$\begin{split}&P_{p}(\omega)\propto\dfrac{1}{\omega^{2}}\cdot\dfrac{(\mu_{-}^{2}+\omega^{2})(\mu_{+}^{2}+\omega^{2})}{(\xi_{-}^{2}+\omega^{2})(\xi_{+}^{2}+\omega^{2})}\\\
&\mu_{\pm}^{2}=\dfrac{1}{2}\left[\,n_{2}\pm\sqrt{n_{2}^{2}-4n_{0}}\,\right]\qquad\xi_{\pm}^{2}=\dfrac{1}{2}\left[\,d_{4}\pm\sqrt{d_{4}^{2}-4d_{2}}\,\right]\,.\end{split}$
(45)

Since $0<n_{0}<d_{2}$ and $0<n_{2}<d_{4}$, which can be seen from the formulae
in (44), we have that the four corner frequencies $\mu_{\pm}$ and $\xi_{\pm}$
are naturally ordered as

$0<\mu_{-}\leq\mu_{+}\qquad\quad 0<\xi_{-}\leq\xi_{+}\qquad\quad
0<\mu_{-}<\xi_{+}$ (46)

Their values define where a spectral turnover should be expected, see section
4.

##  3 Setting the strength of fluctuations

In this section, we discuss a parametrization of the phenomenological
quantities $\sigma_{\infty}$ and $\sigma_{\mathcal{T}}$ in terms of the
physical parameters appearing in the deterministic part of the model, that
have a more transparent hydrodynamic interpretation (Montoli et al., 2020).

To set the amplitude of the fluctuations in both the external and internal
torque in terms of the physical parameters we adopt the following scheme. Let
us consider a generic “force” $F=\delta F+F_{0}$, where $\delta F$ is a white
noise process and $F_{0}$ is deterministic, and its associated physical
timescale $T_{F}$, in the sense that the effect of $F_{0}$ is to drive the
system’s evolution on timescales of order $T_{F}$. We can set the strength of
$\delta F$ by requiring that its effect, over a time interval of order
$T_{F}$, is a fraction $0<\alpha<1$ of the one produced by the deterministic
part $F_{0}$,

$\int_{0}^{T_{F}}\delta F(t)\,dt\approx\alpha\int_{0}^{T_{F}}F_{0}(t)\,dt\,.$
(47)

Strictly speaking, the above equation only serves to convey the general idea
and is not formally correct (the left hand side is a random variable, but on
the right we have a deterministic value). However, the left hand side is
typically realised within a few standard deviations from its statistical
expectation, that is zero. Therefore, the size of one of its “typical”
realisations is given by its standard deviation (this would be the rigorous
interpretation of the above equation). In view of this, our prescription to
set the strength of the noise is

$\text{std. dev.}\left[\int_{0}^{T_{F}}\delta
F(t)\,dt\right]=\alpha\int_{0}^{T_{F}}F_{0}(t)\,dt\,.$ (48)

In practice, we are requiring that, over the relevant timescale $T_{F}$, the
impulse caused by fluctuations is (typically) a fraction $\alpha$ of the
impulse imparted by $F_{0}$. In this way, the noise strength can be
parametrized in terms of $\alpha$, which is a relative measure of fluctuations
with respect to the deterministic part of the process.

In the following, to estimate the standard deviation in (48), we will use the
property that the increments

$W^{i}_{t+\Delta t}-W^{i}_{t}\sim\mathcal{N}\left(0,\sqrt{\Delta
t}\right)\qquad\quad\forall\,i,\,\Delta t>0\,.$ (49)

are normally distributed with standard deviation $\sqrt{\Delta t}$.

###  3.1 Noise in the external torque

In order to set the strength of the external torque $\sigma_{\infty}$, it is
natural to consider the process in (3), namely

$\dot{L}\,dt=-|\dot{\Omega}_{\infty}|+\sigma_{\infty}dW^{\infty}_{t}$ (50)

The integration should be performed over the relevant timescale, namely the
pulsar’s life timescale666 For a “standard” braking index of 3, the pulsar’s
characteristic dipole age would be $T/2$.  $T=\Omega/\dot{\Omega}$, where
$\Omega\approx\Omega_{p}$ and $\dot{\Omega}\approx|\dot{\Omega}_{\infty}|$ are
the observed rotational parameters of the pulsar. Thanks to the fundamental
property in (49), we have

$\int_{0}^{T}\\!\\!\sigma_{\infty}dW^{\infty}_{t}=\sigma_{\infty}(W^{\infty}_{t+T}-W^{\infty}_{t})\sim\mathcal{N}\left(\,0\,,\,\sigma_{\infty}\sqrt{T}\,\right)\,,$
(51)

Now, we require that the standard deviation of this process is a fraction
$0<\alpha_{\infty}<1$ of the observed angular velocity $\Omega$,

$\sigma_{\infty}\sqrt{T}\,=\,\alpha_{\infty}\,\Omega\qquad\qquad\sigma_{\infty}^{2}\,=\,\alpha_{\infty}^{2}\,\Omega\,\dot{\Omega}\,.$
(52)

This choice encodes the fact that random fluctuations in the spin down torque
can result in a total change of the present angular velocity $\Omega$, of the
order of a fraction $\alpha_{\infty}$ of it, when their cumulative effect over
the pulsar’s life is taken into account. In other words, we are imposing that
the typical effect of fluctuations in the external torque piles up to a
fraction of the effect of the deterministic spin-down

$\left|\int_{0}^{T}\mathcal{T}_{\infty}\,dt\,\right|=|\dot{\Omega}_{\infty}|T\approx\Omega\,.$
(53)

###  3.2 Noise in the internal torque

We apply the same reasoning to set the strength of the internal torques. For
$m=1$, the relevant timescale is in terms of the coupling parameter $b$ in
(36). We demand that the “typical” effect of the fluctuation in a time
interval of order $\sim b^{-1}$

$\int_{0}^{b^{-1}}\\!\\!\\!\eta_{\mathcal{T}}\,dt=\int_{0}^{b^{-1}}\\!\\!\\!\sigma_{\mathcal{T}}\,dW^{\mathcal{T}}\sim\mathcal{N}\left(\,0\,,\,\sigma_{\mathcal{T}}/\sqrt{b}\,\right)$
(54)

amounts to a fraction $0<\alpha<1$ of the effect of the steady state internal
torque (30), that is given by

$\left|\int_{0}^{b^{-1}}\\!\\!\\!\mathcal{T}\,dt\,\right|=\frac{|\mathcal{T}|}{b}\,=\,\frac{x_{1}\,|\dot{\Omega}_{\infty}|}{b}\,\approx\,\frac{x_{1}\,\dot{\Omega}}{b}\,.$
(55)

Therefore, we require that,

$\int_{0}^{b^{-1}}\\!\\!\\!\eta_{\mathcal{T}}\,dt\sim\mathcal{N}\left(0\,,\,\alpha\,x_{1}\dot{\Omega}/b\,\right)\,,$
(56)

namely

$\sigma_{\mathcal{T}}/\sqrt{b}=\alpha\,x_{1}\dot{\Omega}/b\qquad\qquad\sigma_{\mathcal{T}}^{2}\,=\,\dfrac{\alpha^{2}\,x_{1}^{2}\,\dot{\Omega}^{2}}{2\,\mathcal{B}\,\Omega}\,.$
(57)

For $m=2$, the parametrization of the fluctuations in the external torque is
still given by (52), while the analogue of (57) is

$\sigma_{\mathcal{T}i}^{2}=\frac{\alpha_{i}^{2}\,x_{i}^{2}\,\dot{\Omega}^{2}}{2\,\mathcal{B}_{i}\,\Omega}\qquad\,\,0<\alpha_{i}<1\qquad\,\,i=1,2\,.$
(58)

##  4 Spectral flattening and corner frequencies

We can now rewrite the PSD in terms of our “natural” choice of the noise
variables $\sigma_{\infty}$ and $\sigma_{\mathcal{T}}$ derived the previous
section. In this way, we can link the presence of a possible spectral turnover
and the corresponding corner frequencies to the physical properties (i.e.,
$x_{p}$, $\tau$, $T$…) of a given pulsar.

###  4.1 Spectral turnover for $m=1$

The PSD for the angular velocity in (37) is defined in terms of two
frequencies $\mu$ and $\xi$. If $\mu\ll\xi$, we can identify 3 power-law
regimes with spectral index $-2$, $0$ and $-2$:

$P_{p}(\omega)\propto\begin{cases}\,\,\mu^{2}/(\xi^{2}\,\omega^{2})\qquad&\omega\ll\mu\\\
\,\,1/\xi^{2}\qquad&\mu\ll\omega\ll\xi\\\
\,\,1/\omega^{2}\qquad&\xi\ll\omega\end{cases}$ (59)

Since the PSD undergoes a spectral turnover around $\omega\approx\xi,\mu$,
those two frequencies play the role of corner frequencies, see e.g. Coles et
al. (2011). By looking at (37), one may also expect a possible regime
$P_{p}\propto 1/\omega^{4}$. However, given our setting, this possibility is
forbidden by the parametrization of the PSD in terms of the physical
parameters of the spin-down model: to obtain a region where $P_{p}\propto
1/\omega^{4}$, we should have $\xi\ll\omega\ll\mu$, that is inconsistent with
the physical ordering $\mu<\xi$. This is provides a possible way to test the
model versus observations: we know that the spectral index $-4$ for the PSD of
the angular velocity can never be realised for any choice of the model’s
parameters in (37). On the other hand, observational evidence for a spectral
index of $-4$ may be taken as an indication that gravitational radiation
(which is not accounted for in the present model, but may be included as done
in Meyers et al. 2021a, see App. A) is contributing to spinning down the
pulsar.

In general, a certain ordering of the corner frequencies gives rise to a
sequence of spectral indexes in the PSD,

$\begin{split}&\text{Possible
case:}\quad&\mu\,<\,\xi\qquad&(\,-2,\,\,\,0,-2\,)\\\ &\text{Inconsistent
case:}\quad&\xi\,<\,\mu\qquad&(\,-2,-4,-2)\end{split}$ (60)

The underlying physical model tells us which spectral index sequence is
consistent or inconsistent with it. The spectral sequences for the phase
residuals are obtained from the ones above by adding $-2$ to each entry, e.g.
$(-4,-2,-4)$ for the first case in (60). We can now use the noise
parametrization in (52) and (57) to write the corner frequencies as

$\mu^{2}\,=\,\frac{1}{1+(x_{1}^{2}/x_{p})(\alpha_{\mathcal{T}}/\alpha_{\infty})^{2}(\tau/T)}\cdot\dfrac{x_{p}^{2}}{\tau^{2}}\,<\,\xi^{2}\,=\,\dfrac{1}{\tau^{2}}\,,$
(61)

which tells us that the region $\omega\in[x_{p}\,\xi\,,\,\xi]$ is always flat,
whatever the choice of the parameters (this is, however, only the minimum
possible extension of the flat region). In any case, to properly test the
presence of this flattened frequency interval, the observed PSD should contain
relevant information on the low-frequency domain, meaning that long
observations (at least as long as a few times $\tau$) are needed. Moreover,
according to our model pulsars with a bigger $\tau/T$ ratio are expected to
display a more extended flat region. Below, we consider in more detail two
interesting extreme cases.

_Pure external noise -_ For $\alpha_{\mathcal{T}}=0$ we have $\xi=1/\tau$ and
$\mu=x_{p}/\tau$: the flat region $\mu<\omega<\xi$ can not be less extended
than this. The PSD in (35) reduces to

$P_{p}(\omega)=\frac{\alpha_{\infty}^{2}\,\Omega\,\dot{\Omega}}{x_{p}^{2}\,\omega^{2}}\cdot\frac{\omega^{2}+x_{p}^{2}\tau^{-2}}{\omega^{2}+\tau^{-2}}\,\approx\frac{\alpha_{\infty}^{2}\Omega\dot{\Omega}}{\omega^{2}}\,,$
(62)

where the approximation is valid only in the limit in which the loosely
coupled superfluid component is restricted to the inner crust (for
$x_{p}\approx 1$, we have $\mu\approx\xi$ and the flattened region
disappears). As expected, in this limit we recover exactly the $m=0$ result in
(26).

_Pure internal noise -_ For $\alpha_{\infty}=0$ we have $\xi=1/\tau$ and
$\mu=0$, meaning that the flat region extends indefinitely at low frequencies.
The PSD in (35) reduces to

$P_{p}(\omega)=\frac{\alpha_{\mathcal{T}}^{2}\,x_{1}^{2}\,\dot{\Omega}^{2}}{2\,x_{p}^{2}\mathcal{B}\,\Omega}\cdot\frac{1}{\omega^{2}+\tau^{-2}}\,.$
(63)

This is the usual Lorentzian form of the PSD of a mean-reverting Ornstein-
Uhlenbeck process, see Appendix B.

###  4.2 Spectral turnover for $m=2$

For the three-component model $m=2$ the analysis of the PSD requires more
care. Given the PSD in (45), we have 5 possibilities, which depend on the
relative ordering of the four corner frequencies: following the scheme already
used in (60), we have that (only the possible cases that are consistent with
(45) are listed)

$\begin{split}\text{case
1:}\quad&\mu_{-}\,<\,\mu_{+}\,<\,\xi_{-}\,<\,\xi_{+}&(\,-2,\,\,\,0,+2,\,\,\,0,-2\,)\\\
\text{case
2:}\quad&\mu_{-}\,<\,\xi_{-}\,<\,\mu_{+}\,<\,\xi_{+}&(\,-2,\,\,\,0,-2,\,\,\,0,-2\,)\end{split}$
(64)

For example, if the physical parameters of the $m=2$ model are such that case
1 is realised, then the expected sequence of spectral indexes is
$(-2,0,+2,0,-2)$, namely

$P_{p}(\omega)\propto\begin{cases}\,\,\mu_{-}^{2}\mu_{+}^{2}/(\xi_{-}^{2}\xi_{+}^{2}\,\omega^{2})\qquad&\omega\ll\mu_{-}\\\
\,\,\mu_{+}^{2}/(\xi_{-}^{2}\xi_{+}^{2})\qquad&\mu_{-}\ll\omega\ll\mu_{+}\\\
\,\,\omega^{2}/(\xi_{-}^{2}\xi_{+}^{2}\,)\qquad&\mu_{+}\ll\omega\ll\xi_{-}\\\
\,\,1/\xi_{+}^{2}\qquad&\xi_{-}\ll\omega\ll\xi_{+}\\\
\,\,1/\omega^{2}\qquad&\xi_{+}\ll\omega\end{cases}$ (65)

Again, the spectral sequences for the phase residuals are obtained from the
ones above by adding $-2$ to each entry. Below, we consider in more detail two
interesting extreme cases. With no loss of generality, we follow Pizzochero et
al. (2020) and set $0<b_{1}<b_{2}$: we will refer to the component 1 as
“loose” superfluid component and to component 2 as “tight” superfluid
component (the loose component is the one that is less coupled to the normal
one).

_Pure external noise -_ For
$\alpha_{\mathcal{T}_{1}}=\alpha_{\mathcal{T}_{2}}=0$ we have that the two
corner frequencies at the denominator $\xi_{\pm}$ are unchanged (they do not
depend on the noise parameters), while

$\mu_{-}=b_{1}\qquad\mu_{+}=b_{2}\,.$ (66)

Moreover, it is possible to show that only case 2 in (64) can be realised,
since

$b_{1}<\xi_{-}<b_{2}<\xi_{+}\,.$ (67)

The PSD in (44) boils down to

$P_{p}(\omega)=\frac{\alpha_{\infty}^{2}\Omega\dot{\Omega}}{x_{p}^{2}\,\omega^{2}}\cdot\frac{(\omega^{2}+b_{1}^{2})(\omega^{2}+b_{2}^{2})}{(\omega^{2}+\xi_{-}^{2})(\omega^{2}+\xi_{+}^{2})}\,.$
(68)

_Pure internal noise in the “loose” component -_ For $b_{1}<b_{2}$, we set
$\alpha_{\infty}=\alpha_{\mathcal{T}_{2}}=0$, so that only the less coupled
part of the superfluid is responsible for the fluctuations of the internal
torque. This process is a Ornstein-Uhlenbeck process (see App. B) and is
asymptotically stationary. Again, only case 2 in (64) can be realised, since

$\mu_{-}=0\qquad\mu_{+}=b_{2}\,.$ (69)

and

$0<\xi_{-}<b_{2}<\xi_{+}\,.$ (70)

The PSD in (44) reads, cf. with (63),

$P_{p}(\omega)=\frac{\alpha_{\mathcal{T}_{1}}^{2}\Omega\dot{\Omega}}{x_{p}^{2}\,\omega^{2}}\cdot\frac{x_{1}^{2}}{b_{1}T}\cdot\frac{\omega^{2}(\omega^{2}+b_{2}^{2})}{(\omega^{2}+\xi_{-}^{2})(\omega^{2}+\xi_{+}^{2})}\,.$
(71)

Recalling that $b_{i}=2\Omega\mathcal{B}_{i}$, the PSD can also be written as

$P_{p}(\omega)=\frac{\alpha_{\mathcal{T}_{1}}^{2}\dot{\Omega}^{2}}{x_{p}^{2}}\cdot\frac{x_{1}^{2}}{2\Omega\mathcal{B}_{1}}\cdot\frac{(\omega^{2}+b_{2}^{2})}{(\omega^{2}+\xi_{-}^{2})(\omega^{2}+\xi_{+}^{2})}\,.$
(72)

_Pure internal noise in the “tight” component -_ For $b_{1}<b_{2}$, we set
$\alpha_{\infty}=\alpha_{\mathcal{T}_{1}}=0$, so that only the more tightly
coupled part of the superfluid is responsible for the fluctuations of the
internal torque. Again, this process is a Ornstein-Uhlenbeck process (see App.
B) and is asymptotically stationary. Now, only case 1 in (64) can be realised,
since

$\mu_{-}=0\qquad\mu_{+}=b_{1}\,.$ (73)

and

$0<b_{1}<\xi_{-}<\xi_{+}\,.$ (74)

The PSD in (44) reads, cf. with (63),

$P_{p}(\omega)=\frac{\alpha_{\mathcal{T}_{2}}^{2}\Omega\dot{\Omega}}{x_{p}^{2}\,\omega^{2}}\cdot\frac{x_{2}^{2}}{b_{2}T}\cdot\frac{\omega^{2}(\omega^{2}+b_{1}^{2})}{(\omega^{2}+\xi_{-}^{2})(\omega^{2}+\xi_{+}^{2})}\,.$
(75)

and we can have a region of the spectrum $b_{1}<\omega<\xi_{-}$ that is blue
(spectral index $+2$ for the angular velocity).

##  5 Timing noise strength

Several metrics have been proposed to quantify the strength of observed timing
noise in a certain pulsar. To extract the timing noise strength from our
simulated signal $\Omega_{p}$, we have to mimic the observational procedure.
First, we introduce the total phase difference of the observable signal
between two arbitrary times $0$ and $t$ in the usual way,

${\phi}_{t}=\int_{0}^{t}{\Omega_{p}}(z)\,dz\,.$ (76)

In timing observations, the phase residuals (i.e., the observational
counterpart of our theoretical process $\delta\phi_{t}$) are measured with
respect to a given timing model with phase $\mathbf{\phi}^{\text{mod}}_{t}$,
typically modelled as a power series in $t$ (Lorimer & Kramer, 2005), namely

$\delta{\phi}_{t}={\phi}_{t}-{\phi}^{\text{mod}}_{t}\,\qquad\dot{\phi}^{\text{mod}}_{t}=\Omega^{\text{mod}}+\dot{\Omega}^{\text{mod}}\,t+...$
(77)

In our case, to extract the theoretical residuals we can set,

$\phi^{\text{mod}}_{t}=\Omega_{p}(0)\,t-|\dot{\Omega}_{\infty}|\,t^{2}/2\,,$
(78)

or, in terms of angular velocity,

$\delta\Omega_{p}=\Omega_{p}(t)-\dot{\phi}^{\text{mod}}_{t}\qquad\qquad\langle\Omega_{p}(t)\,\rangle\,=\,\dot{\phi}^{\text{mod}}_{t}\,.$
(79)

Therefore, the biggest difference between observations and our theoretical
setting is that in real observations the timing model
${\phi}^{\text{mod}}_{t}$ is derived from the prior knowledge of the
rotational ($\Omega$ , $\dot{\Omega}$) and astrometric parameters. We may
follow the same procedure and consider a more general $\phi^{\text{mod}}_{t}$
that does not necessarily coincide with the deterministic solution
($\Omega_{p}(0)$ and $|\dot{\Omega}_{\infty}|$ are just the input of a
mathematical model and can not be extracted directly from observations),
namely $\phi^{\text{mod}}_{t}=\phi_{0}+\Omega\,t-|\dot{\Omega}|\,t^{2}/2$,
with $\Omega\approx\Omega_{p}(0)$ and
$\dot{\Omega}\approx|\dot{\Omega}_{\infty}|$ inferred directly from the
simulated process over a finite time window. However, this complicates
considerably the following calculations, so we stick to the prescriptions in
(78) and (79). Therefore, our analysis refers to the ideal case in which the
long-term spin-down rate $\dot{\Omega}$ extracted from a fit to the data well
approximates the intrinsic (not directly observable) parameter
$|\dot{\Omega}_{\infty}|$ that sets the average spin-down rate of the pulsar
over many decades.

###  5.1 Reconstructing the times of arrival

Numerical simulations of pulsar rotation are based on evolving some
(simplified) dynamical equations for the angular velocity, which is then
obtained directly by numerical integration of the assumed theoretical model.
However, in pulsar timing analysis, the quantity that is measured directly are
the TOAs of a train of pulses, while the phase $\phi(t)$ and the angular
velocity are derived quantities. Therefore, if we want to study the residuals
of the TOAs with respect to a certain timing model $\phi^{\text{mod}}(t)$, we
first have to reconstruct them.

Let us indicate the measured pulsar’s TOAs by $\\{t_{j}\\}$, while the TOAs of
the reference polynomial model are $\\{t^{\text{mod}}_{j}\\}$, meaning that
each $t_{j}$ and $t^{\text{mod}}_{j}$ are the times such that the pulsar’s
phase $\phi\left(t_{j}\right)$ coincides with the one of the model777 From the
observational point of view, a fiducial point in the observed pulse profile is
chosen a priori to keep track of the pulsar’s phase. The pulsar must have spun
around by an integer number of rotations $\Delta N$ between two measured TOAs
$t_{j}$ and $t_{j+1}$, provided there is no change in the assumed timing
model. Therefore, (80) is valid because according to the observational
procedure used to extract the TOAs, $\phi(t_{j+1})-\phi(t_{j})=2\pi\Delta
N=\phi^{\text{mod}}(t^{\text{mod}}_{j+1})-\phi(t^{\text{mod}}_{j})$ and we are
free to set the phase relative to the first TOA equal to zero.

$\phi\left(t_{j}\right)=\phi^{\text{mod}}\left(t^{\text{mod}}_{j}\right)\,.$
(80)

In terms of the timing residuals

$\delta a_{j}=t_{j}-t^{\text{mod}}_{j}\,,$ (81)

that measure how much the observed TOAs advance or lag behind the ones
expected for the timing model, the above equation reads

$\phi(t^{\text{mod}}_{j}+\delta
a_{j})=\phi^{\text{mod}}(t^{\text{mod}}_{j})\,.$ (82)

We can now generalise this concept to a continuous observation in time by
demanding that $\delta a_{t}$ is the solution of the implicit equation

$\phi(t+\delta a_{t})=\phi^{\text{mod}}(t)\,.$ (83)

Implicit differentiation of the above equation gives

$\dfrac{d}{dt}\,\delta
a_{t}=-\dfrac{\Omega_{p}(t)-\dot{{\phi}}^{\text{mod}}_{t}}{\Omega_{p}(t)}-\dfrac{\dot{{\phi}}^{\text{mod}}_{t}\dot{\Omega}_{p}(t)}{{\Omega_{p}(t)}^{2}}\,\delta
a_{t}+O\left(\dfrac{\delta a_{t}^{2}}{T^{2}}\right)\,,$ (84)

where $T$ is of the order of the pulsar’s age. On timescales much smaller than
the pulsar age, we can approximate $\Omega_{p}$ with a suitable constant value
$\Omega$, so that the expression

$\delta a_{t}\,\approx\,-\,\dfrac{\delta\phi(t)}{\Omega}$ (85)

is (for our purposes) a very good approximation to the exact solution of (84)
for time intervals much shorter than $T$: while (84) leads to a more accurate
reconstruction of $\delta a_{t}$, for the sake of studying the timing noise
strength $\sigma^{2}$ (to be introduced in the next subsection), the simple
prescription (85) turns out to be good enough.

###  5.2 Timing noise strength

The variance $\sigma^{2}$ of the timing residuals from a least squares
polynomial fit over an observation interval of length $T_{o}$ provides a
measure of the timing noise strength in real timing data (Cordes & Helfand,
1980; Hobbs et al., 2010; Shannon & Cordes, 2010):

$\sigma^{2}=\frac{1}{N_{o}}\sum_{j=1}^{N_{o}}\,\left|\delta
a_{j}\right|^{2}\,,$ (86)

where $\delta a_{j}$ are the measured timing residuals in (81) and $N_{o}$ is
the number of measurements. For a continuous-time model we can express the
above discrete expression in terms of the variance of the continuous process
$\delta a_{t}$,

$\sigma^{2}=\frac{1}{T_{o}}\int_{t_{0}}^{t_{0}+T_{o}}\\!\\!dz\,\left|\delta
a_{z}\right|^{2}\,\approx\frac{1}{T_{o}\,\Omega^{2}}\int_{t_{0}}^{t_{0}+T_{o}}\\!\\!dz\,\left|\delta\phi_{z}\right|^{2}\,,$
(87)

where the timing residuals $\delta a_{t}$ have been translated into phase
residuals (that are more easily obtained in numerical simulations) via (85).
In this way, $\sigma$ is the root mean square of the observed phase residuals
${\phi}_{j}$ and its physical dimension is that of time (Cordes & Helfand,
1980). To estimate $\sigma$ from our model, we follow the same procedure that
we have used to obtain the PSD for the angular velocity. First, we introduce
the PSD for the phase residuals and for the TOA residuals as

$P_{\phi}(\omega)=\omega^{-2}\,P_{p}(\omega)\qquad\quad
P_{a}(\omega)\approx(\omega\Omega)^{-2}\,P_{p}(\omega)\,,$ (88)

where the expression for the PSD of the TOAs residuals $P_{a}(\omega)$ is an
approximation based on (85). We may estimate $\sigma^{2}$ by considering the
total power contained in a band of frequencies $\omega\in[\nu_{o},\nu_{s}]$,
namely

$\sigma^{2}\approx
2\int_{\nu_{o}}^{\nu_{s}}\dfrac{d\omega}{2\pi}\,P_{a}(\omega)\approx\frac{2}{\Omega^{2}}\int_{\nu_{o}}^{\nu_{s}}\dfrac{d\omega}{2\pi}\,P_{\phi}(\omega)\,.$
(89)

The above integral is band-limited because (in observations as well as in
simulations) we can not access frequencies higher than the sampling one
$\nu_{s}\sim N_{o}/T_{o}$. Moreover, the lower frequency $\nu_{0}\sim 1/T_{0}$
acts as a cut-off parameter that regulates the divergence of
$P_{\phi}(\omega)$ for $\omega\rightarrow 0$. The prefactor of $2$ is needed
since we are integrating over positive frequencies only, and the overall
scaling with $\Omega^{2}$ comes from equation (85).

The expression for $\sigma^{2}$ in (89) is a guess and it is not equivalent to
(87); however, the scaling of $\sigma^{2}$ with the physical and rotational
parameters turn out to be the same, whether we use (89) or (87): the two
estimates of $\sigma^{2}$, can differ by a dimensionless constant of order 1
(not surprisingly, given the arbitrariness of the cut-off $\nu_{0}$). We
analytically check this claim in App. B.

###  5.3 Timing noise strength for $m=1$ and $m=2$

In order to see if our model can reproduce the observed scaling of
$\sigma^{2}$ reported in the literature, e.g. (Hobbs et al., 2010; Shannon &
Cordes, 2010; Parthasarathy et al., 2019) we explicitly calculate the band-
limited integral in (89) for the $m=1$ model.

The timing noise strength is often studied across a sample of pulsars. For
example, Cordes & Helfand (1980) found that, out of a sample of 11 pulsars,
the timing noise strength is correlated with the period derivative and weakly
with the period. More precisely, analysis of the timing noise strength is
based on the study of correlations in the data couples
$(\,\Omega^{a}_{i}\dot{\Omega}^{b}_{i}\,,\,\sigma_{i}\,)$, where $\Omega_{i}$,
$\dot{\Omega}_{i}$ and $\sigma_{i}$ are observed properties of a given $i$-th
pulsar. Following Hobbs et al. (2010), multiple combinations of $(a,b)$ can be
tested to maximise some measure of correlation (like the Pearson coefficient),
which leads them to find that
$\sigma\propto-1.37\log_{10}[\Omega^{0.29}\dot{\Omega}^{-0.55}]$, i.e.
$a\approx-0.39$, $b\approx 0.75$. They also conclude that timing noise is
inversely correlated with characteristic age: in fact, the phenomenological
scaling of $\sigma$ can also be expressed as

$\sigma\propto\Omega^{a}\dot{\Omega}^{b}=\Omega^{a+b}T^{-b}=T^{a}\dot{\Omega}^{a+b}\,,$
(90)

implying that we can equivalently interpret the result of Hobbs and
collaborators as $\sigma\propto P^{-0.36}T^{-0.75}$, where
$P\propto\Omega^{-1}$ is the pulsar’s period.

More recently, Shannon & Cordes (2010) reported an observational scaling
relation of $\sigma\propto\Omega^{-0.9\pm 0.2}\dot{\Omega}^{1\pm 0.05}$ for a
sample of canonical pulsars (CPs) and $a\approx-1.5$, $b\approx 1.2$ for the
whole pulsar population. This relationship has also been investigated by
Parthasarathy et al. (2019); they found that
$\sigma\propto\Omega^{1}\dot{\Omega}^{-0.9\pm 0.2}$ for the sample of young
pulsars, whereas after adding a few millisecond pulsars (MSPs) they found that
$\sigma\propto\Omega^{1}\dot{\Omega}^{-0.6\pm 0.1}$ on the basis of an
analysis based on the Spearman’s rank correlation coefficient. However, this
kind of statistical analysis can only probe the ratio $b/a$ of
$\Omega^{a}\dot{\Omega}^{b}$ (see Jankowski et al. (2017)). In fact, any
rescaling $\Omega^{ca}\dot{\Omega}^{cb}$ with $c>0$ preserves the ranking of
the data: the only thing that the Spearman’s coefficient can probe is the
degree of monotonic behaviour (i.e., the ranking) of the observed data couples
$(\Omega^{1}_{i}\dot{\Omega}^{b/a}_{i},\sigma_{i})$, where $i$ indicates a
particular pulsar. Hence, the relationship found by Parthasarathy et al.
(2019), $b/a\approx-0.6$, is only qualitatively consistent with the previous
result of Shannon & Cordes (2010) for a sample of both canonical pulsars (CP)
and millisecond pulsars (MSP), $b/a\approx-0.8$, see Table 1. Finally, Lower
et al. (2020) have developed a Bayesian framework to estimate the scaling
parameters $a$ and $b$. Their result (see Tab. 6 in their original work), is
also listed in Tab. 1 and is particularly close to the earlier result of
Shannon & Cordes (2010) for non-recycled pulsars (CP).

We now provide a theoretical estimate of the timing noise strength $\sigma$
and its scaling parameters $a$ and $b$ by assuming that timing noise follows
our stochastic model. In principle, it is possible to obtain the exact
analytical expression for (89) by plugging the PSD (35) into (88), but the
final expression is quite convoluted. Therefore, we only consider the two
limiting cases of pure internal and pure external noise. To find the dominant
dependence of $\sigma$ on the rotational parameters, we assume that there is a
hierarchy of timescales, namely

$T_{o}/N_{o}\ll\tau\ll T_{o}\ll
T\quad\quad\left(\,\dot{\Omega}\ll\Omega\,\nu_{o}\ll\mathcal{B}\,\Omega^{2}\ll\Omega\,\nu_{s}\,\right)$
(91)

The explicit calculations for $m=1$ are done by considering the results in
(35) and the formula for $\sigma^{2}$ in (89).

_Pure external noise ( $m=1$) -_ Plugging the exact PSD in (62) into our
definition (89), we obtain the leading terms

$\sigma^{2}\approx\dfrac{\alpha_{\infty}^{2}\dot{\Omega}}{3\pi\nu_{o}^{3}\Omega}\left(1+\dfrac{3(1-x_{p}^{2})\nu_{o}^{2}}{4\mathcal{B}^{2}\Omega^{2}}\right)$
(92)

To obtain this result we have expanded the exact result to implement the
hierarchy of timescales in (91). The leading term tell us that
$\sigma^{2}\propto\Omega^{-1}\dot{\Omega}$.

_Pure internal noise ( $m=1$) - _ We now plug the exact PSD in (63) into (89)
and expand according to the hierarchy of timescales in (91):

$\sigma^{2}\approx\dfrac{\alpha_{\mathcal{T}}^{2}\,x_{1}^{2}\,\dot{\Omega}^{2}}{8\pi\,\nu_{o}\,\mathcal{B}^{3}\,\Omega^{5}}\left(1-\dfrac{\pi\,x_{p}\,\nu_{o}}{4\,\mathcal{B}\,\Omega}\right)$
(93)

Hence, the expected noise strength is expected to scale with the rotational
parameters as $\sigma^{2}\propto\Omega^{-5}\dot{\Omega}^{2}$.

For $m=2$, the calculation proceeds along the same steps as previously, albeit
it is more laborious. Again, we take component 1 to be the loose one, so that
the analogous of (91) is

$\dot{\Omega}\ll\Omega\,\nu_{o}\ll\mathcal{B}_{1}\,\Omega^{2}\ll\mathcal{B}_{2}\,\Omega^{2}\ll\Omega\,\nu_{s}\,.$
(94)

We perform the integral in (89) by using the relations in (44) and then expand
the result thanks to the hierarchy of timescales in the above equation; the
final result reads:

$\begin{split}&\sigma^{2}\approx\dfrac{\alpha_{\infty}^{2}\dot{\Omega}}{3\pi\nu_{o}^{3}\Omega}\left(1+\dfrac{3(2-x_{1})x_{1}\,\nu_{o}^{2}}{4\,\mathcal{B}^{2}_{1}\,\Omega^{2}}\right)\quad&(\alpha_{\infty}\neq
0)\\\
&\sigma^{2}\approx\dfrac{\alpha_{\mathcal{T}1}^{2}\,x_{1}^{2}\,\dot{\Omega}^{2}}{8\pi\,\nu_{o}\,\mathcal{B}_{1}^{3}\,\Omega^{5}}\left(1-\dfrac{\pi\,(1-x_{1})\,\nu_{o}}{4\,\mathcal{B}_{1}\,\Omega}\right)\quad&(\alpha_{\mathcal{T}1}\neq
0)\\\
&\sigma^{2}\approx\dfrac{\alpha_{\mathcal{T}2}^{2}\,x_{2}^{2}\,\dot{\Omega}^{2}}{8\pi\,\nu_{o}\,\mathcal{B}_{2}^{3}\,\Omega^{5}}\left(1+\dfrac{\pi\,(2-x_{1})x_{1}\,\nu_{o}}{4(1-x_{1})\,\mathcal{B}_{1}\,\Omega}\right)\quad&(\alpha_{\mathcal{T}2}\neq
0)\end{split}$ (95)

for the pure external noise ($\alpha_{\mathcal{T}1},\alpha_{\mathcal{T}2}=0$),
noise in the loose component ($\alpha_{\infty},\alpha_{\mathcal{T}2}=0$) and
noise in the tight component ($\alpha_{\mathcal{T}1},\alpha_{\infty}=0$),
respectively. The leading terms are identical to the corresponding ones for
the $m=1$ cases. The sub-leading terms have different dimensionless
parameters, but they still obey the same scaling of the $m=1$ model. In the
end, recalling that $\nu_{0}\sim 1/T_{0}$, we can conclude that, for both
$m=1,2$,

$\begin{split}&\sigma^{2}\approx
a_{e0}\,\dot{\Omega}\Omega^{-1}T_{0}^{3}+a_{e1}\,\dot{\Omega}\Omega^{-3}T_{o}\quad&(\text{pure
external noise})\\\ &\sigma^{2}\approx
a_{i0}\,\dot{\Omega}^{2}\Omega^{-5}T_{o}+a_{i1}\,\dot{\Omega}^{2}\Omega^{-6}\quad&(\text{pure
internal noise})\end{split}$ (96)

where $|a_{e1}|\ll a_{e0}$ and $|a_{i1}|\ll a_{i0}$ are dimensionless
constants that depend on the physical parameters of the pulsar. The results
are summarised in Table 1. We see that the stochastic model for pure external
noise captures the value of the ratio $b/a$, but their individual values are
rather different from the recent Bayesian estimates by Lower et al. (2020).
Therefore, we are left with three possibilities. The first one is that a model
with pure external torque is not good enough to reproduce the scaling and that
some amount of internal fluctuations are necessary (adding internal
fluctuations raises the value of $b$ and decreases $a$, bringing their values
closer to observed ones). The second possibility is that a model with injected
white noise should be abandoned in favour of coloured noise (noise due to
internal turbulence is red, see e.g. Melatos & Link 2014). Finally, it is also
possible that our parameters $\alpha_{\infty,\mathcal{T}}$ depend on the
pulsar’s age, so that they carry extra correlation with $\Omega$ and
$\dot{\Omega}$ that is unaccounted here.

$\sigma\propto\Omega^{a}\dot{\Omega}^{b}$ | $a$ | $b$ | $b/a$
---|---|---|---
“external” noise, eq. (96) | -0.5 | 0.5 | -1
“internal” noise, eq. (96) | -2.5 | 1 | -0.4
Lower et al. (2021) | -0.84 | 0.97 | -1.15
Parthasarathy et al. (2019) | 1 (assumed) | -0.6 | -0.6
Shannon & Cordes (2010) (MSP+CP) | -1.5 | 1.2 | -0.8
Shannon & Cordes (2010) (CP) | -0.9 | 1 | -1.1
Hobbs et al. (2010) | -0.39 | 0.75 | -1.9
Table 1:  Scaling of the timing noise strength $\sigma$ for our model
(considering the leading terms in (96), that are valid for both the $m=1$ and
$m=2$ models) and extracted from pulsar timing observations. The result of
Parthasarathy et al. (2019) is closer to the case of pure internal noise,
while the analysis of Shannon & Cordes (2010) returns a ratio $b/a$ that is
consistent with the one of pure external noise for canonical pulsars (CP). The
values of Hobbs et al. (2010) are obtained from their equation (6) after
multiplying both $a$ and $b$ by their prefactor $-1.37$. For notation
convenience, $\dot{\Omega}>0$ is the absolute value of the spin down rate
throughout the paper.

Keeping in mind the three possibilities above, we stick to our minimal setting
and tentatively explore how our results for the timing noise strength $\sigma$
vary across the known pulsar population. This is done in Fig. 1, where each
pulsar in the $P-\dot{P}$ diagram is plotted according to a colour scheme
representing the theoretical expectation for its timing noise strength: for
each object we use the exact formula (89) for $m=1$ and $\mathcal{B}=2\times
10^{-}9$, $x_{1}=0.1$, $\nu_{o,s}=2\pi/T_{o,s}$, where the observation time is
$T_{o}=20\,$yr and the sampling time is $T_{s}=1$ day. In the case of pure
external noise (right panel) the expected scaling is $\sigma\propto
1/\sqrt{T}$, in accordance with the leading term of the first equation in
(96).

We see that fluctuations in the internal torque tend to define the scaling in
the upper part of the diagram (the one containing younger pulsars). On the
other hand, even a small amount of noise in the external torque dominates in
the region of older pulsars, as can be seen in the central panel.

The use of the exact formula (89) results in corrections that do not modify
the general approximate scaling defined by the leading terms in (96), at least
in the region $P$ and $\dot{P}$ of known pulsars. In fact, equation (96) tells
us that, in the case of pure external noise, the leading term is
$\sigma\propto T^{-0.5}$, which happens to be the case in the right panel of
Fig. 1, where the lines of constant $\sigma$ coincide with the lines of
constant age. On the other hand, for pure internal noise, the approximate
scaling in (96) can be written as $\sigma\propto P^{1/2}\dot{P}$, so that the
lines of constant $\sigma$ should have slope $-0.5$ (to be compared with the
slope $-1$ of the lines of constant magnetic field). Again, this is what can
be seen in the left panel of Fig. 1: the full result (89) is in good agreement
with the approximate scaling defined by the leading term in (89).

Clearly, these results for the pulsar population are biased by the fact that
we have set constant values of $\mathcal{B}$, $x_{1}$, $\alpha_{\infty}$ and
$\alpha_{\mathcal{T}}$ across the whole population. Actually, these physical
parameters are likely to have an intrinsic scaling with the pulsar age (e.g.,
mature pulsars may have a smaller value of $\mathcal{B}$, or a larger $x_{1}$,
than young ones), implying that the coefficients $a_{0,1}$ in (96) are also
functions of $\Omega$ and $\dot{\Omega}$ via their combination $T$. This also
affects the theoretical values of $a$ and $b$ in Tab. 1.

Figure 1:  The $P-\dot{P}$ diagram for all pulsars listed in the ATNF
catalogue (accessed on 2022 April 1). The colour bar represents the
theoretically predicted timing noise strength $\sigma$ for $m=1$ on a
logarithmic scale, normalised over the Crab’s value of $\sigma$. In all three
panels the common parameters are: $\mathcal{B}=2\times 10^{-}9$, $x_{1}=0.1$,
$\nu_{o,s}=2\pi/T_{o,s}$, where $T_{o}=20\,$yr and $T_{s}=1$ day. The cross
marks the Crab pulsar (PSR J0534+2200). Left: pure internal noise for
$\alpha_{\mathcal{T}}=10^{-2}$, the timing noise strength approximately scales
as $\sigma\propto\dot{P}P^{0.5}$. Centre: $\alpha_{\mathcal{T}}=10^{-2}$,
$\alpha_{\infty}=10^{-7}$. Right: pure external noise for
$\alpha_{\infty}=10^{-7}$, the timing noise strength approximately scales as
$\sigma\propto\sqrt{\dot{P}/P}\propto T^{-0.5}$. The lines of constant
magnetic field for $10^{8}-10^{14}\,$G (red, dashed) and the lines of constant
characteristic age $T/2=0.1-10^{9}\,$kyr (grey, dashed) are also indicated to
guide the eye.

##  6 Numerical tests

The Itô’s process in (4) may be integrated forward in time thanks to the
simple Euler–Maruyama scheme, that is often used to solve the
Ornstein–Uhlenbeck process (see App. B); for a finite time step of amplitude
$\Delta t>0$, we have

$\delta\mathbf{\Omega}_{t+\Delta
t}=\delta\mathbf{\Omega}_{t}+B\,\delta\mathbf{\Omega}_{t}\,\Delta
t+M\,\mathbf{X}\,\sqrt{\Delta t}+O(\Delta t^{3/2})$ (97)

where $\mathbf{X}=(\mathbf{W}_{t+\Delta t}-\mathbf{W}_{t})/\sqrt{\Delta t}$ is
a vector of i.i.d. normal random variables sampled at every time step. To give
an explicit example, for the $m=1$ system defined by (32) and (34), the above
scheme reads (the shorthand $\Omega(t)=\Omega^{t}$ is used):

$\begin{split}&\delta{\Omega}_{p}^{t+\Delta
t}=\delta{\Omega}_{p}^{t}+x_{1}\frac{\delta{\Omega}_{1}^{t}-\delta{\Omega}_{p}^{t}}{\tau}\,\Delta
t+\left[X_{p}\frac{\sigma_{\infty}}{x_{p}}-X_{1}\frac{\sigma_{\mathcal{T}}}{x_{p}}\right]\sqrt{\Delta
t}\\\ &\delta{\Omega}_{1}^{t+\Delta
t}=\delta{\Omega}_{1}^{t}-x_{p}\frac{\delta{\Omega}_{1}^{t}-\delta{\Omega}_{p}^{t}}{\tau}\,\Delta
t+X_{1}\frac{\sigma_{\mathcal{T}}}{x_{1}}\sqrt{\Delta t}\\\
&X_{p}\sim\mathcal{N}(0,1)\,,\quad X_{1}\sim\mathcal{N}(0,1)\end{split}$ (98)

where $\sigma_{\infty}$ and $\sigma_{\mathcal{T}}$ are given in (52) and (57),
respectively. In the end, the evolution in (98) will depend on six physical
parameters: $x_{1}$, $\alpha_{\infty}$, $\alpha_{\mathcal{T}}$, $\mathcal{B}$,
$\Omega$, $\dot{\Omega}$. Moreover, the scheme in (98) is consistent with the
evolution of the angular momentum residuals888 Recall that $\dot{\Omega}>0$ is
just a benchmark value for the observed spin-down rate $|\dot{\Omega}_{p}|$,
to be set from observations together with $\Omega>0$.

$\delta{L}_{t+\Delta t}=\delta{L}_{t}+X_{0}\,\sigma_{\infty}\,\sqrt{\Delta
t}=\delta{L}_{t}+X_{0}\,\alpha_{\infty}\sqrt{\Omega\,\dot{\Omega}\,\Delta
t}\,.$ (99)

The $m=2$ case in (42) and (43) is completely analogous, with
$\sigma_{\infty}$ that is still set according to (52), while the two
parameters $\sigma_{\mathcal{T}_{1}}$ and $\sigma_{\mathcal{T}_{2}}$ are given
in (58). In order to capture the features of the PSD and to ensure that the
scheme is linearly stable, the “sampling rate” $\nu_{s}=1/\Delta t$ is always
chosen to be considerably higher than the largest eigenvalue of $-B$. For
example, we choose $\Delta t\ll b^{-1}$ for $m=1$ and $\Delta
t\ll\min(b_{1}^{-1},b_{2}^{-1})$ for $m=2$. Note that $b_{1}$ and $b_{2}$ are
not eigenvalues of $B$ for $m=2$, but this condition automatically ensures
that the linear stability criterion is satisfied, see equation (A.19) of
Montoli et al. 2020.

###  6.1 Evolution of a Vela-like pulsar

In order to build some intuition about the phenomenology associated with the
stochastic process introduced in Section 1, we show the results of some simple
numerical tests for a pulsar that has the rotational parameters of the Vela
(PSR J0835-4510, see e.g. Manchester et al. 2005); we set
$\Omega_{p}(0)=70\,$rad/s, $|\dot{\Omega}_{\infty}|=10^{-10}\,$rad/s2. We
advance the equations in (98) for $20\,$yr, roughly the span of the phase-
coherent timing solution for the Vela that encompasses $\sim 10^{3}$ TOAs
obtained with the Parkes radio telescope (roughly) between 1993 and 2014
(Shannon et al., 2016).

We first consider the case of fluctuations in the external torque
($\alpha_{\mathcal{T}}=0$), see Figs. 2 and 3. In Figure 2 we study the
qualitative effect of varying the extension of the superfluid region that
participates to the internal torque. In the upper panel we have set
$x_{p}=0.1$, meaning that the star is mostly superfluid: the massive
superfluid bulk tends to act as a stabilising reservoir of angular momentum
and the resulting $\delta\Omega_{p}$ process looks almost like a white-noise
process. On the contrary, when the superfluid reservoir is limited (as in the
lower panel of Fig. 2, $x_{p}=0.99$), we are qualitatively closer to the $m=0$
case: in this degenerate limit we know that $\delta\Omega_{p}$ is a Wiener
process (i.e. pure Brownian motion, or integrated white-noise). Therefore, the
two panels show how changing $x_{p}$ results in a transition from an evolution
that visually looks like a white process to a red one. This is consistent with
the fact that the “white region” $\omega\in[\mu,\xi]$ where $P_{p}(\omega)$ is
flat reduces to a point for $x_{p}\approx 1$, see equation (62).

Another expected feature of the model with pure external noise is that the
superfluid component should somehow lag behind the non-superfluid one. In
fact, the fluctuations in the external torque act directly on the observable
component, that, in turn, exchanges angular momentum with the internal one via
the mutual friction. We check this in Fig. 3, that is completely analogous to
Fig. 2 apart from the fact that the dimensionless friction parameter
$\mathcal{B}$ is taken to be one order of magnitude smaller. Again, we can
visually recognise that the evolution reddens by increasing the value of
$x_{p}$ but now the long timescale $1/b\approx 800\,$d (common to both the
upper and lower panels) makes more evident the fact that the evolution of
$\delta\Omega_{1}$ lags behind $\delta\Omega_{p}$. This is nothing but the
stochastic version of the deterministic steady state solution in (30) or (41),
where $\Omega_{1}$ assumes the same value of $\Omega_{p}$ exactly after a time
$b^{-1}=\tau/x_{p}$. Since the superfluid lags behind, its evolution is also
considerably smoother, as it is the result of a sum of many random kicks
imparted by the normal component over the timescale $b^{-1}$.

We briefly comment also on the case with pure internal noise, see Fig. 4,
where we have set $\alpha_{\mathcal{T}}=0.1,\alpha_{\infty}=0$. Again, the
general trend of reddening in the limit $x_{p}\rightarrow 1$ still holds true
but the evolution of the superfluid is qualitatively different, as it does not
lag behind the observable one and it is not smoother. In fact, for pure
internal noise the spectral properties of $\delta\Omega_{p}$ and
$\delta\Omega_{1}$ are identical, since one process is the (rescaled) mirror
image of the other, see (22):

$\delta{\Omega}_{p}(t)=-\dfrac{x_{1}}{x_{p}}\delta{\Omega}_{1}(t)\qquad\quad(\text{i.e.}\,\,\,\delta
L=0)\,.$ (100)

This is particularly evident in the central and right panels of Fig. 4. Let us
stress that in all the three figures 2, 3 and 4, it may look like the lag is
reversed when the pink curve is above the blue one. However,
$\delta\Omega_{1}-\delta\Omega_{p}$, is not the lag, but rather the lag
residual with respect to its deterministic value; given our initial condition,
the deterministic lag is the constant steady state in (30).

Figure 2:  A typical realisation of the processes $\delta\Omega_{p}$ (pink)
and $\delta\Omega_{1}$ (blue) for the two component model ($m=1$) with pure
external noise ($\alpha_{\infty}=10^{-5}$, $\alpha_{\mathcal{T}}=0$), over a
time span of 20 years. The rotational parameters are the ones of the Vela
pulsar, $\Omega\approx 70\,$rad/s, $\dot{\Omega}\approx 10^{-10}\,$rad/s2. In
both cases $\mathcal{B}=10^{-9}$ and the only difference between the two
panels is the value of the superfluid ratio: $x_{p}=0.1$ in the upper panel
(the star is mostly superfluid and the resulting relaxation timescale is
$\tau\approx 8\,$d), $x_{p}=0.99$ in the lower panel (the superfluid is
confined in the inner crust, $\tau\approx 80\,$d). The fluctuations in the
external torque drive the observable component, while the superfluid tends to
follow it with some delay (this is even more evident in Fig. 3).  Figure 3:  A
single realisation of the processes $\delta\Omega_{p}$ (pink) and
$\delta\Omega_{1}$ (dark blue) for the $m=1$ model with fluctuations only in
the external torque ($\alpha_{\mathcal{T}}=0$, $\Omega\approx 70\,$rad/s,
$\dot{\Omega}\approx 10^{-10}\,$rad/s2). The only difference with respect to
Fig. 2 is that now $\mathcal{B}=10^{-10}$. Upper panel: $x_{p}=0.1$ (the star
is mostly superfluid and the resulting relaxation timescale is $\tau\approx
80\,$d). Lower panel: $x_{p}=0.99$ (the superfluid is confined in the inner
crust, $\tau\approx 800\,$d). Since the driving fluctuations are in the
external torque, the evolution of the superfluid angular velocity residual
$\delta\Omega_{1}$ lags behind the one of the observable component and is
considerably smoother, as expected.  Figure 4:  A single realisation of the
processes $\delta\Omega_{p}$ (pink) and $\delta\Omega_{1}$ (blue) for the
$m=1$ model with fluctuations only in the internal torque
($\alpha_{\infty}=0$, $\alpha_{\mathcal{T}}=0.1$, $\Omega\approx 70\,$rad/s,
$\dot{\Omega}\approx 10^{-10}\,$rad/s2). Right: $\mathcal{B}=10^{-9}$,
$x_{p}=0.1$ ($\tau\approx 8\,$d). Centre: $\mathcal{B}=10^{-10}$, $x_{p}=0.1$
($\tau\approx 80\,$d). Left: only in this panel the process $\delta\Omega_{1}$
has been multiplied by $x_{1}=0.01$ in order to fit it within the canvas
($\mathcal{B}=10^{-10}$, $x_{p}=0.99$, $\tau\approx 800\,$d), see equation
(100).

Finally, it is interesting to to see if the stochastic model can reproduce, at
least qualitatively, the 20 years of timing residual data for the Vela pulsar
given in Fig. 1 of Shannon et al. (2016). Since our model does not include the
possible effect of glitches, we have to compare our residuals with the ones
reported in panel (b) therein, where the residual arrival times, measured in
units of time, are fitted for the glitches. This removes the effect of the
glitch recovery from the displayed residuals of Shannon et al. (2016); over
$20\,$yr the quasi-periodic oscillation of the timing residuals has an
amplitude of $\sim 10\,$s. We find that the same features of $\delta a(t)$ can
be reproduced by evolving our model with pure external noise, as shown in Fig.
5.

Following a widespread observational procedure (Cordes & Downs, 1985), the
curves $\delta a(t)$ in Fig. 5 do not start from the origin of the axes
because the best-fit quadratic polynomial has been removed. This only modifies
their spectrum in the inaccessible infrared region, while their spectral
properties for $\omega\gtrsim\nu_{o}$ are left practically unchanged. From a
purely visual point of view, this also eases the comparison of the simulated
residuals with the ones of Shannon et al. (2016), cf. the Vela residuals in
Fig. 2 of Lower et al. (2020). More precise analysis should be done by
contrasting the empirical PSD extracted from timing data with the theoretical
one, which possibly allows falsifying the assumptions of our theoretical model
(e.g., the assumption of uncorrelated white noise terms in the Langevin
equation). For the moment, we limit ourselves to notice that, for a very low
value of $x_{p}$ (in the upper panel we used the unphysical value
$x_{p}=0.01$), the overall evolution is qualitatively similar, but some small
scale decorations (that are absent in the observed timing residuals, the black
dotted line in Fig. 5) are present. This may be a by product of the fact that
observations have a lower resolution than our simulations.

Figure 5:  Four typical realisations of the timing residuals $\delta a(t)$ for
the $m=1$ model with fluctuations only in the external torque
($\alpha_{\mathcal{T}}=0$, $\Omega\approx 70\,$rad/s, $\dot{\Omega}\approx
10^{-10}\,$rad/s2, $\mathcal{B}=10^{-9}$). Upper panel: $x_{p}=0.01$. Lower
panel: $x_{p}=0.99$. The dashed black line represents the timing residuals
from 20 years of timing of the Vela pulsar extracted from Fig. 1 of Shannon et
al. (2016).

###  6.2 Power spectral density

Even though our simple derivation of the PSD $P_{p}(\omega)$ is based on the
well known trick of computing the Fourier transform of the Langevin equation,
this procedure is far from being rigorous in our case (Priestley, 1965) and,
more importantly, it is not obvious how well the final $P_{p}(\omega)$ may
approximate the PSD obtained from a finite sample of the signal
$\delta\Omega_{p}$ over a time span $T_{o}$ by means of the discrete Fourier
transform (DFT).

In fact, the procedure we have used to obtain $P_{p}(\omega)$ is typically
used for mean-reverting multivalued Ornstein-Uhlenbeck processes for which,
however, all eigenvalues of $B$ are strictly negative (Gardiner, 1994; Singh
et al., 2018), see also App. B.

On the other hand, our matrix $B$ is singular (it has at least one null
eigenvector). On top of this theoretical problem, we also have the well known
fact that windowing and sampling of a continuous signal, in this case
$\delta\Omega_{p}$ leads to a distortion of the power spectrum , see e.g. (van
der Klis, 1989). Therefore, it is better to check numerically the validity of
our result for the PSD. It will be sufficient to do so for the $m=1$ model
since, morally, the case $m=2$ is identical: the PSD for both $m=1$ and $m=2$
have the same kind of infrared divergence for $\omega\rightarrow 0$, that
ultimately results from long-lived correlations – the autocovariance is not
integrable – typical of some non-stationary processes, including the Wiener
process, see (Kasdin, 1995).

Given the signal $\delta\Omega_{p}$, sampled at times $t_{l}=lT_{0}/N$ for
$l=0,...,N-1$, our exact definition of DFT is,

$\delta\tilde{\Omega}_{p}(\omega_{k})=\,\sum_{l=0}^{N-1}\,\delta{\Omega}_{p}(t_{i})\,e^{-2\pi
i\,lk/N}\qquad\quad\omega_{k}=k\nu_{s}\,,$ (101)

where the sampling frequency is $\nu_{s}=2\pi N/T_{0}$ and the first non-DC
frequency is $\nu_{0}=\omega_{1}=2\pi/T_{0}$. In this way, the “negative
frequency” part of the spectrum is contained in the second half of the DFT
vector, and we can restrict our attention to
$\delta\tilde{\Omega}_{p}(\omega_{k})$ for $k=0,...,N/2$. Clearly, in contrast
with he continuum transform $\delta{\Omega}_{p}(\omega)$,
$\delta\tilde{\Omega}_{p}$ has the same physical dimensions of an angular
velocity. Therefore, to obtain something that well approximates
$P_{p}(\omega)$ we have to rescale the discrete power spectrum
$p_{k}=|\delta\tilde{\Omega}_{p}(\omega_{k})|^{2}$ as

$P_{p}(\omega=\omega_{k})\approx\dfrac{T_{0}}{N^{2}}\,p_{k}\,.$ (102)

We apply these prescriptions for a benchmark case in Fig. 6, where we plot our
theoretical PSD (35) for the $m=1$ model together with the one obtained from
the DFT of simulated data over a time span $T_{0}=20$yr (the scaling in (102)
is used). The parameters used have been selected to emphasise the presence of
the white region (both internal and external fluctuations are switched on). We
see that the theoretical result $P_{p}(\omega)$ approximates well the spectrum
obtained by averaging $p_{k}$ over an ensemble of $50$ realisations of
$\delta\Omega_{p}(t)$. The PSD extracted from a single realisation (that is
always the case for real timing data) is considerably more noisy, but the
overall shape and scaling is the same. However, the PSD extracted from a
single realisation shows that it can be difficult to identify exactly the
whitened region in real data, especially if its extent is less than one order
of magnitude in $\omega$.

We apply the same strategy to the timing residuals $\delta a$, obtained via
the prescription (85). The corresponding PSD is shown in the lower panel of
Figure 6: inside the whitened region the spectral index is -2, outside -4.

Figure 6:  Upper panel: comparison between the PSD for the angular velocity in
equation (35), the black dashed line, and the PSD obtained from the discrete
Fourier transform of simulated data (orange and blue curves) obtained via
equation (102). The blue curve has been extracted from a single realisation,
the orange one obtained after taking the ensemble average of $p_{k}$ in (102)
over 50 realisations. The parameters used are: $\alpha_{\infty}=10^{-3}$,
$\alpha_{\mathcal{T}}=0.1$, $x_{p}=0.03$, $\mathcal{B}=10^{-9}$,
$\Omega\approx 70\,$rad/s, $\dot{\Omega}\approx 10^{-10}\,$rad/s2. Given these
parameters, the relaxation time is $\tau\approx 2.5\,$d. The two corner
frequencies that define the whitened region $\mu<\omega<\xi$ are indicated by
the green dashed lines. The vertical dotted lines correspond to the the first
non-DC frequency $\nu_{o}=2\pi/T_{o}$, for $T_{0}=20\,$ yr, and to a sampling
frequency of 1 data point per day, $\nu_{s}=2\pi/$day. Lower panel: the
corresponding PSD for the timing residuals $\delta a$.

##  7 Conclusions

We have studied the possibility to describe the timing noise by means of a
simple stochastic model for the fluctuations in both the external and internal
torques in a neutron star. The present framework offers a guideline to
interpret the observed features of timing noise in both single pulsars and
across pulsar population by including an arbitrary number of internal
components. The framework allows one to model the phenomenology associated to
different physical situations: when the source of fluctuations is the mutual
friction associated to the motion of superfluid vortices, i.e. fluctuations in
the internal torques between the components, and when the source of
fluctuations is associated to the external braking torque.

We find that the scaling of the timing noise strength
$\sigma\sim\Omega^{a}\dot{\Omega}^{b}$ is different for these two physical
scenarios, see Table 1. Moreover, assuming that fluctuations in both the
internal and external torques are active, the global scaling of $\sigma$ can
not be written just in terms of the two exponents $a$ and $b$, as it is
evident in Fig. 1. In this case, one should consider the full result in (89)
or, at least, the sum of the two leading terms in (96).

Unfortunately, analysis that aim to extrapolate $a$ and $b$ from observations
across the pulsar population is difficult and results reported in the
literature agree only at the qualitative level, see Tab. 1. However, taking as
a reference the recent estimates of $a$ and $b$ by Lower et al. (2020), we see
that our stochastic model could reproduce the timing noise across the pulsar
population only if at least one of the three possibilities is met (further
work is needed to assess which is the most realistic possibility): injected
noise is coloured rather than white (as in the the model for internal
turbulence of Melatos & Link 2014); fluctuations in the internal and external
torques are both present; the fluctuation strength varies across the
population because it is correlated with the pulsar’s age. Regarding this last
possibility, the physical parameters can indeed evolve as the neutron star
cools down (e.g. the phenomenological friction parameter may be higher in
younger pulsars because of their higher internal temperature), so that the
final scaling $\sigma\propto\Omega^{a}\dot{\Omega}^{b}$ is also affected.
However, in the case of pure external torque we have no dependence on the
friction parameter in the leading term of $\sigma$, meaning that the overall
scaling is less likely to be affected by physical differences due to age.

Similarly, also the power spectra of the angular velocity, phase and TOAs
residuals are dependent on the different physical parameter of the neutron
star. Therefore, using these analytical results to model the timing noise from
the timing data data may help in constraining some global properties of
pulsars, like the body averaged friction coupling or the moment of inertia of
the components. This provides a way to probe the interior of neutron stars by
means of timing noise that is analogous to what is done for pulsar glitches,
where observations of the spin up in tandem with a Bayesian analysis based on
minimal (deterministic) models of the kind discussed here, allows to extract
estimates of the coupling parameters and moments of inertia (Montoli et al.,
2020).

We have also tentatively applied our model to estimate the timing noise
strength $\sigma$ across the population of known pulsars. Under the working
hypothesis that the physical parameters have no dependence the pulsar’s age,
we have found a general trend where the younger pulsars have much larger
timing noise than the older ones (Fig. 1), in accordance with observations
(Arzoumanian et al., 1994; Shannon & Cordes, 2010). A similar pattern has been
observed for the glitch activity, see e.g. (Fuentes et al., 2017; Basu et al.,
2021). Since the modelling proposed here is closely related to the kind of
minimal models used in pulsar glitch studies, it will be interesting to extend
the present framework beyond the hypothesis of pure white noise, to include
also the effect of glitches.

Finally, we have reproduced the qualitative behaviour and amplitude of the
timing residuals extracted from 20 years of almost phase-coherent monitoring
of the Vela pulsar (Shannon et al., 2016). Our limited goal was to start to
probe the parameter space of our Langevin model, and we have found that the
limit in which the external noise is dominant may be a good starting point for
a minimal (in the sense that it contains the least amount of free physical
parameters) model for the timing noise in the Vela (see Fig. 5). Of course,
such a minimal model will have to be contrasted with real data in a detailed
analysis.

##  8 Acknowledgements

Partial support comes from PHAROS, COST Action CA16214. M.A. and B.H.
acknowledge support from the Polish National Science Centre grant SONATA BIS
2015/18/E/ST9/00577. A.B. acknowledges the support from the UK Science and
Technology Facilities Council(STFC). A consolidated grant from STFC supports
the pulsar research at the Jodrell Bank Centre for Astrophysics.

##  9 Data Availability

The numerical results and the code underlying this work are available from
M.A., upon reasonable request.

## References

  * Alpar et al. (1986) Alpar M. A., Nandkumar R., Pines D., 1986, ApJ, 311, 197
  * Antonelli & Pizzochero (2017) Antonelli M., Pizzochero P. M., 2017, MNRAS, 464, 721
  * Antonelli et al. (2018) Antonelli M., Montoli A., Pizzochero P. M., 2018, MNRAS, 475, 5403
  * Arzoumanian et al. (1994) Arzoumanian Z., Nice D. J., Taylor J. H., Thorsett S. E., 1994, ApJ, 422, 671
  * Ashton et al. (2019) Ashton G., Lasky P. D., Graber V., Palfreyman J., 2019, Nature Astronomy, 3, 1143
  * Basu et al. (2021) Basu A., et al., 2021, MNRAS, 510, 4049
  * Baym et al. (1969) Baym G., Pethick C., Pines D., Ruderman M., 1969, Nature, 224, 872
  * Boynton et al. (1972) Boynton P. E., Groth E. J., Hutchinson D. P., Nanos G. P. J., Partridge R. B., Wilkinson D. T., 1972, ApJ, 175, 217
  * Celora et al. (2020) Celora T., Khomenko V., Antonelli M., Haskell B., 2020, MNRAS, 496, 5564
  * Coles et al. (2011) Coles W., Hobbs G., Champion D. J., Manchester R. N., Verbiest J. P. W., 2011, MNRAS, 418, 561
  * Cordes & Downs (1985) Cordes J. M., Downs G. S., 1985, ApJS, 59, 343
  * Cordes & Helfand (1980) Cordes J. M., Helfand D. J., 1980, ApJ, 239, 640
  * D’Alessandro (1996) D’Alessandro F., 1996, Ap&SS, 246, 73
  * Easson (1979) Easson I., 1979, ApJ, 228, 257
  * Fuentes et al. (2017) Fuentes J. R., Espinoza C. M., Reisenegger A., Shaw B., Stappers B. W., Lyne A. G., 2017, A&A, 608, A131
  * Gardiner (1994) Gardiner C. W., 1994, Handbook of stochastic methods for physics, chemistry and the natural sciences
  * Gavassino et al. (2020) Gavassino L., Antonelli M., Pizzochero P. M., Haskell B., 2020, MNRAS, 494, 3562
  * Graber et al. (2018) Graber V., Cumming A., Andersson N., 2018, ApJ, 865, 23
  * Groth (1975) Groth E. J., 1975, ApJS, 29, 453
  * Gügercinoǧlu & Alpar (2017) Gügercinoǧlu E., Alpar M. A., 2017, MNRAS, 471, 4827
  * Haskell et al. (2012) Haskell B., Pizzochero P. M., Sidery T., 2012, MNRAS, 420, 658
  * Heintzmann et al. (1973) Heintzmann H., Kundt W., Schrufer E., 1973, A&A, 27, 45
  * Hobbs et al. (2010) Hobbs G., Lyne A. G., Kramer M., 2010, MNRAS, 402, 1027
  * Jankowski et al. (2017) Jankowski F., van Straten W., Keane E. F., Bailes M., Barr E. D., Johnston S., Kerr M., 2017, MNRAS, 473, 4436
  * Janssen & Stappers (2006) Janssen G. H., Stappers B. W., 2006, A&A, 457, 611
  * Johnston & Galloway (1999) Johnston S., Galloway D., 1999, MNRAS, 306, L50
  * Jones (1990) Jones P. B., 1990, MNRAS, 246, 364
  * Kasdin (1995) Kasdin N., 1995, Proceedings of the IEEE, 83, 802
  * Lam et al. (2019) Lam M. T., et al., 2019, ApJ, 872, 193
  * Liu et al. (2019) Liu X. J., Keith M. J., Bassa C. G., Stappers B. W., 2019, MNRAS, 488, 2190
  * Lorimer & Kramer (2005) Lorimer D. R., Kramer M., 2005, Handbook of Pulsar Astronomy. Cambridge Observing Handbooks for Research Astronomers, Cambridge University Press
  * Lower et al. (2020) Lower M. E., et al., 2020, MNRAS, 494, 228
  * Lower et al. (2021) Lower M. E., et al., 2021, MNRAS, 508, 3251
  * Lyne et al. (2010) Lyne A., Hobbs G., Kramer M., Stairs I., Stappers B., 2010, Science, 329, 408
  * Manchester et al. (2005) Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, Aj, 129, 1993
  * Melatos & Link (2014) Melatos A., Link B., 2014, MNRAS, 437, 21
  * Meyers et al. (2021a) Meyers P. M., Melatos A., O’Neill N. J., 2021a, MNRAS, 502, 3113
  * Meyers et al. (2021b) Meyers P. M., O’Neill N. J., Melatos A., Evans R. J., 2021b, MNRAS, 506, 3349
  * Montoli et al. (2020) Montoli A., Antonelli M., Magistrelli F., Pizzochero P. M., 2020, A&A, 642, A223
  * Parthasarathy et al. (2019) Parthasarathy A., et al., 2019, MNRAS, 489, 3810
  * Pizzochero et al. (2020) Pizzochero P. M., Montoli A., Antonelli M., 2020, A&A, 636, A101
  * Priestley (1965) Priestley M. B., 1965, Journal of the Royal Statistical Society. Series B (Methodological), 27, 204
  * Rice (1944) Rice S. O., 1944, Bell System Technical Journal, 23, 282
  * Rickett (1975) Rickett B. J., 1975, ApJ, 197, 185
  * Shannon & Cordes (2010) Shannon R. M., Cordes J. M., 2010, ApJ, 725, 1607
  * Shannon et al. (2016) Shannon R. M., Lentati L. T., Kerr M., Johnston S., Hobbs G., Manchester R. N., 2016, MNRAS, 459, 3104
  * Shaw et al. (2022) Shaw B., et al., 2022, MNRAS, 513, 5861
  * Singh et al. (2018) Singh R., Ghosh D., Adhikari R., 2018, Phys. Rev. E, 98, 012136
  * Sourie & Chamel (2020) Sourie A., Chamel N., 2020, MNRAS, 493, 382
  * Sourie et al. (2017) Sourie A., Chamel N., Novak J., Oertel M., 2017, MNRAS, 464, 4641
  * Taylor et al. (1975) Taylor J. H., Manchester R. N., Huguenin G. R., 1975, ApJ, 195, 513
  * van Haasteren & Levin (2013) van Haasteren R., Levin Y., 2013, MNRAS, 428, 1147
  * van der Klis (1989) van der Klis M., 1989, Fourier Techniques in X-Ray Timing. Springer Netherlands, Dordrecht, pp 27–69, doi:10.1007/978-94-009-2273-0_3

##  Appendix A Extension to a generic number of noise terms and spin-down
torques

The $m=1$ model in (27) can be slightly generalised to account for the
possibility that part of the total external torque $\mathcal{T}_{\infty}$ acts
directly on the superfluid component,

$\begin{split}&x_{p}\dot{\Omega}_{p}=-\mathcal{T}+\mathcal{T}^{p}_{\infty}\\\
&x_{1}\dot{\Omega}_{1}=\mathcal{T}+\mathcal{T}^{1}_{\infty}\,,\end{split}$
(103)

where
$\mathcal{T}_{\infty}=\mathcal{T}^{p}_{\infty}+\mathcal{T}^{1}_{\infty}$. This
allows to model the possibility of having an extra spin-down torque
$\mathcal{T}^{1}_{\infty}$ associated with the emission of continuous
gravitational waves by oscillation modes in the superfluid component, as
proposed in Meyers et al. (2021a); Meyers et al. (2021b).

The model in Sec. 2 can account for $m$ components and (up to) $m+1$
independent additive noise terms ($m$ fluctuating torques due to the friction
of each superfluid layer with the unique normal component plus the external
torque). However, the matrix form of the equations (6) is more general: if $M$
is a rectangular matrix with $(m+1)$ rows, an arbitrary number of noise terms
can be included and the PSD is still the one in (21). This allows to model the
situation where part of the spin-down is partially due to gravitational wave
emission from the superfluid layers.

Assuming that each of the $m$ superfluid components also have an associated
fluctuating external torque, so that the noise terms are up to $(m+1)+m$, it
may also be convenient to work with two matrices, $M_{\mathrm{e}}$ and
$M_{\mathrm{i}}$ that define the external ($\mathrm{e}$) and internal
($\mathrm{i}$) fluctuations, namely

$\begin{split}&\dot{\mathbf{\Omega}}_{t}=B\,\mathbf{\Omega}_{t}+\mathbf{A}_{t}+M_{\mathrm{e}}\dot{\mathbf{w}}^{\mathrm{e}}_{t}+M_{\mathrm{i}}\dot{\mathbf{w}}^{\mathrm{i}}_{t}\\\
&\langle\dot{\mathbf{w}}^{\mathrm{i}}_{t}\rangle=\langle\dot{\mathbf{w}}^{\mathrm{e}}_{t}\rangle=0\\\
&\langle\dot{\mathbf{w}}^{\mathrm{i}}_{t}\dot{\mathbf{w}}^{\mathrm{i\top}}_{s}\rangle=\langle\dot{\mathbf{w}}^{\mathrm{e}}_{t}\dot{\mathbf{w}}^{\mathrm{e\top}}_{s}\rangle=\delta(t-s)\mathbb{I}\\\
&\langle\dot{\mathbf{w}}^{\mathrm{i}}_{t}\dot{\mathbf{w}}^{\mathrm{e\top}}_{s}\rangle=\langle\dot{\mathbf{w}}^{\mathrm{e}}_{t}\dot{\mathbf{w}}^{\mathrm{i\top}}_{s}\rangle=0\end{split}$
(104)

where $\dot{\mathbf{w}}^{\mathrm{i}}_{t}$ and
$\dot{\mathbf{w}}^{\mathrm{e}}_{t}$ are two independent sets of white noise
processes. Using $M_{\mathrm{e,i}}$ is not strictly necessary, but it allows
us to better keep track of the different nature of the torques in the
analysis, as we briefly outline.

To guarantee that (104) is consistent with the constraint (3), we have to
require the analogous of the action-reaction property (8),

$B^{\top}\mathbf{x}=0\qquad\quad
M_{\mathrm{i}}^{\top}\mathbf{x}=0\qquad\quad\text{(action-reaction property)}$
(105)

as well as the analogous of (9): in this case the fluctuation
$\eta_{t}^{\infty}$ of the total external torque $\mathcal{T}_{\infty}$ is
given by

$\mathbf{x}^{\top}\cdot\mathbf{A}_{t}=\mathcal{T}_{\infty}\qquad\mathbf{x}^{\top}M_{\mathrm{e}}\dot{\mathbf{w}}^{\mathrm{e}}_{t}=\eta_{t}^{\infty}$
(106)

From (104) we have that the Fourier components of the angular velocity
residuals are given by

$\delta\mathbf{\Omega}_{\omega}=(i\omega\,\mathbb{I}-B)^{-1}(M_{\mathrm{e}}\dot{\mathbf{w}}^{\mathrm{e}}_{\omega}+M_{\mathrm{i}}\dot{\mathbf{w}}^{\mathrm{i}}_{\omega})\,,$
(107)

implying that the PSD is

$P_{p}(\omega)\propto\left[(i\omega\,\mathbb{I}-B)^{-1}(M_{\mathrm{i}}M_{\mathrm{i}}^{\top}+M_{\mathrm{e}}M_{\mathrm{e}}^{\top})(-i\omega\,\mathbb{I}-B)^{-1}\right]_{pp}\,.$
(108)

It is interesting to note that the spectrum associated with the total angular
momentum still has a purely Brownian spectrum also in this more general
situation, just its overall power is due to all the external fluctuations:
thanks to (105), the analogous of (24) is

$P_{L}(\omega)=\frac{\mathbf{x}^{\top}M_{\mathrm{e}}M_{\mathrm{e}}^{\top}\mathbf{x}}{\omega^{2}}$
(109)

The Euler–Maruyama scheme for the system in (104) is completely analogous to
the one in (97):

$\begin{split}&\delta\mathbf{\Omega}_{t+\Delta
t}\approx\delta\mathbf{\Omega}_{t}+B\,\delta\mathbf{\Omega}_{t}\,\Delta
t+\left(M_{\mathrm{e}}\,\mathbf{X}^{\mathrm{e}}+M_{\mathrm{i}}\,\mathbf{X}^{\mathrm{i}}\right)\sqrt{\Delta
t}\\\ &{X}^{\mathrm{e}}_{j}\sim\mathcal{N}(0,1)\qquad\forall j=0,...,m\\\
&{X}^{\mathrm{i}}_{j}\sim\mathcal{N}(0,1)\qquad\forall
j=1,...,m\qquad({X}^{\mathrm{i}}_{0}=0)\,\end{split}$ (110)

that is consistent with the evolution of the total angular momentum

$\delta L_{t+\Delta t}=\delta
L_{t}+\mathbf{x}^{\top}M_{\mathrm{e}}\,\mathbf{X}^{\mathrm{e}}\sqrt{\Delta
t}+O(\Delta t^{3/2})\,.$ (111)

###  A.1 Comparison with Meyers et al. (2021a)

The $m=1$ model of Meyers et al. (2021a) is a particular realisation of the
one in (104), where $M_{\mathrm{i}}=0$. In fact, this model implements $2$
fluctuating terms, and has all the properties discussed for the system in
(104), provided that the coupling parameters appearing in $B$ are such that
the action-reaction property (8) is respected. Therefore, it is interesting to
discuss the general $m=1$ case, which contains both our model of Sec. 2.2 and
the one of Meyers et al. (2021a) as special cases. Let us consider a system
defined by

$\mathbf{\Omega}=\begin{pmatrix}\dot{\Omega}_{p}\\\
\dot{\Omega}_{1}\end{pmatrix}\quad B=\begin{pmatrix}-x_{1}/\tau&x_{1}/\tau\\\
x_{p}/\tau&-x_{p}/\tau\end{pmatrix}\quad\mathbf{A}=\begin{pmatrix}-A_{p}\\\
-A_{1}\end{pmatrix}$ (112)

for some spin-down torques $A_{p}>0$ and $A_{1}>0$, namely999 The
deterministic steady-state is
$\dot{\Omega}_{p}=\dot{\Omega}_{1}=-|\dot{\Omega}_{\infty}|=-A_{p}-A_{1}$, and
steady lag is $\Omega_{1}-\Omega_{p}=(x_{1}A_{p}+x_{p}A_{1})/(\tau
x_{1}x_{p})\,$.

$\begin{split}&x_{p}\dot{\Omega}_{p}=\frac{x_{p}x_{1}}{\tau}(\Omega_{1}-\Omega_{p})-A_{p}\\\
&x_{1}\dot{\Omega}_{1}\,=\frac{x_{p}x_{1}}{\tau}(\Omega_{p}-\Omega_{1})-A_{1}\end{split}$
(113)

In this simple case additive noise is conveniently introduced via the
rectangular matrix, cf. (34),

$M=\begin{pmatrix}\sigma_{p}/x_{p}&0&-\sigma_{\mathcal{T}}/x_{p}\\\
0&\sigma_{1}/x_{1}&\sigma_{\mathcal{T}}/x_{1}\,,\end{pmatrix}$ (114)

where $\sigma_{p,1}$ define the strength of the fluctuations associated with
the torques $A_{p,1}$. Consistently with (109), we have that only the external
fluctuations define the strength of $P_{L}$, and it is easy to find that
$P_{L}(\omega)=(\sigma_{p}^{2}+\sigma_{1}^{2})/\omega^{2}$. The PSD reads, cf.
(35),

$P_{p}(\omega)=\dfrac{(\sigma_{\mathcal{T}}^{2}+\sigma_{p}^{2})\tau^{2}\omega^{2}+x_{p}^{2}(\sigma_{p}^{2}+\sigma_{1}^{2})}{x_{p}^{2}(\tau^{2}\omega^{4}+\omega^{2})}\,.$
(115)

The corner frequency $\xi=1/\tau$ is left unchanged by the presence of the
extra fluctuating torque, while for $\mu$ we have, cf. (37),

$\mu^{2}=\dfrac{x_{p}^{2}\,(\sigma_{p}^{2}+\sigma_{1}^{2})}{\tau^{2}\,(\sigma_{p}^{2}+\sigma_{\mathcal{T}}^{2})}\,.$
(116)

For $\sigma_{1}=0$ we recover the model in Sec. 2.2 with
$\sigma_{p}=\sigma_{\infty}$. The ordering $0<\mu<\xi$ is maintained either
when

$\sigma_{1}^{2}\leq\sigma_{\mathcal{T}}^{2}\quad\text{and}\quad 0<x_{p}<1$
(117)

or

$\sigma_{\mathcal{T}}^{2}\leq\sigma_{1}^{2}\quad\text{and}\quad
x_{p}^{2}<\dfrac{\sigma_{p}^{2}+\sigma_{\mathcal{T}}^{2}}{\sigma_{p}^{2}+\sigma_{1}^{2}}\,.$
(118)

Otherwise, namely for

$\sigma_{\mathcal{T}}^{2}\leq\sigma_{1}^{2}\quad\text{and}\quad\dfrac{\sigma_{p}^{2}+\sigma_{\mathcal{T}}^{2}}{\sigma_{p}^{2}+\sigma_{1}^{2}}<x_{p}^{2}<1\,,$
(119)

the ordering $0<\xi<\mu$ is realised and the PSD behaves as
$P_{p}\sim\omega^{-4}$ in the frequency range $\xi<\omega<\mu$. For pure
external noise ($\sigma_{\mathcal{T}}=0$) we recover the model of Meyers et
al. (2021a), see in particular App. A of Meyers et al. 2021b: both orderings
of the corner frequencies are possible depending on the value of $x_{p}$, as
given in (118) and (119). It is interesting to note that in the case in which
the fluctuations in the internal torque are negligible and
$\sigma_{p}\ll\sigma_{1}$ the only possible ordering is $0<\xi<\mu$, whatever
the value of $x_{p}$.

##  Appendix B Pure internal noise as an Ornstein-Uhlenbeck process

For pure internal noise (i.e., $\sigma^{2}_{\infty}=0$), our model for
$\delta\Omega^{p}_{t}$ and generic $m>0$ is a multivariate Ornstein-Uhlenbeck
process (see, e.g., Singh et al. (2018) for a contained introduction to the
multivariate Ornstein-Uhlenbeck process). The Ornstein-Uhlenbeck process is
the only nontrivial process that admits a stationary distribution and is both
Gaussian and Markovian. Over time, the process tends to drift towards its
mean, meaning that the mean acts as an equilibrium for the process, a property
often referred to as mean-reverting (Gardiner, 1994; Singh et al., 2018).

For $\sigma_{\infty}\neq 0$ our process $\delta\Omega_{p}$ contains a Brownian
component that is not mean-reverting and not (wide sense) stationary.
Therefore, the analytical estimation of the PSD is a delicate task (Priestley,
1965): our results are formal and should be checked with the aid of numerical
simulations. In particular, we can not directly use the Wiener-Khinchin
theorem to estimate the PSD (Priestley, 1965). However, we can analytically
check that our formal results are, at least, consistent with the known
properties of a mean-reverting Ornstein-Uhlenbeck process when the limit
$\sigma_{\infty}=0$ is taken.

###  B.1 Angular velocity: Ornstein-Uhlenbeck process

For simplicity, consider the $m=1$ case with $\sigma_{\infty}=0$: by using
(32) and (34) into the general solution (12), we can obtain the explicit
solution for the $p$ component:

$\delta\Omega_{p}(t)=\frac{\sigma_{\mathcal{T}}}{x_{p}}\,e^{-t/\tau}\int_{0}^{t}e^{z/\tau}dW_{z}\,,$
(120)

where $W_{z}$ is a standard zero-mean one-dimensional Wiener process. This
expression is consistent with the mean-reverting zero-drift Ornstein-Uhlenbeck
process defined by the Langevin equation

$\delta\dot{\Omega}_{p}(t)=-\dfrac{1}{\tau}\,\delta\Omega_{p}(t)+\frac{\sigma_{\mathcal{T}}}{x_{p}}\,\dot{W}_{z}\,,$
(121)

where $\dot{W}_{z}$ is the standard delta-correlated white-noise process. The
mean-reversion property is guaranteed by the fact that the relaxation time is
positive, $\tau>0$. The mean and covariance are

$\begin{split}\langle\delta\Omega_{p}(t)\rangle&=\delta\Omega_{p}(0)\,e^{-t/\tau}\\\
c(t,u)&=\langle[\delta\Omega_{p}(t)-\langle\delta\Omega_{p}(t)\rangle][\delta\Omega_{p}(u)-\langle\delta\Omega_{p}(u)\rangle]\rangle\\\
&=\dfrac{\tau\,\sigma_{\mathcal{T}}^{2}}{2x_{p}^{2}}\left[e^{-|t-u|/\tau}-e^{-(t+u)/\tau}\right]\,.\end{split}$
(122)

The first expression is nothing but the mean-reverting property (in the main
text we imposed $\delta\Omega_{p}(0)=0$), while the covariance $c(t,u)$ is
obtained from (14). This process admits an asymptotic stationary state, where
the memory of the initial condition is lost. Setting $u=\Delta t+t$, in the
double limit $t\gg\tau$ and $t\gg\Delta t$ we have that $c(t,t+\Delta t)$ does
not depend on $t$, implying that the process is wide-sense stationary:

$\langle\delta\Omega_{p}(t)\rangle=0\qquad\qquad c(\Delta
t)=\dfrac{\tau\,\sigma_{\mathcal{T}}^{2}}{2x_{p}^{2}}\,e^{-|\Delta
t|/\tau}\,.$ (123)

Therefore, the PSD calculated via the Fourier transform of the Langevin
equation (121) coincides with the one estimated via the Wiener-Khinchin
theorem applied to $c(\Delta t)$:

$P_{p}(\omega)=\int_{-\infty}^{\infty}\\!\\!d\Delta t\,\,e^{-i\omega\Delta
t}\,c(\Delta
t)=\dfrac{\sigma_{\mathcal{T}}^{2}\,\tau^{2}}{x_{p}^{2}(1+\omega^{2}\tau^{2})}\,,$
(124)

that is exactly the $\sigma_{\infty}=0$ limit of the PSD in (35). Note that
the PSD obtained by Fourier transforming the Langevin equation gives the
correct result for the stationary limit of the Ornstein-Uhlenbeck process.
Moreover, the average power in the stationary state is

$\int_{-\infty}^{\infty}\\!\dfrac{d\omega}{2\pi}\,P_{p}(\omega)=c(\Delta
t=0)=\dfrac{\tau\,\sigma_{\mathcal{T}}^{2}}{2x_{p}^{2}}\,,$ (125)

that coincides with the asymptotic value of the process variance. In other
words, we have explicitly checked that the power extracted from the Fourier
transform of the Langevin equation coincides with the one obtained from the
Wiener-Khinchin theorem in the stationary limit. Unfortunately, this does not
hold for the integrated process (namely, for the phase), as discussed below.

###  B.2 Phase: integrated Ornstein-Uhlenbeck process

Since the phase residual $\delta\phi_{t}$ is just the integrated Ornstein-
Uhlenbeck process

$\delta\phi_{t}\,=\,\int_{0}^{t}dz\,\Omega_{p}(z)\,,$ (126)

its associated PSD $P_{\phi}(\omega)$ is expected to be just (124) divided by
$\omega^{2}$, see equation (88). We check to what extent this claim101010 In
principle, a more rigorous introduction of the PDS requires the windowing of
the signal, so that the Fourier transform of the windowed signal converges
(Kasdin, 1995).  about $P_{\phi}(\omega)$ allows us to estimate the average
power of the the signal $\phi_{t}$, that is directly related to the timing
noise strength $\sigma^{2}$ introduced in (87). Applying the stochastic Fubini
theorem directly to (122), it is possible to show that the average and the
covariance are

$\begin{split}\langle\delta\phi^{t}\rangle&=\int_{0}^{t}\\!\\!dz\,\langle\delta\Omega_{p}(z)\rangle=\tau\,\delta\Omega_{p}(0)\,\left(1-e^{-t/\tau}\right)\\\
C(t,u)&=\langle[\delta\phi_{t}-\langle\delta\phi_{t}\rangle][\delta\phi_{u}-\langle\delta\phi_{u}\rangle]\rangle=\dfrac{\tau^{2}\,\sigma_{\mathcal{T}}^{2}}{x_{p}^{2}}\,\min(t,u)\,+\\\
&+\dfrac{\tau^{3}\,\sigma_{\mathcal{T}}^{2}}{2x_{p}^{2}}\left[2e^{-t/\tau}+2e^{-u/\tau}-2-e^{-|t-u|/\tau}-\tau
e^{-(t+u)/\tau}\right]\end{split}$ (127)

Now, the above expression of the covariance $C$ tells us that $\delta\phi_{t}$
has no stationary limit. For example, the variance depends on $t$ even in the
limit $t\gg\tau$,

$C(t,t)\approx\dfrac{\tau^{2}\,\sigma_{\mathcal{T}}^{2}}{x_{p}^{2}}\,\left(t-6\tau\right)+O\left(e^{-t/\tau}\right).$
(128)

However, for the sake of extracting the timing noise strength $\sigma^{2}$, we
can take the average power on a long time interval of duration $T_{o}$ see
equation (87): we write $\sigma^{2}$ in terms of the process’s variance

$\sigma^{2}=\frac{1}{T_{o}\,\Omega^{2}}\int_{0}^{T_{o}}\\!\\!dt\,C(t,t)\approx\dfrac{T_{0}\,\tau^{2}\,\sigma_{\mathcal{T}}^{2}}{2\,x_{p}^{2}\Omega^{2}}-\dfrac{3\,\tau^{2}\,\sigma_{\mathcal{T}}^{2}\,\dot{\Omega}^{2}\,x_{p}}{2\,x_{p}^{2}\Omega^{2}}\,,$
(129)

where we have taken the limit $T_{0}\gg\tau$. This limit is equivalent to
$\nu_{0}\ll\mathcal{B}\Omega$, where $\nu_{0}\sim 1/T_{0}$ is the infrared
cut-off in (89). If we use the noise parametrization in (57), we finally
obtain

$\sigma^{2}\approx\dfrac{T_{0}\,x_{1}^{2}\,\alpha_{\mathcal{T}}^{2}\,\dot{\Omega}^{2}}{16\,\mathcal{B}^{3}\Omega^{5}}\,\left(\,1-\dfrac{3\,x_{p}}{2\,T_{0}\,\mathcal{B}\,\Omega}\,\right)$
(130)

Therefore, apart from constant (dimensionless) numerical factors, the timing
noise strength associated to the variance of the timing residual is expected
to scale as, cf. with (93),

$\sigma^{2}\propto
T_{0}\,\dot{\Omega}^{2}\,\Omega^{-5}\,\propto\,\nu_{0}\,\dot{\Omega}^{2}\,\Omega^{-5}\,\qquad\quad(\text{for}\quad\nu_{0}\ll\mathcal{B}\Omega)$
(131)

This is consistent with the scaling in (96), see also Table 1: we have
explicitly checked that the band-limited integration of the PSD in equation
(89) leads to an estimation of the timing noise strength that has the same
scaling properties of the one estimated directly from the averaged variance of
the signal. This is obvious when the PSD is integrable over the whole real
line, namely when no infrared cutoff is needed $\nu_{0}=0$, as in equation
(124).