Analysing radio pulsar timing noise with a Kalman filter: a demonstration involving PSR J1359$-$6038

In the standard two-component crust-superfluid model of a neutron star, timing noise can arise when the two components are perturbed by stochastic torques. Here it is demonstrated how to analyse fluctuations in radio pulse times of arrival with a Kalman filter to measure physical properties of the two-component model, including the crust-superfluid coupling time-scale and the variances of the crust and superfluid torques. The analysis technique, validated previously on synthetic data, is applied to observations with the Molonglo Observatory Synthesis Telescope of the representative pulsar PSR J1359$-$6038. It is shown that the two-component model is preferred to a one-component model, with log Bayes factor $6.81 \pm 0.02$. The coupling time-scale and the torque variances on the crust and superfluid are measured with $90\%$ confidence to be $10^{7.1^{+0.8}_{-0.5}}$ $\rm{s}$ and $10^{-24.0^{+0.4}_{-5.6}}$ $\rm{rad^2~s^{-3}}$ and $10^{-21.7^{+3.5}_{-0.9}}$ $\rm{rad^2~s^{-3}}$ respectively.


INTRODUCTION
High-precision pulsar timing data are a rich and plentiful resource for probing the properties of neutron star interiors.They show that pulsars decelerate secularly with small, random fluctuations away from the secular trend in the form of occasional glitches (Shemar & Lyne 1996;Hobbs et al. 2010;Espinoza et al. 2011;Haskell & Melatos 2015) and persistent pulsar timing noise (Shannon & Cordes 2010;Lentati et al. 2016;Parthasarathy et al. 2019;Lower et al. 2020).The random fluctuations can be related to internal physics, such as farfrom-equilibrium processes involving superfluid vortices (Anderson & Itoh 1975;Warszawski & Melatos 2011), type II superconductor flux tubes (Ruderman et al. 1998;Drummond & Melatos 2017), and hydrodynamic turbulence (Greenstein 1970;Melatos & Link 2014).One prominent structural theory supported by these observations is the two-component model, in which neutron stars are composed of a rigid outer crust and a superfluid core (Baym et al. 1969).The crust and superfluid are coupled by forces, which act to restore corotation and therefore damp observable fluctuations in the crust's angular velocity, as happens during post-glitch recoveries.With respect to timing noise specifically, the tendency to restore corotation through crust-superfluid coupling should be detectable statistically, e.g. through the auto-correlation time-scale (Price et al. 2012).One advantage is that timing noise occurs all the time, so there is an abundance of data to work with.
The two-component model of neutron stars was originally motivated by observations of pulsar glitches (Baym et al. 1969).Glitches are rare events in which a pulsar spins up impulsively.Following the event, the spin-down rate (which sometimes also changes during the ★ E-mail: noneill@student.unimelb.edu.auevent) returns partially or wholly to the rate before the event, due to some friction-like interaction between the crust and superfluid operating over some time-scale  (Anderson & Itoh 1975;Lamb et al. 1978a;Pines & Alpar 1985;Haskell & Melatos 2015;Antonelli et al. 2022;Zhou et al. 2022;Antonopoulou et al. 2022).As timing noise also involves a perturbation of the crust's angular velocity, it is likely to either drive or be driven by glitch-like differential rotation between the crust and superfluid, which is damped by the same restoring forces that bring the two components back into corotation after a glitch (Lamb et al. 1978b;van Eysden & Melatos 2010;Price et al. 2012;Gügercinoǧlu & Alpar 2017).We aim to measure the coupling between the two components by looking for statistical evidence of its associated relaxation time-scale  in timing noise data.
Many previous timing noise studies focus on calculating ensemble averaged quantities derived from the whole data set such as the noise amplitude (Boynton et al. 1972;Groth 1975;Cordes 1980;Cordes & Helfand 1980;Shannon & Cordes 2010) and power spectrum (van Haasteren et al. 2009;Hobbs et al. 2010;Lentati et al. 2016;Goncharov et al. 2020;Parthasarathy et al. 2019Parthasarathy et al. , 2020)).Although ensemble averages are good for measuring some quantities, they do not make use of the information contained in the specific, timeordered sequence of times of arrival, that is, the specific realization of the random process presented to the observer (Vargas & Melatos 2023).In this paper, we seek to extract this extra information using a Kalman filter, a tool commonly used for signal processing in electrical engineering applications (Kalman 1960).Given a model for how a system evolves with time (here the two-component model), the Kalman filter tracks the specific random fluctuations in the observed time series and estimates the most likely state of the system at every time step, i.e. the most likely values of the dynamical variables (the crust and superfluid angular velocities) in the two-component model.
Moreover, when combined with a nested sampler, the Kalman filter estimates the most likely values of the static parameters in the model, e.g. the crust-superfluid coupling coefficient and the variances of the stochastic torques.Two previous papers (Meyers et al. 2021a,b) explained this approach and validated it with synthetic data.
In this paper we demonstrate in principle that the approach works successfully on real data by applying it to Molonglo Observatory Synthesis Telescope observations of PSR J1359−6038 (Bailes et al. 2017;Lower et al. 2020).This object is chosen, because its timing noise power spectrum is similar to that predicted by the simplest version of the two-component model, as discussed below (see Section 4 and Appendix A).It also offers relatively good quality data, with 429 times of arrival (TOAs) over 3.5 years of monitoring, with an average TOA uncertainty of 5.1 × 10 −5 s.Other pulsars can be analysed but an analysis of a sample of objects is outside the scope of this paper, whose goal is to demonstrate the real-world applicability of the method by way of a worked example, to pave the way for fuller studies in the future.
The paper is structured as follows.In Section 2, the two-component model is introduced.In Section 3, the Kalman filter and its implementation are explained, including the question of parameter identifiability and the choice of priors based on existing observations.In Section 4 we discuss the properties of the data from PSR J1359−6038, including the number of observations and measurement uncertainties, and describe the conversion from TOAs to local frequencies.In Section 5 we present the parameter estimates for PSR J1359−6038.In Section 6 we present a Bayesian comparison between the one-component and two-component models.In Section 7, we conclude by interpreting the results and validating them with Monte Carlo posterior checks to verify their statistical significance.The appendices contain for completeness an analytic derivation of the power spectral density of the two-component model, the explicit formulas for the matrices used in the Kalman filter, the Kalman filter recurrence relations, calculations demonstrating the accuracy of the parameter estimates as a function of data volume and measurement errors, and details of the Bayesian comparison of the two-component model with a simpler, one-component model.

