Radio Variability from a Quiescent Stellar-mass Black Hole Jet

The vast majority of data from the historical VLA were obtained through campaigns led by either Robert Hjellming or Michael Rupen. All observations were taken in continuum mode using two 50 MHz wide spectral windows. We reduced historical VLA observations following standard procedures within the Astronomical Image Processing System (AIPS) version 31DEC14 (Greisen 2003). We set the flux scale (using either 3C 48 or 3C 286) with the task setjy and (time-dependent) Perley & Butler (2013) coefficients. The complex gain solutions were then solved using scans of a secondary point-source calibrator (usually J2025+337, although see Table 2 for exceptions), and the flux scale was bootstrapped from the primary calibrator using the task getjy. Primary flux calibrator scans were not included for 24 observations. For those epochs, we manually set the flux scale to the expected value of the secondary phase calibrator by interpolating the flux densities reported by getjy of time-adjacent observations of the same calibrator at the same frequency. In these cases, we add a systematic uncertainty based on the level of variability of the secondary calibrator in nearby epochs, typically ≈5%. Finally, we add 5% and 10% systematic uncertainties to the flux scale for observations below and above 10 GHz, respectively.

The data were imaged using the task imagr, using Briggs weighting with a robust value of zero to help minimize sidelobes from nearby sources in the field. Of particular note is that a bright Jansky-level source lies 166 southeast of V404 Cygni (J2025+337, which was often used as the phase calibrator). This source is difficult to deconvolve during the imaging process because of bandwidth smearing (given the 50 MHz spectral windows of the historical VLA). We visually inspected a random subset of images of V404 Cygni over a range of frequencies and array configurations. We find that sidelobe artifacts from J2025+337 do not reach V404 Cygni at observing frequencies >8 GHz. At frequencies <8 GHz, there are multiple instances where artifacts from J2025+337 appear to increase the noise level near the location of V404 Cygni, but there are no cases where those artifacts obviously bias the measured flux densities.

Flux densities were measured using the task jmfit by fitting a two-dimensional Gaussian (fixed to the width of the synthesized beam) at the known location of V404 Cygni. Given potential systematics raised by J2025+337 at lower frequencies, we require detections to display peak flux densities at >5σ_rms, where σ_rms is the rms noise measured in a blank region of the sky (with upper limits for nondetections calculated as 5σ_rms).

Our sample also includes five epochs with the upgraded VLA. The first was taken in 2013 (Rana et al. 2016), and the other four were taken at the end of the 2015 outburst after the system reentered quiescence (Plotkin et al. 2017).

2.2.1. 2013 Quiescence

2.2.2. 2015 Post-outburst

We monitored V404 Cygni with the VLBA over 13 (approximately) fortnightly observations between 2014 February 3 and August 22 under program code BM399 (PI: Miller-Jones). Each observation lasted 2 hr, yielding ≈56 minutes on source. We observed with all available antennas, using 256 MHz of bandwidth centered on a frequency of 4.98 GHz. We used J2025+337 (Ma et al. 1998; 166 from V404 Cygni) as both the fringe finder and phase reference calibrator and the somewhat more distant J2023+3153 (Ma et al. 1998; 199 from V404 Cygni) as an astrometric check source. We also identified an additional VLBA observation at 8.4 GHz reported by Miller-Jones et al. (2009) taken on 2008 November 17 under program code BM290 in dual circular polarization with 64 MHz bandwidth per polarization (133 minutes on source), which we rereduced.⁷

We calibrated the data using AIPS (version 31DEC15) following standard procedures. For the 8.4 GHz observation, we used geodetic blocks at the start and end of the observation to correct for unmodeled clock errors and tropospheric delays. For all observations, we applied updated Earth orientation parameters, and we corrected for ionospheric dispersive delays using total electron content maps. We used system temperature information for amplitude calibration, corrected the phases for parallactic angle effects as the antenna feeds rotated with respect to the sky, and iteratively imaged and self-calibrated the phase reference source J2025+337 to derive a model for calibrating the phase, delay, and rate solutions, which were then applied to the target. Here V404 Cygni was too weak for self-calibration. After imaging with natural weighting (for maximum sensitivity), we determined the source flux density by fitting a point source in the image plane. V404 Cygni remains point-like at VLBA resolutions in quiescence (Miller-Jones et al. 2008). Therefore, we do not expect to resolve out any jetted radio emission, allowing a fair flux density comparison between the VLA and VLBA epochs.

