Radio Variability from a Quiescent Stellar Mass Black Hole Jet

TL;DR

This study analyzes 24 years of radio observations of the quiescent black hole binary V404 Cygni, revealing frequent short-term flares and long-term variability, and discusses implications for understanding jet emission and classifying such systems.

Contribution

It provides a comprehensive characterization of radio variability in a quiescent stellar-mass black hole, highlighting the prevalence of rapid flares and long-term fluctuations.

Findings

Radio flux densities follow lognormal distributions.

Frequent radio flares increase flux by factors of 2-4 within minutes.

Long-term variations are consistent with shot noise impulses.

Abstract

Relativistic outflows are believed to be a common feature of black hole X-ray binaries at the lowest accretion rates, when they are in their `quiescent' spectral state. However, we still lack a detailed understanding of how quiescent jet emission varies with time. Here we present 24 years of archival radio observations (from the Very Large Array and the Very Long Baseline Array) of the black hole X-ray binary V404 Cygni in quiescence (totalling 150 observations from 1.4 -- 22 GHz). The observed flux densities follow lognormal distributions with means and standard deviations of (<log f_nu>, sigma) = (-0.53, 0.19) and (-0.53, 0.30) at 4.9 and 8.4 GHz, respectively (where f_nu is the flux density in units of mJy). As expected, the average radio spectrum is flat with a mean and standard deviation of (<alpha_r >, sigma_alpha_r)= (0.02, 0.65) where f_nu \propto nu^alpha_r. We find that radio flares that increase the flux density by factors of 2 -- 4 over timescales as short as <10 min are commonplace, and that long-term variations (over 10--4000 day timescales) are consistent with shot noise impulses that decay to stochastic variations on timescales <10 days (and perhaps as short as tens of minutes to several hours). We briefly compare the variability characteristics of V404 Cygni to jetted active galactic nuclei, and we conclude with recommendations on how to account for variability when placing quiescent black hole X-ray binary candidates with radio luminosities comparable to V404 Cygni (L_r ~ 1e28 erg/s) onto the radio/X-ray luminosity plane.

Tables

Table 1. Table 1: Number of Observations per Frequency

Observing Band	Telescope	$N_{\mathrm{obs}}$	$N_{\det }$
L /1.4 GHz	Historical VLA	4	1
C / 4.9 GHz	Historical VLA	24	14
$\cdots$	Upgraded VLA	5aaFour observations with the upgraded VLA post-2015 outburst were taken in subarray mode, with half of the antennas observing at C-band and the other half at X-band, which we count as separate C- and X-band observations in this table. The fifth upgraded VLA observation was taken at C-band (4-8 GHz) in 2013. We often add a 7.7 GHz flux density measurement from this C-band observation into our analysis of the X-band sample (which is not reflected within this table).	5
$\cdots$	VLBA	13	13
X / 8.4 GHz	Historical VLA	84bbThere are two additional observations in the archive at this frequency for which we could not obtain a calibration solution.	59
$\cdots$	Upgraded VLA	4aaFour observations with the upgraded VLA post-2015 outburst were taken in subarray mode, with half of the antennas observing at C-band and the other half at X-band, which we count as separate C- and X-band observations in this table. The fifth upgraded VLA observation was taken at C-band (4-8 GHz) in 2013. We often add a 7.7 GHz flux density measurement from this C-band observation into our analysis of the X-band sample (which is not reflected within this table).	4
$\cdots$	VLBA	1	1
KU / 14.9 GHz	Historical VLA	11	5
K / 22.5 GHz	Historical VLA	4	0

Table 2. Table 2: Catalog of Radio Observations

Date	MJD	Program ID	Configuration	${\tau }_{\mathrm{source}}$	Frequency	$f_{\nu }$	${\alpha }_r$	Primary	Secondary	PI
				(min)	(GHz)	(mJy)		Calibrator	Calibrator
1991 Sep 25	48525.012	AH424	BnA	36.8	4.9	0.238 ±0.042	1.59 ±0.36	3C286	2025+337	Han
$\cdots$	48525.012	AH424	BnA	38.7	8.4	0.560 ±0.046	$\cdots$	3C286	2025+337	Han
$\cdots$	48524.984	AH424	BnA	55.7	14.9	¡0.695	$\cdots$	3C286	2025+337	Han
1991 Oct 31	48561.022	AH390	BnA	59.8	4.9	0.321 ±0.028	-0.72 ±0.33	$\cdots$	2025+337	Hjellming
$\cdots$	48561.004	AH390	BnA	46.0	8.4	0.218 ±0.034	$\cdots$	$\cdots$	2025+337	Hjellming

Table 3. Table 3: Short-term Variability Statistics

Date	${\left({\chi }_r^2\right)}_{\mathrm{flux}}$	deg of freedom	$\beta$
1998 May 04	11.9	11	$\cdots$
1998 Oct 25	9.0	14	$\cdots$
2000 Jul 21	8.2	38	$\cdots$
2003 Jul 29	2.1	26	$\cdots$
2007 Dec 02	8.5	16	$\cdots$
2008 Nov 17	35.4	6	$\cdots$
2009 Feb 15	1.5	13	$\cdots$
2009 Apr 26	3.7	23	$\cdots$
2013 Dec 02	7.7	52	0.29 ±0.05
2015 Aug 01	2.8	32	0.10 ±0.11
2015 Aug 05	12.3	32	1.08 ±0.08

Equations

$$V(\tau )=⟨{\left[f\left(t+\tau \right)-f\left(t\right)\right]}^2⟩,$$

Full text

accretion, accretion disks — stars: black holes — stars:individual:V404 Cygni — X-rays: binaries

††facilities: VLA, VLBA††software: AIPS (v31DEC2014, v31DEC2015 Greisen, 2003), Astropy (Astropy Collaboration et al., 2013), CASA (v5.1.1 McMullin et al., 2007)

1 Introduction

Black hole X-ray binaries (BHXBs) that accrete at low Eddington ratios ($\lesssim 0.01L_{\rm Edd}$, where $L_{\rm Edd}$ is the Eddington luminosity) are associated with compact radio emission from partially self-absorbed synchrotron jets (Blandford & Königl, 1979; Fender, 2001; Remillard & McClintock, 2006). These relativistic outflows provide channels for accretion flows to shed angular momentum, and to transport energy out to large distances (e.g., Meier, 2001; Fender & Muñoz-Darias, 2016; Romero et al., 2017; Douna et al., 2018). Radio observations of relativistic jets are therefore crucial for understanding how matter is transported through, and away from, accretion flows with low mass accretion rates.

The most weakly accreting black holes reside in the ‘quiescent’ spectral state, which we define here by a soft X-ray spectrum that can be characterised by a power-law photon index111The photon index $\Gamma$ is defined by $N_{E}\propto E^{-\Gamma}$, where $N_{E}$ is the photon number density per unit energy, $E$. of $\Gamma\sim 2.1$ (e.g., Tomsick et al., 2001; Kong et al., 2002; Corbel et al., 2006; Plotkin et al., 2013; Reynolds et al., 2014). Most BHXBs spend the majority of their time in quiescence, which usually corresponds to (0.5-10 keV) X-ray luminosities $L_{\rm X}\lesssim 10^{-6}-10^{-5}L_{\rm Edd}$ (Plotkin et al., 2013, 2017). In quiescence, a larger fraction of the radiative power emitted by the accretion flow/jet system appears to be emitted in the radio waveband (e.g., Fender et al., 2003; Corbel et al., 2013; Gallo et al., 2018, although see Yuan & Cui 2005 for a prediction otherwise). This increased dominance of radio emission implies that the radio domain can be effective for discovering quiescent BHXBs (Maccarone, 2005; Fender et al., 2013), which would produce less biased samples of BHXBs in the Milky Way compared to the traditional method of discovering BHXBs through X-ray emission during an outburst.

Coordinated radio and X-ray surveys are indeed starting to reveal candidate BHXBs in Milky Way globular clusters (Strader et al., 2012; Chomiuk et al., 2013; Miller-Jones et al., 2015; Shishkovsky et al., 2018) and in the field (Tetarenko et al., 2016). Such an approach is highly complementary to other strategies for discovering quiescent BHXBs, through, e.g., H$\alpha$ surveys (Casares, 2018), X-ray surveys (e.g., Agol & Kamionkowski, 2002; Jonker et al., 2014), and optical spectroscopic searches capable of discovering non-accreting black hole candidates in detached binary systems (Giesers et al., 2018; Thompson et al., 2018).

At the moment, even our most sensitive radio facilities are only capable of probing the tip of the quiescent BHXB population. We have currently detected radio emission from only four nearby quiescent BHXBs ($\lesssim$4 kpc), including one of the most luminous known quiescent systems, V404 Cygni ($L_{\rm R}\approx 10^{28}\,{\rm erg~{}s}^{-1}$ at 5 GHz; Hynes et al. 2004; Gallo et al. 2005; Rana et al. 2016), and three of the least luminous known systems (A 0620$-$00, XTE J1118+480, and MWC 656; $L_{\rm R}\approx 10^{26}\,{\rm erg~{}s}^{-1}$; Gallo et al. 2006, 2014; Dzib et al. 2015; Ribó et al. 2017; Dinçer et al. 2018). Thus, we are still establishing the empirical properties of quiescent BHXB radio jets, which is an essential step for eventually defining reliable radio-based selection criteria.

To help inform future radio surveys, this paper focuses on quantifying the radio variability characteristics of quiescent jets. Understanding that variability level is important, as (i) it would establish whether radio variability can discriminate quiescent BHXBs from other classes of radio-emitting objects; and (ii) it would help determine how close in time radio observations must be coordinated with other multiwavelength data.

To open the quiescent radio time domain we focus on V404 Cygni, which represents the luminous end of the quiescent BHXB population. V404 Cygni contains a dynamically confirmed black hole ($9.0^{+0.2}_{-0.6}M_{\odot}$) orbiting a K3 III companion (Khargharia et al., 2010) in a long $6.473\pm 0.001$ day orbit (Casares et al., 1992), and it was the first transient BHXB to have an unambiguous X-ray detection in quiescence ($4\times 10^{33}\,{\rm erg~{}s}^{-1}$ from 0.2-2.4 keV with the ROSAT satellite; Wagner et al. 1994). It also has a well-established distance ($2.39\pm 0.14$ kpc) from radio parallax measurements (Miller-Jones et al., 2009), which has also been measured in the optical (albeit with poorer precision) by Gaia (Gaia Collaboration et al., 2018). Crucially, there already exist $\approx$102 radio observations of V404 Cygni in the Very Large Array (VLA) archive dating back to the 1990s. V404 Cygni is the only quiescent BHXB with such a rich radio dataset (the next best sampled quiescent BHXBs are A0620$-$00 and MWC 656, each of which have 2-3 published radio detections; Gallo et al. 2006; Dzib et al. 2015; Ribó et al. 2017; Dinçer et al. 2018). However, the full VLA archive has yet to be synthesized into a single variability study. The quiescent X-ray variability characteristics of V404 Cygni have already been well-quantified from minute through month timescales (e.g., Wagner et al., 1994; Hynes et al., 2004; Bradley et al., 2007; Bernardini & Cackett, 2014; Rana et al., 2016), making V404 Cygni ripe for radio-to-X-ray comparisons.

In this paper we reanalyse all observations of V404 Cygni in quiescence with the VLA through 2015, and we also consider 14 observations with the Very Long Baseline Array (VLBA). We describe our data reductions in Section 2. In Section 3 we describe the flux density variability characteristics on long (days through decades) and short (minutes to hours) timescales, and we explore radio spectra on long and short timescales. We discuss the radio jet properties in Section 4, and in Section 5 we provide recommendations on how to combat radio variability when coordinating multiwavelength observing campaigns on quiescent BHXBs and on BHXB candidates at luminosities comparable to V404 Cygni.

2 Archival Observations and Data Analysis

V404 Cygni has undergone two major outbursts during the VLA era, the first discovered on 1989 May 22 (Makino, 1989), and the second outburst on 2015 June 15 (Barthelmy et al., 2015; Kuulkers et al., 2015; Negoro et al., 2015; Younes, 2015). After the main 2015 outburst ended in July (see next paragraph), renewed X-ray activity was observed on 2015 December 21 (Malyshev et al., 2015), and V404 Cygni went through a mini-outburst that lasted $\sim$30 days (see, e.g., Kimura et al., 2017; Muñoz-Darias et al., 2017; Tetarenko et al., 2019).

For the 2015 outburst, from Plotkin et al. (2017) we consider that V404 Cygni re-entered quiescence around 2015 July 23, based on when the X-ray spectrum finished softening from $\Gamma\sim 1.6$ to $\Gamma\sim 2.0-2.1$ (see their Figure 1 and Table 3). We therefore exclude VLA observations on or before July 23 in this study, but we include VLA observations after July 23 (four total). Since the X-ray characteristics after July 23 appear similar to those before the outburst, we infer that the physical properties of the underlying accretion flow are no different during these 2015 observations compared to pre-outburst (and by extension, we do not expect physical differences in the radio jet pre- and post-outburst). We therefore include these four observations from 2015 for completeness. However, out of caution, we generally use different symbols/colors to mark the post-outburst observations in figures within this manuscript, and we often remove the post-outburst observations from statistical tests. We do not include any data during or after the 2015 mini-outburst.

The 1989 outburst lasted much longer than the 2015 outburst (e.g., Tetarenko et al., 2019). Similar quality X-ray spectral coverage of the transition into quiescence is not available from 1989, making it difficult to pinpoint when V404 Cygni re-entered quiescence. Han & Hjellming (1992) monitored the 1989 outburst and decay with the VLA for two years, and from their Table 1 the 4.9 and 8.4 GHz flux densities remained brighter than 1 mJy even as late as 1990 September, perhaps indicating elevated accretion (for reference, typical flux densities between outbursts were 0.2–0.4 mJy, e.g., Gallo et al. 2005; Hynes et al. 2009; Rana et al. 2016). We suspect that only the final two epochs in Han & Hjellming (1992) might represent the source in quiescence (taken on 1991 January 31 and 1991 May 31), but we cannot be certain. We therefore conservatively exclude all epochs already reported by Han & Hjellming (1992) from our study, and we begin our dataset with a VLA epoch taken on 1991 September 25. Although, we stress that none of our results are (qualitatively) affected if we were to include the final two epochs from Han & Hjellming (1992).

We found a total of 129 observations over five observing frequencies from 1.4 - 22.5 GHz taken with the historical VLA, i.e., before the VLA was upgraded and re-dedicated as the Karl G. Jansky VLA in 2012. The majority of these data are at 8.4 GHz (86 observations). We also found five observations with the upgraded VLA from 4–8 GHz and four taken from 8–12 GHz. We also include 14 observations with the Very Long Baseline Array (VLBA) at 5.0 GHz (13 observations in 2014) and at 8.4 GHz (1 observation in 2008). A summary of observations is provided in Table 1, and a catalog of flux densities in Table 2. Note that there are two observations at 8.4 GHz for which we could not obtain good calibration solutions. We include entries for those observations in our catalog (Table 2) for completeness, but we have omitted them from the tally of observations in Table 1. In total, we measure flux densities (or limits) for 150 observations.

2.1 Historical VLA

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 version 31DEC14 (aips; 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 $\approx$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 minimise sidelobes from nearby sources in the field. Of particular note is that a bright Jansky-level source lies 16.6 arcmin 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 not 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$\sigma_{\rm rms}$, where $\sigma_{\rm rms}$ is the root-mean-square (rms) noise measured in a blank region of the sky (with upper limits for non-detections calculated as $5\sigma_{\rm rms}$).

2.2 Karl G. Jansky Very Large Array

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 re-entered quiescence (Plotkin et al., 2017).

2.2.1 2013 Quiescence

The 2013 observation lasted $\approx$9 hours in C band (4-8 GHz) in B configuration (maximum baseline $\approx 11$  km), under program code 13B-016 (PI Corbel). Two basebands of bandwidth 1024 MHz were placed at 5.25 and 7.45 GHz. The sources 3C 286 and J2025+337 were used as the primary and secondary calibrators, respectively.

Similar data reduction steps were performed as described in Section 2.1, except for the following: we used the Common Astronomy Software Application v5.1.1 (casa; McMullin et al., 2007), to allow us to account for the larger fractional bandwidth; we Hanning smoothed the data to avoid radio frequency interference from bleeding into nearby frequency channels; we used the primary flux calibrator to solve for delay and complex bandpass solutions; and we imaged the field using the task clean, using Briggs weighting (robust=1.0), and two Taylor terms to model the frequency dependence over the 2 GHz bandwidth. In these observations, the effects of bandwidth smearing on J2025+337 were not significant, on account of the smaller frequency channels (2 MHz vs. 50 MHz), so we were able to adequately deconvolve J2025+337. Flux densities were measured by fitting a two-dimensional Gaussian with the casa task imfit, forcing a point source (fixed to the width of the synthesized beam) at the known location of V404 Cygni.

2.2.2 2015 Post-outburst

Four observations from the end of the 2015 decay are included in this study (from 28 July - 5 August), under program code SG0196 (PI Plotkin). As noted at the beginning of Section 2, V404 Cygni  re-entered quiescence by the time these observations were taken, according to its X-ray spectral behavior (Plotkin et al., 2017). All four observations were taken in the most extended A configuration (maximum baseline $\approx$30 km) using the VLA in subarray mode, where about half of the antennas observed from 4-8 GHz, and the other half observed from 8-12 GHz, yielding strictly simultaneous multi-frequency coverage over four frequencies centered at 5.2, 7.5, 8.6, 11.0 GHz (with 1024 MHz bandwidth at each frequency). These data were reduced using the same procedures as described in Section 2.2.1 (see Plotkin et al. 2017 for details).

2.3 VLBA

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 $\approx$56 min 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, 16.6 arcmin from V404 Cygni) as both the fringe finder and phase reference calibrator, and the somewhat more distant J2023+3153 (Ma et al., 1998, 1.99*∘* from V404 Cygni) as an astrometric check source. We also identfied an additional VLBA observation at 8.4 GHz reported by Miller-Jones et al. (2009) taken on 2008 November 17 under program code BM290, taken in dual circular polarization with 64 MHz bandwidth per polarization (133 min on source), which we re-reduced.222Program BM290 contained four other observations at 8.4 GHz, which were already included in our sample because they used the phased VLA with the VLBA. We used the phased VLA observations in preference to the VLBA, because the phased VLA uses a standard flux calibrator, while the VLBA observations rely on system temperatures for amplitude calibration.

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 unmodelled 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. 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), such that we do not expect to resolve out any jetted radio emission, allowing a fair flux density comparison between VLA and VLBA epochs.

3 Results

3.1 Long-term Variability

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 4 observations). To quantify the distribution of flux densities at each frequency in the presence of non-detections, we perform a survival analysis. Using survfit in R (within the survival package333https://CRAN.R-project.org/package=survival) we calculate the survival function, $S(\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(\log f_{\nu})=1-S(\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$\sigma_{\rm rms}$ level) range between 0.14 – 1.35 mJy bm*-1*.

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 (NADA444https://CRAN.R-project.org/package=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<\log f_{\rm 4.9}/{\rm mJy}\right>=-0.53\pm 0.19$ and $\left<\log f_{\rm 8.4}/{\rm mJy}\right>=-0.53\pm 0.30$, where $f_{\rm 4.9}$ and $f_{\rm 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).555The quoted numbers for the lognormal distributions are not biased by a small number of individual observations in the small- or large flux density tails of the distributions. We show this by bootstrapping the flux density distributions 100 times at each frequency, selecting subsamples of 30 and 50 at 4.9 and 8.4 GHz, respectively, with replacement. After running the same survival analysis on each bootstrapped distribution, we find average values from the 100 distributions to be $\left<\log f_{\rm 4.9}/{\rm mJy}\right>=-0.54\pm 0.19$ and $\left<\log f_{\rm 8.4}/{\rm mJy}\right>=-0.54\pm 0.28,whereerrorsrepresentstandarddeviations$

We find that V404 Cygni shows significant flux density variations that are in excess of statistical fluctuations from measurement errors (which are typically $\pm$0.04 mJy bm*-1* with the historical VLA). Taking 8.4 GHz as example, among 61 detections we find a reduced $\chi^{2}_{r}$ of 46 (for 60 degrees of freedom, as compared to a model with constant flux density). We note that $\chi^{2}_{r}$ is slightly biased by one observation from 2013 with higher signal-to-noise 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_{\rm var}=54\pm 6$ % (Vaughan et al., 2003), which is not biased by the higher signal-to-noise observations (i.e., we also calculate $F_{\rm var}=54\pm 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},\rm var}=\pm f_{\nu}F_{\rm var}$, then propagating errors would yield $\sigma_{\log f_{\nu},{\rm var}}=\pm F_{\rm var}/\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(\tau)$, which characterizes the amount of variability power as a function of time scale in the case of irregularly sampled data,

[TABLE]

where $f(t)$ is a flux density measurement at time $t$, and $\tau$ 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 $\tau_{\rm ij}=t_{j}-t_{i}$. We then bin our set of $V_{\rm ij}$ measures by time difference ($\tau_{\rm ij}$) so that each bin contains 50 data points, and we take the average of those 50 measurements to calculate $V(\tau)$ (where we adopt the midpoint of all time differences within each bin as the value for the time delay $\tau$).666Our error bars represent the error on the mean value of $V(\tau)$ in each time bin. We define our errors in this manner for ease of comparison to Bernardini & Cackett (2014) who provide a structure function for the quiescent X-ray variability of V404 Cygni. The structure function is shown in Figure 3, where we probe long-term variability over timescales of $\sim$10 – 4000 days.

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 summarise below). For example, if variability is characterised as flicker noise (i.e., a power function $P(F)\propto F^{-1}$, where $F$ is the inverse timescale, i.e., frequency, of fluctuations), then the structure function $V(\tau)$ = constant. The constant is expected to be 2$\sigma_{\rm var}^{2}$, with $\sigma_{\rm var}^{2}$ being the variance of the observed flux densities. One could also obtain $V(\tau)=2\sigma_{\rm 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)\propto F^{-2}$) produce first-order structure functions of the form $V(\tau)\propto\tau$, 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 ($\sigma_{\rm err}$) are not meaningful (assuming that $\sigma_{\rm err}$ is dominated by statistical noise). Thus, one expects the structure function to satisfy $V(\tau)\geq 2\sigma_{\rm err}^{2}$ at all timescales.

The structure function in Figure 3 appears flat, with a best-fit slope and normalisation of $\beta=-0.11\pm 0.07$ and $V_{0}=0.16\pm 0.04$ mJy2 (where $V(\tau)=V_{0}\tau^{\beta}$). Measuring the variance directly from the 61 data points yields 2$\sigma_{\rm var}^{2}=0.11$ mJy2. Therefore, the structure function appears consistent with plateauing near 2$\sigma_{\rm var}^{2}$, which signifies either flicker or white noise. In the latter case, we would be probing jet disturbances on timescales ($\gtrsim$10 days) that are longer than characteristic damping timescales.

3.2 Short-term Variability

We also examine variability on sub-hour timescales, focusing on observations long enough (usually $\gtrsim$90 min on source) to produce light curves over multiple time bins (we are able to achieve time resolutions ranging from 3-30 min, 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.

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

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 min timescales on 2015 August 1, and we observe the beginning of a flare in 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 powerlaw to each structure function, and we find powerlaw indices of $\beta=0.29\pm 0.05,0.10\pm 0.11$ and $1.08\pm 0.08$ on 2013 December 2, 2015 August 1, and 2015 August 5, respectively ($V(\tau)\propto\tau^{\beta}$). Note that the structure function on 2015 August 1 is flat, plateauing at $2\sigma_{\rm 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$.

3.2.1 Quiescent Flares

In this subsection we summarise some of the flaring behaviour of V404 Cygni in quiescence. We note that factor of $>$2 changes in flux density appear to be common. However, we do not attempt to define a duty cycle, on account of the limited number of observations (11) with short-term light curves, for which we do not have uniform time resolution.

The 2007 December 2 observation (first reported by Miller-Jones et al., 2008) provides an extreme example of rapid short-term variability, where the 8.4 GHz flux density increased by a factor of 3–4 in $<$10 min, reaching 1.4 mJy, followed by a slower decay that lasted at least 30 min. (The beginning of that light curve also shows a series of three alternating detections and non-detections, suggesting a factor of $>$2.5 flux variability over 10 min). However, not all short-term variability follows a pattern of a fast rise followed by a slower decay. For example, on 1998 May 4 a factor of 2-3 flare rose over $\approx$20 min, on 2009 April 26 a rise and a decay of a factor $\approx$3 in flux density occurred over $\sim$60-90 min in each direction, and on 2015 August 5 V404 Cygni displayed an increase in flux density by a factor $>$2 over $>$60 min.

We also see a potential variety in decay timescales. For example, the 2007 December 2 flare that quickly rose in $<$5-10 min took at least 30 min to decay. At the other extreme, the $\approx$9 hour 2013 December 2 observation appears to be decreasing in flux the entire time (with some smaller-scale variations superposed). If a flare preceded the beginning of that observation, then that implies some flares decay on timescales of at least several hours. From the above we conclude that tens of minutes to hours represent reasonable minimum characteristic timescales for the damping of radio flares. Although flares appear common, V404 Cygni also undergoes quieter periods of time, where either flares are absent or low-amplitude flares occur on short timescales $<$3 min (e.g., 2015 August 1).

3.3 Radio Spectra

For 24 epochs there are observations at multiple frequencies within $\pm$30 min, 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 $\alpha_{r}$ (or place a limit) analytically as $\alpha_{r}=\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 $\nu_{1}$ and $\nu_{2}$, respectively. We assign an error on $\alpha_{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 $\alpha_{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 $\alpha_{r}$ distribution.

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 $\alpha_{r}=-0.26\pm 0.05$, which is consistent with the value $\alpha_{r}=-0.27\pm 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, from $-1.5\lesssim\alpha_{r}\lesssim+1.6$, with the spread of $\alpha_{r}$ being larger at lower flux densities (Figure 7). However, we argue in Section 4.1 that this spread in $\alpha_{r}$ is likely attributed largely to statistical and systematic errors (i.e., large error bars on most $\alpha_{r}$ measurements, especially at lower flux densities, combined with only five epochs having strictly simultaneous multifrequency data).

To quantify the dispersion in radio spectral indices, we fit a Gaussian distribution to the $\alpha_{r}$ measurements (using a Bayesian framework that allows for both upper and lower limits). We find a mean $\left<\alpha_{r}\right>=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 $\alpha_{r}$), we find consistent results: $\left<\alpha_{r}\right>=-0.07\pm 0.25$ and $\sigma_{\alpha_{r}}=0.78\pm 0.21$. Throughout we adopt $\left<\alpha_{r}\right>=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 multi-frequency data allowing us to explore variations in $\alpha_{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 $\approx$15 min (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.777For these light curves, we require $S/N>3$ for the flux densities in each time bin to reduce uncertainties on $\alpha_{r}$. Since the on source integration time should be similar in each time bin, we also remove a small number of time bins ($<$5 across all four sources) where the $\sigma_{\rm rms}$ of the time-resolved images differ by at least a factor of two from the median $\sigma_{\rm rms}$ of all images on each date. Such large variations in $\sigma_{\rm rms}$ might indicate that artifacts remain after the cleaning process, or that there is an unusually large amount of flagged data during that time bin. 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 $\alpha_{r}$ with time; variations of $\alpha_{r}$ are consistent with statistical noise, with reduced $\chi^{2}_{r}$ 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 $\approx$10 min periods on 2013 December 2. While we also see some variations in $\alpha_{r}$ over 10-30 min timescales, they tend to be at the 1–2$\sigma$ 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 $\approx$9 hour observation. Imaging the first 100 min of the observation yields $\alpha_{r}=-0.39\pm 0.10$, while a less steep (and potentially flat) spectral index of $\alpha_{r}=-0.14\pm 0.14$ is measured from images during the last 100 min. The steeper spectral index during the first 100 min appears to be driven by variations at 7.45 GHz that are not mimicked at 5.25 GHz.

4 Discussion

We have presented radio light curves of V404 Cygni in quiescence from 150 VLA and VLBA observations spanning 1991–2015, and we find that factor of 2–4 variations are common on timescales ranging from minutes to decades. Eleven observations are long enough to produce light curves on sub-hour timescales, from which we conclude that radio flares that last from tens of minutes to hours are common. However, there is not a single template to describe quiescent radio variability, either in flare profile, amplitude, or timescale. The observed variety in flare properties could imply that multiple mechanisms control the radio variability, or that a single type of process yields a range of subtly different radiative signatures (e.g., a shock traveling through a steady jet, where the dissipation of energy is highly sensitive to the local conditions at the location of the shock).

Only for a single large flare (2015 August 5) do we have sufficient time resolution to attempt to statistically characterize its properties during the rise. Its structure function displays a slope $\beta\approx 1$, which means that at least some flares are consistent with red noise (i.e., $P(F)\propto F^{-2}$, although we note that we do not have coverage of the entire flare, which could bias our slope measurement). Considering that the structure function of long-term variations is flat (over $\approx$10-4000 day timescales; Figure 3), we suggest that the quiescent jet of V404 Cygni (sometimes) displays flares with ‘random walk’ noise characteristics that dampen to uncorrelated variations on longer timescales. We suspect that the damping time can be as short as tens of minutes to hours, as observed during some flares in Figure 4, but from Figure 3 we can only constrain the damping timescale to $\lesssim$10 days.

We do not assert that all flares have red noise characteristics, since there is only one (large) flare for which we have sufficient time resolution to calculate a structure function, and that observation was taken very shortly after the system returned to quiescence following the 2015 outburst. Nevertheless, the behaviour of red noise flares that dampen to uncorrelated variations on long timescales is reminiscent of decade long radio light curves of BL Lac objects, i.e., low-luminosity active galactic nuclei (AGN) with a jet pointed toward Earth. Hughes et al. (1992) find that BL Lac objects typically show first-order structure functions with slopes $\beta\sim 1$ that become flat at long timescales (i.e., consistent with a damped random walk). They find a broad distribution of characteristic damping timescales from $\approx$1-10 yr for BL Lac objects.888This phenomenology is also common for jet-dominated emission in the optical waveband, where BL Lac object light curves can be well-characterised by a damped random walk (Ruan et al., 2012), and in some cases also in the gamma-ray (where sporadic flaring over a steady flicker/red noise power spectrum is often observed; e.g., Abdo et al. 2010).

To first order, after correcting for beaming effects related to BL Lac orientation, we expect BL Lac objects to be analogous to hard state BHXBs (e.g., Falcke et al., 2004; Körding et al., 2006), and, despite their higher Eddington ratios, we believe that BL Lac objects also make reasonable analogs for comparison to a luminous quiescent BHXB like V404 Cygni (i.e, both types of systems launch compact radio jets from a black hole fed by inefficient accretion). However, to compare BL Lac timescales to V404 Cygni  requires correction for the effects of relativistic beaming from the BL Lac jets, which is difficult given unknown Doppler factors and redshifts for most of the BL Lac objects in Hughes et al. (1992). As an extreme example, we consider a BL Lac object with a Doppler factor $\delta\approx 60$, which would correspond to a very fast jet with bulk Lorentz factor $\Gamma\approx 50$ (e.g., Lister et al., 2009) aligned within only 1 degree to our line of sight. In that case, the intrinsic (i.e., rest-frame) characteristic timescales from Hughes et al. (1992) would be $\lesssim$600 yr (BL Lac objects tend to have low redshifts, e.g., a median $z\sim 0.33$ in Shaw et al. 2013, such that cosmological corrections will be smaller than beaming corrections).

If one were to associate the above timescale with a timescale that scales linearly with black hole mass (e.g., the light travel time across an emitting region that is comparable in size to the radio photosphere of a conical jet), then $\lesssim$600 yr would scale down to $\lesssim$10 min for V404 Cygni (comparing the 9 $M_{\rm Sun}$ mass of V404 Cygni to a typical $3\times 10^{8}M_{\rm Sun}$ mass for a BL Lac object, Plotkin et al. 2011). We suspect that most flares from V404 Cygni decay on timescales of at least several tens of minutes to hours, longer than expectations from the above mass scaling. Therefore, we are either probing physical variations in V404 Cygni that are inaccessible for individual supermassive black holes, or it is possible for some AGN phenomenology to occur on much faster timescales than expected from (simplistic) mass-scaling arguments. We view the latter as a plausible explanation, particularly in light of the existence of changing look AGN, which have accretion disks that appear to operate on much quicker timescales than expected from mass-scaling arguments (e.g., Noda & Done, 2018; Dexter & Begelman, 2018).

Finally, we note that the radio photosphere of the quiescent jet from V404 Cygni (at GHz frequencies) is empirically constrained to be located at a distance $\lesssim$3.4 AU from the black hole (Miller-Jones et al., 2008; Plotkin et al., 2017). The size of the radio emitting region would be even smaller (e.g., by a factor of $\tan\phi$ for a conical jet, where $\phi$ is the jet opening angle; Miller-Jones et al. 2006), such that one should not expect radio variations on timescales longer than $\approx$tens of minutes to be causally connected unless the jet is very slow. Thus, long-term variations likely reflect the jet responding to separate disturbances, like changes in the mass accretion rate through the inner flow, or like changing intrinsic properties of the jet (e.g., how internal energy in the jet is partitioned between particles and magnetic field). That both the radio and X-ray display similar long-term variability characteristics (i.e., flat structure functions, see Bernardini & Cackett 2014) might support that both wavebands release radiative power by tapping into the same energy reservoir, as suggested by, e.g., Malzac et al. (2004).999We note that we are probing longer timescales than considered by Malzac et al. 2004, who consider X-ray variations associated with dynamical timescales at 10-100 Schwarzschild ($\approx 0.1$ s), while in our case the X-ray and radio variations are over timescales longer than several days. However, the general idea that both the radio and X-ray emission regions are responding to a common physical driver seems a reasonable interpretation.

4.1 Spectral Variations

As noted in Section 3.3.1, we measure radio spectral indices ranging from $-1.5<\alpha_{r}<+1.6$, with a mean $\alpha_{r}=0.02$ and a standard deviation $\sigma_{\rm\alpha_{r}}=0.65$. We suspect that the large range in $\alpha_{r}$ can be attributed to the poorer sensitivity of the historical VLA, and multifrequency data that can be offset by $\pm$30 min. For example, a factor of two variability within 30 min would adjust $\alpha_{r}$ by $\pm 1.3$. We are therefore cautious not to over interpret the apparently large range of measured $\alpha_{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 ($\alpha_{r}=-0.26\pm 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 $\alpha_{r}=-0.26\pm 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 optically 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 at approximately 50 and 150 min 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 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) or 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.

5 Coordinating Multiwavelength Observations

Part of our motivation for this work is to quantify the level to which variability-related systematic uncertainties can influence the location of quiescent BHXBs in the radio/X-ray luminosity plane ($L_{\rm R}-L_{\rm X}$). Understanding these uncertainties is important for studies on disk/jet couplings (e.g., fitting $L_{\rm R}-L_{\rm X}$ correlation slopes for individual objects), and also for studies that appeal to $L_{\rm R}-L_{\rm X}$ to identify new BHXB candidates. Below we compare the radio variability of V404 Cygni to its X-ray variability (taken from the literature), and we recommend some guidelines for reducing variability-induced systematics when coordinating and interpreting radio/X-ray observations of quiescent BHXB (candidates) with luminosities comparable to V404 Cygni.

5.1 X-ray Variability in Quiescence

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_{\rm var,X-ray}=57\pm 3$% ($\sim$0.25 dex), and they obtain a flat structure function over timescales of $\approx$5-80 days. These X-ray results are similar to our radio results, where we find $F_{\rm var,radio}=54\pm 6$% ($\sim$0.23 dex) and a flat structure function (although our study extends to longer timescales of of $\approx$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 $>$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 not been 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 hours (2003 July 29 and 2013 December 2), and neither attempt has shown obvious radio/X-ray correlations with light curves on 30–100 min time bins (Hynes et al., 2009; Rana et al., 2016). However, detecting correlated variations may require finer time resolution.

5.2 Recommendations for Coordinating Radio and X-ray Observations

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:

• •

radio and X-ray observations should be scheduled as simultaneously as possible. If the source is bright enough to provide sufficient signal-to-noise on $<$10 min 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 hours).

• •

If data are non-simultaneous, or if one does not have sufficient time resolution to resolve individual flares, then inflate the error bars by $\approx$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.

• •

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 $\alpha_{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).

5.3 A Comparison to the Transitional Millisecond Pulsar PSR J1023+0038

While quiescent BHXBs tend to have higher $L_{\rm R}/L_{\rm 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_{\rm X}\approx 5\times 10^{32}~{}{\rm erg~{}s}^{-1}$) and a high mode ($L_{\rm X}\approx 3\times 10^{33}~{}{\rm erg~{}s}^{-1}$; e.g., Patruno et al. 2014; Stappers et al. 2014; Jaodand et al. 2016; Bogdanov et al. 2015). 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 anti-correlated 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_{\rm R}\approx 2\times 10^{27}\,{\rm erg~{}s}^{-1}$.

In Figure 9 we highlight the location of J1023 in the $L_{\rm R}-L_{\rm X}$ plane101010Data taken from Bahramian et al. (2018), see https://github.com/bersavosh/XRB-LrLx_pub. (large green triangles) compared to V404 Cygni in quiescence and in 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 dataset ($\approx 5\times 10^{27}-5\times 10^{28}~{}{\rm erg~{}s}^{-1}$, which corresponds to flux densities ranging from 0.14-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\times 10^{21}$ cm*-2* (Bernardini & Cackett, 2014), using X-ray count rates extracted from the online Swift-XRT product generator tool111111http://www.swift.ac.uk/user_objects (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_{\rm R}-L_{\rm X}$ plane if considering non-simultaneous data. The mean radio flux density from our study ($\log f_{\nu}/{\rm mJy}=-0.53$) corresponds to $\log L_{\rm R}/{\rm erg~{}s}^{-1}\,\approx 28$, such that V404 Cygni likely spends most of its time toward the bottom left of the shaded region (with the upper right region representing periods of flaring activity).

After accounting for the range of flux density variations exhibited by V404 Cygni, we expect its radio luminosity to always be $\gtrsim$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 $\lesssim$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 hours, 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.121212Although in some cases, e.g., a system being very radio bright or radio faint, one could reasonably disfavor or favor a neutron star from radio and X-ray luminosities. Other multiwavelength data should naturally also be consulted to, e.g., search for an orbital period, and to search for 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:

• •

frequent and relatively short radio flares that both rise and decay on rapid (minute) timescales may favor a neutron star. Bogdanov et al. (2018) find radio flares that last $\lesssim$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 observe as many as 3-4 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;

• •

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) show light curves of J1023 from 13 different epochs (see their Figure 5), and they find a case where the radio flux density steadily increases by an order of magnitude over $\approx$30 min. While V404 Cygni has also been observed to show flares with comparably slow rise times, none of those $>$30 min 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\pm 0.14$ kpc, Miller-Jones et al. 2009, and J1023 is at $1.368^{+0.042}_{-0.039}$ 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 from V404 Cygni unlikely represents the full range of variability that can be displayed by either class of systems (e.g., as evidenced by the nearly 10 hour 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, non-flaring time periods will dilute the radio flux density, thereby moving the tMSP farther away from the quiescent BHXB space in $L_{\rm R}-L_{\rm X}$. Finally, in some cases the X-ray spectrum can provide additional diagnostic power: J1023 has an X-ray photon index of $\Gamma\approx 1.7$ in both low and high X-ray modes (Bogdanov et al., 2018), while quiescent BHXBs are well established to have a softer $\Gamma\approx 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 report a range extending from $-0.5\lesssim\alpha_{r}\lesssim 0.4)$, and that $\alpha_{r}$ is usually not well-measured in low-luminosity sources.

6 Summary

We present archival radio observations of V404 Cygni spanning 24 years (1991–2015), providing the most stringent long-term constraints to date on the radio variability of a synchrotron jet from a quiescent BHXB. We find flux densities that follow a lognormal distribution, with mean and standard deviation (in $\log f_{\nu}/{\rm mJy}$) of $-0.53\pm 0.19$ and $-0.53\pm 0.30$ at 4.9 and 8.4 GHz, respectively. Factor of $>$2-4 variations are common on every observable timescale from minutes to decades. As expected, the radio spectrum of V404 Cygni is flat on average ($\left<\alpha_{r}\right>=0.02\pm 0.65$, where the error represents the standard deviation), but we also find that $\alpha_{r}$ becomes significantly negative on at least one epoch.

Over two decades of observations, V404 Cygni displays a flat structure function, such that the long-term flux density variations (days – years) are consistent with either flicker noise ($P(F)\propto F^{-1}$) or white noise. On epochs when we have sufficient quality data, we observe individual flares that appear to decay within minutes to hours. These results are consistent with an interpretation of shot noise probed over a timescale that is longer than the characteristic damping time of each disturbance. We suspect that typical flare decay timescales are on the order of tens of minutes to hours, but we are formally only capable of constraining that timescale to $<$10 days (with tens of minutes a reasonable lower bound). These properties may be expected from shock instabilities traveling through a steady, compact jet. Given similar characteristics observed in the X-ray band (Bernardini & Cackett, 2014), the radio variability appears consistent with scenarios where the X-ray and radio radiative processes are powered by a common energy source (e.g., Malzac et al., 2004).

Finally, we provide recommendations for combatting variability-induced systematics when attempting to place accreting compact objects onto the radio/X-ray luminosity plane, as is commonly done in surveys for quiescent BHXBs. We recommend that radio and X-ray observations be taken as simultaneously as possible to allow direct detections of flares during each observation. If the data are not of sufficient quality to detect individual flares, or the data are not simultaneous, then we recommend inflating error bars on radio and X-ray luminosities by 0.25 dex. If one must extrapolate radio observations to different frequencies, then ideally one would be able to measure a radio spectrum from strictly simultaneous multi-frequency observations. Otherwise, a flat radio spectrum is a reasonable approximation, except that error bars should also be adjusted according to expectations that the radio spectrum can vary ($\sigma_{\rm\alpha_{r}}=\pm 0.6$ is a reasonable uncertainty). Finally, we repeat warnings by Bogdanov et al. (2018) that some accreting neutron stars (i.e., tMSPs) can obtain radio and X-ray luminosities comparable to those achieved by quiescent BHXBs, such that other types of multiwavelength data (e.g., optical, ultraviolet, and X-ray spectra and timing) should be considered when attempting to identify the nature of accreting compact objects when mass functions are not available.

We are grateful to Robert Hjellming, who obtained many of the observations presented in this paper and who lay the foundations for understanding stellar sources at radio wavelengths. We are also grateful to Michael Rupen for obtaining a large portion of the VLA observations used in this paper. We thank the anonymous referee for helpful comments. The National Radio Astronomy Observatory is a facility of the National Science Foundation (NSF) operated under cooperative agreement by Associated Universities, Inc. This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester. We acknowledge support from NSF grant AST-1308124. R.M.P. acknowledges support from Curtin University through the Peter Curran Memorial Fellowship. J.C.A.M.J. is supported by an Australian Research Council Future Fellowship (FT140101082). L.C. acknowledges support from NSF AST-1412549. J.S. acknowledges support from the Packard Foundation.