TWO-COMPONENT MODEL
The two-component model of a neutron star consists of a rigid crust and superfluid core, which are assumed to rotate uniformly, with angular velocities Ω c and Ω s respectively.The components obey the equations of motion (Baym et al. 1969;Gügercinoǧlu & Alpar 2017) where the subscripts "c" and "s" label the crust and the superfluid respectively,  c and  s are the moments of inertia,  c and  s are constant external torques,  c and  s are stochastic torques, and  c and  s are coupling time-scales.The astrophysical origins of the torques on the right-hand sides of (1) and ( 2) are discussed below.
The stochastic torques obey Gaussian, white noise statistics with where ⟨. ..⟩ denotes the ensemble average, and  c and  s are noise amplitudes.
The coupling between the crust and superfluid is assumed to act like friction between the two components, reducing their relative velocity.In this paper |Ω c − Ω s | is assumed to be small enough that the coupling torque is approximately linear in this difference on the right hand sides of equations ( 1) and ( 2).Over a typical observing time-scale (decades), and over the relaxation time-scales  c and  s (typically weeks) (Price et al. 2012), the torques  c and  s can be approximated as constants. c represents the net secular external torque applied to the crust, including the electromagnetic radiation reaction torque (Goldreich & Julian 1969) and the gravitational radiation reaction torque (Ostriker & Gunn 1969;Ferrari & Ruffini 1969). s represents the net secular torque on the superfluid which includes electromagnetic and gravitational components, if the superfluid is threaded by a corotating magnetic field and possesses a time-varying mass or current quadrupole moment respectively.The torque  c represents the net random torque acting on the crust, which may be negligible in a nonaccreting radio pulsar, except perhaps in the presence of seismic activity (Middleditch et al. 2006;Chugunov & Horowitz 2010;Kerin & Melatos 2022;Giliberti & Cambiotti 2022). s represents the net random torque acting on the superfluid, including the torque from superfluid vortex unpinning (Anderson & Itoh 1975;Warszawski & Melatos 2011).The stochastic torques drive timing noise in the observable Ω c , even in the scenario  c = 0 and  s ≠ 0.
Equations ( 1) and ( 2) have been generalized recently to multiple internal components by Antonelli et al. (2023).We restrict the analysis in this paper to two components, because the available data are insufficient to constrain a more complicated model.Indeed, the data are insufficient even to constrain (1) and (2) fully, as discussed in Section 3.

KALMAN TRACKING AND ESTIMATION
Given an observed time series Ω c ( 1 ), . . ., Ω c (  ), one can use a Kalman filter to infer the most probable sequence of states Ω c (  ) and Ω s (  ) (1 ≤  ≤ ) traversed by the system described by (1) and (2).In this paper, the Kalman filter assumes the measurements are the true state of the system plus some Gaussian measurement noise.The Kalman filter then uses the assumption that the system evolves according to (1) and (2) to separate out the real evolution of the system from the measurement noise on a most probable basis, which gives estimates of Ω c (  ) and Ω s (  ) at each   together with an error on each estimate.This procedure, termed Kalman tracking, is described in Section 3.1.The Kalman filter's ability to recover the model parameters from Ω c data is assessed in Section 3.2.To estimate the parameters, one calculates the Kalman filter likelihood (of the model producing the observed data) as a function of the static parameters to find their most probable values given Ω c (  ).This procedure, termed Kalman estimation, is described in Section 3.3.The results depend on the priors, which are specified in Section 3.4 together with their astrophysical motivations.

Kalman tracking
The state of the system at time   is denoted by The Kalman filter predicts the next state of the system given the current one.Solving (1) and (2) gives a discrete equation for updating the state from time Equation ( 6) multiples  −1 by a transition matrix   , which comes from the coupling terms in (1) and (2); adds a vector   , which comes from the deterministic torques; and adds a random vector   , which comes from the stochastic torques and is drawn from a Gaussian distribution with Explicit expressions for   ,   and   are given in Appendix B.
For general physical systems, the state and measurements at   may be related in a complicated non-linear way.In this paper, however, the relation is simple: the measured and true Ω c (  ) differ by an additive measurement error, and Ω s (  ) is hidden, i.e. it cannot be measured directly at all.Mathematically, we encode this by saying that a measurement   at time   is related to the state of the system   by In ( 9),   is the measurement noise and is drawn from a normal distribution with is the observation matrix which determines which components of the state can be measured.In radio pulsar timing experiments, only the crust can be measured, since the core is hidden from view, implying1 The Kalman filter computes the expected evolution of the system through (6).It updates the expected evolution with new information from the measurement at   through (9), separating the process noise   from the measurement noise   .The detailed implementation of this two-step, predictor-corrector algorithm is discussed in Appendix C to help the interested reader reproduce the results in Sections 5 and 6.