Light curves are shown in Figure 1 at 4.9, 8.4, and 14.9 GHz (we omit our other two frequencies, 1.4 and 22.5 GHz, as they each contain only four observations). To quantify the distribution of flux densities at each frequency in the presence of nondetections, we perform a survival analysis. Using survfit in R (within the survival package⁸ ), we calculate the survival function, $S(\mathrm{log}\,{f}_{\nu })$, via the Kaplan–Meier estimator (see, e.g., Feigelson & Nelson 1985 for a description of the Kaplan–Meier estimator and examples of astrophysical applications). We then estimate the cumulative distribution function as $P(\mathrm{log}\,{f}_{\nu })=1-S(\mathrm{log}\,{f}_{\nu })$, which we display in Figure 2 at 4.9 and 8.4 GHz (omitting the four observations from 2015). In Section 3.3, we image the 2013 observation with the upgraded VLA at four separate frequencies centered at 5.0, 5.5, 7.2, and 7.7 GHz (each using 512 MHz bandwidth) to measure a spectral index. We therefore include the 5.0 GHz flux density from 2013 in the 4.9 GHz distribution here and the 7.7 GHz flux density in the 8.4 GHz distribution. In total, our 4.9 and 8.4 GHz distributions include 38 and 86 data points, respectively (of which 10 and 25 are upper limits, respectively). Detections (at the >5σ_rms level) range between 0.14 and 1.35 mJy bm⁻¹.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We compare the cumulative distribution functions in Figure 2 to lognormal distributions using the Peto & Peto modification of the Gehan–Wilcoxon test, as implemented by cendiff in the R package Nondetects and Data Analysis for Environmental Data (NADA⁹ ). The distributions of flux densities from V404 Cygni are not statistically different from lognormal distributions (p = 0.81 and 0.93 at 4.9 and 8.4 GHz, respectively) with $\left\langle \mathrm{log}\,{f}_{4.9}/\mathrm{mJy}\right\rangle =-0.53\pm 0.19$ and $\left\langle \mathrm{log}\,{f}_{8.4}/\mathrm{mJy}\right\rangle =-0.53\pm 0.30$, where f_4.9 and f_8.4 are flux densities at 4.9 and 8.4 GHz, respectively, in units of mJy. The quoted errors represent standard deviations on the lognormal distribution (i.e., they are not errors on the mean).¹⁰

We find that V404 Cygni shows significant flux density variations that are in excess of statistical fluctuations from measurement errors (which are typically ±0.04 mJy bm⁻¹ with the historical VLA). Taking 8.4 GHz as an example, among 61 detections, we find a reduced ${\chi }_{r}^{2}$ of 46 (for 60 degrees of freedom, as compared to a model with constant flux density). We note that ${\chi }_{r}^{2}$ is slightly biased by one observation from 2013 with a higher signal-to-noise ratio (S/N) using the upgraded VLA. Removing that observation still suggests significant intrinsic variability (${\chi }_{r}^{2}=33$ for 59 degrees of freedom). The fractional rms variability for all 61 observations at 8.4 GHz is F_var = 54% ± 6% (Vaughan et al. 2003), which is not biased by the higher S/N observations (i.e., we also calculate F_var = 54% ± 6% if we exclude the 2013 observation). If we consider the statistical fluctuations induced by variability to have a standard deviation of ${\sigma }_{{f}_{\nu },\mathrm{var}}=\pm {f}_{\nu }{F}_{\mathrm{var}}$, then propagating errors would yield ${\sigma }_{\mathrm{log}{f}_{\nu },\mathrm{var}}=\pm {F}_{\mathrm{var}}/\mathrm{ln}10$, such that a 54% fractional rms variability translates to 0.23 dex in logarithmic space.

