Observations of the Disk/Jet Coupling of MAXI J1820+070 during Its Descent to Quiescence

MAXI J1820+070 was monitored regularly by Swift since its initial discovery. In this work, we use a subset of the observations spanning the date range 2018 October 6 to 2019 April 17 (MJD 58397–58590; Target ID: 10627), covering the decline of the initial outburst of the source and the decline of the first reflare. Swift/XRT operated in windowed timing (WT) mode for all of the observations we analyze in this work. All photon counting (PC) mode observations had count rates that were either too high (such that photon pileup was too great to correct for) or too low (such that not enough photons were collected to construct a spectrum with good signal-to-noise ratio).

Data were reprocessed using the xrtpipeline tool, part of the HEAsoft v6.26.1 software suite for analysis of high-energy astrophysical data.¹⁵ Spectral products were extracted using the xrtproducts tool. Source counts were extracted from a circular region 20 pixels (≈47'') in radius, centered on the source. Background counts were extracted from an annulus centered on the source, with inner and outer radii 80 and 120 pixels, respectively. Response matrices were generated using version 20191017 of the calibration database (caldb).

Individual spectra were grouped such that each spectral bin contained a minimum of 15 counts. Spectral fits were performed with xspec v12.10.1f (Arnaud 1996), using Cash statistics modified for background-subtracted spectra (W-statistic; Cash 1979) as the fitting statistic, which is appropriate for fitting low-count spectra while also tending to χ² in the high-count regime. The majority of spectra were well fit with an absorbed power-law model (powerlaw in xspec), with interstellar absorption accounted for by the tbabs model (Wilms et al. 2000). The first three Swift spectra (MJD 58397–58402) were statistically improved with the addition of a disk blackbody model (diskbb) with an inner disk temperature in the range kT_in = 0.13–0.18 keV. The presence of the disk component was confirmed by a NICER observation on MJD 58400.9. When including a diskbb component for spectral fits at MJD > 58402, we found the disk normalization was consistent with zero. We extracted unabsorbed 1–10 keV fluxes using the cflux model. Uncertainties on the best-fit X-ray parameters are all 90% confidence unless otherwise stated.

The first 10 Swift observations of our campaign (MJD < 58420) were taken primarily through our GTO program, where we compensated for the decrease in flux with time by increasing the length of each subsequent Swift/XRT exposure. In turn, we obtained a relatively steady number of X-ray counts in each spectrum. Later in the decay when MAXI J1820+070 was even fainter, all Swift exposures were relatively short, making it difficult for us to assess degeneracies between best-fit column densities (N_H) and photon indices (Γ). Since a primary goal of our program is to quantify changes in X-ray spectral shape, and since we do not expect N_H to evolve during the decay, we chose to freeze N_H to 1.0 × 10²¹ cm⁻² throughout our entire campaign to allow a more uniform comparison between the higher- and lower-count spectra toward the beginning and end of our campaign, respectively. This value is (a) consistent with the Galactic column in the direction of MAXI J1820+070 (HI4PI Collaboration et al. 2016) and (b) consistent with the measured values of N_H over the course of its outburst (see, e.g., Shidatsu et al. 2018; Kajava et al. 2019; Xu et al. 2020). In addition, during the first 10 Swift observations, where we accounted for the diminishing count rate through increased exposure times, we find a weighted-average N_H = (1.1 ± 0.1) × 10²¹ cm⁻² if we allow N_H to be free.

We utilized Chandra observations of MAXI J1820+070 over the date range 2018 November 13 to 2019 June 11 (MJD 58435–58645), which includes observations during the decline of the initial outburst as well as the final stages of the decline of the first reflare. For all observations, the source was placed at the aim point of the S3 chip on the Advanced CCD Imaging Spectrometer (ACIS; Garmire et al. 2003). Three observations employed the High-Energy Transmission Grating (HETG; Canizares et al. 2005) as a filter in order to avoid photon pileup. Data were reduced using ciao (Chandra Interactive Analysis of Observations) v4.11 and the Chandra caldb v 4.8.4.1 (Fruscione et al. 2006). We reprocessed the data using the chandra_repro script to apply the latest calibrations and bad pixel files. Spectral products were extracted using the specextract script for the majority of the observations. For each of these, source counts were extracted from a circle of radius 2'' centered on the source, while we extracted the background from an annular region centered on the source with inner and outer radii 9'' and 20'', respectively.¹⁶ One observation (ObsID 20207) without the HETG in place suffered from extreme photon pileup, for which we were unable to effectively correct. In this case, we extracted a spectrum from the CCD readout streak that appears in bright observations such as this one.¹⁷ Source counts were extracted using two box regions of size 2'' × 51'' centered on each streak on either side of the piled-up point source, and background counts from two box regions of size 22'' × 51'' at the same location (with the source region excluded). Spectra were extracted with the ciao tool dmextract, and response files were created using mkacisrmf and mkarf. The effective exposure time of the readout streak spectrum was 296 s.

