Measuring the Magnetic Dipole Moment and Magnetospheric Fluctuations of SXP 18.3 with a Kalman Filter

The magnetic dipole moment μ of an accretion-powered pulsar in magnetocentrifugal equilibrium cannot be inferred uniquely from time-averaged pulse period and aperiodic X-ray flux data, because the radiative efficiency η 0 of the accretion is unknown, as are the mass, radius, and distance of the star. The degeneracy associated with the radiative efficiency is circumvented if fluctuations of the pulse period and aperiodic X-ray flux are tracked with a Kalman filter, whereupon μ can be measured uniquely up to the uncertainties in the mass, radius, and distance. Here, the Kalman filter analysis is demonstrated successfully in practice for the first time on Rossi X-ray Timing Explorer observations of the X-ray transient SXP 18.3 in the Small Magellanic Cloud (SMC), which is monitored regularly. The analysis yields μ=8.0−1.2+1.3×1030Gcm3 and η0=0.04−0.01+0.02 , compared to μ=5.0−1.0+1.0×1030Gcm3 as inferred traditionally from time-averaged data assuming η 0 = 1. The analysis also yields time-resolved estimates of two hidden state variables, the mass accretion rate and the Maxwell stress at the disk–magnetosphere boundary. The success of the demonstration confirms that the Kalman filter analysis can be applied in the future to study the magnetic moments and disk–magnetosphere physics of accretion-powered pulsar populations in the SMC and elsewhere.


INTRODUCTION
The magnetic dipole moments µ of neutron stars span several orders of magnitude, from µ ∼ 10 33 G cm 3 for magnetars down to µ ∼ 10 26 G cm 3 for millisecond pulsars (Lyne & Graham-Smith 2012).In rotation-powered objects, µ is inferred from phase-connected pulse timing, assuming magnetic dipole braking (Goldreich & Julian 1969;Ostriker & Gunn 1969).The resulting µ values are accurate up to a factor of order unity stemming from uncertainties in the electrodynamics of the plasma magnetosphere and the angle between the rotation and magnetic axes (Spitkovsky 2006), if the fundamental assumption of magnetic dipole braking does indeed hold.In accretion-powered objects, it is harder to measure µ accurately (Mukherjee et al. 2015).Resonant electron cyclotron lines yield a direct measurement of the magnetic field strength in the line formation region near the stellar surface (Makishima 2003;Caballero & Wilms 2012;Revnivtsev & Mereghetti 2016;Konar 2017;Staubert et al. 2019) but are detected in relatively few objects.Otherwise one usually resorts to inferring µ from standard X-ray timing analysis and time-averaged statistics (Ho et al. 2014;Klus et al. 2014;Mukherjee et al. 2015), coupled with predictions from physical theories of accretion; see Figure 6 in Klus et al. (2014) and Figure 22 in D'Angelo (2017) for example.In some objects, the local cyclotron measurements and global X-ray timing estimates of µ are performed simultaneously and calibrated successfully against each other; see Figure 11 in Makishima et al. (1999).Measurements of µ by the above methods form a key input into population synthesis calculations, whose aim is to present a consistent picture of how µ is distributed at birth and evolves, as a neutron star ages (Faucher-Giguere & Kaspi 2006;Kiel et al. 2008).A popular approach for accretion-powered pulsars is to infer µ from time-averaged observations of the aperiodic X-ray flux L (t) and pulse period P (t), assuming magnetocentrifugal equilibrium (Ho et al. 2014;Klus et al. 2014;Mukherjee et al. 2015).However the aforementioned approach relies on assuming a value for the radiative efficiency, which is unknown at the outset and difficult to measure.Moreover, important time-dependent information is lost, when one analyses time-averaged data only, e.g.(anti)correlations like ⟨P (t) L (t ′ )⟩.Progress has occurred in certain directions, e.g.employing Monte Carlo simulations to estimate ⟨Y (t) L (t ′ )⟩, where Y (t) is the measured pulsed fraction; see Figure 11 in Coe et al. (2015) and Figure 2 in Yang et al. (2018) for example.Further, there is evidence of a break in the X-ray flux autocorrelation function (and hence power spectral density) near the Kepler frequency at the disk-magnetosphere boundary for X-ray sources near spin equilibrium (Revnivtsev et al. 2009;Revnivtsev & Mereghetti 2016).Interpreting all the above data is challenging because of uncertainties in the structure of the disk-magnetosphere boundary, which is complicated, time-dependent, and fundamentally three-dimensional, as revealed by modern numerical simulations (Romanova et al. 2003;Romanova & Owocki 2015).In summary, time-averaged observables such as ⟨L(t)⟩ and ⟨P (t)⟩ do not contain enough independent pieces of information to measure µ uniquely.
In this paper, we generalize previous magnetocentrifugal estimates of µ for accretionpowered pulsars by exploiting the additional information available in the fluctuations of L(t) and P (t) around their time-averaged values, breaking the degeneracy between µ and the radiative efficiency.To do so, we apply the signal processing framework developed by Melatos et al. (2023), which (i) utilizes a Kalman filter to track fluctuations in three hidden state variables associated with magnetocentrifugal accretion (the mass accretion rate, the Maxwell stress at the disk-magnetosphere boundary and the radiative efficiency); and (ii) maximizes the Kalman filter likelihood to estimate the underlying, static physical parameters, including µ.As this is the first time the method has been applied to astronomical data, we focus in this paper on demonstrating its feasibility using a single, regularly monitored object with sufficient samples of L(t) and P (t), namely the X-ray transient SXP 18.3 in the Small Magellanic Cloud (SMC) (Corbet et al. 2003).The results presented herein differ from previous studies in the following two ways.(i) The degeneracy between µ and the radiative efficiency is broken, as noted above, so that µ and the radiative efficiency are measured independently.(ii) The Kalman filter framework exploits the auto-and cross-correlations between P (t) and L(t) to place multiple, simultaneous, interlocking constraints on the hidden state variables describing magnetocentrifugal accretion, breaking the degeneracies between them and producing a unique solution for their temporal evolution, corresponding to the most likely sequence of hidden states consistent with the time-ordered observational data.This approach preserves more statistical information than the alternative, which is to calculate P (t) and L(t) statistics (e.g.power spectral densities) separately and relate the associated ensemble averages to physical models of accretion (Bildsten et al. 1997;Riggio et al. 2008;Klus et al. 2014;Ho et al. 2014;Mukherjee et al. 2015;Serim et al. 2023).It builds on earlier work applying autoregressive moving average models to torque-luminosity correlations (Baykal & Oegelman 1993), and testing for consistency with random walk and shot noise processes (de Kool & Anzer 1993;Baykal 1997;Lazzati & Stella 1997).
The paper is organised as follows.In Section 2 we introduce the stochastic differential equations of motion which govern how the state variables describing accretion evolve in the context of the idealized, canonical, magnetocentrifugal picture (Ghosh et al. 1977), as well as the measurement equations which map the state variables (some of which are hidden) to the measured aperiodic X-ray flux L (t) and pulse period P (t).The equations of motion and measurement equations are linearized about magnetocentrifugal equilibrium and cast into a state-space formulation suitable for Kalman filter analysis.In Section 3 we introduce the SMC X-ray transient SXP 18.3.We briefly describe its properties and the Rossi X-ray Timing Explorer (RXTE) observations used throughout the analysis.We choose to focus our analysis on SXP 18.3 for two practical reasons: (i) the rotational state of the star has been classified by Yang et al. (2017) as being approximately in equilibrium; and (ii) the RXTE Proportional Counter Array (PCA) sensor collected 854 observations over a period of ∼ 16 years, comparable in volume to the representative test case involving synthetic data presented in Melatos et al. (2023).New estimates of the magnetic dipole moment µ and radiative efficiency η 0 associated with SXP 18.3 are presented in Section 4. Astrophysical implications are canvassed briefly in Section 5 together with a note about future plans for generalizing the signal processing framework to population studies.