To further quantify the flux variability, we produce the first-order structure function V(τ), which characterizes the amount of variability power as a function of timescale in the case of irregularly sampled data,

$$\begin{eqnarray}&&V(\tau )=\left\langle {\left[f\left(t+\tau \right)-f\left(t\right)\right]}^{2}\right\rangle ,\end{eqnarray} \tag{ 1 }$$

where f(t) is a flux density measurement at time t, and τ is a time delay. In the following analysis, we include only the 61 detections at 8.4 GHz. For every pair of data points (t_i, t_j) in our light curve, we calculate ${V}_{{ij}}({\tau }_{{ij}})={\left[f\left({t}_{j}\right))-f\left({t}_{i}\right)\right]}^{2}$, where τ_ij = t_j − t_i. We then bin our set of V_ij measures by time difference (τ_ij) so that each bin contains 50 data points, and we take the average of those 50 measurements to calculate V(τ) (where we adopt the midpoint of all time differences within each bin as the value for the time delay τ).¹¹ The structure function is shown in Figure 3, where we probe long-term variability over timescales of ∼10–4000 days.

Zoom In Zoom Out Reset image size

The slope of the structure function provides information on the power distribution of flux variations (see, e.g., Hughes et al. 1992 for details, which we summarize below). For example, if variability is characterized as flicker noise (i.e., a power function P(F) ∝ F⁻¹, where F is the inverse timescale, i.e., frequency, of fluctuations), then the structure function V(τ) = constant. The constant is expected to be $2{\sigma }_{\mathrm{var}}^{2}$, with ${\sigma }_{\mathrm{var}}^{2}$ being the variance of the observed flux densities. One could also obtain $V(\tau )=2{\sigma }_{\mathrm{var}}^{2}$ on long timescales if the structure function probes white noise (e.g., one interpretation would be that the probed timescales are longer than the characteristic timescale on which shot-noise variations dampen). As another example, red-noise variations (P(F) ∝ F⁻²) produce first-order structure functions of the form V(τ) ∝ τ, which can often be interpreted as disturbances that follow a random-walk process (e.g., Kelly et al. 2009). Finally, flux variations smaller than the level of a typical error bar on flux density measurements (σ_err) are not meaningful (assuming that σ_err is dominated by statistical noise). Thus, one expects the structure function to satisfy $V(\tau )\geqslant 2{\sigma }_{\mathrm{err}}^{2}$ at all timescales.

The structure function in Figure 3 appears flat, with a best-fit slope and normalization of β = −0.11 ± 0.07 and V₀ = 0.16 ± 0.04 mJy² (where V(τ) = V₀τ^β). Measuring the variance directly from the 61 data points yields $2{\sigma }_{\mathrm{var}}^{2}\,=0.11\,{\mathrm{mJy}}^{2}$. Therefore, the structure function appears consistent with plateauing near $2{\sigma }_{\mathrm{var}}^{2}$, which signifies either flicker or white noise. In the latter case, we would be probing jet disturbances on timescales (≳10 days) that are longer than the characteristic damping timescales.

We also examine variability on sub-hour timescales, focusing on observations long enough (usually ≳90 minutes on source) to produce light curves over multiple time bins (we are able to achieve time resolutions ranging from 3 to 30 minutes; see Figure 4). We focus on observations near 8 GHz, where we have 11 long observations total, including seven from the historical VLA, one from the VLBA, and three from the upgraded VLA.

Zoom In Zoom Out Reset image size

We characterize the level of variability on each epoch by calculating the reduced ${\chi }_{r}^{2}$ of the flux densities in each light curve. For this calculation, we ignore upper limits (which causes us to underestimate ${\chi }_{r}^{2}$ in some cases). Our measured ${\chi }_{r}^{2}$ values are listed in Table 3. Treating ${\chi }_{r}^{2}\gt 3$ as a crude diagnostic for significant intrinsic variations, we find obvious variability on 8/11 epochs. However, we cannot determine from ${\chi }_{r}^{2}$ alone if we should consider the other three epochs as nonvariable. Rather, lower ${\chi }_{r}^{2}$ values imply that variability on those epochs is less extreme relative to fluctuations from statistical noise, and further tests are required (see next paragraph).