As with the Swift data, we used xspec to perform spectral fitting. For the shorter Chandra observations (exposure time ${t}_{\exp }\lt 10$ ks), we chose to bin the spectra such that each bin contained a minimum of one count. For the remaining, longer observations, we grouped the spectra to a minimum of 15 counts per spectral bin, as with the Swift spectral fits. All Chandra spectra were well fit with an absorbed power-law model and were not statistically improved with the addition of a diskbb model. Unabsorbed 1–10 keV fluxes were extracted in the same way as for the Swift data.

NICER observed MAXI J1820+070 regularly since its initial discovery. In this work, we utilize two observations from 2018 November 19 and 21 (MJD 58441 and 58443; ObsIDs 1200120310 and 1200120312; PI Gendreau) covering the initial decline immediately after the Swift and Chandra coverage ended. Both observations were recalibrated with the nicerl2 task in HEAsoft v6.26.1, using caldb version 20200202. Spectra were extracted from all active detectors, except detectors #14 and #34, which are prone to excessive noise. Background spectra were created using the "3C50_RGv5" model provided by the NICER team. We used response files recommended for caldb version 20200202.

We followed the same spectral fitting procedure as for the Swift data, grouping the data such that each bin contained a minimum of 25 counts and utilizing the Cash (1979) fitting statistic. Spectra were fit in the 0.5–10 keV range, and we included a multiplicative constant of 50/52 in the model to account for the excluded detectors. Both NICER spectra were well fit with an absorbed power-law model.

Our VLA campaign consisted of 19 observations obtained through NRAO programs 18A-277, SK0335, SJ0238, and SK0114; see Table A2. The first eight VLA observations (MJD 58398–58432) were taken in the most compact D configuration (maximum baseline ${B}_{\max }=1.03$ km), with the next five observations (MJD 58441–58517) in C configuration (${B}_{\max }=3.4$ km), and the final six observations (MJD 58603–58645; during the reflare decay) in B configuration (${B}_{\max }=11.1$ km). Observations lasted between ∼30 and 120 minutes (providing ∼1–80 minutes on source), with longer observations generally toward the end of each decay.

All data were taken in the C-band (4–8 GHz) with 4 GHz total bandwidth (2 × 2 GHz basebands centered at 5.0 and 7.0 GHz). For every observation, we used scans on 3C 286 for bandpass calibration and to set the flux density scale (using the Perley & Butler 2017 coefficients), and we interleaved our science observations with scans on the secondary calibrator J1824+1044 to solve for time-dependent complex gain solutions. Data were reduced using the Common Astronomy Software Application v5.6 (CASA; McMullin et al. 2007), and calibrations were performed using the VLA pipeline. A small amount of additional flagging was performed manually prior to imaging the data.

The data were imaged using tclean, using two Taylor terms to model spectral dependences of sources within the field. We used Briggs weighting to reduce sidelobes from other sources in the field, using robust values of 0.0, 0.5, and 1.0 when the array was in D, C, and B configurations, respectively. During the first nine observations, MAXI J1820+070 was >1 mJy, and we performed one to two rounds of phase-only self-calibration (down to 30–60 s solution intervals). We achieved rms noise levels ranging from ≈3 to ≈50 μJy bm⁻¹ across our entire campaign, as measured from source-free regions of our radio images. All flux densities were measured at 6 GHz.

