A new radio census of neutron star X-ray binaries

We report new radio observations of a sample of thirty-six neutron star (NS) X-ray binaries, more than doubling the sample in the literature observed at current-day sensitivities. These sources include thirteen weakly-magnetised ($B<10^{10}$ G) and twenty-three strongly-magnetised ($B\geq10^{10}$ G) NSs. Sixteen of the latter category reside in high-mass X-ray binaries, of which only two systems were radio-detected previously. We detect four weakly and nine strongly-magnetised NSs; the latter are systematically radio fainter than the former and do not exceed $L_R \approx 3\times10^{28}$ erg/s. In turn, we confirm the earlier finding that the weakly-magnetized NSs are typically radio fainter than accreting stellar-mass black holes. While an unambiguous identification of the origin of radio emission in high-mass X-ray binaries is challenging, we find that in all but two detected sources (Vela X-1 and 4U 1700-37) the radio emission appears more likely attributable to a jet than the donor star wind. The strongly-magnetised NS sample does not reveal a global correlation between X-ray and radio luminosity, which may be a result of sensitivity limits. Furthermore, we discuss the effect of NS spin and magnetic field on radio luminosity and jet power in our sample. No current model can account for all observed properties, necessitating the development and refinement of NS jet models to include magnetic field strengths up to $10^{13}$ G. Finally, we discuss jet quenching in soft states of NS low-mass X-ray binaries, the radio non-detections of all observed very-faint X-ray binaries in our sample, and future radio campaigns of accreting NSs.