Table 3.  Short-term Variability Statistics

Date	${\left({\chi }_{r}^{2}\right)}_{\mathrm{flux}}$	Degree of Freedom	β
(1)	(2)	(3)	(4)
1998 May 4	11.9	11	⋯
1998 Oct 25	9.0	14	⋯
2000 Jul 21	8.2	38	⋯
2003 Jul 29	2.1	26	⋯
2007 Dec 2	8.5	16	⋯
2008 Nov 17	35.4	6	⋯
2009 Feb 15	1.5	13	⋯
2009 Apr 26	3.7	23	⋯
2013 Dec 2	7.7	52	0.29 ± 0.05
2015 Aug 1	2.8	32	0.10 ± 0.11
2015 Aug 5	12.3	32	1.08 ± 0.08

Note. Column (1): calendar date of each observation. Columns (2)–(3): reduced ${\chi }_{r}^{2}$ for each light curve (including only detections) and the number of degrees of freedom. Column (4): the best-fit power-law index to the structure function, when available (V(τ) ∝ τ^β).

Download table as:  ASCIITypeset image

The three epochs with the upgraded VLA each display different short-term variability characteristics (i.e., the flux density is slowly decreasing with time on 2013 December 2, there are no obvious flares on >3 minute timescales on 2015 August 1, and we observe the beginning of a flare on 2015 August 5). These three epochs therefore provide useful illustrative examples, and as such, we display their structure functions in Figure 5. We fit a power law to each structure function, and we find power-law indices of β = 0.29 ± 0.05, 0.10 ± 0.11, and 1.08 ± 0.08 on 2013 December 2, 2015 August 1, and 2015 August 5, respectively (V(τ) ∝ τ^β). Note that the structure function on 2015 August 1 is flat, plateauing at $2{\sigma }_{\mathrm{var}}^{2}$, indicating that V404 Cygni indeed shows (uncorrelated) variability in excess of statistical measurement noise fluctuations on 2015 August 1, even though ${\chi }_{r}^{2}=2.8$.

Zoom In Zoom Out Reset image size

3.2.1. Quiescent Flares

For 24 epochs, there are observations at multiple frequencies within ±30 minutes from which we measure spectral indices (${f}_{\nu }\propto {\nu }^{{\alpha }_{r}}$), as shown in Figure 6. For epochs with exactly two frequencies, we calculate α_r (or place a limit) analytically as ${\alpha }_{r}=\mathrm{ln}({f}_{{\nu }_{1}}/{f}_{{\nu }_{2}})/({\nu }_{1}/{\nu }_{2})$, where ${f}_{{\nu }_{1}}$ and ${f}_{{\nu }_{2}}$ refer to flux densities at frequencies ν₁ and ν₂, respectively. We assign an error on α_r by propagating through statistical errors (we ignore errors on frequency for the historical VLA given the relatively small bandwidth at each frequency, and we also ignore systematic errors from short-term variability, since we generally do not know if the spectra were taken during flaring activity). If >2 observing frequencies, then we measure α_r through a least-squares fit to the spectrum (in log space). We estimate error bars through Monte Carlo simulations, where we randomly add noise to each flux density measurement (assuming Gaussian noise with a standard deviation equal to the measurement error on each data point) and then refit the spectral index. We repeat 1000 times and estimate ${\sigma }_{{\alpha }_{r}}$ as the standard deviation on the resulting α_r distribution.

Zoom In Zoom Out Reset image size