MEASURING THE MAGNETIC MOMENT
Disk accretion onto a magnetized, compact star is a complicated, time-dependent process involving hydromagnetic instabilities at the disk-magnetosphere boundary, which imprint a rich, three-dimensional, spatial structure on the system, as seen in magnetohydrodynamic numerical simulations (Romanova et al. 2003(Romanova et al. , 2005;;Kulkarni & Romanova 2008).The spatial structure in the simulations, e.g. the geometry of the magnetic field near the disk-magnetosphere boundary, cannot be inferred uniquely from the observables P (t) and L (t).In this paper, therefore, we choose to describe accretion in terms of the canonical, spatially averaged, magnetocentrifugal model (Ghosh & Lamb 1979), motivated by the promising results presented in Melatos et al. (2023).The reader is directed to the latter reference for a full account of the model, its physical justification, and its implementation in the form of a Kalman filter; we do not repeat here the in-depth discussion in Melatos et al. (2023).Instead, in this section, we introduce briefly the variables, definitions, and dynamical equations needed to implement the Kalman filter, so that the interested reader is equipped to reproduce the key results of this paper, namely the independent measurements of µ and the radiative efficiency for SXP 18.3.In Section 2.1, we relate the observables P (t) and L (t) via nonlinear measurement equations to four time-dependent state variables, three of which are hidden.In Section 2.2, we discuss (i) the canonical magnetocentrifugal torque law; (ii) the degeneracy that exists between µ and the radiative efficiency for systems near magnetocentrifugal equilibrium; and (iii) a phenomenological, idealized model of the stochastic forces associated with instabilities at the disk-magnetosphere boundary, which drive the system away from equilibrium, first introduced in Section 2.4 in Melatos et al. (2023).The nonlinear measurement equations in Section 2.1 and equations of motion in Section 2.2 are linearized about the state of magnetocentrifugal equilibrium in Section 2.3, producing a state-space formulation of the problem which is suitable for Kalman filter analysis.The formulation in Section 2.3 breaks the degeneracy between µ and the radiative efficiency.An explicit formula is presented to calculate µ as a function of the parameters estimated by the Kalman filter from the measurements P (t) and L(t) in Section 2.4.A short review of the key elements of a linear Kalman filter is presented in Appendix A for the convenience of the reader.Details can be found in Melatos et al. (2023).

Observables and state variables
In its simplest form, an X-ray timing experiment returns raw photon counts, which are barycenter corrected and converted to N samples of the pulse period, P (t 1 ) , . . .P (t N ), and aperiodic X-ray luminosity, L (t 1 ) , . . ., L (t N ), using standard X-ray timing analysis; see Section 2.3 in Yang et al. (2017) for practical details about the data reduction pipeline. 1e express the standard magnetocentrifugal model of disk accretion (Ghosh & Lamb 1979) in terms of four time-dependent state variables.The angular velocity Ω (t) of the neutron star is inversely related to the pulse period P (t).The other three state variables, defined as follows, are hidden.Let Q (t) be the rate at which matter flows from the accretion disk into the disk-magnetosphere boundary (units: g s −1 ).Let S (t) denote the Maxwell stress at the disk-magnetosphere boundary (units: g cm −1 s −2 ).Let η (t) denote the dimensionless radiative efficiency with which the gravitational potential energy of material falling onto the star is converted into X-rays.
Hidden state variables cannot be measured directly.Rather, they are related indirectly to the observables P (t) and L(t) through algebraic equations, which are nonlinear in general.The measured pulse period P (t) and aperiodic X-ray luminosity L(t) are given respectively by and where Newton's gravitational constant and the mass and radius of the neutron star are denoted by G, M , and R respectively We assume that the additive noise terms N P (t) and N L (t) are Gaussian and white, with , and ⟨N P (t n ) N L (t n ′ )⟩ = 0, where δ n,n ′ denotes the Kronecker delta.
The radiative efficiency 0 < η(t) < 1 is the product of two factors: (i) the fraction of the specific gravitational potential energy GM/R that is converted into X-rays in the observed band through radiation processes such as bremsstrahlung emission; and (ii) the fraction of Q(t) which strikes the surface of the star, which does not equal unity in general, because some of Q(t) is diverted into polar outflows (for example) by nonconservative hydromagnetic instabilities at the disk-magnetosphere boundary (Matt & Pudritz 2005, 2008;Romanova & Owocki 2015).It is a severe simplification to combine (i) and (ii) into a single scalar variable η(t), as discussed in detail in Appendix C in Melatos et al. (2023).However, more realistic models (D'Angelo & Spruit 2010; D'Angelo & Spruit 2012; D'Angelo 2017) require more data than the N = 854 samples available for SXP 18.3.In the present analysis, therefore, we stick with η (t) and make the additional simplification η (t) = η 0 = constant, i.e. we hold η (t) constant in time but leave its value free to be estimated by the Kalman filter from the data.This approach is appropriate for the volume of data available.Despite its simplicity, it yields new astrophysical insights, namely the first independent measurement of η 0 and hence µ for SXP 18.3, as described in Section 4.

Magnetocentrifugal accretion dynamics
The accretion torque acting on the star is governed by the location of the diskmagnetosphere boundary.In reality, the boundary has a complicated internal structure with nonzero thickness, which affects the torque and other dynamical variables [e.g.η(t)] in different ways (D'Angelo & Spruit 2010;D'Angelo & Spruit 2012;D'Angelo 2017); see also Appendix C in Melatos et al. (2023).In this paper, guided by simplicity and the modest volume of data available, we adopt the canonical magnetocentrifugal model (Ghosh et al. 1977), and approximate the boundary as having zero thickness and occurring at the Alfvén radius R m (t), where the Maxwell stress S(t) balances the disk ram pressure, viz.S ≈ ρv2 /2, where ρ = Q/ (4πR 2 m v) and v = (2 GM/R m ) 1/2 are the mass density and infall speed 2 respectively in free fall (Menou et al. 1999;Frank et al. 2002).The Alfvén radius is given in terms of the state variables by (3) The corotation radius, is located where Keplerian disk matter corotates with the neutron star.Reading Equation (3) at face value, it may seem that R m (t) increases with Q(t), which is the opposite of what one expects physically.However, S(t) = (2π) −1 µ 2 R m (t) −6 contains R m (t) implicitly.Upon substituting the latter expression into Equation (3), we recover the standard expression (up to a factor of order unity) for the Alfvén radius in terms of the mass accretion rate and magnetic dipole moment; see for example Equation (7) in Klus et al. (2014).
The sign of the magnetocentrifugal torque is governed by the ratio R m /R c , known as the fastness parameter.Material strikes the stellar surface for R m < R c , and the neutron star spins up in response to a combination of hydromagnetic and mechanical torques, with dP/dt < 0. The propeller mechanism occurs for R m > R c .Matter is ejected centrifugally by the corotating magnetosphere, causing the star to spin down, with dP/dt > 0. For , and R m = R c are described approximately by the canonical magnetocentrifugal torque law (Ghosh et al. 1977) where the star's moment of inertia is denoted by I. Equations (3)-( 5) feature the hidden variables Q (t) and S (t), but not η (t).For the sake of brevity, the reader is referred to (i) Section 2.2 in Melatos et al. (2023) for a discussion about corrections to Equation ( 5), e.g.due to radiation pressure (Andersson et al. 2005;Haskell et  Assuming a magnetic dipole field, and imposing the zero-torque condition, R m = R c = constant, the star's magnetic moment µ is given by3 While we assume for simplicity that µ is constant, there is observational evidence to suggest that the magnetic dipole moments of recycled neutron stars decrease monotonically with accreted mass (Shibazaki et al. 1989).Several mechanisms have been proposed to explain µ = µ (t) ̸ = constant, including accelerated ohmic decay (Urpin & Geppert 1995;Urpin & Konenkov 1997) and diamagnetic screening or burial (Uchida & Low 1981;Hameury et al. 1983;Brown & Bildsten 1998;Zhang 1998;Melatos & Phinney 2001;Cumming et al. 2001;Rai Choudhuri & Konar 2002;Payne & Melatos 2004).
The zero-torque condition and the two time-averaged observables, ⟨2π/P (t)⟩ = Ω 0 and ⟨L(t)⟩ = L 0 , do not contain enough independent pieces of information to uniquely solve for all four components of the state vector in equilibrium (Ω 0 , Q 0 , S 0 , η 0 ).Therefore one is left with two options: (i) assume a plausible value for one component, e.g.η 0 = 1; or (ii) exploit the time-dependent information in the observations P (t) and L(t), in addition to ⟨P (t)⟩ and ⟨L(t)⟩, and treat one state vector component as an unknown to be estimated from the data.There are at least two approaches available to utilize the time-dependent information in P (t) and L(t), namely χ 2 minimization (Takagi et al. 2016;Yatabe et al. 2018) and Kalman filter analysis (Melatos et al. 2023).Motivated by the results in Melatos et al. (2023), we adopt the latter technique, yielding the first independent estimate of η 0 and hence µ for SXP 18.3.In Appendix C we follow Takagi et al. (2016) and Yatabe et al. (2018) and also employ χ 2 minimization to estimate µ independently for SXP 18.3 and compare the output of both techniques for completeness.
Observations of accretion-powered pulsars reveal pronounced, stochastic variability in the measured aperiodic X-ray flux L (t) and pulse period P (t) (Bildsten et al. 1997).The variability is driven by several mechanisms, discussed in detail in Section 2.4 in Melatos et al. (2023).Here we focus on mean-reverting fluctuations about magnetocentrifugal equilibrium, which are typical of self-limiting instabilities and are observed routinely in the output of three-dimensional magnetohydrodynamic simulations (Romanova et al. 2005(Romanova et al. , 2021)).Accordingly, we adopt an idealized, phenomenological, Ornstein-Uhlenbeck model (Gardiner et al. 1985) for the state variables Q(t) and S(t), and assume that they satisfy the Langevin equations (Melatos et al. 2023) In Equations ( 7) and ( 8), γ Q and γ S are damping constants, and ξ Q (t) and ξ S (t) are fluctuating driving terms satisfying ⟨ξ 7) and ( 8) are supplemented with the deterministic equation where the reader is reminded that η 0 is constant but free to be estimated with the Kalman filter in Section 4. The Langevin equations ( 7) and ( 8) ensure that Q (t) and S (t) wander randomly about their equilibrium values, and do not undergo long-term secular drifts, with characteristic timescales of mean reversion given by γ −1 Q and γ −1 S , and rms fluctuations ∼ γ σ SS respectively.The statistics associated with ξ Q (t) and ξ S (t) are discussed in more detail in Appendix A.
Equations ( 7) and ( 8) are highly idealized.They omit some important pieces of microphysics operating at the disk-magnetosphere boundary (Romanova & Owocki 2015), which are discussed in detail in Section 2.4 and Appendix C of Melatos et al. (2023).They compound the idealizations inherent in the torque law (5), discussed above, and in Equation ( 9) through the definition and constancy of η(t) (see Section 2.1).Altogether, Equations ( 5) and ( 7)-( 9) represent a compromise between realism and the explanatory power of the available volume of data, specifically N = 854 samples for SXP 18.3.

Linear state-space formulation
The Kalman filter is a standard tool for inferring the most probable trajectory of a set of stochastic hidden variables, such as Ω(t n ), Q(t n ), and S(t n ), given a data stream of noisy measurements, such as P (t n ) and L(t n ), a set of measurement equations, such as (1) and (2), and a set of stochastic evolution equations, such as ( 5), (7), and (8) (Gelb et al. 1974;Wan & Van Der Merwe 2000).The key elements of the linear Kalman filter are reviewed in abridged form in Appendix A for the convenience of the reader; details specific to this application can be found in Melatos et al. (2023).In this paper, the measurement equations ( 1) and ( 2) and the magnetocentrifugal torque law (5) are nonlinear.Hence one must choose between two approaches: (i) apply a nonlinear estimation algorithm, e.g. a particle filter (Carpenter et al. 1999;Gustafsson 2010), or an extended or unscented Kalman filter; see Sections 2 and 3 in Wan & Van Der Merwe (2000) and Section 2 of Kandepu et al. (2008); or (ii) linearize Equations ( 1), (2), and (5) about magnetocentrifugal equilibrium, under the assumption that fluctuations are small, i.e. γ −1/2 A σ AA ≪ A 0 , with A ∈ {Q, S}.As this paper reports the first application of the Kalman filter framework to astronomical data, and aims to supply a proof of principle, we adopt approach (ii) to be conservative and defer approach (i) to future work.Nonlinear analyses are straightforward in principle but are more expensive computationally; see Appendices B and C in Melatos et al. (2023) for illustrative examples applied to synthetic data.
Linearizing the equations of motion ( 5), ( 7), and ( 8), and introducing the notation for perturbed quantities, we obtain (Melatos et al. 2023) with Linearizing the measurement equations ( 1) and (2), we obtain and Equations ( 10 13), ( 16), ( 19), and (D11) in Melatos et al. (2023) in order to keep η(t) constant in time but free to be estimated.The simplifications are justified by the limited sample sizes available with modern X-ray timing experiments, i.e. 854 samples for SXP 18.3.Although we elect to track X(t) = [Ω(t), Q(t), S(t), η(t)] with a linear Kalman filter in the present paper, the specific choice of state-space variables is not unique.For example, one may choose to replace Equation ( 8) with an Ornstein-Uhlenbeck process for the Alfvén radius, postulating mean-reverting stochastic fluctuations in R m (t) instead of S(t).In practice, it is unclear from theory or observations what approach better approximates the accretion dynamics in real systems.To support future Kalman filter applications using alternative state-space variables, a linear statespace formulation using R m (t) in lieu of S(t) is presented in Appendix B for the benefit of the reader.

Magnetic dipole moment from the Kalman filter output
The five fundamental model parameters estimated by the Kalman filter from the data are Θ = (γ Ω , γ A , σ AA ), with A ∈ {Q, S}.One of these, namely γ Ω , is combined with the sample means Ω 0 = 2 3/2 π 3/5 (GM ) 1/5 Q −3/5 0 S 3/5 0 and L 0 = GM Q 0 η 0 /R, and the magnetocentrifugal equilibrium condition, to uniquely solve for (Ω 0 , Q 0 , S 0 , η 0 ), adopting approach (ii) in Section 2.2.The explicit expressions for Q 0 , S 0 , and η 0 follow from the step-by-step guide outlined in Appendix A in Melatos et al. (2023) with one minor modification.The Maxwell stress component, used in this paper is 2 −1/2 times its value in Equation (A5) in Appendix A in Melatos et al. (2023); see Equation ( 3) and Footnotes 2 and 3 for details.The expressions for Q 0 and η 0 remain unchanged and are not quoted here for the sake of brevity; see Equations ( A4) and (A6) in Melatos et al. (2023).The magnetic moment µ then follows from Equation (6).

SXP 18.3: A DEMONSTRATION
The SMC is a dwarf irregular galaxy orbiting the Milky Way at a known distance of 62 ± 0.3 kpc (Scowcroft et al. 2016).It hosts more than 120 known high-mass X-ray binaries (HMXBs) (Haberl et al. 2022); see White et al. (1995) and Reig (2011) for overviews of X-ray binaries.The optical counterpart is usually a Be star, with only one system hosting a supergiant donor star, namely SMC X-1 (Haberl & Sturm 2016).The large population of SMC Be-type HMXBs is ascribed to: (i) recent star formation (Galache et al. 2008); (ii) low metallicity (Antoniou et al. 2010); and (iii) recent tidal interactions with the Large Magellanic Cloud (Gardiner & Noguchi 1996).The orbital parameters and stellar properties of SMC HMXBs have been studied extensively (Galache et al. 2008;Townsend et al. 2011;Klus et al. 2014;Coe et al. 2015;Christodoulou et al. 2016Christodoulou et al. , 2017) ) with space-based X-ray observatories, e.g. the RXTE PCA (Jahoda et al. 2006).In Section 3.1 we discuss the population of SMC HMXBs analyzed by Yang et al. (2017), and introduce the X-ray transient SXP 18.3.Details of the RXTE PCA measurements of SXP 18.3 analyzed with the Kalman filter are presented in Section 3.2.

Source properties
Archival X-ray time-of-arrival measurements collected between 1997 and 2014 by the RXTE PCA, XMM-Newton European Photon Imaging Camera (Strüder et al. 2001), and Chandra Advanced CCD Imaging Spectrometer (Garmire et al. 2003) are publicly available via the High Energy Astrophysics Science Archive Research Centre. 4hey are also post-processed into pulse period, pulse amplitude, pulsed fraction, and aperiodic X-ray luminosity time series by Yang et al. (2017).The latter authors estimate the pulse period derivative dP/dt from a linear fit of P (t n ) versus t n for 53 stars with enough data.They characterise the torque distribution of the population in terms of dP/dt divided by its measurement error, denoted by ϵ in Table 3 in Yang et al. (2017).They classify 12 pulsars as spinning up (ϵ ≤ −1.5), 11 as spinning down (ϵ ≥ 1.5), and 30 as being near spin equilibrium (−1.5 < ϵ < 1.5), labelled by the letters "U", "D", and "C" respectively in Table 3 in Yang et al. (2017).
Using the post-processed time series in Yang et al. (2017), we select a subpopulation of accretion-powered pulsars which satisfy two conditions.(i) The number of RXTE PCA samples is ≥ 850, which is sufficient for the Kalman filter to converge accurately, as demonstrated by Melatos et al. (2023); see Table 2 in Yang et al. (2017).(ii) There is no long-term secular drift in the spin period, i.e. each source is near equilibrium; see Table 3 in Yang et al. (2017).Although the Kalman filter framework can be applied to accretion-powered pulsars in magnetocentrifugal disequilibrium, as confirmed for synthetic data in Appendix B of Melatos et al. (2023), we elect to focus on nearequilibrium systems in this paper to be conservative, as this is the first time the framework is applied to real astronomical data.Upon imposing criteria (i) and (ii), we are left with a subsample of nine objects, whose names and timing properties are listed in Table 1.
We focus our analysis on SXP 18.3 for three practical reasons.(i) It has the largest number of significant detections, N det = 73, in Table 1.The concept of significance is defined below.(ii) The standard deviation of the pulse period satisfies σ P = 0.025 s ≪ P 0 (see Table 1), consistent with the assumption of linearity.(iii) It accretes via a disk, consistent with the canonical picture of magnetocentrifugal accretion.With respect to point (iii) above, Klus et al. (2014) investigated whether the X-ray emission is powered by an accretion disk or stellar wind (Shakura et al. 2012) for 42 SMC HMXBs, including SXP 18.3.They estimated the relative velocity V rel between the accreted matter, calculated from the stellar wind velocity at the companion star's circumstellar disk, and the Keplerian orbital velocity of the neutron star.They compared V rel with the critical relative velocity V Crel , below which disk accretion occurs (i.e.R m < R c ), and found V rel /V Crel < 1 for the 42 HMXBs, indicating that all 42 systems accrete via a disk instead of a wind; see Section 3 and Figure 5 in Klus et al. (2014) for further details.With respect to point (i) above, the statistical significance of each timing point is determined by searching the light curve for pulsations and calculating the number of independent frequencies and the Lomb-Scargle power according to Equation (2) in Yang et al. (2017); see Section 2.4.1 and Table 3 2014) estimated the surface magnetic field strength B = 5.0 +1.0 −1.0 × 10 12 G of SXP 18.3 and hence inferred µ = 5.0 +1.0 −1.0 × 10 30 G cm 3 based on Q 0 ∝ ⟨L(t)⟩ and η 0 = 1; see Section 4 and Table 3 in Klus et al. (2014) for further details.

RXTE data
The RXTE PCA is a non-imaging detector (Corbet et al. 2008).The recorded photon counts are converted to light curves using standard X-ray timing analysis; see Van der Klis (1989) for further details, and Sections 1.3 and 4.4.1 in Laycock et al. (2005) and Patruno & Watts (2021) respectively.Once extracted, the light curves are referred to the Solar System barycenter with background noise subtracted (Yang et al. 2017).To estimate the spin period, the light curves are searched for significant pulsations using the Lomb-Scargle periodogram (VanderPlas 2018), and the spin period and pulse amplitude are retrieved as observables.Multiple SMC Xray pulsars regularly occupy the field of view of the RXTE PCA (Yang et al. 2017).In order to avoid source confusion, the aperiodic X-ray flux is not measured directly by the RXTE PCA.Rather, it is estimated empirically from the pulse amplitude; see Section 2.4.2 and Eq. ( 5) in Yang et al. (2017).
In this paper we analyze the data presented in Yang et al. (2017).The time series P 1 (t n ) and L 1 (t n ) are plotted as grey points in the top two panels of Figure 1.Samples whose significance exceeds 99% (defined in Section 3.1) are overplotted as green points.There are a number of important qualitative and quantitative differences between the synthetic data analyzed by Melatos et al. (2023) and the real data in Figure 1 which should be taken into account when comparing the accuracy of the Kalman filter in Section 4 with the results in Melatos et al. (2023).(i) The fractional fluctuations for SXP 18.3 are approximately three orders of magnitude greater than for the synthetic P (t) data and between one and two orders of magnitude greater than for the synthetic L(t) data; compare the first and second panels in Figure 1 with Figure 1 in Melatos et al. (2023).(ii) Only 8% of the measurements in Figure 1 (overplotted as green points) are classified by Yang et al. (2017) as significant pulsations, whereas all 500 synthetic samples in Melatos et al. (2023) are implicitly assumed to be significant.(iii) Between MJD 53750 and MJD 55000, outbursts are visible in the second panel of Figure 1, whereas no outbursts are visible in Figure 1 in Melatos et al. (2023).Outbursts are common in X-ray timing studies and coincide with many significant detections.The increase in L 1 (t) between MJD 53975 and MJD 54175 in Figure 1 is the longest type II outburst recorded for an SMC X-ray pulsar (Schurch et al. 2009).

KALMAN FILTER ANALYSIS
A linear Kalman filter (Kalman 1960;Gelb et al. 1974) combined with a nested sampler (Speagle 2020) can be applied to the measured time series P 1 (t n ) and L 1 (t n ) to estimate the posterior distribution of the five fundamental model parameters Θ = (γ Ω , γ A , σ AA ), with A ∈ {Q, S}, conditional on the linear model described by Equations ( 10), ( 12), and (13).The Kalman filter likelihood, recursion relations, and nested sampling algorithm are discussed in Section 4.1 and Appendix A. We discuss the hidden state evolution reconstructed by the Kalman filter in Section 4.2.New estimates of the magnetic dipole moment µ and radiative efficiency η 0 of SXP 18.3 are presented in Section 4.3.Some consequences of specific systematic uncertainties, e.g.due to the misalignment of the magnetic and rotation axes of the star, are discussed in Section 4.4 for completeness.The random dispersions of the estimated parameters including µ and η 0 are quantified in Section 4.5.

Algorithm
In this paper we employ the dynesty nested sampler (Speagle 2020) with the bilby front-end (Ashton et al. 2019) to evaluate the Kalman filter log-likelihood (Meyers et al. 2021) where is the measurement vector, D Y = 2 denotes the dimension of Y n , and the innovation vector and its covariance are denoted by e n and s n = ⟨e n , e T n ⟩ respectively.In its simplest form, the nested sampling algorithm proceeds as follows.At every step k in the iterative sampling process, the sampler keeps track of N live "live points" within the parameter domain, denoted by The live points are initialized by drawing Θ m for 1 ≤ m ̸ = m ′ ≤ N live .The process repeats until a suitable stopping condition is met, yielding a numerical approximation of the Bayesian evidence and posterior distributions of Θ via weighted histograms or other density estimation techniques.The reader is referred to Sections 5 and 6 in Skilling (2006) for an overview of the nested sampling algorithm and Table 2 in Ashton et al. (2022) for a comparison of nested sampling packages.
The Kalman filter generates the pre-fit innovation vector at every time step t n from where C denotes the 2×3 observation matrix defined implicitly by the noiseless terms on the right-hand sides of Equations ( 12) and ( 13), and A denotes the 3 × 3 state transition matrix defined on the right-hand side of Equation ( 10).Replacing Xn−1 with the updated state estimate Xn in Equation ( 16) yields the post-fit innovation vector, which is an important metric in estimating the accuracy of the Kalman filter state tracking.The state vector is updated recursively via where k n denotes the Kalman gain.Equations ( 16) and ( 17) are standard textbook formulas reproduced here for overall context; the reader is referred to Appendix A for further details regarding their structure and implementation.

Hidden state tracking
In Figure 1 we present the Kalman filter inputs and outputs as functions of time.The time series P 1 (t n ) and L 1 (t n ) are plotted in the top two panels as grey points and the Kalman filter estimates are overplotted as colored, solid curves.The reconstructed measurements are generated at each time step t n from the a posteriori state estimates Xn multiplied by C, defined in Section 4.1 above.We assess the accuracy with which the Kalman filter tracks the state evolution and hence P 1 (t n ) and L 1 (t n ) through the average rms error of the post-fit innovation vector using (i) the entire data set, and (ii) the subset of data corresponding to significant detections.In case (i) the average rms errors are ≈ 0.9σ P 1 and ≈ 0.6σ L 1 , where σ P 1 and σ L 1 denote the standard deviations of the entire P 1 (t n ) and L 1 (t n ) data set, respectively.In case (ii) the average rms errors are ≈ 1.0σ ′ P 1 and ≈ 0.8σ ′ L 1 where σ ′ P 1 and σ ′ L 1 denote the standard deviations of P 1 (t n ) and L 1 (t n ) respectively when restricted to the significant detections only.The results for case (ii) quantify the accuracy with which the Kalman filter tracks the state evolution during outbursts.In the absence of independent measurements of the true hidden state sequence X n , the post-fit innovation is a standard metric used in signal processing applications to assess the accuracy of a Kalman filter; see Section 10.1 in Simon (2006) for an overview on verifying Kalman filter performance.Ultimately, a fuller analysis of accuracy requires N det > 73.However, it is encouraging that the Kalman filter partly recovers the two outbursts visible in the second panel of Figure 1 between MJD 53750 and MJD 55000, even though the idealized state-space model in Section 2.3 does not treat outbursts explicitly.
The Kalman filter state estimates Xn = [ Ω(t n ), Q(t n ), Ŝ(t n )] are plotted as colored, solid curves in the bottom three panels of Figure 1.The shaded regions correspond to the 1σ state estimate uncertainties returned by the Kalman filter, i.e. they are given by the square root of the diagonal elements in Equation (A11).The mean-reverting nature of the hidden state variables is clear upon inspection.The components Ω 0 , Q 0 , and S 0 of the state vector in equilibrium are overplotted as black, dashed lines.We calculate the sample mean Ω 0 from Equation (A1) in Melatos et al. (2023) and infer Q 0 and S 0 from the mode of the posterior of γ Ω , discussed in Section 2.4.The accretion rate fluctuations increase from |Q 1 (t)|Q 0 ≲ 10 −9 M ⊙ yr −1 during quiescence to |Q 1 (t)|Q 0 ≲ 10 −8 M ⊙ yr −1 during outburst which is broadly consistent with other sources; see Figure 3 in Mukherjee et al. (2018) for example.The Maxwell stress fluctuations satisfy |S 1 (t)|S 0 ≲ 10 7 g cm −1 s −2 .They correlate somewhat with the pulse period and the aperiodic X-ray flux in general and during outbursts in particular; we measure the Pearson r-coefficients as r[S 1 (t), L 1 (t)] ≈ 0.19 ± 0.034 and r[S 1 (t), P 1 (t)] ≈ −0.36 ± 0.032, where the uncertainties correspond to the standard error on r.We also find evidence of a moderate correlation between the Kalman filter output for S 1 (t n ) and Q 1 (t n ), with r[S 1 (t), Q 1 (t)] ≈ 0.26 ± 0.033.An instance of the correlation is visible in Figure 1 between MJD 53975 and MJD 54175; the magnetosphere is compressed relative to equilibrium, with Q(t) > Q 0 , causing spikes in S(t).The results about S 1 (t) are examples of the new and physically important information that flows from a Kalman filter analysis.
In the bottom panel of Figure 1, S (t) goes negative near MJD 51700, MJD 53200, MJD 53400, and MJD 54600.Negative excursions of S(t) are unphysical, but they are permitted in principle (albeit rarely) by the Ornstein-Uhlenbeck dynamics in ( 7) and ( 8) in the linear regime; see Section 2.4 in Melatos et al. (2023), where the possibility is discussed.The negative excursions of S(t) are brief, i.e. 0.8 % of the total observation time.Enforcing S(t) > 0 artificially for all t does not change the state estimates appreciably.6

Magnetic moment and radiative efficiency
Figure 2 presents the five-dimensional posterior distribution of the model parameters Θ = (γ Ω , γ A , σ AA ), with A ∈ {Q, S}, returned by the nested sampler and visualised in cross section through a traditional corner plot.The parameter estimates are reported at the top of each one-dimensional posterior (histogram; log 10 scale).We quantify the uncertainty in the estimate using a 68% credible interval, visible as three, vertical, dashed lines in the one-dimensional posteriors, and through the innermost contour in the two-dimensional posteriors in Figure 2. The nominal value corresponds to the posterior median.The mode (peak) of each posterior is contained within the 68% credible interval, and is employed as a point estimate to generate the state sequence Xn plotted in Figure 1.The steps required to reproduce the results in Figures 1   and 2 are outlined in Section 3.1 in Melatos et al. (2023).We assume η (t) = η 0 in the present paper, thereby neglecting the parameters γ η and σ ηη in the Ornstein-Uhlenbeck equation for η(t); see Equation ( 11) in Melatos et al. (2023).
In Figure 2, four parameters are estimated unambiguously, namely γ Ω , γ Q , γ S , and σ QQ .There is evidence of γ Ω -γ S and γ Q -σ QQ correlations, visible as a diagonal tilt in the contours in the γ Ω -γ S and γ Q -σ QQ planes.In contrast, we observe railing in the posterior for σ SS , which cannot be estimated unambiguously.Although the identifiability analysis presented in Appendix D in Melatos et al. (2023) suggests that all five model parameters including σ SS should be identifiable with sufficient data, the estimation problem in this paper using real observational data is more challenging than the synthetic exercise in Melatos et al. (2023), e.g.compare the second panel in Figure 1 in Melatos et al. (2023) with Figure 1 below.Figure 2 suggests that one may require N > 854 for SXP 18.3 to constrain all five model parameters.
Figure 3 presents the posterior of the four inferred parameters (Q 0 , S 0 , η 0 , µ) and is visualised in cross section through a traditional corner plot.To generate the histograms we draw random samples from the marginalized posterior of γ Ω and the normal distribution A ∼ N (A 0 , σ A N −1/2 ), with A ∈ {Ω, L}, where σ A N −1/2 denotes the standard error of the sample means, Ω 0 and L 0 .The strong pairwise covariances (diagonal contours) observed in Figure 3 follow directly from the ratios in Equations ( A4) and (A6) in Melatos et al. (2023) and Equation ( 14) above.The parameter estimates inferred in this paper are reported in the first row in Table 2.We find Q 0 = 2.4 +0.8 −0.7 × 10 17 g s −1 , S 0 = 5.6 +2.0 −1.6 × 10 6 g cm −1 s −2 , η 0 = 0.04 +0.02 −0.01 , and hence µ = 8.0 +1.3 −1.2 × 10 30 G cm 3 , compared to µ = 5.0 +1.0 −1.0 × 10 30 G cm 3 as inferred traditionally from time-averaged data assuming Q 0 ∝ ⟨L(t)⟩ and η 0 = 1 (bottom row of Table 2).The latter, traditional approach implies that 100% of the gravitational potential energy of material falling onto the stellar surface is converted to heat and hence Xrays.It is a standard assumption in X-ray timing analysis.The results presented in this paper suggest otherwise, with η 0 ̸ = 1 in general, i.e. we find that only ≈ 4% of the specific gravitational potential energy of infalling matter in SXP 18.3 is converted to X-rays.As a result, the Kalman filter estimates of Q 0 and hence µ are higher than the estimates based on time-averaged data and η 0 = 1.The present analysis yields independent measurements of η 0 and hence µ in the X-ray pulsar context and is another example of a new and astrophysically important result which stems from the Kalman filter framework.
We draw the reader's attention to the following point of interest.In Klus et al. (2014) andD'Angelo (2017), the authors estimated B and hence µ for 38 SMC Xray pulsars.Assuming magnetocentrifugal equilibrium, Klus et al. (2014) found that ≈ 50% of the stars have B ≥ 4.4 × 10 13 G, i.e. above the Schwinger limit and comparable to magnetar field strengths, based on Q 0 ∝ ⟨L(t)⟩ and η 0 = 1.However, one must be cautious about interpreting B at face value.D'Angelo (2017) argued that the Q 0 estimates in Klus et al. (2014) correspond to the average mass accretion rate during outburst and suggested Q * 0 = (N det /N ) Q 0 as a proxy for the mass accretion rate during quiescence, where the ratio N det /N approximates the fraction of time spent in outburst.Adopting Q * 0 instead of Q 0 yields ≈ 15% of SMC X-ray sources with magnetar field strengths; see Figure 22 in D'Angelo (2017) for a comparison of B estimates using Q 0 and Q * 0 , plotted as green triangles and cyan points respectively.Similarly, the free parameter γ Ω and hence the components Q 0 , S 0 , and η 0 estimated in this paper are influenced by episodes of quiescence and type II outbursts (Reig 2008;Franchini & Martin 2021); see the second panel in Figure 1.Following D'Angelo (2017), we report two values for the magnetic dipole moment in the first and second rows in Table 2.The second value, µ * = 2.3 +0.4  −0.3 × 10 30 G cm 3 , corresponds to the inferred magnetic dipole moment including the N det /N correction for outbursts (D'Angelo 2017).

