SAX J1810.8-2609: An Outbursting Neutron Star X-ray Binary with Persistent Spatially Coincident Radio Emission

Here we report on joint X-ray and radio monitoring of the neutron star low-mass X-ray binary SAX J1810.8-2609. Our monitoring covered the entirety of its ~5 month outburst in 2021, revealing a temporal correlation between its radio and X-ray luminosity and X-ray spectral properties consistent with a `hard-only' outburst. During the outburst, the best-fit radio position shows significant variability, suggesting emission from multiple locations on the sky. Furthermore, our 2023 follow-up observations revealed a persistent, unresolved, steep spectrum radio source ~2 years after SAX J1810.8-2609 returned to X-ray quiescence. We investigated potential origins of the persistent emission, which included an unrelated background source, long-lasting jet ejection(s), and SAX J1810 as a transitional millisecond pulsar. While the chance coincidence probability is low (<0.16%), an unrelated background source remains the most likely scenario. SAX J1810.8-2609 goes into outburst every ~5 years, so monitoring of the source during its next outburst at higher sensitivities and improved spatial resolutions (e.g., with the Karl G. Jansky Very Large Array or Square Kilometre Array) should be able to identify two components (if the persistent emission originates from a background source). If only one source is observed, this would be strong evidence that the persistent emission is local SAX J1810.8-2609, and future monitoring campaigns should focus on understanding the underlying physical mechanisms, as no neutron star X-ray binary has shown a persistent radio signal absent any simultaneous X-ray emission.


INTRODUCTION
Low-mass X-ray binaries (LMXBs) are interacting binary systems that consist of a compact object -a black hole or a neutron star -accreting material from a low-mass companion star (< 1  ⊙ ).The inward-moving accretion flow powers outflows in the form of disk winds and relativistic jets.Many LMXBs are transient systems, spending the majority of their lifetimes in a low-luminosity quiescent state (  ≲ 10 32 erg s −1 ) before sporadically entering into bright transient outbursts (  > 10 35 erg s −1 ) that last weeks to years (e.g., McClintock & Remillard 2006;van der Klis 2006).Since LMXBs rapidly evolve through multiple accretion states during outbursts, LMXBs act as natural laboratories for the study of accretion flows (best measured at X-ray frequencies; e.g., Belloni et al. 1999;Muñoz-Darias et al. 2014;Chakraborty et al. 2021) and relativistic jets (best measured at radio through infrared frequencies; e.g., Corbel & Fender 2002;Russell et al. 2015;Tetarenko et al. 2017).
The standard accretion state nomenclature (i.e., the hard and soft accretion states) was developed to describe the different X-ray spectra ★ E-mail: hughes1@ualberta.caobserved in black hole low-mass X-ray binaries (BHXBs).Moreover, the properties of the relativistic jet(s) are closely correlated with the accretion state (see, Fender et al. 2004;Remillard & McClintock 2006;Belloni 2010;Fender 2010, for detailed reviews).In the hard accretion state, the X-ray emission is dominated by high-energy (i.e., hard) X-ray photons comptonized by an optically thin corona.The X-ray spectra are well described by a power law model with a photon index of Γ∼ 1.7 (where the X-ray flux  X () ∝  −Γ−1 ).Furthermore, in the hard accretion state, the jet adopts a steady, compact structure.The radio spectrum of the compact jet is the result of a superposition of multiple self-absorbed synchrotron spectra originating from different positions along the jet axis (Blandford & Königl 1979).At low frequencies, the jet is best described as an optically thick, partially self-absorbed synchrotron spectrum with an inverted or flat spectral index ( ≳ 0; radio flux density   () ∝   ) up to a break frequency (often at sub-mm wavelengths).Beyond the break frequency, the jet's spectrum becomes optically thin ( ∼ − 0.7; Migliari et al. 2010;Russell et al. 2013;Díaz Trigo et al. 2018).
In the hard state, the X-ray (  ) and radio (  =  , ) luminosities are correlated (henceforth, the   -  relation; Gallo et al. 2003;Corbel et al. 2013).After including a scale for the black hole mass, the   -  relation has been extended to include accreting supermassive black holes (Merloni et al. 2003), thereby spanning 10 orders of magnitude in X-ray luminosity and providing the strongest empirical evidence of the coupling between accretion flows and relativistic jets.Individual BHXBs have exhibited multiple distinct tracks in the   −   plane (e.g., the 'radio-loud ' and 'radio-quiet' tracks;Coriat et al. 2011;Espinasse & Fender 2018;Williams et al. 2020;Carotenuto et al. 2021b) suggesting that the properties of the accretion flow (e.g., geometry and radiative efficiency) may vary significantly in the hard accretion state.Population analyses have both supported (e.g., Gallo et al. 2012) and refuted (e.g., Gallo et al. 2014Gallo et al. , 2018) ) the statistically independent existence of multiple tracks, with the more recent studies not finding any robust statistical evidence for separate tracks, suggesting that, instead, the properties of the 'radioloud' and 'radio-quiet' track sources vary significantly from source to source.
Conversely, in the soft accretion state, low-energy (i.e., soft) thermal emission from a multi-color accretion disk dominates the X-ray spectrum.Furthermore, the compact jet is quenched, decreasing in luminosity by ≳ 3 orders of magnitude (Coriat et al. 2011;Russell et al. 2020).During the hard-to-soft transition, one or more discretized ejection events may be launched.These ejections have been spatially resolved in multiple sources (e.g., Mirabel & Rodríguez 1994;Hjellming & Rupen 1995;Hannikainen et al. 2001;Rushton et al. 2017;Bright et al. 2020).The radio spectra of the ejecta are characterized by a time-variable self-absorbed synchrotron component (sometimes parameterized as the van der Laan (vdL) model; van der Laan 1966; Hjellming & Johnston 1988;Hjellming & Han 1995).As ejecta propagate and expand, they become optically thin at (progressively) lower-frequency emission, steepening the radio spectral index to  ∼ − 0.7.Emission from jet ejections can persist from hours to years (e.g., Miller-Jones et al. 2019;Bahramian et al. 2023), and can exhibit variability that is unrelated to any simultaneous evolution of the accretion flow (e.g., through collision with the surrounding interstellar medium; Carotenuto et al. 2021a).As a result, radio observations of jet ejecta must be excluded from the   -  relation.
For neutron star (low-mass) X-ray binaries (NSXBs), their strong intrinsic magnetic fields and solid surfaces complicate the picture.Historically, radio emission was thought to be exclusive to the weakly-magnetic (< 10 10 G) sub-population, although there have been recent detections of radio emission from strongly-magnetic NSXBs (e.g., van den Eĳnden et al. 2018;van den Eĳnden et al. 2021).The weakly-magnetic NSXBs are most directly analogous to BHXBs; thus, the strongly-magnetic sub-population will not be discussed any further (henceforth, NSXBs only refer to weaklymagnetic NSXBs).NSXBs have two main sub-classes; atoll and Z sources (named for their tracks in colour-colour diagrams, see van der Klis 2006, for a review).Atoll sources tend to be lower luminosity and transient, exhibiting similar hard/soft accretion states as transient BHXBs.In contrast, Z sources are often persistent but show rapid timescale variability.Moreover, although Z sources also transition through multiple accretion states, these states tend to be softer than atoll states (Muno et al. 2002).Some NSXBs have shown transitions from Z to atoll behaviour at lower X-ray luminosities (and thus accretion rates, e.g., XTE J1701-462; Homan et al. 2007;Lin et al. 2009), suggesting that the these may not be unique sub-populations, but instead that Z sources are NSXBs with the largest accretion rates (analogs to the rapidly flaring, semi-persistent BHXBs like GRS 1915+105; Migliari & Fender 2006).
Transient atoll sources more closely follow the evolution of a 'typical' transient BHXBs (see, Migliari & Fender 2006;Muñoz-Darias et al. 2014, for a review).Atoll outbursts exhibit distinct hard (also known as "extreme island") and soft (also known as "banana") accretion states.Atolls (sometimes) exhibit jet quenching in the soft state.State transition-induced jet ejections have been proposed for atolls, although they have only been observed in Z sources (e.g., Fomalont et al. 2001;Spencer et al. 2013).The 'typical' evolution of a transient outburst of an atoll NSXB or BHXB begins with a departure from quiescence through a rapid brightening in the hard state.The source transitions to the soft state following the initial brightening.The system then remains in the soft state for some time (the amount of time varies from system to system) until it begins to dim, eventually returning to the hard state at a lower Xray luminosity.Once back in the hard state, the system dims until it returns to a quiescent state.However, some systems break this paradigm by exhibiting erratic state transitions (e.g., Kajava et al. 2020) or failed (i.e., 'hard-only') outbursts (e.g., Rodriguez et al. 2006;Stiele & Kong 2016;Tarana et al. 2018;Stiele & Kong 2021).Recent analyses have shown that ∼ 40% of outbursts of BHXBs are thought to be 'hard-only' (Tetarenko et al. 2016); this fraction has not been thoroughly explored for NSXBs.
There are several significant differences between the neutron star and black hole X-ray binary sub-populations: (i) NSXBs generally have radio luminosities that are a factor of ∼ 20 lower than BHXBs at comparable X-ray luminosities (the discrepancy cannot be attributed to the difference in compact object mass; Gallo et al. 2018); (ii) NSXBs have shown compact jet radio emission in the soft accretion state (e.g., Migliari et al. 2004;Gusinskaia et al. 2017;van den Eijnden et al. 2021), suggesting the quenching process may not be as extreme as observed in black hole systems or possibly a different jet launching process completely; (iii) all accretion states can have an additional thermal X-ray component (often modeled as a black body component, see Lin et al. 2007) due to emission from the neutron star surface or boundary layer between the accretion disk and surface.Historically, studies of accretion-jet coupling of NSXBs have suffered from their weaker radio emission.Joint X-ray and radio monitoring of NSXBs is critical for understanding the differences between the neutron star and black hole X-ray binary populations and how the presence (or absence) of an event horizon, ergosphere, or solid surface affects the connection between the accretion flow and relativistic jet.In 2021, the NSXB SAX J1810.8−2609exhibited a multi-month outburst that was detected in both X-ray and radio frequencies, allowing for a comprehensive monitoring campaign.
1.1 SAX J1810.8−2609SAX J1810.8−2609(henceforth SAX J1810) is a NSXB that was initially discovered in 1998 by the wide-field X-ray cameras aboard the BeppoSAX satellite (Ubertini et al. 1998).Since its discovery, there have been four subsequent (detected) outbursts that occurred in 2007 (Degenaar et al. 2007), 2012 (Degenaar & Wĳnands 2013), 2018 (Negoro et al. 2018), and 2021(Iwakiri et al. 2021).A Type I X-ray burst (i.e., the runaway thermonuclear detonation of a hotdense surface layer of accreted matter, see Galloway & Keek 2021, for a review) revealed the presence of a solid surface, identifying the accreting object as a neutron star (Natalucci et al. 2000).Furthermore, X-ray modelling of the burst showed a clear signature of photospheric radius expansion (PRE), where the burst luminosity exceeds the local Eddington limit causing a radial expansion of the neutron star photosphere.The PRE X-ray burst was used to estimate the source distance of 4.9 ± 0.3 kpc (see, Kuulkers et al. 2003, for a review of PRE bursts as standard candles).However, we note that the quoted distance error is purely statistical, as it does not take into consideration any systematic effects, such as the potential for the neutron star to deviate from the assumed mass of 1.4 ⊙ or the potential for accreting elements besides hydrogen.Therefore, the error on the distance is likely an underestimation.An analysis of multiple Type I X-ray bursts detected during the 2007 outburst showed timing signals consistent with a neutron star spin frequency of 531.8 Hz (Bilous et al. 2018).These 'millisecond burst oscillations' are thought to be caused by anisotropic X-ray emission (i.e., 'hot spots'; Watts 2012) and allow for the determination of the neutron star spin frequency without the need for consistent pulsations.
The source has not been classified as an atoll or Z source; instead, it has adopted the broader label of neutron star 'soft X-ray transient', which encompasses both sub-classes.However, given its moderate peak X-ray luminosity (  ≤ 4 × 10 36 erg s −1 ) and transient behaviour, it is likely to be an atoll source.The majority of Z sources are persistent and bright, with maximum X-ray luminosities reaching appreciable fractions of the Eddington limit (e.g.,   ∼ 2 × 10 38 erg s −1 ).
On 2021 May 13 (MJD 59347), the gas slit camera (GSC) aboard The Monitor of All-sky X-ray Image (i.e, MAXI; Matsuoka et al. 2009) satellite detected the X-ray brightening of SAX J1810 as it entered its fifth recorded outburst (Iwakiri et al. 2021).Following the X-ray detection, radio observations with the MeerKAT radio telescope on 2021 May 21 (MJD 59356) revealed a spatially coincident radio source, constituting the first radio detection of this source (Motta et al. 2021).Here we present our multi-instrument radio/Xray monitoring campaign of SAX J1810.Our monitoring includes the 2021 outburst and 2023 follow-up that revealed the existence of a spatially coincident, persistent steep spectrum radio source.The remainder of this paper is structured as follows: in Section 2, we introduce our observation and analysis procedure, while in Sections 3 and 4, we present and discuss our results.Finally, we summarize our findings in Section 5.

Weekly Monitoring
We observed SAX J1810 with MeerKAT (a radio interferometer; Camilo 2018) as a part of the large survey project ThunderKAT (Fender et al. 2016).We began a weekly monitoring campaign on 2021 May 22 (MJD 59356), nine days after the outburst's initial detection, and continued until 2021 October 23 (MJD 59508) for a total of 21 observations.Each observation consisted of a single scan of 15 minutes on-source flanked by two 2-minute scans of a nearby gain calibrator (J1830-3602).Each epoch also included a 5minute scan of PKS B1934-638 (J1939-6342) for flux and bandpass calibration.In addition to the weekly monitoring, we observed two deep (1-hour) epochs on 2023 May 22 (MJD 60086) and 2023 August 16 (MJD 60172) when the source was in (X-ray) quiescence.The deep epochs followed the same observing strategy, except the source monitoring was broken into two 30-minute scans.All MeerKAT observations used the L-band receiver, with a central frequency of 1.3 GHz, and a total (un-flagged) bandwidth of 856MHz split evenly into 32768 frequency channels.To decrease the size of each data set, we averaged together every 32 channels (resulting in 1024 total channels) before data reduction and imaging.This averaging will not affect our final results as we are focused on radio continuum emission (as opposed to spectral lines).
We performed flagging, calibration, and imaging using a modified version of the semi-automated routine OxKAT 1 (Heywood 2020), which breaks the process into three steps.Here we will briefly outline the workflow and direct readers to Heywood et al. (2022) for a more comprehensive description.The first step (1GC) uses casa (v5.6;McMullin et al. 2007) to remove data corrupted by radio frequency interference (RFI).After removing RFI, the data is corrected with standard calibration solutions (i.e., flux density, bandpass, and complex gain).The second step (FLAG) applies a second round of flagging using tricolor (Hugo et al. 2022) before creating a preliminary image of the source field using wsclean (v2.9;Offringa et al. 2014).This preliminary image is then used to create an imaging mask.The final step (2GC) begins with a masked deconvolution before using the model image for direction-independent (DI) selfcalibration with CubiCal (Kenyon et al. 2018).Following self-cal, the pipeline ends with a second round of masked deconvolution using the DI self-calibrated visibilities.We adopted the 2GC images as our final data products.We maximize our sensitivity by weighting each image with a Briggs' robustness of 0 (Briggs 1995) 2 .We note that OxKAT has the functionality to solve for direction-dependant (DD) self-calibration solutions if needed (i.e., the 3GC step).However, for SAX J1810, DI self-calibration was sufficient, and thus we omitted the 3GC step.
We measured the source properties in each epoch using the casa task imfit, fitting an elliptical Gaussian component in a small subregion around the source to measure the position and flux density.As the source was unresolved, we set the component shape to be the synthesized beam of each image.We quantified the (1) uncertainty on the flux measurement using the local root-mean-square (rms) noise.We extracted the rms from an annular region for each epoch using the casa task imstat.Each annulus was centered on the position of the Gaussian component.We fixed the inner radius as the major axis of the synthesized beam and scaled the outer radius such that the annular area comprises the area of 100 synthesized beams.We quantified astrometric errors using the method detailed in Appendix A.

Very Large Array
We were approved for a single director's discretionary time observation (Project Code: 23A-417) with the Very Large Array (VLA) as a follow-up of our initial 2023 MeerKAT observation.SAX J1810 was observed on 2023 July 17 (MJD 60142) in the 2-4 GHz (S-band) and 4-8 GHz bands (C-band).For S-band, the observations used the 8-bit sampler comprised of two base-bands, with eight spectral windows of sixty-four 2 MHz channels each, giving a total (unflagged) bandwidth of 2.048 GHz.The 3-bit sampler was used for C-band, which has four base-bands, and thus a 4.096 GHz bandwidth.In each band, we included a single 1-minute scan of the flux calibrator (3C286).For source monitoring the array cycled between SAX J1810, observed for ∼ 8 minutes per cycle in S-band and ∼ 5 minutes in C-band.Each source scan is flanked by ∼ 1 minute observations of a nearby gain calibrator (J1820−2528).The total time on source was ∼ 16 minutes in both bands.We performed flagging, calibration, and imaging using the most recent release of the casa VLA pipeline (v6.4).We imaged the source using wsclean but did not detect the source in either band.As a result, we extract the rms noise from each image to place (3) upper limits on the flux density.We used a circular extraction region (with an area equal to 100 synthesized beams) centered on the archival position of SAX J1810 to measure the rms.The radio flux densities from both MeerKAT and the VLA are presented in Table C1 2.3 Swift-XRT

Weekly Monitoring
We monitored SAX J1810 with the X-ray telescope (XRT; Burrows et al. 2005) aboard the Neil Gehrels Swift Observatory (Gehrels et al. 2004), capturing the quasi-simultaneous evolution of the X-ray flux (i.e., within ∼ 3 days of a MeerKAT observation).During the outburst, we observed 21 epochs (target ID: 32459) between 2021 May 20 (MJD 59364) and 2021 November 6 (MJD 59524) at an approximately weekly cadence.To accompany our deep MeerKAT epochs, we were approved for two Target-of-Opportunity observations on 2023 May 25 (MJD 60089) and 2023 August 16 (MJD 60172).During the initial stages of the outburst, we monitored the source in Windowed Timing (WT) mode, where SAX J1810 exhibited a maximum count rate of ∼ 20 count s −1 during the first epoch.We transitioned to Photon Counting (PC) mode when the sources count rate decayed to ≲ 1 count s −1 on 2021 October 9 (MJD 59496), although there was a single intermittent PC epoch on 2021 September 5 (MJD 59462).
We used the Python API version of the Swift-XRT pipeline, swifttools (Evans et al. 2007(Evans et al. , 2009)), to extract the source and background spectra for all epochs except 2021 August 7 (MJD 59433), where the source exhibited a Type I X-ray burst (see section 2.3.2).We used the HEASOFT package (version 6.25) for our spectral analysis.For observations that had a sufficiently large number of counts (i.e., MJD 59364-59496), we used a modified grppha script to bin the spectra on 25-count intervals and performed spectral fitting using  2 statistics.Towards the end of our 2021 monitoring (i.e., the MJD 59504 and 59511), we used Cash statistics (i.e., cstat; Cash 1979) with single-count binning intervals, due to the small number of counts collected in each observation.The final two epochs of the 2021 monitoring (MJD 59518 and 59524) and the late-time follow-up (MJD 60089 and 60172) were non-detections and thus were omitted from the spectral fitting routine.
Using xspec (Arnaud 1996), we performed our spectral fitting twice, once for the 0.5-10 keV energy range and again for 1-10 keV.As expected, changing the energy range had a negligible effect on the best-fit spectral parameters.We modelled the spectra using an absorbed power law model with an added blackbody component; i.e., tbabs × (pegpwrlw + bbody), where tbabs models the interstellar absorption using an equivalent hydrogen column density (  ) following the abundances from Wilms et al. (2000).The power law accounts for the X-ray emission from the dominant component (i.e., the hard X-ray corona), and the blackbody accounts for any excess soft X-ray emission from a faint accretion disk, neutron star surface, or boundary layer.Initially, we fit each spectrum individually, allowing   to vary epoch by epoch.We then adopted the single epoch fitting as our starting parameters, linking the   values across all epochs and fitting the spectra simultaneously, resulting in a single time-independent value of   .When calculating the degrees of freedom, we treated the linked   as frozen (i.e., each spectrum has four free parameters).The epochs that utilized Cash statistics were omitted from the fitting procedure detailed above.Instead, we fit each of those spectra with a simple absorbed power law model (i.e., tbabs × pegpwrlw), fixing   to our best-fit value of 3.88 × 10 21 cm −2 and the power law photon index (Γ) to the average value of 1.61 from the  2 fitting.As a result, the X-ray flux was the only free parameter in the Cash statistic modelling.The Swift-XRT monitoring and spectral parameters during the 2021 outburst are presented in Table C2.The quoted uncertainties on the X-ray parameters represent the standard 90% confidence intervals.

Type I X-ray Burst
On 2022 August 7 (MJD 59433), SAX J1810 underwent a Type I Xray burst, and, as a result, we performed manual data reduction on the Swift-XRT (WT) observations.First, we ran the task xrtpipeline to produce cleaned event files and exposure maps.Second, using barycorr, we applied the barycentric timing correction.Lastly, we extracted source and background spectra by using xselect.For the pre-burst times, we used a circular source extraction region with a radius of 30 pixels (1 pixel = 2.36 arcsec) and an annular background extraction region with an inner radius of 70 pixels and an outer radius of 130 pixels.The pre-burst spectrum was then processed using  2 statistics and the routine mentioned in 2.3.1.
During the burst, we broke the event file into multiple time bins to analyze the time evolution of the spectral parameters.Due to high count rates during the burst (i.e., maximum count rates ≳ 400 count s −1 ), the observations are affected by systematic effects caused by photon pile-up.As a result, we used an annular source extraction region with an inner (exclusionary) radius that increases with an increasing count rate (ranging from 0 to 3 pixels).Following the Swift-XRT pipeline procedure (see, Evans et al. 2007Evans et al. , 2009)), we choose inner radii that reduce the maximum count rate in a given time bin to < 150 count s −1 .The time ranges were chosen so each bin has ≳ 300 counts corresponding to 21 bins across the 1.5 minute burst.To model the burst parameters in xspec we added a second blackbody component to the pre-burst spectrum, fixing the pre-burst parameters, thereby allowing only the second blackbody to vary.We used the bbodyrad model to directly fit for the normalized radius (i.e., size of the blackbody) and temperature before using the xspec convolution model cflux to calculate the flux.
For the timing analysis, we extracted two light curves.The first light curve was binned on 1 s intervals and was used to model the decay timescales of the burst.We extracted an initial light curve using the circular extraction region.For any time bins with a count rate > 150 count s −1 , we replaced their count rates with the count rate measured by the annular region with a 3-pixel exclusionary inner radius.We corrected for background and annular extraction region effects with lcmath and xrtlccorr, respectively.Following the prescription outlined in (Galloway et al. 2020) we fit an exponential decay function, where  is the time after the burst maximum, () is the count rate at a given ,  0 is the constant background rate,  is the folding decay time, and  is the peak count rate of the bursting component (excluding the contribution from a constant background).We fit for ,  0 , and  with a Markov-Chain Monte Carlo (MCMC) routine using Python's emcee package (Goodman & Weare 2010;Foreman-Mackey et al. 2013), assuming the sampled count rates were independently distributed normal random variables.The number of (sampling) walkers was fixed at five times the number of dimensions (i.e., 15).We chose three flat priors to ensure an unbiased analysis.To ensure convergence, we manually inspected the walkers over many autocorrelation times.Additionally, we analyzed the evolution of the autocorrelation time as a function of the number of MCMC steps following the routine outlined in the emcee documentation3 .The second light curve was extracted using the circular extraction region and binned on 1.8 ms intervals (the minimum bin size possible for WT mode).We used the short timescale light curve to search for millisecond burst oscillations.Given the short timescale binning, no corrections were applied to the 1.8 ms light curves.Appendix B presents the X-ray burst properties.

The WATCHDOG Pipeline
We calculated the X-ray hardness ratio (HR) using a modified version of the pipeline developed for the Whole-sky Alberta Time-resolved Comprehensive black hole Database Of the Galaxy (WATCHDOG; see Tetarenko et al. 2016, for a comprehensive description of the pipeline).The hardness ratio is the ratio between the number of counts in the hard and soft X-ray bands.We used the MAXI/GSC 4-10 keV band as the soft band and 15-50 keV observations from the Burst Alert Telescope (BAT; Barthelmy et al. 2005) aboard Swift as the hard band.Both sets of observations are publicly available4 .We modified the pipeline to average daily observations, ensuring the hard X-ray band had a ≥3 detection.For data where the soft X-ray band detection significance was < 3, we replaced the measured count rate with 3× the noise value to estimate a conservative 3 lower limit.The source appears to have undergone a hard-only outburst, and, as a result, to get meaningful constraints, we needed to measure either a lower limit or detection on the hardness ratio.No further modifications were applied to the WATCHDOG pipeline.
WATCHDOG defined empirical HR limits that corresponded to the different X-ray states: (i)  hard = 0.3204; and (ii)  soft = 0.2846.A hardness ratio is considered consistent with the hard (soft) state if its lower (upper) error bars are above (below) the  hard ( soft ) limits.If neither criterion is met, the source is classified as being in an intermediate state.We note that the values of  hard / soft were calculated for BHXBs; in Section 4.1, we investigate whether it is valid to apply the same standard NSXBs.

Radio Position
In Fig. 1, we show the offset in right ascension and declination between the MeerKAT position and the archival X-ray position of 18h10m44.47s−26 • 09 ′ 01.2 ′′ from (Jonker et al. 2004).The average radio position is 18h10m44.34s−26 • 09 ′ 02.1 ′′ (±0.1 ′′ ).The perepoch declinations are consistent with the average radio position with a reduced  2 = 0.75 (22 degrees of freedom), although the average radio position is offset by ∼ 1 ′′ from the X-ray position.In contrast, the right ascensions show significantly larger offsets ranging from ∼ 1-5 ′′ .Moreover, the measured right ascensions show temporal variability.Adopting the weighted mean offset in right ascension as a model and computing the reduced  2 results in a value of  2 = 4.4 (22 degrees of freedom), suggesting that the variability is not the result of stochastic error fluctuations.We tested the right ascension offsets against a linearly increasing model (i.e., ballistic motion), which resulted in a negligible improvement in the reduced  2 (4.2;21 degrees of freedom), and thus, we found no evidence of ballistic motion.

Outburst Light Curves
In Fig. 2 we show the MeerKAT (1.3 GHz; top panel), Swift-XRT (0.5-10 keV; second panel), MAXI/GSC (4-10 keV; third panel), and Swift-BAT (15-50 keV; bottom panel) outburst light curves.For our MeerKAT observations, 18 (out of 21) epochs were ≥ 5 detections (blue circles).The remaining three epochs (blue diamonds) do not meet the typical reporting threshold of 5, with detection significance of ∼ 4.3-4.9.Given the spatial coincidences between the low (< 5) and high-significance detections (≥ 5), it is likely that we are detecting a source in all of our MeerKAT observations.For the Swift-XRT light curve, we adopted the total fluxes from our spectral fits using the joint power law and blackbody model components (filled black circles).The last two data points (open black circles) correspond to the epochs where the source was too faint for multicomponent spectral modelling; instead, we fit the source with a single power law component.The Swift-BAT and MAXI/GSC light curves display the data at a daily binning frequency.
The observed flux of SAX J1810 displays a common temporal evolution across all observing frequencies.At early times (∼ MJD 59340-59370), all four instruments recorded the brightest signal of the outburst.Following the maxima, the source flux began decreasing, showing a rebrightening between ∼ MJD 59410 and 59440, before the source flux continued to decrease, returning to X-ray qui- escence and plateauing at ∼ 90 Jy in the radio.We find no evidence for additional intra-observation variability beyond the Type I outburst discussed in this paper.
Although the radio and X-ray light curves share a similar evolution in time, the magnitude of the variability is significantly different.In radio, the source exhibits modest variability with a maximum (∼ 230 Jy) and minimum (∼ 80 Jy) flux density separated by a factor of only ∼ 3.In contrast, when only considering the epochs with multi-component spectral modelling, the Swift-XRT fluxes show a factor of ∼ 20 in variability, with a maximum and minimum flux of ∼ 1.6×10 −9 and 6.8×10 −11 erg s −1 cm −2 , respectively.Including the final two Swift-XRT epochs during the source's return to quiescence, the minimum flux is ∼ 5 × 10 −13 erg s −1 cm −2 , which corresponds to a factor of ∼ 2000 decrease from the maximum.The plateauing radio emission at MJD 59463 (and beyond) is consistent with a spatially coincident, persistent radio source (see Section 4.2).

X-ray Spectra
The X-ray modelling parameters are shown in Fig. 3.The best fit equivalent hydrogen column density is   = 3.9 +0.1 −0.2 × 10 21 cm −2 .The Colden: Galactic Neutral Hydrogen Density Calculator5 estimates a value of   ∼ (3.2-4.3)× 10 21 cm −2 along the SAX J1810 line of sight (depending on the choice of neutral hydrogen data set -NRAO or Bell), making the measured   consistent with expectation.
To investigate the relative contributions of each model component, we calculated the power law flux fraction (third panel, Fig. 3); i.e.,  ,PL / ,tot , where  ,PL is the X-ray flux of the power law compo-  We fixed our best fit value to   = 3.88 × 10 21 cm −2 , and left all other parameters free.After including the diskbb component, the pegpwrlw becomes the subdominant component, and the fit becomes insensitive to both the flux and the photon index of the power law component.SAX J1810 may have briefly entered a thermal X-ray dominated accretion state before returning to the hard state.
nent and  ,tot is the total X-ray flux of the model.In all epochs, the power law component is dominant with a flux fraction ranging from ∼ 0.53 to 0.94 with a (variance-weighted) average of 0.72 ± 0.02.The power law photon index (Γ; fourth panel, Fig. 3) shows moderate variability with 0.4 +0.93 −0.43 ≤ Γ ≤ 2.88 +0.18 −0.08 and an average value of 1.61 ± 0.03.The average value is typical of comptonized hard state X-ray emission from (black hole) X-ray binaries (Remillard & McClintock 2006).Moreover, if we exclude the anomalously steep photon index, the maximum photon index becomes Γ = 1.83 +0.10  −0.08 .The blackbody temperature (; third panel, Fig. 3) varied between 0.5 +0.18 −0.08 ≤  ≤ 1.2 +0.18 −0.08 keV, with an average blackbody temperate of  = 0.60 ± 0.01 keV.Black body temperatures ≲ 1 keV are consistent with past analyses of hard state neutron star X-ray binaries (e.g., Lin et al. 2007).The bottom panel of Fig. 3 displays the hardness ratio calculated from the daily Swift-BAT and MAXI/GSC light curves.We observe a moderate degree of variability in hardness ratio, with detections ranging from ∼ 0.5-2.8, and an average value of 1.19 ± 0.06.Including the lower limits increases the maximum hardness ratio to ∼ 4.
The largest single epoch evolution occurs on MJD 59385, where the black body temperature reaches its maximum value of ∼ 1.2 keV, alongside the extreme softening of the power law component (Γ ∼ 2.9).During this epoch, the two-component fit had a reduced  2 value of ∼ 1.17 (216.5/186).To investigate whether we were observing a transition to an intermediate or soft state, we added a multicolour disk to the two-component model; i.e., tbabs × (pegpwrlw + bbody + diskbb).The inclusion of the third component moderately reduces the  2 to ∼ 1.12 (206.1/184) and decreases both the power law photon index and blackbody temperature to levels consistent with the other epochs (See Table 1 for the full model parameters).Moreover, the power law component becomes sub-dominant, suggesting that the source may have briefly transitioned into an intermediate or soft state.The observations on MJD 59413 and 59462 show similarly large reduced  2 values of ∼ 1.22 (237/194) and ∼ 1.52 (50/33), respectively.As a result, we attempted to fit these spectra with the same three-component model.However, the fitting resulted in a negligible improvement of the  2 statistic.We note that, for the latter epochs, both have reduced  2 deviations that are consistent (at the < 3 level) with the expected value of 1.Therefore, the poor fits may result from statistical effects rather than a physical change in the X-ray spectrum.(MJD 60169).The best-fit positions of both 2023 detections are consistent with the 2021 outburst (see Fig. 1).Therefore, we confidently detect a persistent radio source spatially coincident with SAX J1810.We calculated an (intra-band) spectral index of the persistent source using the brighter of the two MeerKAT follow-up observations (MJD 60086).We broke our observations into four evenly spaced sub-bands, ensuring a ≥ 5 detection in each sub-band.Applying a simple linear least squares fit, we measured a spectral index of  = −0.7 ± 0.5.In addition to the large statistical error, we note that intra-band spectral indexes are known to bias towards flatness ( ∼ 0) at detection significances ≲ 35 (Heywood et al. 2016).Given our source was only detected at ∼ 10 and the relatively large error bar, we do not apply any strong physical inference based on this intra-band spectral index During the last seven epochs of 2021 monitoring (MJD 59463 to 59511) -after the radio flux density had plateaued -the average radio flux density is 93 ± 7Jy.This value is consistent with our 2023 observations (at the ∼ 2 level), suggesting the persistent emission is, at most, weakly variable with a ∼ 20% excess variance.Combining the late-time 2021 and 2023 observations results in a (weighted) average flux density of 89 ± 5Jy.The quasi-simultaneous Swift-XRT follow-up on MJD 60089 and 60172 did not detect any spatially coincident X-ray source in either epoch setting 3 upper limits on the 1-10 keV X-ray flux of < 1.3 × 10 −13 erg s −1 cm −2 and < 3.0 × 10 −13 erg s −1 cm −2 , respectively.Furthermore, our scheduled VLA follow-up at 3 GHz and 6 GHz, taken between our two MeerKAT observations on 2023 July 17 (MJD 60142), did not detect the source.The 3 upper limits on the 3 GHz and 6 GHz were 30 Jy and 18 Jy, respectively.Adopting a 1.3 GHz flux density of 78 Jy (conservatively assuming a 3 drop in flux caused by intrinsic variability), we use the 3 GHz non-detection to calculate a conservative upper limit of  < − 1.1.2022), an updated version of the Bahramian & Rushton (2022) catalog.As our Swift-XRT and MeerKAT observations were quasisimultaneous, we applied a one-dimensional linear interpolation to map the radio observations onto the X-ray times for our 2021 observations.We did not apply any interpolation for our 2023 follow-up observations.Instead, we grouped the MeerKAT observations with the nearest Swift-XRT follow-up.We present the   -  relation from the 2021 outburst as red circles.Fitting the 2021 results with a simple power law results in a shallow exponent of  = 0.09 ± 0.03 (for   ∝    ).If we assume that the 2023 MeerKAT detections originate from a persistent hard state jet (purple stars on Fig. 4) and thus should follow the   -  relation, the measured power index becomes an upper limit (due to the X-ray non-detections) adopting a value of  < 0.06.Given that our results strongly suggest the existence of a persistent radio source that is unrelated to the hard state jet of SAX J1810, we present a secondary set of   -  data points (green squares) after subtracting off 93 Jy from each of the radio flux densities from our 2021 outburst.Post-subtraction, there are only four epochs (MJD 59364,59378,59413,and 59437) that show a > 3 excess flux density when compared to the persistent level.For the rest of the epochs, we set the radio flux density to be 3× the rms noise and displayed them as upper limits.The subtracted values are unconstraining but consistent with the broader population of NSXBs.The implications of SAX J1810   -  evolution and the origin of the persistent radio source are discussed in Section 4.2

DISCUSSION
We monitored the NSXB SAX J1810 during its 2021 outburst.The X-ray and radio properties suggest that the source underwent a 'hard-only' outburst, never fully transitioning to a soft accretion state.Moreover, the late-time plateau of radio flux density in 2021, combined with our follow-up in 2023, suggests the existence of a persistent radio source.In the following subsections, we present the evidence of a 'hard-only' outburst and discuss the possible origins of the persistent radio emission.

Hard-Only Outburst
Our observations suggest that SAX J1810 exhibited a 'hard-only' outburst in 2021.We justify this claim with three points of evidence: (i) The hardness ratio between the Swift-BAT and MAXI/GSC observations is above the hard state limit throughout the monitoring.Although the limit was empirically defined using outbursting BHXBs, we expect that the persistent source of thermal X-ray photons (from the neutron star surface or boundary layer) would make all X-ray states softer, thereby decreasing the hard state limit for NSXBs.We investigate this proposition by analyzing the best-studied outbursting (atoll) NSXB, Aql X-1.In Fig. 5, we have plotted a sample light curve of Aql X-1 during its 2016 outburst.The source exhibits a rapid transition of its hardness ratio, with a large fraction of the outburst remaining at a steady value of ∼ 0.05 well below the soft state limit derived for BHXBs.Díaz Trigo et al. (2018) performed an X-ray spectral analysis of four separate observations; the authors identified that the source was in the hard accretion state on 2016 Aug 3 (MJD 57603) and 2016 Sep 19 (MJD 57650) and in the soft accretion on 2016 Aug 5 (MJD 57605) and 2016 Aug 7 (MJD 57607).The hard and soft state epochs are shown with the dashed and dashed-dotted lines in Fig. 5.As expected, the soft and hard state epochs are temporally consistent with small and large hardness ratios.The final (Sep 19) hard state epoch shows a hardness ratio below the BHXB hard state limit, consistent with our prediction that the thermal photons from neutron stars will lower the hard state limits.We note that other outbursts of Aql X-1 (e.g., the 2009 outburst;Miller-Jones et al. 2010) show a similar 'softening' of the hard state limit.Therefore, we are confident that the Swift-BAT and MAXI/GSC hardness ratio for SAX J1810 is consistent with hard state emission throughout the 2021 outburst, and our adoption of the WATCHDOG limits is most likely appropriate (if not a conservative approximation).
(ii) Our Swift-XRT spectral modelling is consistent with hard state emission in nearly all epochs.The X-ray photon indexes (Γ avg ∼ 1.6) and low-energy black body temperatures ( avg ∼ 0.6) are typical of hard state X-ray emission from an NSXB (Lin et al. 2007).Moreover, the power law component is the dominant flux component in all epochs (i.e., power law flux fraction ≥ 50%).Although some epochs show approximately equal contributions between the blackbody and power law components, the narrow (0.5 − 10.0 keV) energy range favors the black body component when calculating band limit flux, as the power law component will dominate at higher energies (≥ 10 keV).The bolometric X-ray flux is more strongly dominated (> 90%) by the power law component than our observations would suggest, consistent with hard state emission.The anomalous epoch (MJD 59385; Table 1) that shows a clear softening of the X-ray spectrum suggests the source may have exhibited a brief deviation from a hard accretion state.Assuming a successful transition to the soft state, and given the cadence of our observations and the bracketed hard sate epochs, the source would have gone through a full cycle (i.e., hard → soft → hard) in ≤ 14 days before remaining in the hard state for the remaining ∼ 120 days of outburst (atypical behaviour for an outbursting NSXB, see, Muñoz-Darias et al. 2014, for a review of outburst timescales).We find it more likely that the source briefly entered an intermediate state, failed to complete a transition to the soft state, and transitioned back to the hard state.
(iii) The evolution of our radio observations is consistent with the hard state.First, the radio and X-ray light curves show a correlated temporal evolution characteristic of hard state emission.Second, we do not detect any significant jet-quenching.Although radio emission from NSXBs has been observed in the soft state, when both hard and soft state (compact jet) radio emission has been detected, the jet emission is brighter in the hard state (at a fixed X-ray flux, e.g., Gusinskaia et al. 2017).Therefore, without a significant increase in the X-ray flux (which was never observed), we would expect a decrease in the radio flux after a transition to the soft state.We recognize that the spatially coincident, persistent radio source contaminates our ability to detect jet-quenching.However, the persistent source can not explain the joint radio-X-ray time evolution, as we would expect the radio flux to drop to the persistent level (∼ 90 Jy) without a similar decrease in X-ray flux.Whenever we observed an increasing X-ray flux, we observed a simultaneous increase in the radio flux density.
Comprehensive monitoring campaigns of future outbursts of SAX J1810 will be critical for confirming whether the source consistently exhibits 'hard-only' outbursts or shows a broader outburst phenomenology that sometimes results in successful transitions to the soft state (as observed in some BHXBs, e.g., H1743-322;Coriat et al. 2011;Williams et al. 2020).

The Origin of the Persistent Radio Emission
Our observations strongly support the existence of an unresolved, persistent, steep-spectrum radio source spatially coincident with the position of SAX J1810 (± 3 ′′ ).Considering the source exhibited a 'hard-only' outburst in 2021, we expect the radio emission to (partially) originate from a hard state jet (i.e., compact jet).The temporal coincidence between the flares at X-ray and radio frequencies is strong evidence for the existence of a steady jet.Moreover, the persistent source is weakly variable with an average flux density of ∼ 90Jy.Considering that we have multiple detections at ≳ 200Jy, we have clearly detected radio emission from the compact jet.
However, a hard state jet associated with SAX J1810 cannot be the source of the persistent radio emission.Hard state jets are stationary and, therefore, would not exhibit the proper motion that we have observed (Fig. 1) Moreover, the locations of its luminosities on the   -  plane (red circles Fig. 4) are inconsistent with a hard state jet.At early times and high X-ray luminosities, the radio/X-ray luminosities are positively correlated, as expected from a compact, steady jet.Towards the end of the outburst (at   ≲ 5 × 10 35 ergs s −1 ), there is a clear flattening of the correlation resulting in a < 0.06 due to the radio luminosity remaining approximately constant while the X-ray luminosity decreased by over three orders of magnitude.The 2023 follow-up, in particular, would make SAX J1810 exceptionally radioloud for a NSXB, consistent with the population of BHXBs.Recent analyses estimate a value of  = 0.44 +0.05 −0.04 for the total population NSXBs, with the atoll sub-population (which SAX J1810 is likely a member of) having  = 0.71 +0.11  −0.09 (Gallo et al. 2018).Both values of  reject our measurements at the > 3 level.Therefore, the observed radio emission likely originates from two components, with the most likely candidates of the persistent emission being either a discrete jet ejection or an unrelated, spatially coincident source.
We disfavor an origin due to jet ejection(s).First, the average decay timescale of an ejection event is ≪ 1 year, and thus a jet ejection persisting for ∼ 2 years and showing no significant decrease in the measured flux density is, in itself, unlikely.Long-lasting jet ejecta have been observed from BHXBs and are thought to be the result of jet-ISM interactions driving in situ particle acceleration and longterm synchrotron emission (e.g., Corbel et al. 2005;Bright et al. 2020;Carotenuto et al. 2021a;Bahramian et al. 2023).However, such long-lasting ejecta have never been observed in NSXB (likely due to their weaker, lower-luminosity jets being unable to power such long-term emission), and when observed in BHXBs, the radio emission of long-lived ejecta is strongly variable.Second, our VLA follow-up observations suggest a 3 upper limit on the radio spectral index of  < −1.1, significantly steeper than expected from opticallythin synchrotron emission from a jet ejection ( ∼ − 0.7).Lastly, our observations show no evidence of ballistic motion despite the source The number of sources within the inclusion radius.We include data for both the unresolved (black line) and for unresolved + extended source populations (red line).We note that the sharp increase and peaks close to SAX J1810 correspond to a regime susceptible to low-count statistics.Regardless we adopt the peak of the red curve as the most conservative estimate of the chance coincidence probability.
persisting for ∼ 2 years, which would be the strongest evidence for a jet ejecta origin of the persistent emission.If the persistent emission originated from jet ejecta, we would have had to observe a longlasting, non-variable, spectrally steep ejecta showing no motion on the sky.Therefore, we can rule out a jet ejecta origin with high confidence.
To estimate the probability of a spurious spatial coincidence with an unrelated source in the field we used the Python Blob Detector and Source Finder (PyBDSF; Mohan & Rafferty 2015) to make a catalog of all sources (in each image) with a flux density > 74 Jy (3 lower than the average persistent radio flux density).We use the deep 2023 observations as their lower rms noise (10 Jy vs. 20 Jy in 2021) makes PyBSDF less prone to mistaking spurious noise spikes as real sources.Due to flux variability, each image catalog has a different number of sources.As a result, we conservatively use the 2023 May 22 image as it has more sources than the August observation and, therefore, a larger source density.We calculate the source density and then convert it to the expected number of sources within a 3 ′′ radius.The choice of 3 ′′ was motivated by the scatter of our best-fit positions.Using the expected number of sources, we then calculate the Poissonian probability of a chance coincidence of one or more unrelated background sources.The instrument's sensitivity decreases as a function of radial distance from the phase center of the array, and thus, there is a progressively smaller number of sources cataloged at larger separations from the phase centre (decreasing the source density).We applied a cut when calculating the probability to investigate this potential bias, only including sources within a certain distance from the phase centre in our calculations.In Fig. 6, we show the chance coincidence probability as a function of the aforementioned 'inclusion radius' for only unresolved sources (following the criteria from Appendix A) and for both unresolved and extended sources (all sources).We adopt the peak value for all sources as our conservative estimate of the chance coincidence probability (i.e., ∼ 0.6%).
Radio-bright active galactic nuclei (AGN) are the dominant population of unresolved background sources.However, background AGN have an average spectral index of  ∼ − 0.7.We use two recent surveys of background AGN spectral indexes to estimate the probability of finding a steep spectrum AGN.Randall et al. (2012) calculated the spectral index of 166 AGN using 325, 610, and 1400 MHz flux densities.Only 43 sources had an  < − 1.1 corresponding to a probability of ∼ 26%.In a more recent, larger sample size survey, de Gasperin et al. ( 2018) measured the spectral indexes of ∼ 540000 radio sources (using 147 and 1400 MHz flux densities), with only a subset of ∼ 32000 having an appropriately steep .The corresponding probability is ∼ 6%.Adopting the older catalog probability as a conservative estimate, we calculate the total probability of finding a spurious radio AGN with a sufficiently steep spectral index as ∼ 0.16% (a ∼ 3.2 event).Alternatively, the spectral index could suggest an origin from a class of sources known to have steep spectral indexes.The most common steep spectrum source is pulsars, with average spectral indexes of ∼ − 1.6 (Jankowski et al. 2018).We searched the Australian Telescope National Facility pulsar catalog (Manchester et al. 2005) for any nearby known radio pulsars but found no pulsars within a radius of 0.6 • .Given that there are only 3000 known radio pulsars (corresponding to an expectation value of ∼ 2 × 10 −7 pulsars within a 3 ′′ radius), there is a chance coincidence probability of ∼ 0.002%.When considering that pulsars tend to be distributed in the Galactic plane (∼ 20% of the sky), and SAX J1810 is also in the galactic plane, the chance coincidence probability would increase by a factor of ∼ 5 but is still less likely than the AGN scenario.We note that the persistent emission would correspond to a time-averaged flux of a pulsar; as a result, recent surveys that looked at this part of the sky would have detected a pulsed source (e.g., Keith et al. 2010).Moreover, MeerKAT's pulsar timing backend (i.e., MeerTRAP Sanidas et al. 2018) was operational during all of our observations but did not detect any pulsed emission from the source.Therefore, our estimated coincidence probability between SAX J1810 and an unknown pulsar is most likely an overestimate.
There is a small possibility that the persistent radio-emission is local to SAX J1810.Transitional millisecond pulsars (tMSPs)accreting neutron stars that transition between accretion-powered (i.e., NSXB-like) and radio pulsar behaviour -have shown anomalously bright radio emission while actively accreting.For instance, the tMSP, 3FGL J0427.9−6704, was measured at a point on the   -  relation that was also more consistent with the population of black hole X-ray binaries; however, its X-ray luminosities were a factor of ≳ 3 larger than our upper limits on MJD 60086 (e.g., Li et al. 2020).Other tMSPs (i.e., PSR J1023+0038) have even exhibited anti-correlations between radio and X-ray luminosities, which could allow for bright radio emission absent any X-ray detections (Bogdanov et al. 2018).
However, the properties of SAX J1810 are inconsistent with what is expected from tMSPs.Firstly, SAX J1810 does not show radio pulsations during X-ray quiescence (although eclipses or highly compact, elliptical binary orbits can prevent the detection of pulsations from tMSPs Lorimer & Kramer 2004;Papitto et al. 2013).Second, at X-ray luminosities ≤ 10 33 erg s −1 , tMSPs spectra are non-thermal (Γ ≤ 1.7 Linares 2014; Bogdanov et al. 2018;Li et al. 2020), whereas SAX J1810 is thermally dominated (Γ ≥ 3 Jonker et al. 2004;Allen et al. 2018).Lastly, SAX J1810 does not exhibit any of the rapid X-ray variability that results from switching between different accretion modes (during outburst), showing, at most, modest variability (Allen et al. 2018).Although it cannot be conclusively ruled out, we find it unlikely that the persistent radio emission results from SAX J1810 being a tMSP.
Local emission, tMSP or otherwise, is difficult to reconcile with the variability in the position, as the source is spatially unresolved.Using the scatter in the measured position (∼ 3 ′′ ) as a proxy for the expected separation of the two-source scenarios (i.e., the persistent emission is non-local), then observations by an instrument with sufficient angular resolution and sensitivity (e.g., the VLA in A-configuration or the Square Kilometer Array) during future outbursts when the compact jet is 'on' should be able to spatially resolve two distinct components.If only a single source is observed, and there continues to be temporally correlated evolution in the radio/Xray light curves, this would strongly support the scenario where the persistent radio emission is local to SAX J1810.

SUMMARY AND CONCLUSIONS
We have presented our ∼ 2 year joint radio and X-ray monitoring of the neutron star X-ray binary SAX J1810.8−2609.Our observations include dense (i.e., weekly cadence) observations during the source's 2021 outburst and a collection of late-time observations in 2023.The X-ray spectral properties suggested that the source remained in the hard state throughout the entire 2021 outburst.Moreover, the radio and X-ray luminosities show a temporally correlated evolution, characteristic of a hard state radio jet.We discovered a spatially coincident, persistent steep-spectrum radio source that shows no correlation with the simultaneous X-ray flux.Therefore, during the outburst, the radio emission originated from a superposition of two components: a variable hard state compact jet (≲ 100Jy), and the unknown persistent source (∼ 90Jy).The spectral index and evolution of the persistent source are inconsistent with jet ejecta.We conservatively estimated the probability of a chance coincidence with an unrelated spectrally steep background source, and although low (∼ 0.16%), a background AGN seems to be the most plausible scenario.SAX J1810.8−2609 is known to go into outburst every ∼ 5 years, and future outbursts should focus on identifying the source of the persistent emission.Of the current generation of radio telescopes, the VLA (A-configuration) and the Very Long Baseline Array (VLBA) both have sufficient angular resolution and sensitivity to resolve two ∼ 100 Jy sources (assuming a separation of ∼ 3 ′′ ).Moreover, nextgeneration radio interferometers, such as the Square Kilometer Array (SKA; of which MeerKAT is a pathfinder), would be able to reach the desired sensitivity with a fraction of the observing time (i.e., ∼ 10 Jy rms for ≲ 3 minutes on source; Braun et al. 2019).During the next outburst, if a second unrelated source is ruled out, follow-up observations should focus on understanding what physical mechanism is driving the persistent radio emission, whether the source is a tMSP or otherwise.tinuous region of pixels with a flux value above a user-defined threshold and at least one pixel has a flux larger than a higher (also user-defined) threshold.For large islands (i.e., extended emission), PyBDSF will fit multiple sources to a single island.For our fitting, we used 3 and 4 for our thresholding.solve for  and  independently along the RA and Dec directions.Below, we outline our fitting routine: (i) For each source, calculate an average SNR and an average position.Calculate the RA/Dec offset from the average position for every source in each epoch using the average position.
(ii) Estimate the error in the astrometric precision of each source by bootstrapping the offsets, adopting the median value of the bootstrapped sample as an initial guess for  and the ranges between the median and the 15 th / 85 th percentiles as the 1 (−)/(+) uncertainties (Δ  ).
(iii) Using the  estimates and the average SNR, solve for the scaling parameters  and  (i.e., the uncorrected fit).The fit implements an MCMC routine and follows the same approach detailed in 2.3.2.
(iv) Solve for the (inverse-variance weighted) average offset of all sources in each epoch (i.e., the epoch-to-epoch correction) weighting each offset using the uncorrected fit.
(v) Correct the source offsets with the epoch-to-epoch correction and re-solve for  and  with the updated -corrected -offsets.
(vi) Repeat (ii)→(v) until the fitting converges on solutions for  and .We defined a convergence parameter  = (  −  −1 )/Δ  ; i.e., the difference between the astrometric error of a source for the current () and previous ( − 1) iterations in units of Δ  .The fit is said to have converged after three consecutive iterations with a mean value of  < 0.1.The post-convergence fit is the corrected fit.Record the final epoch corrections.
The relative astrometric fitting is shown in Fig. A1 and the best-fit parameters are tabulated in Table A1.The uncorrected fits have reduced  2 values of ∼ 1.3 (123 degrees of freedom) in both RA and Dec. Applying the epoch-to-epoch corrections (i.e., the corrected fit) shows a significant worsening of the fit quality with a reduced  2 > 2, suggesting that a single per-epoch correction is not accurately capturing the time-dependent systematics in our observations, and a more complex epoch correction may be appropriate (e.g., one that accounts for distance and direction with respect to the phase center).We intend to expand upon this preliminary work to investigate whether the relative astrometric error is similar across a range of ThunderKAT fields.
The fits show that (for MeerKAT), the systematic threshold of the relative error is significantly lower than the commonly assumed limit of 10% the size of the synthesized beam.Moreover, the signal-to- The sources used for fitting (i.e., 125 sources within 0.3 • of the phase center) are given by the solid blue circles, and the total population of isolated point-like sources (612 sources) is shown as hollow black circles.The uncorrected fits are marginally acceptable for the fitted population (reduced  2 ∼ 1.3 for 123 degrees of freedom), although the corrected fits are poor (reduced  2 > 2.0 for 610 degrees of freedom).Furthermore, the fits (uncorrected and corrected) are poor matches to the total population of isolated point-like sources, suggesting that far-field effects (especially at high signal-to-noise ratios) are significant.Overall, the fits show that the systematic limit is well below 10% of the synthesized beam and that at SNR > 20, the global epoch affects (i.e., affecting every source in a given epoch) are the dominant astrometric error.
noise dependency is similar to the commonly assumed 1/(2 • SNR) scaling.Due to the residual issues in our modeling, for our SAX J1810 analysis, we conservatively rounded our uncorrected fit values, adopting  = 0.5 and  = 0.02 to quantify the relative astrometric errors.
To correct for absolute astrometry effects, we identified nine sources8 within our field of view that are used as phase calibrators for very long baseline interferometry (i.e., with positions measured at < 10 milliarcsecond precision).Eight of the nine sources met our unresolved and isolated requirement, and we used this sub-sample for absolute astrometric corrections.After applying the epoch-to-epoch correction from the relative astrometric fitting, we measured the offsets of the eight calibrators with respect to their known positions.We then calculate each epoch's weighted mean (weighting each source by their relative astrometric errors).Lastly, we calculated a single time-independent absolute astrometric correction (see Fig. A2).The epoch-to-epoch correction removed any (substantial) temporal variability, and, as a result, the per-epoch average offsets are consistent with a single (time-independent) RA/Dec offset.
The final astrometric error ( tot ) was calculated by adding (in quadrature) the relative astrometric precision (), the error on the epoch-correction ( epoch ), and the error on the absolute offset ( abs ), These are the errors shown in Fig. 1.We note that given the signalto-noise ratio of our SAX J1810 detections (SNR ≲ 10), the relative astrometry term, , dominates the quoted errors.where  is the source radius in units of km and  is the distance to the source in units of 10 kpc.
The burst began its rise at 14:14:12 on 2021 August 7 (MJD Table B1.X-ray burst timing parameters. 59433.59319), reaching a peak count rate of ∼ 400 counts s −1 with a rapid 7 ± 1 s rise time before decaying for the remainder of our observations.The timing fit converged on an e-folding decay time of  = 15.8 ± 0.2 s (full fit parameters in Table B1).The burst parameters are consistent with the MINBAR burst catalog Galloway et al. (2020) in both rise (3.4 +5.6 −2.4 s) and e-folding decay times (8 +21 −4 s).During its 2007 outburst, SAX J1810 exhibited 531.8 Hz oscillations in the light curves of a Type I X-ray burst, likely the result of the spin frequency of the neutron star (Bilous et al. 2018).Following the prescription outlined in Bilous et al. (2018) we searched for burst oscillations in our (1.8 ms resolution) light curves by calculating the power spectrum in sliding windows with widths of 0.5, 1, 2, and 4 s, where each subsequent window is offset by 0.5  from the previous one.We found no evidence of burst oscillations.However, the temporal resolution of Swift-XRT WT mode (1.8 ms) makes our power spectra insensitive to frequencies above ∼ 280 Hz.Assuming the oscillations result from the spin period of the neutron star, we do not expect the oscillation frequency to evolve drastically between the 2007 and 2021 outbursts.
Furthermore, SAX J1810 is known to exhibit PRE (i.e., during the 1998 outburst a PRE signature provided the current distance constraint of 4.9 ± 0.3 kpc; Natalucci et al. 2000).Therefore, we per-formed time-resolved intra-epoch spectral modelling to search for evidence of PRE.We observe some evolution of the radius and temperature, although the large errors greatly reduce their significance.Assuming a distance of 4.9 kpc, the radius of the blackbody component ranges from 3.2 +3.5 −2.6 to 6.7 +5.0 −4.2 km (i.e., from ∼ 5-to-100% of the neutron stars surface assuming a 10 km stellar radius).However, the evolution of the radius and temperature does not occur alongside a period of (approximately) constant X-ray flux; thus, we do not detect PRE.

APPENDIX C: DATA TABLES
This paper has been typeset from a T E X/L A T E X file prepared by the author.139.7 +7.9 −7.9

Figure 1 .
Figure 1.The right ascension (top panel) and declination (bottom panel) offsets for the best-fit SAX J1810 positions.The filled blue circles are the offsets of the source.The purple dotted line and cyan dashed line are the 2023 May 22 and 2023 August 13 offsets, respectively.The dashed-dotted black line is the archival X-ray position from Jonker et al. (2004), and the grey shaded area is the error on the archival position (± 0.6 ′′ ).Note the clear offset and temporal variability in the right ascension of the source.

Figure 2 .
Figure 2. Multi-instrument light curves of the 2021 outburst of SAX J1810.The top panel is the MeerKAT 1.3 GHz radio light curves showing both ≥ 5 (blue circles) and 4-5 (blue diamonds) detections during the 2021 outburst.The horizontal lines show the 2023 May 22 (purple dotted) and 2023 August 13 (cyan dashed) flux densities.The second panel is the Swift-XRT (0.5-10.0 keV) light curves.The filled and open circles correspond to the epochs fit with  2 and Cash statistics, respectively.The bottom two panels show the MAXI/GSC (third panel) and Swift-BAT (bottom panel) daily-binned light curves.All four instruments show a common temporal evolution characteristic of the correlation between radio and X-ray emission in the hard accretion state.

Figure 3 .
Figure 3.A summary of the spectral properties of SAX J1810.The top panel shows the MeerKAT radio flux density.The next two panels show the Swift-XRT X-ray flux (second panel) and the power law flux fraction (third panel) in the 0.5-10.0keV (filled circles) and 1.0-10.0keV (open circles) energy bands.The fourth panel shows the temperature of the black body component, and the fifth shows the power law photon index.The bottom panel shows the hardness ratio between the MAXI/GSC (4.0-10.0keV) and Swift-BAT (15.0-50.0keV) energy bands.The upper ( hard ) and lower ( soft ) dotted lines show the empirically defined state boundaries from WATCHDOG (Tetarenko et al. 2016).These spectral properties are characteristic of the hard accretion state MNRAS 000, 1-16 (2023) -  relation Our 2023 follow-up MeerKAT observations revealed a 112 ± 12 Jy radio (point) source on 2023 May 22 (MJD 60086) and another 75 ± 11 Jy radio source three months later on 2023 August 13

Figure 4
Figure 4 presents the   -  relation.The plot includes archival hard state BHXBs (grey circles), hard state NSXBs (blue squares), and accreting millisecond X-ray pulsars (AMXPs; orange triangles).The archival sources were adapted from Fig. 4 of van den Eĳnden et al. (2022), an updated version of theBahramian & Rushton (2022) catalog.As our Swift-XRT and MeerKAT observations were quasisimultaneous, we applied a one-dimensional linear interpolation to map the radio observations onto the X-ray times for our 2021 observations.We did not apply any interpolation for our 2023 follow-up observations.Instead, we grouped the MeerKAT observations with the nearest Swift-XRT follow-up.We present the   -  relation from the 2021 outburst as red circles.Fitting the 2021 results with a simple power law results in a shallow exponent of  = 0.09 ± 0.03 (for   ∝

Figure 4 .
Figure 4.The radio (5 GHz) and hard state X-ray (1-10 keV) luminosity (  -  ) relation.The archival values for black holes (grey circles), neutron stars (blue squares), and accreting millisecond X-ray pulsars (orange triangles) were taken from van den Eĳnden et al. (2022), which is based on the Bahramian & Rushton (2022) catalog.The luminosities for SAX J1810 during the 2021 outburst are represented as red circles.The purple stars are the 2023 values if we assume that all radio emission originates from a hard state jet.The late-time plateau in the radio flux density strongly suggested the existence of a second radio source is uncorrelated with the X-ray emission.The green squares show the 2021 outburst values after subtracting 89Jy from each of the radio flux densities (i.e., the average contribution from the persistent source).Even after subtracting off the persistent source, the two 3 radio detections (large green squares) of SAX J1810 remain consistent with the general population of hard-state NSXB jets.

Figure 5 .
Figure5.The X-ray evolution of the NSXB Aql X-1's 2016 outburst as seen by MAXI/GSC (top panel), Swift-BAT (middle panel), and the hardness ratio between the two instruments (third panel).The horizontal dotted lines adopt the same definition as BHXBs in Fig.3.The vertical dashed lines and dasheddotted lines show the times when the source was independently identified to be in the soft and hard accretion states, respectively(Tasse et al. 2018).Both soft accretion states occur at an HR ∼ 0.05, well below the empirically defined transition values.This suggests one can use the BHXB transition hardness ratio to conservatively estimate if a NSXB undergoes a 'hard-only' outburst.

Figure 6 .
Figure 6.(top panel) The probability of a chance spatial coincidence between SAX J1810 and an unrelated background source as a function of the 'inclusion radius'.(bottom panel)The number of sources within the inclusion radius.We include data for both the unresolved (black line) and for unresolved + extended source populations (red line).We note that the sharp increase and peaks close to SAX J1810 correspond to a regime susceptible to low-count statistics.Regardless we adopt the peak of the red curve as the most conservative estimate of the chance coincidence probability.

Figure A1 .
Figure A1.The relative astrometric fits for the population of isolated point-like sources: (top left) uncorrected Dec; (top right) corrected Dec; (bottom left) uncorrected RA; (bottom right) corrected RA.The sources used for fitting (i.e., 125 sources within 0.3 • of the phase center) are given by the solid blue circles, and the total population of isolated point-like sources (612 sources) is shown as hollow black circles.The uncorrected fits are marginally acceptable for the fitted population (reduced  2 ∼ 1.3 for 123 degrees of freedom), although the corrected fits are poor (reduced  2 > 2.0 for 610 degrees of freedom).Furthermore, the fits (uncorrected and corrected) are poor matches to the total population of isolated point-like sources, suggesting that far-field effects (especially at high signal-to-noise ratios) are significant.Overall, the fits show that the systematic limit is well below 10% of the synthesized beam and that at SNR > 20, the global epoch affects (i.e., affecting every source in a given epoch) are the dominant astrometric error.

Figure
FigureB1shows the parameters of the 2022 August 7 (MJD 59433) Type I X-ray burst.The top panel shows the 1s-binned light curves and the timing fits; the second panel shows the bolometric X-ray flux of the blackbody component; the third panel shows the temperature of the blackbody component; and the bottom panel shows the normalized radius of the blackbody component, defined as  2 / 2 , where  is the source radius in units of km and  is the distance to the source in units of 10 kpc.The burst began its rise at 14:14:12 on 2021 August 7 (MJD

Figure A2 .
Figure A2.The absolute astrometric corrections from the very long baseline interferometry calibrators in the SAX J1810 field-of-view.The open black circles are the offsets of each calibrator in each epoch.The closed blue circles are the average offset in each epoch.The dashed black line is the timeindependent average offset across all sources and all epochs.The blue region shows the 1 errors on the average offset.The per-epoch RA (Dec) average offsets are consistent with the time-independent value at a reduced  2 of ∼ 0.65 (∼ 0.24) for 22 degrees of freedom.These low  2 values suggest that we may be overestimating the relative astrometric error.

Figure B1 .
Figure B1.Spectral and timing fits from the Type I X-ray burst observed on 2022 August 7.The top panel shows the 1s-binned fight curves.We overlayed the best-fit exponential decay (dashed line) and the constant (pre-burst) count rate(dotted line); the timing fit parameters are tabulated in TableB1.The bolometric X-ray flux (top panel), temperature (2nd panel), and normalized radius (bottom panel) of the black body component do not show conclusive evidence of PRE.

Table 1 .
Three component fit of the Swift-XRT observation on MJD 59385.

Table A1 .
Relative Astrometry Parameters:  This column indicates whether the fitting omitted (uncorrected) or used (corrected) the epoch-to-epoch astrometry correction;  The fitting parameters  and  are expressed as a fraction of the synthesized beam FWHM for both the RA and Dec directions;  This column indicates the population of sources used for the corresponding  2 calculations.Population 1 is the nearby (< 0.3 • ) isolated point sources used in the fittings.Population 2 includes all isolated point sources, regardless of distance from the phase centre.The contrast between the  2 values of Population 1 and 2 highlights the effects of far-field errors.

Table C1 .
Radio properties of SAX J1810 MJD Date Instrument Central Frequency [GHz] Radio Flux Density [ Jy]