For the 2013 epoch from Rana et al. (2016) with the upgraded VLA, we create four (strictly simultaneous) images by splitting the bandwidth into four sub-bands with 512 MHz bandwidth (centered at 5.0, 5.5, 7.2, and 7.7 GHz), and we allow the frequency to randomly vary across each 512 MHz when running our Monte Carlo simulations to estimate error bars. We obtain α_r = −0.26 ± 0.05, which is consistent with the value α_r = −0.27 ± 0.03 obtained by Rana et al. (2016). For the four VLA observations from 2015, we adopt the spectral indices reported by Plotkin et al. (2017), which were calculated using the same least-squares fitting method as described above (over 4–12 GHz).

3.3.1. Long-term Spectral Variations

We measure a large range of spectral indices, $-1.5\,\lesssim {\alpha }_{r}\lesssim +1.6$, with the spread of α_r being larger at lower flux densities (Figure 7). However, we argue in Section 4.1 that this spread in α_r is likely attributed largely to statistical and systematic errors (i.e., large error bars on most α_r measurements, especially at lower flux densities, combined with only five epochs having strictly simultaneous multifrequency data).

Zoom In Zoom Out Reset image size

To quantify the dispersion in radio spectral indices, we fit a Gaussian distribution to the α_r measurements (using a Bayesian framework that allows for both upper and lower limits). We find a mean $\left\langle {\alpha }_{r}\right\rangle =0.02\pm 0.17$ and a standard deviation ${\sigma }_{{\alpha }_{r}}=0.65\pm 0.15$, where the errors represent 68% confidence intervals of the posterior distributions. If we exclude the five spectra obtained from the upgraded VLA (which have significantly smaller error bars on α_r), we find consistent results: $\left\langle {\alpha }_{r}\right\rangle =-0.07\pm 0.25$ and ${\sigma }_{{\alpha }_{r}}=0.78\pm 0.21$. Throughout, we adopt $\left\langle {\alpha }_{r}\right\rangle =0.02\pm 0.65$, which is consistent with a flat radio spectral index, on average, as expected for a partially self-absorbed compact synchrotron jet.

3.3.2. Short-term Spectral Variations

Among the 11 observations from which we produce sub-hour light curves in Section 3.2, four contain multifrequency data allowing us to explore variations in α_r over sub-hour timescales. These epochs include 2013 December 2 with the upgraded VLA (Rana et al. 2016, 4–8 GHz), two epochs at the end of the 2015 outburst on August 1 and 5 (Plotkin et al. 2017, 4–12 GHz), and a historical VLA observation on 2003 July 29 that interleaved observations at 4.9 and 8.4 GHz every ≈15 minutes (Hynes et al. 2004, 2009).

The radio spectra from these four epochs are displayed in Figure 8, where we also display the corresponding light curves for reference.¹² The spectral constraints during the 2003 epoch are not very meaningful, but we include that epoch in the figure for completeness. For the other three epochs, we do not see any obvious changes of α_r with time; variations of α_r are consistent with statistical noise, with reduced ${\chi }_{r}^{2}$ of 1.5, 2.5, and 1.5 (for 52, 32, and 32 degrees of freedom) on 2013 December 2, 2015 August 1, and 2015 August 5, respectively. We note that Rana et al. (2016) reported that the radio spectrum of V404 Cygni switched from optically thick to optically thin over ≈10 minute periods on 2013 December 2. While we also see some variations in α_r over 10–30 minute timescales, they tend to be at the 1σ–2σ level, such that we do not consider those variations to be highly significant relative to the measurement error. We tentatively see marginal indications for a long-term evolution of the spectrum over the ≈9 hr observation. Imaging the first 100 minutes of the observation yields α_r = −0.39 ± 0.10, while a less steep (and potentially flat) spectral index of α_r = −0.14 ± 0.14 is measured from images during the last 100 minutes. The steeper spectral index during the first 100 minutes appears to be driven by variations at 7.45 GHz that are not mimicked at 5.25 GHz.

Zoom In Zoom Out Reset image size

As noted in Section 3.3.1, we measure radio spectral indices in the range −1.5 < α_r < +1.6, with a mean α_r = 0.02 and a standard deviation ${\sigma }_{{\alpha }_{r}}=0.65$. We suspect that the large range in α_r can be attributed to the poorer sensitivity of the historical VLA and multifrequency data that can be offset by ±30 minutes. For example, a factor of 2 variability within 30 minutes would adjust α_r by ±1.3. We are therefore cautious not to overinterpret the apparently large range of measured α_r, which might not require a physical explanation.