Impact of model assumptions
Systematic model uncertainties affect the inferred values of µ and η 0 in complicated ways, which are hard to quantify with the data at hand.Yatabe et al. (2018) discussed systematic uncertainties thoroughly in Section 4 of the latter reference, with an emphasis on the canonical picture of magnetocentrifugal accretion.For example, disk warping and precession induced by misalignment of the magnetic and rotation axes (Foucart & Lai 2011;Romanova et al. 2021) cause fluctuations in the mass accretion rate Q(t) (Romanova et al. 2021).If Q(t) and η(t) are anticorrelated, spikes in Q(t) produce dips in η(t), affecting the inferred value of η 0 .Such effects have been investigated recently in Figures 6 and 11 in Romanova et al. (2021).Misalignment of the rotation and magnetic axes leads to complicated three-dimensional accretion flows (Romanova et al. 2021;Melatos et al. 2023), requiring sophisticated, threedimensional magnetohydrodynamic simulations, which lie outside the scope of the current analysis.Other η 0 and µ biases neglected in this paper arise due to uncertainties in the location of the magnetospheric radius R m , and X-ray beaming, among other factors.The accretion dynamics in Section 2.2 are presented as a compromise between realism and the explanatory power of the limited volume of data that is currently available, i.e. ∼ 10 3 samples.Disentangling these and other complicated uncertainties will require larger datasets from future X-ray timing experiments.We emphasize that the Kalman filter approach complements (and does not supplant) other approaches to measuring µ (Klus et al. 2014;Ho et al. 2014;Shi et al. 2015;Christodoulou et al. 2016), which also neglect some of the above systematics but are likely to depend on them in different ways.Applying multiple approaches helps to assemble a complete picture.
In the current single-object analysis, it is difficult to ascertain what properties of SXP 18.3 dictate the inferred value η 0 = 0.04 +0.02 −0.01 .For example, one may ask: does the information that determines η 0 come mainly from P 1 (t) or L 1 (t) fluctuations?In rough terms, the answer is both: P 1 (t) implies higher Q 1 (t) than would be inferred assuming η 0 = 1 traditionally, implying η 0 ≪ 1 to match L 1 (t) as observed.
To get a preliminary sense of how P 1 (t) and L 1 (t) interact to produce the inferred η 0 value, we boost L 1 (t) artificially to test how η 0 changes.The results of the test are presented in Figure 4, where we plot η 0 (vertical axis), inferred using a different set of L 1 (t) fluctuations artificially enhanced by a factor L, versus L (horizontal axis).The numerical experiment follows the same procedure outlined in Section 4.1.It reveals that η 0 increases, as L increases, i.e. higher L 1 (t) fluctuations yield higher η 0 , with η 0 ∝ L1/2 for L ≳ 1.The opposite applies to P 1 (t): η 0 decreases, as P increases (not plotted for brevity).A more elaborate physical understanding in terms of the "moving parts" of magnetocentrifugal accretion is likely to emerge upon extending the present analysis to a larger population of X-ray pulsars.In particular, instantaneous (anti)correlations between the state variables and measurements, such as those presented in Section 4.2 for SXP 18.3, are likely to shed light on physical causes and their effects at the disk-magnetosphere boundary once measured at the population level.We defer this substantial project to a forthcoming paper.
One possible astrophysical interpretation of η 0 = 0.04 +0.02 −0.01 (out of many) is that disk accretion, as modelled in Section 2.2, is supplemented by outflows, as revealed by three-dimensional magnetohydrodynamic simulations (Romanova & Owocki 2015), e.g.slingshot-driven by the corotating magnetosphere in the regime R c ≲ R m .Some qualitative features of outflows are captured fortuitously by the mean-reverting disk dynamics in Equations ( 7) and ( 8); see the fourth and fifth panels of Figure 1 for timeresolved histories of Q(t) and S(t), revealing episodes of quiescence punctuated by accretion outbursts, as well as Appendix C.4 in Melatos et al. (2023) for an abridged review of outflow-launching mechanisms.It is plausible, in this context, that outflows are stronger (and hence η 0 is lower) when Q 0 is higher, as pressure builds up near the disk-magnetosphere boundary.It will be interesting to test for the existence of an inverse relation between η 0 and Q 0 in a future Kalman filter analysis of the SMC HMXB population.

