{'fraction_non_alphanumeric': 0.03991484199383913, 'fraction_numerical': 0.038428144958907694, 'mean_word_length': 4.3326568148664935, 'pattern_counts': {'":': 1, '<': 10, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The NICER"Reverberation Machine": A Systematic Study of Time Lags in Black Hole X-Ray Binaries

We perform the first systematic search of all NICER archival observations of black hole (and candidate) low-mass X-ray binaries for signatures of reverberation. Reverberation lags result from the light travel time difference between the direct coronal emission and the reflected disk component, and therefore their properties are a useful probe of the disk-corona geometry. We detect new signatures of reverberation lags in 8 sources, increasing the total sample from 3 to 11, and study the evolution of reverberation lag properties as the sources evolve in outbursts. We find that in all of the 9 sources with more than 1 reverberation lag detection, the reverberation lags become longer and dominate at lower Fourier frequencies during the hard-to-soft state transition. This result shows that the evolution in reverberation lags is a global property of the state transitions of black hole low-mass X-ray binaries, which is valuable in constraining models of such state transitions. The reverberation lag evolution suggests that the corona is the base of a jet which vertically expands and/or gets ejected during state transition. We also discover that in the hard state, the reverberation lags get shorter, just as the QPOs move to higher frequencies, but then in the state transition, while the QPOs continue to higher frequencies, the lags get longer. We discuss implications for the coronal geometry and physical models of QPOs in light of this new finding.


INTRODUCTION
Black holes are the most extreme laboratories to study accretion and ejection physics, and to eventually test theories of gravity. In addition to the supermassive black holes at the centers of galaxies, there are also stellar-mass black holes, which are mostly discovered when they are X-ray bright in systems called X-ray binaries (XRBs). These are binary systems comprised of a compact object accreting material from a companion star (see Remillard & McClintock 2006 for a review and references therein). If the mass of the companion star is lower than roughly one solar mass, it is called a lowmass XRB (LMXB; e.g., Bradt & McClintock 1983). Depending on the nature of the central compact object, the source is categorized as either a black hole XRB (BHXB) or a neutron star XRB (NSXB).
Typically, LMXBs are transient sources which go through cycles of activity called outbursts Remillard & McClintock 2006;Tetarenko et al. 2016; and see Fig. 1). During the full outburst of a given black hole low-mass X-ray binary (BHLMXB), the system starts from the quiescent state with a very low mass accretion rate, and increases in luminosity through the "hard state", when X-ray emission is dominated by a hard X-ray component that can be approximately described by a power-law with photon index Γ in the range of 1.4-1.8 and is usually attributed to the non-thermal coronal emission. After several weeks, the spectrum can change dramatically, and as the BHXB transitions, it progresses through the "hard intermediate state" (HIMS), to the "soft intermediate state" (SIMS), before switching to the "soft state" dominated by the thermal emission of the accretion disk (see Fig. 1). The fractional rms can reach 30%-40% in the hard state, but drops to only a few percent in the soft state. After spending weeks to months in the soft state, the luminosity steadily drops, and the source makes a transition back to the hard state (at a lower luminosity than the hard-to-soft transition) and then eventually decays to quiescence. There is a remarkable coupling between accretion and ejection in BHLMXBs, as the properties of the radio jets strongly depend on the accretion states (Fender et al. 2004(Fender et al. , 2009). In the hard state, a compact and mildly relativistic jet is almost always observed (e.g. Dhawan et al. 2000); and there is a correlation between the radio and X-ray luminosities (Hannikainen et al. 1998;Corbel et al. 2000;Gallo et al. 2003). During the hard-to-soft transition, the radio luminosity decreases, suggesting that the compact jet is slowly quenched, and at the so-called "jet line" (which loosely corresponds to the boundary between the HIMS and SIMS), a highly relativistic, discrete and ballistic jet is launched . During the soft-to-hard transition and the decay in the hard state, when the compact jet is reformed is not well understood (e.g., Miller-Jones et al. 2012;Kalemci et al. 2013).
Though we have observed the ubiquitous non-thermal X-ray emission from XRBs for several decades, our understanding of its nature is still limited (Done et al. 2007;Nowak et al. 2011). The widely accepted idea is that the X-ray corona, a hot plasma located close to the black hole with electron temperature of the order of ∼ 100 keV and optical depth τ ≈ 1, produces the power-law-like X-rays via thermal Comptonisation of the disk photons. However, how the X-ray corona is generated in the first place, and how it evolves throughout an outburst, are not well understood. The most frequently used coronal models include: (1) an inner hot flow (e.g., advectiondominated accretion flow; Esin et al. 1997;Narayan & McClintock 2008) interior to a truncated thin accretion disk (Mahmoud et al. 2019;Zdziarski et al. 2021), and (2) the corona being the base of a jet (Markoff et al. 2001(Markoff et al. , 2005 or an aborted jet (Ghisellini et al. 2004), the simplest realization of which is the lamppost geometry (Matt et al. 1991(Matt et al. , 1992Dauser et al. 2013).
It then becomes advantageous to study the spectraltiming evolution of XRBs to be able to probe the inner disk-corona region. With a technique called reflection spectroscopy, we model the time-averaged flux-energy spectrum of the disk material that is irradiated by the coronal photons. The most prominent reflection features include the Fe K emission complex (6 − 7 keV), the Fe K-edge (∼ 7 − 10 keV), and the Compton hump (∼ 20 − 30 keV), all broadened by relativistic effects (Fabian et al. 1989;Matt et al. 1993;García et al. 2014). However, in modeling the flux-energy spectrum alone, limitations arise from degeneracies such as the model choice for the coronal emission (Dzie lak et al. 2019), and the positive correlation between the iron abundance and the BH spin parameter (García et al. 2015). These prevent us from unequivocally answering important questions, such as how truncated is the disk in the bright hard state (e.g., Wang et al. 2020b) or how quickly is the BH spinning (Dauser et al. 2013).
X-ray timing analyses can help break some of these degeneracies and build a more complete picture by providing insights about the nature of variability. Fourier timing techniques allow us to separate variability on different timescales that are believed to result from distinct physical processes (Uttley et al. 2014). Besides its unique signature in the flux-energy spectrum, the reflection process can also be seen in the form of a time lag feature, which is commonly referred to as reverberation. It arises because of the light travel time delay between the direct coronal emission and the reflected emission. Reverberation lags have been detected in both AGN and BHXBs at relatively high Fourier frequencies, in the forms of variations in the soft-excess emission and/or the Fe K emission line lagging behind those in the energy band dominated by the coronal emission (e.g., Fabian et al. 2009;Kara et al. 2016;De Marco et al. 2017). At relatively low Fourier frequencies (when probing longer timescale variability), there is a ubiquitous hard lag whose amplitude is too large to result from light travel time difference; the most popular model is mass accretion rate fluctuations propagating through the disk, moving from the colder to the hotter regions, causing soft photons to respond before hard photons (Lyubarskii 1997).
In a typical lag-frequency spectrum, the low-frequency hard lag (positive by convention) has a power-law-like shape and decreases with frequency, making it possible to see the soft lags (negative by convention) caused by reverberation at high frequencies. The soft lag flips into positive above some frequency ∼ 1/(2τ 0 ) due to phasewrapping, where τ 0 is the intrinsic time delay between the coronal and reflected emissions. At even higher fre-quencies, the lag oscillates around zero, again as a result of phase-wrapping (see e.g., Fig. 7 in Emmanoulopoulos et al. 2014, and also Fig. 3e). The observed lag amplitude depends both on τ 0 and the level of "dilution". The lag-frequency spectrum is measured between two energy bands, each dominated by the coronal or reflected emission. However, in reality, both energy bands contain contributions from coronal and reflected emission. This dilutes the time delay causing the observed lag amplitude to be smaller than the intrinsic lag (see e.g., Kara et al. 2013;Uttley et al. 2014;De Marco et al. 2015).
In addition to the time lags, another observable that may encode information about the inner disk-corona geometry is quasi periodic oscillations (QPOs; see the review Ingram & Motta 2019, and references therein). Type-C QPOs are strong and narrow peaks in the power spectrum. They are prominent in the bright hard state and HIMS, with centroid frequency rising from a few mHz to ∼ 10 Hz as the sources evolve. Type-B QPOs are characterized by a relatively strong and narrow peak with centroid frequencies of 5-6 Hz. The most rare Type-A QPOs are weak and broad in the power spectrum. Type-B and Type-A QPOs appear in the SIMS . Therefore, a common definition of the SIMS is the presence of the Type-B and/or Type-A QPOs. One promising interpretation of the Type-C QPOs invokes the relativistic precession model, in which the fundamental frequency is simply assumed to be the Lense-Thirring precession frequency at some characteristic radius (e.g., the disk truncation radius within which the inner hot flow is precessing; Ingram et al. 2009). As the Type-B QPOs are observed to be closely related to the launch of discrete jet ejections (Homan et al. 2020), they are likely to have a jet-based origin, e.g., the jet is launched from a precessing flow and precesses with it (Motta et al. 2015;Stevens & Uttley 2016;de Ruiter et al. 2019).
A systematic study of reverberation lags has been performed for both Seyfert galaxies (De Marco et al. 2013;Kara et al. 2016) and BHXBs (De Marco et al. 2015). Intrinsically, as the number of counts within the variability timescale over which we detect reverberation lags is much lower in BHXBs than in Seyfert galaxies, reverberation lag detection is more challenging for BHXBs. With the XMM-Newton archival data of 10 BHXBs, soft reverberation lags were detected in two sources in the hard state: GX 339-4 andH1743-322 (De Marco et al. 2015).
Now with the Neutron Star Interior Composition Interior Explorer (NICER; Gendreau et al. 2016) we have the opportunity to systematically search for reverberation signatures in BHXBs, and to globally study their evolution. NICER is an ideal instrument for spectraltiming analysis because of its large effective area in the 0.3-12 keV range, superior timing capabilities (i.e., high resolution and very little dead-time), with negligible pileup and low background. In addition, NICER has fast slew capabilities, which allows us to monitor BHXB systems at high cadence, throughout their outbursts. NICER has observed several outbursts, some from newly discovered black holes, covering all spectral states. Therefore, with the rich NICER archive, we can perform a systematic search for reverberation lags.
Thus far, reverberation lags have been reported in two BHXBs with NICER. In the hard state, reverberation lags are detected in MAXI J1820+070 (Kara et al. 2019) and in GX 339-4 (during its 2017 and 2019 outbursts; Wang et al. 2020a). In Wang et al. (2021a), we were able to trace reverberation lags into the HIMS for the first time, thanks to the unprecedented data quality of MAXI J1820+070. Our main finding was that the Fourier frequency range of the reverberation signature decreases during the hard-to-soft state transition with an increasing lag amplitude, possibly due to a vertically expanding corona. We successfully tested this hypothesis with the publicly-available reverberation model reltrans Mastroserio et al. 2021), whereby the longer time lags during the state transition are explained with an increasing height of the lamppost corona over the accretion disk. This X-ray coronal expansion preceded a radio flare by ∼ 5 days, hinting that the corona is the base of the jet, and the hard-to-soft transition is marked by the corona expanding vertically and launching a jet knot that propagates along the jet stream at relativistic velocities. These X-ray reverberation results were confirmed independently by De Marco et al. (2021).
Motivated by these reverberation results on individual sources, in this paper we perform the first systematic search for reverberation lags in BH and BH candidate (BHC) LMXBs observed with NICER. We perform a detailed study on their evolution throughout the hard and intermediate states, to learn whether the behavior seen in MAXI J1820+070 is common to other sources. The lag properties we focus on are which Fourier frequency range is dominated by reverberation and what the reverberation lag amplitudes are. This paper is organized as follows. Section 2 describes the observations (sample and data reduction) and our method to efficiently search for reverberation lags. Section 3 focuses on the search results and the evolution of reverberation lags. We present the discussion in Section 4, and summarize the results in Section 5.

The sample
The systems of interest in our study are the BH and BHC LMXBs. We cross-matched the sample in Tetarenko et al. 2016 (see Tables 8 and 12 therein) and included newly discovered systems since then (private communications with John Tomsick). In total, since its launch in 2017 until June 2021, NICER has observed 26 BH/BHC LMXBs (see Table 1 and A1 for the full list).

Data reduction
The NICER data are processed with the data-analysis software NICERDAS version v2020-04-23 V007a, and energy scale (gain) release "optmv10". We use the following filtering criteria: the pointing offset is less than 60 , the pointing direction is more than 30 • away from the bright Earth limb, and more than 15 • away from the dark Earth limb, and the spacecraft is outside the South Atlantic Anomaly (SAA). Data are required to be collected at either a sun-angle > 60 • or else collected in shadow (as indicated by the "sunshine" flag). We filter out commonly-noisy detectors FPMs #14, 34, and 54. In addition, we flag any "hot detectors" in which X-ray or undershoot (detector resets triggered by accumulated charge) rates are far out of line with the others (∼ 10σ) and exclude those detectors for the GTI in question. We select events that are not flagged as "overshoot" (typically caused by charged particle passing through the detector and depositing energy) or "undershoot" resets (EVENT FLAGS=bxxxx00), or forced triggers (EVENT FLAGS=bx1x000), and require an event trigger on the slow chain which is optimized for measuring the energy of the event (i.e., excluding fast-chain-only events where the fast chain is optimized for more precise timing). A "trumpet" filter on the "PI-ratio" is also applied to remove particle events from the detector periphery (Bogdanov 2019). The resulting cleaned events are barycenter corrected using the FTOOL barycorr. The background spectrum is estimated using the 3C50 background model (Remillard et al. 2021). GTIs with overshoot rate > 2 FPM −1 s −1 are excluded to avoid unreliable background estimation. We use the RMF version "rmf6s" and ARF version "consim135p", which are both a part of the CALDB xti20200722.

The "Reverberation machine" methodology
To maximize the statistical quality of the lag spectrum, data need to be grouped. As the evolution of BHLMXBs happens on typical timescales of days to months, the data grouping needs some criteria to group only the data with similar spectral-timing properties. Therefore, we start from a NICER pilot pipeline (Steiner et al. in preparation) which processes public NICER data and extracts the data products for each good time interval (GTI), including the flux-energy spectrum (both source and background), the light curves (LCs) and the The count rates are for 0.3 − 12 keV, and are normalized for 52 FPMs. The rms is calculated also in 0.3 − 12 keV, with the Fourier frequency range of 1-10 Hz. 1 These values are calculated among the data in the good groups. As we do not have good groups in the soft state (for their low variability) where count rate is usually higher than the hard state, and also with great NICER coverage, the peak count rate and exposure could be underestimated compared to the values taking into account all the archival data. The peak rms could also be underestimated when there is no good group in the faint hard state (for their low count rate). 2 The hard state, intermediate state, and the soft state are represented by H, I, and S respectively; the bold indicates where the soft lag is detected. 3 The first three values correspond to the fractions within which the relative changes of the hardness ratio, the count rate, and the fractional variance are in each group (by default, 0.3, 0.3, and 0.5 respectively); the last value is the minimal SNR ∆t,ind we count as a good group. See text in Section 2.3 for more details.
power spectral densities (PSDs) in several different energy bands.
The GTI-based library of data products enables efficient data grouping. We group data such that the relative changes compared to the first GTI in each group in the hardness ratio, the count rate, and the fractional variance (F var ) are within a fixed fraction (by default, 30%, 30%, and 50% respectively, but the exact values used for each source in our analysis are quoted in Table 1); and the remaining detector combination after filtering out the noisy and hot detectors needs to be the same. Specifically, the hardness ratio is defined as the count ratio of the hard band (4-12 keV) and the soft band (2-4 keV), and the count rate is in 0.3-12 keV. The fractional variance is calculated in 0.3-12 keV and within 1-10 Hz, the range where reverberation lags are expected for BHLMXBs.
As the NICER BHLMXB archive is large, and lag measurements for each data group would be timeconsuming, we introduce a "lag significance indicator" to filter only groups which are promising to generate good quality lag spectra, i.e., those that will result in small lag uncertainties. The phase lag uncertainty is ∆φ(ν j ) = P X,noise P X,signal + P Y,noise P Y,signal + P X,noise P Y,noise P X,signal P Y,signal /2M , where X, Y are the two energy bands between which lag is measured, "noise" denotes the Poisson noise level, and "signal" denotes the Poisson-noise-subtracted intrinsic PSD, and M is the number of LC segments. The last term dominates over the former two terms in the reverberation-dominated frequencies of BHLMXBs (Uttley et al. 2014). Therefore, if the intrinsic lag amplitude is similar and P X,signal /P X,noise ∝ P Y,signal /P Y,noise , the signal-to-noise ratio (SNR) of a lag measurement is where P signal and P noise are the intrinsic PSD and Poisson noise in one energy band. Using the "fractional rmssquared" normalization, the integrated intrinsic PSD over a frequency range is the fractional variance F var due to variability in that range and the above estimate becomes where x is the mean count rate. We therefore define the lag significance indicator as which serves as a quantitative measurement of the expected quality of the lag spectrum. The longer the exposure, the more variable the LC is, the brighter the source gets, the larger the lag significance indicator SNR ∆t,ind is, resulting in a larger SNR of the lag measurement. For instance, a data group of a fairly bright BHLMXB (count rate ∼200 counts/sec) with total exposure of 10 ks (M = 1000 when segment length is 10 s) and a typical fractional rms of 20% in the hard state (F var = 0.04) will have SNR ∆t,ind ∼ 250. In the bright soft state, the fractional rms drops sharply to ∼ 0.5%, so even if the  count rate is as high as 10,000 counts/sec with a very long total exposure of 100 ks, SNR ∆t,ind ∼ 25. From our experience working on previous datasets of GX 339-4 (Wang et al. 2020a) and MAXI J1820+070 (Wang et al. 2021a), we find roughly that an indicator value of ∼ 250 leads to a good-quality lag-frequency spectrum. This is not a strict limit, but rather a tool to help automate the pipeline. Therefore, for each data group, we calculate SNR ∆t,ind and filter only the "good groups" reaching SNR ∆t,ind > 250. In three sources (MAXI J1803-298, MAXI J1727-203, and EXO 1846-031), there was no good group with the default grouping criteria. By relaxing the criteria on the relative changes in hardness ratio, fractional variance, and the count rate, there were a few groups that had indicator values slightly less than 250. We checked the lag-frequency spectra of these sources by hand, and determined that significant soft lags were present despite the slightly lower indicator value, and so we also include these in the sample (see Table 1 for the data grouping criteria used after trials).
We classify the spectral state using the "power-color" (PC) methodology of Heil et al. (2015b). PCs enable the shape of power spectrum to be characterised regardless of its normalisation in the same manner that spectral colors characterise the shape of flux-energy spectrum. We define PC1 as the integrated variance between 0.25-2.0 Hz divided by that between 0.0039-0.031 Hz, and PC2 is defined identically for the 0.031-0.25 Hz and 2.0-16.0 Hz bands. Fig. 2 (b) and (c) shows the power-color diagrams (PC2 vs PC1) of MAXI J1820+070 and all the 10 sources for which we detect soft lag (see Section 3.1 for more details). We see, in line with the results of Heil et al. (2015b), that each outburst follows a circular path on the PC diagram and crucially, each source follows the same path. This is known as the PC "wheel", and provides a very efficient way to define spectral states during an outburst. In a full outburst, the source follows the PC wheel in the clockwise direction in the hard-to-soft state transition (i.e., from the hard state, through the HIMS and SIMS, to the soft state), and follows the same path but in the counter-clockwise direction in the softto-hard state transition (i.e., from the soft state, through the SIMS and HIMS, to the hard state). In particular, the hue is defined as the clockwise angle from the axis defined at a 45 • angle to the x and y-axes. The center of wheel is chosen to be (PC1, PC2)=(4.5, 0.45). The hue ranges for different spectral states are: hard state with hue from -20 • to 140 • , HIMS with hue from 140 • to 220 • , SIMS with hue from 220 • to 300 • , and the soft state with hue from 300 • to 20 • (see Table 2 in Heil et al. 2015b). We note that the original investigation of hue uses the RXTE data in 2-13 keV, and in our analysis, we use NICER data in 4.8-9.6 keV to match the hard X-ray bandpass. We checked the data for MAXI J1820+070 for which the MJD corresponding to the spectral states are known from previous spectral-timing analysis (Buisson et al. 2019;Homan et al. 2020;Buisson et al. 2021), and found that the hue defined in this manner is indeed able to distinguish between hard and hard intermediate states (see Fig. 2 a-b).
For the good groups, we calculate the PSDs, the power-spectral hue, and the lag-frequency spectrum measured between 0.5-1 keV and 2-5 keV. These two energy bands are chosen such that the soft band shows the soft excess, and the hard band is dominated by the coronal emission. We take light curve segment lengths of 10 s for the lag-frequency spectrum and 256 s for the PSD (to match the definition of hue), both with 0.001 s bins. If none of the GTIs in the group is longer than 256 s, the segment length for the PSD is also chosen to be 10 s; in this case, we still count the group as a good group if it meets the before-mentioned criteria because the segment length for the PSD does not affect the lag spectrum, and the only difference lies in the lack of a hue measurement. The logarithmic frequency rebinning factor is 0.4 for the lag-frequency spectrum, and 0.1 for the PSD. The former is chosen to increase the SNR of a lag measurement, and the latter is to to measure the QPO frequency more precisely. These choices do not change our results. Both the source and background flux-energy spectra for each GTI in the data group are combined. The energy resolution of the source spectrum is oversampled by a factor of 3, and is then binned with a minimum count of 1 per channel. The response matrices for different data groups are combined according to the detectors used.
To keep the reverberation machine as automated as possible, we define criteria for determining if reverberation is present in the lag-frequency spectrum. At its most basic, we determine that reverberation is present if a significant soft lag (i.e., a negative lag) is seen in the lag-frequency spectrum. In detail, reverberation is present if at least two frequency bins show a negative lag (i.e., a soft lag) with its absolute amplitude larger than its uncertainty (> 2σ confidence equivalent).
It is further necessary to determine the frequency range over which a soft reverberation lag is present, and that task is performed as follows. As sometimes low statistics can lead to some points having lags consistent with zero, to avoid fragmenting the reverberation frequency range, we require that the range starts from a frequency bin with a negative lag, and continues to higher frequencies as long as the associated lag is consistent with being negative within 1σ confidence. This means the range ends when a positive (hard) lag appears with its amplitude larger than the uncertainty (ν 1 ). In addition, in order to avoid biases due to phase wrapping, we define a "phase wrapping frequency cut-off" (ν 2 ) as the frequency when the absolute lag amplitude becomes more than half of the time lag corresponding to the phase limit ±π. The reverberation lag range ends at either ν 1 or ν 2 , whichever frequency is lower.
The lags at the QPO frequencies can behave differently from their neighboring frequencies dominated by broadband noise, likely due to a different mechanism producing the QPOs (Van den Eijnden et al. 2016; Zhang et al. 2020b). Therefore, as in our earlier work (Wang et al. 2021a), we fit the PSD with the sum of several Lorentzians, and calculate the lag-energy spectra outside of the central frequency ± half-width at half-maximum (HWHM) of those fitted QPOs 1 and their (sub)harmonics. The reference band for the lag-energy spectra is the full band 0.3-10 keV.
All spectral fitting is done with XSPEC 12.11.1 (Arnaud 1996). We use the PG-statistics, for Poisson data with a Gaussian background. We use the wilm set of abundances , and vern photoelectric cross sections (Verner et al. 1996), to account for the galactic absorption with the model tbabs.

Reverberation lag detections
Following the procedures in Section 2.3, we run our autonomous reverberation search algorithm ("reverberation machine") over the sample. We find 10 sources with good data groups (see Table 1), and 16 sources with no good data group (see Table A1). Moreover, all the 10 sources with good groups have soft lag detections (see the detection criteria in Section 2.3). For each source with a reverberation lag detection, we provide a brief introduction to the source properties, the outburst timeline, and the reverberation machine results including the number of good groups, the states covered, and the reverberation lag detections. We also present a summary figure that demonstrates what the data look like in the luminous hard state, and in the state transition (if available). As a case study, we show in this section the results for GX 339-4 which recently went on its first full outburst since NICER's launch. The results for the other 9 sources are detailed in the Appendix.
GX 339-4 is an archetypal BHLMXB which outbursts every 2-3 years (e.g., Belloni et al. 2005;Motta et al. 2009; see also Fig. 1 in García et al. 2019). The mass of the BH is constrained to be between 2.3 and 9.5 M , and our distance towards the system is d > 5 kpc (Heida et al. 2017). It is a mildly absorbed system with the column density of N H ∼ 6 × 10 21 cm −2 (Zdziarski et al. 1998). The inclination via reflection spectroscopy is estimated to lie in the range of 30-48 degrees (e.g., García et al. 2015;Wang-Ji et al. 2018), which makes it a famous low-inclination source. GX 339-4 is significant because it was the first to show evidence for a reverberation feature in the hard state with XMM-Newton (Uttley et al. 2011;De Marco et al. 2015. In the previous two hard-only outbursts in 2017 and 2019, Wang et al. (2020a) detected, for the first time with NICER, a soft reverberation lag with amplitude of 9 ± 3 ms in the Fourier frequency range 2-7 Hz, in the brightest epoch available at that time.
In January 2021, GX 339-4 underwent its first Following the procedures described in Section 2.3, we identify 14 good groups, and a soft reverberation lag is detected in 9 of them, in both the hard and intermediate states (covering data until obsID 4133010105 on 2021-03-28). An example of the reverberation machine summary figure is shown in Fig. 3 for GX 339-4 which includes (a) the hardness-intensity diagram (HID) with a representative hard state group and an intermediate state group shown in blue and red, respectively, (b) the full-band PSDs and any possible QPO feature, (c) the unfolded flux-energy spectra to a constant model, (d) the data-to-model ratio when the 3-10 keV flux-energy spectra are fitted with a (diskbb+powerlaw) model to see phenomenologically the iron line profile, (e)-(f) the lag-frequency spectra, and (g)-(h) the lag-energy spectra in the frequency ranges we detect soft lags. We see that the intermediate state exhibits a much larger lag amplitude that dominates at a lower frequency compared with the hard state; similar to what was previously observed for MAXI J1820+070 (Wang et al. 2021a). We note that for clarity in the summary plots, we show only the lag-energy spectra in the full Fourier frequency band in which we detect soft lags (i.e., not outside of the QPO frequencies), but an example of the lag-energy spectra in the discontinuous frequency bands outside of QPO frequencies are shown in Fig. A1. We also exclude the measured lag within QPO frequencies when we estimate the reverberation lag amplitude for each data group in Section 3.2. The non-detections in the other 5 good groups of GX 339-4 are due to tentative soft lags in the lagfrequency spectra that do not meet the 2σ confidence criteria. Examples of these low-significance lags can be found in Wang et al. (2020a). We also note the detection of a soft lag depends on the chosen energy bands between which we measure the lags. The analysis presented in this systematic study adopts a uniform choice of a hard band of 2-5 keV (with soft band 0.5-1 keV). We have also tested a hard band of 1.5-3 keV as used in Wang et al. (2020aWang et al. ( , 2021a, and in that case, the soft lags in the lag-frequency spectra become more significant for the groups in the hard state, and less significant in the intermediate state. Fundamentally, this difference results from the different shapes of the lag-energy spectra, i.e., the corresponding energy of the "valley" in which the lag is negative is approximately 1.5-3 keV in the hard state, and becomes closer to 2-5 keV during state transition (see panels g and h in Fig. 3).

The evolution of reverberation lag properties throughout the outbursts
With reverberation lags detected in 10 different BH/BHC LMXBs (see Table 1), we can perform a systematic comparison of the reverberation lag properties from source to source, which could help us to understand the disk-corona geometry across spectral states.
As mentioned in Section 2.3, we use the power-spectral hue as an indicator of the state the sources are in during their outbursts (see Section 2.3). The reverberation lag amplitude is approximated with a variance-weighted average of the soft lag amplitude in the lag-frequency spectrum (i.e., the soft lag amplitude is weighted based on its uncertainty σ in each frequency bin we detect a soft lag, by a factor of 1/σ 2 ). We exclude a frequency bin if its centroid lies within the QPO (and their harmonics or sub-harmonics) centroid frequency ± HWHM.
The resulting soft lag amplitude of the 10 BH/BHC LMXBs as a function of the hue is shown in Fig. 4. Our 10 sources show a remarkably similar reverberation lag behaviour as a function of spectral state. We discuss this in more detail later in Section 4.1. As the averaged lag amplitude depends on the level of dilution, a useful quantity to measure the reverberation lag properties is the upper bound of the frequency band where we detect a soft lag, which should not be affected by dilution (while the lower bound can be affected by a change in the lowfrequency hard lags). The evolution of the upper bound of the soft lag frequency band across values of hue is shown in Fig. 5a, together with the fundamental QPO centroid frequency in Fig. 5b. Except for two outliers in the hard state (EXO 1846-031 and MAXI J1803-298, see Section 4.3 for the discussion), all sources appear to follow the same evolution of soft lag amplitude with hue. We therefore average over all sources except for the two outliers by binning on hue. The resulting rebinned evolution of the frequency and amplitude of the soft lag as functions of hue for the remaining 8 sources are shown in Fig. 6, and the soft lag amplitude as a function of QPO frequency is shown in Fig. 7. We summarize the findings as follows: • In the hard state (hue < 140 • ), as hue increases, we see an increase in soft lag amplitude at the beginning (when hue is 20 • , see Section 4.4 for the discussion), and then all the 3 sources with more than 1 reverberation lag detection in this state (GX 339-4, MAXI J1820+070, and MAXI J1348-630) agree on a steadily decreasing trend of lag amplitude (∼ 2 ms to 1 ms, see Fig. 4 and Fig. 6a). Meanwhile, we see a tentative increase of soft lag frequency with hue from tens of Hz to ∼ 100 Hz (see Fig. 5a and Fig. 6b). The evolution of both quantities indicates that the size of the region causing reverberation (i.e., the inner disk-corona region) reduces as the hard state evolves. This is consistent with what has previously been reported for MAXI J1820+070 and interpreted as a contracting corona in the luminous hard state (Kara et al. 2019;Buisson et al. 2019). The QPO frequency increases from ∼ 0.1 Hz to ∼ 1 Hz (Fig. 5b). This can explain why both the Type-C QPO and reverberation evolve to shorter timescales (higher frequencies) before the transition to HIMS. The correlation between QPO and reverberation can also be seen in Fig. 7 as long as the QPO frequency is < 1 Hz, as a higher upper bound of the frequency band over which reverberation occurs generally comes with a shorter lag.
• In the HIMS (140 • < hue < 220 • ), as the hue increases, the reverberation lag amplitude shows a sharp increase by a factor of 10, from 1 to ∼ 10 ms (see Fig. 4 and Fig. 6a). Of the 9 sources with more than 1 reverberation detection, all agree with this increase. In EXO 1846-031, the only hard state detection has an anomalously high lag amplitude and correspondingly low upper bound frequency for which reverberation dominates (see Section 4.3 for the discussion). The reverberation lag amplitudes are mostly consistent (within uncertainties) with MAXI J1820+070, which has the best coverage and the most detections. At the same time, the upper bound of the frequency where reverberation is observed drops from ∼ 100 Hz at the end of hard state to ∼ 10 Hz (see Fig. 5a). In Fig. 6 (b), we see it is not only the upper bound of the reverberation lag frequency that decreases, but also the lower bound, i.e., we could safely say that the reverberation signal evolves to lower frequencies. This result agrees with our MAXI J1820+070 result (Wang et al. 2021a) and provides further evidence for a longer reverberation lag in the HIMS, and is consistent with a growing disk-corona distance, possibly due to a vertically expanding corona sitting at the base of a jet. By comparing the timescale of the reverberation lag to that of the Type-C QPO, we see that unlike in the hard state, during the state transition, the Type-C QPOs continue to evolve to higher frequencies (from ∼ 1 Hz to ∼ 10 Hz), but the reverberation lags, instead, get longer (see Fig. 7 when the QPO frequency is > 1 Hz). We discuss implications of these trends between QPO and reverberation in Section 4.2.
• In the SIMS (hue > 220 • ), out of the 10 sources  in our sample, we have only 2 good groups (in MAXI J1820+070 and MAXI J1348-630) because of the sudden decrease in the variability level compared to the HIMS. A reverberation lag is detected in both data groups, and the mean lag amplitude is even larger (∼ 0.03 s) compared to the HIMS.
To compare the reverberation lag amplitude and the QPO frequency in a more intuitive way, we attempt to unify the HIDs for the 9 sources (except for GRS 1915+105, due to its nonstandard spectral evolution). Simply overplotting the raw HIDs (Fig. 8a) can be misleading because each source is at a different distance, with a different associated galactic absorption, etc. This means, for example, that a hardness of 0.3 for MAXI J1820+070 corresponds to its hard state (according to the power-spectral hue), but the same hardness for MAXI J1535-571 corresponds to the HIMS because MAXI J1535-571 is so heavily obscured. We attempt to scale all sources so that each state (de- The larger galactic absorption diminishes the count rate in the 2-4 keV band when calculating the hardness ratio, and therefore shifts the hardness ratio to higher values, which is exactly what we observe in Fig. 8 (a).
In the 8 sources except for the benchmark MAXI J1820+070, we make linear transformations for both axes, so that observations with similar hues lie in the same region of the HID. The result of this transformation can be seen in Fig. 8 (b). The final result is shown in Fig. 9, where the soft lag amplitudes and the QPO frequencies are presented on the combined HID. One can clearly see that the reverberation lags become much longer in the state transition, while the QPO frequencies continue to increase, which is also shown in Fig. 7. We note that the linear transformation of HID is only for the purpose of presenting the reverberation lag amplitudes and QPO frequencies in a canonical HID. The values of lag amplitudes and QPO frequencies are independent of the linear transformation of HID.

DISCUSSION
With our reverberation machine designed for systematic reverberation lag searches, we have detected soft reverberation lags in 10 BH/BHC LMXBs in the NICER archive. This work has not only increased the sample size by a factor of 5, but also expanded our horizons from the hard state to the intermediate state. We compare the evolution of the soft lag amplitudes in different outbursts in a systematic way, and find that the reverberation lags are short in the hard state (lag amplitudes of ∼ 1 ms) and increase by a factor of ∼ 10 during the state transition (e.g., see Fig. 4).

The disk/corona/jet geometry in the state transition
In the most naive picture, the reverberation lag, τ , relates to the distance between a lamppost corona and reflector, d, as τ = (d/c)[1 + cos i] where i is the inclination angle. Therefore, for a given intrinsic disk-corona geometry in units of R g , a higher BH mass (i.e., a larger R g ) leads to a longer reverberation lag. Alternatively, for a given BH mass, the lag increases with increasing d/R g , which is related to intrinsic change of the diskcorona geometry in units of R g .
Among our sample, MAXI J1820+070 has the most 10 1 10 0 10 1 QPO frequency (Hz) Assuming these values, a reverberation lag amplitude of 1 ms corresponds to a light travel time of ∼ 16 R g /c in the lamppost geometry. Then, a decrease of reverberation lag amplitude from 2 ms to 1 ms in the hard state corresponds to a decrease in the distance between the corona and the inner disk of 32 R g to 16 R g . For MAXI J1820+070, we adopted the reverberation model reltrans to fit the lag-energy spectra in several frequency bands, and measured the coronal height in the lamppost geometry to decrease from 40 to 27 R g in the hard state (Wang et al. 2021a). The reltrans model computes the reverberation lags from the light travel times of the photons using general relativistic ray tracing, and accounts self-consistently for the dilution. The estimate from the proper modeling is close to the naive estimate, and the discrepancy within a factor of 2 is mostly due to dilution. Then, during the state transition, the increase of reverberation lag amplitude to 10 ms leads to a diskcorona distance of ∼140 R g . This is of the same order of magnitude as the fitted coronal height for MAXI J1820+070 with reltrans (starting from ∼ 100− 200 R g and eventually to 300 R g ; Wang et al. 2021a). Moreover, since MAXI J1820+070 has the closest-intime radio coverage to the X-ray state transition, we were able to determine that the radio flare caused by the ballistic transient jet lagged behind the X-ray corona expansion by ∼ 5 days. This finding suggests that the state transition is marked by the corona expanding and ejecting a jet knot that propagates along the jet stream at relativistic velocities, supporting the case for the corona being the base of the jet (see the schematic referred to as "intermediate state" in Fig. 10).
Remarkably, we find that in all the 9 sources with more than 1 reverberation lag detection, as previously reported for MAXI J1820+070 (Wang et al. 2021a), the reverberation lag gets longer in the state transition. This result shows this trend is a generic feature of state transitions of BHLMXBs, not depending on, e.g., inclination, spin, or mass. Unfortunately, we cannot test for a correlation between lag amplitude and BH mass, inclination, or spin, as most of our sample consists of new transients lacking robust measurements of these quantities. For different BHLMXBs, the dependency of the lag amplitude on the BH mass and inclination can be seen from the naive picture at the beginning of this sec-tion. The lag amplitude could also depend on the BH spin because it could change the disk-corona distance in units of R g in two ways: (1) the innermost stable circular orbit (ISCO) radius monotonically decreases with the BH spin; (2) as the jet could be powered by the BH spin (see e.g., Narayan & McClintock 2012;Steiner et al. 2012), the geometry of the corona (as base of a jet) could thus be related to the BH spin. In addition, the inclination could only cause a difference in lag amplitude by a factor of 2 in the naive picture, and the BH mass is usually the same within a factor of 2. That is to say, the reverberation lag amplitude is expected to have a narrow dynamic range, making it very challenging to test for its correlation with physical quantities.
The generic trend of reverberation lag getting longer is highly suggestive of what happens in such state transitions. The reverberation lag evolution during the state transition cannot result merely from dilution, because the upper bound of the Fourier frequency where reverberation is detected also evolves to lower values (see Fig. 5a and 6b). This upper frequency bound is independent of dilution. Therefore, a growing size of the disk-corona region that causes reverberation is necessary. In our modeling of MAXI J1820+070 data, we find an increase in the coronal height possibly due to a vertically expanding corona.
We also note that a new time lag model for accreting BHs could explain the hard lags and soft (reverberation) lags with a relatively compact corona (Uttley & Malzac 2021, in preparation; see the poster). Uttley & Malzac suggested that the power-law pivoting delay which produces the hard power-law continuum lags, modelled em-pirically in reltrans, is due to the delay expected between seed photon variations and coronal heating as the driving mass accretion fluctuations propagate from the disk (providing the seed photons) to the corona (where coronal heating occurs). In this scenario, hard lags are produced at lower frequencies and soft lags are produced at higher frequencies, due to accretion propagation effects. The soft lags result from the lag of the reverberation signal (which responds primarily to the power-law emission in the harder X-rays, driven mainly by coronal heating) relative to the seed photons from fluctuations propagating through the disk (which cause the powerlaw flux in the soft band to rise before the power-law flux in the hard band). However, these lags also depend on the coronal geometry, so that more vertically extended coronae show larger soft lags, i.e., the inferred geometric change is the same as for the reverberation lags being due to light-travel times only. However, in the case where the soft lags are linked to propagation, the vertical scales required to produce the lags are reduced.
The idea of the corona expanding in state transition was previously suggested on different grounds (Pottschmidt et al. 2000;Nowak et al. 2002). The authors used RXTE data of Cygnus X-1 and GX 339-4 and found longer hard lags during state transition. The interpretation was that these prolonged lags were produced in spatially very extended synchrotron-emitting outflows. More recently, this result was generalized into the photon-index-time-lag correlation in BHXBs (e.g., Grinberg et al. 2013;Reig et al. 2018). With NICER's soft X-ray bandpass and large effective area, we are now able to find prolonged soft reverberation lags in state transitions, which provide more direct evidence of the expanding corona since the mechanism that produces hard lags is not well understood. In addition to longer lags, recent broadband spectral modeling for several BHLMXBs finds that the radius of the corona which is assumed to be the base of the jet increases from ∼ 10 R g to ∼ 40 R g through the hard state and the HIMS; and the reflection fraction occasionally becomes very low in the HIMS when approaching the SIMS, suggesting an increase in also the vertical extent of the corona and/or its bulk velocity Cao et al. 2021, but see Steiner et al. 2017).

The QPO/reverberation connection
We compare the timescale of the reverberation lag to that of the Type-C QPO in Fig. 5. In the hard state, the reverberation lags evolve to shorter timescales (higher frequencies) before the transition. Similarly, the Type-C QPOs also evolve to higher frequencies. In fact, all characteristic frequencies in the power spectrum increase during the outburst, not just the QPO frequencies (Wij-nands & van der Klis 1999; Psaltis et al. 1999). However, during the state transition, the Type-C QPOs continue to evolve to higher frequencies, but the reverberation lags, instead, get longer (Fig. 7). In short, reverberation and Type-C QPO are correlated in the hard state, but they anticorrelate during the state transition. By modeling the reflection spectrum with a dual-lamppost model (two lampposts at different coronal heights) to mimic a vertically extended corona, Buisson et al. (2019) also found the upper coronal height to decrease first in the bright hard state, and then increase when approaching the state transition in MAXI J1820+070, while the characteristic features of the temporal variability kept increasing in frequency.
Type-C QPOs' amplitudes show great evidence for an inclination dependence (e.g., Motta et al. 2015;Heil et al. 2015a), which suggests that they arise from some geometric effect, e.g., Lense-Thirring precession of an inner hot flow or corona (Ingram et al. 2009). The evolution of the QPOs to higher frequencies has been proposed as evidence for a smaller inner hot flow or central corona (along with a less truncated disk) before the source reaches the soft state. This could also explain the increasing disk blackbody temperature, the softening coronal emission with a decreasing coronal temperature (Done et al. 2007), and also the shorter reverberation lags in the bright hard state (Kara et al. 2019). However, the anticorrelation of the QPO frequency and reverberation lag in the state transition breaks this simple picture.
If both the QPOs and reverberation lags depend mostly on the geometry of the corona-inner-disk system as discussed before, the anticorrelation seems to suggest that the corona is both radially and vertically extended. While the QPO frequency is mostly coupled to the radial extent, the reverberation lag is coupled to both the radial and vertical extents. The idea that the corona has both radial and vertical extents is very natural, and is previously suggested by reverberation features in AGN, which require both a radially extended corona over the disk surface and a vertically extended one sitting at the launching point of a jet (Wilkins et al. 2016(Wilkins et al. , 2017. If the QPO is produced by the Lense-Thirring precession of an inner hot flow and/or the jet (Ma et al. 2021;Liska et al. 2021), which shrinks in the radial dimension before the source reaches the soft state, the longer reverberation lags require an additional component of the corona on the vertical axis of the BH that expands during the state transition. In the narrow line Seyfert 1 galaxy Mrk 335, the emissivity profile of the reflection spectrum indicates the radially extended corona becomes gathered up in the center and is then launched away from the BH, reminiscent of a jet, during a flaring event (Wilkins & Gallo 2015). This scenario is very similar to the physical picture we propose here to explain the evolution in the QPO frequency and reverberation lag in the hard state and the state transition of BHLMXBs. More detailed spectral-timing modelings assuming extended coronal geometries and incorporating QPO mechanisms are required to understand the QPO/reverberation connection.

Outliers with long reverberation lags
In the hard state, there are two outliers that show long soft lags of ∼ 10 ms compared to the bestsampled sources (MAXI J1820+070, GX 339-4, and MAXI J1348-630) with ∼ 1 ms soft lags (see Fig. 4). At the same time, the upper bound of the soft lag frequency is only ∼ 10 Hz, which is 10 times smaller than that of MAXI J1820+070 (see Fig. 5a). These are EXO 1846-031 (Fig. A7) and MAXI J1803-298 (Fig. A10). The reason for the much longer reverberation lag amplitude in EXO 1846-031 is likely that it has the largest column density among the 10 sources in our sample (N H ∼ 10 23 cm −2 ). The absorption diminishes significantly the photons at soft X-rays (< 1 keV), and therefore the lag spectral quality is not adequate to confidently detect the reverberation lag and measure its amplitude, especially when we search for reverberation in the soft band of 0.5-1 keV. The limit of lag spectral quality can be seen from the large uncertainties in the lag-frequency spectrum where soft lag is detected, and the lag-energy spectrum below 1 keV in Fig. A7. The soft lag amplitude is measured to be 16 ± 11 ms, also with a large uncertainty. In fact, we find the 0.5-1 keV band is dominated by the shelf of the response matrix with such high absorption. This means that the photons registered in these energy channels are actually the photons redistributed downward from higher energies, most likely ∼ 2-3 keV which is the peak area band with remaining signal. However, we note this effect does not make our soft reverberation lag detection artificial because an intrinsic hard lag (i.e., the hard band when measuring the lag-frequency spectrum, 2-5 keV, lags behind the ∼ 2-3 keV band) remains a hard lag. Encouragingly, we observe a very promising iron lag whose energy is not affected by the high galactic absorption in EXO 1846-031. If, instead, we use the tentative iron lag amplitude, this source is no longer an outlier in our sample.
The other outlier, MAXI J1803-298, is not affected by large galactic absorption (N H ∼ 3 × 10 21 cm −2 ; Homan et al. 2021), and indeed in this source, the measured lag is highly significant. The longer reverberation lag is more likely to be intrinsic to the system. Candidate explanations include that MAXI J1803-298 contains a more massive black hole or a different disk-corona geometry. Perhaps this system is unique compared to the rest of the sample because it showed a long-lived excursion to the "steep power-law state" (Ubach et al. 2021), which may suggest a different coronal geometry relative to the other sources.

A closer look at reverberation lags in the hard state
In the hard state, we see an increase in soft lag amplitude as a function of hue when hue is 20 • (see Figs. 4 and 6a). To explore this further, we plot the soft lag amplitude in the hard state as a function of the luminosity in Eddington units 2 for the three sources with more than one soft lag detection in the hard state (MAXI J1820+070, GX 339-4, andMAXI J1348-630) in Fig. 11. In the initial rise of the hard state in MAXI J1820+070's 2018 outburst, the source luminosity in NICER's obsIDs 1200120101-1200120105 increased from 1% to 5% L Edd . During this stage, the hue decreased from −3 • to −6 • , and then increased to 2 • . At the same time, the reverberation lag amplitude decreased from 25 to 6 R g /c, and then increased to 15 R g /c. This correlation between hue and reverberation lag amplitude results in the increase of soft lag amplitude with hue at 20 • . After the initial rise, MAXI J1820+070 was in a plateau where the luminosity slightly decreased from 20% to 18% L Edd , hue gradually increased from 20 • to 120 • , and the soft lag amplitude decreased from 30 to 10 R g /c. After the plateau, the source went through a decay and then a rise, before finally making the hard-to-soft state transition (see e.g., Fig. 2a). During this stage, the hue spanned the range of 110 • to 140 • and correlated with luminosity (in the range of 6% to 11% L Edd ), while soft lag amplitude (in the range of 7 to 10 R g /c) anticorrelated with luminosity. Therefore, after the initial rise, the soft lag amplitude kept decreasing with hue in the hard state as seen at hue 20 • in Figs. 4 and 6(a). In GX 339-4 and MAXI J1348-630, the outburst history is not as complicated, and the data groups in the hard state are all during the initial rise.
In Fig. 11, we can see that for L/L Edd 2%, the reverberation lag amplitude decreases with increasing 2 MAXI J1820+070 has the most reliable measurements of its BH mass (8.48 +0.79 −0.72 M ; Torres et al. 2020), distance (2.96±0.33 kpc; Atri et al. 2020), and inclination (63±3 degrees; Bright et al. 2020). From near-infrared study, the latest measurement of the mass function of GX 339-4 is f (M ) = 1.92±0.08 M , and the mass ratio of companion star and the BH is q = 0.18±0.05; the distance to GX 339-4 has only a lower limit of ∼ 5 kpc, and a larger distance is favored (Heida et al. 2017). If a distance of 10 kpc and an inclination of ∼ 40 • (García et al. 2015;Wang et al. 2020a) are assumed, the mass becomes 10 M from the mass function. As for MAXI J1348-630, the distance is assumed to be 3.4 kpc, the mass to be 10 M (Lamer et al. 2021)  ) and distance (2.96±0.33 kpc) of MAXI J1820+070 as it is the benchmark whose HID is not linearly transformed. (Lower ) The evolution and geometry of the disk-corona-jet system we infer in the faint hard state, bright hard state, and the intermediate state. In the faint hard state (< 1 − 2% L Edd ), the disk is truncated at several tens of Rg and the inner disk fills in, causing the reverberation lag to become shorter (see Section 4.4 and Fig. 11). In the bright hard state (> 4% L Edd ), the disk is no longer truncated, and the X-ray corona vertically contracts from several tens of Rg to ∼ 10 Rg (see Section 4.1 and Fig. 11), while the steady jet fades. In the intermediate state, the X-ray corona expands to ∼ 100 Rg, while the steady jet is quenched and the transient jet is ejected (see the discussion in Section 4.1). In the soft state, there is no good data group (due to the low variability) or soft lag detection, but we note that there is a great amount of evidence that the disk extends to the ISCO radius (e.g., Steiner et al. 2010); the disk emission (with a disk temperature of ∼ 1 keV) dominates over the coronal emission (Done et al. 2007), and the radio jet is off (Fender et al. 2004). During the soft-to-hard state transition and the following decay in the hard state, the soft lag amplitudes are consistent with those during the rise at earlier time in the outbursts. However, as the number of detection at this stage is limited, future data and more detailed modeling are needed to test for any hysteresis effect during the rise and decay in the hard state, or in the upper and lower branches of the intermediate state.
luminosity, as expected from previous results. However, in the 2% L/L Edd 4% range, the lag amplitude surprisingly appears to increase with increasing luminosity, before once again decreasing with increasing luminosity for 4% L/L Edd 10%. Among the 3 sources, only MAXI J1820+070 shows the plateau stage at ∼ 20% L Edd , where both the luminosity and the lag amplitude steadily decrease. We also note that this finding has no conflict with previous reverberation lag measurements in De Marco et al. (2015, as the increase of lag amplitude we see in Fig. 11 occurs in a luminosity range not covered in those papers. Although the increase of reverberation lag between 2% and 4% L Edd is subtle because of the large uncertainties of the lag amplitudes, it is interesting that we see this behaviour in all 3 sources with more than 1 soft lag detection in the hard state. A possible interpretation is that below 1-2% L Edd , the disk is truncated at the order of several tens of R g , and above that, the disk is no longer truncated (the inner edge of the disk is only could only be seen if we use the soft lags detected between 0.5-1 keV and 1.5-3 keV because the soft lag is not significant enough to detect when using 2-5 keV as the hard band. The values from previous reverberation lag measurements (De Marco et al. 2015 are the "observed" lag amplitudes (without considering the dilution as shown in Table 3 in De Marco et al. 2017), and the L/L Edd values are converted accordingly with our assumed BH mass and distance. See Section 4.4 for the assumed BH mass, distance, and inclination, and how the dip results in the correlation of soft lag amplitude and hue in Fig. 4 and 6(a) when hue is < 20 degrees.
several R g ), which is consistent with the spectral analysis of GX 339-4 (García et al. 2015;Wang-Ji et al. 2018; see the schematic referred to as "faint hard state" in Fig. 10). The levels of disk truncation also match the soft lag amplitudes in light travel time (in units of R g /c, see Fig. 11 below 2% L Edd ). The change from a truncated disk to a non-truncated one is very likely to be due to the radiative efficiency of the inner disk increasing above ∼ 1% L Edd , which is theoretically predicted (Yuan & Narayan 2014) and supported by the change of slope in the radio/X-ray correlation (e.g., Koljonen & Russell 2019). Instead, above 2% L Edd , the increase and then decrease of soft lag amplitude is more likely to be driven by changes in the coronal geometry, which would be consistent with an evolving base of the compact jet. Recent GRMHD simulations (Liska et al. 2019) show that, as the inner disk is more radiatively efficient and thus thinner, the inflowing material exerts less pressure on the jet boundary, which becomes less collimated, resulting in the base of the jet being more extended compared to lower luminosities Cao et al. 2021). The decreasing soft lag amplitude above 4% L Edd could then be caused by a contracting X-ray corona before the quenching of the compact radio jet (see the schematic referred to as "bright hard state" in Fig. 10). This is similar to the expanding X-ray corona that happens before the transient radio jet ejection in the state transition. However, the interpretation should be taken with some caveats. First, the evolution of soft lag amplitudes is tentative because of the large uncertainties of the measured lags. Second, the uncertainties in BH mass, distance, and the bolometric correction factor preserve the trend, but shift the values of luminosities in Eddington units. Further detailed reflection modeling of the data in the hard state (especially measurements of the coronal height, reflection fraction, and the inner edge of the disk), and comparison with infrared/radio data are needed.

SUMMARY
We have constructed a NICER reverberation machine aiming to systematically search for reverberation signatures in the BH/BHC LMXBs in the NICER archive, and to study the reverberation lag evolution throughout the outbursts. Our major findings are as follows: 1. Out of 26 BH/BHC LMXBs NICER has observed, we find data groups with high enough signal to noise to enable a soft lag detection in 10 sources, and soft reverberation lags are detected in all 10 of these sources. Detections have previously only been reported for MAXI J1820+070 and GX 339-4, and therefore for the other 8 we report the first ever reverberation detections; increasing the sample size of BHLMXBs with reverberation lag detections from 2 to 10 (3 to 11 when accounting for detections with XMM-Newton in H1743-322).
2. We adopt the power-spectral hue to quantify what stages the sources are in during their outbursts, and find a consistent reverberation lag evolution in the hard state and state transition for the 10 sources.
3. In all the 9 sources with more than 1 reverberation detection, the reverberation lag is longer and dominates at lower Fourier frequencies during the hard-to-soft transition. This reverberation lag evolution has been reported only for MAXI J1820+070, and now with our systematic searches, our result suggests this trend is a global property of state transitions of BHLMXBs that does not depend on inclination, BH mass or spin. This reverberation evolution suggests a growing Xray emitting region, possibly due to an expanding corona.
4. The frequency range where we see reverberation correlates with the Type-C QPO frequency in the hard state, but they become anticorrelated during the state transition. This result hints at a corona that is both radially and vertically extended, and while the QPO frequency is mostly coupled to the radial extent, the reverberation lag is coupled to both extents.
JW, GM, EK and JAG acknowledge support from NASA ADAP grant 80NSSC17K0515. JAG thanks support from the Alexander von Humboldt Foundation. AI and DA acknowledge support from the Royal Society. RAR acknowledges support from NICER grant 80NSSC19K1287.

A. REVERBERATION LAG DETECTIONS IN INDIVIDUAL SOURCES
There are 26 BH/BHC LMXBs observed by NICER. Among those, the reverberation machine finds good data groups and detects soft reverberation lags in 10 sources (see Table 1). The reverberation machine summary plots include the spectral-timing results for a group in the bright hard state and a group in the intermediate state (if available). We note that an example of the lag-energy spectra of GX 339-4 in the discontinuous frequency bands outside of QPO frequencies is shown in Fig. A1. In the 10 reverberation machine summary plots (Figs. 3, A2-A10), we show only the lag-energy spectra in the full Fourier frequency band in which we detect soft lags for clarity. The results for GX 339-4 are detailed in Section 3.1, and those for the other 9 sources are presented below. The reverberation machine finds no good data group in the other 16 sources, and the peak count rate, peak rms, total exposure, the states achieved, and the reasons for no good groups are shown in Table A1.

A.1. MAXI J1820+070
MAXI J1820+070 is one of the brightest BHLMXBs detected in history as it reached 5.8 Crab at the peak luminosity in the 2018 outburst. The brightness and great exhibition of multi-wavelength activity across spectral states make it an extensively studied new BHLMXB (e.g., Homan et al. 2020;Tetarenko et al. 2021). The BH mass is measured to be 8.48 +0.79 −0.72 M (Torres et al. 2020), and the distance to be 2.96±0.33 kpc (Atri et al. 2020). The inclination constrained from the jet motion is 63±3 degrees .
Thanks to the high brightness and near-daily cadence coverage of this source, we find 88 good groups (making use of data with obsID 1200120101 in the initial rise of hard state to obsID 2200120309 in the quiescence), and detect soft reverberation lags in 84 of them. The mean lag significance indicator is 3115, which is the highest among all 10 BHLMXBs with detected reverberation lags. The non-detection groups lie in the lower branch of HIMS (i.e., soft-to-hard transition) approaching the hard state for the significance of promising soft lag being lower than 2σ. We note that with NICER, we have performed detailed time lag analysis for MAXI J1820+070 with a conventional method of manually grouping the data (Kara et al. 2019;Wang et al. 2021a). With the new automatic grouping, we verify in Fig. A2 that the major result of Wang et al. (2021a), i.e., that the reverberation lag frequency decreases during transition with growing lag amplitude, is well reproduced.

A.2. MAXI J1348-630
MAXI J1348-630 is a bright X-ray transient reaching 4.5 Crab at its peak luminosity in its Jan 2019 outburst (Tominaga et al. 2020), after which 4 hard-only mini-outbursts were observed (Zhang et al. 2020a). With H I absorption, the distance of MAXI J1348-630 is measured to be 2.2 +0.5 −0.6 kpc (Chauhan et al. 2021), and the jet motion constrains the ejecta inclination to our line of sight to be i ≤ 46 • (Carotenuto et al. 2021). Moreover, MAXI J1348-630 was believed to host a relatively massive black hole (Tominaga et al. 2020), but with the more recent distance measurement, the evidence has become weaker because the lower limit of M BH is only 3.7 M assuming a * = 0 and i = 0 • (Chauhan et al. 2021).
We find 34 good groups and detect soft reverberation lags in 15 groups. Among them, 14 groups are during the major outburst covering both the hard and intermediate states, and the remaining 1 group is in the first hard-only mini-outburst after the major outburst. The non-detections all lie in the hard state because of the intrinsic shape difference of the lag-energy spectrum compared to the HIMS. The evolution in the soft reverberation lag frequency and amplitude is consistent with what has been detected in MAXI J1820+070 and GX 339-4 (see Fig. A3). We note that NICER had a visibility gap for this source between MJD 58513.2 (obsID 1200530105) when it was in the bright hard state and MJD 58521.9 (obsID 1200530106) when the source was already in the HIMS, during which time a relativistic transient jet was estimated to be ejected on MJD 58518.9 ± 2.4 (Carotenuto et al. 2021). Therefore, we cannot conclude when the state transition starts relative to the transient jet launching, nor can we precisely determine when the X-ray reverberation lag becomes longer during the state transition.
A.3. MAXI J1535-571 MAXI J1535-571 is a LMXB that was first detected by MAXI on 2017 September 2 (Negoro et al. 2017). By fitting the soft state spectrum with a relativistic thin disk component kerrbb, the black hole mass is estimated to be 10.39 +0.61 −0.62 M with a distance of 5.4 +1.8 −1.1 kpc (Sridhar et al. 2019). Reflection spectroscopy suggests an inner disk inclination in the range of 55-77 degrees (Xu et al. 2018;Miller et al. 2018); but the jet inclination inferred from the motion of the superluminal knot is ≤45 degrees . Notice that MAXI J1535-571 has a relatively high interstellar absorption with N H = (3 − 8) × 10 22 cm −2 (Kennea et al. 2017;Gendreau et al. 2017;Xu et al. 2018), which makes the soft lag more challenging to measure than in less obscured sources. In the radio, a transient jet knot was estimated to be ejected on MJD 58010 +2.65 −2.5 , and the compact jet was rapidly quenching on MJD 58013.6 (Russell et al. 2020).
NICER started observing this transient from 2017 September 9 (MJD 58005.3), corresponding to the observation (obsID 1050360103) that composes the first good group we find with our reverberation machine (see Fig. A4). In this data group, no GTI is longer than 256 s to measure the power down to 0.0039 Hz. As one ingredient in measuring the hue is the power in 0.0039-0.031 Hz, we cannot use the hue to determine the spectral state for this data group. Instead, we confirm its hard state nature from spectral-timing features. The QPO centroid frequency is relatively low at 0.9 Hz. The photon index is ∼ 1.8, and the disk is cool at T in ∼ 0.18 keV when the flux-energy spectrum is fitted in 0.5-10 keV range with a basic reflection model of (tbabs*diskbb+reltransDCp) following Wang et al. (2021a).
The other 13 good groups are all in the HIMS covering MJD 58008.5-58038.4 (obsID 1050360104 to 1130360115) 3 , and the 2 non-detections happen along with the two lowest lag significance indicators. A representative group in the HIMS is also shown in Fig. A4, and we observe again that the frequency over which soft reverberation lag dominates decreases during state transition. Because of the gap between MJD 58005.3 and 58008.5, and the large uncertainty in the transient jet ejection time (MJD 58010 +2.65 −2.5 ), we cannot determine the time between the X-ray corona expansion and the launch of the transient jet.
A.4. MAXI J1631-479 MAXI J1631-479 was discovered as a bright hard X-ray source on 2018 December 21 by MAXI/GSC (Kobayashi et al. 2018), and was in the hard state until December 26 (MJD 58478;Fiocchi et al. 2020). MAXI J1631-479 is likely to be a low-inclination source based upon several arguments including from reflection spectroscopy with NuSTAR (i = 29±1 • ; Xu et al. 2020), the shape of the HID (Monageng et al. 2021), and the phase lag of its Type-C QPO (van den Eijnden et al. 2019). Both the black hole mass and distance are currently unknown. A radio flare was observed to peak on MJD 58560 to 58582 (though notably, the HIMS-to-soft-state transition happened earlier, between MJD 58545 and 58560). The delay of the radio flare hints at the possibility of multiple unresolved flares and/or jet-ISM interactions (Monageng et al. 2021). NICER started monitoring this source on MJD 58498 (obsID 1200500101) when it was already in the soft state, and caught a transition back to the HIMS on MJD 58506 (obsID 1200500110; van den Eijnden et al. 2019). The source stayed in the HIMS until MJD 58545 (obsID 2200500102), after which NICER had a visibility gap and the source was already back to the soft state when the monitoring resumed on MJD 58560 (obsID 2200500108; Rout et al. 2021). As no NICER data is available in the hard state, all the good groups we find are in the HIMS, and the results for one representative group is shown in Fig. A5. Soft reverberation lags are detected in 9 out of 12 good groups. These lags are always found at low Fourier frequencies (i.e., even at 0.1 Hz, no prominent hard lag is observed), consistent with results from previous sources in the HIMS.
3 It was suggested that the source entered the SIMS between MJD 58014.2 and 58015.4 (Tao et al. 2018). Meanwhile, all the good groups we find in the intermediate state are determined to be in the HIMS based on the hue. There is no controversy because by comparing with the observations classified to be SIMS with the NICER data, our good groups do not cover any of those data segments (see Fig. 1 in Stevens et al. 2018). This is owing to the criteria of the relative change in fractional variance when grouping the data, and the requirement on the lag significance indicator that is most sensitive to the rms.

A.5. MAXI J1727-203
MAXI J1727-203 was discovered by MAXI/GSC on 2018 June 5 (Yoneyama et al. 2018). NICER observed this source since the day it was discovered, and caught the initial HIMS followed by the soft state, the soft-to-hard transition, and the final hard state (see Fig. 1 in Alabarta et al. 2020). Spectral-timing analysis with NICER data has been reported in Alabarta et al. (2020), and the nature of the compact object in this system is by far the most uncertain.  For the first time, we show here the lag analysis for this source. There are only 3 good groups found with our reverberation machine: one in the initial HIMS (MJD 58275.0-58276.5, obsID 1200220102 to 1200220103), one in the HIMS during the transition back to the hard state (MJD 58298.6-58301.4, obsID 1200220120-1200220123), and one in the hard state at the end of the outburst (MJD 58305.6-58327.8, obsID 1200220127-1200220140). The group in the hard state has only prominent low-frequency hard lags, and has no strong evidence for a high frequency soft lag, probably for the same reason as discussed in Section 3.1 (see the blue data group in Fig. A6). The soft reverberation lag is detected in both HIMS groups (e.g., Fig. A6). Although the soft lag < 1 keV is not as significant as in previous sources, we note that the iron lag is very significant in the group shown. This confirms that reverberation is present in this source. Interestingly, the relativistic iron line in the flux-energy spectrum is not as prominent as the iron lag in the lag-energy spectrum. This type of discrepancy was reported previously for several AGN, and was thought to suggest the inadequacy of the lamppost corona assumption (Zoghbi et al. 2020;Mastroserio et al. 2020). The detection of its reverberation lag (especially the iron lag) hints at a black hole nature for MAXI J1727-203.
A.6. EXO 1846-031 EXO 1846-031 was first discovered in 1985 (Parmar & White 1985). After 25 years in quiescence state, it underwent an outburst in 2019, and was detected first by MAXI (Negoro et al. 2019). Draghis et al. (2020) performed reflection spectroscopy analysis with NuSTAR data, and found that EXO 1846-031 has a very high spin (a * ∼ 0.997), a high inclination (i ∼ 73 • ), and a high column density of N H ∼ (10−12)×10 22 cm −2 , which makes it difficult to dynamically measure the black hole mass in the optical. From the model kerrbb, the authors estimate a black hole mass of 9±5 M . Both the extremely high spin and its position in the radio/X-ray plane ) hint at a black hole nature for this system.
After its renewed activity, NICER quickly followed up on 2019 July 31 (MJD 58695, obsID 2200760101) and confirmed the source was in the hard state (Bult et al. 2019). With our reverberation machine, we have found 3 good groups, 1 in the initial hard state, and 2 in the intermediate state. A soft reverberation lag is detected in the first two groups (see Fig. A7). We note that although the soft X-ray data quality is largely affected by the large column density towards this source (see also Section 4.3), an iron lag is present in the lag-energy spectrum in the hard state data group. On the other hand, the shape of the lag-energy spectrum for HIMS data group is certainly not a typical hard lag (which  would monotonically increase with energy without any feature). Therefore, we confirm that reverberation is present in these two data groups. By comparing results from these two groups, we see the frequency over which soft reverberation lag is detected decreases during the hard-to-soft transition, and the iron line seems to get blueshifted, exactly as seen previously in MAXI J1820+070 (Wang et al. 2021a). The nature of the blueshifting is not well understood, and physically-driven models that are preliminarily promising include the emission from the plunging region (Fabian et al. 2020) and returning radiation (Connors et al. 2020). Other possibilities include extra Comptonization of the reflection spectrum in the corona (Steiner et al. 2017), a change in the inclination (i.e., warps) in the reflecting surface of the inner disk, a change in the disk thickness as the mass accretion rate rises, or a change in the irradiation profile due to changing coronal geometry.  Figure A8. Reverberation machine summary plot for GRS 1915+105 (see Section A.7 for more details). All the good data groups are in the HIMS determined from the power-spectral hue.
GRS 1915+105 is the first source of detected superluminal jets in our galaxy, and therefore is the first "microquasar" (Mirabel & Rodriguez 1994). It displays at least a dozen of variability states that are believed to result from disk-jet coupling (e.g., Belloni et al. 2000). The BH mass is measured to be 12 ± 2 M , with an orbital period of 33.5 days (Reid et al. 2014). GRS 1915+105 is a peculiar BHLMXB as it remained in a bright outburst for 26 years since its discovery in 1992 (Castro-Tirado et al. 1992). Since its launch, NICER frequently observed GRS 1915+105 in 2017(until 2017. This dataset was analysed in Neilsen et al. (2018) and the well-known disk wind evidence including Fe XXV and Fe XXVI absorption was detected. From 2018 April to 2020 March, GRS 1915+105 surprisingly declined in count rate, and eventually settled in a very faint and hard state Homan et al. 2019). Very recently, in 2021 July, the source was detected to be rebrightening (Neilsen et al. 2021).
Our reverberation machine results in 3 good data groups all during the declining stage (making use of ob-sID 1103010147 to 2596010901, corresponding to 2018-05-27 to 2019-05-01). A soft reverberation lag is detected in the first two groups. The results of the first group is shown in Fig. A8. The lag spectral quality is limited, and the soft lag amplitude is measured to be 8 ± 5 ms, with a relatively large uncertainty. However, we note that this measurement is within uncertainties with MAXI J1820+070 (see Fig. 4). We find a broad emission line at 6.4 keV and a prominent narrow absorption feature at 6.97 keV which corresponds to Fe XXVI line from highly ionized wind, consistent with previous NICER result in the more luminous state (Neilsen et al. 2018). We also note that the hue for these groups suggests that GRS 1915+105 was in the HIMS in the decay stage.
A.8. AT2019wey AT2019wey was initially discovered as an optical transient in December 2019 (Tonry et al. 2019), and later was classified as a galactic XRB for its strong X-ray emission and hydrogen lines detected at redshift z = 0 (Yao et al. 2020c). The source is expected to be a LMXB because the mass of the companion star was constrained to be < 0.8 M (Yao et al. 2020b). The multi-wavelength (radio to X-ray) spectral-timing features are consistent with a BHLMXB nature for AT2019wey (Yao et al. 2020b;Mereminskiy et al. 2021). A resolved radio source was detected and could be interpreted as the steady compact jet common in the BHLMXB hard state (Yadlapalli et al. 2021). Currently, we lack a mass estimate for this source, and the inclination is constrained to be i ≤ 30 • from reflection spectroscopy (Yao et al. 2020a).
During the NICER coverage, AT2019wey was in the hard state at the beginning from MJD 59070 to MJD 59082 (obsID 3201710105-3201710117), and then transitioned into the HIMS as reported by Yao et al. (2020a), which is in great agreement with spectral states classified by the hue in our reverberation machine. Though we have found in total 22 good groups (6 in the hard state and 16 in the HIMS), the fraction of good groups with soft reverberation lag detected is by far the lowest: only 3 out of 22, all in the HIMS. In addition, the lag-energy spectra found in the soft lag frequencies this way are not representative of a soft reverberation lag (i.e., there is no strong evidence for reverberation in either the soft X-ray band or the iron line energies, see the upper panel in Fig. A9). If we use 1.5-3 keV instead of 2-5 keV as the hard band in calculating the lag-frequency spectrum, the lag-energy spectra show quite promising reverberation lag features (see the lower panel in Fig. A9). The detection of reverberation lag may suggest that this source is a BHLMXB.
The reason for a very low reverberation lag detection fraction and a better-suited hard band of 1.5-3 keV instead of 2-5 keV is again intrinsically the shape of the lag-energy spectrum. This would mean that even in the HIMS where we observe a prominent reverberation lag in the other sources when we match the hard band of lag-frequency measurement to the energy of the "valley" in the lag-energy spectrum (2-5 keV), AT2019wey still has a lag-energy spectrum shape similar to the other sources in the hard state (valley at 1.5-3 keV). This would imply a lower disk temperature, a lower ionization, or a lower disk density. We note that this system has a very short orbital period (< 8.6 hours), and a low mass accretion rate that could result in the "hard-only" outburst implied in Yao et al. (2020a).
A.9. MAXI J1803-298 MAXI J1803-298 was first discovered by MAXI/GSC in May 2021 and quickly followed up by NICER (Serino et al. 2021;Gendreau et al. 2021), making it the most newly discovered XRB in our sample. It was initially detected in the hard state between MJD 59336 and 59340 (obsID 4202130101-4202130104), and was already in the intermediate state on MJD 59352.7 (obsID 4202130105) when NICER came out of a visibility gap Ubach et al. 2021). With the promising nearly-periodic X-ray absorption dips seen by NICER and NuSTAR, and the observed disk wind, it is likely to be a high-inclination source (i > 60 • ), and the orbital period is between 7 and 8 hours Xu & Harrison 2021;Miller & Reynolds 2021).
There are 2 good groups in our reverberation machine: one in the hard state (obsID 4202130103-4202130104), and the other in the HIMS (obsID 4202130107, very close to the SIMS judging by the value of hue), both are shown in  Figure A9. Reverberation machine summary plot for AT2019wey. The soft lag is detected in 3 out of 22 good data groups, but the reverberation nature is difficult to confirm because of a lack of the iron line feature in the lag-energy spectrum (upper ). If we use 1.5-3 keV instead of 2-5 keV as the hard band in calculating the lag-frequency spectrum, the lag-energy spectra show more promising reverberation lag features (lower, see Section A.8 for more details). Fig. A10. The reverberation lag is detected only in the hard state group, while prominent low-frequency hard lags are observed in the intermediate state group. Notice that the intermediate state group ) has a type-B QPO detected at 6.7 Hz with a sub-harmonic at 3.3 Hz, which is consistent with what was reported in Ubach et al. (2021). It was surprising at first to see that only in MAXI J1803-298, the low-frequency hard lag remains evident in the intermediate state, where for previous sources, we see the reverberation lag dominates over the hard lags at the low frequencies. One possibility is that the intermediate state here is different from the canonical HIMS, and rather, it is in the steep power-law (SPL) state (Remillard & McClintock 2006), as suggested by Ubach et al. (2021) based on a relatively high scattering fraction. In Fig. A10, besides the blue-shifted iron line in the intermediate state, we also find evidence for absorption feature in the Fe band at 6.97 keV, suggesting the system is observed closer to edge-on as the disk wind is expected to be confined to the equatorial plane, in agreement with what was seen with Swift (Miller & Reynolds 2021).   Table A1. The BH/BHC LMXB samples observed by NICER without any good reverberation data from our reverberation machine, sorted by the peak count rate.