The above variability-related concerns, however, do not preclude the possibility of less extreme spectral variations. The upgraded VLA provides strictly simultaneous multifrequency observations, and the 2013 December 2 epoch (α_r = −0.26 ± 0.05) provides evidence that the spectral index can indeed stray negative at times (as originally pointed out by Rana et al. 2016). We stress that α_r = −0.26 ± 0.05 is inconsistent with a purely optically thin spectrum, as expected from transient ejecta (e.g., Fender et al. 1999; Corbel et al. 2004). Possible explanations for a mildly negative spectral index, as observed on 2013 December 2, could include the combination of optically thick and thin emission (e.g., the fading optically thin stage of a flare superposed over a steady, flat-spectrum jet), a jet that expands more slowly than a conical jet, or a decelerating jet. Also contributing could be that lower-frequency variations are "smeared," since they are emitted over a larger volume (farther from the black hole) compared to higher-frequency emission. Indeed, the 2013 December 2 light curve (Figure 8) displays dips at 7.45 GHz approximately 50 and 150 minutes after the start of the observation that are not mimicked at 5.25 GHz, which likely contributes to the overall negative spectral index.

To our knowledge, the 2013 December 2 observation of V404 Cygni is the only observation of a quiescent BHXB to show a (well-measured) negative spectral index. However, comparably negative spectral indices have been measured from compact jets in the hard state (see, e.g., Espinasse & Fender 2018 and references therein, although we note that many of those observations also suffer from low sensitivity). Reasonable explanations for different spectral indices in the hard state include intrinsic differences in jet properties (Espinasse & Fender 2018), differences in inclination (Motta et al. 2018), or both. One difference in our work, however, is that we are seeing both positive and negative spectral indices from the same source. Therefore, we cannot explain a varying radio spectral index from V404 Cygni in quiescence as an inclination effect.

Bernardini & Cackett (2014) obtained dense X-ray monitoring of V404 Cygni with the Neil Gehrels Swift Observatory, taking 33 observations over 75 days. Their study provides the most appropriate comparison to our long-term radio light curve(s) in Figure 1. They find a fractional rms variability in the X-ray of ${F}_{\mathrm{var},{\rm{X}} \mbox{-} \mathrm{ray}}=57 \% \pm 3$% (∼0.25 dex), and they obtain a flat structure function over timescales of ≈5–80 days. These X-ray results are similar to our radio results, where we find F_var,radio = 54% ± 6% (∼0.23 dex) and a flat structure function (although our study extends to longer timescales of of ≈10–4000 days).

While the level of long-term X-ray variability is comparable to that in the radio, short-term X-ray variability can be larger in amplitude. For example, on hour timescales, factor of 4–8 X-ray variations are typical (e.g., Bradley et al. 2007; Bernardini & Cackett 2014; Rana et al. 2016); at the extreme end, Hynes et al. (2004) observed a flare that increased the Chandra count rate by a factor of >20 (also see Wagner et al. 1994 for a factor of 10 variation in <0.5 days). For comparison, in our work, we commonly observe flares that change the flux density by factors of 2–4 in the radio.

Whether or not one should expect correlated variations between the radio and X-ray bands on short timescales is an open question. To our knowledge, there have been no multiwavelength campaigns that simultaneously take radio and X-ray observations in quiescence over multiple epochs separated by only 1–2 days, thereby making it impossible to empirically test if radio/X-ray variations are correlated on sub-week timescales. On minute-through-hour timescales, only two attempts have been published so far that searched for coordinated radio/X-ray variability over observations lasting several hr (2003 July 29 and 2013 December 2), and neither attempt has shown obvious radio/X-ray correlations with light curves on 30–100 minute time bins (Hynes et al. 2009; Rana et al. 2016). However, detecting correlated variations may require finer time resolution.

