Quasi-simultaneous radio and X-ray observations of Aql X-1 : probing low luminosities

ABSTRACT

Aql X-1 is one of the best-studied neutron star low-mass X-ray binaries. It was previously targeted using quasi-simultaneous radio and X-ray observations during at least seven different accretion outbursts. Such observations allow us to probe the interplay between accretion inflow (X-ray) and jet outflow (radio). Thus far, these combined observations have only covered one order of magnitude in radio and X-ray luminosity range; this means that any potential radio–X-ray luminosity correlation, L_R ∝ L_X^β, is not well constrained (β ≈ 0.4–0.9, based on various studies) or understood. Here we present quasi-simultaneous Very Large Array and Swift-XRT observations of Aql X-1’s 2016 outburst, with which we probe one order of magnitude fainter in radio and X-ray luminosity compared to previous studies (6 × 10³⁴ erg s⁻¹ < L_X <3 × 10³⁵ erg s⁻¹, i.e. the intermediate to low-luminosity regime between outburst peak and quiescence). The resulting radio non-detections indicate that Aql X-1’s radio emission decays more rapidly at low X-ray luminosities than previously assumed – at least during the 2016 outburst. Assuming similar behaviour between outbursts, and combining all available data in the hard X-ray state, this can be modelled as a steep β =|$1.17^{+0.30}_{-0.21}$| power-law index or as a sharp radio cut-off at L_X ≲ 5 × 10³⁵ erg s⁻¹ (given our deep radio upper limits at X-ray luminosities below this value). We discuss these results in the context of other similar studies.

1 INTRODUCTION

Astrophysical jets (collimated outflows) are an important yet not-well-understood component in the physics of accretion. Despite the development of numerous theoretical models (e.g. Blandford & Znajek 1977; Blandford & Payne 1982; Parfrey, Spitkovsky & Beloborodov 2016), their formation mechanism is still unclear. Observations of jets in various types of accreting systems can help to elucidate this problem. A particularly useful class of objects in the study of jet formation are Galactic transient low-mass X-ray binaries (LMXBs; Lewin, van Paradijs & van den Heuvel 1997) – systems consisting of a compact object (which can be either a black hole, BH, or a neutron star, NS), accreting from a low-mass companion star (⁠|$\lesssim 1\, \mathrm{M_{\odot }}$| for NS systems, and less than a few solar masses for BH systems). During their characteristic episodes of active accretion towards the central object (outbursts), the accretion rate of LMXBs evolves on time-scales of days to months, allowing us to observe the formation and disappearance of jets on practical time-scales.

When a transient LMXB goes into outburst, its multiwavelength luminosity increases by a few orders of magnitude. As the system’s accretion rate increases, accumulated in-falling matter is heated up and emits in X-rays (e.g. Shakura & Sunyaev 1973). Not all of the in-falling material ends up on the central object, however, some fraction is expelled from the system. In some cases, a jet is formed: outflowing accelerated particles produce synchrotron radiation that can be visible at radio wavelengths. Thus, studying the relationship between the radio luminosity (L_R) and the X-ray luminosity (L_X) of such a system allows one to trace the connection between the inflow and outflow of material in LMXBs – as well as various other types of accreting systems (e.g. active galactic nuclei; Plotkin et al. 2012).

Observations of X-ray binaries in the hard X-ray state have shown that L_R and L_X often follow a power-law relation L_R ∝ L_X^β (Corbel et al. 2002; Migliari & Fender 2006; Gallo, Degenaar & van den Eijnden 2018). This demonstrates that both L_R and L_X are functions of mass accretion rate, and (under the assumption that the jet power is a fixed fraction of the accretion power; e.g. Markoff et al. 2003) the power-law index β mostly reflects the properties of the accretion inflow (e.g. how efficient it is at producing radiation; Narayan & Yi 1995). In such scale-invariant coupled disc–jet models, β ∼ 1.4 is usually associated with radiatively efficient accretion and β ∼ 0.7 with radiatively inefficient accretion. In the case of BH accretors, such radiative inefficiency can be explained by the loss of part of the accreting material within the event horizon (Narayan & Yi 1995). In the case of NS accretors, a radiatively inefficient accretion flow could be explained, e.g. by a strong outflow (due to, e.g. propeller effect D’Angelo et al. 2015).

It has been argued previously that NS-LMXBs and BH-LMXBs show different relations between their L_R and L_X, when modelling observations in the hard X-ray state as a power-law relation (Migliari & Fender 2006). More recently, however, Gallo et al. (2018) find that NS- and BH-LMXBs show roughly consistent power-law slopes at the 2.5σ level. They find |$\beta = 0.44^{+0.05}_{-0.04}$| and β = 0.59 ± 0.02 for NS and BH systems, respectively (Gallo et al. 2018).

