Disk, Corona, Jet Connection in the Intermediate State of MAXI J1820+070 Revealed by NICER Spectral-timing Analysis

There are two populations of astrophysical black holes with strong observational evidence: stellar-mass black holes and supermassive black holes. Despite the orders-of-magnitude black hole mass difference, the two populations share many spectral and timing features, particularly in the X-ray band (Merloni et al. 2003; Falcke et al. 2004; Uttley et al. 2014). On the other hand, they also have distinct astrophysical significance. In the case of stellar-mass black holes, the distribution of spins and masses measured in black hole binary (BHB) systems observed in X-rays, and in those of mergers involving black holes detected via gravitational waves by LIGO/Virgo, will shed light on their formation channels (Reynolds 2020). The supermassive black holes in active galactic nuclei (AGNs) are believed to profoundly influence galaxy evolution via the radiation and outflows, which is referred to as AGN feedback (Fabian 2012).

BHBs exhibit a variety of accretion states with rich multi-wavelength characteristics in a single system, on a human timescale (Remillard & McClintock 2006). The outburst of BHBs can be described by a hysteresis pattern in the hardness-intensity diagram (HID; e.g., Belloni & Motta 2016, see also Figure 1). In a complete outburst, the BHB rises from the quiescent state when the accretion rate is very low, usually ≲10⁻⁵ of the Eddington limit, to the hard state when the X-ray emission is dominated by the non-thermal coronal emission. The luminosity keeps rising until it makes a transition into the soft state dominated by the thermal disk emission, and the transitional state is called the intermediate state. Then, with a steadily decreasing mass accretion rate, the BHB eventually transitions back to hard state and fades into quiescence. In the radio band, a compact and mildly relativistic jet is commonly observed in the hard state (e.g., Dhawan et al. 2000). During the hard-to-soft transition, the core jet is slowly quenched, and at the so-called "jet line," which loosely corresponds to the boundary of the hard and soft intermediate states, a highly relativistic, discrete, and ballistic jet is launched with unclear mechanism (Fender et al. 2004).

Zoom In Zoom Out Reset image size

The X-ray emission from a BHB system consists of the thermal emission from the accretion disk, the ubiquitous non-thermal continuum emission caused by Compton upscattering of disk photons intercepted by hot (of the order of ∼100 keV) electrons in the "corona," and a reflection component as the disk material reprocesses the illuminating coronal photons. The local reflection spectrum is composed of prominent features such as the Fe K emission complex (6–7 keV), the Fe K edge (∼7–10 keV), and the Compton hump (∼20–30 keV), and it is smeared by relativistic effects when produced close to the central object (Fabian et al. 1989). Within this paradigm, a type of time-lag feature called a reverberation lag is naturally anticipated because of the light travel time delay between the direct coronal and the reflected emission. The reverberation lag has been detected in both AGN and BHBs, in the forms of the soft excess emission and/or the Fe K emission line lagging behind the continuum-dominated band (e.g., Fabian et al. 2009; Kara et al. 2016; De Marco et al. 2017).

There are many remarkable questions waiting to be answered that are significant to understanding accretion and ejection physics in the strong gravity regime, such as the nature of the corona, and the state transition mechanism. Despite the ubiquity of the non-thermal X-ray emission, the fundamental properties of the corona are not well understood (Done et al. 2007). The difficulty lies in the fact that the innermost region of BHBs are too compact to be resolved via imaging, and it is challenging for X-ray spectroscopy alone to distinguish the underlying continuum models (e.g., single or double Comptonization regions; Dziełak et al. 2019) from which the geometry and formation of the corona could then be inferred. There are three different types of coronal models: radiatively inefficient accretion flow (RIAF) within the truncated accretion disk (advection-dominated accretion flow or other types of coronal flow; Ferreira et al. 2006; Narayan & McClintock 2008), static corona with different configurations (e.g., slab; Haardt et al. 1994), and the corona being the base of a jet (Markoff et al. 2001, 2005). One commonly adopted geometry is the lamppost geometry, which could very likely be due to the base of a jet, is also motivated by relativistic reflection itself in terms of the steep emissivity profile (Wilkins & Fabian 2012; Dauser et al. 2013).

The hypothesis that the corona is the base of a jet is motivated by observational evidence. The connection between the X-ray corona and the radio-emitting jet was first observed in the BHB GX 399−4 (Hannikainen et al. 1998; Corbel et al. 2000), and now is seen universally in the hard state of low-mass X-ray binaries (LMXBs; Gallo et al. 2003; Fender et al. 2004; Gallo et al. 2012). Strong correlations between optical/near-infrared (NIR) and X-ray observations also indicate a connection between the corona and the jet (Homan et al. 2005; Russell et al. 2006; Gandhi et al. 2008). This relation can even be extended to an AGN, in which it is called the "fundamental plane of black hole activity," and shows that there is a connection between the inner accretion flow (as probed by X-rays) and large-scale jets seen in the radio (Merloni et al. 2003; Falcke et al. 2004; Plotkin et al. 2012). Such observations led to models that suggest that the corona is in fact the base of a relativistic jet, where the hard X-rays are produced by inverse Compton scattering of thermal disk photons and synchrotron photons in the jet base (Fender et al. 1999; Markoff et al. 2005). Such a model, fitted to the broadband spectral energy distribution (SED) from radio to hard X-ray, suggests a corona/jet-base outflow (Markoff & Nowak 2004). Similar models have been extended to higher luminosities and larger masses to explain the reflection and reverberation-constrained coronal geometry in AGNs (King et al. 2017).

