A Bayesian approach for torque modelling of BeXRB pulsars with application to super-Eddington accretors

In this study we present a method to estimate posterior distributions for standard accretion torque model parameters and binary orbital parameters for X-ray binaries using a nested sampling algorithm for Bayesian Parameter Estimation. We study the spin evolution of two Be X-ray binary systems in the Magellanic Clouds, RX J0520.5-6932 and RX J0209-7427, during major outbursts, in which they surpassed the Eddington-limit. Moreover, we apply our method to the recently discovered Swift J0243.6+6124; the only known Galactic pulsating ultra-luminous X-ray source. This is an excellent candidate for studying the disc evolution at super-Eddington accretion rates, for its luminosity span several orders of magnitude during its outburst, with a maximum $L_{\rm X}$ that exceeded the Eddington limit by a factor of $\sim 10$. Our method, when applied to RX J0520.5-6932 and RX J0209-7427, is able to identify the more favourable torque model for each system, while yielding meaningful ranges for the NS and orbital parameters. Our analysis for Swift J0243.6+6124 illustrates that, contrary to the standard torque model predictions, the magnetospheric radius and the Alfv\'en radius are not proportional to each other when surpassing the Eddington limit. Reported distance estimates of this source range between 5 and 7 kpc. Smaller distances require non-typical neutron star properties (i.e. mass and radius) and possibly lower radiative efficiency of the accretion column.