Identifiability
In general, there is no guarantee that all six of the static parameters  c ,  s ,  c ,  s ,  c ,  s in (1) and (2) can be identified from the time series Ω c ( 1 ), . . ., Ω c (  ), even if the data volume is arbitrarily large ( → ∞).Certain parameters cannot be separated from the rest and can be estimated only in combination.This issue, known as identifiability, is widely recognized in electrical engineering applications.A number of formal techniques have been developed to handle it, as summarized by Bellman & Åström (1970) for example.In this section, we analyze what combinations of the static parameters in (1) and (2) are identifiable, in order to interpret the posterior distribution.
As a starting point, it is instructive to express (1) and (2) as an equation for Ω c on its own, viz.
Equation ( 13) determines Ω c () fully, supplemented by initial conditions Ω c ( 1 ) and Ω c ( 1 ).Therefore the static parameter combinations that appear in (13) are identifiable -that is, they can be estimated from the data -but they cannot be decomposed into their elements.For example,  c and  s cannot be estimated separately, because they appear in the irreducible combination  −1 c +  −1 s in the first term on the right-hand side of (13), and no additional, independent equation of motion exists beyond (13), in which  c and  s appear separately.
The first term on the right-hand side of (13) contains the parameter which is a combined relaxation time-scale for the system as a whole.
The second term (when multiplied by ) contains the parameter which is the ensemble-averaged frequency derivative of the pulsar.
The noiseless evolution of Ω c is governed solely by  and ⟨ Ω c ⟩ so these parameters are easiest to estimate.For notational convenience, and to assist with physical interpretation, we also introduce the complementary parameter combinations, In ( 16) and ( 17),  is the ratio of relaxation time-scales, which equals  s / c when the coupling torques form an action-reaction pair, and ⟨Ω c − Ω s ⟩ is the ensemble-averaged angular velocity lag between the crust and superfluid.
It is less straightforward to read off by sight whether the noise amplitudes  c =  2 c / 2 c and  s =  2 s / 2 s can be estimated by the Kalman filter.However, in Appendix D, the question is answered empirically using synthetic data.It turns out that the Kalman filter can estimate  c reliably (in addition to  and ⟨ Ω c ⟩) but usually not  s .An analytic treatment of the identifiability of  c and  s is given in Appendix E. The analysis suggests that it is easiest to recover the parameter combinations , ,  c ,  s , ⟨Ω c − Ω s ⟩, ⟨ Ω c ⟩, and we switch to them in what follows.

Kalman filter likelihood
We construct a likelihood function ( 1: |) with  1: = { 1 ,  2 , ...,   } for the data given a choice of parameters.We calculate ( 1: |) from the optimal state sequence output by the Kalman filter.Intuitively, if the parameters are chosen well (i.e.near their true values), the model's trajectory in time matches the data closely, and ( 1: |) is relatively large.
Let   be the difference between the measurement and the prediction of the state at   , and let   be the covariance of   .Then Bayes's rule gives the posterior distribution of  as ln ( where   is the dimension of   , which in our case is always unity, () is the prior distribution of , and  is the evidence, which acts as a normalisation constant.Appendix C supplies more details on the derivation of ( 18).The dynesty sampler (Speagle 2020) is used to efficiently sample the posterior distribution.

Prior distribution
In the absence of other information, we choose log-uniform priors for the parameters that are positive definite and plausibly span several decades, namely , ,  c and  s .We choose uniform priors for the parameters with small ranges, namely ⟨ Ω c ⟩ and ⟨Ω c − Ω s ⟩.As more pulsars are analysed in the future, and a picture of the distribution of (say)  across the population emerges, it may become appropriate to use more informative priors in future work.We can use previous measurements of glitch relaxation time-scales to infer reasonable bounds on , motivated by ( 14).The data in Yu et al. (2013) imply 10 5 ≤ /(1 s) ≤ 10 8 .Models of stellar structure suggest 10 −2 ≤  =  s / c =  s / c ≤ 10 2 (Link et al. 1999;Lyne et al. 2000;Espinoza et al. 2011;Chamel 2012).Reasonable bounds on  c and  s follow from previous, independent measurements of timing noise amplitude.For the crust one measures typically 10 −30 ≤  c /(1 rad 2 s −3 ) ≤ 10 −16 (Cordes & Downs 1985;Çerri-Serim et al. 2019;Lower et al. 2020;Vargas & Melatos 2023).In the absence of additional information, we assume the same range for  s , noting that radio pulsar timing experiments cannot track the core to measure  s directly.For ⟨ Ω c ⟩, a central estimate ⟨ Ω c ⟩ cent and an uncertainty  can be obtained by fitting a linear trend to the time series Ω c (  ).We adopt a uniform prior for where the subscripts "max" and "min" denote the maximum and minimum values of the prior range given earlier in this section.For most pulsars, one finds −10 −11 ≤ ⟨ Ω c ⟩ ≤ 0 (Helfand et al. 1980) and hence −10 −3 ≤ ⟨Ω c − Ω s ⟩ ≤ 10 −3 .For more details on justifying the prior distribution, the reader is referred to Meyers et al. (2021a,b).