The state transition mechanism also remains a puzzle, with challenges arising from both theoretical and observational perspectives. On the theory side, it is computationally expensive to simulate long-term behaviors of accretion flows to track through accretion states; while in the observational arena the hard-to-soft transition happens on a fairly fast timescale of usually a week, so the data sets are valuable. It has been proposed that the inner hot flow collapses and transitions into a standard optically thick accretion disk during the hard-to-soft transition. However, this point of view is challenged by measurements of the inner disk radius, R_in, being within ∼ 10R_g (R_g = GM_BH/c²) in the bright hard state (above 1% of the Eddington limit; García et al. 2015; Wang-Ji et al. 2018), and is persistent at < 9 R_g during the transition (Sridhar et al. 2020).

MAXI J1820+070 is one of the brightest BHBs detected so far, and has exhibited a great amount of activity in radio (Bright et al. 2020; Homan et al. 2020), NIR (Sánchez-Sierras & Muñoz-Darias 2020), optical (Muñoz-Darias et al. 2019; Paice et al. 2019; Shidatsu et al. 2019; Veledina et al. 2019), ultraviolet (UV; Kajava et al. 2019), and X-ray (Buisson et al. 2019; Kara et al. 2019; Espinasse et al. 2020; Wang et al. 2020b; Buisson et al. 2021; You et al. 2021; Zdziarski et al. 2021). In the hard state, reverberation analysis using the Neutron Star Interior Composition Interior Explorer (NICER; Gendreau et al. 2016) showed that while the broad iron line profile remained constant (and only the narrow component decreased), the frequency range where the (soft) reverberation lag dominated increased (Kara et al. 2019). This result was interpreted as a contraction in the X-ray corona. This hypothesis has been later confirmed by reflection spectroscopy in the context of a dual-lamppost, where the height of the upper lamppost decreased while the lower one stayed constant (Buisson et al. 2019), and by the increase of the peak frequency of the quasi-periodic oscillations (QPOs; Ingram & Motta 2019). During the hard-to-soft transition, an extremely powerful superluminal ejection was observed in the radio band (Bright et al. 2020) very close in time with a small flare in the 7–12 keV band and a transition from type-C to type-B QPO (Homan et al. 2020).

While reflection spectroscopy and reverberation mapping fundamentally probe the same region around the black hole, until very recently, the spectral and timing modeling were done independently. Spectral-timing modeling, together with the vast improvement in high-throughput data, mean that we have the potential for studying the inner accretion flow with unprecedented detail. One model with this capability is reltrans (Ingram et al. 2019). In reltrans, a lamppost corona and razor-thin disk are assumed, where the lamppost is the simplest prescription of the corona as the base of the jet. There are two main kinds of time-lags: the low-frequency hard lags, which are fairly ubiquitous in both AGN and BHBs (e.g., Van der Klis et al. 1987; Nowak et al. 1999; Papadakis et al. 2001), and are described as fluctuations in the photon index Γ of the illuminating coronal emission; and the high-frequency soft/reverberation lags, which are caused by the light-crossing time delay between the continuum and reflected emission. Non-linear changes in the reflection spectrum due to the variable continuum spectrum are approximated via first-order Taylor expansion (Mastroserio et al. 2018). The model can be used to fit directly and simultaneously to both the real and imaginary parts of the cross-spectrum as a function of energy in different frequency bands, accounting correctly for the instrumental response, and self-consistently produces the time-averaged energy spectrum. It has led to successes in fitting to observational data of both AGN (Mrk 335; Mastroserio et al. 2020) and BHB systems (Cyg X-1; Mastroserio et al. 2019).

The NICER data set of MAXI J1820+070 is exceptional because it covers the bright state transition that happens on a short timescale. It was so bright that we can extend previous reverberation lag analysis into the state transition, which is usually difficult because of the decreased variability and limited signal-to-noise ratios of lag spectra.

In this Letter, we first qualitatively explore the evolution of iron line profile and the reverberation lag frequency and amplitude, and then fit the lag-energy spectra and time-averaged energy spectrum with the physically motivated reverberation model reltrans. This Letter is organized as follows. Section 2 describes the observations and data reduction, Section 3 provides the details of results, including a phenomenological first look (Section 3.1), the lag-energy spectral fitting (Section 3.2), and the time-averaged spectral fitting (Section 3.3). We present the discussion in Section 4, and summarize the results in Section 5.