INTRODUCTION
X-ray pulsars (XRPs) are astronomical objects powered by accretion that display periodic variations in X-ray intensity. They are formed when highly magnetized (B > 10 11 G) neutron stars (NSs) are found in close binary systems, allowing material to be transferred by the donor star onto the NS surface. The spin period evolution of a NS is indicative of the type of its accretion mechanism, the accretion disc structure and the neutron star magnetic field B. XRPs can also be useful in developing our understanding of the evolution process of binary systems with neutron-star members (Bildsten et al. 1997).
The majority of XRPs are found in Be X-ray binaries (BeXRBs) (see Reig 2011, for a review on BeXRBs). In this case material escapes the massive donor through a slow equatorial outflow, which is usually known as the decretion disc or the Be disc (e.g. Krtička et al. 2011). The mechanism behind the formation and depletion of the Be disc is still a matter of debate, however its transient nature results in highly variable mass transfer, and causes outbursts in BeXRBs. Observations of BeXRBs point to transient activity that is manifested in the form of two types of outbursts (e.g. Stella et al. 1986; Bildsten E-mail: georgios.vasilopoulos@astro.unistra.fr et al. 1997). Type I outbursts (L X ∼10 36−37 erg s −1 ) may occur during a close passage of the NS to the decretion disc and thus show a correlation with the binary orbital period. Giant or type II outbursts (L X ≥10 38 erg s −1 ) that can last for several orbits are associated with warped Be-discs (Okazaki et al. 2013).
During outbursts an accretion disc is formed around the NS, resulting in angular momentum transfer to the NS and a change of its spin. At a zeroth order approximation the NS spin changes due to mass accretion. However, the consensus is that torque acts through what is called a "magnetically threaded disc model" (first introduced by Ghosh & Lamb 1979) that describes the coupling of the NS magnetic field lines and the accretion disk resulting in torques acting on the NS (see overview by Parfrey et al. 2016).
Studies of the brightest type II outbursts (L X > 10 38 erg s −1 ) became especially relevant in the advent of the recent discoveries of pulsating ultra-luminous X-ray sources (PULXs, e.g. Bachetti et al. 2014;Israel et al. 2017;Carpano et al. 2018). ULXs are extragalactic point sources with an apparent isotropic luminosity above the Eddington limit for a 20 M black hole (see Kaaret et al. 2017, for a recent review). The discovery and subsequent study of PULXs confirmed that at least some ULXs are powered by highly magnetised NSs. Indeed, an increasing number of authors have put forward the hypothesis that a large fraction of ULXs may actually be powered by strongly magnetized NSs (see e.g. Koliopanos et al. 2017;King et al. 2017;Walton et al. 2018), building upon the early models for XRPs (Basko & Sunyaev 1976) but also on more recent theoretical considerations (Mushtukov et al. 2015a).
Estimates of the NS magnetic field can be made directly through the detection of cyclotron emission lines in XRP spectra (Staubert et al. 2019). These lines may be directly formed in the accretion column (e.g. Basko & Sunyaev 1976;Mushtukov et al. 2015b) or through reflection onto the NS atmosphere (e.g. Poutanen et al. 2013;Kylafis et al. 2021). However, these direct measurements are hindered by the spectral resolution and energy range of our instruments that make difficult the detection of lines corresponding to magnetic field strengths B 10 13 G. Note however that INTEGRAL has pushed this limit for nearby bright systems (Winkler et al. 2003). Alternatively, indirect measurements of B may be derived from the spin evolution of the NS during major outbursts. Such calculations require the use of torque models and proper corrections for the orbital motion of the binary. This method may be applied to systems with a wide range of magnetic field strengths, including PULXs.
Given that PULXs host magnetized NSs several authors have invoked standard torque models (i.e. Ghosh & Lamb 1979;Wang 1995) to estimate the magnetic field of the NS (e.g. Vasilopoulos et al. 2018Vasilopoulos et al. , 2019Vasilopoulos et al. , 2020aBachetti et al. 2020;Erkut et al. 2020;Chen et al. 2021). At the same time theoretical studies have demonstrated that it is required to adjust these standard torque models due to change in the disc structure when exceeding the Eddington limit (e.g. Bozzo et al. 2009). Moreover, according to Chashkina et al. (2017Chashkina et al. ( , 2019, and their numerical calculations, the radius of the magnetosphere should not be regarded as being to scale with the Alfvén radius for all mass accretion rates as suggested by the standard models (Ghosh et al. 1977;Koenigl 1991;Wang 1996;Kluźniak & Rappaport 2007). Instead, the ratio of the magnetospheric and Alfvén radii was found to depend on the mass accretion rate in a way that leads to an almost constant magnetospheric radius for super-Eddington mass accretion rates (see also Mushtukov et al. 2019).
In this work we study the spin evolution of accreting NSs during major outbursts of BeXRBs that reached or exceeded the Eddington limit using torque models that are widely used in the literature (i.e. Ghosh & Lamb 1979;Wang 1995;Ho et al. 2014). We also implement a nested sampling algorithm for Bayesian parameter estimation and apply it to our sample of sources to simultaneously estimate posterior distributions for the parameters of standard accretion torque models and binary orbital parameters.
We first test our approach against RX J0520.5-6932 (RX J0520 hereafter) and data obtained during a major outburst in 2014 that lasted for several orbits. Then we apply our method to two of the most energetic systems monitored by the Fermi Gamma-ray Burst Monitor (GBM, Meegan et al. 2009), namely RX J0209-7427 (RX J0209 hereafter) and Swift J0243.6+6124 (J0243 hereafter). For RX J0520 and RX J0209 we found that our method converges to a solution with almost no fine-tuning of the parameter space. In addition it provides more realistic uncertainties to the model parameters than typical methods based on least square fitting, and also enables investigation of degeneracies between parameters. The challenge was the modelling of J0243, a system with data that cover a large dynamic range in luminosity, and with maximum luminosity exceeding the Eddington limit by a factor of ∼ 10 considering a distance of ∼ 5-7 kpc (Reig et al. 2020;Doroshenko et al. 2018). The big variation of the bolometric luminosity of J0243 provides us with an excellent test case to examine the relation of the magnetospheric radius with the accretion rate. Standard accretion models for J0243 should be modified to account for the change in the magnetospheric radius at super-Eddington accretion rates, as demonstrated in recent theoretical and observational studies (e.g. Chashkina et al. 2017Chashkina et al. , 2019Mushtukov et al. 2019;Mönkkönen et al. 2019;Doroshenko et al. 2020). For this purpose, we find a parametric expression of the coefficient ξ, which is defined as the ratio of the magnetospheric radius to the Alfvén radius, as a function of the accretion rate. In other words, we move beyond the assumption of a constant ξ, usually made in the study of accreting pulsars. Our empirical approach would be applied for the first time in observational data of systems above the Eddington limit, but we refer the reader to Bozzo et al. (2009) for a parametric study of the disc-magnetosphere interaction models in lower luminosity accreting systems. This paper is structured as follows. In Secs. 2 and 3 we outline respectively the torque models and the observational data that will be used in our study. In Sec. 4 we present our methodology for modelling the spin evolution of the XRPs in our sample, and describe the Bayesian approach we implemented for the latter. In Sec. 5 we introduce the three systems in our sample and present the results of our analysis for each source in Sec. 6. We continue with a discussion of our results in Sec. 7 and finish with our conclusions in Sec. 8.

ACCRETION TORQUE MODELS
The problem of mass and torque transfer in accreting NS has been investigated by several studies in the past 50 years (e.g. see Frank et al. 2002;Parfrey et al. 2016, and references within). In the following paragraphs we will introduce the basic equations that we invoked in our work.
Assuming spherical accretion, the gas will stop at the so-called Alfvén radius, which is estimated by equating the magnetic pressure from the stellar dipole to the ram pressure of gas free-falling from infinity (Elsner & Lamb 1977;Davidson & Ostriker 1973): where M is the NS mass, µ = BR 3 is the magnetic dipole moment, with R the NS radius and B the NS magnetic field strength at the equator 1 ,Ṁ is the accretion rate, and G is the gravitational constant. We define the truncation radius of a thin Keplerian disc as the magnetospheric radius where ξ ∼ 0.5 − 1 for all kinds of magnetic stars (see Campana et al. 2018). After material gets halted at R m it may continue flowing towards the NS if its angular momentum is high enough to penetrate the centrifugal barrier set by the rotating magnetosphere. The radius where a particle attached to a field line would rotate at the Keplerian rate is defined as the corotation radius and is expressed as: where Ω is the NS angular velocity. Since matter inside the corotation radius flows along the field lines, for steady accretion to occur, the Keplerian angular velocity at R m has to be larger than the angular velocity of the star (and of the field lines). Following this rational Elsner & Lamb (1977) defined the fastness parameter as: where Ω K (r) is the Keplerian angular velocity at distance r.
A major consequence of accretion and the general interaction of the disc with the NS through its field lines is that the NS spin can change as a result of the induced torques (e.g. Ghosh & Lamb 1979). On the one hand, there is the torque applied to the star by the accretion, N acc , defined as On the other hand, there is a torque, N field that tends to spin-down the pulsar, and is applied by the dragging of the field lines by the disc and the sweeping of the open field lines due to the effective inertia of the electromagnetic field (see e.g. Bogovalov 1997). The total torque N tot is the sum of the two terms and it is usually expressed as where n(ω fast ) is a function of the dimensionless fastness parameter that incorporates the details of N field (Parfrey et al. 2016).
In the literature several torque models have been developed to explain the coupling of the disc with the magnetosphere and to estimate the induced torque onto the NS (e.g. Ghosh & Lamb 1979;Wang 1995;Kluźniak & Rappaport 2007;Lovelace et al. 1995;Rappaport et al. 2004). In our study we will focus on the Ghosh & Lamb (1979) model (hereafter GL79) and the Wang (1995) model (hereafter W95), as they are the most commonly used in the literature for accreting pulsars.
These models, expressed as seen in Eqs. (6)-(8) are only applicable when ω fast < 1. When we study bright systems during outbursts, we tend to ignore the ω fast terms in these expressions, as generally ω fast 1; in other words, the systems are away from equilibrium -a state in which the NS rotation frequency is constant. However, during the low luminosity phases the assumption of ω fast 1 might not hold and the full version of Eqs. (7) and (8) should be used for the treatment of the torque evolution.
When we have transitions from ω fast < 1 to ω fast > 1 (i.e. the system goes through equilibrium) the GL79 and W95 models can no longer be applied. J0243 is such an example. To model the spin evolution in such systems we can use an approximate expression for the total torque that reads Even though several studies have used the approximation of Eq. (9) (e.g. Menou et al. 1999;Wang & Tong 2020), Ho et al. (2014) have perhaps presented the first extended application to accreting XRPs for the study of their equilibrium state. In what follows, we therefore refer to Eq. (9) as the H14 model. Finally, the equation describing the spin up of the NS is given bẏ where I is the NS moment of inertia and N tot may be derived by Eqs. (5) and (7), (8) or (9). Using the above prescription one may indirectly estimate one of the fundamental parameters of the NS, its magnetic field strength. This is made possible because theν anḋ M can be inferred from observations, while parameters like the NS mass and radius are well determined.

OBSERVATIONAL DATA
Our methodology requires measurements of the spin period and mass accretion rates during major outbursts. It is crucial to obtain a baseline of measurements that would allow an estimation of the orbital parameters and the intrinsic spin-up due to accretion.
Outbursts are daily monitored in the X-rays by all-sky surveys like the Swift Burst Alert Telescope (Barthelmy et al. 2005) (BAT, 15-150 keV), the Fermi Gamma-ray Burst Monitor (Meegan et al. 2009) (GBM, and the Monitor of All-sky X-ray Image (Mihara et al. 2000) (MAXI, 0.5-30 keV). Moreover pointing observations may be performed by various observatories. In particular, the Neutron star Interior Composition Explorer (NICER) (Gendreau et al. 2016) and the Swift X-Ray Telescope (XRT) (Burrows et al. 2005) can perform multiple short observations (i.e. 1-2 ks) over weeks or months; thus, they are ideal for monitoring systems in the soft Xrays (i.e. 0.2-10 keV). Target of opportunity observations may also be performed by the Nuclear Spectroscopic Telescope Array NuS-TAR (Harrison et al. 2010) (3-79 keV) or AstroSat (Singh et al. 2014) (0.3-100 keV). These triggered observations last typically over 20 ks and are not repeated more than a couple times over the course of a major outburst. However, they are crucial as they deliver broadband spectra with high energy resolution and enable proper characterization of the spectral shape, the bolometric luminosity and the mass accretion rate.
In our study we will mainly use results that are available in the literature (through repositories), and perform limited data reduction of Swift/XRT data. For the latter case, we retrieved and analyzed the data from the UK Swift science data centre 2 using standard procedures as outlined in Evans et al. (2007Evans et al. ( , 2009.

Spin-period monitoring
While spin period measurements may be obtained from monitoring observations by NICER or Swift/XRT, they usually have larger uncertainties than the Fermi/GBM measurements. Hence, in this work, we will use Fermi/GBM data products from the GBM accreting pulsar project 3 to study the spin evolution of the NS (for details see Malacaria et al. 2020). These products contain spin measurements of data chunks that are typically binned every one to three days depending on the source luminosity.

Mass accretion rate estimation
As it was discussed in Sec. 2, measuring the mass accretion rate through monitoring observations is crucial for the torque modelling.
There are several ways one can deductṀ from observational proxies.
The Fermi/GBM products contain pulsed fluxes for epochs where a spin period could be obtained. Pulsed fluxes have been known to correlate with the luminosity of the pulsar, and evenν (e.g. for 2S 1417-624 Finger et al. 1996). However, pulsed fluxes are affected by changes in the pulse profile and the pulse shape. Nevertheless, in some systems, like RX J0209-7427, changes in the pulse profile are minimal; Vasilopoulos et al. (2020b) showed that the pulse shape remained almost constant during the evolution of its 2019 outburst. Thus, in certain cases the pulsed flux could be a good proxy of the accretion rate.
Alternatively one could use the Swift/BAT transient monitor results provided by the Swift/BAT team 4 that delivers daily binned data products in the form of count rates (Krimm et al. 2013). The advantage of BAT over GBM is that it provides intensities that are not tied to pulsation searches. Nevertheless, BAT count rates have larger uncertainties and more scatter (i.e. day-to-day) compared to the GBM detection of the same source. Thus, it is often required to bin the data over longer intervals, perform some smoothing of the overall light curve and even exclude outliers that often appear as flaring or dipping points 5 .
Upon selecting a proxy for the intensity, the next crucial step is the conversion to bolometric X-ray luminosity and finallyṀ. Ideally, for this purpose one should use broadband spectra (i.e. 0.3-70 keV). These can be obtained using a combination of instruments like NuSTAR and AstroSat for the hard band (∼1-100 keV) with XMM-Newton, Swift or NICER for the soft band (∼0.2-10 keV). The absorption-corrected flux can be estimated through spectral fitting, and can be converted to bolometric X-ray luminosity, L X . The latter can be then translated toṀ adopting some efficiency η eff under which gravitational energy is converted to radiation (typically assumed to be 100 per cent, Campana et al. 2018), namely L X = GṀM/R ≈ 0.2Ṁc 2 (M/1.4M )(R/10 km) −1 (Frank et al. 2002).

METHODOLOGY AND IMPLEMENTATION
In this section we will discuss the methodology we used when applying our model to the observed data of accreting pulsars. We will present the model parameters and the steps followed to construct the model. For the application to the data we will employ a Bayesian approach to ensure accurately estimated model parameters and their associated uncertainties.

Modelling the intrinsic spin-up
The first step is to estimate the intrinsic spin-up based on the selected torque model. The free parameters of the model are: • the magnetic field strength at the NS equator B, • the ratio of the magnetospheric radius to the Alfvén radius (i.e. the ξ parameter), • the NS spin frequency at some reference time (i.e. v 0 ) and, • the distance d to the source.
In general there is a degeneracy between ξ, B and d. For example, if the system is away from equilibrium with ω fast 1 there is a power-law dependence 6 of B on d (i.e. B ∝ d −6/7 ). This scaling can be easily understood as follows. For a larger distance d (and the same observed flux), the derived L X andṀ are higher, thus a lower B field is needed to explain the measuredν. However, a lower magnetic field means that the spin-equilibrium is going to be reached at lower fluxes. If the transition to spin-equilibrium is covered by monitoring observations, the degeneracy may be partially broken, since we now need to make assumptions for only one of the three free parameters ξ, B and d. The above discussion motivates studies of extragalactic XRBs where the distance of the host galaxies is well determined, such as the Magellanic Clouds. A similar power-law dependence exists between B and ξ. For thin accretion discs one may restrict ξ within a small range of values and in most cases it is safe to assume ξ = 0.5 (Ghosh et al. 1977). Nevertheless, we will also consider a mass-accretion dependent ξ parameter whenever relevant (see 4.2).
To calculate the intrinsic spin-up of the NS as a function of time, ν(t), for the three torque models described in Sec. 2, we used Eq. (10) and the inferredṀ(t) from one of the proxies described in Sec. 3.2. To minimize any sawtooth-like effects in the derived time series of the mass accretion rate, we re-sampled it with a finer resolution (i.e. 10-20 steps per day).
In all calculations we assumed a typical value for the NS moment of inertia, i.e. I = 1.3 × 10 45 g cm 2 unless stated otherwise, and considered that R co is constant in time. While the corotation radius can in principle evolve during an outburst, the expected change is very small given the minimal change in v and the large dynamical range of L X that drives theν evolution during an outburst (see Sec. 5).

A mass accretion dependent ξ
The value of the ξ parameter is motivated by theory of disc accretion (e.g. Ghosh et al. 1977). The inner region of a geometrically thin disc is gas-pressure dominated and ξ≈0.5 (Campana et al. 2018). However, as the accretion rate increases, radiation pressure becomes increasingly more important and the disc structure changes (becoming geometrically thick). Significant outflows from the disc play also an important role in this picture. It has been shown (e.g. Eq. 61 of Chashkina et al. 2017) that at a given mass accretion rate the magnetospheric radius becomes almost independent of the accretion rate, if the radiation pressure dominates at the inner parts of a disc (see also Mushtukov et al. 2019). In this case the ξ parameter gradually changes from ∼0.5 to ∼1.0 as a function of mass accretion rate. For even higher accretion rates the advection of viscously generated heat in the inner disc plays a more important role. Because of heat advection the radiation energy flux transported by diffusion in the vertical direction is less than the one released locally in the disc. As a result, the advection process effectively leads to a reduced mass loss from the disc (e.g. Mushtukov et al. 2019) and the relation between the magnetospheric radius and the accretion rate is closer to that of the standard models (i.e. R m ∼Ṁ −2/7 ) (see Fig. 12 of Chashkina et al. 2019). This change in disc structure is supported also by observational evidence in the power density spectra of pulsars (e.g. Mönkkönen et al. 2019;Doroshenko et al. 2020).
When modelling the spin evolution of XRPs a constant ξ is usually assumed, as it is rare to observe a transition through the above Figure 1. The magnetospheric radius given by Eq. (11) plotted against the dimensionless accretion rate for different choices of the free parameters a 1 , a 2 and a 3 (from left to right). Other parameters used are: a 0 = 3.459, corresponding to B = 10 13 G for M = 1.4M and R = 1.2 × 10 6 cm. Each a 1 value translates to a different ξ range and we expect ξ min 0.5 for a 1 0.14. The a 2 parameter indicates the characteristic value ofṁ where the outflows become significant and therefore the ξ parameter deviates from the standard ∼ 0.5 value. The a 3 parameter determines how wide is the area where ξ deviates from a constant value. mentioned accretion regimes during an outburst. Although extragalactic BeXRBs have been known to exceed the Eddington limit, it is difficult to obtain quality data at lower luminosity levels. Thus, J0243 offers a unique case-study with quality monitoring observations spanning over a large dynamical range around the Eddington limit. Therefore, we will implement an accretion-dependent ξ parameter in the modelling of this source, as described below.
Taking the logarithm of Eq.
(2) and using Eq. (1) we may write where a 0 is defined as In the equations above R g = GM/c 2 is the NS gravitational radius, andṁ is the accretion rate normalized to the Eddington accretion rate for a NS,Ṁ Edd . This is defined from the relation GMṀ Edd /R = L Edd = 4πcGM/κ, namelyṀ Edd = 4πcR/κ 1.3 × 10 18 (R/12 km) g s −1 , where κ ≈ 0.2(1 + X) cm 2 g −1 is the Thomson opacity and X = 0.7 is the hydrogen abundance. For typical NS parameter values, i.e. B = 10 12−13 G, R = 12 km, and M = 1.4M we find a 0 3 − 3.5. Motivated by the results of Chashkina et al. (2019) we developed a functional form for the ξ parameter, namely where a 1 , a 2 , and a 3 are parameters to be determined by the fit to the data (see Sec. 6). Here, a 1 describes the range of ξ values, a 2 describes the range ofṁ values where R m deviates from the standard scaling relation (i.e. R m ∝ṁ −2/7 ), and a 3 describes how fast ξ changes. To better illustrate the dependence of R m onṁ we plot Eq. (11) in Fig. 1 for different choices of the parameters a 1 , a 2 , and a 3 .

Modelling orbital spin evolution
The Doppler shifts induced by the orbital motion in an XRP are described with five orbital parameters: orbital period (P orb ), orbital eccentricity (hereafter e), the epoch of a mean longitude of 90 degrees of the star's orbit (T π/2 ), the semi-projected binary separation (a sin i) and the orbital phase that is commonly expressed as the angle of periapse (ω). Synthesizing NS radial velocities for a set of orbital parameters and at given times involves solving Kepler's equations, which can be done by an iterative method (e.g. Danby 1988;Fulton et al. 2018). Upon computing the radial velocities, one can combine it with the intrinsic spin evolution to derive a complete model for the evolution of the NS frequency in time.
where V r refers to the radial velocity of the NS due to the Keplerian orbit with parameters contained in the vector θ Kep .

Bayesian Inferenceultranest
Our goal is to infer the posterior probability density p given a dataset (D) and priors from the Bayes' Theorem for a model with a set of parameters contained in the vector θ. Having calculated the model we construct a likelihood function. Given the nature of the physical problem we added a term ln f to account for the systematic scatter and noise of our data not included in the statistical uncertainties of the measurements. This term results in an excess variance compared to statistical uncertainties, i.e.
where σ i are the GBM frequency errors, σ tot,i are the errors after accounting for the systematic scatter and noise not included in the statistical uncertainties of the measurements and i runs over the times of measurements. The likelihood function for a dataset D j can be then written as: where ν data are the measured spin frequencies. In principle, different datasets ( j = 1, · · · , N) can be combined to construct the total likelihood function of the model. To derive the posterior probability distributions and the Bayesian evidence we used the nested sampling Monte Carlo algorithm ML- .5 -58497.5 F, SB, Nu † Observatories whose data we used in this study: Fermi/GBM (F), NICER (N), NuSTAR (Nu), Swift/XRT (SX) and BAT (SB).
Friends (Skilling 2004;Buchner 2019) that employs the ultranest 7 package (Buchner 2021). The overall procedure is similar to methods used to derive Keplerian orbits from the time series of radial velocities (e.g. Fulton et al. 2018) that use Markov Chain Monte Carlo (MCMC) methods (Foreman-Mackey et al. 2013). The advantage of using ultranest lies in its overall strengths that are the unsupervised navigation of complex, potentially multi-modal posteriors until a well-defined termination point. Thus, no initial optimization is needed and minimal adjustment of the priors is necessary.

APPLICATION TO SYSTEMS
We present the systems that will be used as test beds of our methodology. We selected three BeXRBs that underwent outbursts exceeding the Eddington limit and are listed in Table 1. RX J0520 is a BeXRB located in the Large Magellanic Cloud (LMC) hosting a 8.04 s pulsating NS (i.e. LXP 8.04 Vasilopoulos et al. 2014a). In 2014 the system went through a major outburst that exceeded the Eddington limit (Vasilopoulos et al. 2014b;Tendulkar et al. 2014). The 2014 major outburst lasted for several months and was monitored by Fermi/GBM, Swift/XRT and Swift/BAT. Fermi/GBM monitoring resulted in determination of orbital parameters of the system (Malacaria et al. 2020). The major outburst was monitored by Fermi/GBM and Swift/BAT all sky detectors and by pointed Swift/XRT observations for more than seven orbital periods. Given that GBM detected pulsations for about 80 consecutive days, this makes the system an ideal test-case for our method.
RX J0209 is a BeXRB located in the outer wing of the Small Magellanic Cloud (SMC) hosting a 9.3 s pulsating NS (Vasilopoulos et al. 2020b). In November 2019 it exhibited a particularly bright outburst, among the brightest we have observed from a BeXRB in the Magellanic Clouds, reaching super-Eddington luminosity, that was detected by MAXI. During the outburst, the system was monitored by NICER, Fermi/GBM, AstroSat and Swift/BAT. Furthermore Fermi/GBM monitoring resulted in determination of preliminary orbital parameters of the system 8 .
Swift J0243 is the first and only known Galactic PULX (Wilson-Hodge et al. 2018). It was first detected by Swift/BAT on October 3, 2017 (Kennea et al. 2017) during an outburst that lasted until 2018. This source is characterized by a spin period of ∼9.86 s (Jenke & Wilson-Hodge 2017) and at its peak the L X is well above the Eddington limit (∼ 10 39 erg/s) (Doroshenko et al. 2018). During the period over which we have observations by both Fermi/GBM and Swift/BAT its luminosity varied over several orders of magnitude, making it an excellent test case for studying the evolution of the 7 https://johannesbuchner.github.io/UltraNest/ 8 GBM Accreting Pulsars project: https://gammaray.msfc.nasa.gov/gbm/science/pulsars.html magnetospheric radius with accretion rate. Upon the initial discovery of the system, the Gaia Data Release 2 estimated the distance of the NS as 6.8 +1.5 −1.1 kpc (Bailer-Jones et al. 2018), that was adopted by most follow up studies (e.g. van den Eijnden et al. 2018). However, analysis of data from the Gaia Data Release 3 (DR3) yielded a distance of 5.2 ± 0.3 kpc (Bailer-Jones et al. 2021). Finally, Reig et al. (2020) computed a distance of 5.5 ± 1.7 kpc based on BVRI photometric measurements to estimate the interstellar absorption. Assumptions about the distance on the source play an important role in the modelling of the NS spin evolution as we will see in the next sections.

RESULTS
We present the results of our analysis for the three systems introduced in the previous section.
6.1 RX J0520.5-6932 (LXP 8.04) To estimate the mass accretion rate as a function of time we used the standard methods described in Sec. 3. We first obtained the XRT and BAT count rates and pulsed flux from GBM. Then we used the bolometric X-ray luminosity obtained by NuSTAR data (i.e. L X = 4×10 38 erg s −1 at MJD 56682, Tendulkar et al. 2014) to scale GBM data. This resulted in the light curves shown in the upper panel of Fig.  2. The conversion factors we used are 6.55 × 10 40 erg s −1 count −1 s, 3.2 × 10 39 erg s −1 count −1 s and 4.2 × 10 37 erg s −1 count −1 s for BAT count rates, GBM pulsed fractions and XRT count rates respectively. The inferred light curves from all three instruments agree for the bright luminosity state. However, when the luminosity drops below about 2.5×10 38 erg s −1 , estimates based on GBM overshoot both XRT and BAT measurements. Another interesting feature is that BAT and XRT estimates are in good agreement between them when data from both instruments are available. Given that GBM pulsed fractions can be affected by changes in the pulse profile we opted to use the BAT data as proxy for the bolometric luminosity and the inferred mass accretion rate.
Having an estimate for the mass accretion we applied our recipe and fitted the data using the W95 and GL79 models. As an example, we show in Fig. 2 (b) the fitting result to the measured spin frequencies (including modulation because of the orbital motion) using the GL79 model. The evidence of the W95 model is ln Z =300.3 versus 302.9 for the GL79 (see Table 2), which translates to the latter being ∼ 10 times more probable than the W95 model, assuming the models are equally probable a priori. The posterior distributions of the orbital solution and the GL79 model parameters are presented in Fig. A1. The orbital parameters we recovered by both models are close to the estimates from previous works related to this system (see Table 2). Our orbital solution is also in agreement with the one presented in Malacaria et al. (2020) where the GBM pulse profiles phase offsets were modelled to refine the orbital solution of the source. Most importantly our method enables estimation of each model parameter and their uncertainties more accurately than the standard least-square minimization method (e.g. Sugizaki et al. 2017).
Based on the W95 model we estimated a polar magnetic field of ∼ 1.6 × 10 12 G for the NS. This is in agreement with other direct measurements of B. In particular, the study of the broadband spectrum of RX J0520 by NuSTAR also revealed the presence of a cyclotron resonance scattering feature at ∼ 31.5 keV yielding a direct measurement of B ∼ 2 × 10 12 G (Tendulkar et al. 2014).

RX J0209.6-7427
To estimate the mass accretion rate as a function of time we used the methods introduced in Sec. 4. First, we obtained NICER luminosity measurements and pulsed fluxes from GBM which we also used as a proxy for the L X . We modelled the luminosity using two methods, first using both the NICER and GBM data and then using only the GBM pulsed flux. As evident by comparing panels (c) and (d) of Fig. 3 the latter method yields better results. In fact, comparing the bolometric luminosities scaled using the GBM or NICER energy ranges alone (see Fig. 3 upper panel) we see that the NICER L X is systematically higher for the brightest phase of the outburst. This could be a result of a contribution to the NICER band from disc soft X-ray radiation, thus leading to an overestimation of the peak value of L X . A similar excess due to contribution from a soft component has also been reported in other super-Eddington accreting systems and has been proposed to be related to the hot accretion disc and/or outflows (Tao et al. 2019;Doroshenko et al. 2020).
Comparing the GL79 and W95 models, we find that the latter yields better results, with ln Z = 504.0 as compared to 492.3 for the GL79 model (see Table 2). Thus the W95 model is ∼ 2 × 10 5 times more probable than the GL79 model, assuming the models are equally probable a priori. The model parameters are listed in Table  2, while in Fig. A2 we show the corner plot of posterior distributions for the better model. The orbital parameters we recovered with the two models are similar to each other and to the estimates from previous works related to this system, with the exception of the semiprojected binary separation (a sin i), which is found to be smaller than the previously published values.

Swift J0243.6+6124
For the estimation of the system's bolometric X-ray flux, F X , we considered a linear relation with the Swift/BAT count-rates CR BAT , i.e. F X = A · CR BAT , where A is determined as follows. Tao   ing NuSTAR observations (see Table 3). Assuming that these are a good proxy of the bolometric flux, we performed a linear fit to those fluxes and the Swift/BAT count rates on the same days. The slope 9 of the linear fit was found to be 1.47 ± 0.13 × 10 −7 erg cm −2 s −1 (counts/s) −1 . Contrary to the other systems we studied, which lie in the Magellanic Clouds, the distance to J0243 is more uncertain despite the very accurate parallax measurements by Gaia (see Sec. 5 for more details). Therefore, we treated d as a free parameter, allowing it to take values between 4 kpc and 8 kpc. For the estimation of the mass accretion rate from the L X we used the same method as in the previous systems (for more details, see Sec. 3). Initially, we tried to fit our complete data set (MJD 58027.5-58497.5) following a similar procedure as for RX J0209 and RX J0520. This approach was hampered mainly by two issues: (i) the appearance of gaps in the GBM frequency monitoring due to the source entering a faint state; (ii) in the low flux states the source intensity is not properly characterised by Swift/BAT or other all-sky monitoring programs. The typical methodology used is to assume a steady spin-down term during these epochs (e.g. Sugizaki et al. 2017). However, given that this spin-down term is a function of accretion rate, we adopted a brute-force approach where we introduce a "jump" in frequency for every large gap in the available data. In the Fermi/GBM data we identified 6 large gaps. We therefore added 6 extra free parameters in our model. The first results obtained by  our MCMC approach (not shown here) revealed that a model with constant ξ for the whole duration of the outburst yields large residuals and cannot explain the data. Moreover, due to the extra free parameters, it required excessive computational time to run tests with different torque models and expressions of ξ(ṁ). For these reasons, we opted to model the orbital modulation and intrinsic torque separately. First, to derive the orbital parameters we used a chunk of observational data (MJD 58260 to 58460), where no large fluctuations of the flux were apparent, to calculate and remove the orbital effects from our problem, considering a constant ξ ≈ 0.5. The orbital parameters we recovered fall within 1σ from the results obtained from previous works on that system (see Table 4). Having removed the orbital modulation from our data, we proceeded with the modelling of the intrinsic NS spin up. Instead of estimating the spin-up rate us- ing a torque model and then fitting it to GBM observed frequencies (as we did in the previous systems), we calculated the derivative of the GBM frequency data and fitted the theoretical spin-up rate predicted by our model to the F X −ν space (we remind that distance is a free parameter). This method is more efficient, because multiple intermediate steps from our process can be eliminated from every run of our algorithm. Nevertheless, the results of the fit performed to ν(t) (similar to RX J0520 and RX J0209) or to the F X −ν space are consistent to each other.
We fitted the F X −ν data using the H14 torque model, since the other models cannot describe both the low and high L X regimes (for more details, see Sec. 2). We then considered two cases, one with ξ = 0.5 and another one with a physically motivated parametrization of ξ on the accretion rate, as described in Sec. 4 (see Eq. 13). The corresponding results are presented in Table 5 and in panels (b) to (d) of Fig. 4, where we compare the Fermi/GBM frequencies with the best-fit model and the residuals of both models for ξ. Inspection of panels (c) and (d) shows that a model with constant ξ yields larger residuals than the parametric ξ model of Eq. (13). The ln Z factor of the latter model was estimated to be 5506.5 as compared to 5417.7 for a constant ξ (see Table 5), meaning that the ξ(ṁ) model is ∼ 10 38 times more probable, assuming the both are equally probable a priori. In our model the magnetic field is calculated indirectly from the a 0 parameter, as described in Eq. (12). The result from the fit was log(B/G) = 13.061 ± 0.017, or a polar magnetic field of ∼ 2.3 × 10 13 G. Fig. 5 shows the best fit-solution for the constant ξ For clarity, the bolometric flux is plotted instead of the luminosity, since the distance to the system is a free parameter. Panel (b): Fitting results to the observational data (after removing the orbital effects) using the H14 torque model and the massaccretion dependent ξ parameter of Eq. (13). The orbital effects were removed by independently modelling the data obtained within the gray shaded region as described in the text. Panel (c): Residuals of the fit shown in panel (b). Panel (d): Residuals of the fit using the standard ξ =const. approach and the H14 torque model. and ξ(ṁ) models in the F X −ν space, with theν measurements of Fermi/GBM overplotted for comparison (see Fig. A3 for posterior distributions). We note that the variable ξ model can explain quite well the transition between the high and low luminosity regimes 10 .
At this point we should comment that the updated distance by Gaia (i.e. Gaia DR3, d = 5.2 ± 0.3 kpc) does not fall within the 3σ range of values derived by torque modelling (see Table 5). In fact, if we fix the distance at 5 or 6 kpc the data cannot be fitted by our model using standard NS parameters (see further discussion in Sec. 7.2). However, if we treat the NS mass, radius and moment of inertia as free parameters and set a hard limit on the distance at 6 kpc we are able to get acceptable fits to the data (see Table 5 and Figure  A4). 10 During the revision of this manuscript we became aware of an independent study (i.e. Liu et al. 2022a) that also noted a so-called flattening in thė ν − L X space of J0243. Same flattening effect would be evident in out Fig. 5 if plotted in linear scale. We note that our results quantitatively match their findings. Figure 5. Plot of the absolute value of the frequency derivative |ν| versus the X-ray flux F 3−79keV for J0243. Symbols indicate the observational values from Fermi/GBM, the solid magenta line shows the prediction of our bestfit model with ξ(ṁ) given by Eq. (13) and the dashed line shows the best-fit model using the standard approach with a constant ξ = 0.5.

Application to systems in the Magellanic Clouds
We have studied the properties of two BeXRB systems in the Magellanic Clouds using different torque models (i.e. Ghosh & Lamb 1979;Wang 1995), based on X-ray data collected during their outbursts. We used data from NICER, Fermi/GBM, XRT and Swift/BAT as a proxy for the luminosity and then scaled our time-series using NuSTAR observations to retrieve the bolometric L X . With that data we were able to get an expression for the theoretically predicted spin evolution of our sources, and to obtain posterior distributions for the magnetic field and the orbital parameters of our systems, using ultranest to fit to the Fermi/GBM frequency data.
For RX J0520 and RX J0209 we were able to simultaneously obtain orbital solutions for the binary and estimates on the B field of the NS. Most importantly, by using the full expression of the torque models, we are not limited to the asymptotic behaviour during ω fast 1, which is commonly used in the literature (e.g. Sugizaki et al. 2017;Weng et al. 2017), but we also explore behaviours where the slope changes in a log |ν| versus log L X diagram at lower accretion rates as we approach equilibrium. This approach enables to test which of the GL79 and W95 torque models can better explain the data. Thus we can favor one model over the other for RX J0209 and somewhat less prominently in the case of RX J0520 (see also Table 2). In Fig. 6 we plot |ν| versusṀ together with both models with their best fit parameters. However, the question about which torque model is favorable over the other is not entirely solved. There are other systematic uncertainties, like a luminosity dependent bolometric correction factor (e.g. Anastasopoulou et al. 2022), which could reduce the amount (i.e. ln(Z)) one model is favored over another.
The problem of accretion disc threading by stellar magnetic field still lacks a comprehensive solution as demonstrated by theoretical and observational studies (e.g. Bozzo et al. 2009;Filippova et al. 2017;Malacaria et al. 2020). Nevertheless, regardless of the torque model, it is possible to relax some of the underlying assumptions or   introduce extra terms to create wider or narrower cusp-like transitions around equilibrium. For example by assuming a misalignment between magnetic and rotation axes it is possible to induce a sharper transition near equilibrium (see Benli 2020, for an application to 4U 1626-67).

Implications from MC depth estimations
In the calculations concerning sources in the Magellanic Clouds a possible source of uncertainty is distance. While the average distances of the SMC and LMC are well determined, the depth of each galaxy is of the order of 4-8 kpc (Subramanian & Subramaniam 2009). Thus, it is prudent to at least explore if the fit to the data sets of the two systems can improve by treating the distance to each source as a free parameter. For simplicity we fitted the model to thė M −ν parameter space, although performing the fit on a similar manner as above yielded the same results and trends in the corner plots.
The results are shown in Fig. 7. For RX J0520 we found that the data favour a somewhat smaller distance than the average one of LMC, which is however still consistent with the depth of the galaxy. Interestingly, for RX J0209 regardless of the torque model used we find is a degeneracy between the distance and the magnetic field, meaning that we cannot put any constraints on the position of the system compared to the average distance of the nearby galaxies. The application to RX J0209 also demonstrates how the uncertainty in distance affects the B-field estimates. For example, using a uniform prior for the distance between 50-70 kpc, the 1 σ uncertainty in the derived magnetic field strength is about 3 or dex of 0.5.
Recently, an independent study of the spin-up of RX J0209 with the use of a generalised torque model (i.e.ν ∝ L a ) revealed a somehow steep dependence (i.e. a = 1) of spin-up rate on luminosity (Liu et al. 2022b). Considering that the authors used a fixed distance of 55 kpc, this would be consistent with our results.

What did we learn for the first Galactic PULX?
For J0243 we adapted the method used for the other two pulsars to tackle the complexity of the dataset. Given the dynamic range in the observed luminosity and the transitions between spin-up and spindown phases we opted for using a more generalised form for the H14 torque model.
One of the main difficulties arising when trying to create an empirical model for ξ(ṁ) is the lack of a way to attain direct measurements of the magnetospheric radius. As a result, we have to rely on comparing how well different methods describe our observational data, namely the luminosity and the NS spin. A way to overcome this difficulty is using a torque model and solving backwards for R m . This is possible using the H14 model, which can be successfully solved for R m without making any assumptions for the form of the ξ(ṁ) or the magnetic field, provided that we havev andṁ measurements.
We can create a set of observationalv andṁ values for J0243 by calculating the gradient of the GBM frequency data and assuming a distance of 6.95 kpc which is the result we got from our fit (see  Table A1) in order to be able to compare the data we generate to our model. Now, after interpolating the two sets of data to the same dates, using a linear interpolation method, we can solve the equation for R m . This equation has two solutions for the spin-up phase, which correspond to the two intersection points of a horizontal line (at a given value ofν) with a dashed curve shown in Fig. 8). The non acceptable solutions in the spin-up phase can be easily identified after plotting them on a R m −ṁ graph, since they fall away from the standard R m (ṁ) = ξR A (ṁ) (with ξ = 0.5 − 1) solutions by several orders of magnitude. Using this method we generated a set of R m data points derived from the observational data, making no assumptions for the dependency of the magnetospheric radius on the accretion rate or the magnetic field other than the ones inherent in the torque model we used (H14). The results of this method are portrayed in Fig. 9. The most intriguing result for J0243 is that we found evidence of an evolving disc in qualitative agreement with theoretical predictions (i.e. Chashkina et al. 2019). We can clearly see from our results that the magnetospheric radius is not to-scale with the Alfvén radius at high accretion rates, in contradiction with the standard torque models' predictions. Our model describes very well the magnetospheric radius evolution at super-Eddington accretion rates (see Fig.  9). However, we should stress the degeneracy between the range of ξ values and B in our approach (see Eq. 13). This degeneracy is the same to the constant ξ approximation that is evident if we leave both ξ and B free parameters (see Fig. A5). This introduces extra uncertainty in our estimation of B apart from the statistical uncertainty derived from the fit. For our estimations we opted for setting ξ max = 1 based on the upper limits on ξ reported in the literature. We could instead set a lower limit on ξ used for standard disc accretion, i.e. ξ min = 0.5. Fixing the lower ξ value to 0.5 yields log B(G) = 13.43 ± 0.08. Thus, the range of the accepted B values is (2 − 5) × 10 13 G. This estimate value is in agreement with the detection of a cyclotron resonance scattering feature between 120-146 keV in the insight-HMXT spectra (Kong et al. 2022).
Since we argue that our findings for an evolved ξ may be a result of changes in the disc, it is interesting to compare the transitions found here with independent studies. In particular, a sharp state transition in the spectral and temporal properties of the system has been reported based on insight-HMXT observations (Doroshenko et al. 2020). Based on the study of power-spectra and quasi periodic oscillations (QPOs) a transition in the pulse profile was found to occur at about two times the Eddington limit. This transition was proposed to mark the border between gas-pressure dominated and radiation-pressure dominated regions of the disc. Another interesting transition marks the change of the pulse profile of the pulsar from single peaked to double peaked (see Wilson-Hodge et al. 2018;Doroshenko et al. 2020). This critical transition has been attributed to the formation of the accretion column (see Becker et al. 2012) and has been used as a proxy for an indirect estimate of the magnetic field. For J0243 this transition was found atṁ ∼ 1 or L X ∼ L Edd for the distance of 7 kpc that we computed from our model (see Wilson-Hodge et al. 2018;Doroshenko et al. 2020). For comparison purposes with mark these transitions with vertical lines in Fig. 9.
As we mentioned earlier the updated Gaia distance of J0243 introduces difficulties in finding a torque model that can fit the observed Fermi/GBM data (see Sec. 6.3). This is because a smaller distance yields a lower maximum L X and lower mass accretion rate estimates. Thus, to explain the highest observed spin-up rates at the peak of the outburst a larger magnetic field strength is required (see Fig. 8). This increases further the magnetospheric radius and, in our case, pushes it very close to the corotation radius, prohibiting essentially any further spin up. A way around this problem was to let the NS parameters free. Indeed, a good fit was found for a NS with larger radius and smaller mass than the typically assumed values (see Table 5 and Fig. A4). Searching the literature for NSs with reliable mass estimates (Özel & Freire 2016), the double NS system J0453+1559 hosts the NS with the smallest measured mass of Figure 9. Plot of the magnetospheric radius versus the dimensionless accretion rate for J0243. Symbols indicate the values inferred by solving Eq. (17). The most probable model obtained by the Bayesian method is overplotted with a dashed magenta line. The vertical dash-dotted lines represent the change of the pulse profile of the pulsar from single peaked to double peaked and the transition from gas-pressure dominated to radiation-pressure dominated disc respectively from left to right (see Sec. 7.2).
1.174(4) M (Martinez et al. 2015). Moreover, a recent study of the isolated NS in the center of supernova remnant HESS J1731-347 (Doroshenko et al. 2022) indicated that the NS may be extremely light with having a mass of 0.77 +0.20 −0.16 M (1 σ errors). Thus, J0243 could potentially host a very low mass NS. However, this approach yielded a magnetic field strength smaller (by a factor of 2) than the one inferred by the reported cyclotron line (Kong et al. 2022). This would potentially mean that the cyclotron line is formed in regions with multi-polar magnetic field components (e.g. see evidence of such configuration Riley et al. 2019;Chen et al. 2020), compared to the torques that are associated to the dipole component. Alternatively, one needs to revise the assumptions of our model and in particular the radiative efficiency of the accretion column. More specifically, to reconcile the observed spin evolution of J0243 with the updated Gaia DR3 distance, the radiative efficiency should be lower by a factor of ∼ (7/5.2) 2 ∼ 1.8 (assuming standard NS parameters). State-of-the-art physical models about the emission of the accretion column generally assume that all gravitational energy is transformed to radiation (e.g. Wolff et al. 2016;West et al. 2017b,a). As new fitting strategies are implemented into these models (Thalhammer et al. 2021) one could further test the radiative efficiency in the accretion column in systems like J0243.

Further application
Our approach demonstrates that self-consistent modelling of the intrinsic and orbital spin-up of the system is essential for major outbursts. Coupling the torque models with a Bayesian interface delivers much more realistic uncertainties. Moreover, implementation of nested sampling may allow fitting data while using a wide parameter space for priors enabling better investigation of degeneracies and possible multimodal solutions. The Bayesian modelling is also useful for exploring orbital modulation in systems with low quality of data monitored with Swift/XRT or NICER, as seen in a recent application we made in SXP 15.6 (Vasilopoulos et al. 2022). In terms of the physical problem, inclusion of other torque models and extra terms would be the next step so the code can be applied to a wider family of accreting XRPs. Finally, we plan to build upon our current code, and provide a user friendly version to the community with parallelization capabilities.

CONCLUSION
We have used a nested sampling algorithm for Bayesian Parameter Estimation to study the spin evolution of nearby super-Eddington accreting pulsars. By coupling torque and orbital models for systems with well determined distance we were able to simultaneously estimate the orbital parameters and the magnetic field of the NS. A similar application to J0243, the closest known PULX, revealed a transition that may be quantitatively linked to changes in the accretion disc structure close to the Eddington luminosity. According to the most recent Gaia parallax measurements J0243 seems to be closer than previously thought. The study of the NS spin up demonstrates that typical NS parameters cannot be used to explain the NS spin evolution using the updated distance. A possible solution is to assume a low-mass NS (M ≈ 1.1M ) or assume a lower accretion column radiative efficiency (by a factor of 2) than typically assumed.

DATA AVAILABILITY
X-ray data are available through the High Energy Astrophysics Science Archive Research Center: heasarc.gsfc.nasa.gov. Swift/BAT data are available through Swift transient monitoring project: https://swift.gsfc.nasa.gov/results/ transients/weak/. Fermi/GBM data are available through the GBM Accreting Pulsars project: https://gammaray.msfc.nasa.gov/gbm/science/pulsars. html. Figure A1. Corner plot for RX J0520 using the GL79 model. We plot the logarithm of the magnetic field strength B in G. The eccentricity (e), orbital period in days, the longitude of periastron in degrees (ω) and the projected semi-major axis in light-sec (a sin i). For clarity T π/2 is given relative to a reference MJD of 56666.91, while the reference frequency (F 0 ) is given relative to 124.3927 mHz. Finally, ln( f ) is the systematic scatter that is used to estimate the excess variance of the model. Figure A2. Corner plot for RX J0209 using the GBM pulsed flux as a proxy of accretion rate and the W95 model. We plot the logarithm of the magnetic field strength B in G. The eccentricity (e), orbital period in days (P orb ), the longitude of periastron in degrees (Per) and the projected semi-major axis in light-sec (a sin i). For clarity, T π/2 is given relative to a reference MJD of 58793.32, while the reference frequency (F 0 ) is given relative to 107.4909 mHz. Finally, ln( f ) is the systematic scatter that is used to estimate the excess variance of the model. Figure A3. Corner plot for J0243 using the Swift/BAT count rates as a proxy of accretion rate, the H14 model, our ξ(ṁ) model and fitting to the frequency derivative instead of the frequency. Instead of a 0 we plot the magnetic field logarithm log(B) directly, in G. Finally, ln( f ) is the systematic scatter that is used to estimate the excess variance of the model. The upper right contour is calculated for the standard ξ = 0.5 approach. Figure A4. Corner plot for J0243 using the Swift/BAT count rates as a proxy of accretion rate, the H14 model, our ξ(ṁ) model, with the NS mass M, radius R and moment of inertia I ass free parameters and fitting to the frequency derivative instead of the frequency. Instead of a 0 we plot the magnetic field logarithm log(B) directly, in G. The mass M is given in M , the radius R in 10 6 cm and the moment of inertia I in 10 45 g cm 2 . Finally, ln( f ) is the systematic scatter that is used to estimate the excess variance of the model. Figure A5. Corner plot for J0243 using the Swift/BAT count rates as a proxy of accretion rate, the H14 model and the standard ξ =const. approach. We plot the logarithm of the magnetic field strength B in G. For clarity the reference frequencies (F 0 to F 6 ) are given relative to the result of the fit displayed in table A1. Finally, ln( f ) is the systematic scatter that is used to estimate the excess variance of the model. Table A1. Results of modelling J0243 in the ν(t) parameter space.