In comparison to BH-LMXBs, NS-LMXBs are typically radio fainter, by roughly an order-of-magnitude for the same X-ray luminosity (Fender & Kuulkers 2001; Migliari & Fender 2006; Gallo et al. 2018), and hence more challenging to study. A particularly interesting NS-LMXB for studying jet activity is Aql X-1, one of the first Galactic X-ray sources to be discovered (Friedman, Byram & Chubb 1967; Kunte et al. 1973). It is a binary system (orbital period ∼19 h) composed of an NS (Czerny, Czerny & Grindlay 1985) and a low-mass (⁠|$\sim 0.4\, \mathrm{M_{\odot }}$|⁠) main-sequence-like companion (Chevalier et al. 1999). Aql X-1 goes into outburst almost every year (e.g. Campana, Coti Zelati & D’Avanzo 2013; Güngör, Güver & Ekşi 2014), and hence offers regular opportunities to study its L_R–L_X relation.

Based on the X-ray spectral and timing properties, different states are identified during LMXB outbursts (e.g. van der Klis 1994). During its outbursts, Aql X-1 typically evolves from a low-hard X-ray state, peaking in the high-soft X-ray state, and returning back to the hard state at the end of the outburst (Muñoz-Darias et al. 2014). The majority of outbursts of Aql X-1 can be characterized by fast rise and exponential decay of the X-ray luminosity (‘FRED-type’ outburst; Campana et al. 2013). Some of Aql X-1’s outbursts have symmetrical rise and decay (‘Gaussian-like’ outburst; Campana et al. 2013) or have multiple X-ray luminosity peaks. An alternative classification of outbursts of Aql X-1 is based on their duration and peak luminosity (Güngör et al. 2014): low-luminosity outbursts reach a maximum X-ray luminosity of L_X ∼ 10³⁷ erg s⁻¹ before decaying, while high-luminosity outbursts can reach up to a factor of 5 higher in L_X and tend to be longer in duration.

Previous radio monitoring of Aql X-1 was performed for at least seven different outbursts (all of the low-luminosity class). These studies demonstrated that its radio jet behaviour changes depending on the X-ray spectral state and luminosity (Tudose et al. 2009; Miller-Jones et al. 2010a). Flat-spectrum¹ radio emission (indicative of an optically thick jet) was observed during the rise in the low-hard X-ray state or transition between low-hard and high-soft X-ray states (Miller-Jones et al. 2010a; Díaz Trigo et al. 2018), and ‘quenches’ (i.e. disappears) during the peak of the high-soft X-ray state (L_X ∼10³⁷ erg s⁻¹; Miller-Jones et al. 2010a).

Aql X-1 is one of the only three NS-LMXBs in which the correlation between radio and X-ray luminosities has been investigated in detail using a relatively large number of observations. The power-law index was fitted in different studies, using different data selections (see Section 3.2.2), which found β = 0.40 ± 0.07 (Tudose et al. 2009) and |$\beta =0.76^{+0.14}_{-0.15}$| (Tetarenko et al. 2016). This is similar to what has been found in NS- and BH-LMXBs when fitting the populations globally and accepting more scatter in the correlation (Gallo et al. 2018). For two other NS-LMXBs, 4U 1728−34 (Migliari et al. 2003) and EXO 1745−248 (Tetarenko et al. 2016), steeper power-law indices were found: β = 1.5 ± 0.2 and |$\beta =1.68^{+0.10}_{-0.09}$|⁠, respectively. This could suggest a radiatively efficient accretion regime (Migliari & Fender 2006). Different power-law indices have also been observed in individual BH-LMXBs (e.g. the ‘outliers track’ described by Coriat et al. 2011); however, they have been observed to return to the same canonical track for L_X <10³⁵ erg s⁻¹ (Maccarone 2012; Corbel et al. 2013).

The situation can be clarified by observing such systems over the widest-possible range of radio and X-ray luminosities. While the radio–X-ray correlation was probed for a few BH-LMXBs over 8 dex (orders of magnitude) in L_X and 5 dex in L_R, in 4U 1728−34 and Aql X-1 this correlation has previously only been studied over ∼1 dex in radio and X-ray luminosities.

Since NS-LMXBs tend to be more radio faint compared to BH-LMXBs (Fender & Kuulkers 2001; Gallo et al. 2018), radio observations below the peak of their outburst (L_X <10³⁶ erg s⁻¹) have previously been hampered by the limited sensitivity of available radio telescopes. With the upgraded Karl J. Jansky Very Large Array (VLA; Perley et al. 2011), an order-of-magnitude deeper radio observations are now possible. However, such observations are still challenging to schedule because NS-LMXBs often pass through the 10³⁴⁻³⁶ erg s⁻¹ X-ray luminosity range only during the rise and decay of the outburst, which often lasts no longer than ∼3–7 d (Campana et al. 2013; Güngör et al. 2017).

