Extreme Jet Ejections from the Black Hole X-ray Binary V404 Cygni
A.J. Tetarenko, G.R. Sivakoff, J.C.A Miller-Jones, E.W. Rosolowsky, G., Petitpas, M. Gurwell, J. Wouterloot, R. Fender, S. Heinz, D. Maitra, S.B., Markoff, S. Migliari, M. P. Rupen, A.P. Rushton, D.M. Russell, T.D. Russell,, C.L. Sarazin
TL;DR
This study presents a detailed multi-frequency analysis of extreme jet ejections from the black hole X-ray binary V404 Cygni during its 2015 outburst, revealing rapid flaring activity and modeling jet dynamics.
Contribution
It introduces a novel multi-frequency observational approach combined with a customized relativistic jet model to analyze bi-polar ejections in V404 Cygni.
Findings
Detected multiple rapid flares reaching Jy levels across all frequencies.
Modeled 8 bi-polar jet ejections explaining observed emission.
Highlighted the critical role of mm/sub-mm bands in jet studies.
Abstract
We present simultaneous radio through sub-mm observations of the black hole X-ray binary (BHXB) V404 Cygni during the most active phase of its June 2015 outburst. Our hour long set of overlapping observations with the Very Large Array, the Sub-millimeter Array, and the James Clerk Maxwell Telescope (SCUBA-2), covers 8 different frequency bands (including the first detection of a BHXB jet at ), providing an unprecedented multi-frequency view of the extraordinary flaring activity seen during this period of the outburst. In particular, we detect multiple rapidly evolving flares, which reach Jy-level fluxes across all of our frequency bands. With this rich data set we performed detailed MCMC modeling of the repeated flaring events. Our custom model adapts the van der Laan synchrotron bubble model to include twin bi-polar ejections, propagating away from the…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4
Figure 5
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11| Parameterb | Description | Prior Distributionc | Prior | Prior |
| Minimum | Maximum | |||
| amplitude of compact jet component at in mJy | uniform | 10 | 1000 | |
| spectral index | truncated normal (, ) | -1 | 0 | |
| ∗ | ejection time of bi-polar components in secondsd,e | uniform | ||
| ∗ | inclination angle of the jet axis in degrees | truncated normal (, ) | 0 | 90 |
| ∗ | observed opening angle in degrees | uniform | 1 | 20 |
| ∗ | optical depth at the reference frequency | truncated normal (, ) | 1 | 3 |
| ∗ | peak flux density at the reference frequency in mJy | truncated normal (, ) | 0 | 6000 |
| ∗ | bulk jet speed in units of c | inverse gamma ()f | 0 | 1 |
| Compact Jet Parametersa | |||||||||
| (mJy) | |||||||||
| Individual Jet Ejecta Parameters | |||||||||
| Ejection | (HH:MM:SS.S) | (MJD) | (degrees) | (degrees) | b | (mJy) | (v/c) | (v/c)c | |
| 1 | |||||||||
| 2 | |||||||||
| 3 | |||||||||
| 4 | |||||||||
| 5d | |||||||||
| 6 | |||||||||
| 7 | |||||||||
| 8 | |||||||||
| Compact Jet Parameters | |||||||
| (mJy) | |||||||
| Individual Jet Ejecta Parameters | |||||||
| Ejection | (sec) | (MJD) | (degrees) | (degrees) | (mJy) | (v/c) | |
| 1 | |||||||
| 2 | |||||||
| 3 | |||||||
| 4 | |||||||
| 5 | |||||||
| 6 | |||||||
| 7 | |||||||
| 8 | |||||||
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
\settimeformat
hhmmsstime
Extreme Jet Ejections from the Black Hole X-ray Binary V404 Cygni
A.J. Tetarenko,1 G.R. Sivakoff,1 J.C.A. Miller-Jones,2 E.W. Rosolowsky,1 G. Petitpas,3 M. Gurwell,3 J. Wouterloot,4 R. Fender,5 S. Heinz,6 D. Maitra,7 S.B. Markoff,8 S. Migliari,9,10 M.P. Rupen,11,12 A.P. Rushton,5,13 D.M. Russell,14 T.D. Russell,8,2 and C.L. Sarazin15
1Department of Physics, University of Alberta, CCIS 4-181, Edmonton, AB T6G 2E1, Canada
2International Centre for Radio Astronomy Research- Curtin University, GPO Box U1987, Perth, WA 6845, Australia
3Harvard-Smithsonian Center for Astrophysics, Cambridge, MA 02138
4East Asian Observatory, 660 North Aohoku Place, University Park, Hilo, Hawaii 96720, USA
5Department of Physics, Astrophysics, University of Oxford, Keble Road, Oxford OX1 3RH, UK
6Astronomy Department, University of Wisconsin-Madison, 475. N. Charter St., Madison, WI 53706, USA
7Department of Physics and Astronomy, Wheaton College, Norton, MA 02766, USA
8Anton Pannekoek Institute for Astronomy, University of Amsterdam, P.O. Box 94249, NL-1090 GE Amsterdam, the Netherlands
9Department of Astronomy and Meteorology, University of Barcelona, Martí i Franquès 1, 08028 Barcelona, Spain
10XMM-Newton Science Operations Centre, ESAC/ESA, PO Box 78, 28691 Villanueva de la Cañada, Madrid, Spain
11National Research Council, Herzberg Astronomy and Astrophysics, 717 White Lake Road, PO Box 248, Penticton, BC V2A 6J9, Canada
12National Radio Astronomy Observatory, P.O. Box 0, Socorro, NM 87801, USA
13School of Physics and Astronomy, University of Southampton, Highfield, Southampton SO17 1BJ, UK
14New York University Abu Dhabi, P.O. Box 129188, Abu Dhabi, United Arab Emirates
15Department of Astronomy, University of Virginia, P.O. Box 400325, Charlottesville, VA 22904, USA E-mail: [email protected]
(Accepted XXX. Received YYY; in original form ZZZ)
Abstract
We present simultaneous radio through sub-mm observations of the black hole X-ray binary (BHXB) V404 Cygni during the most active phase of its June 2015 outburst. Our hour long set of overlapping observations with the Very Large Array, the Sub-millimeter Array, and the James Clerk Maxwell Telescope (SCUBA-2), covers 8 different frequency bands (including the first detection of a BHXB jet at ), providing an unprecedented multi-frequency view of the extraordinary flaring activity seen during this period of the outburst. In particular, we detect multiple rapidly evolving flares, which reach Jy-level fluxes across all of our frequency bands. With this rich data set we performed detailed MCMC modeling of the repeated flaring events. Our custom model adapts the van der Laan synchrotron bubble model to include twin bi-polar ejections, propagating away from the black hole at bulk relativistic velocities, along a jet axis that is inclined to the line of sight. The emission predicted by our model accounts for projection effects, relativistic beaming, and the geometric time delay between the approaching and receding ejecta in each ejection event. We find that a total of 8 bi-polar, discrete jet ejection events can reproduce the emission that we observe in all of our frequency bands remarkably well. With our best fit model, we provide detailed probes of jet speed, structure, energetics, and geometry. Our analysis demonstrates the paramount importance of the mm/sub-mm bands, which offer a unique, more detailed view of the jet than can be provided by radio frequencies alone.
keywords:
black hole physics — ISM: jets and outflows — radio continuum: stars — stars: individual (V404 Cygni, GS 2023+338) — submillimetre: stars — X-rays: binaries
††pubyear: 2017††pagerange: Extreme Jet Ejections from the Black Hole X-ray Binary V404 Cygni–E
1 Introduction
Black hole X-ray binaries (BHXBs), the rapidly evolving, stellar-mass counterparts of active galactic nuclei, are ideal candidates with which to study accretion and accretion-fed outflows, such as relativistic jets. These transient binary systems, containing a black hole accreting mass from a companion star, occasionally enter into bright outburst phases lasting days to weeks, providing a real time view of the evolving relativistic jets (probed by radio through IR frequencies) and accretion flow (probed at X-ray frequencies).
BHXBs display two different types of relativistic jets, dependent on the mass accretion rate in the system (Fender et al., 2004). At lower mass accretion rates ()111The Eddington luminosity is the theoretical limit where, assuming ionized hydrogen in a spherical geometry, radiation pressure balances gravity. This limit corresponds to , where is the black hole mass., during the hard accretion state (see Remillard & McClintock 2006 and Belloni 2010 for a review of accretion states in BHXBs), a steady, compact synchrotron-emitting jet is believed to be present in all BHXBs. It has also been shown that the compact jet is not only present during outburst phases, but can persist down into quiescence, at (Gallo et al., 2005; Plotkin et al., 2013; Plotkin et al., 2015; Plotkin et al., 2016). At higher mass accretion rates, during the transition between accretion states, discrete jet ejecta are launched (e.g., Mirabel & Rodríguez 1994; Hjellming & Rupen 1995; Corbel et al. 2002; Miller-Jones et al. 2012a), and the compact jet may become quenched (Fender et al., 1999d; Corbel et al., 2001; Russell et al., 2011; Coriat et al., 2011; Rushton et al., 2016). A small number of BHXBs have been observed to display multiple jet ejection events within a single outburst (e.g., Mirabel & Rodríguez 1994; Hjellming & Rupen 1995; Tingay et al. 1995; Fender et al. 1999a; Kuulkers et al. 1999; Brocksopp et al. 2002; Brocksopp et al. 2013).
Compact jets are characterized by a flat to slightly inverted optically thick spectrum (; where ; Fender 2001), extending from radio up to sub-mm or even infrared frequencies (Corbel & Fender, 2002; Casella et al., 2010; Tetarenko et al., 2015a). Around infrared frequencies the jet emission becomes optically thin (; Russell et al. 2013a), resulting in a spectral break. Each frequency below this break probes emission (from the optical depth, surface) coming from a narrow range of distances downstream in the jet, where higher frequencies originate from regions along the jet axis that are closer to where the jet is launched (Blandford & Königl, 1979; Falcke & Biermann, 1995). The exact spectral shape (i.e., spectral index, location of the spectral break) is believed to evolve with changing jet properties such as geometry, magnetic field structure, and particle density profiles (Heinz & Sunyaev, 2003; Markoff et al., 2005; Casella & Pe’er, 2009; Russell et al., 2013b; van der Horst et al., 2013; Russell et al., 2014), as well as the plasma conditions in the region where the jet is first accelerated (Polko et al., 2010, 2013, 2014; Koljonen et al., 2015).
In contrast to the compact jets, jet ejecta are characterized by an optically thin spectrum (), give rise to bright flaring activity, and can be routinely resolved with Very Long Baseline Interferometry (VLBI; e.g., Fender 2006). The accompanying flares typically have well defined rise and decay phases, where the flares are usually optically thick in the rise phase, until the self-absorption turnover in the spectrum has passed through the observing band. These jet ejection events are believed to be the result of the injection of energy and particles to create an adiabatically expanding synchrotron emitting plasma, threaded by a magnetic field (i.e., van der Laan synchrotron bubble model, hereafter referred to as the vdL model; van der Laan 1966; Hjellming & Johnson 1988; Hjellming & Han 1995). In this model, as the source expands the evolving optical depth results in the distinct observational signature of the lower frequency emission being a smoothed, delayed version of the higher frequency emission. The ejection events have been linked to both X-ray spectral and timing signatures (e.g., Fender et al. 2009; Miller-Jones et al. 2012a; Russell et al. 2014; Kalemci et al. 2016), although a definitive mechanism or sequence of events leading to jet ejection has not yet been identified.
Additionally, an extremely rare jet phenomenon, so called jet oscillation events, has also been observed in two BHXBs, GRS 1915105 (radio, mm, IR; Pooley & Fender 1997) and V4641 Sgr (optical band; Uemura et al. 2004). Such rare events seem to occur only when the accretion rate is at very high fractions of the Eddington rate. These quasi-periodic oscillations (see Fender & Belloni 2004 for a review) show lower frequency emission peaking at later times (consistent with the vdL model for expanding discrete jet ejecta), rise and decay times of the repeated flares that are similar at all frequencies, and time lags between frequencies that vary within a factor of two. Moreover, no discrete moving components were resolved with VLBI during these oscillation events (although we note this could very well be due to sensitivity limits or the difficulty of synthesis imaging of fast-moving, time-variable components). As such, the exact nature of these events remains unclear, with theories including discrete plasma ejections, internal shocks in a steady flow, or variations in the jet power in a self-absorbed, conical outflow (e.g., Fender & Pooley 1998, 2000; Collins et al. 2003). In GRS 1915105 , these oscillations have also been clearly associated with dips in hard X-ray emission, possibly linking the launching of jet ejecta to the ejection and refilling of the inner accretion disc or coronal flow (Mirabel et al., 1998; Belloni et al., 1997; Vadawale et al., 2001).
While several transient BHXBs may undergo an outburst period in a given year, in which the jet emission becomes bright enough for detailed multi-wavelength studies, only rare (e.g., once per decade) outbursts probe the process of accretion and the physics of accretion-fed outflows near (or above) the Eddington limit. Observing the brightest and most extreme phases of accretion during these outbursts presents us with a unique opportunity to study jet and accretion physics in unprecedented detail. On 2015 June 15, the BHXB V404 Cygni entered into one of these rare near-Eddington outbursts. In this paper we report on our simultaneous radio through sub-mm observations of V404 Cygni during the most active phase of this outburst.
1.1 V404 Cygni
V404 Cygni (aka GS 2023338; hereafter referred to as V404 Cyg) is a well studied BHXB that has been in a low-luminosity quiescent state since its discovery with the Ginga satellite in 1989 (Makino, 1989). This source has been observed to undergo a total of three outbursts prior to 2015; most recently in 1989 (Han & Hjellming, 1992; Terada et al., 1994; Oosterbroek et al., 1997), and two prior to 1989 which were recorded on photographic plates (Richter, 1989). V404 Cyg is known to display bright X-ray luminosities and high levels of multi-wavelength variability, both in outburst and quiescence (Tanaka & Lewin, 1995; Hynes et al., 2002; Hjellming & Han, 1989; Kitamoto et al., 1989). The prolonged quiescent period of V404 Cyg, and high quiescent luminosity (; Corbel et al. 2008), has allowed the complete characterization of the system. The optical extinction is low, with , enabling the study of the optical counterpart, and the determination of the mass function as (Casares et al., 1992; Casares & Charles, 1994). Subsequent modelling determined the black hole mass to be , with an inclination angle of , and an orbital period of 6.5 days (Khargharia et al., 2010; Shahbaz et al., 1994). However, we note that this inclination angle estimate is dependent on the assumed level of accretion disc contamination in the optical light curves being modelled. Khargharia et al. (2010) assumed accretion disc contamination, but given that V404 Cyg is known to be variable in quiescence in the optical, it is plausible that the accretion disc contamination may be larger (Zurita et al., 2003; Bernardini et al., 2016a), which would imply a larger inclination angle. Further, the faint, unresolved radio emission from the quiescent jets was used to determine a model-independent parallax distance of kpc (Miller-Jones et al., 2009b), making V404 Cyg one of the closest known BHXBs in the Galaxy. The close proximity, well-determined system parameters, and bright multi-wavelength activity make this system an ideal target for jet and accretion studies.
On 2015 June 15222Bernardini et al. (2016a) serendipitously detected an optical precursor to this outburst on June 8/9, approximately one week prior to the first X-ray detection., V404 Cyg entered into its fourth recorded outburst period. The source began exhibiting bright multi-wavelength flaring activity (e.g., Ferrigno et al. 2015a; Gandhi et al. 2015; Gazeas et al. 2015; Mooley et al. 2015; Motta et al. 2015a, b; Tetarenko et al. 2015b, c) immediately following the initial detection of the outburst in X-rays (Barthelmy et al., 2015; Negoro et al., 2015; Kuulkers et al., 2015), and swiftly became the brightest BHXB outburst seen in the past decade. This flaring behaviour was strikingly similar to that seen in the previous 1989 outburst (Terada et al., 1994; Oosterbroek et al., 1997; Zycki et al., 1999). Towards the end of June the flaring activity began to diminish across all wavelengths (e.g., Ferrigno et al. 2015b; Martin-Carrillo et al. 2015; Oates et al. 2015; Scarpaci et al. 2015; Tetarenko et al. 2015d; Tsubono et al. 2015), and the source began to decay (Sivakoff et al., 2015a, b), reaching X-ray quiescence333V404 Cyg entered optical quiescence in mid October 2015 (Bernardini et al., 2016b). in early to mid August (Sivakoff et al., 2015c; Plotkin et al., 2017). V404 Cyg also showed a short period of renewed activity from late December 2015 to early January 2016 (e.g., Lipunov et al. 2015; Trushkin et al. 2015; Beardmore et al. 2015; Malyshev et al. 2015; Tetarenko et al. 2016; Motta et al. 2016a), and Munoz-Darias et al. (2016b) present radio, optical, and X-ray monitoring during this period.
We organized simultaneous observations with the Karl G. Jansky Very Large Array (VLA), the Sub-millimeter Array (SMA), and the James Clerk Maxwell Telescope (JCMT) on 2015 June 22 (approximately one week following the initial detection of the outburst), during which time some of the brightest flaring activity seen in the entire outburst was observed. This comprehensive data set gives us an unprecedented multi-frequency view of V404 Cyg, in turn allowing us to perform detailed multi-frequency light curve modelling of the flaring events. In §2 we describe the data collection and data reduction processes. §3 describes the custom procedures our team developed to extract high time resolution measurements from our data. In §4 we present our multi-frequency light curves, outline our model, and describe the modelling process. A discussion of our best fit model is presented in §5, and a summary of our work is presented in §6.
2 Observations and Data Analysis
2.1 VLA Radio Observations
We observed V404 Cyg with the VLA (Project Code: 15A-504) on 2015 June 22, with scans on source from 10:37:24–14:38:39 UTC () in both C () and K () band. The array was in its most extended A configuration, where we split the array into 2 sub-arrays of 14 (sub-array A) and 13 (sub-array B) antennas. Sub-array A observed the sequence C-K-C, while sub-array B observed the sequence K-C-K, with an 80 second on target and 40 second on calibrator cycle, in order to obtain truly simultaneous observations across both bands. All observations were made with an 8-bit sampler, comprised of 2 base-bands, with 8 spectral windows of 64 2 MHz channels each, giving a total bandwidth of 1.024 GHz per base-band. Flagging, calibration, and imaging of the data were carried out within the Common Astronomy Software Application package (CASA; McMullin et al. 2007) using standard procedures. We used 3C48 (0137331) as a flux calibrator, and J20253343 as a phase calibrator for both sub-arrays. No self-calibration was performed. Due to the rapidly changing flux density of the source, we imaged the source (with natural weighting; see the Appendix for details on our choice of weighting scheme) on timescales as short as the correlator dump time (2 seconds) using our custom CASA timing scripts (see §3.1 for details).
2.2 SMA (Sub)-Millimetre Observations
We observed V404 Cyg with the SMA (Project Code: 2015A-S026) on 2015 June 22, with scans on source from 10:16:17–18:20:47 UTC (), and the correlator tuned to an LO frequency of 224 GHz. The array was in the sub-compact configuration with a total of 7 antennas (out of a possible 8 antennas). These observations were made with both the ASIC and SWARM (Primiani et al., 2016) correlators active, to yield 2 side-bands, with 48 spectral windows of 128 0.8125 MHz channels (ASIC) and an additional GHz spectral windows (SWARM), giving a total bandwidth of 8.32 GHz per side-band. The SWARM correlator had a fixed resolution of 101.6 kHz per channel, and thus originally 16383 channels for each SWARM spectral window. Given the continuum nature of these observations, we performed spectral averaging, to yield 128 13 MHz channels in both SWARM spectral windows, to match the number of channels in the ASIC spectral windows, and in turn make it easier to combine ASIC and SWARM data. We used 3C454.3 (J22531608) as a bandpass calibrator, MWC349a and J20153710 as phase calibrators, and Neptune and Titan as flux calibrators444The SMA calibrator list can be found at http://sma1.sma.hawaii.edu/callist/callist.html.. We note that only the second IF (spectral windows 25-50) was used for flux calibration in the upper side-band due to a CO line that was present in both flux calibrators at 230.55 GHz. Our observing sequence consisted of a cycle of 15 min on target and 2.5 min on each of the two phase calibrators. As CASA is unable to handle SMA data in its original format, prior to any data reduction we used the SMA scripts, sma2casa.py and smaImportFix.py, to convert the data into CASA MS format, perform the corrections, and spectrally average the two SWARM spectral windows. Flagging, calibration, and imaging of the data were then performed in CASA using procedures outlined in the CASA Guides for SMA data reduction555Links to the SMA CASA Guides and these scripts are publicly available at https://www.cfa.harvard.edu/sma/casa.. Due to the rapidly changing flux density of the source, we imaged the source (with natural weighting; see the Appendix for details on our choice of weighting scheme) on timescales as short as the correlator dump time (30 seconds) using our custom CASA timing scripts (see §3.1 for details).
2.3 JCMT SCUBA-2 (Sub)-Millimetre Observations
We observed V404 Cyg with the JCMT (Project Code: M15AI54) on 2015 June 22 from 10:49:33–15:12:40 UTC (), in the (350 GHz) and (666 GHz) bands. The observation consisted of eight min scans on target with the SCUBA-2 detector (Chapin et al., 2013; Holland et al., 2013). To perform absolute flux calibration, observations of the calibrator CRL2688 were used to derive a flux conversion factor (Dempsey et al., 2012). The daisy configuration was used to produce 3 arcmin maps of the target source region. During the observations we were in the Grade 3 weather band with a 225 GHz opacity of 0.095–0.11. Data were reduced in the StarLink package using both standard procedures outlined in the SCUBA-2 cookbook666http://starlink.eao.hawaii.edu/devdocs/sc21.htx/sc21.html and SCUBA-2 Quickguide777https://www.eaobservatory.org/jcmt/instrumentation/continuum
/scuba-2/data-reduction/reducing-scuba2-data, as well as a custom procedure to create short timescale maps (timescales shorter than the 30 minute scan timescale) to extract high time resolution flux density measurements of the rapidly evolving source (see §3.2 for details).
3 High Time Resolution Measurements
3.1 VLA and SMA
To obtain high time resolution flux density measurements of V404 Cyg from our interferometric data sets (VLA and SMA) we developed a series of custom scripts that run within CASA. A detailed account of the development and use of these scripts will be presented in Tetarenko & Koch et al. 2017, in prep., although we provide a brief overview of the capabilities here.
Our scripts split an input calibrated CASA Measurement Set into specified time intervals for analysis in the image plane or the uv plane. In the image plane analysis, each time interval is cleaned and the flux density of the target source is measured by fitting a point source in the image plane with the native CASA task imfit. All imaging parameters (e.g., image size, pixel size, number of CLEAN iterations, CLEAN threshold) can be fully specified. In the uv plane analysis, the uvmultfit package (Marti-Vidal et al., 2014) is used to measure flux density of the target source. In either case, an output data file and plot of the resulting light curve are produced. These scripts are publicly available on github888https://github.com/Astroua/AstroCompute_Scripts, and are being implemented as a part of an interactive service our team is developing to run on Amazon Web Services Cloud Resources.
All VLA and SMA flux density measurements output from this procedure (fitting only in the image plane) are provided in a machine readable table online, which accompanies this paper. Additionally, to check that the variability we observed in V404 Cyg is dominated by intrinsic variations in the source and not due to atmospheric or instrumental effects, we also ran our calibrator sources through these scripts (see the Appendix for details).
3.2 JCMT SCUBA-2
To obtain high time resolution flux density measurements of V404 Cyg from our JCMT SCUBA-2 data we developed a custom procedure to produce a data cube, containing multiple maps of the target source region, at different time intervals throughout our observation.
We run the StarLink Dynamic Iterative Mapmaker tool on each of the target scans, using the bright compact recipe, with the addition of the shortmap parameter. The shortmap parameter allows the Mapmaker to create a series of maps, each of which will include data from a group of adjacent time slices. The number of time slices included in each map is equivalent to the shortmap parameter value. At we use to produce 362 time slices for a 32 minute scan, resulting in 5 second time bins. At would produce the same number of time slices, where a factor of 2 is applied as the default pixel size is 2 arcsec at and 4 arcsec at . However, as the noise is higher at , we use to produce 32 time slices for a 32 minute scan, resulting in 60 second time bins. The stackframes task is then used to combine all of the short maps into a cube for each scan. The and parameters ensure the maps are ordered chronologically in time, with the resulting cube having the dimensions, position X (pixels), position Y (pixels), time (MJD). Using the wcsmosaic task we then combined the cubes from all the scans. We calibrated the combined cube into units of Jy using the scuba2checkcal and cmult tasks. Finally, the combined cube can be viewed in Gaia, and converted to FITS format with the ndf2fits task.
To extract flux densities from each time slice in the combined cube, we fit a 2D gaussian999The python package gaussfitter is used in the gaussian fitting; https://github.com/keflavich/gaussfitter with the size of the beam (FWHM of 15.35 arcsec at and 10.21 arcsec at ; derived using the task scuba2checkcal) to each slice of the cube. All JCMT SCUBA-2 flux density measurements output from this procedure are provided in a machine readable table online, which accompanies this paper. As with our interferometric data sets, to check that the variability we observed in V404 Cyg is dominated by intrinsic variations in the source and not due to atmospheric or instrumental effects, we also ran this procedure on our calibrator source scans (see the Appendix for details).
4 Results
4.1 Multi-frequency Light Curves
A composite light curve of all of our VLA, SMA and JCMT observations from June 22 is presented in Figure 1. We observe rapid multi-frequency variability in the form of multiple large scale flares, reaching Jy flux levels. In the SMA data, the largest flare (at 13:15 UTC) rose from mJy to a peak of Jy on a timescale of min. The JCMT SCUBA-2 data appear to track the SMA data closely, with the largest flare at GHz rising from mJy to a peak of Jy on a timescale of min. This is the largest mm/sub-mm flare ever observed from a BHXB, far surpassing even the brightest events in GRS 1915105 (Fender & Pooley, 2000). The VLA radio data lag the mm/sub-mm (where the lag appears to be variable among the flares; 20–45 min & 40–75 min between 350 GHz and the 18–26 GHz & 4–8 GHz bands, respectively), with flares in the 18–26 GHz band rising to a peak of Jy on a timescale of min, and flares in the 4–8 GHz band rising to a peak of mJy on a timescale of min.
Upon comparing the multi-frequency emission, it is clear that the mm/sub-mm data provide a much more extreme view of the flaring activity than the radio emission. In particular, there is more structure present in the mm/sub-mm light curves when compared to the radio light curves. As such, while not immediately apparent in the radio light curves, the mm/sub-mm data suggest that each of the three main flares in the light curves is actually the result of the superposition of emission from multiple flaring components. Additionally, the lower frequency emission in the light curves appears to be a smoothed, delayed version of the high frequency emission (with the flares showing longer rise times at lower frequencies). This emission pattern is consistent with an expanding outflow structure, where the mm/sub-mm emission originates in a region (with a smaller cross-section) closer to the black hole, and has thus not been smoothed out to as high a degree as the radio emission, as the material expands and propagates outwards. Therefore, all of these observations suggest that the emission in our light curves could be dominated by emission from multiple, expanding, discrete jet ejection events (van der Laan, 1966).
Further, we notice that the baseline flux level at which the flaring begins at each frequency in our light curves appears to vary. This suggests that there is an additional frequency-dependent component contributing to our light curves, on top of the discrete jet ejecta. In an effort to determine the origin of this extra emission, assuming that the baseline emission is constant in time, we create a spectrum of this emission by estimating the baseline flux level at each frequency (we performed iterative sigma clipping and take the minimum of the resulting sigma clipped data). This spectrum101010We note that these are only empirical initial estimates of the baseline flux at each frequency, and do not necessarily represent the flux of the compact jet in our model presented in §4.2. is presented in Figure 2, where it appears as though the baseline emission could be described by a broken power-law or a single power-law (with higher frequency emission displaying a lower baseline level than lower frequency emission). This spectral shape, combined with the fact that we observe a strong compact core component (in addition to resolved ejecta components) within simultaneous high resolution radio imaging (Miller-Jones, et al. 2017, in prep), suggests that the baseline emission originates from an underlying compact jet that was not fully quenched.
4.2 V404 Cyg Jet Model
Given the morphology of our light curves outlined in the previous section, we have constructed a jet model for V404 Cyg that is capable of reproducing emission from multiple, repeated, discrete jet ejection events, on top of an underlying compact jet component. We define two coordinate frames, the observer frame and the source frame (at rest with respect to the ejecta components). We will compute our model primarily in the source frames, and then transform back to the observer frame. All variables with the subscript obs are defined in the observer frame. Schematics displaying the geometry of our model from different viewpoints are displayed in Figures 3 & 4.
In our model, the underlying compact jet is characterized by a broken power-law spectrum, where the flux density is independent of time and varies only with frequency according to,
[TABLE]
Here represents the frequency of the spectral break, represents the amplitude of the compact jet at the spectral break frequency, represents the spectral index at frequencies below the break, and represents the spectral index at frequencies above the break. In the case where the spectral break frequency is located below the lowest sampled frequency band, or above the highest sampled frequency band, the underlying compact jet can be characterized by a single power-law spectrum, where . Here represents the amplitude of the compact jet at , and represents the spectral index.
On top of the compact jet, we define a discrete ejection event as the simultaneous launching of two identical, bi-polar plasma clouds (an approaching and receding component). Each of these clouds evolve according to the vdL model (van der Laan, 1966). In this model, a population of relativistic electrons, with a power-law energy distribution (), is injected into a spherical cloud threaded by a magnetic field. The cloud is then allowed to expand adiabatically, while the electrons and magnetic field are assumed to be kept in equipartition. As a result of the expansion, this model predicts the flux density of each cloud will scale as,
[TABLE]
Here indicates the time-dependent radius of the cloud, and the synchrotron optical depth, , at a frequency, , scales as,
[TABLE]
Note that the subscript 0 in all our equations indicates values at the reference frequency111111We defined our reference frequency as the upper-sideband in our SMA data (230 GHz)., at the time (or radius) of the peak flux of the component.
Taking the derivative of Equation 2 with respect to time121212Our expression in Equation 4 differs from that of van der Laan (1966), as he takes the derivative with respect to instead of time, yielding . (or radius), allows us to relate the optical depth at which the flux density of the reference frequency reaches a maximum, , to the power-law index of the electron energy distribution, ,
[TABLE]
Equation 4 has no analytic solution and thus must be solved numerically. Therefore, we choose to leave our model in terms of , and solve for after the fitting process.
To describe the time-dependence of the cloud radius, a linear expansion model is used, according to,
[TABLE]
Here represents the expansion velocity of the cloud, while can be expressed in terms of the distance to the source, , peak flux, , and optical depth, , of the cloud at the reference frequency (Yusef-Zadeh et al., 2008),
[TABLE]
At the same time that the clouds are expanding, they are also propagating away from the black hole at bulk relativistic velocities, along a jet axis that is inclined to the observer’s line of sight (see Figure 3). As such, the emission we observe will have been affected by projection effects, relativistic beaming, and a geometric time delay between the approaching and receding clouds in each ejection event.
To account for these effects, we first assume that the clouds are travelling at a constant bulk velocity, , and that the jet has a conical geometry (with an observed opening angle, ). In turn, the apparent observed velocity across the sky (derived via the transverse Doppler effect) is represented as (Mirabel & Rodríguez, 1999),
[TABLE]
where is the distance travelled by the cloud away from the black hole, represents the ejection time, represents the speed of light, and represents the inclination angle of the jet axis to the line of sight.
Equation 7 can be simplified by substituting in our expression for to yield,
[TABLE]
where the Doppler factor and bulk Lorentz factor are given by and , respectively. The sign convention in the Doppler factor indicates that a should be used for the approaching cloud and a should be used for the receding cloud.
[TABLE]
Combining Equations 8 & 9, and solving for the bulk Lorentz factor, , yields,
[TABLE]
Rearranging Equation 10 (and substituting in ) gives the expansion velocity (to be input into Equation 5) in terms of only the bulk velocity and jet geometry (inclination and opening angle), such that,
[TABLE]
Further, we wish to write our model in terms of only the ejection time (), rather than the time of the peak flux at the reference frequency (), without introducing any additional parameters. Using our definition that at the instant , the two timescales are related by,
[TABLE]
Lastly, we correct for relativistic beaming by applying a factor of (Longair, 2011) to our flux density in Equation 2, according to,
[TABLE]
Here indicates the flux density of the cloud in the observer frame, at the observing frequency , at the observed times since the zero point of our observations, , while indicates the flux density of the clouds in the source frame, at the frequency, , at the times, .
All of the ejection events we model are not correlated, and thus evolve independently of each other. The total observed flux density in our model is represented as,
[TABLE]
4.2.1 Jet Precession
In addition to our VLA, SMA, and JCMT observations, we also obtained simultaneous high angular resolution radio observations with the Very Long Baseline Array (VLBA). Through imaging the VLBA data set in short 2 minute time bins, we resolve multiple discrete ejecta. Our analysis of these VLBA images has shown clear evidence of jet precession, where the position angle of the resolved ejecta change by up to 40 degrees on an hourly timescale (this result will be be reported in detail in Miller-Jones et al. 2017, in prep.). As the emission predicted by our model is highly dependent on the inclination angle of the jet axis, we account for the effect of this rapid, large scale jet precession in our model by allowing our inclination parameter, , to vary between ejection events.
4.2.2 Accelerated Motion
While we have assumed that the jet ejecta are travelling at constant bulk velocities, it is possible that they undergo some form of accelerated motion. To test this hypothesis we generalized our model to allow the input of a custom bulk velocity profile, where we implemented simple velocity profiles to mimic a finite acceleration period where the cloud would approach a terminal velocity (e.g., a linear ramp function, a body subject to a quadratic drag force). However, in all cases, our best fit model either tended towards a constant velocity profile, or would not converge. This result, while not ruling out the possibility of accelerated motion, suggests that any potential acceleration period may have only lasted for a short enough period of time that we are not able to discern the difference between the resulting light curves for the accelerated and constant bulk motion.
4.2.3 Sub-Conical Jet Geometry
While we have assumed that the jet in our model is conical (constant opening angle), it is possible that the jet geometry could deviate from a strictly conical shape (especially on the AU size scales we are probing), where the opening angle (and in turn the expansion speed of the ejecta) could change with time. In particular, if we assume that the jet confinement mechanism is external, then the jet geometry will depend on the adiabatic indices of the two media (i.e., the jet and its surrounding medium). A relativistic plasmon confined by the internal pressure of a terminal spherical wind (made up of a gas) will expand sub-conically, according to . To test this scenario, we modified our model to use the above sub-conical expansion expression in place of Equation 5. In doing this we find that our best fit model still tends toward constant expansion speed/opening angle profiles for all the ejecta. This result, while not ruling out a non-conical jet geometry, could suggest that any deviations from a conical jet shape only occur on sub-AU size scales, probing timescales before the sub-mm emission peaks, and thus we are not able to discern the difference between the resulting light curves for conical/sub-conical jet geometry.
4.2.4 Bi-polar vs. Single-Sided Ejections
Our jet model assumes that each ejection event takes the form of two identical, oppositely directed plasmons. However, in principle our light curves could also be fit with a collection of single-sided ejections. These unpaired components could occur as a result of Doppler boosting of highly relativistic plasmons causing us to observe only the approaching component of an ejecta pair, or intrinsically unpaired ejecta. Our simultaneous VLBA imaging may help distinguish between these two scenarios. We resolve both paired and (possibly131313Given the rapid timescales of the ejections, multiple ejecta can become blended together in these images, making it difficult at times to conclusively identify and track individual components.) unpaired ejecta components in our VLBA images, which could suggest that the emission in our light curves is produced by a combination of bi-polar and single-sided ejection events. Using these VLBA results to include stricter constraints within our model on ejecta numbers, type (single/bi-polar), and ejection times, is beyond the scope of this work, but will be considered in a future iteration of the model.
4.3 Modelling Process and Best Fit Model
Due to the large number of free parameters in our model, we use a Bayesian approach for parameter estimation. In particular, we apply a Markov-Chain Monte Carlo algorithm (MCMC), implemented with the emcee package (Foreman-Mackey et al., 2013), to fit our light curves with our jet model. This package is a pure-Python implementation of Goodman & Weare’s Affine Invariant MCMC Ensemble Sampler (Goodman & Weare, 2010), running a modified version of the commonly used Metropolis-Hastings Algorithm, whereby it simultaneously evolves an ensemble of “walkers" through the parameter space. We use 500 walkers (10 the number of dimensions in our model) for our MCMC runs.
Prior distributions used for all of our parameters are listed in Table 1. We choose physically informative priors that reflect our knowledge of V404 Cyg (or commonly assumed values for BHXBs) where possible, and wide uninformative uniform priors when we have no pre-defined expectation for a specific parameter. For instance, the prior for the inclination angle is set as a truncated normal distribution, centered on 67 degrees (the measured inclination angle of the system), with boundaries of 0 and 90 degrees (allowed values of the inclination angle). On the other hand, the prior for the ejection time is simply a uniform distribution, sampling a wide range of possible times around our best initial guess.
Before running the MCMC, the initial position of the walkers in the parameter space needs to be defined. As the performance of the emcee algorithm tends to benefit heavily from well defined initial conditions, we do an initial exploration of the parameter space using a harmony search global optimization algorithm141414Implemented in the python package, pyHarmonySearch; https://github.com/gfairchild/pyHarmonySearch . This metaheuristic algorithm, that is similar to, but much more efficient than a brute force grid search method (which would not be computationally feasible in this case), yields a reasonable initial guess for our model, and we place our walkers in a tight ball around this initial guess in the parameter space.
As our jet model can predict emission at multiple frequencies, to reduce the degeneracy in our model, we choose to simultaneously fit all of our multi-frequency data sets, except for the JCMT 666 GHz data set, due to its sparser sampling and larger uncertainty in flux calibration (see the Appendix for details). To do this, we use an iterative process whereby we start with our reference frequency data set, run the MCMC (the walkers are evolved over a series of steps, where the first 500 step “burn in" period is not retained) until convergence is reached, and use the final position of the walkers for the first run as the initial guess for the next run of the MCMC, which will include increasingly more data sets in the fit. To monitor the progress of the MCMC and ensure that correct sampling was occurring, we checked that the acceptance fraction stayed within the suggested bounds (between 0.25 and 0.75). Our criteria for convergence requires that the positions of the walkers are no longer significantly evolving. We determine whether this criteria is met by monitoring the chains of each of the walkers through the parameter space, and ensuring that, for each parameter, the intra-chain variance across samples is consistent with the inter-chain variance at a given sample.
Using the multi-dimensional posterior distribution output from the converged MCMC solution, we create one dimensional histograms for each parameter. The best fit result is taken as the median of these distributions, and the uncertainties are reported as the range between the median and the 15th percentile (-), and the 85th percentile and the median (+), corresponding approximately to errors. All of the best-fit parameters and their uncertainties are reported in Table 2. Figure 5 & 6 show the best fit model overlaid on our multi-frequency light curves. Additionally, with our multi-dimensional posterior distribution we can explore the possible two-parameter correlations for our model, where a significant correlation between a pair of parameters can indicate a model degeneracy or a physical relationship between the parameters. In the Appendix section we show correlation plots (Figure 11), along with the one-dimensional histograms, for pairs of parameters for which we find a correlation, and discuss the significance of such a correlation.
Within the Bayesian formalism, the uncertainties reported in Table 2 are purely statistical, and only represent the credible ranges of the model parameters under the assumption that our model is correct. Given the residuals with respect to the optimal model (Figure 5 bottom panel), the observations contain physical or observational effects not completely accounted for in our model. To factor in how well our chosen model represents the data, we estimated an additional systematic error for our parameters (displayed in Table 3 of the Appendix). To do this we rerun our MCMC, starting from the best fit solution, with an extra variance parameter (effectively modelling all the physical/observational effects not included in our model) in our log probability for each frequency band. This variance is equivalent to the square of the mean absolute deviation of the residuals with respect to our optimal model at each frequency (difference between the best fit model and the data). The resulting uncertainties in the parameters after this extra MCMC run will reflect the full (statistical + systematic) uncertainties.
Our broad frequency coverage, in particular the high sub-mm frequencies, is crucial to the success of our modelling. Detailed substructure detected in the sub-mm bands can be used to separate out emission from different ejections, where their lower frequency counterparts are smoothed out and blended together. As such, modelling the lower frequency emission would not be possible without the critical information the high frequency sub-mm emission provides and vice versa.
5 Discussion of the Best Fit Model
Our jet model for V404 Cyg, with a total of 8 bi-polar ejection events on top of an underlying compact jet, is able to reproduce the emission in all of our observed frequency bands, matching the flux levels, time lags between frequencies, and the overall morphology remarkably well. With such a large sample of jet ejecta, we can probe the intrinsic ejecta properties, and the distribution of these properties between the different ejection events. In particular, our model characterizes the bulk speeds, peak fluxes, the electron population injected during each event, and the jet geometry, all of which we find can vary between events, with bulk speeds of c, peak fluxes of mJy, electron energy distribution indices of (corresponding to ), and observed opening angles of ∘. In the following sections we discuss these ejecta parameters and what they can tell us about jet speeds, energetics, mass loss, and geometry. Additionally, we draw comparisons between the V404 Cyg ejection events and the jet oscillation events in GRS 1915+105, as well as other multi-wavelength observations of V404 Cyg.
5.1 Jet Speeds
The bulk speeds of jet ejecta measured in BHXBs151515An important caveat when considering the value of the bulk Lorentz factor (), estimated using proper motions of discrete jet ejecta, is that depends strongly on the assumed distance to the source (Fender, 2003). While the distance is well known for V404 Cyg, this is not the case for the majority of BHXBs, and as a result constraints on in these systems typically represent lower limits. can vary from system to system (e.g., in SS 433; Hjellming & Johnson 1981, in V4641 Sgr; Hjellming et al. 2000), where some systems that are known to enter high luminosity states, like V404 Cyg, have been shown to launch jet ejecta with (e.g., GRO J1655-40; Hjellming & Rupen 1995). However, in V404 Cyg we find that the bulk speeds of our modelled ejecta are quite low, with bulk Lorentz factors of only (excluding ejection 5; see footnote c in Table 2 for details).
Moreover, V404 Cyg shows bulk speeds that vary substantially between ejection events, on timescales as short as minutes to hours. There is some evidence in the literature that jet speeds can vary within a BHXB161616There is also evidence of jet speeds varying in neutron star XBs, most notably, Sco X-1 (Fomalont et al., 2001a, b) and Cir X-1 (Tudose et al., 2008). source. For example, Blundell & Bowler (2005) find small variations in jet speed up to 10% in SS 433, jet speeds have been reported to vary between outbursts of H1743-322 (Corbel et al. 2005; Miller-Jones et al. 2012b), and varying proper motions have been measured in GRS 1915+105 (Miller-Jones et al., 2007). However, no other source has shown variations as large, or on as rapid timescales as V404 Cyg.
Performing a Monte Carlo Spearman’s rank correlation test, we find no correlation between jet speed and ejection time, where, for instance, the bulk speed of the ejections (i.e., ) increased or decreased throughout our observation period. However, we find a potential correlation (Spearman coefficient of with a p-value of ) between bulk speed and peak flux of our modelled ejecta, where brighter ejecta tend to have higher speeds. This correlation is consistent with what was seen in H1743-322, where higher bulk ejecta speeds corresponded to higher radio luminosity measurements (Corbel et al. 2005; Miller-Jones et al. 2012b).
The factors that govern jet speed in BHXBs are not well understood, but our measurements of surprisingly slow speeds, which can vary between sequential jet ejection events, suggest that the properties of the compact object (i.e., black hole mass) or peak luminosity of the outburst are likely not the dominant factors that affect jet speed.
Additionally, given the varying bulk speeds between the ejection events, it is plausible that later, faster ejections could catch up to earlier, slower ejections. Such a collision between ejecta may result in a shock that could be as bright or even brighter than the initial ejections, and in turn produce a flaring profile that could mimic a new ejection event. While including ejecta collisions in our model is beyond the scope of this work, we briefly consider the possibility here by examining the bulk motion of all of the ejections. We find that a collision between ejection 3 and ejection 2 would occur at 11:30 (if they were ejected at the same PA), which is very close to the predicted ejection time of ejection 4. Moreover, ejection 4 has a bulk speed which is in between the bulk speeds of ejection 2 and 3, as we might expect for the bulk motion of the plasmon after such a collision. However, given that the jet appears to be rapidly precessing in V404 Cyg (Miller-Jones et al. 2017, in prep), ejection 2 and ejection 3 are launched at very different inclination angles, which would prevent such a collision from occurring. Therefore, given the precessing jet, we find this collision scenario unlikely.
5.2 Jet Energetics, Mass Loss, and Particle Acceleration
In our model we assumed that the radiating electrons follow a power-law energy distribution. The power-law index of this distribution, , informs us about the population of accelerated electrons initially injected into each discrete jet component, where the value of this energy index is governed by the electron acceleration mechanism. Fermi acceleration by a single shock can produce values of , which are typically found in XRBs (Blandford & Eichler, 1987; Bell, 1978; Markoff et al., 2001). However, the energy index can take on a wider range of values under certain conditions, where for example, lower values of (which result in a more asymmetric flare profile) can be produced if the acceleration occurs in multiple shocks (Melrose & Pope, 1993), or if the electrons carry away kinetic power from the shock (Drury & Volk, 1981), and higher values of could be produced in the presence of oblique shocks (although this case requires highly relativistic shocks to produce large values; Ballard & Heavens 1992). Distributions with values of are nearly indistinguishable from a thermal (Maxwellian) distribution, which in the shock acceleration paradigm, implies very little acceleration has occurred (a shock essentially takes an input thermal distribution of electrons and builds a power-law distribution up over time). Magnetic reconnection in a relativistic plasma is another viable mechanism that can accelerate electrons into distributions with similar values to shock acceleration. In this case, smaller values can be produced in the case of a strongly magnetized plasma (; where represents the magnetization parameter), and larger values can be produced in a weakly magnetized plasma (; Sironi & Spitkovsky 2014; Guo et al. 2014; Sironi et al. 2016). In either theory of particle acceleration, we would expect a link between the speed (for shock acceleration) or magnetization (for magnetic reconnection), and the energy index, .
The energy indices of our modelled ejecta appear to vary between sequential ejection events, with (where we find no clear correlation between values and jet speed). These values could be produced by shock acceleration or magnetic reconnection (under the right conditions), although we would need to invoke different mechanisms to produce distributions in both the very low and very high regimes (e.g., 1.4 in ejection 7, and 5.7 in ejection 2), which is not entirely physical for a single source. Further, this significant range seen in our energy indices suggests that our model may not be capturing all of the complexities of these ejection events, where the more extreme values of the energy index could be mimicking the effect of physics that has not been included in our model. For instance, the vdL model assumes equipartition, but as the plasmons expand they must do work, which will result in some of the magnetic field dissipating into kinetic or thermal pressure, and in turn, the assumption of equipartition may break down. Simplifications in our model such as this could also explain the lack of expected correlation between our energy indices and the speed of the ejecta. A more rigorous treatment, which, for example, calculates the full synchrotron flux (and does not rely on the equipartition assumption), is beyond the scope of this work, but will be considered in future iterations of this model.
For synchrotron emitting clouds of plasma injected with our measured electron distributions, we estimate that the minimum energy171717In our minimum energy calculations, we perform the full calculations outlined in Longair (2011), where we integrate the electron energy distribution from to . The minimum frequency represents the lowest radio detection with LOFAR on June 23 & 24 (Broderick et al., 2015), and the maximum frequency represents our highest frequency sub-mm detection. When we consider an electron-proton plasma, we assume the ratio of the energy in the protons over that of the electrons is . needed to produce each of our modelled ejection events range from , with minimum energy magnetic fields181818We note that while these calculation assume equipartition, the system could be far from equipartition. In this case the magnetic field would not necessarily be equivalent to the minimum energy field, but rather could be either much higher or much lower. on the order of a few Gauss ( G). Taking into account the duration of each event, these energies correspond to a mean power into each event ranging from . Due to the slow bulk speeds of the ejecta, including the kinetic energy from the bulk motion (in an electron-positron plasma ) yields only slightly higher values of . The minimum energy and power released within each of our modelled ejection events is comparatively lower than estimated for other major ejection events in BHXBs (; e.g., Fender et al. 1999b and ; e.g., Fender et al. 1999c; Brocksopp et al. 2007; Curran et al. 2014). This difference is dominated by the difference in the estimated size of the emitting region, where the radii of our modelled ejecta are smaller than is normally estimated for major ejection events, and the low bulk speeds, which result in a much smaller kinetic energy contribution. Considering that the flaring activity in V404 Cyg lasted weeks (and assuming our observations to be representative of this entire period), we estimate that the the total (minimum) energy (radiative + kinetic) released into jet ejections over the full flaring period is , which is more on par with typical energies estimated for major ejection events in BHXBs. This total energy is also comparable to that carried by the accretion disc wind ()191919A rough estimate of the energy lost in the accretion disc wind is equivalent to . Using and (Munoz-Darias et al., 2016a), we estimate ..
If we assume that the jet ejecta contain some baryonic content, in the form of one cold proton for every electron, we calculate that the mean power into each event (including the kinetic energy from bulk motion) ranges from . In this baryonic case, we estimate a total mass lost through the jet in our observation period of (corresponding to over the 2 week flaring period). To compare this jet mass loss to the mass accreted onto the black hole, we follow a procedure similar to Munoz-Darias et al. (2016a), using simultaneous INTEGRAL X-ray observations (only including the harder ISGRI bands, ranging from 25-200 keV) to calculate the total energy radiated (integrated X-ray luminosity) during our observations. To do this we convert the count rate into flux in the 10–1000 keV band using a power law model with photon index , and approximating the integral as a sum (, where is the weighted mean, is size of the time bins, and is the total observation time). Assuming an accretion efficiency of 0.1, we calculate a total mass accreted during our observations of . Therefore, the mass lost in the jet is a small fraction of the total mass accreted, , and much less than the mass estimated to be lost in the accretion disc wind (; Munoz-Darias et al. 2016a).
5.3 Jet Geometry and Ejecta Size Scale
Measurements of jet geometry in BHXBs, in particular the observed opening angle, only exist for a handful of systems, where all but one are upper limits (e.g., see Table 1 in Miller-Jones et al. 2006, as well as Yang et al. 2010 and Rushton et al. 2017 for recent measurements in XTE J1752-223 and XTE J1908+094). Our simultaneous light curve modelling technique allows us to directly derive the first measurements of the jet geometry in V404 Cyg, where we model observed jet opening angles of ∘. These measurements are consistent with the opening angle estimates for the other BHXB systems with constraints, where the majority show .
With the opening angles, we can estimate the level of confinement of the jets in V404 Cyg by solving for the intrinsic expansion speed (using Equation 11; see last column of Table 2) of our modelled ejecta ( indicates freely expanding components, where represents the speed of sound in a relativistic gas). We find intrinsic expansion speeds of c, indicating a highly confined jet in V404 Cyg. There are many possible mechanisms that could be responsible for confining the jets in V404 Cyg. In particular: the jet could be inertially confined (Icke et al., 1992), where the ram pressure of the strong accretion disc wind detected in V404 Cyg (Munoz-Darias et al., 2016a) could inhibit the jet ejecta expansion202020Although, we note that if the confinement is external, this would suggest that a very large amount of pressure surrounds the ejections. If this is supplied solely by the ram pressure from an accretion disc wind, then the mass-loss rate (proportional to the velocity ratio of the ejections to the wind) would be unrealistically large (i.e., greater then the mass accretion rate).; the jet could be magnetically confined by a toroidal magnetic field (Eichler, 1993); the jet could contain cold protons, which may impede the jet ejecta expansion (Miller-Jones et al., 2006); or a combination of these different mechanisms could be at work.
Further, as we alluded to in the previous section, the initial radii of the jet ejecta (i.e., the radii of the ejecta at the time the sub-mm emission peaks) estimated by our model are noticeably smaller than those typically estimated for major ejection events in other BHXBs. This is likely a result of the much slower expansion velocities () we find for the V404 Cyg ejecta. In particular, we infer a range of initial radii for our ejecta ranging from . These radii appear to remain similar (to within a factor of 2) between ejection events.
5.4 Underlying Compact Jet
In addition to the jet ejecta component, we observe an extra constant flux component in our light curves, which varies with frequency. Due to the shape of our estimated spectrum of this emission (see Figure 2) and the strong compact core jet present throughout the span of our simultaneous VLBA imaging (Miller-Jones et al. 2017, in prep.), we interpret this extra flux term as emission from an underlying compact jet. We believe that this compact jet was switched on during the launching of the multiple discrete ejection events. In our best fit model, this compact jet emission is characterized by a single power-law spectrum, with an optically thin spectral index of .
Our suggestion of a compact jet, that has not been fully quenched, is in agreement with the findings of Sanchez-Fernandez et al. (2016), who show that V404 Cyg never fully reached a soft accretion state (where we would likely expect strong quenching of the compact jet; e.g., Rushton et al. 2016), but rather remained in either a harder intermediate or very high state during our observations. Under this interpretation, based on our lowest radio frequency measurement, we place limits on the optically thick to thin jet spectral break frequency of GHz, and flux at the spectral break of mJy. However, we note that simultaneous VLITE observations (Kassim et al., 2015) detect V404 Cyg at a total time-averaged flux density of mJy at 341 MHz. Given that our best fit model predicts a maximum jet ejecta flux component of mJy at 341 MHz, it is clear that the 341 MHz compact jet component cannot lie along the single power-law stated above. As such, the spectral break would therefore occur within the range of GHz, which is significantly lower than previous estimates for V404 Cyg made during the hard accretion state (; Russell et al. 2013a, b). This evolution in the location of the spectral break is consistent with the pattern suggested by recent observations (e.g., Corbel et al. 2013; van der Horst et al. 2013; Russell et al. 2014) and MHD simulations (Polko et al., 2014), where, as the mass accretion rate increases during softer accretion states of BHXB outbursts (which usually occur at high luminosities; Koljonen et al. 2015), the jet spectral break moves toward lower radio frequencies prior to the jet switching off (or at least fading below our detection limits).
Up to now we have only considered the compact jet and the ejection events as separate entities. In the presence of explosive, energetic ejection events, we might expect a compact jet to be disrupted. In particular, as the ejecta collide with the pre-existing compact jet, a shock would likely develop, due to the difference in bulk speeds between the two. In this situation, if the compact jet rapidly re-establishes itself after being destroyed by ejecta (before the ejecta propagate far enough away from the black hole to be resolved), we would observe a compact core jet which appears to never shut off. Therefore, we believe it is plausible that a compact jet is being repeatedly destroyed and re-established (on rapid timescales) following ejection events in V404 Cyg. Further, the emission from such a shock interaction could display an optically thin spectrum (similar to the interaction between the discrete ejecta and the surrounding ISM; e.g., Corbel et al. 2004), like the one we observe for our baseline emission component. Thus, while we interpreted the baseline emission in our light curves as originating only from a compact jet, emission from a possible interaction of the jet ejecta with this compact jet, and/or continuous lower-level, fainter jet ejecta that never get resolved, could also be contributing to the baseline flux level we observe.
Moreover, in our model we have assumed that the compact jet flux component is constant in time. However, as the accretion rate (and in turn the jet power) changes, the flux of a compact jet is expected to change as well (Russell et al., 2014). If we consider the erratic X-ray behaviour observed in the source (which presumably traces a rapidly changing accretion rate), it is plausible that the compact jet component could in fact be variable as well. Exploring the possibility of a variable compact jet component in our model is left for future work.
5.5 Ejecta Time Lags
Our model predicts that the intrinsic time lag (in the source frame) between a certain frequency () and the reference frequency (), is represented by,
[TABLE]
where the observed time lag can be obtained through the transformation, .
Figure 7 shows the observed time lags, predicted by our model, between each frequency band and the ejection time, for the approaching (top panel) and receding (bottom panel) components. The time lags are clearly variable between different ejection events (e.g., min between the ejection and our reference frequency, 230 GHz), which is a result of the varying ejecta properties (i.e., , , ).
Further, it is interesting to note that, for a different flaring event that occurred days after our data set, Shahbaz et al. (2016) measured a time lag of 2.0 hours & 3.8 hours between the predicted ejection time (indicated by an r’-band polarization flare, which these authors suggest could be the signature of the launching of major jet ejection event) and the flare peaks at 16 GHz & 5 GHz, respectively. These lags are slightly higher than predicted for the approaching components of our modelled ejection events, but appear to share a similar slope across frequencies.
5.6 Comparison to GRS 1915+105
GRS 1915+105 is the only other BHXB in which a similar multi-frequency variability pattern to that seen in V404 Cyg has been reported. While flaring activity has been seen in other systems, it is often only detected in one frequency band (e.g., V4641 Sgr in optical; Uemura et al. 2004), or the flares in question evolved over much longer (days rather than minutes/hours) timescales (e.g., 4U 1630-47 in radio/X-ray; Hjellming et al. 1999). GRS 1915+105 has displayed some correlated radio, sub-mm, and IR flares (with lower frequency emission delayed from higher frequency emission), which repeated every minutes for a day period (Fender & Pooley, 2000). While no discrete components were resolved with VLBI during the events, the similar rise and decay times of flares at different frequencies suggest that adiabatic energy losses, likely during the expansion of discrete components, played a key role in determining the flaring profiles of these events. In fact, Mirabel et al. (1998) found that the timing of the radio emission during these events was consistent with synchrotron emission from adiabatically expanding plasma clouds, where each event required an energy input of , and carried an estimated mass of . Many studies of these jet ejection events suggest that they occur as a result of instabilities causing the repeated ejection and refilling of the inner accretion disc or coronal flow (e.g., Belloni et al. 1997; Nandi et al. 2001; Vadawale et al. 2001).
In V404 Cyg our modelled ejection times appear to occur on a similar rapid timescale as seen in GRS 1915+105, where we observe groups of 2-4 ejections (separated by at most minutes), followed by longer periods of up to hour between groups (see Figure 8). Each group of ejections seems to correspond to a large flaring event in the light curves. Our estimates of the energetics and mass-loss of the V404 Cyg events (§5.2) are also similar to those estimated for the oscillation events in GRS 1915+105, where both are consistent with being smaller-scale analogues of major ejection events seen in other BHXBs. Further, Naik et al. (2001) suggested that multiple ejections in GRS 1915+105 could manifest as a single radio flare, similar to the ejection groupings we see in V404 Cyg. However, a noticeable difference in the timing of the V404 Cyg and GRS 1915+105 events is that the V404 Cyg events are not as quasi-periodic (i.e., they do not occur on as regular intervals) when compared to the GRS 1915+105 events, which occurred every min (Fender & Pooley, 2000). The absence of quasi-periodicity in the V404 Cyg events could indicate that the jet production process is not as stable in V404 Cyg as it was during the GRS 1915+105 events.
The similarity between the morphology, duration and energetics of the rapid ejection periods in V404 Cyg and GRS 1915+105 suggests that the events may have a common origin, possibly in the repeated ejection and refilling of some reservoir in the inner accretion flow. This hypothesis is consistent with the recent finding of Radhika et al. (2016), who report the non-detection of the disc component in the X-ray spectra following major radio flares in V404 Cyg. Although, given the large intrinsic absorption (Motta et al., 2016b) seen in V404 Cyg during the outburst, it is conceivable that we may not have been able to detect the soft disc component, even if it was present. Nevertheless, as both V404 Cyg and GRS 1915+105 are long period systems, with large accretion discs, a key ingredient in fuelling rapid, repeated ejection events may be a large accretion disc (as suggested by Kimura et al. 2016; Munoz-Darias et al. 2016a).
5.7 Alternative Emission Models
Other than the vdL model, an alternative emission model that has been used to reproduce flaring light curves in XRBs is the shock-in-jet model (Marscher & Gear, 1985; Bjornsson & Aslaksen, 2000; Turler et al., 2000). This analytical model, while traditionally favoured for extragalactic sources, has been successfully applied to flaring events in Cyg X–3 (Lindfors et al., 2007; Miller-Jones et al., 2009a), (Stevens et al., 2003), and GRS 1915+105 (Turler et al., 2004). The shock-in-jet model considers a shock wave travelling downstream in a jet flow as the source of each flare in the light curve. Each shock wave will evolve through three different phases; (1) Compton losses dominate, (2) synchrotron losses dominate, and (3) adiabatic losses dominate. The main differences between the shock-in-jet model and the vdL model are that the shock-in-jet model considers an initial phase where Compton losses dominate over adiabatic losses, all shock wave events are self-similar, and the electron energy scales differently when compared to the vdL model (shock-in-jet flow expands in 2D, ; vdL cloud expands in 3D, ). These differences will result in a different flare profile between models, where for the same electron population (i.e., same value), the shock-in-jet model flares will show a much shallower decay, and the peak fluxes at frequencies that are initially optically thin (likely IR and above) will be smaller than predicted by the vdL model (which will over predict peak fluxes at these frequencies).
As our adapted vdL model is able to reproduce all our light curves (at 7 different frequencies) remarkably well, and simultaneous VLBA imaging resolves multiple, discrete components (Miller-Jones et al. 2017, in prep), we favour the expanding plasmon model over the shock-in-jet model for the V404 Cyg events (although we can not rule out the shock-in-jet model).
However, for the GRS 1915+105 oscillation events, the emission has been shown to be consistent with both an expanding plasmon model (Mirabel et al. 1998; although we note that these authors only model a single flaring event, and did not include any relativistic/projection effects in their model) and a shock-in-jet model (Turler et al., 2004). If the GRS 1915+105 events are in fact a result of shock waves rather than expanding plasmons, this could explain the notable differences to the V404 Cyg events, namely the quasi-periodicity and the lack of VLBI resolved components212121Although, we note that these GRS 1915+105 oscillation events were only observed with MERLIN (Fender et al., 1999c), which does not have the resolution to see ejection events of a few mas in size (like those of V404 Cyg).. Additionally, as Turler et al. (2004) point out, the shock-in-jet model is still consistent with the scenario that these oscillation events originated with instabilities in the inner accretion disc, as these instabilities could be the catalyst that leads to an increased injection rate of material at the base of the jet, and in turn a downstream shock wave.
5.8 Connection to X-ray & OIR
If the jet ejection events in V404 Cyg are linked to processes occurring in the accretion flow, we might expect our predicted ejection times to correlate well with X-ray/OIR emission. For instance, in GRS 1915+105, IR and radio flares (which are presumably tracers of the ejection events) followed an X-ray peak and occurred during a period of spectral softening (dips in hard X-ray emission). However, the connection is not as clear in V404 Cyg. Figure 8 displays our predicted ejection times on top of simultaneous X-ray222222All X-ray data presented in this paper are taken from the INTEGRAL public data products available at http://www.isdc.unige.ch/integral/analysis#QLAsources. (Rodriguez et al., 2015) and OIR (Kimura et al., 2016) emission.
Flares in the OIR light curve appear to coincide with flares in the X-ray light curves. However, an unfortunate gap in the OIR coverage makes it difficult to confirm that such a pattern holds for the final X-ray flare. In terms of our modelled ejection times, we may be able to tentatively match groups of ejections with specific X-ray/OIR peaks, and possibly local dips in hardness (where the start/end of a steep gradient in hardness appears to correspond to ejections). But it is puzzling that the group which contains the largest number of ejections and produces the largest sub-mm flares appears to be connected to the X-ray flare with the smallest amplitude (although, if an X-ray flare is indicative of a strong dissipative process, more energy dissipated in the X-ray implies less energy would be available to the jets, and vice versa). Further, the second X-ray flare appears to have no jet ejecta counterpart.
Given the extremely high intrinsic absorption during this time period (Motta et al., 2016b), it is entirely possible that the flaring in the X-ray light curves is not always dominated by intrinsic source variation, but rather dependent on how much of the inner accretion flow is obscured. This effect was seen in the 1989 outburst, where large changes in column density were determined to be the origin of some of the extreme X-ray variability observed (Oosterbroek et al., 1997; Zycki et al., 1999). Thus, even if the jet ejections are linked to processes in the accretion flow, we may not expect to see a clear correlation between our jet ejections and the X-ray/OIR emission. On the other hand, if the high absorption reduced the X-ray flux artificially, we would expect the high energy bands (60-200 keV) to be less affected than the lower energy bands (5-10 keV), which does not seem to be the case here. Therefore, the nature of the connection (if any) between our jet ejections and the X-ray/OIR emission is still not fully understood.
5.9 The Critical Sub-mm Perspective
Traditionally XRB jet studies have been dominated by radio frequency observations, such that there only exists a limited set of XRBs that have been observed at mm/sub-mm frequencies (e.g., Paredes et al. 2000; van der Horst et al. 2013; Russell et al. 2013b; Fender et al. 2001; Tetarenko et al. 2015a). When considering time-resolved ( day cadence) mm/sub-mm observations this number decreases to two (i.e., GRS 1915+105; Pooley & Fender 1997, Cygnus X-3; Baars et al. 1986; Fender et al. 1995). However, our work in this paper has clearly shown the vital importance of high time resolution mm/sub-mm data in XRB jet studies. In particular, the mm/sub-mm bands can be used to isolate emission from different flaring events in the light curves, while the lower frequency counterparts of these events tend to be smoothed out and blended together. As such, we find that radio frequency observations alone can often be misleading, especially in terms of identifying and pinpointing the timing of individual rapidly variable flaring events. Including mm/sub-mm monitoring during future XRB outbursts will continue to add key insight to our understanding of jet behaviour.
6 Summary
In this paper we present the results of our simultaneous radio through sub-mm observations of the BHXB V404 Cyg during its June 2015 outburst, with the VLA, SMA and JCMT. Our comprehensive data set, taken on 2015 June 22 ( week following the initial detection of the outburst), extends across 8 different frequency bands (5, 7, 21, 26, 220, 230, 350, and 666 GHz). Using custom procedures developed by our team, we created high time resolution light curves of V404 Cyg in all of our sampled frequency bands. In these light curves, we detect extraordinary multi-frequency variability in the form of multiple large amplitude flaring events, reaching Jy level fluxes.
Based on the overall morphology, we postulate that our light curves were dominated by emission from a relativistic jet. To understand the source of the emission we constructed a detailed jet model for V404 Cyg. Our model is capable of reproducing emission from multiple, discrete, bi-polar plasma ejection events, which travel at bulk relativistic speeds (along a jet axis inclined to the line of sight), and evolve according to the van der Laan synchrotron bubble model (van der Laan, 1966), on top of an underlying compact jet. Through implementing a Bayesian MCMC technique to simultaneously fit all of our multi-frequency light curves with our jet model, we find that a total of 8 bi-polar ejection events can reproduce the emission we observe in all of our sampled frequency bands.
Using our best fit model to probe the intrinsic properties of the jet ejecta, we draw the following conclusions about the ejection events in V404 Cyg:
- •
The intrinsic properties of the jet ejecta (i.e., speeds, peak fluxes, electron energy distribution indices, opening angles) vary between different ejection events. This results in varying time lags between the flares produced by each ejection at different frequencies.
- •
The ejecta require (minimum) energies on the order of . When taking into account the duration of each event, these energies correspond to a mean power into the ejection events of .
- •
The ejecta carry very little mass (% ), especially when compared to that carried by the other form of outflow detected in V404 Cyg, the accretion disc wind (). However, despite carrying much less mass, we estimate that the ejecta carry similar energy to that of the accretion disc winds.
- •
We place the first constraints on jet geometry in V404 Cyg, where we find that V404 Cyg contains a highly confined jet, with observed opening angles of the ejecta ranging from . While we can not pin down the main jet confinement mechanism in V404 Cyg, it is possible that the ram pressure of the strong accretion disc wind detected in V404 Cyg (Munoz-Darias et al., 2016a) could contribute to inhibiting the jet ejecta expansion, and thus be a key cause of the highly confined jet in this system.
- •
The ejecta travel at reasonably slow bulk speeds, that can vary substantially between events, on timescales as short as minutes to hours ().
- •
Brighter ejections tend to travel at faster bulk speeds.
- •
Our modelled ejection events appear to occur in groups of 2-4 ejections (separated by at most minutes), followed by longer periods of up to hour between groups.
- •
The rapid timescale of the ejections is similar to the jet oscillation events observed in GRS 1915+105. Although the V404 Cyg events do not occur on as regular intervals as the GRS 1915+105 events, possibly suggesting the jet production process is not as stable in V404 Cyg.
- •
We can tentatively match groups of ejections with peaks in simultaneous X-ray/OIR emission. However, the nature of the connection (if any) between our modelled ejection events and X-ray/OIR emission is still not completely clear.
Based on these conclusions, it appears as though the V404 Cyg ejection events are smaller-scale analogues of major ejection events, typically seen during the hard to soft accretion state transition in BHXBs. Given the similarity between these rapid ejection events in V404 Cyg and those seen in GRS 1915+105, we postulate that the ejection events in both systems may have a common origin, in the repeated ejection and refilling of some reservoir in the inner accretion flow. This suggests that, in agreement with the findings of Kimura et al. (2016) & Munoz-Darias et al. (2016a), the presence of a large accretion disc in both systems may be a key ingredient in producing these rare, rapid ejection events.
Overall, the success of our modelling has shown that, multiple expanding plasmons, on top of a compact jet, is a good match to the emission we observe from V404 Cyg in multiple frequency bands. However, it is also apparent from our results that some simplifications within our model may not fully capture all of the physics of these ejection events (e.g., assuming equipartition, assuming a constant flux from the compact jet), and future iterations of this model will work to address these assumptions and explore their effect on the ejecta properties.
In this work we have demonstrated that simultaneous multi-band photometry of outbursting BHXBs can provide a powerful probe of jet speed, structure, energetics, and geometry. Additionally, our analysis has revealed that the mm/sub-mm bands provide a critical new perspective on BHXB jets (especially in the time-domain) that can not be achieved with radio frequency observations alone. Future high time resolution, multi-band observations of more systems, including the mm/sub-mm bands, have the potential to provide invaluable insights into the underlying physics that drives jet behaviour, not only in BHXBs but across the black hole mass and power scale.
Acknowledgements
We extend our sincere thanks to all of the NRAO, SMA, and JCMT staff involved in the scheduling and execution of these observations. Without their tireless hard work and constant support during this observing campaign we would never have obtained such an extraordinary data set. We offer a special thanks to Iain Coulson for continuing to share his JCMT expertise. We thank M. Kimura et al. for sharing their OIR data. AJT thanks Eric Koch for many helpful discussions on MCMC implementation and cloud computing, and Patrick Crumley for his helpful comments and feedback on particle acceleration mechanisms. Also, many thanks to Christian Knigge for creating the V404 mailing list and Tom Marsh for creating the V404 observations website. Both of your efforts made it much easier to organize coordinated multi-frequency observations of this outburst. AJT is supported by an NSERC Post-Graduate Doctoral Scholarship (PGSD2-490318-2016). AJT, GRS, and EWR are supported by NSERC Discovery Grants. JCAMJ is the recipient of an Australian Research Council Future Fellowship (FT140101082). SM acknowledges support from VICI grant Nr. 639.043.513/520, funded by the Netherlands Organisation for Scientific Research (NWO). TDR acknowledges support from the Netherlands Organisation for Scientific Research (NWO) Veni Fellowship, grant number 639.041.646. Cloud computing time on Amazon Web Services, used for the development and testing of our CASA timing scripts, was provided by the SKA/AWS Astro-Compute in the Cloud Program. Additionally, we acknowledge the use of Cybera Rapid Access Cloud Computing Resources, and Compute Canada WestGrid Cloud Services for this work. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Sub-millimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics, and is funded by the Smithsonian Institution and the Academia Sinica. The James Clerk Maxwell Telescope is operated by the East Asian Observatory on behalf of The National Astronomical Observatory of Japan, Academia Sinica Institute of Astronomy and Astrophysics, the Korea Astronomy and Space Science Institute, the National Astronomical Observatories of China and the Chinese Academy of Sciences (Grant No. XDB09000000), with additional funding support from the Science and Technology Facilities Council of the United Kingdom and participating universities in the United Kingdom and Canada. The authors also wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.
Appendix A Image Weighting Scheme
As we are imaging the source on very short timescales, the uv-coverage in each time-bin will be limited. While we do not need to worry about the lack of uv-coverage affecting the fidelity of the images, as the source is point-like at the VLA and SMA resolutions, the side-lobe levels may be a concern. In particular, if the amplitude is changing significantly in each time bin, this implies that we cannot deconvolve the side-lobes properly. As such, the choice of weighting scheme used while imaging could affect the quality of the images, and in turn the flux density measurements for each time bin. While the side-lobe level is not much of a concern for the VLA, which has reasonably good instantaneous uv-coverage, the SMA is only an 8-element interferometer. In this case, imaging the source with a more uniform weighting scheme minimizes the side-lobe level, and could improve the quality of the images in each time bin. On the other hand, imaging with a natural weighting scheme would maximize sensitivity, leading to lower rms noise levels. After testing different weighting schemes we find that the choice of weighting had very little effect on the output SMA light curves, where any differences in the flux measurements in each time bin were well within the rms noise. We find that the natural weighting scheme led to lower rms noise and slightly higher dynamic range in the majority of the time bin images. Therefore, we opted to use natural weighting, as the side-lobe level/rms noise trade-off appears to be optimized for natural weighting.
Appendix B Calibrator Light Curves
Given the large flux variations we detected in our data of V404 Cyg, we wished to check the flux calibration accuracy of all of our observations on short time scales, and ensure that the variations we observed in V404 Cyg are dominated by intrinsic variations and not atmospheric or instrumental effects. Therefore, we ran our custom procedures to extract high time resolution measurements from our data (see §3 for details) on all of our calibrator sources. Figure 9 displays target & calibrator light curves at all frequencies.
We find that all of our interferometric calibrator sources and our JCMT 350 GHz calibrator display relatively constant fluxes throughout our observations, with any variations (/ of the average flux density at radio/(sub)-mm frequencies) being a very small fraction of the variations we see in V404 Cyg. However, our JCMT 666 GHz calibrator scan shows noticeably larger scale variations ( of its average flux level). While these larger variations are not unexpected at this high frequency, as the atmosphere is much more opaque, when combined with the fact that higher noise levels at this frequency prevent us from sampling timescales shorter than 60 seconds, we choose to not include the 666 GHz data set in our modelling (although see Appendix B below for a discussion of how well our best fit model agrees with the 666 GHz data).
Overall, based on these results, we are confident that the high time resolution light curves of V404 Cyg used in our modelling are an accurate representation of the rapidly changing intrinsic flux of the source.
Appendix C JCMT SCUBA-2 666 GHz Model Comparison
While we did not include the JCMT SCUBA-2 666 GHz data in our model fitting, it is still of interest to compare our best fit model prediction for the 666GHz band to the data (see Figure 10). While our best fit model appears to match the timing of the flares in the 666 GHz data quite well, the model tends to over predict flux in some areas when compared to our data. It is possible that the deviations between the best fit model and the data are dominated by the higher flux calibration uncertainty in this band, especially when considering such short timescales. On the other hand, our model (and the vdL model) are only capable of predicting emission at frequencies which are initially self-absorbed (i.e. optically thick). Thus the deviations between the best fit model and the data could also suggest that the emission we observe from the jet ejecta in the 666 GHz band is initially optically thin.
Appendix D Systematic Errors
As described in §4.3, we estimated additional uncertainties on our best fit parameters, to factor in how well our chosen model represents the data. Table 3 displays these uncertainties (+ for upper confidence interval, for lower confidence interval) for each fitted parameter.
Appendix E Two-Parameter Correlations
With the multi-dimensional posterior distribution output from our MCMC runs, we explored possible two-parameter correlations for our model. A significant correlation between a pair of parameters, that is common to all of the ejecta, could indicate a model degeneracy or a physical relationship between the two parameters. Out of the possible two-parameter pairs, we find interesting correlations involving the , , , and parameters. Figure 11 displays the correlation plots, along with the one-dimensional histograms of the parameters232323We make use of the corner python module to make these correlation plots; https://github.com/dfm/corner.py. The correlation between and (first column) indicates a known degeneracy in the vdL model. The correlation between and (second column) likely indicates a physical relationship between the parameters, where faster ejecta tend to have brighter fluxes. We find the same relationship when we look at the distribution of bulk speeds and fluxes across all the ejecta, and this relationship has been seen in other sources (see §5.1 for details). The final four correlations (columns 3 through 6) seem to indicate a degeneracy between all four parameters (or at least a sub-set of them), where different combinations of the four parameters could potentially produce similar flaring profiles.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Baars et al. (1986) Baars J., Altenhoff W., Hein H., Steppe H., 1986, Nature, 324, 39
- 2Ballard & Heavens (1992) Ballard K., Heavens A., 1992, MNRAS, 259, 89
- 3Barthelmy et al. (2015) Barthelmy S., D’Ai’ A., D’Avanzo P., Krimm H., Lien A. Y., Marshall F. E., Maselli A., Siegel M. H., 2015, GCN Circular, 17929, 1
- 4Beardmore et al. (2015) Beardmore A. P., Page K. L., Kuulkers E., 2015, ATEL, 8455, 1
- 5Bell (1978) Bell A., 1978, MNRAS, 182, 147
- 6Belloni (2010) Belloni T., 2010, in The Jet Paradigm - From Microquasars to Quasars, ed. Belloni, T.M., Lecture Notes in Physics 794, Springer
- 7Belloni et al. (1997) Belloni T., Mendez M., King A. R., van der Klis M., van Paradijs J., 1997, Ap J, 479, 145
- 8Bernardini et al. (2016 a) Bernardini F., Russell D. M., Shaw A. W., Lewis F., Charles P. A., Koljonen K. I. I., Lasota J., Casares J., 2016 a, Ap JL, 818, L 5