Long-lived ejecta that were launched by MAXI J1820+070 during its soft state transition (see Section 2) were present during our coverage of the initial decay. During our D-configuration observations, we found it challenging to resolve the core from these ejecta in the image plane, which reached distances up to 13'' from the compact core. We therefore measured flux densities by fitting the blended core and ejecta in the uv plane using uvmultifit¹⁸ (Martí-Vidal et al. 2014), as described below.

For each D-configuration observation, we first produced an image in tclean, as described above. We then took the sky model produced by tclean, masked out the MAXI J1820+070 complex, and subtracted remaining "field" sources from the visibility set using the task uvsub. Depending on the date of observation and the time on source, we expected to detect anywhere from zero to two relativistic ejecta in each observation in addition to the compact core. We therefore ran three iterations of uvmultifit, requiring one, two, and three point sources. For the core, we left the radio spectral index α (f_ν ∝ ν^α, where f_ν is the flux density at frequency ν) as a free parameter, and for the relativistic ejecta we fixed the spectral index to α = −0.7 (Bright et al. 2020). We examined the resulting fit statistics to guide our decision on the best fit (i.e., the one-, two-, or three-component model); we also examined residual images after subtracting each model (i.e., we used uvsub to subtract each model from the visibility set and then created "dirty images" using tclean). The peak flux densities and spectral indices of only the core (from our preferred model) are reported in Table A2. Uncertainties incorporate both rms noise and errors related to the fitting process. Flux densities (and positions) for the relativistic ejecta were reported in Bright et al. (2020).

When the VLA was in its C configuration, we could resolve the compact core from the ejecta in the image plane. We generally found for these observations, when the core was fainter, that flux density measurements appeared more reliable in the image plane compared to the uv plane. We therefore measured the peak flux density of the compact core using the task imfit for our C-configuration observations, but using a model that required one, two, or three point-source components. As above, we examined residual images with each model subtracted to decide on the number of required components. The flux densities reported in Table A2 (for our C-configuration observations) are the peak flux densities reported by imfit only for the core. All of our B-configuration observations are during the reflare decay, by which time radio emission from the relativistic ejecta had completely faded.

For the C- and B-configuration observations, we measured radio spectral indices by splitting our bandwidth into 4 × 1 GHz basebands centered at 4.5, 5.5, 6.5, and 7.5 GHz. We then imaged each baseband separately and measured the flux density of the core in imfit (using multicomponent models in C configuration, as described above). The spectral index was then measured by performing a least-squares fit, with the uncertainty on α estimated by Monte Carlo simulations that randomized (and refit) each spectrum (see Plotkin et al. 2017b). Note that for the D-configuration observations, the spectral index was fit as a free parameter when running uvmultifit. Throughout the text, we adopt radio luminosities at 5 GHz as calculated from the measured radio flux density (at 6 GHz), the radio spectral index, and the source distance.

The X-ray light curve of MAXI J1820+070 is shown in the upper panel of Figure 1, highlighting the end of the initial outburst and the decay of the first reflare in 2019 March. We incorporate the uncertainty on the distance into our luminosity calculation. The decline of both the initial outburst and the first reflare cover a similar dynamic range in Figure 1, with each covering ∼3 orders of magnitude in L_X. It is unclear if MAXI J1820+070 reached its minimum quiescent L_X between the end of the initial outburst and the onset of the first reflare, because X-ray observations of the source at an L_X lower than that measured by our final Chandra epoch do not exist at the time of writing. However, considering the P_orb of MAXI J1820+070, we suspect that it may eventually settle to an L_X up to 1–2 orders magnitude lower than the minimum observed during our campaign (see, e.g., Figure 4 of Reynolds & Miller 2011).

Zoom In Zoom Out Reset image size

The lower panel of Figure 1 shows the evolution of the X-ray spectrum, parameterized by the photon index Γ. The hard X-ray spectrum (Γ ≲ 1.7), coupled with the presence of a compact jet, indicates that MAXI J1820+070 had transitioned to the hard state by the time we commenced our L_R–L_X monitoring program. This is confirmed by the detection of significant variability (fractional rms ≳35%) and quasi-periodic oscillations in the Swift and NICER light curves (Stiele & Kong 2020), typical of BH-LMXBs in the hard state.