Due to these time constraints and their relative radio faintness, only a few NS-LMXBs have been observed in the luminosity range L_X < 10³⁶ erg s⁻¹. In the case of non-pulsating NS-LMXBs, these observations led exclusively to non-detections (Tetarenko et al. 2016; Gusinskaia et al. 2017; Tudor et al. 2017). However, there are two classes that potentially stand out in this low X-ray luminosity regime (Gallo et al. 2018): accreting millisecond X-ray pulsars (AMXPs; systems that exhibit coherent millisecond X-ray pulsations as a result of accreting material channelling on to the magnetic polar caps of the NS; Wijnands & van der Klis 1998; Patruno & Watts 2012) and transitional millisecond pulsars (tMSPs; NSs that switch between rotation-powered radio millisecond pulsars and an accretion-disc state; Papitto et al. 2013; Jaodand, Hessels & Archibald 2018). Both these classes show comparatively bright radio emission at low L_X (Deller et al. 2015; Tudor et al. 2017). There were previously tentative suggestions that tMSPs might follow their own correlation in the radio – X-ray luminosity plane (Deller et al. 2015). However, recent strictly simultaneous radio and X-ray observations of tMSP PSR J1023+0038 in the low-luminosity accretion-disc state showed anticorrelated radio–X-ray brightness variations that argue against an origin from a collimated jet (Bogdanov et al. 2018). Thus, the accretion-disc state of low-luminosity tMSPs is likely quite different compared to NS-LMXBs in outburst. Consequently, care is needed when comparing different NS-LMXB systems in the L_R–L_X plane because they could be in fundamentally different accretion regimes.

Ultimately, many more systems need to be studied in detail before we can come to any robust conclusions on whether different classes of NS-LMXBs have different L_R–L_X correlations, because of, e.g. the influence that different magnetic field strengths, geometry, or spin rate may have on disc–jet coupling.

Aql X-1 is identified as an intermittent AMXP because of the discovery of 550.27 Hz coherent X-ray pulsations (Casella et al. 2008) that lasted for only 150 s out of the 10⁶ s that the source was observed. Two other AMXPs are also classified as intermittent and have shown pulsations that appear and disappear on time-scales of days to years: HETE J1900.1−2455 (Galloway et al. 2007; Patruno 2012) and SAX J1748.9−2021 (Altamirano et al. 2008). Since Aql X-1 is the most extreme case of intermittency, it is also possible that the origin of the pulsations is different in nature.

In this paper we present sensitive radio observations of Aql X-1 in the low X-ray luminosity range (<10³⁶ erg s⁻¹), taken while the source was decaying from its 2016 outburst. We compare these with archival measurements from this and previous outbursts to study the L_R–L_X correlation of Aql X-1 – under the assumption that the radio – X-ray correlation behaves similarly between outbursts. In Section 2 we present the observations and data analysis. In Section 3 we present the results, and discuss their interpretation in Section 4.

2 OBSERVATIONS AND DATA REDUCTION

2.1 2016 outburst

Aql X-1 entered into outburst at the end of 2016 July (Sanna et al. 2016). The outburst lasted for approximately 2 months, reaching a peak X-ray luminosity of ∼8 × 10³⁷ erg s⁻¹ (for an assumed distance of 4.5 kpc; Güngör et al. 2017), making it the brightest outburst observed from Aql X-1 to date.

Using the VLA, we performed three radio observations of Aql X-1 while it was fading in brightness at the end of its 2016 outburst. These observations were spaced daily (on the 26th, 27th, and 28th of September 2016), and were coordinated quasi-simultaneously with Neil Gehrels Swift-XRT daily monitoring. Table 1 provides a log of all 2016 outburst observations used in the present study.

2.1.1 Radio data and analysis

Aql X-1 was observed with the VLA at three epochs (project ID: 15B-239), each at X-band (8–12 GHz). The first two epochs had total duration ∼1 h and the last was ∼3 h including calibration scans, with ∼0.5 and ∼2 h on source, respectively. During all three epochs the array was in A configuration (synthesized beam ∼0.3 arcsec). Data were initially calibrated using the VLA Common Astronomy Software Application (casa;² McMullin et al. 2007) Calibration Pipeline (version 5.4.1). We used J1331+305 and J1407+2807 as flux and polarization calibrators, respectively, and J1907+0127 as phase and amplitude calibrator. The data were additionally flagged and calibrated following standard procedures within casa.

In our analysis below, we also incorporate results from Australia Telescope Compact Array (ATCA) observations of Aql X-1’s 2016 outburst that were previously presented in Díaz Trigo et al. (2018).

2.1.2 X-ray data and analysis

Based on the three observational epochs obtained using the VLA we selected four quasi-simultaneous Swift-XRT observations (target ID: 34719), all done in photon counting (PC) mode. As such, each of our radio observations is framed by a pair of X-ray observations. Additionally, we selected seven Swift-XRT observations (target ID: 33665) – two done in PC mode and five done in window timing (WT) mode – that are spaced quasi-simultaneously between the aforementioned ATCA radio observations (Díaz Trigo et al. 2018). For more details see Table 1 and Fig. 1. Due to the high count rate of Aql X-1, all observations were corrupted by pile-up, and consequently two observations that were performed in PC mode (observation IDs: 33665073 and 33665088) during the peak of the outburst were excluded.