MAXI J1820+070 is also known as ASASSN-18ey as it was discovered in optical by the All-sky Automated survey for supernovae (ASAS-SN; Shappee et al. 2014) on 2018 March 6 (Tucker et al. 2018), and was only discovered in X-ray by the Monitor of All-sky X-ray Image (MAXI; Matsuoka et al. 2009) 5 days later (Kawamuro et al. 2018). NICER monitored it on a daily cadence, from 2018 March 6 to 2018 November 21, generating an extraordinary data set covering the bright hard-to-soft transition that usually happens quickly (∼a week) and thus is difficult to catch. The NICER HID is shown in Figure 1, where Epochs 1–5 were chosen to cover the hard-to-soft transition are shown, along with Epoch 0 serving as a reference point in the bright hard state (the same as the first epoch investigated in Kara et al. 2019). In Epochs 1 and 2, we combine several observations that are close in time and hardness ratio, to boost the signal-to-noise ratios (see Table 1 for details in other epochs). We note that Epoch 5 (obs.ID 1200120197) corresponds to the X-ray flare seen in Homan et al. (2020), when type-C QPOs transitions into type-B, and a strong radio flare was observed ∼2 hr after that (see the "jet line" in Figure 1). To avoid the type-B QPO, we select the photons with arrival time before MJD 58,305.67740 (TIME < 142359302) in that observation.

The NICER data are processed with the data-analysis software NICERDAS version 2019-05-21_V006, and CALDB version xti20200722 with the energy scale (gain) version "optmv10." We use the standard filtering criteria: the pointing offset is less than 54'', the pointing direction is more than 40° away from the bright Earth limb, more than 30° away from the dark Earth limb, outside the South Atlantic Anomaly. We also filter out the events from the two noisy detectors, FPMs 34 and 14. In addition, we select events that are not flagged as "overshoot" or "undershoot" resets (EVENT_FLAGS = bxxxx00), or forced triggers (EVENT_FLAGS = bx1x000). A "trumpet" filter is also applied to filter out known background events (Bogdanov 2019). The cleaned events are barycenter corrected using the FTOOL barycorr. The spectra are then binned with a minimum count of 25 per channel, and the energy resolution was oversampled by a factor of 3. We do not include background in the analysis because of the very high count rate of the source, and the fitted energy range is 0.5–10 keV. We use the RMF version "rmf6s" and ARF version "consim135p," which are both a part of the CALDB xti20200722, and we combine 50 modules excluding FPMs 34 and 14. We note that some detectors were turned off starting from ObsID 1200120196 due to the high count rate, so the response matrices for those observations are combined differently and accordingly. For the purpose of timing analysis, the light curve segment length is 10 s with 0.001 s bins, which covers frequencies from 0.1 to 500 Hz.

All of the uncertainties quoted in this Letter are for a 90% confidence range, unless otherwise stated. All spectral fitting is done with XSPEC 12.11.1 (Arnaud 1996). In all of the fits, we use the wilm set of abundances (Wilms et al. 2000), and vern photoelectric cross sections (Verner et al. 1996). We used the χ² fit statistics.

3.1. A Phenomenological First Look

We first fit the time-averaged energy spectra in 3–10 keV with a model of (diskbb+powerlaw), and the data-to-model ratios are shown in Figure 2(a), in order to qualitatively asses the excess flux in the iron Kα region. A line-like feature exists in all epochs, which clearly evolves through the epochs. In the hard state (Epochs 0–1), it remains stable, while the strength of the narrow component decreases, as shown in Kara et al. (2019) and Buisson et al. (2019). During the hard-to-soft transition, the line feature remains very similar until Epoch 5, when the line gets blueshifted. When an extra Gaussian component is added to approximate the iron K emission, the line energy increases from 6.52 ± 0.05 keV (Epoch 1) to 6.92 ± 0.05 keV (Epoch 5), with tied line width. We explore models for the evolution of the reflection features in Section 3.2.

Zoom In Zoom Out Reset image size

From the timing perspective, we calculate the cross-spectrum in each epoch following standard Fourier timing techniques (Uttley et al. 2014), between the energy bands of 0.5–1 keV and 1.5–3 keV. We show in Figure 2(b) the lags as a function of frequency for Epoch 0 and 1, representing the beginning and end of the hard state, and also for the combined Epochs 2–5 due to the significantly reduced variability during the transition. All epochs show hard lags at low frequencies and soft lags at high frequencies, but the frequency range where the soft lag dominates is evolving. In the hard state, we confirm the results from Kara et al. (2019) that the soft lag frequency increases, along with a decreasing lag amplitude; while in the transition, we find the soft lag shows an opposite trend, when its frequency decreases and the amplitude increases. The behavior in the transition suggests a growing emitting region, possibly due to an expansion of the corona. We show the combined lag-frequency spectrum (Epochs 2–5) here for clarity; the individual epochs (see Figure A1 in the Appendix A) exhibit the trend in more detail that is also confirmed in the lag-energy spectral fits (see Section 3.3).

3.2. Time-averaged Flux Spectral Fitting

As shown in Section 3.1 and Figure 2, a reflection feature is present in the time-averaged energy spectra of all epochs. We first perform a simultaneous fit to Epochs 1–5 with the model (tbabs∗diskbb+reltransDCp) where reltransDCp ¹² is a flavor of reltrans, which adopts a more physically motivated thermal Comptonization model nthComp (Zdziarski et al. 1996) as the illuminating coronal emission, and can also probe high electron density (N_e) up to 10²⁰ cm⁻³ rather than a fixed value of 10¹⁵ cm⁻³. The high density model has been shown to contribute to soft excess emission below 2 keV (García et al. 2016) and may have the effect of decreasing the fitted iron abundance, which is usually found to be super-solar (García et al. 2018; Jiang et al. 2019).