We see a mild hardening of the X-ray spectrum of MAXI J1820+070 for the first ≈15 days of our program, with Γ evolving from ≈1.7 to 1.6 (reaching its hardest values as the source declined from L_X ≈ 10^36.5 to 10^35.6 erg s⁻¹, equivalent to ≈10^−2.5 to 10^−3.4 L_Edd). During the first three observations (MJD 58397–58402), we detect a cool, thin disk (kT = 0.13–0.18 keV) in the hard state that slowly fades and is not detectable at later epochs with lower-count spectra (L_X ≲ 10⁻³ L_Edd). An accretion disk component is not uncommon in the spectra of BH-LMXBs in the canonical hard state, particularly at these modest Eddington ratios (e.g., Miller et al. 2006; Reis et al. 2010; Reynolds & Miller 2013). During the decline of the second reflare of MAXI J1820+070, Xu et al. (2020) found strong evidence for a truncated disk at L_X ≈ 10^−2.6 L_Edd, similar to the luminosity at which we require a disk component in the spectra during the decline of the main outburst.

At lower X-ray luminosities during the initial decay, we observe Γ evolve from ∼1.6 to 2.0. The initial hardening and then softening of Γ with decreasing X-ray luminosity is a well-known trend (one suggestion is that it is driven by a change in the source of seed photons for inverse Comptonization; see, e.g., Sobolewska et al. 2011; Kajava et al. 2016). Based on two Chandra observations taken on MJD 58519, the X-ray spectrum did not appear to continue to soften indefinitely. Unfortunately, we cannot empirically confirm this statement, since we lack X-ray coverage for MJDs 58444–58519 because the source was too close to the Sun for X-ray observations. Nevertheless, even though we cannot determine exactly when the X-ray spectrum "saturated" to its maximum Γ during the initial decay, we can place a limit that Γ saturated on or after MJD 58444, when L_X ≈ 10³⁴ erg s⁻¹ (≈10⁻⁵ L_Edd).

Our coverage of the decay of the reflare started 1–2 days after its peak, in contrast with the observations of the decline of the initial outburst, which we commenced ∼200 days postpeak, after the source had undergone state transitions. During the decay of the reflare, we only observed the stage of the decay when Γ increases. Here, it is more apparent that once Γ ≈ 2 was reached, the X-ray spectrum plateaued to that value, which occurred sometime between the last Swift observation and the first Chandra observation during the reflare (58591 < MJD < 58603, 10⁻⁴ ≳ L_X/L_Edd ≳ 10⁻⁵). The X-ray spectrum then remains soft (Γ ≈ 2) over >1.5 decades in X-ray luminosity (10⁻⁵ ≳ L_X/L_Edd ≳ 10^−6.4).

4.1.1. Comparisons to Other Quiescent Sources

In Figure 2 we show Γ as a function of Eddington fraction for both the initial decay and reflare. This figure demonstrates the softening of the X-ray spectrum (i.e., increasing Γ) as X-ray luminosity decreases, with Γ eventually "saturating" to ≈2 at L_X ≈ 10⁻⁵ L_Edd. The bottom panel of Figure 2 shows the X-ray spectrum binned by X-ray luminosity, to ease comparisons to the ensemble average of other BH-LXMBs in the Chandra archive from Plotkin et al. (2013), combined with additional data from V404 Cygni (Plotkin et al. 2017b).¹⁹ Among all of the quiescent BH-LMXBs in the Chandra archive, the only individual source (besides MAXI J1820+070) where the X-ray softening has been tracked through the "plateau" stage is V404 Cygni (during its 2015 outburst decay; Plotkin et al. 2017b).²⁰ For V404 Cygni, the softening proceeded quickly (on a timescale of <3 days) and saturated at a slightly lower luminosity L_X ≈ 10^−5.5–10^−5.6 L_Edd.

Zoom In Zoom Out Reset image size

Comparing MAXI J1820+070, V404 Cygni, and the ensemble of other quiescent BH-LMXBs, we are starting to uncover potential variations in the X-ray softening between different sources. The general qualitative trend—that the softening completes below ≈10⁻⁵ L_Edd and then plateaus—holds. However, the softening appears to occur gradually in MAXI J1820+070, at least according to the initial decay where Γ increased over ∼30 days (we do not have sufficient coverage during the reflare decay; see Figure 1). Conversely, the active X-ray softening of V404 Cygni was much more rapid, occurring over <3 days (Plotkin et al. 2017b). Considering the strong relation between L_X and Γ, this softening timescale can be considered to be a proxy for the timescale of the luminosity decay.