DATA
The data used in this paper come from the UTMOST pulsar observing program 2 (Lower et al. 2020).They are in the form of barycentered pulse times of arrival (TOAs).To assist parameter estimation, we select a test object with a relatively large number of TOAs with relatively small error bars.We also avoid objects exhibiting glitches, because glitches are not included in the dynamical model ( 1) and ( 2).An object which satisfies these criteria is PSR J1359−6038.The UTMOST data release for PSR J1359−6038 contains 429 TOAs, one of which we discard as an outlier 3 .The average TOA error is 51 s.The TOAs span 1263 days.In Lower et al. (2020) the timing noise in PSR J1359−6038 is measured to have a spectral index of  = −5.1 +0.8 −0.9 , i.e. the power spectral density of the phase residuals scales ∝   at  ≳ 10 −8 Hz.This is comparable to the theoretical prediction for the two-component model ( 1) and (2), of  = −4, given in (A13).
The Kalman filter in Section 3 ingests pulse frequency data.Hence the phase information in the TOAs must be converted to local frequencies Ω c ( 1 ), . . ., Ω c (  ).This is done using the standard pulsar timing software package, tempo2 (Hobbs et al. 2006;Edwards et al. 2006).To generate a local frequency estimate, we feed a set,   , of consecutive TOAs into tempo2 and extract Ω(  ) and its associated error, where   is the average of the TOAs in the set   .We construct the disjoint sets,  1 ,  2 , ...,   , from the TOAs,  1 ,  2 , ...,   ( ≥ ), starting with  1 , according to the following two rules: (i)   contains   and all   >   with |  −   | < 10 days; and (ii) there must be at least three TOAs per set, and the 10-day window is lengthened when necessary to achieve this.
The fitting process yields 62 frequency data points spanning 1220 days (the   values span a slightly shorter time than the   values).We find Ω c (  ) = Ω c ( 1 ) + (  −  1 ) Ω c ( 1 ) to an excellent approximation, with Ω c ( 1 ) = 49.2767rad s −1 and Ω c ( 1 ) = −2.4472× 10 −12 rad s −2 .Subtracting the linear, long-term trend yields the frequency residuals plotted in Fig. 1.The standard deviation of the residuals is 1.7 × 10 −8 rad s −1 .The mean uncertainty is 7.8 × 10 −9 rad s −1 .However the size of the error bars varies with the errors of the original TOAs, the number of TOAs used in a fit, and the time span of a fit.The interval from MJD 57400 to MJD 57900 features TOAs with relatively large error bars and spacings; the median uncertainty on Ω(  ) therein is 1.3 × 10 −8 rad s −1 compared to 2.6 × 10 −9 rad s −1 for the measurements outside that interval.Some fitted frequencies deviate significantly from the linear trend, possibly due to an issue with fitting to TOAs with large errors.We remove two outliers with residuals from the trend below −10 −7 rad s −1 .
The conversion from TOAs to frequencies is imperfect.Fitting frequencies filters the timing noise on the time-scale over which the fit is done, so some information is lost.The effect of filtering is tested on simulated TOAs in Appendix F. We find that it causes  c to be slightly underestimated and makes  harder to recover.

ESTIMATED PARAMETERS OF THE TWO-COMPONENT MODEL
In this section, we apply the Kalman tracking and estimation procedure described in Section 3 to the UTMOST data described in Section 4. Marginalized posterior distributions for the static parameters that can be identified reliably, namely ,  c ,  s and ⟨ Ω c ⟩, are presented in Section 5.1.The four parameters are discussed individually in greater detail in Sections 5.2-5.4 and interpreted astrophysically.The two parameters that cannot be identified reliably, namely  and ⟨Ω c −Ω s ⟩, are discussed in Section 5.5.

Joint posterior distribution
The six-dimensional joint posterior distribution calculated from the PSR J1359−6038 frequency data in Fig. 1 using equation ( 18) is displayed in the traditional format of a corner plot in Appendix G. Four parameters can be estimated reliably and exhibit well-formed peaks: ,  c ,  s and ⟨ Ω c ⟩. Two parameters,  and ⟨Ω c − Ω s ⟩, cannot be estimated reliably and rail against the prior bounds.One-dimensional cross-sections of the posterior, marginalized over all but one parameter, are displayed as histograms in Fig. 2 for the four parameters that can be estimated reliably.The marginalized posteriors in Fig. 2 are unimodal, with the peaks falling comfortably within the respective prior ranges in Table 1.There are distinct peaks for ,  c ,  s and ⟨ Ω c ⟩.They occur at  = 1.3 × 10 7 s,  c = 1.05 × 10 −24 rad 2 s −3 ,  s = 1.82 × 10 −22 rad 2 s −3 and ⟨ Ω c ⟩ = −2.4475× 10 −12 rad s −2 .
The widths of the 90% confidence intervals of the histograms in Fig. 2 are 1.3 dex, 6.0 dex, 4.4 dex and 2.4 × 10 −15 rad s −2 for ,  c ,  s , ⟨ Ω c ⟩ respectively.The 90% confidence intervals cover 43%, 43%, 31%, and 1.4% of their respective prior ranges, confirming that the four parameters are identified reliably.Numerical values for the parameter estimates for PSR J1359−6038 conditional on the two-component model given by ( 1) and ( 2) are summarized for completeness in Table 2.

Coupling time-scale
In Fig. 2 we find that the marginalised posterior for  has a well defined peak.With 90% confidence we recover 3.8 × 10 6 ≤ /(1 s) ≤ 7.6 × 10 7 .This result is broadly consistent with previous measurements of timing noise and post-glitch relaxation time-scales.Price et al. ( 2012) computed the auto-correlation time-scale of timing residuals in the objects PSR B1133+16 and PSR B1933+16 and obtained ≈ 10 days and ≈ 20 days respectively.Post-glitch relaxation time-scales have been measured previously by many authors (Mc-Culloch et al. 1990;Alpar et al. 1993;Shemar & Lyne 1996;Lyne et al. 2000;Wang et al. 2000;Wong et al. 2001;Dodson et al. 2002;Yuan et al. 2010;van Eysden & Melatos 2010;Espinoza et al. 2011;Yu et al. 2013;Espinoza et al. 2021;Lower et al. 2021;Gügercinoğlu et al. 2022).They mostly range from ≈ 10 days to ≈ 3 × 10 2 days and occasionally reach as high as ≈ 10 3 days.The 90% confidence interval in Fig. 2 overlaps this range.The prior on  in Table 1 is motivated by glitch observations, so the overlap is expected, but it is significant that the posterior peaks comfortably within the prior range.
The detailed shapes of the peaks in  c and  s and to a lesser extent their positions are influenced by the tempo2 fitting process which constructs local frequencies from TOAs, as described in Section 4. To calibrate for this effect, we conduct Monte Carlo simulations with synthetic data comparing the  c and  s estimates inferred from local tempo2 TOA fits versus direct frequency data.The results are presented in Appendix F. When calculating local frequencies from TOAs, tempo2 fits a linear frequency model ∝ Ω c (  −  −1 ), which does not make any allowance for a random walk in the interval  −1 ≤  ≤   .Consequently tempo2 acts as a low-pass filter, which smooths out features on inter-TOA time-scales.If TOAs are spaced widely enough, so that one has  ≪   −  −1 , then tempo2 smooths out some evidence for  as well.Smoothing may also bias  c and  s (Deeter 1984;Meyers et al. 2021a).Because (1) and ( 2) are linear differential equations, the timing noise in Ω c can be written as a linear combination of the timing noise from  c and  s .These two contributions are separated and plotted in Fig. 3 for simulated data.The contribution from  c is rougher than from  s .This is because  c appears directly in equation (1) for Ω c /, whereas  s only affects Ω c indirectly through its time integral in Ω s in the coupling term.Hence smoothing the data nullifies  c fluctuations more than  s fluctuations, overestimating  s and underestimating  c .