The global parameters including the column density, inclination, and iron abundance are all tied among epochs, and a maximal spin is assumed (a_* = 0.998) to fit for R_in. We also add a 1% systematics to energies below 3 keV to reduce the effects from calibration uncertainties at low energies, and at the same time raise the statistical importance of the iron line band. The resulting χ²/d.o.f. is 2548.7/1112 = 2.29, with strong residual features at low energies <3 keV (see the lower panel of Figure B1 in the Appendix B for the contributions to χ² in Epoch 3 as a representative case). The most noticeable features include edge-like shapes near ∼0.5 keV and ∼2.4 keV, and Gaussian-like absorption around 1.6 keV. These remaining features are likely from NICER's calibration systematics because they could be attributed to oxygen (∼0.5 keV), aluminum K edge (∼1.56 keV), and gold M edge (2.1–4.5 keV complex), and thus we model them empirically.

Similar to Wang et al. (2020a), we adopt phenomenological models (edge and gabs) to model those features empirically. This leads to the final model in this work:

$$\begin{eqnarray*}&&\begin{array}{l}({\mathtt{tbabs}}\ast {\mathtt{diskbb}}+{\mathtt{reltransDCp}})\\ \ast {\mathtt{edge}}\ast {\mathtt{edge}}\ast {\mathtt{gabs}}.\end{array}\end{eqnarray*}$$

In the phenomenological models, only their energies are tied through epochs, allowing for possible attitude-dependent variations that are not accounted for in the current response matrix. The fit with the phenomenological calibration model has χ²/d. o. f. = 1349.2/1089 = 1.24, decreasing significantly by 1200 with 23 fewer d.o.f. compared to the model without calibration features (see the Appendix for more information on the phenomenological models). The contributions from model components, and the data-to-model ratios, are shown in Figure 3. The best-fit parameters are presented in Table 2.

Zoom In Zoom Out Reset image size

Table 2. Best-fit Parameters when Fitting the Time-averaged Spectra with the Model (tbabs∗diskbb+reltransDCp)∗edge∗edge∗gabs

Parameter	Epoch 1	Epoch 2	Epoch 3	Epoch 4	Epoch 5
NH	$({1.73}_{-0.07}^{+0.10})\times {10}^{21}$ cm−2
a*	0.998
i (degrees)	${81.0}_{-0.1}^{(p)}$
AFe	${1.2}_{-0.2}^{+0.1}$
Eedge,1 (keV)	${0.5000}_{(p)}^{+0.0002}$
Eedge,2 (keV)	${2.40}_{-0.01}^{+0.02}$
Egabs (keV)	${1.62}_{-0.02}^{+0.01}$
Tin (keV)	${0.318}_{-0.004}^{+0.005}$	${0.314}_{-0.005}^{+0.006}$	${0.311}_{-0.005}^{+0.006}$	0.406 ± 0.004	${0.601}_{-0.006}^{+0.004}$
Ndiskbb (104)	${3.0}_{-0.3}^{+0.2}$	4.7 ± 0.4	${9.6}_{-0.8}^{+0.7}$	3.5 ± 0.2	${1.12}_{-0.04}^{+0.06}$
Rin (RISCO)	<1.2	1.3 ± 0.2	<1.1	${1.6}_{-0.2}^{+0.1}$	${7.5}_{-0.9}^{+0.5}$
$\mathrm{log}{N}_{{\rm{e}}}$ (cm−3)	15.9 ± 0.1	${15.8}_{-0.3}^{+0.2}$	${16.8}_{-0.5}^{+0.3}$	16.2 ± 0.2	<15.5
Γ	1.84 ± 0.01	1.91 ± 0.02	${2.03}_{-0.02}^{+0.03}$	2.21 ± 0.01	${2.42}_{-0.01}^{+0.03}$
kTe (keV)	<74	<217	>31	>136	<34
h (Rg)	25 ± 2	25 ± 2	${22}_{-2}^{+4}$	20 ± 2	${9}_{-1}^{+2}$
$\mathrm{log}\xi$ (erg cm s−1)	${2.48}_{-0.17}^{+0.10}$	${2.47}_{-0.25}^{+0.23}$	${3.07}_{-0.80}^{+0.22}$	${3.48}_{-0.11}^{+0.09}$	${4.54}_{-0.14}^{+0.07}$
boost	1.8 ± 0.2	${2.1}_{-0.3}^{+0.2}$	${1.7}_{-0.2}^{+0.4}$	1.7 ± 0.1	3.0 ± 1.5
Nr (10−2)	${7.0}_{-0.2}^{+0.5}$	${8.8}_{-0.6}^{+1.9}$	${12.2}_{-0.6}^{+0.4}$	${11.0}_{-0.8}^{+0.4}$	${19.7}_{-1.8}^{+1.2}$
${\tau }_{\max ,1}$	$-{0.44}_{-0.06}^{+0.09}$	$-{0.37}_{-0.09}^{+0.06}$	−0.0 ± 0.08	${0.24}_{-0.09}^{+0.06}$	${0.56}_{-0.08}^{+0.07}$
${\tau }_{\max ,2}$	$-{0.055}_{-0.005}^{+0.006}$	−0.049 ± 0.007	−0.04 ± 0.01	$-{0.059}_{-0.007}^{+0.009}$	−0.019 ± 0.006
σ (keV)	>0.09	>0.09	>0.09	>0.09	${0.08}_{-0.03}^{(p)}$
Strength	0.012 ± 0.002	0.011 ± 0.002	0.011 ± 0.002	0.012 ± 0.002	0.004 ± 0.002
χ2/d.o.f.	1349.2/1089 = 1.24