INTRODUCTION
Accretion is a fundamental process occurring across the Universe in a plethora of objects and circumstances.These range from accreting supermassive black holes in Active Galactic Nuclei (AGN) and accreting white dwarfs, via stellar mass black holes and neue-mail: jakob.vandeneijnden@physics.ox.ac.uk tron stars in X-ray binary systems, to forming stars surrounded by proto-planetary discs.All these systems show states where the accretion of matter is observed to be accompanied by a coupled outflow of material, either in the form of wide-angled, relatively slow winds, and/or collimated and often relativistic jets.These outflows influence both the accreting systems, for instance contributing to the angular momentum loss in the accretion flow and reducing the effective accretion rate (e.g.Tetarenko et al. 2018b), and the sur-rounding medium (e.g.Gallo et al. 2005;Fabian 2012).Feedback from X-ray binaries and AGN is thought to contribute to the stellar feedback regulating star formation and ionising the early Universe (e.g.Fender et al. 2005;Mirabel et al. 2011;Justham & Schawinski 2012;Fragos et al. 2013a,b).
The formation of jets from an accretion flow is often fundamentally attributed to either one or a combination of two mechanisms.In both those models, twisted magnetic field lines close to the compact object launch material away from the accretion flow, but an important difference lies in the origin of the twisting of the magnetic field lines.In the Blandford & Znajek (1977) mechanism, these field lines are spun up as they thread the rotating ergosphere of a black hole.A key prediction of this model, that has remained difficult to test unambiguously, is the dependence of jet power on the black hole spin (e.g.King et al. 2013;McClintock et al. 2014;Russell et al. 2013b).Alternatively, for the many jet-launching systems do not contain a black hole, the Blandford & Payne (1982) mechanism proposes that the differential rotation of the accretion flow itself tangles up the magnetic field (note that, naturally, the Blandford & Payne (1982) mechanism can also occur in black hole systems).While this model applies to neutron stars, it predicts an upper limit on the magnetic field of neutron stars that can be spun up by an accretion flow.Therefore, it predicts that neutron stars with magnetic fields above a certain threshold should not launch jets (e.g.Massi & Kaufman Bernadó 2008;Migliari 2011).
While the accretion flow in AGN and X-ray binaries typically emits strongly in the X-ray band, the jet dominates at low frequencies through the emission of synchrotron radiation.This emission results from free electrons in the jet spiralling around magnetic field lines, producing a synchroton spectrum (Rybicki & Lightman 1979;Longair 1992).The observed spectrum of an unresolved jet depends on the jet type: discrete ejecta typically have steep spectra in the radio band, defined as radio spectral index α ≈ −0.7 (where Sν ∝ ν α ) as they emit as a single, optically thin population.A compact, steadily-outflowing jet is instead observed as the superposition of synchrotron spectra from different distances downstream, resulting in a flat (α = 0) to inverted (α > 0) spectral shape up to the jet break frequency.The highest frequency jet emission originates from the highest energy electrons, located near the base of the jet (Markoff et al. 2001;Corbel & Fender 2002;Markoff et al. 2005;Romero et al. 2017;Malzac 2013Malzac , 2014)), while lower frequencies are emitted further down the jet (e.g.Blandford & Königl 1979).The jet synchrotron emission can extend into the sub-mm (Russell et al. 2014;Tetarenko et al. 2015;Díaz Trigo et al. 2018), nIR, and optical (Russell et al. 2006(Russell et al. , 2007(Russell et al. , 2013b,a;,a;Gandhi et al. 2017;Baglio et al. 2018), and might contribute up to the X-ray band via synchrotron-self-Compton emission (e.g.Markoff et al. 2005).With few confusing radiative processes and many sensitive observatories, the radio band is particularly suitable for jet studies.
Accreting black hole systems in their hard spectral state show a correlation between their X-ray and radio luminosity, that holds over orders of magnitudes in black hole mass from X-ray binaries to AGN, after incorporating a mass normalization term (the fundamental plane of black hole activity) and indicates a coupling between the inflow and outflow of matter (Hannikainen et al. 1998;Corbel et al. 2000Corbel et al. , 2003;;Falcke et al. 2004;Merloni et al. 2003;Gallo et al. 2003;Plotkin et al. 2013;Gallo et al. 2014).While this sample of stellar-mass black holes is dominated by binaries with a low-mass ( 1M ) donor -the low-mass X-ray binaries (LMXBs) -it includes two high mass X-ray binaries (HMXBs) hosting black holes (Cyg X-1 and MWC 656; the candidate black hole HMXB Cyg X-3 is not included).Together, these sources fol-low as similar correlation between X-ray and radio luminosity, suggesting that this coupling might be independent of mass transfer / donor type for black hole systems (Ribó et al. 2017).Looking more closely at the X-ray -radio correlation for stellar-mass black holes, there is evidence for both a radio-loud and radio-quiet track (Soleri & Fender 2011;Dinc ¸er et al. 2014;Meyer-Hofmeister & Meyer 2014;Drappeau et al. 2015).Despite several possible explanations, including inclination (Motta et al. 2018), variable jet Lorentz factors (Soleri & Fender 2011;Russell et al. 2015), and X-ray (Koljonen & Russell 2019) and radio (Espinasse & Fender 2018) spectral shapes, these tracks remain not fully understood.Moreover, the statistical evidence of their existence remains debated (Gallo et al. 2018).The correlation between X-ray and radio luminosity for black hole LMXBs disappears as the system transitions via the intermediate states into the soft state: during this transition, the compact jet quenches while fast ejecta can be launched (Fender et al. 2004;Russell et al. 2020); any remaining radio emission during the soft state is typically attributed to those ejecta, either unresolved or tracked as they move away from the LMXB (e.g.Bright et al. 2020).

A brief history of neutron star jet observations
The story is more complicated for neutron stars.Weaklymagnetized neutron stars accreting above ∼ 1% of the Eddington luminosity (L Edd ) can be divided into two classes based on their tracks in the X-ray color-color diagram: the Z and atoll sources (Hasinger & van der Klis 1989).The difference between these source classes and their various sub-classes is thought to be driven by instantaneous mass accretion rate (Lin et al. 2009;Homan et al. 2010).Z sources, accreting around the Eddington limit, are the radio brightest class of accreting neutron stars and their jets were therefore characterised first (Ables 1969;Lampton et al. 1971;Penninx et al. 1988;Penninx 1989;Hjellming et al. 1990b,a).These studies found different jet types along the different branches in the X-ray color-color diagram, changing in tandem with changes in accretion flow properties: from a compact jet to discrete ejecta and finally quenching at the highest mass accretion rates (qualitatively similar to the black hole behaviour detailed above; Migliari & Fender 2006).Jet studies of the X-ray and radio fainter atolls came later, with Migliari et al. (2003) presenting the first multiepoch X-ray and radio campaign for such a source (4U 1728-34).
Using the enhanced sensitivity of current day radio telescopes, neutron star LMXBs have now been studied extensively down to LX ∼ 10 36 erg/s (i.e.∼ 1% of the Eddington limit for a 1.4M neutron star).Observations at lower X-ray luminosities are dominated by radio non-detections (Tudor et al. 2017;Gallo et al. 2018;Gusinskaia et al. 2020a), although some neutron stars have been detected down in this regime as well (e.g.SAX J1808.4-3658 and IGR J00291+5934, down to ∼ 2 × 10 34 erg/s; Tudor et al. 2017).The behaviour of neutron star jets at low mass acccretion rates remains poorly explored.Another open question, for jets in atolls, regards the presence and mechanism of the jet quenching seen in black hole systems (Fender & Muñoz-Darias 2016).Atolls can change (although not all do) between thermaland Comptonisation-dominated accretion flow states: their soft and hard states, respectively.Several sources show a quenched jet in their soft spectral state (Migliari et al. 2003;Miller-Jones et al. 2010;Gusinskaia et al. 2017;Díaz Trigo et al. 2018;Gusinskaia et al. 2020a), as seen in black holes (e.g.Fender et al. 2004).Others, however, do not (Rutledge et al. 1998;Migliari et al. 2004).Finally, coordinated X-ray and radio studies have often focused on transient neutron star LMXBs, in order to probe different accretion rates -leaving persistently accreting sources more poorly explored.
The first comprehensive investigations of the radio properties of accreting neutron stars were presented by Fender & Hendry (2000) and Migliari & Fender (2006).Since then, many individual accreting neutron stars have been added to the X-ray -radio luminosity plane; see the compilations by Tetarenko et al. (2016) and Gallo et al. (2018), and, for the most up to date database, Bahramian et al. (2018) 1 .All of these sources -both atoll and Z -contain weakly-magnetised neutron stars (e.g.B < 10 10 G), where the Blandford & Payne (1982) mechanism can be applied.As a sample, these weakly-magnetised neutron stars are radio-fainter by a factor ∼ 20 than stellar-mass black holes -a difference that cannot simply be accounted for by the difference in accretor mass, bolometric X-ray corrections, or the presence of a boundary layer around the neutron star (Fender & Kuulkers 2001;Gallo et al. 2018).In the statistical analysis by Gallo et al. (2018), these sources show a scatter similar to the black hole population, assuming the latter follow a single track.
The key difference between accreting black holes and neutron stars -with respect to jet formation -is the presence of a solid stellar surface in the latter.This comes with additional differences, such as anchored magnetic fields and a compact object spin that can be measured via pulsations for neutron stars.As shown by Xray pulsations detected in a subset of accreting neutron stars, their magnetic fields can dynamically alter the geometry of the inner accretion flow, where jets are launched from.The neutron star magnetic fields can be measured directly through the detection of the cyclotron resonance scattering feature (hereafter cyclotron lines; see Staubert et al. 2019, for a recent review).Alternatively, and more indirectly, one can measure the inner disc radius and use that to constrain the magnetic field (Cackett et al. 2008;Degenaar et al. 2017;Ludlam et al. 2019), or constrain the magnetic field strength from the relation between the spin evolution and mass accretion rate (e.g.Ghosh & Lamb 1978;Campana et al. 2002;Strohmayer et al. 2018b).Neutron star spins are measured directly through Xray pulsations or nearly-coherent oscillations during thermonuclear bursts on their surface (Patruno & Watts 2012;Staubert et al. 2019).
Despite decades of radio studies, jets were until recently not detected in strongly-magnetised neutron stars.A number of sample studies between the 1970s and 2000s reported only non-detections (Duldig et al. 1979;Nelson & Spencer 1988;Fender & Hendry 2000;Migliari & Fender 2006;Migliari et al. 2011b), while the radio detection of the HMXB X-ray pulsar GX 301-2 by Pestalozzi et al. (2009) could be attributed to the radio emission from the stellar wind.The young neutron star X-ray binary, Cir X-1 (Heinz et al. 2013), is known to launch strong jets (Stewart et al. 1993;Fender et al. 1998;Tudose et al. 2006;Heinz et al. 2007;Soleri et al. 2009;Coriat et al. 2019).However, despite claims of a strong magnetic field (see, e.g.Schulz et al. 2020), its field strength has not been measured directly.The series of radio non-detections for stronglymagnetised neutron stars also formed the basis for the theoretical reasoning from Massi & Kaufman Bernadó (2008), explaining why the Blandford & Payne (1982) mechanism should not operate in this magnetic field regime.
Recently, the first jet from a strongly-magnetised neutron star was observed, contrary to this theoretical expectation (van den Eijnden et al. 2018a).The slow (spin period exceeding 1 second) X-ray pulsar Swift J0243.6+6124, that accretes from a Be star (Kouroubatzakis et al. 2017), launched a jet both during its so-called giant outburst in 2017/2018, and during X-ray rebrightenings in the outburst decay (van den Eijnden et al. 2019).While this jet is not expected to be launched via the Blandford & Payne (1982) mechanism, it remains unknown what alternative process could be responsible.Additionally, the strongly-magnetised neutron stars Her X-1 and GX 1+4 were detected in radio (van den Eijnden et al. 2018b,c), although the origin of this emission was not conclusively attributed to a jet.All three sources were detected at faint radio flux densities below ∼ 100 µJy, which explains the radio non-detections of this source class in earlier decades.

An extended parameter space for neutron star jets
The inclusion of strongly-magnetised neutron stars into the class of jet-launching sources, greatly expands the parameter space to study (neutron star) jets.Firstly, in addition to the larger range in neutron star magnetic field, a much greater spin range can now be accessed: while weakly-magnetised neutron stars show spins in the millisecond range (Patruno et al. 2017), their strongly-magnetised counterparts reach spins up to thousands of seconds (Staubert et al. 2019).Secondly, all confirmed HMXB neutron stars are, when measured, strongly magnetised2 , with typical magnetic fields of the order > 10 12 G.Therefore, these systems probe a much wider range of donor types, binary periods, eccentricities, and mass transfer mechanisms, than is accessible through only weakly-magnetised neutron stars.As the vast majority of HMXBs contains a neutron star, probing the effect of these binary and donor properties on jet formation was also barely possible with black hole systems.
At the same time, the origins of radio emission in HMXBs can be more complicated to untangle.Unresolved radio emission from a stellar-mass black hole or weakly-magnetised neutron star with a low mass donor is often automatically assumed to originate from a jet.In systems with a high-mass donor, the donor's stellar wind can also contribute to the radio emission.Studies of isolated O/B giants and Be-stars have shown them to be radio emitters (Lamers 1998a,b;Güdel 2002), implying that the companion's wind cannot be ignored when interpreting the radio emission from HMXBs: depending on properties such as mass-loss rate, clumpiness, and density, it can show either thermal (Bremsstrahlung) or non-thermal (shock) radio emission (Wright & Barlow 1975;Blomme & Runacres 1997;Dougherty & Williams 2000).These processes can lead to radio luminosities up to several times 10 27 erg/s in C-band, where jets are often studied, depending on the exact wind properties (see Section 4.2.1 for more details and discussion).Estimating the required wind properties to explain a HMXB radio detection, and having (multiple) coordinated radio and X-ray observations, preferably with spectral or polarization constraints, are therefore important tools to distinguish the stellar winds from the jet.
While a larger parameter space for jet studies is now accessible in terms of both neutron star magnetic field strength and spin frequency, it remains unclear what model can explain the existence of jets such as observed in Swift J0243.6+6124.Possible models, including jets powered by opened neutron star field lines (Parfrey et al. 2016) or magnetic propellers (Romanova et al. 2009), typically predict a dependence of jet power on magnetic field and spin frequency.Searches for such dependencies have been performed earlier for samples of weakly-magnetised neutron stars (Migliari et al. 2011b), or samples including both neutron stars and black holes (King et al. 2013).However, these studies were typically inconclusive or unable to find strong evidence for these relations, partially due to the limited spin range covered by weakly-magnetised neutron stars (Patruno et al. 2017).With a sample of neutron stars including a sufficient number of strongly-magnetised, slowly spinning sources, such studies can be revisited and extended -assuming, fundamentally, a common jet launching mechanism across the entire sample.
In this paper, we present a detailed study of a large set of new radio and X-ray observations of accreting neutron stars.Our sample includes 13 weakly-magnetised and 23 strongly-magnetised targets, more than doubling the total number of neutron stars in the radio/X-ray-luminosity plane observed at current radio sensitivities.Preliminary results from the weakly-magnetised sample were reported in Gallo et al. (2018, which included, in total, 41 neutron stars), while here we provide the full data analysis and the most up-to-date results.With this large set of new observations, we present the first systematic study of differences between the jets of weakly-and strongly-magnetised neutron stars where sources from both categories are detected, and aim to observationally constrain the neutron star jet formation mechanism.

OBSERVATIONS AND DATA ANALYSIS
In this section, we will first introduce the sample of observed neutron star X-ray binaries and the observing campaigns in radio and X-rays.In Sections 2.2 and 2.3, we will give a general introduction into the reduction and analysis of the radio and X-ray observations, respectively.Given the large number of analysed sources and the wide variety in setup of both the radio and X-ray campaigns, we discuss all details per source in the Online Supplementary Materials.

Targets and observing campaigns
The sample of targets presented in this paper consists of 13 weakmagnetic field neutron stars and 23 strong-magnetic field neutron stars, where we define weak and strong magnetic fields as B 10 10 G and B > 10 10 G, respectively.For most strongmagnetic field sources, the detection of a cyclotron line provides direct, robust magnetic field measurements (e.g.Staubert et al. 2019).For the remaining targets in this category, from a combination of the neutron star spin (evolution) and comparisons with slow pulsars with measured field strengths, their strong magnetic fields have been inferred (although therefore no actual measurement is performed).For weakly-magnetised accreting neutron stars, on the other hand, direct magnetic field measurements are not available.Instead, indirect measurements based on reflection spectroscopy (Cackett et al. 2008;Degenaar et al. 2017;van den Eijnden et al. 2018d;Ludlam et al. 2019), modeling of the evolution of the spin frequency (Patruno 2012), and magnetic propeller states (e.g.Mukherjee et al. 2015) imply typical magnetic fields between ∼ 10 7 and ∼ 10 9 G.In Table 2 of the Online Supplementary Materials, we list magnetic field measurements and estimates for all sources, alongside their spin and orbital period, when known.We want to stress again that these B-field measurements can be rather uncertain and affected by systematic effects -hence our broad classification in weak and strong magnetic fields.
The observations of all but one of the weak-magnetic field sources -IGR J17379-3747 -were already included in the statistical analysis of neutron star LMXBs by Gallo et al. (2018).However, that work only focused on the overall sample properties without focusing on individual systems.Moreover, it did not include details on the radio and X-ray data analysis, which are presented in the current paper.To the comparison literature sample, we have also added the recently published observations of two weakly-magnetised neutron stars: IGR J16597-3704 (Tetarenko et al. 2018a) and IGR J17591-2342 (Russell et al. 2018;Gusinskaia et al. 2020b), for the latter assuming the 7.6 kpc distance reported by Kuiper et al. (2020).Given the novelty of radio detections of the strongly-magnetised neutron stars, in the studied sample we include the few recently published observations of such sources: the radio detections of GX 1+4 (van den Eijnden et al. 2018b) and Her X-1 (van den Eijnden et al. 2018c), and the multi-epoch monitoring of Swift J0243.6+6124 (van den Eijnden et al. 2019).
An overview of all sources, divided based on neutron star magnetic field, is shown in Tables 1 and 2. These sources cover different, and overlapping, source classes: the thirteen weaklymagnetised neutron stars include (i) atolls; (ii) accreting millisecond X-ray pulsars (AMXPs), hosting a neutron star whose X-ray pulsations reveal a spin of several hundreds of Hz (Patruno et al. 2017); (iii) ultra-compact X-ray binaries (UCXBs), systems with an orbital period less than one hour; and (iv) very-faint X-ray binaries (VFXBs), where the neutron star persistently emits below ∼ 10 36 erg/s (Wijnands & Degenaar 2016).In the VFXB category, we include the quasi-persistent source XMMU J174716.1-281048, that was likely continuously active between 2003 and 2011 (see e.g.Degenaar et al. 2011) but was not detected in X-rays during the 2014 radio observations reported in this work.These four source classes combined tackle different, relatively unexplored science cases: the jet properties of neutron star LMXBs that are persistently accreting (i and iii), are in the soft state (i), or are X-ray faint (ii, iii, and iv).We note, finally, that our source sample does not include any Z-sources: given their radio brightness, these sources have been studied in detail previously, and this work's approach of single or a small numbers of observations per source, would not significantly contribute to their understanding.Hence, we decided not to perform observations of such targets and, lacking new observational results, will not further discuss this source class in much detail in this work.
The strongly-magnetised neutron stars (Table 2) include one candidate and two confirmed symbiotic X-ray binaries, where the neutron star accretes from the stellar wind of an evolved low-mass donor in a wide orbit; three LMXBs (of which one is a UCXB); one intermediate-mass X-ray binary (IMXB); and sixteen HMXBs.Based on the donor type, the latter are categorised either as Be/Xray binaries (BeXRBs), or as Super-giant X-ray binaries (SgXBs).For an overview of the differences between these source classes, see Reig (2011).Finally, 3A 1239-599 is simply denoted as HMXB as it is unknown in what subcategory it falls.The sources in this class were mainly targeted to study the poorly-understood jet properties of strongly-magnetised neutron stars and explore the effect of other radio emission mechanisms, such as their stellar winds.The four new BeXRBs (i.e.all but Swift J0243.6+6124) were targeted specifically to probe their radio properties at very low accretion rates, close to or in their propeller regimes.
The latter class also includes GRO J1744-28, known as the Bursting Pulsar: a LMXB where the neutron star has an intermedi-ate spin frequency of 2.1 Hz (Cui 1997), causing it to fall somewhat between the two source categories.Its magnetic field is likely lower than typically seen in strongly-magnetised systems (e.g.B 10 12 G): using different methods, it has been claimed to lie between 2 × 10 10 and 7 × 10 11 G (Cui 1997;Rappaport & Joss 1997;Bildsten et al. 1997;Degenaar et al. 2014;Younes et al. 2015).For this paper, we include it in the strong-magnetic field class, although ultimately, this classification does not affect our conclusions significantly: as shown by the preliminary results published by Russell et al. (2017), the Bursting Pulsar is not detected at radio frequencies, with a relatively unconstraining upper limit due its proximity to the Galactic centre.
Where possible, we used parallax measurements from Gaia Early Data Release 3 (Gaia Collaboration et al. 2020) to constraint the distance.This was done for sources where (i) the detected Gaia counterpart has (ii) a positive parallax measurement with (iii) a signal-to-noise-ratio larger than three.To assess the affect of choice of prior in converting the parallax into distance, we then compared inferred distances from the priors in Atri et al. (2019) (for Galactic LMXBs), Bailer-Jones et al. (2018), andBailer-Jones et al. (2020).We found these to be consistent within errors and use the Atri et al. (2019) prior in the analysis.For sources where Gaia was not used, we searched the literature for distance measurements instead.Finally, we note that using literature (non-Gaia) distances for all sources does not alter the main conclusions of this work.

Overview of radio data analysis
The radio observations of our neutron star sample were performed with the Karl G. Jansky Very Large Array (hereafter VLA) for most sources with declinations above ∼ −40 • , and the Australia Telescope Compact Array (ATCA) for the remaining Southern targets (negative declinations).ATCA was in the most extended 6-km configurations during all campaigns, while the VLA changed configuration between observations.All raw data sets are publicly available under VLA programme codes SD0134,, and ATCA programmes C3108, C3184, C3243, and CX379 (see table 1 in the Online Supplementary Materials).
To flag, calibrate, and image the observations, we used the COMMON ASTRONOMY SOFTWARE APPLICATION (CASA; Mc-Mullin et al. 2007) package v4.7.2.We removed RFI using a combination of automatic flagging routines and visual inspection.We performed imaging using the multi-scale multi-frequency CLEAN task, with a Briggs robust parameter adjusted to the target field in order to balance sensitivity and confusion.We then fitted an elliptical Gaussian with Full-Width Half Maxima equal to the synthesized beam size of the observation using the CASA-task IMFIT.We measured the root-mean-square variability over a nearby region devoid of sources in case the target was detected, or over the target region for a non-detection.In the latter case, we set the 3σ upper limit to three times this RMS measurement.Target-specific details, such as VLA configuration, primary and secondary calibrator, and beamsizes, are listed per source in Table 1 of the Online Supplementary Materials.
For a subset ofd neutron stars, data was recorded at two frequencies.In those cases, we calculated the radio spectral index α, where Sν ∝ ν α , between the two bands.The error on α is estimated through a propagation of the uncertainties on the radio flux densities measured at each frequency, using Monte-Carlo simulations.In case radio emission is detected in only the lower frequency band (2 sources), we follow the Monte-Carlo approach of van den Eijnden et al. (2019) to estimate an upper limit on the spectral index.For those sources, we show the diagnostic figures of this method in Section 4 of the Online Supplementary Materials.In this work, we will refer to negative and positive spectral indices as steep and inverted spectra, respectively.

Overview of X-ray data analysis
We measured unabsorbed X-ray fluxes of our targets using pointed observations of the Neil Gehrels Swift Observatory (Swift Gehrels et al. 2004) or monitoring observations with the Monitor of All-Sky X-ray Image/Gas Slit Camera (MAXI/GSC Matsuoka et al. 2009).Details about the X-ray analysis per source, such as observations used and spectral fit parameters, can be found in the second section of the Online Supplementary Materials.We aimed to use observations taken on the same day as the radio observation; in the cases where such observations did not exist, we used the closest Xray observation in time.In those cases, we used longer-term X-ray monitoring, combined with typical time-scales of state changes, to ensure that the source did not change its state or X-ray flux significantly around the radio observation.We preferentially used pointed Swift X-ray Telescope observations, either measuring the flux directly from the spectrum, or for very faint sources converting the count rate or count rate upper limit into the flux.These Swift analyses where permormed in the 0.5-10 keV range.When no pointed observations were available, we measured the flux from the (multiday) MAXI spectrum (fitted between 2 and 10 keV) or converted the detected MAXI 2-10 keV count rate into a flux.For all sources that are systematically undetected in MAXI, we ensured that Swift observations were performed.
In order to convert either Swift or MAXI count rates into unabsorbed fluxes, we used the WEBPIMMS tool3 .When the source was in a state with a typical and known X-ray spectrum, we used the literature to model the spectrum used in the count rate conversion.Otherwise, we used the Crab conversion following the NuS-TAR measurement of the Crab spectrum by Madsen et al. (2017).Note that, whether a spectral model was assumed/fitted (see below) or the Crab spectrum was used, we use the full flux and do not attempt to distinguish between different spectral components; see (Miller et al. 2012) for a study where such effects are taking into account.
All measured fluxes were calculated in the 0.5-10 keV range; hence, we note that we extrapolated the model fitted to the MAXI spectra down to lower energies, thereby possibly introducing extra uncertainy in the flux measurement due to interstellar absorption.While differences exist in the shape of the X-ray spectrum between LMXBs and HMXBs, we used the same energy band to enable direct comparison between the source classes and with the literature.We fitted X-ray spectra using XSPEC v.12.10.1 (Arnaud 1996), setting the ISM abundances and cross-sections to Wilms et al. (2000) and Verner et al. (1996), respectively.As some spectra contain few counts, we used C-statistics to find the best fit (Cash 1979).All spectra were modelled with three models, combining interstellar absorption (TBABS) with power law or/and blackbody models: TBABS*POWERLAW, TBABS*BBODYRAD, and TBABS*(POWERLAW + BBODYRAD).We picked the best-fitting model of the former two, based on the lowest test statistic, and then compared this with the latter, combined model using a f-test.We selected the combined model if the f-test preferred it at a > 5σ probability over the best-fitting single-component model.Finally, we measured the flux and its uncertainty by convoluting the selected model with CFLUX.While this approach is phenomenological, it suffices to measure the X-ray flux even for observations with low numbers of total counts.
For two radio observations, no X-ray information from either monitoring or pointed observations was available sufficiently close in time, compared to the source's typical time scale of variability and state transitions.These observations were the second radio observation of GX 1+4 and the second radio observation of 4U 1954+31.The former was detected in radio during this epoch, while the latter was not.Given the lack of X-ray data, we do not include these observations in the X-ray -radio luminosity diagrams later in this paper.

RESULTS
We list all sources, radio flux densities and X-ray fluxes in Tables 1 (the 13 weakly-magnetised sources) and 2 (the 23 stronglymagnetised sources).Out of the weakly-magnetised category, a radio counterpart is detected for four targets for the first time: the persistent atolls GX 3+1, GS 1826-24, and 4U 1702-429, and the AMXP IGR J17379-3747.From the 23 strongly-magnetised neutron stars, nine sources are detected: the symbiotic X-ray binaries GX 1+4 and 4U 1954+31, the intermediate-mass X-ray binary Her X-1, the Supergiant X-ray binaries 1E 1145.1-6141,4U 1700-37, Vela X-1, IGR J16318-4848 and IGR J16320-4751, and the Be/Xray binary Swift J0243.6+6124.The radio detections of the latter and Her X-1 were already reported (van den Eijnden et al. 2018c,a) but we include these as they were not compared with a larger sample of strongly-magnetised neutron stars yet.GX 1+4 was also presented before (van den Eijnden et al. 2018b), but here we add three more detections in new observations at different observing frequencies.
In Figure 1, we show the X-ray -radio luminosity plane for black hole and neutron star X-ray binaries, in order to search for coupling between the X-ray emitting accretion flow and radioemitting jets.For now, we do not yet attempt to distinguish between sources where the radio emission is clearly attributable to a jet, or other processes that might contribute; see Section 4 for an extensive discussion on this topic.The newly added sources from our sample are shown per magnetic-field class as the blue (weak magnetic field) and red (strong magnetic field) data points.We originally compiled the comparison sample, shown with grey crosses for black holes and black circles for neutron stars, for the statistical study in Gallo et al. (2018), and complemented it by two neutron stars discovered since (see Section 2.1).We stress that this comparison sample does not include soft state atolls, while our weakly-magnetised sample does.We plot the 6-GHz radio luminosity, which we calculated by first estimating the 6-GHz flux density Sν using the spectral shape where known or otherwise assuming a flat radio spectrum.Then we calculated LR = 4πνSν D 2 , where ν = 6 GHz and D is the distance to the source.All X-ray luminosities LX are calculated in a similar fashion in the 0.5-10 keV range using LX = 4πFX D 2 , where FX is the measured, unabsorbed X-ray flux.
We briefly note that in the X-ray -radio luminosity plane, we only plot statistical errors on both luminosities.Underlying assumptions, for instance a flat spectrum for single-frequency observations, and issues such as non-simultaneity of observations, uncertainties on distances, and errors in absolute flux calibration, cause systematic uncertainties in these comparisons.We do not include those in Figures 1 and 2 as the literature sample similarly uses statistical errors only and these systematics are challenging to constrain accurately.However, one should keep their existence in mind when interpreting X-ray -radio luminosity diagrams.
From Figure 1, it is immediately apparent that only few neutron stars are detected in the radio band below an X-ray luminosity of LX ≈ 5×10 34 erg/s.Conversely, above LX ≈ 10 36 erg/s, most of the radio observations of accreting neutron stars yield a detection at current sensitivities.As discussed in detail in Section 5.1, none of the strongly-magnetised neutron stars reach above LR ≈ 3 × 10 28 erg/s, independent of the X-ray luminosity.This faint apparent maximum radio luminosity -corresponding to ∼ 170 µJy for a typical distance of 5 kpc -explains why these sources remained undetected in previous observing campaigns at lower radio sensitivity (Duldig et al. 1979;Nelson & Spencer 1988;Fender & Hendry 2000).
In the sample of sources in Figure 2, two subtle points might be easily missed.Therefore, we point them out explicitly here: firstly, two strongly-magnetised neutron stars, the SyXRBs 4U 1954+31 and the HMXB IGR J16318-4848, were only detected in radio and not in X-rays.Secondly, we plot the detected X-ray Apart from GX 1+4 and Her X-1, we do not show previous observations of strongly-magnetised neutron stars, as these are all upper limits at typically unconstraining levels (although we discuss several exceptions in Section 5.3.4).Radio luminosities have been converted to 6 GHz, using a measured spectral index if available, or assuming a flat spectrum otherwise.
a Kolmogorov-Smirnov test comparing the strongly-and weaklymagnetized neutron star samples returns a p-value of p = 7 −11 for the hypothesis that both are drawn from the same underlying distribution.Alternatively, the Anderson-Darling test consistently finds p < 10 −3 (we note that this limit is due to the implementation of this test in PYTHON/SCIPY: the measured test statistic greatly exceeds the critical value corresponding to p = 10 −3 ).
We note that a bias can be introduced by differences in the X-ray luminosity distribution, as the three different source classes might be dominated by different X-ray luminosities, translating in different dominant radio luminosities if a coupling exists between both (see e.g.Gallo et al. 2018, for a detailed discussion).To assess this effect, we also show the X-ray luminosity histograms of each source class in the right panel.While there are differences in the distributions, these are minor compared to the differences in radio luminosity.In fact, the only major difference is the peak at super-Eddington X-ray luminosities for strongly-magnetised sources, attributable completely to Swift J0243.6+6124.However, this difference only emphasises the striking radio faintness of strongly-magnetised neutron stars: the super-Eddington peak in X-rays should shift the radio luminosity distribution of strongly-magnetised sources to higher values, if these sources would show a black hole-like X-ray-radio coupling (without an apparent maximum radio luminosity).In other words, the strongly-magnetic neutron stars are radio faint, despite the domi-nance of the X-ray bright Swift J0243.6+6124 in the HMXB sample.

Comparing low-mass and high-mass X-ray binaries
Radio emission observed in Roche-lobe overflowing LMXBswhether the primary is a black hole or a weakly-magnetised neutron star -is typically attributed to synchrotron processes in a relativistic jet (e.g.Corbel et al. 2000;Dhawan et al. 2000;Stirling et al. 2001;Fender et al. 2004;Migliari & Fender 2006;Gallo et al. 2018).
Our samples of weakly and strongly-magnetised neutron stars contain 16 LMXBs: all weakly-magnetised neutron stars, plus the slow pulsars 2A 1822-371, 4U 1626-67, and GRO J1744-28 (here, we ignore the SyXRBs, which have a low-mass donor but accrete from the stellar wind; see Section 4.4).Out of these sixteen, four sources are detected (e.g.Table 1): the Atolls GX 3+1, GS 1826-24, and 4U 1702-429, and the AMXP IGR J17379-3747.Their radio spectral shapes and positions on the X-ray -radio luminosity diagram fit with a jet identification for accreting neutron stars, being consistent with the larger sample of sources in the literature (Russell et al. 2013b;Gallo et al. 2018).In the undetected sources, the nondetections can typically be attributed to either their spectral state   Figure 3. Histograms of the radio (left) and X-ray (right) luminosity for the full (archival plus this work) samples of black hole, weakly-magnetised, and strongly-magnetised X-ray binaries.Note that we plot all observations, including when multiple observations of the same source are available.Coloring is black, blue, and red, respectively, following Figure 1.We only show detected sources and do not take upper limits into account.The histograms are normalised by their peak values to ease comparison of their shapes.The dashed line in the right panel indicates the Eddington limit for a 1.4 M neutron star.
Table 1.Radio and X-ray observations of 13 weakly-magnetised accreting neutron stars.For source type acronyms, see Section 2.1 -hard and soft for atolls refer to their state during the radio observation.*IGR J17062-6143 is also an AMXP and VFXB.Upper limits in radio are listed at the 3σ level, while X-ray upper limits are shown at the 90% level.All other uncertainties are 1σ.The X-ray fluxes are unabsorbed fluxes in the 0.5-10 keV band.Distance references: a poorly constrained, see Savolainen et al. (2009)  or faint X-ray luminosity.For a detailed discussion, we refer the reader to Section 5.3, where we will comment on individual (detected and undetected) sources.
The identification of the radio emission origin is less straightforward for high-mass and symbiotic X-ray binaries.As alluded to in the introduction, these sources are more complicated for two reasons: firstly, their donors launch stellar winds that could contribute to the detected radio flux via various mechanisms.Secondly, the jet properties of strongly-magnetised neutron stars are more poorly explored and therefore existing literature offers few comparison studies.In the remainder of this section, we will review six options for the origin of the radio emission: emission from the donor star itself, emission from stellar winds and their interaction with other components of the binary system, coherent processes, emission from a slow, wide-open outflow from a propeller-type mechanism, and relativistic jets.
Three options can be excluded directly.While stars do emit in the radio band, such emission is typically seen in coronally-active low-mass stars (type F or later).In addition, these stars reach maximum X-ray luminosities of ∼ 10 32 erg/s, orders of magnitude below the X-ray luminosities of our radio-detected targets (Guedel & Benz 1993;Güdel 2002;Kurapati et al. 2017).Secondly, a wideopen, slow gas outflow driven by a magnetic propeller-type mechanism, as for instance seen in simulations of accretion onto magnetised (neutron) stars by Romanova et al. (2009) and Parfrey et al. (2017), is similarly unlikely: as discussed also in Section 5.2, none of our radio-detected targets resided in the propeller regime during the observations.Thirdly, coherent processes have been inferred the radio emission of types of accretion magnetic white dwarfs (e.g.Barrett et al. 2017, see also Section 4.3).These processes, namely an electroncyclotron maser or gyrosynchrotron emission, are unlikely to operate in the systems considered here: firstly, electron-cyclotron maser emission is associated with high (∼ 100%) levels of circular polarization.While the ATCA observations studied here were not set up to measure circular polarization, the earlier VLA studies of Her X-1, GX 1+4 and Sw J0243, did not show any circular polarization (e.g.van den Eijnden et al. 2018b).Secondly, gyrosynchrotron emission is expected to show lower levels of polarization.Such emission would originate from the magnetosphere of the neutron star.However, in actively accreting HMXBs, the magnetospheric radius does not, for reasonable neutron star parameters and accretion rate, exceed ∼ 10 4 Rg (Tsygankov et al. 2017).In comparison, the typical minimum emission size of the detected radio emission is roughly two orders of magnitude higher (see Section 4.3 and equation 3).
4.2 Winds from massive donor stars and their interactions?

Stellar winds
Stellar winds and their emission properties have been studied extensively for both single and binary stars in the past decades (see Negative spectral indices have been measured in several stellar winds from high-mass stars, however.In those cases, it is thought that these systems are (massive) binary systems, and the shock between the two stellar winds (that does not move outwards as it would for a single star) is observed (Moran et al. 1989;Churchwell et al. 1992;Dougherty et al. 1996;Dougherty & Williams 2000;Williams et al. 1997;Chapman et al. 1999;Contreras et al. 1997;Ortiz-León et al. 2011;Blomme & Volpi 2014;Blomme et al. 2017).At the most extreme end, where both stars are early-type / Wolf-Rayet stars with strong winds, these sources are observed as colliding wind binaries, which are known sources of X-ray and non-thermal radio emission (e.g.Dubus 2013).
For thermal wind emission, the expected radio flux densities can be estimated using the formalism derived by Wright & Barlow (1975, note how similar results were derived around the same time by Olnon 1975 andPanagia &Felli 1975).Based on the observing frequency ν, electron temperature Te, wind mass loss rate, the mean atomic weight per electron µe, the terminal wind velocity v∞, and the distance d, we can estimate the flux density as: Winds from massive stars are known to be clumped, which affects their observational appearance.However, the clumpiness decreases with distance from the star, and hence, the radio emission is the least affected by such effects.As a result, we can apply the above formalism, which was formulated for smooth plasmas (Puls et al. 2006).The brightness temperature of these thermal winds is typically 10 4 -10 5 K (Longair 1992); however, even with VLBI-like resolution, none of the targeted HMXBs are bright enough to reject the thermal wind hypothesis based on a minimum brightness temperature argument (see also Section 6).For the non-thermal, shocked emission, the brightness is more difficult to predict, but observational constraints exist.
Such observational constraints on the radio emission of stellar winds (in single and binary stars) have been obtained through many studies.For our comparison, the most relevant are the studies of OB supergiants.The most constraining and sensitive study of this kind was recently performed with ATCA and ALMA, observing the stellar cluster Westerlund 1 (Fenech et al. 2018;Andrews et al. 2019).Located at approximately 5 kpc (Clark et al. 2019), Westerlund 1 hosts at least 100 OB supergiants (Negueruela et al. 2010).However, at current sensitivities, only seven of those stars (i.e. at most a few percent) are detected at GHz frequencies.Combining the ATCA and ALMA observations reveals that all radio detected OB supergiants have steep spectra, indicating that these are likely binary star systems (Andrews et al. 2019).Therefore, Andrews et al. ( 2019) state that at current sensitivities, they 'would not expect to detect any emission from purely thermal stellar wind emitters in radio'.Importantly, these results imply that the radio emission from single OB supergiants is difficult to detect and might be overpredicted by, for instance, the Wright & Barlow (1975) formalism.
The studies of Westerlund 1, and in fact all radio studies on massive stars in the literature, focus on single and binary massive, nondegenerate stars.Similarly, the Wright & Barlow (1975) model was not developed for massive stars in X-ray binaries, but for single massive stars.Therefore, this model does not include the effect of X-rays emitted by the accretion flow on the stellar wind properties,   and its resulting radio luminosity.However, the picture painted by the comparison between Westerlund 1 and the model predictions, fits with our results on OB supergiants with a neutron star companion.We show this graphically in Figure 4: we plot the measured radio luminosity (upper limits) as a function of the wind mass loss rate and velocity as parameterized in Equation 1, alongside their relation according to that equation.In our calculations, we ignore the clumpy nature of winds in HMXBs (Martínez-Núñez et al. 2017;Grinberg et al. 2015Grinberg et al. , 2017)), for the reasons discussed above.
As is apparent from Figure 4, four sources out of the nondetected HMXBs (Cen X-3, 3A 1239-599, OAO 1657-41, IGR J16207-5129) should have been detected given the ATCA sensitivity, stellar wind properties, distance, and equations for thermal wind radio emission.A fifth source, EXO 1722-363, would be above the detection limit, while the remaining non-detected SgXBs are too distant for the wind to be detected.The non-detection of the five named sources above supports the notion that OB supergiants are radio fainter than in the simple Wright & Barlow estimate and challenging to detect in radio, unless they are in stellar (not X-ray) binary systems.
One can consider a scenario wherein the apparent theoretical over-prediction of the radio flux density for the single stars in Westerlund 1 would result from incomplete ionisation.However, given the temperatures of O/B stars, such an explanation appears unlikely.Moreover, in the X-ray binary estimates above, we assumed that the literature mass loss rate corresponds fully to ionised material.Comparing the radio non-detected with detected HMXBs in our sample, we do not find a systematic difference in X-ray luminosity or orbital period; therefore, we do not necessarily expect a systematic difference in the degree of X-ray ionisation of the wind.Hence, we expect that even if a partial ionisation of the wind could explain their radio non-detections, this possibility should also exist for the radio-detected systems discussed below.
In our sample, five Supergiant X-ray binaries are detected with ATCA.For two of those sources (1E 1145.1-6141,IGR J16318-4848), their donor (wind) properties predict a radio flux density substantially below the detected levels, assuming the wind is fully ionised (e.g Figure 4).Of those, only IGR J16318-4848 has a spectral shape that may be consistent with such stellar wind emission (α = 1.06 ± 0.35).4U 1700-37 is detected at slightly higher radio luminosity than predicted, although systematic uncertainties in this comparison might account for that difference.For the fourth source, IGR J16320-4751, no wind characteristics are known, preventing a flux density estimate.For Vela X-1, both the flux density estimates and spectral index fit with the Wright & Barlow (1975) description.Therefore, we conclude that for Vela X-1 and 4U 1700-37, the stellar wind might have a substantial contribution to the observed emission (again with the caveat that thermal winds appear difficult to detect for OB supergiants), while it likely contributes less in the other detected sources.All details on these estimates can be found in the third section of the Online Supplementary Materials.

Intrabinary wind shocks
Could the stellar wind interact with the pulsar wind, causing the observed radio emission?It is commonly assumed that, in accreting systems, the radio pulsar mechanism is suppressed as the magnetosphere is filled with ionised material, implying that no pulsar wind is launched.However, with that in mind, we can briefly compare the radio properties of the strongly-magnetised sources in our sample with the known types of systems where the pulsar and stellar wind interact.For instance, shocks between the stellar and pulsar winds, similar to those in stellar binaries resulting in negative spectral indices in the wind radio emission, can occur in binary systems where no accretion takes place.In such systems, this shock creates non-thermal emission that dominates the spectrum from radio to gamma-rays (Dubus 2013).As the spectral energy distribution of these systems peaks above 1 MeV, they are referred to as γ-ray binaries.Only a small number of γ-ray binaries are known to date (Paredes & Bordas 2019;Corbet et al. 2019).None of our detected sources are γ-ray binary candidates; indeed, the typical radio luminosities of γ-ray binaries are ∼two orders of magnitude larger (see e.g.Dubus 2013) than the apparent maximum radio luminosity of neutron star HMXB of ∼ 3 × 10 28 erg/s we find in our results.
γ-ray binaries are one example of shock interaction between the pulsar wind and surrounding material; similar shocks can be seen in pulsar wind nebulae (PWNe), systems where the pulsar wind of an isolated neutron star interacts with surrounding supernova matter or the interstellar medium.The magnetic wind of relativistic electrons and positrons can carry away the majority of energy lost in the spin down of the pulsar (Rees & Gunn 1974;Michel 1982;Kennel & Coroniti 1984).The shock can be detected through radio synchrotron emission as the electrons gyrate around the magnetic field lines.If we again ignore, for the sake of the comparison, the common assumption that accreting neutron stars do not launch a pulsar wind, we can ask the following question: could a similar interaction be at play in HMXBs -creating a dialed-up version of a PWN, or a dialed-down version of the extreme γ-ray binariesto explain our radio detections?The most quantitative constraints follow from the radio spectral shape.PWNe typically show powerlaw radio spectra with an index α ≈ −0.3 (Weiler & Panagia 1978;Gaensler et al. 2000).The best spectral constraints in our sample of HMXBs (i.e.Vela X-1, 4U 1700-37, IGR J16318-4848) show α > 0 instead.
The energetics of the shocks can also provide constraints.The radio luminosity of an individual system's shock is challenging to predict, as the spin down energy of the isolated pulsar population spans a wide range (e.g. 10 28 -10 39 erg/s; Gaensler et al. 2000), and the fraction transferred into radio emission can vary between PWNe.Furthermore, in accreting strongly-magnetised neutron stars, the pulse frequency evolution is regulated by the interactions with the accretion flow.Therefore, a measurement of the spin period and its derivative do not translate into a spin-down energy estimate, as it does for isolated pulsars.However, using magnetic field and spin measurements, one can estimate what the corresponding spin-down energy would be in the absence of accretion (ignoring the effects of accretion on the strength of the magnetic field trough, for instance, burial).Assuming a typical NS mass of 1.4 M and radius of 10 km, one can combine the spin down energy Ėspin = 4π 2 I Ṗ /P 3 and the estimate of the magnetic field (B/10 12 G) 3.2 × 10 7 (P Ṗ /1s) 1/2 to derive: For the sources in our sample, we find a wide range of values between ∼ 3 × 10 14 erg/s to 10 32 erg/s, with just six sources overlapping with the low end of the distribution of the PWNe population.More importantly, there are no systematic differences between radio-detected and non-detected targets, and for five out of the seven detected sources where both B and P are known, Ėspin is significantly lower than the observed radio luminosity.Therefore, the energetics argue against contributions from intrabinary shocks.Also, we again stress that this comparison only holds if accreting neutron stars do launch pulsar winds, which is not usually assumed.
We are left with four detected, strongly-magnetised sources that we have not yet discussed in this context; for Swift J0243.6+6124, a wind was excluded by van den Eijnden et al. (2018aEijnden et al. ( , 2019)), based on its luminosity, spectral change, and variations in those two properties throughout the source's giant outburst in 2017/2018.Secondly, Her X-1 has an intermediate-mass, Rochelobe overflowing donor, which does not launch strong stellar winds.Therefore, it is unlikely that the radio emission in this system is attributable to such a mechanism.Finally, we observed and detected two SyXRBs, where the neutron star accretes from the stellar wind of an evolved low-mass donor.Given the different wind properties of these sources compared to HMXBs, we separately review them in Section 4.4.

Relativistic jets from strongly-magnetized neutron stars?
Relativistic jets, as observed in weakly-magnetised neutron star and black hole X-ray binaries, emit radio emission through synchrotron processes.Given the broad range of observed spectral and brightness properties of jets in X-ray binaries, those observations pose few stringent constraints on what the emission can look like.However, we can make the comparison with these properties.Depending on the type of jet, this emission is either optically thick (α 0) or thin (α < 0, typically α ≈ −0.7; Fender et al. 2004;Russell et al. 2013b) for a compact, steady jet or discrete ejecta, respectively.In terms of the observed spectrum, the detected stronglymagnetised neutron stars fit within the expectation for a compact jet.However, with the wide range of indices possibly generated by jet synchrotron processes, this hardly implies the emission necessarily originates from a jet.
Secondly, we can consider the observed radio luminosity, given the X-ray luminosity.No unique prediction, from either theory or observations, exists for the radio luminosity based on source state, mass accretion rate, or X-ray luminosity.However, it is known observationally that weakly-magnetised neutron star Xray binaries are typically substantially radio fainter than accreting black holes (Fender & Kuulkers 2001;Migliari & Fender 2006;Gallo et al. 2018).All radio-detected, strongly-magnetised neutron stars are radio-fainter than the black hole tracks in the X-ray -radio luminosity plane, fitting with our expectations for relativistic jets.
Finally, we can consider the constraints from the Compton limit on the brightness temperature for synchrotron radiation of T = 10 12 K.The angular size of the emitting region θ can be calculated as (e.g.Longair 1992): For a typical flux density of 100 µJy at 5 GHz (i.e.∼ 6 cm), setting T < 10 12 K yields θ > 2 µas.At a typical distance of 5 kpc, this angular size sets a minimal physical size of the emitting region of > 10 11 cm, or > 7 × 10 5 Rg for a neutron star.Therefore, the observed flux densities are consistent with jet synchrotron emission, as the radio-emitting regions of the jet typically lie further out (i.e. 10 7 -10 9 Rg).
In this comparison with other source classes, we can also briefly compare the strongly-magnetised neutron stars with jetlaunching white dwarfs.On the one hand, these systems are far apart in their fundamental physical properties such as magnetic field, accretor size.However, they share an important property: a 10 12 G, 1.4 M neutron star, accreting from a disc at LX ≈ 10 36 erg/s, has a magnetospheric radius of Rm ∼ 6600 km (3.2 × 10 3 gravitational radii).This scale is comparable to the typical size of a white dwarf.Therefore, jet formation mechanisms at play in the accretion discs of accreting white dwarfs, might play a role beyond Rm in the disc of strongly-magnetised neutron stars as well.
What types of accreting white dwarfs launch jets?Similar to strongly-magnetised neutron stars, accreting white dwarfs were long thought not to launch jets (Livio 1997(Livio , 1999)).However, radio observations of SS Cyg in the past two decades reveal a jet launched by this famous accreting white dwarf, whose magnetic field is weak enough that the accretion flow extends to its surface (Körding et al. 2008;Russell et al. 2016;Fender et al. 2019).SS Cyg is an example of a non-magnetic Cataclysmic Variable (CV), a white dwarf accreting from a low-mass donor via Roche-lobe overflow4 .Subsequent observations of persistent (i.e.nova-likes) and transient CVs in outburst (i.e.dwarf novae) mostly resulted in radio detections as well, consistent with jet emission (Coppejans et al. 2015(Coppejans et al. , 2016)).However, in those systems, alternative mechanisms could not be ruled out as confidently as in SS Cyg.Therefore, while it appears that non-magnetic CVs are capable of launching jets, the increase in detected sources has introduced many remaining questions (Coppejans & Knigge 2020), a development that repeats in this study of strongly-magnetised neutron stars.
Comparing non-magnetic CVs with our sample more quantitatively, one finds that the radio luminosities of non-magnetic CVs lie significantly below those reported here for strongly-magnetised neutron stars: 4 × 10 26 erg/s (Russell et al. 2016;Coppejans et al. 2015Coppejans et al. , 2016;; Coppejans & Knigge 2020) versus ∼ 2 × 10 27 -3 × 10 28 erg/s, respectively.However, the former are detected at lower X-ray luminosities as well (LX < 10 33 erg/s), which are observable due to the smaller distances to the observed targets: the accretion flux is not dominated by the X-ray band to the extent that it is in X-ray binaries.To finish, we note that this comparison assumes the strongly-magnetised neutron star accretes from a disc -a separate comparison with the wind-accreting white dwarfs in symbiotic stars is made below (Section 4.4).
The recent radio campaigns on the 2017/2018 giant outburst of Swift J0243.6+6124 have demonstrated that strongly-magnetised neutron stars can launch jets (van den Eijnden et al. 2018aEijnden et al. , 2019)).Combined with the above issues with explaining the observed radio properties purely through stellar winds, we therefore conclude that it is possible that the radio emission observed from the neutron star HMXBs in our sample is dominated by synchrotron emission from a relativistic jet, and will assume so in the following section of the discussion.The exception to this interpretation are Vela X-1 and 4U 1700-37, which are the sources where the theoretical stellar wind predictions are similar to the observed radio emission.We will also assume such a jet origin for the emission in Her X-1, although we stress that the origin of its radio emission could not be unambiguously identified (van den Eijnden et al. 2018c).

The case of symbiotic X-ray binaries
Finally, we turn to the SyXRBs.All but one (Shaw et al. 2020) of known SyXRBs (nine confirmed and three candidate systems; e.g.Bahramian et al. 2014Bahramian et al. , 2017;;Qiu et al. 2017;Kennea et al. 2017;Bozzo et al. 2018) host a strongly-magnetised neutron star.Therefore, they are interesting analogues to (wide) wind-fed HMXBs.Within the wind capture radius, their accretion flows and possible jets could be similar, while the donor wind itself can be quite different.Particularly, the wind velocities of late-type giants are significantly lower, at typically ∼ 100 km/s, while their mass loss rates tend to be lower and span a large possible range of ∼ 10 −10 -10 −5 M /yr (Espey & Crowley 2008;Enoto et al. 2014).This combination of lower velocity, boosting possible thermal wind radio emission, and wide range in mass loss rate, complicates the comparison of the jet and wind scenario for the three SyXRBs in this study: GX 1+4, 4U 1954+31, and the candidate 3XMM J181923.7-170616.
We detect two out of the three SyXRBs systems in the radio, namely GX 1+4 and 4U 1954+31.The former is detected during four epochs at levels between ∼ 40-100 µJy, while the latter is detected in one of two observations at 21 ± 5 µJy.During its second observation, no X-ray information was available, although it is likely that 4U 1954+31 remained in the same faint X-ray state as during the first observation.In that scenario, the radio nondetection is likely due to a poorer sensitivity during the second observation (see also Table 1 in the Online Supplementary Materials).Finally, 3XMM J181923.7-170616 was not detected, with a 3-σ upper limit of 36 µJy at 5.5 GHz.
The wind velocity in GX 1+4 is inferred to be ∼ 100 km/s (Chakrabarty et al. 1997;Hinkle et al. 2006), while no direct measurements are available for the other two sources.Assuming all three have similar wind velocities, we can invert Equation 1 to estimate the required mass loss rates to explain their observed radio properties.This yields Ṁwind ≈ (1.4-2.7)× 10 −7 M /yr for GX 1+4, Ṁwind ≈ 1.1 × 10 −8 M /yr for 4U 1954+31, and Ṁwind 4.5 × 10 −7 M /yr for 3XMM J181923.7-170616.All these estimates are consistent with the large range observed for latetype giants.However, via independent methods based on the donor star properties, van den Eijnden et al. (2018b) estimate an upper limit of Ṁwind 7 × 10 −8 M /yr on the mass loss rate in GX 1+4.While estimating wind mass loss rates is challenging, if true, that estimate would rule out a wind origin of the radio emission of GX 1+4.
During the first (and likely during the second) radio observation, 4U 1954+31 was observed to be in a faint X-ray state by Swift/BAT and MAXI monitoring.This source shows X-ray variability on a wide range of time scales, from flares lasting hundreds of seconds pointing to a very clumpy, inhomogeneous wind in density and ionisation, to slower variations on 200-400 day time scales (Masetti et al. 2007;Enoto et al. 2014).The origin of the latter, slow evolution is unknown, but its long time scale makes it unlikely to be related to local inhomogeneities in the wind.Instead, it might relate to more large-scale changes in wind velocity, mass loss rate, or ionisation, due to the orbital phase or changes in the donor.In such a scenario, the low-flux state where we caught 4U 1954+31 might result from lower ionisation or mass loss rate, or high wind velocity.All three of those factors decrease the expected radio luminosity of the wind.
No radio spectral information was available for 4U 1954+31 and 3XMM J181923.7-170616, while it was for the final three observations of GX 1+4.Of those, only the second observation shows a spectral shape (α = −0.6 ± 0.6) inconsistent with the expected thermal wind spectrum at 1σ, although the uncertainty is large and α varies by more than 1σ between the observations.We can also briefly compare the observed radio properties of the three SyXRBs with a jet scenario.The measured spectral indices for GX 1+4 are consistent with a synchrotron-emitting compact radio jet, while the radio luminosities of the two detected SyXRBs fit with the distribution seen from other stronglymagnetised neutron stars.Therefore, a jet could also plausibly explain the radio emission.Without a better theoretical understanding of jets from strongly-magnetised neutron stars (see Section 5.1), there are no further tests of this scenario that we can perform with the current data.However, once the wind has been captured by the neutron star, we do not expect significant differences with jetlaunching HMXBs.
The non-detection of 3XMM J181923.7-170616 is not surprising in either the wind or jet interpretations.Firstly, its 8-kpc distance yields an upper limit of LR 2 × 10 28 erg/s, while we find that no strongly-magnetised neutron stars appear above a similar maximum luminosity.In addition to its large distance, we do not find a strong constraint on the wind mass loss rate corresponding to the upper limit: it is realistic that the stellar wind indeed follows Ṁwind 4.5 × 10 −7 M /yr (Espey & Crowley 2008).Finally, a large orbital size might imply a low wind capture rate and ionisation, possibly decreasing the wind and jet luminosities.However, GX 1+4 likely has a 1161 day period (Hinkle et al. 2006), and for this explanation, 3XMM J181923.7-170616 would require an even larger orbit.
The white dwarf analogues of SyXRBs, symbiotic stars, are much more numerous and have been characterised in detail in radio over the past decades, allowing for an interesting comparison.Seminal early studies by, e.g., Seaquist et al. (1984), Seaquist &Taylor (1990), andSeaquist et al. (1993) detected unresolved radio emission from a large fraction of symbiotic stars at radio luminosities reaching up to 10 30 erg/s (converted to 10 GHz).This emission is attributed typically to thermal wind radiation, requiring mass loss rates of the order of 10 −7 -10 −10 M /yr.More recent studies have confirmed that those symbiotic stars that show hydrogren shell burning on their surface show unresolved thermal wind emission at luminosities between 10 28 -10 30 erg/s, while the nonshell-burning systems tend to be a factor 10-100 fainter in radio (Weston et al. 2016a,b).Possibly, this difference might arise from a difference in mass accretion rate, which could correspond to lower mass loss rate and therefore wind radio luminosity.Alternatively, the X-ray photons from shell burning might enhance the ionisation of the stellar wind (see Sokoloski et al. 2017, for a recent review).
Symbiotic stars also launch jets: MWC 560 launches an unresolved radio jet (Lucy et al. 2019), while in several close-by symbiotic stars, resolved jets have also been detected (Padin et al. 1985;Dougherty et al. 1995;Ogley et al. 2002;Brocksopp et al. 2003Brocksopp et al. , 2004;;Karovska et al. 2010).The radio luminosities of these jets lie in the range of 10 27 -10 30 erg/s, overlapping with, but up to higher luminosities than, the range for strongly-magnetised neutron stars.A comparison of stellar wind properties between symbiotic stars and SyXRB might be obvious due to the shared donor star types, but the inner accretion flow has similarities too: as argued before, these non-magnetic white dwarfs have sizes similar to the magnetospheric radius of a 10 12 G, 1.4 M neutron star accreting at LX ≈ 10 36 erg s −1 .Therefore, a jet comparison between SyXRBs and symbiotic stars can be valid as well.
What does a radio comparison with symbiotic stars reveal about SyXRBs?Firstly, the wind loss rates required to explain the SyXRB radio emission are consistent with those seen in symbiotic stars.More interestingly, despite similar donor wind properties, symbiotic stars reach up to two orders of magnitude higher radio luminosities.What could explain such a difference?Those high radio luminosities are seen in shell-burning systems, which have a continuous, additional source of ionising photons.This can maintain a higher degree of ionisation in the wind and thereby increase its radio emission.Similar thermonuclear burning is possible on neutron stars, but is not sustained for similar lengths of time and does not occur on strongly-magnetised neutron stars (e.g.Galloway et al. 2008).In addition, the accretion flux is significantly softer and therefore more ionising in symbiotic stars comared to SyXRBs.As an alternative to ionisation, we might simply observe an effect of small number statistics: the number of known SyXRB is of the order of ten, while the number of known and candidate symbiotic stars is of the order of thousands (Belczyński et al. 2000;Corradi et al. 2008Corradi et al. , 2010)).Therefore, we might simply not know of any SyXRBs with high enough mass loss rates to reach higher radio luminosity.Finally, the difference in maximum radio luminosity might arise from a maximum jet luminosity for stronglymagnetised neutron stars (see Section 5.1), if a jet dominates the radio emission from SyXRBs.
Concluding our discussion of SyXRBs, we find that their radio properties fit either a jet or stellar wind interpretation.Symbiotic stars offer an interesting comparison, but do not make allow us to disentangle these scenarios for SyXRBs.Following our ap-proach for Vela X-1 and 4U 1700-37, we will in the remainder of this discussion assume that a jet contributes to the observed radio emission, in order to derive as complete as set of constraints on jet physics from strongly-magnetised neutron stars as possible.