Unidentifable parameters
The marginalized posteriors for  and ⟨Ω c − Ω s ⟩ are not sharply peaked, as is apparent from columns 1 and 5 in the corner plot in Fig. G1 in Appendix G.There is no evidence of railing against the prior bounds, but the marginalized posterior is flat and therefore uninformative; the probability at the extremes of the prior range is ≈ 80% and ≈ 90% of the peak probability for  and ⟨Ω c − Ω s ⟩ respectively.Table 2. Static parameters inferred by the Kalman tracker and nested sampler for the two-component model, extracted from the six-dimensional joint posterior in Fig. G1.The ,  c ,  s and ⟨ Ω c ⟩ rows correspond to the four identifiable parameters, whose marginalized posteriors are plotted in Fig. 2. The  and ⟨Ω c − Ω s ⟩ rows correspond to the unidentifiable parameters discussed in Section 5.5.The fourth and fifth columns list two complementary measures of the dispersion in the posterior, viz. the full-width half-maximum (FWHM) and 90% confidence intervals respectively.

Parameter
The difficulty in recovering  and ⟨Ω c − Ω s ⟩ is predicted by the identifiability analysis in Section 3.2.The deterministic part of equation (13) does not feature  and ⟨Ω c − Ω s ⟩.Relatedly, the simulations in Appendix D show the recovery of the six model parameters on simulated Ω c data and confirm that  and ⟨Ω c − Ω s ⟩ are poorly recovered.The average distances of the recovered log 10  and ⟨Ω c − Ω s ⟩ parameters in Fig. D1 from their true values as a fraction of their prior widths are 41% and 32% respectively.