Considering the above X-ray characteristics, and that correlated (short-term) radio and X-ray variability might only be detectable with minute time resolution, we suggest the following strategies when coordinating multiwavelength campaigns on quiescent X-ray binaries.

• 1.  
  Radio and X-ray observations should be scheduled as simultaneously as possible. If the source is bright enough to provide sufficient S/N on <10 minute time bins in both the radio and X-ray, then (some) individual flares should be resolvable, and one can directly control for multiwavelength variability. Ideally, the observations would be long enough to observe the entire flare rise and decay (which can last several hr).
• 2.  
  If data are nonsimultaneous, or if one does not have sufficient time resolution to resolve individual flares, then inflate the error bars by ≈0.25 dex in both the radio and the X-ray. However, one should still bear in mind that variations as large as factors of 2–4 in the radio and 4–8 in the X-ray are common.
• 3.  
  Radio observations should be taken as close as possible to the frequency required to achieve one's science goals. If spectral constraints from strictly simultaneous multifrequency data are lacking, then we support the standard assumption of a flat radio spectrum when extrapolating radio observations to other frequencies. However, we find evidence that α_r is not always strictly zero, and we recommend propagating errors on flux densities by assuming that the radio spectrum could vary by ${\sigma }_{{\alpha }_{r}}\approx \pm 0.6$ (see Section 3.3.1).

While quiescent BHXBs tend to have higher L_R/L_X ratios than other classes of Galactic compact accreting objects (e.g., Migliari & Fender 2006; Tudor et al. 2017; Gallo et al. 2018), the transitional millisecond pulsar (tMSP) PSR J1023+0038 (hereafter J1023) was recently shown to sometimes be an exception to that trend (Bogdanov et al. 2018). Considering periods when J1023 does not show radio pulsations as accretion-powered states, J1023 exhibits aperiodic and rapid switching between two distinct X-ray flux levels, a low mode (L_X ≈ 5 × 10³² erg s⁻¹) and a high mode (L_X ≈ 3 × 10³³ erg s⁻¹; e.g., Patruno et al. 2014; Stappers et al. 2014; Bogdanov et al. 2015; Jaodand et al. 2016). Bogdanov et al. (2018) discovered that the low modes are nearly always accompanied by radio flares that both rise and decay on minute timescales. The radio flux density appears anticorrelated with the X-ray flux, although some radio flaring is also observed when the X-ray flux remains in a high mode. During these low X-ray mode states, the increased radio flaring tends to reach L_R ≈ 2 × 10²⁷ erg s⁻¹.

In Figure 9, we highlight the location of J1023 in the L_R − L_X plane¹⁵ (large green triangles) compared to V404 Cygni in quiescence and the hard state (large blue circles). We also shade in a region that represents the minimum and maximum luminosities displayed by V404 Cygni within our radio data set (≈5 × 10²⁷–5 × 10²⁸ erg s⁻¹, which corresponds to flux densities ranging from 0.14 to 1.35 mJy if one assumes a flat radio spectrum) and the minimum and maximum X-ray luminosities displayed during the Swift campaign from Bernardini & Cackett (2014). We estimate 1–10 keV X-ray luminosities assuming a photon index of 2.1 and a column density of 9 × 10²¹ cm⁻² (Bernardini & Cackett 2014) using X-ray count rates extracted from the online Swift-XRT product generator tool¹⁶ (Evans et al. 2007, 2009). The shaded region in Figure 9 assumes that radio and X-ray luminosities are not correlated in quiescence and therefore represents a conservative range for where one might find V404 Cygni in the L_R − L_X plane if considering nonsimultaneous data. The mean radio flux density from our study ($\mathrm{log}\,{f}_{\nu }/\mathrm{mJy}=-0.53$) corresponds to $\mathrm{log}{L}_{{\rm{R}}}/\mathrm{erg}\,{{\rm{s}}}^{-1}\,\approx 28$, such that V404 Cygni likely spends most of its time toward the bottom left of the shaded region (with the top right region representing periods of flaring activity).

Zoom In Zoom Out Reset image size