IMPLICATIONS FOR NEUTRON STAR JET PHYSICS
After discussing the possible physical origins of the detected radio emission in the previous section, we will now turn to the implications for jet physics.We stress that for this section, we assume that all detected radio emission originates in a jet in order to use the most information possible.However, as detailed in the previous section, it is not certain that all detected sources indeed launch a jet.

The radio brightness of neutron stars' jets
As shown in Figure 3, despite some overlap, the typical radio luminosities of accreting BHs, weakly, and strongly magnetized NSs differ substantially: black holes are, as a sample, radio brighter than weakly-magnetised neutron stars, which are in turn radio brighter than strongly-magnetised neutron stars.This difference was known already between the first two source classes (e.g.Fender & Kuulkers 2001;Migliari & Fender 2006), and is further confirmed by our new observations (e.g.Gallo et al. 2018, which already include some of our preliminary results).We have extended this comparison to stronger magnetic field sources, which also possess systematically slower spin periods (e.g.ν 1 Hz) than their weaklymagnetised counterparts (e.g.ν > 100 Hz for AMXPs).We will investigate the role of magnetic field and spin in the next section, but here we will first discuss the radio luminosity trends in broader terms.
Comparing the black hole and neutron star X-ray binaries, the former appear to be radio brighter over the entire explored range in X-ray luminosity (e.g.∼ 10 33 -10 39 erg/s).However, comparing the two classes of accreting neutron stars in more detail, the story is more complicated.While generally, strongly-magnetised neutron stars are fainter, this effect establishes itself most clearly above X-ray luminosities of ∼ 10 37 erg/s (see Figure 1).Between ∼ 10 35 and 10 37 erg/s (i.e.∼ 0.1 -10% of the Eddington luminosity for neutron stars), the strongly-magnetised neutron stars overlap with the fainter detections of their weakly-magnetised counterparts.In the highest X-ray luminosity ranges, the comparison is dominated by Swift J0243.6+6124, as that parameter space is poorly explored for HMXBs.However, we do not expect that these sources -or, at least, the persistent sources -can be much radio brighter than Swift J0243.6+6124, as otherwise they would have been detected more easily in earlier decades.Therefore, it appears as if the strongly-magnetised neutron stars reach a ceiling radio luminosity of ∼ 3 × 10 28 erg/s, while the weakly-magnetised neutron stars continue a scattered but broadly-coupled increase between X-ray and radio luminosity.This difference becomes apparent above 10 37 erg/s.
One can then ask whether the overlap between the two neutron star samples between LX ∼ 10 35 and 10 37 erg/s is merely a coincidence or hints towards similarities in the launching process.For instance, at lower X-ray luminosity, the characteristic radii for the jet launching process might be located further away from the neutron star -whether these radii are set by the size of the Comptonising medium, the inner radius of the accretion disc, or the location and height of the jet base.At such launching radii, less gravitational potential energy is available to tap for the launch of an outflow.As the accretion rate and X-ray luminosity increase, and the characteristic radii move inwards, more energy becomes available.However, in the strongly-magnetised neutron star case, the magnetospheric radius, where accretion flow and magnetosphere pressures are equal, will be significantly larger and could, in this scenario, create a maximum power available to launch a jet (where we assume that radio luminosity and power are directly related).This description is purposefully vague, as we discuss more detailed models in the next section.However, while we cannot exclude a coincidental overlap in radio luminosity and this scenario ignores the significant differences in accretion flow structure between Roche-lobe and windaccreting systems, this scenario might qualitatively account for the observed divergence at higher accretion rates.