The basic data reduction and calibration was done using the xrt_pipeline task within heasoft-6.24³. For all nine epochs, the spectrum of the source was extracted from a 30-pixel-wide area (circular in the case of PC mode observations and square in the case of WT mode observations), centred on the source position (taken from Miller-Jones et al. 2010a), while the background was extracted from a surrounding annulus with 30-pixel inner radius and 60-pixel outer radius. The extraction of spectra was done using xselect. To correct for known CCD artefacts we used ancillary arf-files (produced using the xrtmkarf script and exposure map) together with the ‘swxpc0to12s6_20130101v014.rmf’ response file. Each spectrum was grouped into 20-photon bins. In order to account for pile-up, we used the technique described in Romano et al. (2006) and excluded the necessary amount of pixels from the central part of the images (or events in the case of WT mode) before performing spectral model fitting. We used xspec (version 12.10.0c) to perform the spectral model fit and to extract the source flux. For spectral fitting we used photons from the 0.4–10 and 0.7–10 keV ranges in the case of PC and WT mode observations, respectively. We achieved sufficiently good fits (⁠|$\chi ^2_{\nu } \lt 1.1$|⁠) using a simple absorbed power-law model (TBabs * powerlaw) or with an additional blackbody component (TBabs * (powerlaw + bb)) in xspec. By fitting all spectra together with the N_H parameter tied between data sets, we found N_H = 0.39 ± 0.02 × 10²² cm⁻², which is consistent with the value reported by Rutledge et al. (2001). Thus we simply used N_H = 0.4 × 10²² cm⁻² as a fixed parameter while performing fits to individual observations. We added a cflux component to the spectral model in order to determine unabsorbed 1–10 keV fluxes and their associated errors for each epoch. The resulting parameter values, spectral models, |$\chi ^2_{\nu }$| values and number of degrees of freedom (dof) are listed in Table 1.

In order to place our radio observations in context, we produced a hardness-intensity diagram (HID) using daily-averaged count rates of all-sky monitors MAXI⁴ and Swift-BAT⁵. Both count rates were converted to mCrab using the average (over the time period of Aql X-1’s 2016 outburst) count rate of the Crab nebula with the corresponding instrument. Here we define hardness ratio as 15–50 keV Swift-BAT flux / MAXI 2–10 keV flux and intensity as 2–10 keV MAXI flux. The obtained HID is shown in Fig. 1 (right). The colour scale represents the hardness ratio, where yellow is more soft and blue is more hard. For illustration purposes we inferred the hardness ratio for each Swift-XRT observation by interpolating the results from MAXI/Swift-BAT observations taken around the time of the Swift-XRT observations.

In order to directly compare the 2016 outburst with the other seven outbursts (which were observed using RXTE-PCA) we additionally produced HIDs where the intensity and the hardness ratio are defined as 2–16 keV flux and ratio of 9.7–16 and 6.0–9.7 keV fluxes, respectively. Analogous to the above, we converted MAXI count rates into mCrab units using the average (during the 2016 outburst) count rate of the Crab nebula in the corresponding energy bands. The mCrab units were then converted to |$\mathrm{erg\, s^{-1}\, cm^{-2}}$| flux (see Fig. 2).

As shown in Fig. 1, Aql X-1 was in the hard X-ray state during all VLA observations. The first VLA observation was performed just 6 d after the last ATCA detection, which was performed a few days after Aql X-1 had started the spectral transition to the hard state. As such, two ATCA observations of Aql X-1’s 2016 outburst were performed when the source was still relatively hard and two when it was in the soft X-ray state.

2.2 Archival data of Aql X-1 from previous outbursts

2.2.1 Radio data

In our analysis below (Section 3), we combine data from the 2016 outburst with previous radio–X-ray monitoring campaigns of Aql X-1. We used the results of all published radio continuum observations that were taken before 2011 (table 1 from Tudose et al. 2009 and the entire⁶ table 1 from Miller-Jones et al. 2010a). From this set we selected only radio observations that were taken within 1 d of an RXTE-PCA X-ray observation.

2.2.2 X-ray data

From the RXTE-PCA data archive we selected the closest observation to each selected radio point. In most cases it was possible to bracket the radio observation with an additional X-ray observation that was within 2 d of the radio measurement, and which allowed us to interpolate the X-ray flux (see Section 3.1).

Selected X-ray observations were analysed using the standard heasoft tools. Each spectrum was produced by grouping 20 photons per bin, subtracting background, applying telescope response corrections and – in the case of seven observations – removing Type-I X-ray bursts. In order to correct for the systematic uncertainty of the RXTE-PCA (Jahoda et al. 2006) energy bins, we increased the uncertainties of each spectral bin by 0.5 per cent (following e.g. Nowak et al. 2011). We performed basic spectral fitting with xspec using three different possible models: a power-law, a power-law plus a Gaussian (if the iron Kα line was present), and a power-law plus a black-body component (if an additional thermal component was needed to produce an acceptable fit). For this spectral fitting we only used the 3–16 keV band; thus we kept the N_H value fixed to the value quoted above. For the neutral iron line modelling we restricted the mean of the Gaussian to the range 6.3–6.5 keV (since the iron line appears at 6.4 keV). We achieved acceptable fits for all observations (⁠|$\chi ^2_{\nu } \lt 1.3$| and p-value >0.08). We then used the cflux model to obtain Aql X-1’s unabsorbed 1–10 keV flux. We also obtained the unabsorbed 2–16, 6–9.7, and 9.7–16 keV fluxes from each spectrum in order to calculate the intensity and hardness ratio (HR) in the same way that it was defined in Miller-Jones et al. (2010a). The resulting fluxes and hardness ratios for each observation are presented in Table A1. We repeated this analysis for all RXTE-PCA observations in all seven previous outbursts discussed in this paper. The associated light-curves and HIDs for all previous outbursts are presented in the Appendix (Figs A1–A7).