After accounting for the range of flux density variations exhibited by V404 Cygni, we expect its radio luminosity to always be ≳2.5 times larger than the radio luminosity of J1023, even when J1023 is in a low X-ray flux (i.e., radio-flaring) state. Still, J1023 ventures very close to the parameter space expected for quiescent BHXBs, and it is easy to envision that it could overlap with a BHXB that has similar variability characteristics as V404 Cygni but at a slightly lower Eddington ratio. Although J1023 only ventures toward the BHXB parameter space on short timescales, ≲500 s (i.e., over longer time-averaged observations, J1023 will generally appear radio fainter), other tMSPs could venture toward the BHXB parameter space for longer periods of time (e.g., the tMSP IGR J18245−2452 has displayed at least one low mode that lasted for nearly 10 hr; e.g., Linares et al. 2014). We therefore reiterate the conclusion from Bogdanov et al. (2018) that radio and X-ray luminosities alone are not always sufficient to confidently assert that a quiescent X-ray binary contains a black hole instead of a neutron star.¹⁷ Other multiwavelength data should naturally also be consulted, e.g., to search for an orbital period, emission lines, or other properties to give clues on the nature of the companion star and accreting compact object (e.g., Bahramian et al. 2017; Shishkovsky et al. 2018; Tudor et al. 2018).

Lacking detailed multiwavelength data, we assert that radio light curves could have some diagnostic power to determine if one is observing a black hole or a neutron star.

• 1.  
  Frequent and relatively short radio flares that both rise and decay on rapid (minute) timescales may favor a neutron star. Bogdanov et al. (2018) found radio flares that last ≲500 s, and they rise and decay rapidly (see their Figures 1, 4, and 5). These flares can occur quite frequently at times (e.g., Bogdanov et al. 2018 observed as many as three to four radio flares per hour during some portions of their observation). For V404 Cygni, we observe less frequent flaring, and we only observe rapid flare rises. The decays tend to proceed on longer timescales of tens of minutes to hours.
• 2.  
  An order-of-magnitude change in the quiescent radio flux density over minute-to-hour timescales might argue for a neutron star. Deller et al. (2015) showed light curves of J1023 from 13 different epochs (see their Figure 5), and they found a case where the radio flux density steadily increases by an order of magnitude over ≈30 minutes. While V404 Cygni has also been observed to show flares with comparably slow rise times, none of those >30 minute flares have increased in flux density by an order of magnitude.

We stress that the above assertions are predicated on observing a system that is close enough to produce radio light curves on minute timescales (e.g., V404 Cygni is at 2.39 ± 0.14 kpc, Miller-Jones et al. 2009; and J1023 is at ${1.368}_{-0.039}^{+0.042}$ kpc, Deller et al. 2012). Also, in some cases, radio light curves will not hold any diagnostic power, since, particularly with shorter observations, J1023 and V404 Cygni can both show time periods of low activity. Furthermore, the variability observed so far from J1023 and V404 Cygni is unlikely to represent the full range of variability that can be displayed by either class of systems (e.g., as evidenced by the nearly 10 hr X-ray low mode displayed by the tMSP IGR J18245−2452).

If one does not have access to radio light curves with minute time resolution, then taking a radio flux density integrated over long observations could provide a less ambiguous interpretation: over long time intervals, nonflaring time periods will dilute the radio flux density, thereby moving the tMSP farther away from the quiescent BHXB space in L_R − L_X. Finally, in some cases, the X-ray spectrum can provide additional diagnostic power: J1023 has an X-ray photon index of Γ ≈ 1.7 in both low and high X-ray modes (Bogdanov et al. 2018), while quiescent BHXBs are well established to have a softer Γ ≈ 2.1 in quiescence (Plotkin et al. 2013; Reynolds et al. 2014). The radio spectral index, however, is a poor discriminant, given that both V404 Cygni and J1023 exhibit a range of mildly positive and negative spectral indices (Bogdanov et al. 2018 reported a range of −0.5 ≲ α_r ≲ 0.4) and that α_r is usually not well measured in low-luminosity sources.