Note. Errors are at 90% confidence level and statistical only. The inclination is limited to be <81° given by the binary inclination (Torres et al. 2020). For the phenomenological model gabs, the width is set to have an upper limit of 0.1 keV, and it is pegged at that upper limit in all epochs. Parameters without uncertainties were fixed at the quoted value. The index (p) means that the uncertainty is pegged at the bound allowed for the parameter.

Download table as:  ASCIITypeset image

In Figure 3, we see that this basic model is insufficient, and the shapes of residuals at the iron line are different among the epochs. At the end of hard state (Epoch 1) and the beginning of transition (Epoch 2), the residual is very narrow and peaks precisely at 6.4 keV, which has been seen previously with reflection spectroscopy in the hard state, and is usually interpreted as a distant reflector (e.g., García et al. 2015). Moreover, the narrow line has been seen in MAXI J1820+070 (Buisson et al. 2019). If an additional xillverDCp component is added in Epochs 1 and 2, with the xillverDCp parameters tied to the corresponding ones in reltransDCp except for the ionization parameter and the normalization, χ² decreases by 38 (Epoch 1) and 10 (Epoch 2) each with 2 fewer d.o.f., and the key parameter values shown in Table 2 stay similar.

However, when the source approaches the jet line, the residuals in Epochs 4 and 5 are more blueshifted (i.e., excess above 6.4 keV, instead of a narrow peak at 6.4 keV), as already indicated by Figure 2. This behavior makes the modeling more complicated. In the standard reflection spectroscopy regime, we find either an additional xillverDCp or a second reltransDCp could not fully describe the spectra, suggesting a different origin rather than the distant reflection. We preliminarily tested out two physically driven models, which turned out to both be promising: emission from the plunging region modeled by bbody (Fabian et al. 2020), and returning radiation approximated by relxillNS (Connors et al. 2020). Other possibilities include extra Comptonization of the reflection spectrum in the corona (Wilkins & Gallo 2015; Steiner et al. 2017), or a change in the disk thickness as the mass accretion rate rises.

The best-fit parameter values are presented in Table 2. The iron abundance is close to the solar value at ${1.24}_{-0.15}^{+0.10}$. The disk temperature increases from ∼0.3 keV to ∼0.6 keV, which is expected in a hard-to-soft transition because the soft state spectrum becomes dominated by the thermal disk emission. The inner edge of the disk remains < 2 R_ISCO in Epochs 1–4, but increases to ${7.5}_{-0.9}^{+0.5}\,{R}_{\mathrm{ISCO}}$ in Epoch 5. The boost parameter enables the model to account for deviations from the model assumptions of an isotropically radiating source and point-like geometry, and it is ∼2 in Epochs 1–4, and becomes 3.0 ± 1.5 in Epoch 5, which causes the continuum and reflection components to overlap in Figure 3(e). However, we note that even though the best-fit value is larger than previous epochs, it is consistent within 90% confidence. We notice several behaviors worth discussing: the inclination is pegged at the higher limit given by the binary inclination i < 81° (Torres et al. 2020), the ionization increases from ∼2.5 to ∼4.5. In fact, in another LMXB XTE J1550–564, the inclination was also measured to be high when the source progressed to softer states (Connors et al. 2019). The reason behind these behaviors is very likely to be the same: the fit tries to accommodate the blueshifted iron line that is prominent in Epochs 4–5, because a blueshifted line could result from either a higher ionization or a higher inclination angle.

3.3. Fitting the Lag-energy Spectra with the Reverberation Model

For each epoch, we fit simultaneously the lag-energy spectra in multiple frequency bands chosen to be outside of the central frequency ± half-width at half-maximum (HWHM) of QPOs and their (sub)harmonics fitted with the sum of several Lorentzians (Olive et al. 1998, see Figure 4 for an example). The QPO frequencies and the frequency bands in which we extract the lag-energy spectra in all epochs are shown in Figure C2 in the Appendix C. The reason to avoid QPO frequencies in the fit is that we observe harder lags at the QPO frequencies in the lag-frequency spectra compared to adjacent frequency bins, which could be due to a more complicated mechanism producing the lags (Ingram et al. 2016) rather than a combination of the pivoting power-law and the reverberation lag assumed in reltrans (see the Appendix for more details).

Zoom In Zoom Out Reset image size