The role of neutron star magnetic field and spin
In this section, we will discuss in more detail whether the neutron star magnetic field and/or spin can affect the jet power; or, phrased more accurately, whether our observations contain evidence for such an effect of the neutron star's properties on the radio jet.The brief answer is that we find no such evidence, beyond the maximum radio luminosity of strongly-magnetised neutron stars that might be attributable to their strong magnetic fields, as hypothesised in the previous section.

Fast-rotating, weakly-magnetized versus slowly-rotating, strongly-magnetized neutron stars
The weakly-and strongly-magnetised neutron star samples overlap in radio luminosity in the X-ray luminosity range between ∼ 10 35 -10 37 erg/s.We cannot exclude that this similarity in spectral shape and radio brightness is completely due to coincidence, while the jets of strongly and weakly magnetised neutron stars would be launched by completely independent mechanisms.For instance, such a scenario might be at work if the radio observing band probes regions down the jet where the emission properties are uncoupled from the exact launching mechanism.However, assuming a relation between launch mechanism, jet power, and jet luminosity (but see below), this scenario still requires a similar jet power by coincidence.
Alternatively, the overlap in radio luminosity could result from a similar jet launching mechanism for both classes of accreting neutron stars.In that case, the enormous differences in magnetic field of three to four orders of magnitude, in neutron star spin of up to six orders of magnitude, and in accretion flow geometry between Roche-lobe-overflowing and wind-accreting systems, do not have a large effect on the radio luminosity of the jet (beyond the apparent maximum for strongly-magnetised systems).This lack of effect might imply that the jet launching mechanism does not depend, in any significant and detectable manner, on the neutron star properties or mass transfer type.On the other hand, these properties might have an effect, but degeneracies between their influence could mean that our sample remains too small to disentangle these.Finally, there is the recurring challenge in radio observations of jets that these might not directly probe the jet power.What we measure in this work is the radio flux at maximally two frequencies, while the conversion to jet luminosity depends for instance on the complete, unmeasured, broadband jet spectrum.Our limited information about the actual jet power therefore further washes out any effects that magnetic field and spin might have -not only in our sample, but also in the archival comparison data.