3 ANALYSIS AND RESULTS

In the following analysis, we assume a distance of 4.5 kpc – as was done in Campana et al. (2014) based on the results of Galloway et al. (2008), where the peak luminosity of Type-I X-ray bursts was used to estimate distance.

3.1 Radio and X-ray luminosities

All three of our VLA observations of Aql X-1 resulted in non-detections. Using the imfit task in casa and the exact position of the source (taken from Miller-Jones et al. 2010a), we measured flux densities of 2.5 ± 4.2, −1.9 ± 4.3, and −0.2 ± 2.0 μJy beam⁻¹ at the three epochs, respectively, where the error corresponds to the off-source local rms image noise. These correspond to 3σ upper limits of 15, 11, and 6 μJy beam⁻¹, respectively (Table 1).

Aql X-1 was detected in all Swift-XRT observations that are used in our study (Table 1). During the 2016 outburst, Aql X-1 had an X-ray luminosity typical of NS-LMXBs in outburst (10³⁶⁻³⁸ erg s⁻¹). None the less, we note that the 2016 outburst was a particularly high-luminosity outburst – in fact the brightest observed from Aql X-1 to date – reaching a peak X-ray luminosity of almost 10³⁸ erg s⁻¹ (Güngör et al. 2017). Our VLA observations were performed as the source was decaying in brightness, with the first VLA observation taken when Aql X-1 was at an X-ray luminosity of 3 × 10³⁵ erg s⁻¹ and the last one when it was at 6 × 10³⁴ erg s⁻¹.

Since radio and X-ray observations were not strictly simultaneous, we obtained the X-ray flux values corresponding to the 2016 VLA and ATCA radio observations by logarithmically interpolating between pairs of closely spaced Swift-XRT observations (see Fig. 1). In the same manner, we obtained approximate X-ray fluxes and hardness ratios (9.7–16 keV flux/6–9.7 keV flux) at the time of each archival radio observation (Section 2.2.1) using RXTE-PCA data (Section 2.2.2). We increased the uncertainties on the interpolated X-ray fluxes at the times of the radio epochs in order to account for potential short-timescale (<1 d) X-ray variability. Using RXTE-PCA data from all 7 archival outbursts of Aql X-1 (190 observations), we found that the average maximum fractional difference in fluxes of two adjacent X-ray observations (Δ_fl = |fl₁ − fl₂|/fl₁) is proportional to the time between observations (Δ_t in days) in the form: Δ_fl = 0.2 × Δ_t + 0.005. We therefore increased the uncertainty on each interpolated X-ray flux (typically by ∼6 times compared to the formal uncertainties) using this formula and the separation in time between the radio epoch and closest available X-ray observation. When the source was dim, the formal errors on the X-ray fluxes dominated.

The 1–10 keV X-ray fluxes were converted into luminosities using L_X = 4πD²S, where S is the observed flux and D is distance to the source. The 5 GHz radio luminosities were calculated using L_R = 4πνD²S_ν, where ν is the central frequency and S_ν is the observed flux density. To incorporate radio measurements at different frequencies (∼4–12 GHz) we assumed a flat spectral index, which is consistent with previous multiband radio measurements of Aql X-1 in the hard X-ray state (Miller-Jones et al. 2010a; Díaz Trigo et al. 2018).

A summary of all 2016 and archival radio–X-ray measurements used in this study is provided in Table A1. The figures of the Appendix (Figs A1–A7) show light curves, hardness-intensity diagrams and L_R–L_X diagrams for each of the archival outbursts individually.

3.2 Radio–X-ray correlation

3.2.1 Fitting method

As in Gallo et al. (2018), we calculated the median values of β, ξ, and σ₀ from 10 000 draws of the posterior distribution, and determined 1σ confidence intervals using the 16–84th percentiles of the posterior distributions. Additionally, we performed 500 realizations of this method in order to confirm that the results were robust. The quoted values in Table 2 are the average of these realizations.

3.2.2 Selection criteria and fit results

Differences in the criteria used to select what data points are included in the fitting procedure can lead to significant differences in the resulting inferred power-law index. Such selection effects are, in fact, responsible for the difference in power law indices previously reported in the literature: β = 0.40 ± 0.07 (Tudose et al. 2009) and |$\beta = 0.76^{+0.14}_{-0.15}$| (Tetarenko et al. 2016). Tudose et al. (2009) used radio–X-ray observations of Aql X-1 performed during its 2002, 2004, and 2005 outbursts, including all radio observations, regardless of X-ray state. Tetarenko et al. (2016) used data points from the 2004, 2006, 2007, and 2009 outbursts of Aql X-1, but excluded radio observations that were performed in the soft X-ray state (HR <0.75 according to their definition, which is the same as used here). In both cases the authors did not include radio upper limits when performing their fits.