We employ the exact same model, reltransDCp, used in the time-averaged spectral fit (Section 3.2). ¹³ As it is computationally expensive to fit all the epochs simultaneously, we need to fit epochs individually. For consistency, we fix some global parameters to reasonable values for all epochs. Using the proper motion and the distance to the system from radio parallax measurement, 2.96 ± 0.33 kpc, the inclination of the jet to the line of sight is estimated to be i = 63 ± 3° (Atri et al. 2020). Optical spectroscopy reveals a mass function of $f{(M)=({M}_{1}\sin i)}^{3}/{({M}_{1}+{M}_{2})}^{2}\,=\,5.18\pm 0.15{M}_{\odot }$, where M₁ and M₂ are the masses of the black hole and its companion (Torres et al. 2019); assuming the conservative range of mass ratio 0.03 < q ≡ M₂/M₁ < 0.4, the black hole mass is 10 ± 2 M_⊙. Therefore, we fix the inclination to i = 63°, and the mass to M = 12 M_⊙. The spin parameter is fixed at the maximal value of 0.998 to allow the inner disk radius to fit to any physically allowed value. The iron abundance is fixed at the solar value (A_Fe = 1). In the source direction, the Galactic column density is N_H = 1.3–1.5 in units of 10²¹ cm⁻² (Bekhti et al. 2016; Kalberla et al. 2005), so we fix it at 1.3 × 10²¹ cm⁻². In addition, we fix the coronal electron temperature kT_e = 100 keV because NICER's energy coverage is insufficient to constrain it, and the boost parameter at 1, corresponding to a static, isotropically emitting lamppost corona. Guided by the prior expectation that the inner edge of the accretion disk extends to the innermost stable circular orbit (ISCO) by the time the source enters the soft state (Remillard & McClintock 2006), and to simplify our model, we set the inner edge of the disk, R_in to have an upper bound of 10 R_ISCO. As it is difficult to constrain the photon index Γ with time-lag fits alone, we only allow Γ to vary between ±0.1 around the corresponding best-fit value in the time-averaged energy spectral fit (Section 3.2).

The resulting parameter values in all epochs are shown in Table 3. The power spectral density, the frequency ranges used to extract the lag-energy spectra and then fitted, the best-fitted reltransDCp model and lag-energy spectra for Epoch 1, are shown in Figure 4, as an example. We find that reltransDCp can describe the two types of time-lags very well: the low-frequency hard lag, and the high-frequency soft/reverberation lag. At low frequencies, the reverberation lag contributes a dip at the iron line energy (Mastroserio et al. 2018). The inner edge of the disk, R_in is consistent at ∼6 R_ISCO ∼ 7 R_g (when a_* = 0.998) where R_ISCO is the radius of the ISCO, whose radius is a function of the spin, decreasing from 6 R_g to 1.23 R_g when the spin increases from a_* = 0 to 0.998. With reflection spectroscopy in the hard state (R_in ∼ 5.3 R_g; Buisson et al. 2019) and continuum fitting method in the soft state (Zhao et al. 2020; Guan et al. 2020), MAXI J1820+070 is expected to be a slow-spinning black hole with a_* ∼ 0 − 0.2. Therefore, our measured disk truncation level is consistent with a disk extending very close to the ISCO radius during the state transition. We also note that the disk electron density tends to peg at the upper limit of the model, 10²⁰ cm⁻³, confirming the necessity for BHB modeling to go beyond the usually tabulated value.

Table 3. Best-fit Parameters when Fitting the Lag-energy Spectra in Several Frequency Ranges (Outside of the QPO Frequencies) with the Model reltransDCp

Parameter	Epoch 0	Epoch 1	Epoch 2	Epoch 3	Epoch 4	Epoch 5
i (degrees)	63
AFe	1
kTe (keV)	100
a*	0.998
M (M⊙)	12
NH (1021cm−2)	1.3
boost	1
h (Rg)	${42}_{-2}^{+1}$	${27}_{-3}^{+2}$	33 ± 5	${162}_{-45}^{+62}$	${114}_{-27}^{+34}$	>335
Rin (RISCO)	4.0 ± 0.4	6 ± 1	${6}_{-1}^{+2}$	>5	${6}_{-3}^{(p)}$	${5}_{-1}^{+2}$
Γ	${1.57}_{-0.02}^{+0.02}$	${1.802}_{-0.002}^{+0.007}$	${1.800}_{(p)}^{+0.004}$	${2.01}_{-0.03}^{+0.07}$	${2.26}_{-0.04}^{(p)}$	${2.41}_{-0.01}^{+0.05}$
$\mathrm{log}\xi$ (erg cm s−1)	${3.02}_{-0.06}^{+0.10}$	2.09 ± 0.06	2.22 ± 0.08	${1.40}_{-0.22}^{+0.39}$	<3.1	${2.26}_{-0.52}^{+0.20}$
$\mathrm{log}{N}_{{\rm{e}}}$ (cm−3)	${19.74}_{-0.08}^{+0.03}$	>19.99	>19.93	${19.46}_{-0.36}^{+0.20}$	>18.9	${18.59}_{-0.15}^{+0.28}$
χ2/d.o.f.	857.0/205	270.2/163	350.9/205	223.4/205	228.8/163	129.1/79

Note. Errors are at 90% confidence level and statistical only. Parameters without uncertainties were fixed at the quoted value. The index (p) means that the uncertainty is pegged at the bound allowed for the parameter. Guided by the prior that the inner edge of the accretion disk extends to the ISCO by the time the source enters the soft state (Remillard & McClintock 2006), and to simplify our model, we set the inner edge of the disk, R_in to have an upper bound of 10 R_ISCO. It is pegged at that limit in Epochs 3–4. The photon index, Γ, is set to vary between ±0.1 around the best-fit value in the time-averaged energy spectral fit in Section 3.2 (1.8, 1.9, 2.0, 2.2, 2.4 for Epochs 1–5; 1.6 for Epoch 0, which corresponds to Epoch 2 in Buisson et al. 2019), and is pegged at the lower and upper bounds in Epochs 2 and 4, respectively.