Among slowly-rotating, strongly-magnetized neutron stars
Prior searches of spin and magnetic field effects, such as by Migliari et al. (2011a), have always focused on the class of weaklymagnetised neutron stars.With the first sample of radio detections from strongly-magnetised neutron stars, we can search for such effects within this source class as well.Focusing on this subgroup circumvents the question whether different jet launch mechanisms exist for the weakly and strongly-magnetised types, as we will assume that within this subgroup, only a single mechanism is responsible for jet launching.
While neutron stars in HMXBs tend to have similar magnetic fields -or at least within an order of magnitude, between 10 12 -10 13 G -they show a spread in their radio luminosities at similar LX (e.g. Figure 3).Most notably, there is the difference between detected and non-detected targets that is not merely due to differences in distance, despite their similar neutron star properties and X-ray luminosity.What could explain these variations between sources?Firstly, it does not clearly correlate with subclass, as a BeXRB, SgXBs, SyXRBs, and an IMXB are all detected.Secondly, it could simply be another example of the spread in radio luminosity that is also observed in LMXB jets, regardless of the accretor nature (Gallo et al. 2018): both black hole and neutron star X-ray binaries show such spread at a similar level.While this spread remains poorly understood (Tudor et al. 2017;Espinasse & Fender 2018;Motta & Fender 2019) both for individual sources and the sample as a whole, it shows that jet launching mechanisms do not create a one-to-one relation between X-ray and radio luminosity.Given the small difference between maximum and sensitivitylimited radio luminosity for strongly-magnetised neutron stars (i.e.3×10 27 to 3×10 28 erg/s), a measurement of the scatter in the radio luminosities would be heavily biased by the sensitivity limit, artificially posing an upper limit on the amount of observable scatter in the detected sources.
Could the neutron star magnetic field or spin contribute to or explain the spread in observed radio luminosity (e.g.Migliari et al. 2011a)?For this, we consider several scenarios (where we, again, simplistically use radio luminosity as a direct proxy for jet power).Firstly, the radio luminosity does not scale with the magnetic field in our sample: all strongly-magnetised accreting neutron stars have similar magnetic fields (when measured), which cannot account for the variation between sources that we observe.The targets in our sample span a much larger range in spin frequencies, between < 10 −3 and > 1 Hz.However, plotting the radio luminosity at 6 GHz versus neutron star spin in Figure 5 (top left panel) does not reveal a clear relation -in fact, the radio brightest object, 1E 1145.1-6141, is amongst the most slowly rotating neutron stars in our sample (3.367 × 10 −3 Hz), while the radio detected but faint Vela X-1 has a similar spin of 3.5 × 10 −3 Hz.Of course, searching for a scaling with only spin and not taking the effect of mass accretion rate differences into account, is a naive approach.However, no scaling between X-ray and radio luminosity is apparent from the data, nor have we yet considered any model of these effects.
Such models, that take into account magnetic field, spin, and mass accretion rate, do exist.For instance, we can consider whether the launching mechanism might be related to a magnetic propeller.In the magnetic propeller regime, the radius where the accretion flow's and magnetosphere's pressures are equal (the magnetospheric radius Rm), lies far outside the radius where the neu-tron star's and accretion flow's rotational velocities are equal (the co-rotation radius Rco; Illarionov & Sunyaev 1975; D'Angelo & Spruit 2012).As the Keplerian velocity of the disc decreases with distance to the neutron star, this implies that the magnetosphere rotates faster than the flow and expels material in an outflow, instead of channeling it to the magnetic poles.Whether a source sits in this regime depends on a careful balance of magnetic field strength and mass accretion rate, which set Rm, with the neutron star spin, which sets Rco.
During our radio observations, none of the detected stronglymagnetised neutron stars were at luminosities which would have been likely to place them in the propeller regime, regardless of the calculation assumptions about the accretion flow structure (i.e.disc or spherical accretion).This supports the hypothesis put forward in van den Eijnden et al. (2019) to explain the sudden turn on of the jet of Swift J0243.6+6124 during an X-ray re-brightening: depending on the exact magnetic field strength, this turn on could be matched to the transition out of the propeller regime, allowing material to reach closer to the neutron star into a possible jet launching region.Therefore, in such a scenario, no jet is expected for sources residing in their propeller regimes.
A more extensive model for the launch of neutron star jets, regardless of magnetic field strength, is the model proposed by Parfrey et al. (2016).Both analytically and through simulations, it investigates the possibility that the accretion flow opens up magnetic field lines, releasing energy to launch outflows.The strength of the outflow in this scenario scales with magnetic moment to the power 6/7 and with the spin frequency squared.Therefore, we plot the radio luminosity as a function of this scaling in Figure 5 (bottom left panel) in arbitrary units: (B/10 12 G) 6/7 × (ν/1Hz) 2 .For sources where the magnetic field is not measured from the cyclotron line, we plot a range corresponding to B = 10 12 -10 13 G, as is typical for these slow pulsars.As in the top left panel in the same figure, we again highlight Swift J0243.6+6124 and GX 1+4, which have multiple observations.However, this exercise again reveals no obvious scaling or trends.
So far, we have neglected the effects of variations in mass accretion rate.While the global X-ray -radio luminosity diagram does not clearly reveal such a correlation in the stronglymagnetised sample, we can consider it for the Parfrey et al. (2016) model: this model also includes a scaling with X-ray luminosity, as a probe of the mass accretion rate, with a power of 4/7.Therefore, we add this scaling in the two right panels of Figure 5.We have removed the two sources that were only detected in radio and therefore do not have a measured X-ray luminosity (i.e.4U 1954+31 and IGR J16318-4848).However, including this correction, we again find no clear trend or correlation.So, to conclude, our radio sample study of stronglymagnetised accreting neutron stars does not show any obvious evidence for a coupling between radio luminosity and neutron star spin, magnetic field, or mass accreting rate.The only possible effect of magnetic field strength, revealed through the comparison with weakly-magnetised neutron stars, is the apparent presence of a ceiling in jet luminosity for strong-magnetic field neutron stars.Beyond that, the scatter in radio luminosities is similar in the strongly and weakly-magnetised neutron stars, as is it in the black hole LMXBs.For three source classes, we do not obtain any new radio detections: UCXBs, VFXBs, and BeXRBs.The former two are both classes of LMXBs, that overlap: UCXBs are binary systems with orbital periods shorter than 1 hour, while the VFXBs are persistent or transient X-ray binaries with maximum X-ray luminosities of ∼ 10 36 erg/s (Wijnands & Degenaar 2016).A possible explanation for their faintness is that the VFXBs are UCXBs with small accretion disk.Alternatively, in neutron star systems, a dynamically active magnetic field could also inhibit efficient accretion and high X-ray luminosities (e.g.Starting with the UCXBs, the deepest constraints on radio emission are found for 2S 0918-549 and 4U 1246-588: both are detected below 10 36 erg/s in X-rays while their radio luminosity does not exceed ∼ 3 × 10 27 erg/s -an order of magnitude lower than typical radio-detected weakly-magnetised neutron stars at similar X-ray luminosities.From the perspective of jet formation, UCXBs might offer an interesting view given their non-standard disc composition: while ordinary LMXBs have a hydrogen-dominated disc, the discs of UCXBs are typically hydrogen-deficient (Nelemans & Jonker 2010;Hernández Santisteban et al. 2019).In the case of a white dwarf donor, the disc might instead be formed predominantly of helium, with high abundances of carbon, oxygen, neon, and/or magnesium (Juett et al. 2001).Such a different composition changes the ionisation properties of the disc (Ludlam et al. 2019) as well as the charge-to-mass ratio, possibly affecting the interaction between the disc and magnetic field and the formation of outflows.
Reported radio observations of UCXBs typically detect a jet, assuming that the source's X-ray luminosity exceeds LX ∼ 10 36 erg/s.Comparing the donor types of these reported UCXBs (e.g.Rappaport et al. 1987;Nelemans et al. 2004;Madej et al. 2013;Homer et al. 1996;Dieball et al. 2005;Galloway et al. 2002;Sanna et al. 2017Sanna et al. , 2018) ) and the four in our sample (Hemphill et al. 2019;in't Zand et al. 2005;van den Eijnden et al. 2018d;Hernández Santisteban et al. 2019) reveals no systematic difference between the two classes: both groups are dominated by white dwarf donors, of various types when identified.Therefore, differences in accreted material or donor type do not appear to explain the difference in radio brightness.
The radio behaviour of VFXBs is interesting for a different reason: several authors have previously suggested that these systems might be transitional millisecond pulsars, or tMSPs (Heinke et al. 2015;Degenaar et al. 2017).These systems are neutron star binaries switching between an X-ray binary and radio pulsar state (Archibald et al. 2009;Papitto et al. 2013), of which only three are known and which are thought to form an evolutionary link between X-ray binaries and radio millisecond pulsars.At higher X-ray luminosities than we targeted for the VFXBs XMMU J174716.1-281048,1RXH J173523.7-354013, and AX J1754.2-2754, the tMSP M28i appears relatively radio bright within the neutron star sample (e.g Deller et al. 2015;Jaodand et al. 2018).Our radio upper limits on the three VFXBs lie firmly below the radio luminosities of the tMSPs at low X-ray luminosities.That argues that either these three VFXBs are not tMSPs5 , or that no single tMSP X-ray -radio correlation exists.
The final undetected class of sources in our new radio observations are the BeXRBs.The specific aim of observing these targets was to search for radio signatures of the propeller regime: with their low accretion rates and strong magnetic fields, BeXRBs would be prime targets to observe this state.Out of the four sources, V 0332+53 clearly resided in the propeller state during the radio observations: we measure an X-ray luminosity upper limit of LX < 3 × 10 34 erg/s, while the maximum X-ray luminosity of the propeller state is estimated to be L lim = 3 × 10 35 erg/s (fol-lowing Tsygankov et al. 2018, and assuming disk accretion, i.e. k ≡ 0.5.See Table 2 in the Online Supplementary Materials for the used source parameters).Therefore, we can place an upper limit of LR < 3 × 10 27 erg/s on the radio emission of any type of propeller outflow that might be launched in this system.
The two X-ray detected BeXRBs, V*V490 Cep and MXB 0656-072, did not reside in their propeller states during out observations: we detect these sources at respectively LX = 8.9 × 10 34 erg/s and 7.6 × 10 33 erg/s, above their respective propeller limits of L lim = 4.6 × 10 33 erg/s and 5.0 × 10 32 erg/s.Finally, SAX J2239.3+6116 was not detected at LX < 8 × 10 33 erg/s, while its propeller regime is expected below L lim = 6 × 10 29 erg/s -significantly below the observed upper limit.However, at the order of magnitude of the propeller regime limit, the source is likely not actively accreting.In addition, both MXB 0656-072 and SAX J2239.3+6116spin at frequencies below 0.01 Hz.In the cold disk model by Tsygankov et al. (2017), only neutron stars spinning faster than 0.01 Hz are capable of entering the propeller state; otherwise, the disk recombines into its un-ionised state before decreasing in X-ray luminosity below the propeller limit.Therefore, we conclude that these three BeXRBs were not in their propeller states.

Soft state neutron stars: jet quenching and spectral shapes
The previous subsection discussed sources that were radio nondetected at low X-ray luminosity.However, radio non-detections can also occur at higher X-ray luminosities, in particular in the soft state, where the inner disc is expected to be geometrically thin.This phenomenon, where the radio emission from the compact jet quenches as the source transitions between hard to soft state, is observed in all black hole LMXBs with good monitoring across the transition or in both states (Tananbaum et al. 1972;Harmon et al. 1995;Fender et al. 1999;Gallo et al. 2003;Fender et al. 2009;Miller-Jones et al. 2012).This state transition is also often accompanied with optically thin radio flaring, likely associated with the launch of discrete ejecta (Fender et al. 2004).Once the black hole has fully entered the soft state, no compact core radio emission is detected anymore (Coriat et al. 2011); most extremely, upper limits on the compact core emission have been measured more than 3.5 order of magnitude below the hard state flux in the black hole LMXB MAXI J1535-571 (Russell et al. 2019).Once black hole LMXBs later return to the hard state, at lower X-ray luminosity, the compact jet re-establishes (Fender et al. 2004).
For neutron star X-ray binaries, this picture is more complex (note: we only focus on LMXBs here, where hard and soft states equivalent to black hole systems can be identified).In three neutron star X-ray binaries, the radio emission has been studied in both the hard and the soft state.Comparing these two states, Migliari et al. (2003) find only marginal evidence of radio quenching in the persistent X-ray binary 4U 1728-34.For the transient Aql X-1, Miller-Jones et al. (2010) report quenching by at least an order of magnitude in the soft state (although a more recent campaign also caught a much more radio bright soft state of Aql X-1; see Díaz Trigo et al. 2018, and below).Finally, Gusinskaia et al. (2017) observed how the transient 1RXS J180408.9-342058quenched as well during its transition to the soft state.Considering these three sources, one notices that the difference between hard state detections and soft state limits is not as large as in black hole LMXBs, as a result of the faintness of the neutron star systems.Four more atoll sources have been observed in radio only in their soft state: three of those (4U 1820-30 and Ser X-1 in Migliari et al. 2004, and MXB 1730-355 in Rutledge et al. 1998) were detected while GX 9+9 was not (Migliari 2011).Whithout a hard state observation, however, it is unclear whether the three detected sources significantly differ in radio luminosity between states.All combined, the picture emerges that some neutron stars quench while others do not, and the origin of this difference remains debated (e.g.Fender 2016;Gusinskaia et al. 2017).
The sample of 13 neutron star LMXBs presented in this paper contains several sources in their soft state: GX 9+1, GX 9+9, and GX 3+1 all persistently reside in the soft state.The persistent atoll 4U 1702-429 switches between hard and soft states.While no simultaneous, high-quality X-ray data was taken during its radio observation, we think it is likely that 4U 1702-429 also resided in a soft state: during a Swift observation taken 8 days before the radio epoch, the source showed an X-ray luminosity of ∼ 1.4 × 10 37 erg/s, or close to 10% of the Eddington limit for neutron stars.For comparison, when the source was observed to be in its hard state with NuSTAR by Ludlam et al. (2019), its X-ray luminosity was a factor 10 lower.The XRT spectrum also requires the addition of a soft thermal component, not seen in this source's hard state, at a significance of 3.7σ (p = 0.0001).This interpretation is supported by Swift/BAT monitoring, shown in Figure 6, wherein 4U 1702-429 switches between brighter (∼ 0.015 ct/s/cm 2 ) and fainter states (< 0.01 ct/s/cm 2 ).As Swift/BAT is sensitive at hard X-ray energies, i.e. 15-50 keV, the brighter state corresponds to the source's hard state.At the time of the radio observation, 4U 1702-429 was in the faint Swift/BAT state, that had already started during the Swift/XRT observation.Finally, as discussed more extensively later on in this section, it shows a steep radio spectral index, which is not typically seen in neutron star hard states (compare with, for instance, IGR J17379-3747 in our sample).
What do these four soft state atolls tell us?Firstly, we do not detect radio emission from GX 9+9, confirming the earlier results from Migliari (2011), and GX 9+1, adding an extra example of a radio-undetected soft-state-only neutron star.However, to complicate matters, GX 3+1, which is very similar to the other two sources in terms of its X-ray properties, is detected in radio.While black holes appear to always quench in the soft state, neutron stars show two types of behaviour, as highlighted by these three sources.Finally, 4U 1702-429 is detected as well, with a steep radio spectrum between the two ATCA bands at 5.5 and 9 GHz.
As none of the four sources above has been observed in radio during a hard state, we stress that a radio detection does not automatically exclude quenching: if GX 3+1 were to be a factor > 10 times radio brighter in the hard state (i.e.> 7 × 10 28 erg/s), this would be consistent with other hard state atolls (see Figure 2).4U 1702-429, on the other hand, would be among the radio brightest neutron stars in its hard state (i.e.LR > 4 × 10 29 erg/s) in the case for a similar soft state radio quenching factor of 10.
Neutron stars, both in the literature and our sample, thus show a range in soft state radio behaviour.To better understand this variance, we can consider the nature of the radio emission during the soft state, if present.Does it originate from the compact jet, or do we observe discrete ejecta?When the spectral index has changed, do we observe a change in jet break frequency?To answer such questions, we need to consider sources with multi-band observations, such as 4U 1702-429: how does its negative spectral index compare to other (soft state) neutron star X-ray binaries?
The best jet spectral shape constraints for neutron stars are the four cases where the jet spectral break is measured: 4U 0614+091, a persistent atoll and UCXB observed in its hard state with spectral break frequency ν b ≈ (1 − 4) × 10 13 Hz (Migliari et al. 2010); 1RXS J180408.9-342058, the aforementioned transient with a hard state ν b ∼ 5 × 10 14 Hz (Baglio et al. 2016, see also Díaz Trigo et al. 2017); the persistent source 4U 1728-34, where a ν b of (1.3 − 11) × 10 13 Hz was measured during the transition from its hard to soft state (Díaz Trigo et al. 2017); and the transient Aql X-1 (Díaz Trigo et al. 2018).In the latter, the break frequency was constrained in four epochs, with a break frequency in the range of 30-100 GHz during the hard state but a break frequency below the ATCA band (ν b < 5.5 GHz) during the soft state.Finally, deep ATCA+ALMA observations of the persistent UCXB 4U 1820-30 in its soft state did not detect a spectral break, implying a inverted radio-to-submm spectrum.
In our sample, we cannot infer a spectral shape for the two non-detected soft-state sources, and GX 3+1 is too faint for a meaningful spectral shape measurement (26.6 ± 5.1 µJy).4U 1702-429, we measure a spectral index of α = −0.8± 0.3 between 5.5 and 9 GHz.The source likely transitioned into its soft state 15 days prior (e.g. Figure 6), suggesting relatively steady soft state emission: while it is possible we caught an optically thin radio flare associated with the launch of individual ejecta, that would be highly coincidental.Instead, we deem it more likely that we are observing a more steady radio source with a spectral break frequency below the ATCA bands.
Prior to this work, Díaz Trigo et al. (2018) suggested that the soft state spectral shapes (up to the ALMA band) of neutron stars might depend on their transient or persistent nature: the transient Aql X-1 shows a steep soft-state radio spectrum, while the persistent 4U 1820-30 shows an inverted spectrum in that state.This also fits with the results for 4U 1728-34, which showed a spectral break in the nIR, although it was only transitioning into the soft state (Díaz Trigo et al. 2017).However, the spectral index measurement in 4U 1702-429 counters this suggestion, showing once more that the differences and similarities in soft state radio behaviour of neutron stars remain challenging to summarise.If the negative spectral index has a similar origin as that in Aql X-1, it suggests that the spectral break might be located at higher frequencies during its hard state.More observations of 4U 1702-429, particularly in its hard state, would be valuable to see how the spectrum and flux evolves between states.Finally, we remind the reader of an important caveat: during one of its hard state epochs with an inverted broad-band spectrum between the ATCA and ALMA bands, Aql X-1 also showed a lower radio flux at 9 than at 5.5 GHz (Díaz Trigo et al. 2018).The spectral index measurement for 4U 1702-429 using only ATCA might suffer the same effect, although our data does not allow us to test this option.

5.3.3
The nature of the compact object in 4U 1700-37 4U 1700-37 is a well-known HMXB system that was discovered by the Uhuru mission (Jones et al. 1973).It is located relatively nearby at ∼ 1.9 kpc (Ankay et al. 2001) and its companion OB supergiant HD 153919 is the hottest and most luminous HMXB donor star known (van der Meer et al. 2004).The nature of the compact object -neutron star or black hole -remains unconfirmed.While its X-ray spectrum resembles that of accreting X-ray pulsars, no pulsations or a cyclotron line have been unambiguously detected (although see Reynolds et al. 1997, for a possible cyclotron line detection at ∼ 37 keV).The mass of the compact object is possibly larger than 2M , consistent with a heavy neutron star, similar to Vela X-1, or a low-mass black hole (Rubin et al. 1996;Clark et al. 2002;Falanga et al. 2015).
In terms of radio flux density, 4U 1700-37 is the brightest HMXB in our sample.However, as it is also the most proximate, its position in the radio -X-ray luminosity plane is consistent with the neutron star sources (see Figure 2).This position is completely inconsistent with the radio-bright track for black hole X-ray binaries, and roughly an order of magnitude below the radio-faint track of these sources.While not many black hole HMXBs are known, all radio observations of these systems currently place them at the radio-brightest part of the correlation seen in LMXBs (Ribó et al. 2017).Therefore, while we stress that it does not amount to direct proof, the radio luminosity of 4U 1700-37 is strongly suggestive of a neutron star primary.

Deep radio constraints in the literature: 4U 2206+54, A 0535+262, and X Per
Finally, we briefly turn to the deepest radio non-detections of strongly-magnetised neutron stars in the literature.With our survey, we have probed down radio luminosities of approximately 10 27 erg/s.Combined with the apparent maximum radio luminosity of these sources of LR ≈ 3 × 10 28 erg/s, this limited range implies that we merely scratched the surface in terms of the full population.
What can be expected in future studies probing deeper in radio luminosity?Before the VLA and ATCA sensitivity upgrades in the past decade, after which all our observations were taken, deep radio observations were made of the neutron star HMXB 4U 2206+54 (Blay et al. 2005).A non-detection of FR < 39 µJy at 8.4 GHz, given a source distance of ∼ 3 kpc, translates into an upper limit of LR < 3.5 × 10 27 erg/s -still comparable with our sensitivity.While no simultaneous X-ray detection was obtained, two INTE-GRAL observations within a few weeks places the source at a 4-12 keV luminosity of 4.5-8.5 × 10 35 erg/s.Westerbork Synthesis Radio Telescope observations of A 0535+262, at a distance of only 2 kpc, similarly resulted in a non-detection of LR < 3.8 × 10 27 erg/s (Migliari et al. 2011b).While no X-ray luminosity is reported, the MAXI monitoring count rate suggests a 2-20 keV luminosity around 6 × 10 36 erg/s, assuming a Crab spectrum for simplicity.
The deepest existing radio observation of a HMXB neutron star is, to our knowledge, that of the persistent source X Per, with a spin period of 835 seconds (White et al. 1976).The Westerbork Synthesis Radio Telescope did not detecting the source down to LR 8 × 10 25 erg/s (Nelson & Spencer 1988).This exceptionally deep upper limit, due to the small distance (∼ 800 pc) to X Per, shows that our radio non-detections are not necessarily due to limited sensitivity alone.Unfortunately, no X-ray information is reported for the time of the radio non-detection.Therefore, it remains unclear what accretion state and rate X Per showed around the radio observation.However, assuming X Per was active, its radio upper limit shows that other accreting, strongly-magnetised neutron stars might also not be detected with future instruments, such as the next-generation Very Large Array (see Section 6 Selina et al. 2018;Coppejans et al. 2018), because their jets are either incredibly weak or simply do not exist.

CONCLUSIONS AND FUTURE OUTLOOK
We have presented radio observations of a sample of 36 accreting neutron stars.Based on these observations, we draw the following eight conclusions: (i) Strongly-magnetised accreting neutron stars can be detected at radio frequencies; the previously detected source Swift J0243.6+6124 is no outlier (Section 3).
(ii) Strongly-magnetised accreting neutron stars are, as a sample, radio fainter than their weakly-magnetised counterparts, which are in turn fainter than black holes.This difference is established especially at high X-ray luminosities, e.g.LX 10 37 erg/s, as the radio luminosity of strongly-magnetised systems appears to reach a ceiling value around 3 × 10 28 erg/s (Sections 3 and 5.1).
(iii) No new very-faint X-ray binaries, ultracompact X-ray binaries, or faint Be/X-ray binaries were detected in radio (Sections 3 and 5.3.1).
(iv) Strongly-magnetised accreting neutron stars, as a sample, do not show a clear correlation between X-ray and radio luminosity, although this might be caused by sensitivity limits (Section 3).
(v) Accretion-powered jets explain the radio detections of weakly-magnetised neutron stars, and are the likely explanation for the radio emission from strongly-magnetised neutron stars.The exceptions to the latter statement are Vela X-1, 4U 1700-37, and the symbiotic X-ray binaries, where the stellar wind possibly contributes (Section 4.1).
(vi) Strongly-magnetised accreting neutron stars do not show any obvious scaling between their radio luminosity and their spin, magnetic field strength, and accretion rate, or combination thereof, as predicted by existing models (Section 5.2).
(viii) Like their transient counterparts, persistent neutron star LMXBs can show steep soft-state radio spectra (Section 5.3.2).
Finally, we will end this work by briefly turning to the future.As this work is an observational study, we will focus mainly on the observational developments.Firstly, since the first systematic comparison of accreting black holes and neutron stars in the X-ray -radio luminosity plane (Migliari & Fender 2006), tens of neutron stars have been observed or monitored in radio, amounting to a sample of over sixty X-ray binaries (e.g.Gallo et al. 2018, and this work).It has however proven difficult to obtain monitoring over a large range in X-ray luminosity, due to the inherent faintness of neutron star jets and difficulty in coordinating X-ray and radio observations.Such monitoring would be especially valuable for strongly-magnetised neutron stars, where the only monitored outburst (of Swift J0243.6+6124) has revealed surprising behaviour that will be interesting to compare to other binaries.
The 23 strongly-magnetised neutron stars with observations at current-day sensitivities have merely scratched the surface; these constitute only a small fraction of all such sources.Observing a larger sample, especially down to lower radio luminosities, will be essential to probe whether or not a scattered, but global correlation between X-ray and radio luminosity exists, that is currently cut off by the detection limit.Future arrays, such as the next-generation VLA (ngVLA) with its planned sensitivity of ∼ 0.23 µJy at 8 GHz in one hour of observing time (Selina et al. 2018), will probe 50 times deeper than the current VLA, bringing X-ray binaries in M31 above sensitivity limits.In addition, monitoring of transient sources will complement the persistent sources that currently dominate the sample, while probing the jet's turn-on and fade.
What are the prospects for high spatial resolution observations of neutron star jets?Among the newly-added neutron star X-ray binaries -and the strongly-magnetised neutron stars in generalmany persistent sources might have launched jets long enough to create feedback structures, despite their low radio luminosity.At even higher resolution, very-long baseline observations mightresolve the compact jets of neutron stars depending on the size of the radio emitting region: at the high end of the typical range, ∼ 10 7 -10 9 Rg, this size corresponds to a resolvable ∼ 7 mas at 2 kpc.However, if the radio emission regions are located towards 10 7 Rg, these will not be resolvable.Discrete ejecta in transient sources, however, might be resolved using the newly-employed MeerKAT radio telescope, following such detections in black hole systems (e.g.Russell et al. 2019;Bright et al. 2020).
For a sufficiently radio bright neutron star, VLBI observations could help distinguish between radio emission from a stellar wind or jet through brightness temperature estimates.Stellar winds typically have brightness temperatures of T b ≈ 10 4 -10 5 K (Rybicki & Lightman 1979;Longair 1992;Russell et al. 2016).At 6 cm, the Very Long Baseline Array can reach a spatial resolution of 1.4 mas.If an accreting neutron star, with a flux density 0.5 mJy at 5 GHz, is unresolved at that resolution, it implies a minimum brightness temperature of T b 10 7 K, arguing against a stellar wind origin.Beyond the flux density required, this constraint does not depend on distance.Assuming the apparent maximum radio luminosity of ∼ 3 × 10 28 erg/s holds for all neutron star HMXBs, such a measurement is therefore feasible for sources within ∼ 3 kpc.
Two promising avenues that remain poorly explored for neutron star jets are eclipse mapping and polarimetry.Recently, Maccarone et al. (2020) presented model calculations for the eclipse of the jet by the donor star, showcasing how such measurements can constrain the geomtrical and radiative properties of the jet.Due to the large extend and density profile of the stellar wind in high-mass systems, this technique is most suitable for neutron star LMXBs.Polarimetry could help constrain the origin of the radio emission in HMXBs, as radio jets can be linearly polarized while stellar winds are not expected to be.However, the level of polarization, especially in compact jets, can be lower (i.e.few per cent) than current sensitivities allow to detect (van den Eijnden et al. 2018a).
From a multi-wavelength perspective, the recent detection of a jet spectral break in an accreting neutron star with ALMA (Díaz Trigo et al. 2018) has shown such studies are possible for the brightest neutron star X-ray binaries.Strongly-magnetised neutron stars would be an interesting target for ALMA or ngVLA, going up to ∼ 100 GHz, as their jets are likely launched further out, possibly resulting in lower jet break frequencies.The low luminosity of these sources would pose a challenge for ALMA, but are realistically detectable with the ngVLA.At higher frequencies, in for instance the infrared band, studies of strongly-magnetised neutron stars will often be limited by the bright emission of their high-mass companion stars.However, in sources with low or intermediatemass donors, such as Her X-1, nIR observations could also probe the jet launch regions.
As a final remark, all above suggestions for future studies are observational.These studies will add to the broad set of observational constraints on jet formation already existing.From the side of new jet launch models, testable predictions, such as scaling relations between radio luminosity and magnetic field, spin, and mass accretion rate, are essential to match to the observational constraints and guide the design of new observing campaigns.Similarly, predictions from either analytical theory or GRMHD simulations for the effect of binary orbit properties, mass transfer type, and donor star type on the jets might in the future be testable given the large variety in types of HMXBs.
Table A1.Background information about the radio observations presented in this paper.An * in the ATCA primary calibrator indicates that 0823-500 was also observed but not used the calibration.The observing frequencies are listed in Tables 1 and 2 in the main paper.When two robust parameters are listed, these correspond to the two (increasing) observing frequencies.The synthesized beam parameters correspond to the lowest observing frequency when observations were performed at two frequencies simultaneously.Shortened source names: XMMU J174716.1 ≡ XMMU J174716.1-281048,1RXH J173523.7 ≡ 1RXH J173523.7-354013, and 3XMM J181923.7 ≡ 3XMM J181923.7-170616.When multiple dates are listed for one epoch, the data were combined in the analysis.**Not included due to the presence of bright, proximate background source, yielding any flux density upper limit physically unconstraining.A2.Magnetic field strengths, spin frequencies, and orbital periods of the X-ray binaries in this work.Sources where the magnetic field has not been measured/inferred, are classified as weakly or strongly magnetized based on their (combination of) bursting properties, pulse evolution, X-ray behavior, and/or evolutionary state.We note that all parameters, especially spin frequencies, evolve over time, and therefore refer to the listed references for more details.

Source
Question marks indicate that a measurement might be debated.Again, more details are found in the referenced work(s).does not show evidence for stellar wind emission, supporting the notion that the wind mass loss rate does not greatly exceed typical values.

Source
Finally, we turn to the BeXRBs, of which only Swift J0243.6+6124 was detected in radio.For individual BeXBRs, the stellar wind parameters are typically poorly constrained.However, the Be star winds are significantly weaker than those of the supergiant O/B donors in SgXB: typical mass loss rates are of the order 10 −9 M yr −1 (Snow 1981).At such mass loss rates, even a low terminal velocity (yielding the highest luminosity) of 100 km s −1 leads to expected 5.5 GHz flux densities below 1 µJy for distances larger than 1 kpc.As all observed BeXRBs are located at larger distances (at least 3.6 kpc), their stellar wind are not expected to be observable at current radio sensitivities.

APPENDIX D: SPECTRAL INDEX UPPER LIMITS
In Figures B1 and B2, we show the estimate of the 3σ upper limit on the radio spectral index for the two sources with a detection only at 5.5 GHz but not at 9 GHz: IGR J16320-4751 and 1E 1145.1-6141.For the details of the calculation, see van den Eijnden et al. ( 2019).
Table A3.Radio positions as measured with the CASA-task IMFIT.The uncertainty is calculated as the beamsize divided by the S/N ratio, limited by 10% of the beam size for the maximum accuracy as suggested by the VLA.The position with the best accuracy, out of all frequency bands and observations of the source, is shown.Table A6.Input stellar wind parameters and references for the 5.5-GHz wind radio flux density estimate.For values denoted with a *, no literature measurement was known to the authors -therefore, we assumed typical values of 10 −6 M yr −1 and 400 km s −1 for the wind mass loss rate and terminal velocity, respectively.The range in predicted 5.5 GHz rado flux density values for EXO 1722-363 corresponds to the range in terminal wind velocities.For IGR J16320-4751, no estimates were found in the literature and therefore no estimate is made.SyXRBs are discussed in the main text, BeXRBs in the appendix.(1974) 6.9 ∼ 42 See also Krzeminski (1973)   would not be detected at 9 GHz with a 3σ significance, as a function of spectral index, assuming the observed 5.5-GHz flux density and 9-GHz sensitivity.

Figure 1 .
Figure1.The X-ray -radio luminosity plane for X-ray binaries.Archival observations of black holes are shown as the light grey crosses, while those of weakly-magnetised neutron stars are shown as dark grey circles.The observations presented in this work are shown in blue (weakly-magnetised) and red (strongly-magnetised).Upper limits are shown in open markers, while detections are plotted as filled markers.Apart from GX 1+4 and Her X-1, we do not show previous observations of strongly-magnetised neutron stars, as these are all upper limits at typically unconstraining levels (although we discuss several exceptions in Section 5.3.4).Radio luminosities have been converted to 6 GHz, using a measured spectral index if available, or assuming a flat spectrum otherwise.

Figure 2 .
Figure 2. The same as Figure 1, now labelling individual sources in our sample.Filled markers indicate strongly-magnetised neutron star, while the open ones show weakly-magnetised sources.

Figure 4 .
Figure 4.A visual comparison of the measured radio luminosities of the SgXBs in our sample and their predicted Wright & Barlow (1975) stellar wind radio luminosity.The dashed line indicates the theoretical wind prediction, as a function of the parameterization of wind velocity and mass loss rate in Equation 1.Most undetected target lie below the predicted relation, while three our of four detected sources are too radio-bright, suggesting a second/other dominant emission mechanism.

Figure 5 .
Figure5.The relation between radio luminosity and neutron star parameters for strongly-magnetised neutron stars.In all panels, Swift J0243.6+6124(referred to as Sw J0243) and GX 1+4 are highlighted separately, as these are the only two sources with multiple data points.Top left: radio luminosity plotted versus neutron star spin ν.Top right: radio luminosity plotted versus neutron star spin, corrected for mass accretion rate, parameterised by X-ray luminosity L X following the prescription inParfrey et al. (2016).Bottom left: radio luminosity versus the jet power scaling with magnetic field B and spin of the model byParfrey et al. (2016).Bottom right: same as bottom left, again correcting for mass accretion rate followingParfrey et al. (2016).

Figure 6 .
Figure 6.Swift/BAT 15-50 keV monitoring light curve of 4U 1702-429.The red dashed line show the time of the radio observation.Two distinct states, with count rates around 0.015 and below 0.01 ct/s/cm 2 can be distinguished: the hard and soft state, respectively.During the radio observation, 4U 1702-429 was in the soft state.Data between MJD 58061 and 58130 have been re-binned by a factor 5 to increase signal-to-noise.

Figure B1 .Figure B2 .
Figure B1.The fraction of simulated data sets where IGR J16320-4751 would not be detected at 9 GHz with a 3σ significance, as a function of spectral index, assuming the observed 5.5-GHz flux density and 9-GHz sensitivity.
(Dougherty & Williams 2000;van Loo 2007)007, and references  therein).Through radio and IR observations, fundamental wind properties such as mass loss rate and terminal velocity can be inferred, probing stellar feedback and the effect on stellar evolution.Wind radio emission is typically attributed to one of two emission processes: either thermal Bremsstrahlung emission from the ionised gas in the wind, or non-thermal emission from shocks in the outer wind.As shown byWright & Barlow (1975), the thermal emission is expected to have a positive (α ≈ 0.6) radio spectral index, which is indeed often observed in isolated O stars (e.g.van Loo 2007).Non-thermal shocks locally create steep spectra, i.e. with negative spectral indices α.But, taking into account that the shocks both weaken and peak at lower frequencies as they move away from the star, their cumulative spectrum for single stars also has a positive spectral index(Dougherty & Williams 2000;van Loo 2007).