Besides MAXI J1820+070 and V404 Cygni, we identify one other source in the literature with comparable coverage of the X-ray softening and plateau of Γ: 4U 1543−47, which softened from Γ ≈ 1.6 to 2.0 over ∼5 days and then plateaued during the decay of its 2002 outburst (Kalemci et al. 2005). Another BH-LMXB, Swift J1357.2−0933, also has excellent spectral constraints of the X-ray softening during its 2011 and 2017 decays, which took ∼90 days (Armas Padilla et al. 2013; Beri et al. 2019), but became too faint to obtain useful constraints during the plateau stage. Considering the ensemble of these four systems, BH transitions into quiescence appear to fall into categories with gradual softenings (i.e., MAXI J1820+070 and Swift J1357.2−0933) and fast softenings (i.e., V404 Cygni and 4U 1543−47). One might be tempted to link the softening timescale to the P_orb of the system (and therefore the disk size for Roche Lobe overflow systems), since Swift J1357.2−0933 with the slowest softening (∼90 days) has the shortest orbital period (P_orb = 2.8 hr; Corral-Santana et al. 2013) and V404 Cygni with the fastest softening (<3 days) has the longest orbital period (P_orb = 6.5 days; Casares et al. 1992), with MAXI J1820+070 (P_orb = 16.5 hr; Torres et al. 2019) and 4U 1543−47 (P_orb = 1.1 days; Orosz et al. 1998) falling in between. Nevertheless, the apparent differences in the softening timescales between these four systems further highlight the need for improved X-ray spectral coverage during decays of systems spanning a wide range of properties like P_orb, inclination, and donor mass, in order to attach a physical scenario to the timescale of the transition into quiescence.

In the upper panel of Figure 3, we show the 5 GHz radio light curve of the compact core of MAXI J1820+070, and we show the time evolution of the radio spectral index α in the lower panel. The radio spectrum is inverted in nearly all of our observations (typically α > 0.2; the weighted-average spectral index is $\bar{\alpha }=0.24\pm 0.06$, where the quoted uncertainty represents the standard deviation about the weighted average).

Zoom In Zoom Out Reset image size

We see some variations in α with time. To investigate a potential correlation between log₁₀ L_R and α, we calculate the Spearman correlation coefficient, ρ, and incorporate the uncertainties on the input parameters by simulating 10,000 data sets based on the original sample, allowing L_R and α to vary within their error bars according to a (log)normal distribution. We find ρ = −0.2 ± 0.2, where this value represents the median of the 10,000 values of ρ calculated from the simulated data, and the 1σ uncertainties are calculated as the 16th and 84th percentiles. The value of ρ is consistent with no correlation between α and log₁₀ L_R. To test the significance, we adopt a Monte Carlo method by randomly shuffling the data, creating 10,000 new log₁₀ L_R, α pairs and calculating ρ each time. We find that 17% of the time we measure an anticorrelation that is as strong as, or stronger than, ρ = −0.2, implying that our measured value of ρ is not significant.

Perhaps the most interesting variation in α occurs on MJD 58603, where α veers slightly negative (α = −0.11 ± 0.04), and possibly also on MJDs 58485 and 58517 (at much lower statistical significance). Such negative deviations, however, are not exceptional and do not approach becoming optically thin (where typically α ≈ −0.7), and similarly, mild negative deviations in α have been observed for other low-luminosity hard-state and quiescent BH-LMXBs (e.g., Rana et al. 2016; Espinasse & Fender 2018; Plotkin et al. 2019). A slightly negative value of α such as the instances discussed here could perhaps indicate the jet break temporarily moving below gigahertz radio frequencies. Alternatively, we may be seeing deviations from a simple conical jet, or we may be seeing an optically thin emission region (perhaps related to a small flare) superposed over a flat/inverted radio spectrum from the compact jet.