ONE-COMPONENT MODEL
One may ask whether the two-component model described by Equations (1) and ( 2 challenges of identification arise because the model is more complicated than it needs to be. To test the above hypothesis, we repeat the Kalman filter analysis in Sections 3-5 for a one-component model.The single component corresponds to the crust (subscript 'c'), which is phase-locked to the radio pulsations.The one-component equation of motion takes the form  c Ω c / =  c +  c (), analogous to Equation (1) but without the crust-superfluid coupling, where  c () is a Langevin driving torque; see Appendix H for details.Parameter estimates for the model parameters, namely  c and ⟨ Ω c ⟩, are presented in the format of a traditional corner plot in Appendix H; see Fig. H1 and Table H1.The posterior peaks at log 10  c /(1 rad 2 s −3 ) = −22.73+0.23  −0.19 and ⟨ Ω c ⟩ = −2.4473× 10 −12 ± 7 × 10 −16 rad s −2 (90% confidence intervals).
The ⟨ Ω c ⟩ value recovered above is similar to that recovered for the two-component model, because it is insensitive to how the timing noise is modelled.However, the recovered  c for the one-component model is larger than for the two-component model by a factor of 10 1.3 ≈ 20.This is likely because  c and  s combine through the crust-superfluid coupling to generate the noise in Ω c () in the two-component model, whereas only  c is responsible in the onecomponent model, and we find  s ∼ 10 2  c in Fig. 2 and Table 2 for the two-component model.
We can compare  c for the one-component model to the PSD normalisation  measured by other researchers as in Section 5.3.The one-component estimate of  c in Table H1 converts to log 10  = −9.87+0.12  −0.09 .This is consistent with the estimate for PSR J1359−6038 in Lower et al. (2020) and with estimates for ordinary and millisecond pulsars more broadly (Parthasarathy et al. 2019;Lower et al. 2020;Keith & Niţu 2023).
A Bayesian model comparison shows that the two-component model is strongly preferred in a Bayesian sense, with a log 10 Bayes factor of 6.81 ± 0.02 relative to the one-component model.Further details can be found in Appendix H.

CONCLUSION
Traditionally, timing noise studies proceed by comparing the measured variance in the phase residuals to that predicted by a microphysical or phenomenological model (Boynton et al. 1972;Groth 1975;Cordes 1980;Cordes & Helfand 1980;Shannon & Cordes 2010;Melatos & Link 2014), fitting a microphysical or phenomenological model to the power spectral density of the phase residuals (van Haasteren et al. 2009;Hobbs et al. 2010;Lentati et al. 2016;Parthasarathy et al. 2019Parthasarathy et al. , 2020;;Goncharov et al. 2020), or measuring a relaxation time-scale using the autocorrelation function of the phase residuals (Price et al. 2012).These approaches raise interesting questions about the type of variance measured, e.g.Allan variance (Matsakis et al. 1997;Hobbs et al. 2010;Shannon & Cordes 2010;Melatos & Link 2014), or biases in constructing the power spectral density (Coles et al. 2011;van Haasteren & Levin 2013;Keith & Niţu 2023).In this paper we fit a model directly to the time-ordered Ω c () data without averaging implicitly over an ensemble, i.e. without analyzing the phase residual PSD.The extra information made available by the analysis of a unique, time-ordered, random realization makes it feasible to estimate reliably four of the six static parameters in the classic, two-component, crust-superfluid model of a neutron star, and to distinguish statistically between one-and two-component models, with interesting astrophysical implications.
Bayesian model selection shows that the classic two-component model in Section 2 fits the data better than the representative onecomponent model given by equations (H1)-(H3), with a log 10 Bayes factor of 6.81 ± 0.02.The Kalman filter's ability to recover successfully ,  c , and  s for the two-component model, noting that  and  s do not feature in the one-component model, adds support to the conclusion of the model selection exercise.
The results in this paper exemplify the usefulness in astrophysics of parameter estimation methods based on Kalman filtering and similar algorithms (Meyers et al. 2021a).In the future, more data and more sophisticated models will give more accurate parameter estimates.The Kalman filter can be rewritten easily to track using physical models other than equations ( 1) and ( 2).A nonlinear torque can be incorporated to calculate the braking index, if the data volume is sufficient (Vargas & Melatos 2023).The model currently assumes that  c and  s are white noise torques, which implies that the PSD of Ω c fluctuations scales as  −4 at large  , a property which is satisfied approximately by some but not all pulsars (Lower et al. 2020;Antonelli et al. 2023;Keith & Niţu 2023).Generalizing  c and  s to colored noise is a standard procedure in electrical engineering (Gelb 1974).Bayesian model selection between these physically motivated alternatives is straightforward too, because the Kalman tracker and nested sampler in this paper together generate the Bayesian evidence as a by-product of the analysis.
The next step in this investigation is to extend the Kalman tracker so that it operates on TOAs directly instead of converting them first to local frequencies Ω c (  ) using tempo2.Monte Carlo simulations in Appendix F show that local tempo2 computation of Ω c (  ) biases the static parameter estimation results, underestimating  c and making  harder to infer.This is not a fault with tempo2; it is a general consequence of the low-pass filtering introduced by any local frequency fitting process.Once the Kalman tracker is extended to operate on TOAs directly, it will be appropriate to apply the method to a wider selection of pulsars beyond PSR J1359−6038 and conduct population studies of the recoverable quantities ,  c , and  s , which are important physically and are measured in only a few objects to date.We note in closing that the Kalman tracker and nested sampler in this paper are easy to implement and quick to run.By way of calibration, the PSR J1359−6038 analysis in this paper takes of order one hour to run on ∼ 10 3 TOAs.More generally, the run time scales in direct proportion to the number of TOAs.

APPENDIX C: KALMAN FILTER
In this appendix, for the sake of completeness and reproducibility, we summarize the structure and operation of the Kalman filter used to generate the results in Sections 5 and 6.The discussion includes a short justification of the form of the Kalman likelihood in equation ( 18).The Kalman filter is a predictor-corrector algorithm.Given a sequence of  − 1 measurements,  1 ,  2 , ...,  −1 , the Kalman filter predicts the value of the state   at the next time step   .The estimate is denoted by  −1  , and its error is given by the covariance matrix  −1  .When the measurement   is made, the estimate  −1  is updated to give the new estimate of the current state,    , and its error,    .This procedure is carried out from  1 to   .The following update formulas (C1)-(C9) implement the prediction and correction steps described above.The symbols are explained in the text immediately following. Initialisation: State prediction: State correction: Equations ( C1) and (C2) estimate the initial state and its error.Equations (C3) and (C4) predict the next state from the previous data.Equations ( C5) through (C9) combine the new measurement with the prediction to get a new estimate of the current state.In (C5), we define   =   −  (  | 1:−1 ).In (C6), we define   = var(  ).The matrix   in (C7)-( C9) is usually called the Kalman gain.
To get the likelihood, ( =  1: |), from the Kalman filter estimates, we apply standard rules for conditional probability to get and hence recursively is Gaussian if all the errors are assumed to be Gaussian.The mean and covariance matrix of the ( 1: |) distribution are which imply where N ( ; , ) denotes a Gaussian with mean  and covariance matrix .Equivalently, writing   =   −   −1  , we obtain The full formula for the log-likelihood then becomes log ( where   is the dimension of .Equation (C18) follows from (C12) and (C19) follows from (C17).Equation (C19) is the same as equation ( 18) in Section 3.3.Once ( 1: |) is known, Bayes's theorem can be used to get the probability of the parameters in terms of the data,

APPENDIX D: PARAMETER ESTIMATION WITH SIMULATED FREQUENCIES D1 Simulation under ideal conditions
In this appendix our parameter estimation method is tested on simulated frequency data.To simulate the data, we must first select values for the model parameters  c ,  s ,  c / c ,  s / s ,  c / c ,  s / s .We then chose an initial value Ω c (0) and set Ω s (0) to be so the pulsar is initially in equilibrium.We choose  obs random times from  = 0 to  =  obs , where simulated measurements occur.
The differential equations ( 1) and ( 2) are then integrated from the initial time to each measurement time, giving Ω c (  ) and Ω s (  ) values at each 1 ≤  ≤  obs .Realistic timing experiments only yield Ω c observations, so the Ω s values are discarded.To simulate measurement errors we add a number drawn from a Gaussian with variance  to each Ω c value.For simplicity in these simulations we assume every data point has the same measurement uncertainty, unlike in real data.Once the data are simulated, the parameter estimation algorithm is executed, and the posterior distribution is compared to the true parameters.Fig. D1 shows the results for 200 simulations with randomly chosen values for , ,  c ,  s , ⟨Ω c − Ω s ⟩, ⟨ Ω c ⟩ and with Ω c (0) = 10 rad s −1 .The simulations are under idealised conditions (i.e.low noise, many samples) with  obs = 1000 days,  obs = 1000 and  = 10 −30 rad 2 s −2 .The six panels plot the recovered values for each of the six parameters against the injected parameters.The diagonal red lines mark where the recovered and injected parameters are equal, indicating a successful simulation.The colours indicate the width of the marginalised posterior for that parameter (width here meaning the inter-quartile range).The closeness of the points to the red line in the ,  c and ⟨ Ω c ⟩ panels means those parameters are generally well recovered, while the large vertical spread of the points in the ,  s and ⟨Ω c − Ω s ⟩ panels means that they are usually harder to recover.This agrees with the identifiability analyses in Section 3.2 and Appendix E. Interestingly, small  values tend to be better recovered, which is unsurprising since these correspond to strong damping.
As a further verification of the method, in Fig. D2 we show a standard PP-plot (probability-probability plot) for the same 200 simulations as in Fig. D1.The curves for each parameter remain close to the diagonal line, indicating that the posteriors give unbiased estimates.For more details on interpretation of PP-plots see Meyers et al. (2021b).

D2 Measurement noise
In this appendix we test the Kalman filter's ability to recover parameters with different levels of measurement error.Large measurement errors make it difficult for the Kalman filter to isolate the true random process, making the parameters harder to recover.To show the effect of large measurement noise we run three sets of 200 simulations, with  = 10 −24 rad 2 s −2 , 10 −20 rad 2 s −2 and 10 −16 rad 2 s −2 .The results are shown in the three columns of Fig. D3.It is easiest to see the effect of varying  by looking at  and  c because they are the best recovered (except ⟨ Ω c ⟩ which is trivial).We can see in each test that  c is poorly recovered below some critical value, when the measurement noise makes up most of the total noise.If  c is above the critical value, parameter estimation is usually successful.

D3 Impact of incorrect measurement uncertainties
In this section we investigate the effect of feeding incorrect measurement uncertainties into the Kalman filter and introduce a modified parameter estimation algorithm to correct for the effect.
The frequency uncertainties for the real data in this paper are calculated from the TOA uncertainties provided in the UTMOST data release (Lower et al. 2020).Frequency uncertainties may be wrong, if the raw TOA errors or their conversion to frequency errors are incorrect.The quoted TOA uncertainties depend on many factors such as dispersion measure (Verbiest & Shaifullah 2018).
To simulate the problem, we generate data with a measurement error variance  true but use a different measurement error variance,  KF , in the Kalman filter.Results with  KF = 10 −23 rad 2 s −2 , R true = 10 −27 rad 2 s −2 and  KF = 10 −27 rad 2 s −2 , R true = 10 −23 rad 2 s −2 are shown in the left hand columns of Fig. D4 and Fig. D5 respectively.The results suggest that when  KF is too large (column 1 of Fig. D4) the  c values are underestimated and when  KF is too small (column 1 of Fig. D5) the  c values are overestimated.The effect is most pronounced for  c ≲ 10 −26 rad 2 s −3 .
The Kalman filter separates measurement noise from the underlying random process guided by  KF .If it is guided to remove too much noise ( KF ≥  true ), then some process noise is removed too, and the inferred strength of the process noise is weaker than it should be.The reverse is true for  KF ≤  true .If the measurement noise is not properly separated then  is also difficult to recover.
We adjust for  KF being input incorrectly by sampling  as well as the six dynamical model parameters , ,  c ,  s , ⟨Ω c − Ω s ⟩, ⟨ Ω c ⟩.The results are shown in the right columns of Fig. D4 and Fig. D5.The recovered  c values are now close to their correct values and the spread on the recovered  values is lower, which is an encouraging outcome.When analysing the real data in Section 5, sampling  makes little difference to the posterior, so the result is not displayed for brevity.We present this modified estimation method in this appendix in preparation for analysing more objects in the future.

APPENDIX E: IDENTIFIABILITY OF NOISE PARAMETERS
In this appendix we consider a simplified version of the parameter estimation problem to get a heuristic for the identifiability of the noise parameters  c and  s , as we did in Section 3.2 for the other parameters.We aim to separate the measurement noise from  c and  s by exploiting their different behaviours over time.Specifically,  c acts directly on the crust so it affects the frequency faster than  s , whose effect is delayed by the coupling time-scale.The measurement noise does not grow with time, nor does it evolve according to the equations of motion.The relative strengths of   ,  s and the measurement noise affect the size and shape of the observed random walk in Ω c .By looking at the rate that the random walk in Ω c grows over different time-scales, the strengths of the three noise types can be estimated.In this section we set  c =  s = 0 to analyse just the random behaviour.
Let Ω c be observed at some instant and then again a time Δ later, such that Ω c changes by an amount ΔΩ c .Because it is a random process, ΔΩ c is drawn from a random distribution with a variance we can calculate as a function of Δ.By comparing the observed distribution of ΔΩ c to the calculated distribution we can fit for the unknown parameters to infer  c ,  s and the size of the measurement noise.
The variance of a jump in Ω c over a time Δ is given by equation (B3), viz. with and For Δ ≫  one has and for Δ ≪  one has Hence the influence of  s through  s on short time-scales is negligible.
The variance for a step in the measurement can be calculated from (E3).Denote the measurement at time   by   and the measurement noise by   .Then one has In this appendix we assume for simplicity that all the measurement errors have variance  even though in reality different data points have different measurement uncertainties.Then the change in  between   and   has variance Let us define the auxiliary function   from the many recovered  values lying near the bottom of the prior range.
We also simulate data with the same TOA spacings and measure- The PP-plot in Fig. F3 shows that the parameter curves deviate significantly from the diagonal.However, the estimates do not seem to have a significant bias in either direction.So while the parameter estimates in Section 5 are uncertain due to the fitting procedure, they have no significant systematic bias.The large spread of the recovered parameters is due in large part to the small number of frequency data points.More data would give a greater level of convergence to the true values as seen in Fig. F1.

APPENDIX G: FULL POSTERIOR DISTRIBUTION FOR PSR J1359−6038
In this appendix, for the sake of completeness, we present in Fig. G1 a visualisation of the posterior distribution ( |) for all six parameters  = {, ,  c ,  s , ⟨Ω c − Ω s ⟩, ⟨ Ω c ⟩}.The figure is formatted as a traditional corner plot: the panels with contours display ( |) marginalised over four out of six parameters, and the panels with histograms display ( |) marginalised over five out of six parameters.The distinctly peaked histograms in columns 2, 3, 4 and 6 are reproduced in Fig. 2. The flatter histograms in columns 1 and 5 correspond to parameters that cannot be inferred reliably from the data, as discussed in Section 5.5.

APPENDIX H: ONE-COMPONENT MODEL
In this appendix, Bayesian model selection is performed to determine whether the two-component model outperforms a simpler onecomponent model when modelling PSR J1359−6038.
The model given by ( 1) and ( 2) assumes there is a second component of the star, hidden from view, which couples to the crust, viz. the superfluid with angular velocity Ω s .We could instead construct a simpler model where the pulsar is one rigid body described by the equation of motion where Ω is the angular velocity of the crust (to which the pulses are tied),  is the constant spin-down torque,  is the pulsar's moment of inertia, and  is the stochastic torque responsible for timing noise modelled as a white noise process with strength .
We can make a Kalman filter for the one-component model and apply the same method of parameter estimation as in Section 3.For the one-component model, the update and measurement equations are still given by ( 6)-( 9) with The dynesty sampler generates samples from the posterior distribution to create the posterior plots and computes (H8) as a by-product.
Let  1 and  2 denote the one-and two-component models respectively.We assert no prior preference for either model so set ( 1 ) = ( 2 ), yielding the evidence ratio (Bayes factor)  H1.The contour plot in Fig. H1 shows no evidence for correlations between ⟨ Ω c ⟩ and  c .
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .Figure 2 .
Figure 1.Residuals of frequencies fitted to 428 TOAs versus time (units: MJD) for the PSR J1359−6038 after subtracting a linear trend (Ω( ) [rad s −1 ] = 49.2888− 2.4472 × 10 −12  [s]) and removing two outliers due to poor fitting.The vertical scale is in units of 10 −8 rad s −1 .The frequencies and error bars are calculated using tempo2 from the TOAs and their uncertainties in the UTMOST data release.

Figure 3 .
Figure 3. Residuals of simulated pulsar frequencies to illustrate the qualitative difference between the timing noise produced by  c and  s .The first panel shows the fluctuations in Ω c for a simulated pulsar with  c ≠ 0 but  s = 0, the second panel shows Ω c produced for  s ≠ 0 and  c = 0 and the third panel is for a simulated pulsar, where  c and  s are calculated by adding the fluctuations from the first and second panels.

Figure D1 .
Figure D1.Test of the parameter estimation scheme on simulated frequency data showing the injected (horizontal axis) and recovered (vertical axis)  , ,  c ,  s , ⟨Ω c − Ω s ⟩ and ⟨ Ω c ⟩ values for 200 simulations.The injected parameters are those used to simulate the data and the recovered parameters are the peak of the posterior distribution obtained using the Kalman filter and nested sampler.The colour of each point can be compared to the colour bar next to its panel to get the width (interquartile range) of the marginalised posterior for that parameter.

Figure D3 .
Figure D3.Test of the parameter estimation scheme on simulated frequency data showing the injected (horizontal axis) and recovered (vertical axis)  and  c values (top and bottom rows respectively).Each column shows the results of 200 simulations with a different level of measurement error.Column 1 has  = 10 −24 rad 2 s −2 , column 2 has  = 10 −20 rad 2 s −2 , and column 3 has  = 10 −16 rad 2 s −2 .The colour of each point can be compared to the colour bar next to its panel to get the width (interquartile range) of the marginalised posterior for that parameter.

Figure D4 .
Figure D4.Test of the parameter estimation scheme on simulated frequency data showing the injected (horizontal axis) and recovered (vertical axis)  and  c values (top and bottom rows respectively) for 200 simulations with  true = 10 −27 rad 2 s −2 , R KF = 10 −23 rad 2 s −2 .The colour of each point can be compared to the colour bar next to its panel to get the width (interquartile range) of the marginalised posterior for that parameter.The two plots in the left-hand column show the results with the normal parameter estimation algorithm.The plots in the right-hand column show the results on the same data supplemented by sampling of  KF .

Figure F1 .
Figure F1.Accuracy of tempo2 fits to local frequencies.Plots of the injected (horizontal axis) and recovered (vertical axis)  and  c values (top and bottom rows respectively) for 200 simulations.The two plots in the left-hand column show the results when the method is applied to frequencies fitted to simulated TOAs with random spacings.The two plots in the right-hand column show the results when the method is applied to simulated frequencies without any tempo2 fitting.The colour of each point can be compared to the colour bar next to its panel to get the width (interquartile range) of the marginalised posterior for that parameter.

Figure F2 .Figure F3 .
Figure F2.Accuracy of tempo2 fits to local frequencies for synthetic data mimicking UTMOST observations of PSR J1359−6038.Plots of the injected (horizontal axis) and recovered (vertical axis)  and  c values (top and bottom rows respectively) for 200 simulations.The two plots in the left-hand column show the results when the method is applied to frequencies fitted to simulated TOAs with the same spacing as the real data from PSR J1359−6038 used in Section 5.The two plots in the right-hand column show the results when the method is applied to simulated frequencies with the same spacing and uncertainties as on the left, but without any tempo2 fitting.The colour of each point can be compared to the colour bar next to its panel to get the width (interquartile range) of the marginalised posterior for that parameter.

Figure G1 .
Figure G1.Six-dimensional joint posterior distribution for the two-component model for PSR J1359−6038.Marginalised posteriors for each parameter are shown as histograms on the diagonal.The lower half triangle of panels shows two-dimensional contour plots (with dark blue representing high probability) marginalised over all but two variables.
and Δ  =   −  −1 .We compare the two models' ability to explain the PSR J1359−6038 data by calculating the Bayesian evidence, , for each model given the data.We calculate the Bayesian evidence   for a model  from the posterior distribution via the formula  = ∫  ( |, ) ( |).(H8)

Table 1 .
Prior distributions introduced in Section 3.4 for the parameter estimation carried out in Figs.G1 and H1.The rightmost column indicates the prior distribution,  (•), on each parameter.U (, ) indicates a uniform distribution between  and .The priors of the first four (last two) parameters are log-uniform (uniform).