Accuracy
Without independent measurements of the true hidden state sequence X n = [Ω(t n ), Q(t n ), S(t n ), η(t n )], it is challenging to verify how accurately the Kalman filter state estimates Xn reproduce X n in an astrophysical setting, e.g.SXP 18.3.Melatos et al. (2023) conducted Monte Carlo validation tests on the linear problem defined by Equations ( 13)-( 16) in Melatos et al. (2023), using synthetic time series with N = 500, i.e. comparable in volume to the SXP 18.3 time series.The tests return an error in the peak of the γ Ω posterior of ≈ 0.15 dex relative to the injected value, comparable to the uncertainty in the posterior for γ Ω in Figure 2.
The dispersion introduced by randomness in the noise realizations also limits the accuracy of the parameter estimates, because one cannot know where the actual, astrophysical noise realization lies within the ensemble of possibilities.Melatos et al. (2023) found that the full width half maximum (FWHM) of the γ Ω posterior in validation tests with synthetic data equals ≈ 0.60 dex, which agrees approximately with the FWHM in Figure 2 (≈ 0.31 dex).Monte Carlo validation tests also yield FWHMs of the posteriors of γ Q , γ S , σ QQ , and σ SS in line with the results in Figure 2, e.g. the FWHM of the γ Q posterior is ≈ 0.60 dex and ≈ 0.24 dex in the validation tests in Melatos et al. (2023) and Figure 2, respectively.
The analysis yields new estimates of the mean-reversion time-scales γ Q and γ S , as well the rms noise amplitude of the accretion rate.We find γ Q = 2.1 +0.5 −0.4 × 10 −7 s −1 , γ S = 2.5 +0.9 −0.7 × 10 −7 s −1 , and σ QQ = 8.0 +0.8 −0.8 × 10 13 g s −3/2 .Independent estimates of γ A and σ AA with A ∈ {Q, S} are not available at the population level, e.g. for all SMC X-ray pulsars, because traditional analysis methods based on time-averaged data do not estimate γ A and σ AA , as discussed in Section 1.However, the results in Figure 2 should be regarded as reasonable for two qualitative reasons.First, the characteristic timescale of mean reversion in the bottom two panels in Figure 1 and the second and third columns in Figure 2 is γ −1 A ∼ 10 7 s, with A ∈ {Q, S}.This is consistent with the Lomb-Scargle PSD computed from the light curves of other accreting systems such as Scorpius X−1 and 2S 1417−624, which roll over at ∼ 10 −7 s −1 (Mukherjee et al. 2018;Serim et al. 2022).Second, the inferred accretion rate fractional fluctuations satisfy 3 in Mukherjee et al. (2018).