Figure 4 shows MAXI J1820+070 in the L_R–L_X plane. The parameters L_R and L_X, along with their corresponding spectral indices, are also tabulated in Table 1. MAXI J1820+070 displays a nonlinear correlation over ≈4 decades in X-ray luminosity extending from the hard state through quiescence. To measure the slope of the correlation, we adopt the Bayesian modeling package linmix,²¹ the python port of the linmix_err idl package (Kelly 2007), to measure the dependence of L_R on L_X. We assume a linear correlation of the form ℓ_R = b + mℓ_X, where ℓ_R and ℓ_X are the logarithms of L_R and L_X, respectively, and b and m are the ℓ_R-intercept and slope, respectively. Following Gallo et al. (2012, 2014, 2018), we include additional uncertainties of 0.3 dex on ℓ_R and ℓ_X to account for the lack of strict simultaneity between our radio and X-ray observations. MAXI J1820+070 follows a relation with slope m = 0.52 ± 0.07, with an intrinsic random scatter ${\sigma }_{0}={0.14}_{-0.07}^{+0.10}$ dex, where the best-fit values represent the median of 10,000 draws from the posterior distributions and the 1σ uncertainties are calculated as the 16th and 84th percentiles. The correlation slope is consistent with that derived for the larger hard-state BH population (m = 0.59 ± 0.02 from 36 BHs; Gallo et al. 2018).

Zoom In Zoom Out Reset image size

Table 1.  Summary of Results from the Quasi-simultaneous X-Ray and Radio Observations

MJDX	${\mathrm{log}}_{10}\,{L}_{{\rm{X}}}$(1–10 keV)a	Γ	MJDR	${\mathrm{log}}_{10}\,{L}_{{\rm{R}}}$(5 GHz)	α
	(erg s−1)			(erg s−1)
(1)	(2)	(3)	(4)	(5)	(6)
58397.09	36.5 ± 0.1	1.66 ± 0.02	58398.04	30.0 ± 0.1	0.23 ± 0.00
58400.01	36.2 ± 0.1	1.63 ± 0.03	58399.99	29.6 ± 0.1	0.27 ± 0.01
58402.28	36.0 ± 0.1	1.59 ± 0.04	58402.85	29.5 ± 0.1	0.42 ± 0.02
58404.26	35.9 ± 0.1	1.63 ± 0.02	58403.91	29.4 ± 0.1	0.51 ± 0.04
58406.73	35.8 ± 0.1	1.60 ± 0.02	58405.90	29.3 ± 0.1	0.42 ± 0.05
58410.30	35.6 ± 0.1	1.65 ± 0.02	58409.83	29.3 ± 0.1	0.23 ± 0.05
58419.61	35.3 ± 0.1	1.67 ± 0.01	58418.85	29.2 ± 0.1	0.28 ± 0.03
58432.55	34.7 ± 0.1	1.74 ± 0.14	58431.85	28.9 ± 0.1	0.48 ± 0.07
58441.11	34.2 ± 0.1	1.84 ± 0.07	58440.90	28.8 ± 0.1	0.29 ± 0.05
58518.61	32.9 ± 0.1	1.84 ± 0.10	58516.70	28.0 ± 0.2	−0.12 ± 0.17
58603.13	34.0 ± 0.1	1.97 ± 0.14	58603.32	28.8 ± 0.1	−0.11 ± 0.04
58609.33	33.2 ± 0.1	2.19 ± 0.23	58608.35	28.4 ± 0.1	0.30 ± 0.07
58618.67	32.7 ± 0.1	${1.99}_{-0.23}^{+0.24}$	58618.34	27.8 ± 0.1	0.55 ± 0.18
58627.42	32.8 ± 0.1	1.95 ± 0.15	58627.34	27.8 ± 0.1	0.44 ± 0.18
58636.67	32.6 ± 0.1	2.11 ± 0.12	58636.26	27.6 ± 0.1	0.56 ± 0.27
58645.08	32.6 ± 0.1	2.04 ± 0.07	58645.24	27.6 ± 0.1	1.06 ± 0.34

Note.

^aUnabsorbed.

Download table as:  ASCIITypeset image