Download table as:  ASCIITypeset image

The key result is that the fitted coronal height decreases from ∼40 R_g at the beginning of bright hard state (Epoch 0) to ∼30 R_g at the end of hard state (Epochs 1–2), and then increases to ∼ 100–200R_g when the transition begins (Epochs 3–4), and eventually to ≳ 300R_g in Epoch 5 just before the jet line. This trend suggests that the corona is contracting in the hard state, and then is expanding in the hard-to-soft transition, which is consistent with the qualitative argument by tracking the evolution of reverberation lag frequency in Section 3.1. The increase in coronal height is strongly suggested by the lag spectra in our reverberation model because the statistics are much worse if the coronal height is fixed at the value found in Epoch 2 at the end of hard state (33R_g). The χ² increases by 11, 31, and 25 with 1 more d.o.f. in Epochs 3–5, respectively. The largest difference in fitted models lies in both low energies below ∼0.8 keV, and also energies above ∼3 keV (e.g., see Figure 5). We also want to emphasize that even though some parameters are fixed, the trend in the coronal height persists in all the fitting trials with different setups, which is reasonable because within the reverberation model, the reverberation lag frequency fundamentally correlates with the size of the corona-disk emitting region. We note that the fit statistics are worse in Epochs 0–2 compared to Epochs 3–5, likely because the signal-to-noise ratios of lags are higher in the early epochs, while the variability decreases rapidly when the source approaches the jet line/soft intermediate state. There are common residuals at low energies <1 keV in all epochs, which is similar to what is shown in Figure 4, and could suggest some other contributions to the lags (e.g., an additional thermalization time-lag; see the discussion in Section 4).

Zoom In Zoom Out Reset image size

3.4. Comparison Between Time-lags and Flux Spectra, and Lessons Learned

With the extraordinary NICER data set of MAXI J1820+070 in the bright hard-to-soft transition, we are able to find that the reverberation lag frequency decreases in the transition, with an increasing amplitude, which indicates a larger emitting region that causes reverberation, likely due to an expanding corona. With the reverberation model reltransDCp, we are able to quantitatively measure the coronal height, by fitting both the lag-energy spectrum and the time-averaged energy spectrum. In the time-lag fits, the X-ray coronal height increases from several tens to several hundred R_g several days before the radio flare, consistent with the qualitative trend seen in reverberation lags. An expanding or ejected corona during a state transition may explain earlier Rossi X-Ray Timing Explorer observations of the high-mass X-ray binary Cyg X-1, where the low-frequency hard lags were seen to increase during the transition (Pottschmidt et al. 2000; Grinberg et al. 2013). We want to stress that to our knowledge, this is the first BHB in which the increase of reverberation/soft lags is clearly observed.

Next, we will discuss this result in the context of the radio jet. Previous work suggested that the corona contracts during the hard state (Kara et al. 2019), when the steady jet becomes weaker. This may suggest that the jet luminosity is related to the vertical extent of the X-ray corona, and would favor the hypothesis that the corona is the base of a jet. Now, we have extended this methodology to the hard-to-soft state transition. Our measured X-ray coronal heights from fitting the time-lags, and the radio flux density from the Arcminute Microkelvin Imager (AMI; Zwart et al. 2008; Hickish et al. 2018) during the transition are presented in Figure 6(a). The radio flux density indicates a weakened steady jet, following the trend in the hard state, followed by a powerful transient jet in the transition that features ballistic ejections (Bright et al. 2020), before the source transitions to the soft state.

Zoom In Zoom Out Reset image size

We propose that the increase in coronal height inferred from our reverberation analysis coincides with the launch of the blob that ∼5 days later causes a radio flare when it reaches a large enough scale to become radio bright (Figures 6(b), (c)). We note that due to the low-cadence radio monitoring between MJD 58302 and 58304, we cannot rule out that there was another earlier flare that was missed, and thus the time delay between corona expansion and radio flare would be less than 5 days. If the blob is launched when the corona is most compact (MJD ∼58,298) and travels at ∼0.9 c (Bright et al. 2020), then a delay of ∼5 days until a radio flare would suggest that the radio-emitting blob is at 10⁹ R_g. Unfortunately, there have been no published radio measurements of the synchrotron break to independently estimate the radius and distance of the jet ejecta, but this distance is reasonable, when compared to the BHB MAXI J1535–571 where the emission height above the black hole is ∼ 10⁸ R_g (Russell et al. 2020). We also note that recent Chandra imaging showed a transient X-ray jet ejecta at the time of the radio flare, inferring that the ejecta was a cone of 1.5 × 10⁴ au by 7.7 × 10³ au (Espinasse et al. 2020). As 1 au ∼ 10⁷ R_g for a 10 M_⊙ black hole, this would lead to a vertical extent of the transient jet emitting in X-ray of ∼10¹¹ R_g, and the light travel distance of several days is always below this upper limit. Therefore, it is plausible that the expansion of X-ray corona is accompanied with a jet knot that propagates downstream (moves higher up) with the jet material when approaching the time when the radio flare was observed. This interpretation is also suggested to explain the increased hard lag amplitude during the transition in the high-mass X-ray binary Cyg X-1 (Pottschmidt et al. 2000). It is consistent the physical picture presented in Russell et al. (2019), where the soft X-ray flux is observed to increase just before the radio flare in another LMXB MAXI J1535–571, because the X-ray flare could result from inverse Compton scattering of the synchrotron radiation from the knot after it was ejected (Russell et al. 2019). Finally, this behavior may be analogous to observations in some high-luminosity AGNs, where X-ray flares and iron line narrowing precede radio ejection events (Marscher et al. 2008; Lohfink et al. 2013; King et al. 2016). In this scenario, the source of the hard power-law emission is more likely to arise from synchrotron emission in addition to thermal Comptonization because the optical depth along the jet axis (τ ∝ N_j /r where N_j is the power injected in the jet, and r is the radius of the jet) drops quickly (Lucchini et al. 2021), while efficient thermal Comptonization requires a mildly optically thin medium (τ ∼ 0.1–1).