CONCLUSION
Observations of accretion-powered pulsars in the SMC (Galache et al. 2008;Coe et al. 2015;Yang et al. 2017) and elsewhere reveal stochastic variability in the measured aperiodic X-ray flux L(t) and pulse period P (t) (Bildsten et al. 1997), driven by several complicated mechanisms such as hydromagnetic instabilities at the diskmagnetosphere boundary (Romanova et al. 2003;Kulkarni & Romanova 2008;Das et al. 2022) and flicker noise due to propagating fluctuations in the disk α parameter (Lyubarskii 1997).The time series L(t n ) and P (t n ) can be combined in several ways to shed light on the accretion physics at the disk-magnetosphere boundary and impor-    Table 2. Equilibrium quantities Q 0 (second column), S 0 (third column), and η 0 (fourth column), together with the magnetic dipole moment µ (fifth column).The first row corresponds to the linear Kalman filter analysis in Sections 2 and 4 above.The second row corresponds to Q * 0 and µ * , with a correction made for outbursts (D'Angelo 2017).The third row corresponds to using time-averaged data assuming Q 0 ∝ ⟨L(t)⟩ and η 0 = 1 (Klus et al. 2014).ple using the linear Kalman filter and state and measurement models presented in Melatos et al. (2023) and discussed in Sections 2 and 4 above.The latter subsample, however, poses a new challenge and necessitates three adjustments to the current framework.First, we replace Equations ( 12) and ( 13) with their nonlinear counterparts, Equations ( 1) and (2), respectively.Second, we replace Equation ( 10) with the nonlinear magnetocentrifugal torque law (5), and a system of modified Langevin equations; see Equations ( B1)-(B3) in Melatos et al. (2023).Third, we replace the linear Kalman filter with a nonlinear filter, e.g. an unscented or extended Kalman filter.Estimating µ, η 0 , and instantaneous (anti)correlations between time series at a population level may shed light on accretion physics in the SMC in general.We will explore this opportunity in a forthcoming paper.
where s n denotes the covariance of the innovation vector.In practice we work with log p({Y n } N n=1 |Θ), which is equivalent to Equation (15) in the main body of the text.The Kalman filter likelihood is combined with the prior distribution p (Θ) to estimate the posterior on the parameters Θ according to Bayes' rule ) where the denominator on the right-hand side of Equation ( A13) is known as the evidence and is defined according to The primary goal of nested sampling is to numerically approximate Equation (A14).Posterior samples are an important by-product.The specific steps of the algorithm are sketched out in Section 4.1.

B. ALTERNATIVE STATE-SPACE VARIABLES
It is possible in principle to employ a different set of state-space variables to replace those introduced in Section 2, i.e.X = [Ω(t), Q(t), S(t), η(t)].For example, one may elect to replace S(t) with R m (t) and Equation ( 8) with a mean-reverting Langevin equation for R m (t), viz.
where γ Rm and ξ Rm (t) are interpreted analogously to the damping constants and fluctuating driving terms in Equations ( 7) and ( 8).The Kalman filter analysis proceeds as in Section 2 with the following modifications: (i) S 0 is replaced with (2π) −1 µ 2 R −6 m0 ; and (ii) the linearized equations ( 10) are replaced with One subtle astrophysical implication of replacing Equation (8) with Equation (B15) is as follows.Equation (B15) implies that fluctuations in R m (t) are driven stochastically and independently from fluctuations in Q(t) and S(t), e.g.due to fluctuations as the companion star evolves or instabilities in the magnetic geometry at the diskmagnetosphere boundary, whereas Equation (3) implies that R m (t) is slaved to Q(t) and S(t).Hence the two approaches are different physically and lead to different results in general.It is impossible to say at present, theoretically or observationally, whether real HMXBs are approximated better by Equation ( 8) or (B15); even three-dimensional magnetohydrodynamic simulations are not yet equipped to offer a clear-cut answer (Romanova & Owocki 2015), and data volumes are too small to return a clear-cut preference from Bayesian model selection at present.
By selecting X = [Ω(t), Q(t), S(t), η(t)] as the state-space variables, as in this paper, it is possible to see how R m (t) and R c (t) evolve with time by substituting their timeresolved histories from the Kalman filter into Equations ( 3) and ( 4), respectively.One potential avenue for future work, once more data become available, is Bayesian model selection between the state-space framework presented in Section 2 and one involving alternative state-space variables; see Section 3 in Thrane & Talbot (2019) for an overview on Bayesian model selection in the context of gravitational wave astronomy.
We focus on the Maxwell stress in this paper because it is interesting astrophysically and notoriously difficult to measure.For example, the present analysis reveals important information about how S(t) evolves in time and its relationship with (for example) the measured pulse period, e.g.we find moderate evidence for an anticorrelation with r[S 1 (t), P 1 (t)] ≈ −0.36±0.032 in Section 4.2.Extending the present analysis to a larger population of SMC HMXBs has the potential to reveal important information about the accretion physics and S(t) in particular at the disk-magnetosphere boundary.

C. DYNAMICAL MODE OF TRADITIONAL GHOSH-LAMB ANALYSIS
The Kalman filter analysis presented in Section 4 is not the only way to extract time-dependent information from P (t) and L(t) fluctuations and hence infer µ.In this appendix, we discuss one alternative studied previously by other authors, namely the dynamical mode of traditional Ghosh-Lamb analysis (Yatabe et al. 2018).
The traditional technique to infer µ using the time-averaged observables ⟨P (t)⟩ = P 0 and ⟨L(t)⟩ = L 0 was generalized by Yatabe et al. (2018) to infer the magnetic field strength B and mass M of the HMXB X Persei.The latter authors numerically estimated a dP/dt time series, calculated from dP/dt = [P (t n+1 ) − P (t n )]/(t n+1 − t n ), and inferred the underlying, static parameters by χ 2 minimization.The technique, when applied to Equation (1) in Yatabe et al. (2018), yields a range of B and M estimates which return the same χ 2 minimum, with 5.0 ≤ B/(10 13 G) ≤ 23.
The dynamical Ghosh-Lamb analysis performed by Yatabe et al. (2018) complements time-averaged techniques and Kalman filter analyses.In this appendix, for the sake of completeness and by way of comparison, we apply the dynamical Ghosh-Lamb analysis to infer µ for the X-ray transient SXP 18.3.In what follows we assume fixed canonical values for the mass M and radius R of the star and generate the dP/dt time series numerically from P (t n ) measurements in the top panel of Figure 1.We encourage the reader to consult Equation (1) and Section 4 in Yatabe et al. (2018) for fuller details of the analysis.
In Figure 5 we graph χ 2 as a function of µ.The minimum χ 2 occurs at µ ≈ 4.4 × 10 30 G cm 3 , to be compared with µ = 8.0 +1.3  −1.2 × 10 30 G cm 3 inferred in Section 4.3 with a Kalman filter, and µ = 5.0 +1.0 −1.0 × 10 30 G cm 3 inferred by Klus et al. (2014) using time-averaged data and assuming η 0 = 1.It is encouraging that the three µ estimates are comparable for two reasons.(i) SXP 18.3 is classified in Table 3 in Yang et al. (2017) as being near magnetocentrifugal equilibrium, with a pulse period standard deviation satisfying σ P = 0.025 s ≪ P 0 (i.e.linear regime).(ii) The three estimates are derived from similar, Ghosh-Lamb-type models of magnetocentrifugal disk accretion.Klus et al. (2014) set dP/dt = 0 to infer µ; likewise, point (i) above implies dP/dt ≈ 0, when minimizing χ 2 in Figure 5.We do not expect the output from both techniques to coincide in general.For example, the Kalman filter features temporal autocorrelations between successive samples Q(t n ), Q(t n+1 ), . . .[and indeed S(t n ), S(t n+1 ), . ..], caused by the factors γ Q and γ S , whereas χ 2 minimization does not; it treats every sample independently.Although detailed comparisons between χ 2 minimization and the Kalman filter are beyond the scope of the current paper, we refer the interested reader to Sections 2-4 in Zarchan (2005) for such details.In the future, it will be interesting to compare the output from the nonlinear Kalman filter analysis in Appendix B in Melatos et al. (2023) with the output from χ 2 minimization and the dynamical Ghosh-Lamb analysis for more of the SMC HMXBs classified by Yang et al. (2017).

N
live randomly from the prior distribution p(Θ).For each k, and for all 1 ≤ m ≤ N live , we calculate L i.e. we run the Kalman filter while fixing Θ (k) m and evaluate Eq. (15).The sampler replaces the live point whose likelihood is lowest, namely Θ Figure1.Kalman filter state tracking applied to SXP 18.3.Inputs: observational data of fractional pulse period fluctuations P 1 (t) (first panel) and aperiodic X-ray luminosity fluctuations L 1 (t) (second panel).The RXTE data are plotted as grey points, with significant pulsations overplotted as green points (first and second panels).Outputs: Kalman filter estimates P1 (t) (top panel, blue curve) and L1 (t) (second panel, orange curve), and the absolute (not fractional) state variables Ω (t) (third panel, blue curve), Q (t) (fourth panel, cyan curve) and S (t) (fifth panel, green curve).In the bottom three panels, the state vector equilibrium parameters are plotted as horizontal, dashed lines.The shaded regions correspond to the 1σ state estimate uncertainties, returned by the Kalman filter.Ω 0 and Q 0 are inferred from Equations (A1) and (A4) in Melatos et al. (2023) and S 0 is inferred from Equation (14) above, discussed in Section 2.4.The time units on the horizontal axis are MJD.

Figure 2 .
Figure 2. Corner plot of the posterior distribution of the five fundamental model parameters Θ = (γ Ω , γ A , σ AA ), with A ∈ {Q, S}, for SXP 18.3.All distributions are plotted on a log 10 scale.Contour plots depict the posterior marginalized over three out of five parameters.Histograms depict the posterior marginalized over four out of five parameters.We estimate four parameters unambiguously, namely γ Ω , γ Q , γ S , and σ QQ .In contrast, σ SS rails against the upper bound of the prior and cannot be recovered.

Figure 3 .
Figure 3. Corner plot of the posterior distribution of the four induced model parameters (Q 0 , S 0 , η 0 , µ) for SXP 18.3.All distributions are plotted on a log 10 scale.Contour plots depict the posterior marginalized over two out of four parameters.Histograms depict the posterior marginalized over three out of four parameters.

Figure 4 .
Figure 4. Radiative efficiency η 0 inferred using a different set of L 1 (t) fluctuations artificially enhanced by a factor L, versus L. The plotted points correspond to the η 0 onedimensional posterior median values.The horizontal axis is plotted on a log 10 scale.

Table 1 .
Yang et al. (2017)3)), as well asVanderPlas (2018)for how to detect periodic signals with Lomb-Scargle periodograms.In X-ray timing analysis, it is standard practice to preference observations with statistical significance ≥ 99%.SXP 18.3 was discovered byCorbet et al. (2003), who detected pulsations with periods 18.37 ± 0.01 s, 18.36 ± 0.01 s and 18.37 ± 0.1 s on 2003 November 26, 2003 Candidate SMC HMXBs for Kalman filter analysis.The pulse period derivative (column 2), proximity to spin equilibrium (column 3), pulse period standard deviation (column 4), total number of RXTE PCA observations (column 5), and number of significant detections (column 6) are denoted by Ṗ , ϵ, σ P , N , and N det respectively.The subsample in Table1is defined by |ϵ| ≤ 1.5 and N ≥ 850.December 4, and 2003 December 11 respectively, using the RXTE PCA.SXP 18.3 is categorised in Table3inYang et al. (2017)as being near magnetocentrifugal equilibrium; a linear fit of P (t n ) versus t n yields dP/dt ≈ (3 ± 4) × 10 −6 s day −1 and ϵ = 0.75.Klus et al. ( Kalman filter state tracking applied to SXP 18.3.Inputs: observational data of fractional pulse period fluctuations P 1 (t) (first panel) and aperiodic X-ray luminosity fluctuations L 1 (t) (second panel).The RXTE data are plotted as grey points, with significant pulsations overplotted as green points (first and second panels).Outputs: Kalman filter estimates P1 (t) (top panel, blue curve) and L1 (t) (second panel, orange curve), and the absolute (not fractional) state variables Ω (t) (third panel, blue curve), Q (t) (fourth panel, cyan curve) and S (t) (fifth panel, green curve).In the bottom three panels, the state vector equilibrium parameters are plotted as horizontal, dashed lines.The shaded regions correspond to the 1σ state estimate uncertainties, returned by the Kalman filter.Ω 0 and Q 0 are inferred from Equations (A1) and (A4) in Melatos et al. (2023) and S 0 is inferred from Equation (14) above, discussed in Section 2.4.The time units on the horizontal axis are MJD.