MAXI J1820+070 is only the fourth individual source for which an unbroken L_R–L_X correlation has been tracked over four (or more) decades in X-ray luminosity, joining V404 Cygni (m = 0.54 ± 0.03; Corbel et al. 2008; Plotkin et al. 2017b), GX 339−4 (m = 0.62 ± 0.01; Corbel et al. 2013; see also Tremou et al. 2020), and XTE J1118+480 (m = 0.72 ± 0.09; Gallo et al. 2014). Note that we are aware of two other sources with L_R–L_X coverage over a comparable range of luminosities, H1743−322 and Swift J1753.5−0127; however, both sources appear to follow different correlation slopes above and below L_X ≈ 10³⁴ erg s⁻¹ (Coriat et al. 2011; Plotkin et al. 2017a), such that we do not compare MAXI J1820+070 to those two sources here. MAXI J1820+070 also shows the same L_R–L_X correlation over multiple decays, similar to V404 Cygni (Plotkin et al. 2017b) and GX 339−4 (Corbel et al. 2013). This suggests that the coupling between the disk and the jet is robust between outburst decays.

We can combine the slope of the L_R–L_X correlation with radio spectral information to draw some conclusions about the dominant emission mechanisms at work as the source decays. If the radio luminosity is responding to changes in the mass accretion rate ($\dot{M}$), then L_R ∝ ${\dot{M}}^{\tfrac{17}{12}-\tfrac{2}{3}\alpha }$ (Heinz & Sunyaev 2003; Markoff et al. 2003).²² The above assumes that jets are emitting synchrotron radiation from nonthermal particles accelerated into a power-law distribution, ${dn}/d\gamma \propto {\gamma }^{-p}$ (where n is the number density of particles, γ is the Lorentz factor, and p is the power-law index that describes the energy distribution of the particles; following Heinz & Sunyaev 2003, we adopt p = 2). For an arbitrary X-ray emission process, we can assume L_X ∝ ${\dot{M}}^{q}$, where q is a number that parameterizes the radiative efficiency of the X-ray process²³ (i.e., q ≠ 1 indicates a nonthermal or inefficient mechanism). Thus, the measured slope of the L_R–L_X correlation can be written as

$$\begin{eqnarray}&&m=\displaystyle \frac{\tfrac{17}{12}-\tfrac{2}{3}\alpha }{q}.\end{eqnarray} \tag{ 1 }$$

For m = 0.52 ± 0.07 and $\bar{\alpha }=0.24\pm 0.06$ (i.e., the weighted-average radio spectral index during our campaign), we find on average during the two decays that q = 2.4 ± 0.3. Thus, we require a radiatively inefficient source of X-ray emission.

Given our empirical constraint on the radiative efficiency (through the q parameter), we can make some inferences on the most likely source of X-ray emission. Inverse Comptonization or SSC processes in an RIAF can readily produce the observed X-ray signatures. As the accretion rate drops, one expects the density of the inner region of the accretion flow to decrease. In turn, the optical depth to inverse-Compton scattering may then decrease and yield a gradual increase in Γ (e.g., Esin et al. 1997).

In terms of jet origins for X-ray emission, we describe below that most scenarios are unlikely, except perhaps in some instances of jet-related SSC emission, which is difficult to distinguish from SSC originating from an RIAF. As expanded upon below, different mechanisms for (jet) X-ray emission are expected to emit with different radiative efficiencies. Thus, in the following, we walk through different scenarios for jet emission, and in order for a scenario to be deemed viable, its predictions must be consistent with our observations on the average radio spectral index, with the observed evolution in Γ, and with the behavior in L_R–L_X.

Since we have a precise measurement of the average radio spectral index, our combined constraints on the X-ray spectral softening and L_R–L_X allow us to exclude other scenarios for jet-dominated X-rays. However, the following discussion comes with the caveat that we assume that changes in radio/X-ray luminosities are driven primarily by changes in $\dot{M}$. If we assume that X-rays are always dominated by optically thin synchrotron radiation (again, emitted by a nonthermal electron population described by ${dn}/d\gamma \propto {\gamma }^{-p}$), then the X-ray spectral softening from Γ = 1.5 to 2 (as seen during the reflare) could be explained as optically thin synchrotron emission emanating from particles evolving to a steeper nonthermal distribution (p = 2Γ −1 evolving from 2 to 3). Then, according to Equation 17(c) of Heinz & Sunyaev (2003), we would expect the slope of the L_R–L_X correlation to evolve from 0.7 to 0.6 in coordination with the X-ray spectrum softening from Γ = 1.5 to 2.0 (again, adopting a radio spectral index of $\bar{\alpha }=0.24\pm 0.06$). Therefore, if p were truly evolving to a steeper distribution, then we would expect to see the L_R–L_X slope become shallower as the source decays to quiescence, an effect that does not appear to be seen in Figure 4.