For our analysis, we used data from all eight outbursts for which quasi-simultaneous radio – X-ray data are available (see Fig. 2, Table A1, and Figs A1–A7). According to the classification of Campana et al. (2014), the 2002, 2004, 2005, 2006, and 2016 outbursts are FRED-like; the May 2007 outburst is multipeaked; and the October 2007 and 2009 outbursts are Gaussian-like. Alternatively, Güngör et al. (2017) classifies all these outbursts as low-luminosity, with the exception of the 2016 outburst. As in previous studies, we assume that – regardless of the outburst classification – the L_R– L_X correlation behaves similarly between outbursts. We do this because of the relatively low number of available observations, and argue that this is justified by the fact that these measurements occupy the same region in the L_R–L_X diagram (within the scatter), without obvious systematic offsets between outbursts (Fig. 2).

We performed fits using three different HR thresholds (HR >0.54, 0.75, and 0.96), both including and excluding upper limits, as well as with and without the 2016 outburst – see Table 2 for a summary. The data points corresponding to the various selection criteria are shown in Fig. 3.

The inferred slope changed with HR threshold, but within the uncertainties (Fig. 3 and Table 2). To allow direct comparison with previous studies, we adopt the threshold value of HR >0.75. The inclusion of upper limits produces a steeper slope compared to that previously reported in the literature, even if the 2016 outburst is left out. The inclusion of the 2016 outburst also leads to a consistent steep power-law index, which is consistent with the strong VLA radio upper limits at relatively low X-ray luminosity. In summary, our preferred fit uses HR >0.75, and includes upper limits as well as the 2016 outburst data.

One can also consider the constraints on the L_R– L_X correlation that can be derived using the 2016 outburst measurements alone. Unfortunately, given the relatively small number of hard-state observations (four in total with HR >0.75), the linmix_err method does not converge. Alternatively, we fit the 2016 outburst data points using a simple least-squares method, in which the upper limits were treated as 3σ detections. This resulted in a best-fitting slope β = 0.82, which is roughly speaking a 3σ lower limit on the steepness of the slope. This indicates that the 2016 data are consistent with the slope inferred from modelling all outbursts together.

4 DISCUSSION

Using the VLA, we have performed high-sensitivity, quasi-simultaneous radio – X-ray observations of Aql X-1 that probe an order of magnitude deeper in radio luminosity compared to previous studies. This allows us to meaningfully constrain the L_R–L_X correlation of Aql X-1 over 2 dex, down to low X-ray luminosities (6 × 10³⁴ erg s⁻¹ < L_X < 3 × 10³⁵ erg s⁻¹), as the source was fading from outburst.

ATCA detections of radio emission from Aql X-1 during the same outburst indicate that the jet was present throughout the outburst when L_X >5 × 10³⁶ erg s⁻¹ (Díaz Trigo et al. 2018). The last ATCA detection (with HR =0.64) was performed just 6 d before our first VLA observation. Thus our VLA non-detections are indeed indicative of a fading jet.

Here we discuss the derived constraints on the L_R– L_X correlation for Aql X-1, compare its behaviour to other NS-LMXBs, and discuss the motivation for future radio–X-ray observations.

4.1 L_R–L_X correlation in Aql X-1

Aql X-1’s jet behaviour during the 2016 outburst was possibly somewhat unusual: unlike what was previously observed (Miller-Jones et al. 2010a), the radio emission was not seen to quench at X-ray luminosities above 10³⁷ erg s⁻¹ (in the soft X-ray state), and rather reached a record brightness of 810 µJy at 5.5 GHz. However, the sparse observational sampling may mean that the quenching period was simply missed. Noteworthy is that the 2016 outburst is the only high-luminosity outburst (reaching about factor of 5 higher X-ray luminosity; see Güngör et al. 2014) of Aql X-1 for which both radio and X-ray observations were performed, and thus they probed unprecedentedly high X-ray luminosities. None the less, and as discussed before, the radio luminosities of the two observations performed during the (relatively) hard state (at L_X∼ 10³⁶⁻³⁷ erg s⁻¹) are consistent with previous observations from past outbursts. This suggests that it is appropriate to combine the 2016 outburst measurements with those of previous outbursts to study the L_R–L_X correlation of Aql X-1 using the largest available data set.

4.1.1 Modelling L_R–L_X with a single power law

We collected all archival radio and X-ray measurements prior to 2016 and fit these with a slope β = |$0.94^{+0.38}_{-0.28}$|⁠. As we discuss in Section 3 and summarize in Table 2, the slope is highly dependent on the inclusion of radio upper limits. While we exclude soft-state measurements, the slope depends only slightly on the chosen hardness ratio threshold. Ultimately, we find a slope that is consistent within uncertainties with the |$\beta = 0.76^{+0.14}_{-0.15}$| of Tetarenko et al. (2016), despite the fact they used data from fewer outbursts and only included detections in their fitting. Based on this slope, Tetarenko et al. (2016) suggested that Aql X-1 is consistent with being radiatively inefficient, and similar to BH-LMXBs. However, it is hard to draw robust conclusions from the previously available data because of the limited number of measurements and the fact that they span only ∼1 dex in radio and X-ray luminosity.