We also notice that the corona is found to be physically and radiatively compact within several tens of R_g from microlensing in lensed quasars (Chartas et al. 2016), X-ray reverberation mapping in Seyfert galaxies (Kara et al. 2016), and X-ray eclipses of the corona (e.g., Gallo et al. 2021). Therefore, the exact values of the coronal height (especially the ones at hundreds R_g in time-lag fits) need to be taken with caveats, such as possible contribution from thermalization time in the disk (García et al. 2013) which we will implement in reltrans the near future, or even the effects from returning radiation that are shown to increase the reverberation lags by ∼50% (Wilkins et al. 2020). With all this in mind, we stress that within the reverberation model, the relative trend of increasing coronal height during the state transition comes fundamentally from the change in the reverberation lag frequency.

We have analyzed five epochs of MAXI J1820+070 in the hard-to-soft transition using NICER, to study the spectral-timing evolution of the system. Our major findings are as follows.

• 1.  
  During the hard-to-soft transition, the frequency range where the reverberation lag dominates decreases, and reverberation lag amplitude increases, suggesting a growing X-ray emitting region, possibly due to an expanding corona.
• 2.  
  From the flux spectral fit, the inner edge of disk does not decrease. This result challenges the state transition model of a collapsed hot inner flow.
• 3.  
  By fitting the lag-energy spectra in a set of Fourier frequency ranges with the reverberation model reltransDCp, we measure the coronal height to be anomalously large. The height of the corona increases ∼5 days before a radio jet ejection event.
• 4.  
  We propose that the correlation between the X-ray corona height and the radio jet behavior can be explained with a model where the X-ray corona is the base of the jet thatis ejejcted in its final moments before transitioning to the soft state.

While this work was being finalized, we became aware of De Marco et al. (2021), showing the decrease in the observed soft lag frequency during the state transition, corroborating this aspect of our work. J.W., G.M., E.K., and J.A.G. acknowledge support from NASA grant 80NSSC17K0515, and thank the International Space Science Institute (ISSI) and participants of the ISSI Workshop "Sombreros and lampposts: The Geometry of Accretion onto Black Holes" for fruitful discussions. J.A.G. thanks the Alexander von Humboldt Foundation. A.I. and D.B. acknowledge support from the Royal Society.

Due to the reduced variability during the transition, in Figure 2, the decrease in the frequency range where the soft lag dominates during the transition is shown by comparing lag-frequency spectra of Epoch 1 and combined Epochs 2–5. The individual epochs exhibit the same decrease of soft lag frequency compared to the end of hard state, as shown in Figure A1.

Zoom In Zoom Out Reset image size

In the flux-energy spectral fit (Section 3.2), we see strong residual features at low energies $\lt 3\,\mathrm{keV}$. Using Epoch 3 as a representative case, the lower panel of Figure B1 shows the contributions to ${\chi }^{2}$. We adopt phenomenological models (edge and gabs) to model those features empirically. The effective area of NICER with the vertical lines indicating the fitted characteristic energies of phenomenological models, in this work and in a previous NICER spectral analysis of another BHB GX 339–4 (Wang et al. 2020a), are shown in Figure B1 (top panel). The features in MAXI J1820+070 could be attributed to oxygen (∼0.5 keV), aluminum K edge (∼1.56 keV), and gold M edge (2.1–4.5 keV complex). The characteristic energies are not identical in the two cases, and that could be due to the updated gain solution and response matrices; also, in previous work, a Crab correction was applied. We suggest that readers be aware of such residuals and follow the most recent NICER analysis guidelines.

Zoom In Zoom Out Reset image size

For the epochs during the transition, we find "harder" lags at the frequencies of type-C QPOs and their (sub)harmonics than at neighboring frequencies (Figure C1). This suggests a very likely different, or at least a more complicated, mechanism to produce the time-lags involving geometric effects such as Lense-Thirring procession (Stella & Vietri 1997), rather than a combination of the pivoting power law and reverberation as considered in reltrans. Therefore, we extract the lag-energy spectra and fit those with reltrans in the frequency bands avoiding the QPO frequencies (see Figure C2).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The phase lags at the type-C QPO frequencies have been systematically investigated in Van den Eijnden et al. (2017) and Zhang et al. (2020), but a similar feature has not been reported. Perhaps this is because those observations were mostly of type-C QPOs in the hard state, while we find that such hard lags are specific to type-C QPOs observed during the relatively short-lived transition and are not seen as clearly in the hard state.