From L_R–L_X we can also exclude the possibility that the X-ray spectral softening is caused by a synchrotron-cooled jet dominating the X-rays. If this were the case, we would expect the transition to occur at L_X ≈ 10³⁴ erg s⁻¹ (i.e., the luminosity at which Γ reaches 2), at which point the radiative efficiency would be ${q}_{\mathrm{cool}}=p+2-\left(\tfrac{3}{2}\right){\rm{\Gamma }}$ (Heinz 2004; see also Section 4.2 of Plotkin et al. 2017b). For p = 2 as expected for a distribution of synchrotron-emitting particles for which we measure Γ = 1.5 at its hardest (p = 2Γ −1; see also, e.g., Heinz & Sunyaev 2003), the radiative efficiency would be q_cool = 1. Adopting the weighted-average value of $\bar{\alpha }=0.24\pm 0.06$ and the weighted average of all measured values of Γ at luminosities L_X ≤ 10³⁴ erg s⁻¹ (Γ = 1.99 ± 0.09), we would therefore expect to see the slope of the L_R–L_X correlation increase to m = 1.23 ± 0.12 at L_X = 10³⁴ erg s⁻¹ (Yuan & Cui 2005). In Figure 4 we show that all of the lowest radio luminosity data points lie above the expected synchrotron cooling decay, so we can likely rule it out.

Though we have ruled out the majority of pure jet models for the X-ray emission, we still identify a jet origin for SSC emission as a possibility. Jet SSC emission could plausibly come from either a thermal or a nonthermal particle distribution. A softening might then be expected if the Comptonized spectrum is produced by single scatterings of synchrotron-cooled seed photons (see, e.g., Corbel et al. 2008; Plotkin et al. 2017b, for discussions). However, confirming or refuting this scenario requires numerical modeling that is outside the scope of this paper. This would allow us to understand how the extra source of cooling would affect particle energies and therefore the power and spectrum of the seed synchrotron photons, as well as the average number of scatterings before inverse-Comptonized photons escape, and so on. In lieu of such modeling, for the time being we believe that an RIAF origin for the X-ray emission is the most reasonable.

4.3.1. Deviations from L_R–L_X

We stress that our best-fit slope only describes the average path taken through the L_R–L_X plane, since we observe temporary deviations about our best-fit line (which our regression technique models as intrinsic scatter). Similar deviations have also been seen for GX 339–4, and they illustrate the importance of measuring L_R–L_X slopes over a wide dynamic range in luminosity (see Section 4.2.3 of Corbel et al. 2013 for a detailed discussion). In Figure 5 we plot the (logarithmic) radio luminosity residuals against α and find a marginal possibility of an anticorrelation. We calculate a Spearman correlation coefficient ρ = −0.5 ± 0.2, which we find has an ∼2% probability of occurring by chance in randomly shuffled data (see Section 4.2 for our methodology).

Zoom In Zoom Out Reset image size

Although we do not consider this to be a statistically significant anticorrelation, it is an intriguing result that warrants further observational scrutiny (on future outburst decays). Deviations from a simple Blandford & Königl (1979) compact jet can result in changes in α, which would then cause an anticorrelation between α and m (e.g., according to Equation (1), if α increases, then m decreases; also see Corbel et al. 2013 for a discussion on how α can influence a source's position in the L_R–L_X plane). Thus, if the anticorrelation is indeed real, it may not be unexpected, and it could be manifesting itself as the (observed) deviations from the global L_R–L_X correlation. Unfortunately, we cannot quantify if α really is the driver because we are not able to make precise-enough measurements of the slopes of individual deviations (since they cover too little dynamic range in luminosity) to directly compare changes in m to changes in α. Nevertheless, Figure 5 motivates the importance of obtaining meaningful constraints on α during high-cadence multiwavelength monitoring of future outburst decays.