The 2016 measurements probe ∼2 dex in luminosity and provide significant new constraints. Including the 2016 outburst and all archival data, the resulting slope is β = |$1.17^{+0.30}_{-0.21}$| for HR >0.75. This is steeper than what has been reported in the previous studies mentioned above, but also consistent with what we obtained using only archival data. Given the small number of available observations, it is not possible to perform a robust separate fit to the 2016 outburst data alone, but these data are consistent with the global fit. In the framework of models in which jet power is a fixed fraction of the accretion power, this newly inferred slope β = |$1.17^{+0.30}_{-0.21}$| indicates that Aql X-1 has a more radiatively efficient accretion inflow than what was previously reported.

Recently, Qiao & Liu (2019) also used pre-2016 archival data of Aql X-1 and found a shallow (β ∼ 0.4) correlation between the radio and X-ray luminosities. These results are comparable to ours if we omit radio upper limits in our fitting (see Table 2). However, as we demonstrate, the inclusion of upper limits results in a steeper β =|$0.94^{+0.38}_{-0.28}$|⁠. In their work, Qiao & Liu (2019) attempt to explain their fitted correlation of β ∼ 0.4 with a coupled advection-dominated accretion (ADAF)-jet model (Yuan & Cui 2005; Done, Gierliński & Kubota 2007), where they assume that the accretion flow is always radiatively efficient for NSs, but that the fraction of the jet power to the accretion power is not constant. They assume that β ∼ 0.4 correlation holds down to low X-ray luminosities (L_X ∼10³⁴⁻³⁶ erg s⁻¹) and predict the corresponding radio luminosities to be L_R ∼1–4 × 10²⁸ erg s⁻¹. Our VLA upper limits on radio luminosity are about an order of magnitude fainter than was predicted by Qiao & Liu (2019). Therefore, such a scenario is inconsistent with our new 2016 observations. Moreover, the β = |$0.94^{+0.38}_{-0.28}$| that we find for archival data (including upper limits) suggests that the scenario is also inconsistent with the observations from previous outbursts.

4.1.2 Modelling L_R–L_X with a low-luminosity radio cut-off

Instead of fitting all data with an assumed single power-law correlation, we can also consider a model in which the radio brightness suddenly drops below a critical X-ray luminosity. This could be the case if a minimal accretion rate is necessary to sustain a steady jet.

The position of our new VLA observations of Aql X-1 in the L_R–L_X diagram may indicate a sudden step in radio luminosity, rather than a smooth decay. Indeed, during 2016 outburst, the X-ray luminosity at the time of the last ATCA detection and the first VLA upper-limit is different by a factor of 14, while the radio luminosity drops by more than a factor of 24. Furthermore, our first VLA non-detection (with L_X ∼ 3 × 10³⁵ erg s⁻¹, L_R ≤ 1.5 × 10²⁷ erg s⁻¹) is different from the lowest X-ray luminosity archival data point (L_X ∼ 5.5 × 10³⁵ erg s⁻¹, L_R ≤ 1.4 × 10²⁸ erg s⁻¹) by approximately a factor of 2 in X-ray luminosity, but by more than an order of magnitude in radio luminosity. This again assumes that we can compare separate outbursts.

Given the above, it may be questionable to impose a single power law over the entire L_X regime between L_X ∼ 6 × 10³⁴ erg s⁻¹ and L_X ∼10³⁷ erg s⁻¹. Instead, it could be that Aql X-1 is switching between two physically different accretion modes around some critical X-ray luminosity, above which L_R–L_X follows a shallow β ∼ 0.8 correlation, and below which the jet switches off (or becomes undetectable at 5–10 GHz) or does not follow a power-law correlation at all over this broad luminosity range. Dense sampling of the outburst decay is needed to explore this possibility.

4.2 Comparison to other NS-LMXBs

We now compare Aql X-1 to other NS-LMXBs with quasi-simultaneous radio and X-ray observations (Fig. 4, and associated references).

The best-studied regime is at high X-ray luminosities (L_X > 10³⁶ erg s⁻¹), where over a dozen systems have been observed. These show a broad range of radio luminosities for comparable X-ray luminosity. Aql X-1 is one of the fainter systems, being e.g. ∼15× less bright compared to IGR J17591−2342 (Russell et al. 2018). The origin of this large spread in radio luminosity is unclear. Presumably it could depend on viewing geometry, the magnetic field strength of the neutron star, its spin rate, and the inclination of the magnetic axis with respect to the orbital angular momentum. Determining which of these (or other) physical parameters are most relevant requires additional observations and numerical simulations. To date, the larger the number of NS-LMXBs that have been studied, the larger the diversity in their L_R–L_X behaviour. Aql X-1 has shown intermittent pulsations for a period of only 150 s (Casella et al. 2008). In the (perhaps unlikely) situation that it shows pulsations for a longer period during a future outburst, simultaneous radio observations can determine whether this has a demonstrable effect on the observed L_R–L_X behaviour. The same argument applies to other intermittent AMXPs.

The situation is even more unclear at intermediate X-ray luminosities (10³⁴ erg s⁻¹ < L_X < 10³⁶ erg s⁻¹), where only five systems have been observed, and only two detected. Both detected systems are AMXPs (SAX J1808.4−3658 and IGR J00291−5934), which led Tudor et al. (2017) to contemplate why these are surprisingly bright and apparently launch jets at relatively low X-ray luminosity. Our strong radio upper limits in this X-ray luminosity range for Aql X-1 show that it is certainly dimmer than the two detected AMXPs and consistent with the non-detections of two non-pulsating NS-LMXB systems (1RXS J180408.9−342058 and EXO 1745−248). However, given the large spread in radio luminosity at higher X-ray luminosity, the interpretation is unclear. It is too early to state that AMXPs are brighter on average compared to non-pulsating NS-LMXBs. Indeed, Gallo et al. (2018) concluded that there is no strong evidence for a systematic difference between AMXPs and non-pulsating NS-LMXBs. Here too, many more NS-LMXBs need to be observed, though this is challenging given their low radio brightness and will require using the most sensitive recently commissioned (MeerKAT) and future (Square Kilometre Array) radio telescopes.

Finally, at low X-ray luminosities (L_X < 10³⁴ erg s⁻¹), typically referred to as ‘quiescent’ state, the only radio-detected NS systems are the tMSPs PSR J1023+0038 and XSS J12270−4859. As previously discussed (Section 1), their radio–X-ray behaviour is demonstrably different compared to NS-LMXBs observed at higher accretion rate, and they are likely in a fundamentally different accretion regime. It would be very interesting to know whether Aql X-1 enters a similar regime during quiescence; however, given that it is at least three times further away compared to PSR J1023+0038, its expected radio flux density may be only ∼10 µJy. PSR J1023+0038 showed anticorrelated radio–X-ray brightness variations in which the radio flux density changed by a factor of ∼3 between X-ray high and low modes. A similar, strictly simultaneous radio–X-ray campaign on Aql X-1 during quiescence would likely require the sensitivity of future instruments like the Square Kilometre Array and Athena in order to be successful. However, Cen X-4 is similarly nearby to PSR J1023+0038 and remains undetected in radio (Tudor et al. 2017). This suggests that not all NS-LMXBs enter a low-luminosity accretion regime similar to that observed in the tMSPs. Possible explanations include the neutron star magnetic field strength or spin rate.

Despite the complications described above, there have been attempts to describe the population of NS-LMXBs with a single L_R–L_X relation. Gallo et al. (2018) studied the L_R–L_X correlation of all NS-LMXBs combined as well as for different NS-LMXB classes. They found a correlation for the combined sample of NS-LMXBs with a slope |$\beta = 0.44^{+0.05}_{-0.04}$| and |$\beta = 0.71^{+0.11}_{-0.09}$| for only atoll-type NS-LMXBs, both suggestive of radiatively inefficient accretion. Surprisingly, this result does not match what was so far observed in individual atoll-type NS-LMXBS for which the L_R–L_X correlation was studied in detail. As mentioned earlier, power-law slopes of β = 1.5 ± 0.2, |$\beta = 1.68^{+0.10}_{-0.09}$|⁠, and β =|$1.17^{+0.30}_{-0.21}$| were found for 4U 1728−34, EXO 1745−248, and Aql X-1, respectively – all consistent with radiatively efficient accretion. This argues that differences in the behaviour of individual sources – such as inconsistent L_R– L_X relationships (Tudor et al. 2017) or radio luminosity offset (Tetarenko et al. 2016) – are significant.

4.3 Conclusions

In summary, we have acquired quasi-simultaneous radio–X-ray observations that allow us to provide the strongest-ever constraints on the radio brightness of Aql X-1 at low X-ray luminosities (L_X < 10³⁶ erg s⁻¹). We combine these new measurements with an exhaustive re-analysis of all available radio–X-ray data from seven previous outbursts of Aql X-1. When modelling this combined data set with a single power law, we find a slope of |$1.17^{+0.30}_{-0.21}$|⁠, which is steeper than previously reported results in the literature and is consistent with radiatively efficient accretion. The inclusion of radio upper limits, which were ignored in the fits of previous studies, already lead to a steep inferred slope (⁠|$0.94^{+0.38}_{-0.28}$|⁠) even if the 2016 outburst is excluded.

ACKNOWLEDGEMENTS

NVG acknowledges funding from Netherlands Research School for Astronomy (NOVA). JWTH acknowledges funding from an Netherlands Organisation for Scientific Research (NWO) Vidi fellowship and from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013) / ERC Starting Grant agreement nr. 337062 (‘DRAGNET’). The Australia Telescope Compact Array is part of the Australia Telescope National Facility, which is funded by the Australian Government for operation as a National Facility managed by Commonwealth Scientific and Industrial Research (CSIRO). JCAM-J is the recipient of an Australian Research Council Future Fellowship (FT140101082), funded by the Australian government. ND is supported by a Vidi grant from NWO.

Footnotes

REFERENCES

APPENDIX A: