© ESO 2021

1 Introduction

Black hole X-ray binaries (BHXRBs) are ideal laboratories for studying in detail the evolution of the accretion flow as a function of the accretion state during an outburst (e.g. Done et al. 2007). During the soft state, thermal emission with a peak temperature of ~1 keV is observed. The low or undetected rapid X-ray variability (on timescales between milliseconds and hundreds of seconds; e.g. Gierliński & Zdziarski 2005; Belloni et al. 2005) suggests the accretion flow is stable against fast changes in the mass accretion rate in this state (Lyubarskii 1997; Churazov et al. 2001). A geometrically thin, optically thick accretion disc (Shakura & Sunyaev 1973; Novikov & Thorne 1973) extending close to the innermost stable circular orbit (ISCO) of the black hole (BH) is the accepted paradigm to explain the thermally dominated soft state (e.g. Gilfanov 2010). However, this standard model is unable to explain the hard spectral states, which are instead dominated by rapidly variable (Gierliński & Zdziarski 2005; Muñoz-Darias et al. 2011; Heil et al. 2012) emission from a hot plasma with peak temperatures of a few tens to hundreds of keV (e.g. Joinet et al. 2008; Motta et al. 2009; Buisson et al. 2019; Zdziarski et al. 2021b). In the hard state the primary broadband X-ray continuum more likely results from Compton up-scattering processes of less energetic photons (thermal and/or synchrotron) by the electrons in the hot plasma (Poutanen & Coppi 1998; Gierliński et al. 1999; Poutanen & Vurm 2009). On the other hand, the standard disc is observed to be cooler (kT ≲0.3 keV; e.g. Kolehmainen et al. 2014; Zdziarski & De Marco 2020) and more variable (Wilkinson & Uttley 2009; De Marco et al. 2015a) than in the soft state. Transitions between hard and soft states are characterised by abrupt changes in the X-ray spectral and timing properties of the source (e.g. Belloni 2010; Bogensberger et al. 2020). The intermediate states during which the transition occurs have been further classified into hard-intermediate and soft-intermediate states, based on their significantly different X-ray timing properties (Homan & Belloni 2005; Belloni et al. 2005). The observed changes characterising the intermediate states are likely a consequence of major, but still not well-understood, variations in the physical properties of the accretion flow.

All the models proposed to explain the hard state postulate that the hot phase of the accretion flow is located at relatively small distances from the BH; however, its exact geometry is unknown, and a number of viable possibilities currently exist. For example, the hot flow might fill the inner regions of a truncated optically thick and geometrically thin accretion disc, with the transition region characterised by an overlap between the two phases (i.e. a corona sandwiching the inner disc and/or disrupted disc clumps embedded in the hot flow; e.g. Done & Zycki 1999; Yuan & Narayan 2014; Poutanen et al. 2018). Alternatively, the standard disc might reach close to the BH and the hard X-ray source may be a compact region, possibly located on the rotationaxis of the BH (the so-called lamppost; e.g. García et al. 2014; Niedźwiecki et al. 2016) and linked to the base of a jet (Miller et al. 2006). As a matter of fact, radio observations indicate the presence of a bipolar relativistically expanding and outflowing plasma associated with this state (Fender et al. 2004; 2006; Corbel et al. 2001; 2004; Brocksopp et al. 2002; Miller-Jones et al. 2011). However, the exact relation (if any) between the X-ray source and the jet is still unclear (e.g. Yuan et al. 2019).

The main challenge for these models is explaining the high level of X-ray spectral and timing complexity that characterises the hard state, including: complex broadband spectral curvature (Makishima et al. 1986; Nowak et al. 2011); evolution of characteristic variability timescales (Done et al. 2007); detection of hard X-ray lags (hard band flux variations lagging behind soft band variations) in Comptonised emission-dominated regions of the spectrum (Nowak et al. 1999; Grinberg et al. 2014); and the recurrent appearance of low-frequency quasi-periodic oscillations (QPOs) during high luminosity hard and hard-intermediate states (e.g. Rodriguez et al. 2002; Casella et al. 2004; Sobolewska & Życki 2006; Ingram & Motta 2019). Truncated disc models can (at least qualitatively) explain most of the observed X-ray spectral and timing properties of BHXRBs and their evolution during an outburst, provided some spectral stratification of the hot flow is postulated for the hard state (Done et al. 2007; Axelsson & Done 2018; Mahmoud & Done 2018a). Within this scenario the key ingredient is an evolving inner disc radius, moving closer to the BH as the source evolves from the hard to the soft state, and receding again as the source goes back to quiescence (Esin et al. 1997; Yuan & Narayan 2014). On the other hand, a compact X-ray source with a disc inner radius at ~ R_ISCO already at moderate Eddington ratios is usually invoked to explain the broad profile of Fe K emission lines in the disc reflection component, which are observed particularly in bright hard states (Miller et al. 2006; García et al. 2015; Fürst et al. 2015). Within this model, variations in the vertical extent of the primary hard X-ray source have been proposed to explain the evolution of some of the X-ray spectral and timing properties throughout the hard state (Kara et al. 2019; Buisson et al. 2019).

Theoretical arguments and supporting observational evidence exist for the optically thick and geometrically thin disc being highly truncated during quiescence and in the early stages of an outburst (R_in ~ 10³⁻⁴ R_g for M_BH ~ 10 M_⊙, where R_g = GM_BH∕c² is the gravitational radius; e.g. Lasota et al. 1996; Narayan et al. 1997; Dubus et al. 2001; McClintock et al. 2001; Esin et al. 2001; Bernardini et al. 2016). Nonetheless, the question of when the disc reaches the ISCO (whether early in the hard state, in the bright hard state, or during the intermediate states characterising the hard-to-soft state transition) remains unanswered. Attempts to settle this issue have so far yielded controversial results (see Bambi et al. 2021 for a review). Among the main problems is the difficulty in obtaining a long, almost uninterrupted monitoring of an outburst with high throughput X-ray detectors covering energy ranges with significant contribution from all the main spectral components. In particular, while the rise in luminosity throughout the hard state can be quite slow (of the order of a few weeks to months), the hard-intermediate states preceding the transition to the soft-intermediate states are spanned in a very short amount of time (a few days), and thus they are very difficult to catch (Belloni et al. 2005; Dunn et al. 2010). Yet, these phases are of the utmost importance for understanding how the accretion flow evolves.

The Neutron star Interior Composition Explorer (NICER; Gendreau et al. 2016) is making significant inroads in this context by performing high cadence, high quality X-ray observations of Galactic transient sources. Among the sources observed by NICER is the very bright BHXRB MAXI J1820+070 (ASASSN-18ey), which was discovered in March 2018 (Tucker et al. 2018). During the 2018 outburst the source rapidly brightened, reaching extremely high X-ray luminosities, with a peak X-ray flux of about 4 Crab. Such a high brightness and low absorbing column density to the source (N_H ~ 10²¹cm⁻², Uttley et al. 2018) make MAXI J1820+070 an excellent target for detailed studies of BHXRB outbursts and state transitions. Intensive monitoring at radio wavelengths, followed by high spatial resolution Chandra observations, revealed the ejection of long-lived discrete relativistic jets (Bright et al. 2020; Espinasse et al. 2020). Major changes in the X-ray variability properties of the source preceded these ejections by a few hours (Homan et al. 2020), suggesting that these events marked the transition from the hard to the soft state.

Since its discovery, NICER has been observing the source regularly, providing an incredibly rich archive of data covering all the phases of the outburst. Here we focus on the first part of the 2018 outburst, covering the hard and hard-intermediate states before transition to the soft-intermediate states. We employ two independent methods to constrain the inner accretion flow geometry and its evolution. Both methods rely on the application of X-ray spectral-timing analysis techniques (Uttley et al. 2014) to study the quasi-thermal soft component produced by the hard X-ray irradiation of the inner disc (e.g. Zdziarski & De Marco 2020). Since the hard X-ray flux impinging on the cold disc medium is variable, it will induce a temporal response in the disc quasi-thermal emission, which is the focus of this work.

The first method (Sect. 4.1) is based on studying the evolution of high-frequency (ν ≳1 Hz) soft X-ray lags, namely soft band flux responding to rapid hard band flux variations with a time delay. These lags are thought to be the signature of X-ray thermal reverberation, namely the time-delayed response of the re-emitted quasi-thermal disc emission to variations in the hard X-ray irradiating flux (e.g. Uttley et al. 2011). Being strongly dependent on the light path between the X-ray source and the reprocessing region, X-ray reverberation lags can be used to map the geometry of the inner accretion flow and its evolution. To date, such lags have been detected in some of the best monitored BHXRBs (Uttley et al. 2011; De Marco et al. 2015b, 2017; De Marco & Ponti 2016; Wang et al. 2020a), including MAXI J1820+070 (Kara et al. 2019).

The second method (Sect. 4.2) is based on studying the evolution of covariance spectra (Wilkinson & Uttley 2009; Uttley et al. 2014), which are used to single out the spectral components responsible for the observed thermal reverberation lags. The peak temperature of the quasi-thermal component responding to the hard X-ray irradiation of the inner disc is linked to the disc inner radius, thus providing an alternative way of constraining the geometry of the inner accretion flow (Zdziarski & De Marco 2020; Zdziarski et al. 2021a).

For our computations we make use of the latest results from the optical spectroscopy of MAXI J1820+070, which yield a BH mass of M_BH = (5.95 ± 0.22)sin⁻³i M_⊙ (Torres et al. 2020). Hereafter we assume a fiducial value of M_BH = 7.6 ± 1.7 M_⊙. This encompasses the range of values inferred from current constraints on the orbital (66° < i < 81°; Torres et al. 2020), and jet inclination (63 ± 3°; Atri et al. 2020; 64 ± 5°, Wood et al. 2021). We use the distance to the source of d = 2.96 ± 0.33 kpc (where the uncertainty is for 68 per cent confidence; Atri et al. 2020).

2 Observations and data reduction

We analysed data from NICER X-ray Timing Instrument (XTI) observations of MAXI J1820+070 between MJD 58 189 and MJD 58 306. These observations cover the first part of the 2018 outburst, namely the hard and hard-intermediate accretion states, up to the transition to the soft-intermediate states (the observations analysed are listed in Table 1, the observation identification number, ObsID, increases progressively with the time of the observation; we note that our sample also includes the observations previously analysed in Kara et al. 2019). For brevity, the different observations will be hereafter referred to as Oxxx, with xxx being the last three digits of their ObsID.

The data were reprocessed with the NICERDAS tools in HEASOFT v6.28 and NICER calibration files as of July 27, 2020. Calibrated, unfiltered, all Measurement/Power Unit (MPU) merged files (ufa) were created using NICERL2 task. We applied standard filtering criteria (e.g. Stevens et al. 2018), using NIMAKETIME and NICERCLEAN routines. We extracted light curves in the energy range 13− 15 keV (corresponding to the range where NICER XTI effective area quickly drops) to check for periods of high particles background. Time intervals showing background events with rate >2 counts s^-1 (e.g. Ludlam et al. 2018) were discarded.

Of the 56 focal plane modules (FPMs) of NICER XTI, FPMs 11, 20, 22, and 60 are not operational. Following ‘NICER Analysis Tips & Caveats¹ ’, FPMs 14 and 34 were removed from the analysis as they can exhibit increased detector noise. The single observations were screened to identify additional FPMs showing anomalous behaviour, which, for example, can manifest as the total number of registered counts significantly exceeding the average counts from all detectors. However, none was found during the analysed observations. Occasionally (particularly during the brightest phases of the outburst before the transition to the soft-intermediate states), the number of active FPMs was reduced to avoid telemetry saturation and an increase in dead time. Once accounting for this reduction and the outlined filtering criteria, the net number of FPMs used for this analysis ranges between 50 and 27, as reported in Table 1.

The resulting event lists were used to extract light curves in different energy bands with a time bin of 0.4 ms. Custom ancillary response files (ARFs) and redistribution matrix files (RMFs) were computed for the specific subset of FPMs of each observation. This was done by combining the corresponding per-module, publicly distributed ARFs and RMFs. Fits were performed using Xspec v.12.11.1 and errors are hereafter reported at the 90 per cent confidence level.

In order to speed-up the reduction process, short good time intervals (GTIs) with length < 10 s were removed. This led us to the exclusion of observations O128, O129, O154, O181, O191, O192, O193, because of their resulting very short net exposure (≲200 s). The effective, on-source, exposure times of analysed observations after screening and filtering of short GTIs are listed in Table 1.

3 Hardness-intensity diagram

We built the hardness-intensity diagram (HID; Fig. 1) of the source in order to select observations corresponding to the accretion states of interest. To this aim we extracted NICER count rates in the energy ranges 2–4 keV and 4–12 keV. The total 2–12 keV count rate was normalised for the number of FPMs used for the analysis (N_FPM, Table 1). In some cases N_FPM changed during a single observation. In these cases the observation was split so as to have datasets characterised by a constant number of FPMs (each point in the HID corresponds to either a single NICER observation or the part of an observation characterised by constant N_FPM). We note that the HID of Fig. 1 reports also observations excluded from the analysis because too short (Sect. 2).

During the first part of the outburst, we identified four main phases: an initial brightening at almost constant hardness (which lasted for about 12 days, ‘rise’ in Fig. 1), a persistent (lasting about 2 months) bright hard state (‘plateau’ in Fig. 1), a sharp dimming phase (lasting about one month) at slightly lower, and almost constant hardness (‘bright decline’ in Fig. 1), and the quick (lasting about two weeks) transition to the soft-intermediate state (‘hard-soft transition’ in Fig. 1). These phases are the focus of our analysis. We note that the last analysed dataset includes only the first part (~64 ks) of O197, corresponding to an on-source effective exposure of ~12 ks. This spans the phases immediately before the fast transition from the hard-intermediate to the soft-intermediate state (i.e. the ~ 47 ks preceding the transition, of which ~8 ks of effective time) and the transition itself. The excluded part of O197 corresponds to the soft-intermediate state. In this phase X-ray spectral-timing measurements are noisier due to a drop in the intrinsic variability of the source. As shown by Homan et al. (2020), apart from the drop of broadband X-ray variability strength, the transition from the hard-intermediate to the soft-intermediate state was also marked by an intense X-ray flare and the switch from type-C to type-B QPO at MJD~58 305.7. This epoch is highlighted in the HID of Fig. 1. Also shown in Fig. 1 is the estimated epoch of the launch of discrete relativistic ejecta, detected first at radio wavebands (Bright et al. 2020) and later on, at larger scales, in the X-rays (Espinasse et al. 2020). An earlier ejection has been recently reported (Wood et al. 2021), coinciding with O197, and preceding the switch to a type-B QPO. The empty circles in Fig. 1 show the spectral changes of the source during the remaining part of the 2018 outburst, till MJD 58 404.

Fig. 1 HID of MAXI J1820+070 during its 2018 outburst (from MJD 58 189 to MJD 58 305). The hardness is computed as the ratio between the 4–12 keV and 2–4 keV NICER count rates. On the y axis we report the total 2–12 keV NICER count rate re-normalised for the number of used FPMs. Highlighted in different colours and symbols are the phases of the outburst that are the focus of our work. The grey empty circles show the subsequent spectral evolution of the source up until MJD 58 404. The epochs of the appearance of a type-B QPO (MJD ~ 58 305.7; Homan et al. 2020) and the launch of relativistic discrete radio ejecta (MJD 58 305.6 ± 0.04; Wood et al. 2021, and ~MJD 58 306, Bright et al. 2020; Espinasse et al. 2020) are marked, respectively, by the bigger star and orange arrows.

4 X-ray spectral-timing analysis

We measured the Fourier cross-spectrum to extract frequency-dependent time lags and covariance spectra (e.g. Nowak et al. 1999; Uttley et al. 2014). Given that the on-source exposure time of each NICER observation is usually quite short (a few hundred to a few thousand seconds; see Table 1), we adopted the approach of combining consecutive observations in order to increase the signal-to-noise ratio (S/N) of spectral-timing products. Before combining consecutive observations we verified thatthey do not show deviations from stationarity (Vaughan et al. 2003). The main reason for this is that non-stationary behaviour is likely associated with substantial changes in the physical properties of the accretion flow, possibly including geometry, whose evolution we aim to test here. Details on this selection are reported in Appendix A. In addition, we combined only datasets with the same number of FPMs. This is relevant for the fit of covariance spectra (Sect. 4.2), since the response of the instrument depends on the specific FPMs used for the analysis (Sect. 2). The final groups of combined observations can be inferred from Table 1.

At the beginning of the outburst, the power spectral density (PSD) changed relatively slowly. In this phase of the outburst we split large groups of observations having compatible PSDs into smaller groups each characterised by sufficiently long net exposure time (~10–15 ks). On the other hand, before the hard-soft transition the PSD changed quickly, thus preventing us from combining consecutive observations and obtaining longer exposures (e.g. O191, O192, and O193 were discarded because of the short exposures and the inability to combine them Sect. 2).

4.1 The evolution of high-frequency soft reverberation lags

Soft X-ray lags (soft variations lagging behind hard variations) in BHXRBs are interpreted as a signature of thermal reprocessing in the disc. Light travel time delays due to the extra path between the X-ray source and the disc are thought to contribute significantly to the observed soft lags (thermal reverberation; e.g. Uttley et al. 2011; De Marco et al. 2015b; 2017). As a consequence, these lags can be used as a diagnostic of the geometry of the innermost accretion flow. Soft X-ray lags ascribable to thermal reverberation were consistently observed in the plateau phase of MAXI J1820+070, with clear indications of an evolution of lag amplitude and frequency (Kara et al. 2019).

We measured soft X-ray lags throughout the first part of the outburst of MAXI J1820+070. Following results reported in Zdziarski et al. (2021b) we selected a soft (0.5–1 keV) and a hard (2–5 keV) energy band. The selected soft band shows a soft excess that includes some contribution from disc emission (Fig. 6 in Zdziarski et al. 2021b), the hard band is instead dominated by a hard X-ray Comptonisation component. Although the analysis of Zdziarski et al. (2021b) does not cover all the phases considered here, we verified that these choices remain appropriate even in the latest phases. We did not include photons with energies > 5 keV to avoid strong contribution from Fe K line photons. We computed frequency-dependent time lags between the selected bands as τ(ν) = ϕ(ν)∕2πν, where ϕ(ν) is the phase of the cross spectrum as a function of frequency (e.g. Nowak et al. 1999). We corrected the lags for dead-time-induced cross-channel talk (van der Klis et al. 1987; Vaughan et al. 1999), which introduces instrumental ± π-phase lags. However, we verified that up to very high frequencies (~100 Hz) the contribution from the source dominates, thus making cross-channel talk negligible.

A number of examples among the obtained lag-frequency spectra are shown in Fig. 2 (which reports the broad-frequency band ~0.03−100 Hz lag spectrum) and Fig. 3 (which focuses on the high-frequency ~ 2− 100 Hz range behaviour). The lags are re-binned using a multiplicative re-binning factor of 1.2, with each frequency bin containing at least 500 points (Ingram 2019). At low frequencies (≲1–2 Hz) the lag-frequency spectra are dominated by (positive) hard X-ray lags (variations in hard X-ray photons follow variations in soft X-ray photons; see Fig. 2). These lags are commonly observed in BHXRBs at frequencies ≲ 1 Hz and are thought to be intrinsic to the primary hard X-ray continuum (e.g. Nowak et al. 1999; Pottschmidt et al. 2000; Grinberg et al. 2014; Rapisarda et al. 2016).

At higher frequencies a (negative) soft X-ray lag is clearly observed in almost all the analysed observations. Only during observations O101 and O102 the detection is marginal². The soft X-ray lag typically starts to be observed at frequencies ≳2 Hz (inset inFigs. 2 and 3), after getting to its maximum amplitude (in absolute value) it starts decreasing, reaching zero-amplitude at a given very high frequency (hereafter ν₀).

As clear from Fig. 3, soft X-ray lags appear to vary both in amplitude and in characteristic frequencies throughout the analysed observations. The maximum measured amplitudes range between ~ 3 ms and ~ 0.2 ms³. However, the most compelling observed phenomenon is the net increase in lag amplitude in the last analysed dataset (O197), as compared to previous ones (Fig. 3, bottom-right panel). Here the soft X-ray lag amplitude suddenly increases to ~15–20 ms and dominates the entire sampled frequency range (see bottom panel of Fig. 3), with no clear evidence of hard X-ray lags (down to ~0.05 Hz). This dataset corresponds to the phases immediately before and during the transition as identified by the appearance of a type-B QPO (Homan et al. 2020). However, the disappearance of the hard lags/appearance of a long ~15–20 ms lag precedes the appearance of the type-B QPO. Indeed, we verified that the long, dominant soft X-ray lag is already present when considering only the ~ 50 ks preceding the transition itself (i.e. the first ~8 ks of effective on-source time during O197, Sect. 3). Such a long reverberation lag is observed also right after the transition from the soft to the hard state (from O262 to O265). This will be studied in more detail in a follow-up paper. It is worth noting that a much longer lag (~20 s) has been reported by Buisson et al. (2021) in a NuSTAR observation performed on MJD 58 306.35 (quasi-simultaneous to NICER observation O198), that is, right after the transition to the soft-intermediate state.

The evolution of soft X-ray reverberation lags informs us about the evolution of the geometry of the system. From a visual inspection of the lag-frequency spectra, the observed maximum amplitude of the reverberation lag clearly shows a complex evolutionary pattern (the trend of observed maximum reverberation lag amplitude as a function of spectral hardness is shown and discussed in more detail in Appendix B). However, it is well known that the presence of both primary continuum and reprocessed photons in the energy bands used for the computation of cross-spectra has the effect of diluting the amplitude of soft lags (e.g. Uttley et al. 2014). In particular, it can be easily shown that the amount of dilution also depends on the amplitude of the lags intrinsic to the primary X-ray continuum, being minimum when the continuum has no intrinsic interband lags. Nonetheless,the primary continuum in BHXRBs is characterised by (positive) hard X-ray lags, which show a complex evolution during the outburst (e.g. Pottschmidt et al. 2000; Grinberg et al. 2014; Altamirano & Méndez 2015; Reig et al. 2018). As clear from Fig. 2, this complexity characterises also MAXI J1820+070 (see also Wang et al. 2020b). As a consequence, the hard X-ray lags in the continuum can bias the maximum observed amplitude of the reverberation lag, and their behaviour throughout the outburst can potentially distort the observed evolution of soft reverberation lags.

Nonetheless, the frequency at which the soft lag goes to zero, ν₀, has been shown to be unaffected by dilution of continuum photons (e.g. Fig. 8 of Wilkins & Fabian 2013; Uttley et al. 2014; Mizumoto et al. 2018). This frequency is also unaffected by the presence of hard X-ray lags associated with the primary continuum, because these lags are commonly observed to steadily decrease with frequency, eventually becoming negligible (e.g. Nowak et al. 1999; De Marco et al. 2015b). This behaviour is also seen in MAXI J1820+070, as shown in Fig. 2. Thanks to the suppression of hard X-ray lags towards high frequencies, ν₀ directly depends on the reverberation lag intrinsic (and unbiased) amplitude as . We note that ν₀ is unbiased by the presence of hard lags even when these lags are not totally suppressed at high frequencies (see Appendix C). Therefore,ν₀ can be used to study, in a model-independent way, the evolution of soft reverberation lags.

In order to constrain ν₀ in an objective way we modelled the frequency-dependent lag profile. Because of the marginal detection of soft X-ray lags during O101 and O102, and given that no significant differences in the lag are observed, only for this part of the analysis these two observationswere combined. For the fits we used the RELTRANS spectral-timing model (Ingram et al. 2019) in the version presented inMastroserio et al. (2019; 2020). This model assumes a lamppost geometry for the source of Comptonised hard X-ray photons anddescribes the hard X-ray lags as due to spectral pivoting of the inverse-Compton spectrum (Mastroserio et al. 2018, 2019). Hard X-ray lags generally show hump-like features at specific frequencies which might be associated with a more structured Comptonising region (e.g. Mahmoud & Done 2018a; 2018b). These complexities are clearly observed in the lag-frequency spectra of MAXI J1820+070 (see Fig. 2) and did not allow us to always obtain good fits of both the hard and soft lags on a broad range of frequencies. Therefore, we excluded most of the hard lags complex profile from the fits, restricting the modelling of the lag-frequency spectra to frequencies ranging between ~1 and 2 Hz and 100 Hz. For O197 we used the frequency range 6−100 Hz because of the anomalous behaviour observed at lower frequencies (lack of visible hard lags) that did not allow us to obtain good fits over a broader frequency range. Since this observation catches the source right before and during the transition, it is likely to present additional spectral-timing complexities not included in current models.

We note that, for a given value of the BH mass, the parameters of the model influencing the position of ν₀ are the X-ray source height, h, and the disc inner radius, R_in. However, these parameters are degenerate when fitting lag-frequency spectra alone. For this reason, and because the geometry of the Comptonising region might be more complex than assumed in the model, we do not discuss the inferred best-fit values for these parameters, and use the best-fit models only to obtain a phenomenological description of the soft X-ray lag profile. Additional technical details of these fits are reported in Appendix D.

The derived best-fit models provide a satisfactory description of the high-frequency lag spectrum at each epoch (see the examples in Fig. 3, where the red solid curve represents the best-fit model). We estimated ν₀ as the frequency at which the best-fit model approaches zero at high frequencies⁴. The errors on ν₀ are obtained by accounting for the uncertainty in the parameters of the model (as described in Appendix D). The derived values of ν₀ are listed in Table E.1. These values are plotted as a function of spectral hardness in Fig. 4.

For most of the analysed observations (more specifically those corresponding to the rise, plateau, and bright decline phases), the data show a steady trend of increasing ν₀ with decreasing spectral hardness (Fig. 4 left panel). A linear model in log-log space to describe this trend in the rise, plateau, and bright decline is preferred to a constant model at > 99.99 per cent confidence level (i.e. corresponding to a Δχ² ~ 154 for a difference of 1 degree of freedom; Fig. 4 left panel).

Notably this trend appears to break during the hard-intermediate states preceding the transition (hard-soft transition phase in Fig. 4, right panel). Indeed, ν₀ is observed to decrease by a factor of ~2–3 at spectral hardness ≲0.28, remaining systematically lower than observed in most of the bright hard state. During observation O197 (right before and during the transition) we estimate Hz (Table E.1).

An increase in ν₀ (as characteristic of most of the analysed part of the outburst, Fig. 4, left panel) indicates a decrease in the intrinsic lag amplitude, and thus a decrease in the relative distances between the Comptonising plasma and the reprocessing region in the disc. This might be due to the inner disc truncation radius moving inwards or the Comptonising region decreasing in its height orvertical extent. The behaviour close to transition (Fig. 4, right panel) is instead more puzzling and might be associated with the presence of a jet (Bright et al. 2020; Homan et al. 2020; Espinasse et al. 2020, see Sect. 5.1).

The intrinsic amplitudes that can be inferred from ν₀ correspond to Euclidean distances of the order of several tens-to-hundreds of R_g (see also discussion in Sect. 5.2). However, these distances represent the weighted mean over all the light travel paths of the Comptonised photons to the disc (e.g. Wilkins & Fabian 2013), and can be significantly larger than the minimum travelled distance set by the X-ray source height and the inner disc radius (Mahmoud et al. 2019; see also Sect. 5.2). The distribution of light travel paths depends on the geometrical parameters of the disc and the Comptonising plasma, and is encoded in the impulse-response function (e.g. Gilfanov et al. 2000; Poutanen 2002; Wilkins & Fabian 2013). However, simultaneous fits of cross-spectra and time-averaged spectra, possibly with models that account for more complicated geometries for the X-ray source, are needed in order to constrain the impulse-response function. In the following section we employ an alternative technique to obtain tighter constraints on the geometry of the inner accretion flow responsible for the observed soft reverberation lags in MAXI J1820+070.

Fig. 2 0.5–1 keV vs. 2–5 keV lag-frequency spectrum of some of the analysed observations of MAXI J1820+070 during the rise, plateau, bright decline, and hard-soft transition phases (with colour codes and symbols as in Fig. 1). The positive lags dominating at low frequencies (≲2 Hz) are hard X-ray lags. The plot illustrates the hump-like structures characterising these low-frequency hard lags and their complex evolution throughout the first part of the outburst. At high frequencies (≳2 Hz) the (negative) soft X-ray lag becomes dominant (see inset, with the y axis in linear scale).

Fig. 3 0.5–1 vs. 2–5 keV lag-frequency spectra of selected observations of MAXI J1820+070. The observations illustrate the diverse behaviour of soft (negative) lags during different phases of the outburst (as also reported in the insets). The solid red curve is the best-fit RELTRANS model. The blue arrow marks the frequency range used for the extraction of covariance spectra.

Fig. 4 Evolution of ν0 in lag-frequency spectra as a function of spectral hardness. Left panel: results from the rise, plateau, and bright decline phases, which are observed to follow the same trend. The dotted black line is the best-fit model. Right panel: results from the hard-soft transition phase as compared to the other phases (grey symbols), which is observed to break from the observed trend at higher spectral hardness.The arrows are 90 per cent lower and upper limits.

4.2 Constraining the inner radius of the disc from covariance spectra

The hard X-ray flux that irradiates the inner parts of the accretion disc is partly reflected (i.e. producing absorption features and emission lines plus a Compton scattering continuum) and partly reprocessed as a quasi-thermal continuum in the soft X-ray band. The latter adds up to the disc black body component due to internal dissipation (e.g. Gierliński et al. 2008; Zdziarski & De Marco 2020) and is expected to give a significant contribution in accretion states energetically dominated by the hot Comptonising plasma. Therefore, the observed disc black body temperature at the inner radius of the disc is T_in = f_colT_eff, where T_eff is the effective temperature due to the internal dissipation and irradiation processes (such as ), and f_col > 1 is the colour-correction factor accounting for spectral hardening (e.g. Ebisawa et al. 1993; Shimura & Takahara 1995; Davis et al. 2005; Salvesen & Miller 2021).

As discussed in Zdziarski & De Marco (2020) and Zdziarski et al. (2021a), T_in can be estimated from the Stefan-Boltzmann law, as: (1)

where a is the disc albedo, (1 − a) is the fraction of irradiating flux reprocessed as a quasi-thermal component, F_irr is the total irradiating hard X-ray flux impinging on the inner radius of the disc R_in, and σ is the Stefan-Boltzmann constant.

Assuming a standard emissivity profile for the X-ray source (i.e. ∝ R^−q, with the emissivity index q = 3) and a small scale height compared to the disc radius (such as in a coronal geometry or a lamppost relatively close to the BH horizon, this assumption is discussed in Sect. 5.1), the flux irradiating the inner edge of the disc is (Zdziarski & De Marco 2020). Therefore, following Zdziarski & De Marco (2020), Eq. (1) translates into a lower limit on R_in: (2)

where l_irr = 4πd²F_irr∕L_Edd (assuming isotropy) is the total (i.e. including contribution from both variable and constant components) Eddington-scaled (with L_Edd for pure hydrogen) luminosity of the irradiating hard X-ray source, and k is the Boltzmann constant. The factor is the reflection fraction and was introduced by Zdziarski & De Marco (2020) to account for the possible reduction of hard X-ray flux impinging on the disc (e.g. as a consequence of anisotropic emission; see discussion in Sect. 5.2). The approximate equality in Eq. (2) holds for the case of null disc internal dissipation (). Under this condition an estimate of the disc truncation radius can be obtained. We note that this additional quasi-thermal component is incorporated in state-of-the-art reflection spectra (e.g. Tomsick et al. 2018), but for the remainder of the paper we fit it as a separate component in order to constrain its temperature.

The quasi-thermal emission due to irradiation responds to the fast variability of the hard X-ray flux (thermal reverberation), as inferred from the detection of high-frequency soft X-ray lags (Sect. 4.1). Previous analyses showed that the thermal response of the disc to variable irradiation also produces a soft, disc black body-like component on top of the hard X-ray Comptonisation continuum in high-frequency covariance spectra (Wilkinson & Uttley 2009; Uttley et al. 2011). These spectra single out components varying in a linearly correlated way (coherently) and produced in the innermost regions of the accretion flow, thus can be used to study the quasi-thermal emission produced via irradiation of the inner disc.

We estimated the temperature of the quasi-thermal component responding to hard X-ray irradiation of the inner disc, from the fits of high-frequency covariance spectra (kT_in,covar). We emphasise that these are the frequencies where the soft X-ray lags ascribed to thermal reverberation are observed, thus allowing us to constrain the geometry of the region of the disc producing these delays. Indeed, the temperature inferred for this region was used together with Eq. (2) in order to obtain constraints on the inner radius of the irradiated disc. In order to extract the covariance spectra we chose a broad energy band as reference (0.5–10 keV) and computed cross-spectra between the reference band and a series of adjacent energy bins. We followed the standard approach of removing the contribution of each energy bin from the reference band before computing each cross spectrum (see Uttley et al. 2014 for details, and Ingram 2019 for a different approach). The cross-spectra were averaged over a frequency interval of interest. Since we are interested inthe thermal reprocessing component, we considered the frequency intervals where soft X-ray reverberation lags are detected (Sect. 4.1). In order to exclude as much contribution as possible from the process responsible for the hard X-ray lags (which produces correlated disc and hard X-ray continuum variability at low frequencies; e.g. Wilkinson & Uttley 2009; Uttley et al. 2011; De Marco et al. 2015a), we focused only on the highest sampled frequencies. Specifically, we averaged the cross-spectra over the frequencies that range between the maximum observed soft lag amplitude and τ ~ τ(ν₀) (the frequency intervals for the lags displayed in Fig. 3 are indicated by the blue arrows). In general, the tested frequencies range between ~6 and ~90 Hz. Since kT_in,covar is expected to increase with frequency (as higher frequencies test regions of the disc closer to the BH), we verified that even limiting our analysis to the highest-frequency end of the chosen frequency intervals does not influence significantly the results reported here.

We first fit the covariance spectra of each observation at energies E ≥ 3 keV with a Comptonisation model, absorbed by a column of cold gas (TBabs*nthComp in Xspec, Zdziarski et al. 1996; Życki et al. 1999; Wilms et al. 2000). The seed photons temperature was fixed at 0.2 keV, the photon index, and normalisation were left free to vary. Since NICER data are not sensitive to the Comptonisation high energy cut-off, the temperature of the Comptonising plasma was fixed at kT_e = 30 keV (Zdziarski et al. 2021b). The column density was fixed at N_H = 1.4 × 10²¹cm^-2 (Kajava et al. 2019; Dziełak et al. 2021). The best-fit models were then extrapolated down to E = 0.3 keV, revealing significant residuals at soft energies. In order to fit such residuals, we added a quasi-thermal contribution from a disc component. However, since we are fitting covariance spectra (which are essentially obtained from the difference between the time-dependent spectrum and the time-averaged spectrum of the source, Wilkinson & Uttley 2009) we need to account for the fact that the observed thermal emission in such spectra deviates from that of a disc blackbody (van Paradijs & Lewin 1986). In other words, fitting covariance spectra with a disc blackbody model such as diskbb in Xspec (Makishima et al. 1986) would lead to an overestimate of the intrinsic disc temperature. Therefore, following van Paradijs & Lewin (1986), within Xspec we defined the model deltadisk⁵ that represents the derivative of the diskbb model over T_in. This corresponds to computing the variable spectrum as a small perturbation of the time-averaged one. In this model, the temperature of the disc in covariance spectra T_in,covar was left free to vary, and the seed photons temperature of nthComp was tied to T_in,covar. The fits were performed in the 0.3–10 keV band, simultaneously for groups of consecutive observations belonging to the same phase of the outburst (as reported in Table 1 and Fig. 1), thus resulting in a total of four groups. Due to its short exposure, O168 did not allow us to constrain the parameters of the model well, so it was discarded. Results of the fits are reported in Table E.1 for each fitting group. As an example, the best-fit model to high-frequency (7− 30 Hz) covariance spectra of combined observations O104-O105 is shown in Fig. 5.

The TBabs[deltadisk + nthComp] model is statistically preferred to the simple Comptonisation model in the joint fits for each phase of the outburst (for the four fitting groups the χ² improves by 128, 167, 74, and 71, respectively, for a difference of 7, 18, 12, and 8 in the number of degrees of freedom). These results point to the presence of linearly correlated variability between the primary source of hard X-ray photons and the disc on the short timescales corresponding to the sampled frequencies. This is expected if the excess soft emission is produced via hard X-ray irradiation of the inner disc, and the process is linear at least to first order (Sect. 5). This is also in agreement with the interpretation of the soft X-ray lags detected in MAXI J1820+070 as due to thermal reprocessing (Sect. 4.1).

In Fig. 6 we report the best-fit values of the Γ_in,covar and kT_in,covar parameters as a function of spectral hardness. The slope of the hard X-ray Comptonisation component does not show significant changes throughout this part of the outburst (left panel of Fig. 6). This indicates that the spectral structure of the innermost parts of the hot plasma does not vary significantly throughout the hard and hard-intermediate states: it remains quite hard, although the errors become large towards the transition and a spectral softening cannot be excluded during the last observation.

On the other hand, the temperature of the disc black body component in high-frequency covariance spectra is observed to steadily vary with hardness (right panel of Fig. 6). We checked that letting the N_H free to vary in covariance spectra results in a slight, systematic shift of kT_in,covar towards lower values, which does not significantly influence the obtained results. The measured disc black body component in covariance spectra becomes hotter, by a factor of ~5, as the time-averaged spectrum softens. Notably, during the rise, plateau, and bright decline, the characteristic frequency ν₀ of the soft X-ray lag (and thus its intrinsic amplitude) is consistent with an increase (decrease in intrinsic lag amplitude) by a similar or slightly smaller factor (of ~3−5; Fig. 4, left panel). This strongly hints at a connection between the two in these phases of the outburst. Indeed, the steady increase in the characteristic frequency of soft X-ray lags (Fig. 4) in these phases of the outburst suggests a change in the geometry of the innermost accretion flow. This might ultimately lead to a change in the irradiation of the disc, causing (according to Eq. (1)) the observed increase in the inner disc temperature kT_in,covar (Fig. 6, right panel).

Therefore, we assumed that the scale height of the X-ray source is relatively small, and used the best-fit kT_{in, covar} values from high-frequency covariance spectra to put constraints on the inner radius of the irradiated disc via Eq. (2). The effective temperature of the disc is influenced by the variable hard X-ray flux, but intrinsic dissipation as well irradiation from constant (or variable on longer timescales) hard X-ray emission may also contribute to it. Thus, in general Eq. (2) yields a lower limit on the inner radius, unless heating from irradiation dominates in these regions (kT_in∕f_col = kT_eff ~ kT_eff,irr).

We assumed f_col = 1.7 (Shimura & Takahara 1995; Kubota et al. 1998; Gierliński et al. 1999; Davis et al. 2005). In order to estimate the total irradiating luminosity L_irr we performed systematic fits of time-averaged spectra in the energy band 3–10 keV. Details of these fits are reported in Appendix E. The Eddington-scaled irradiating luminosity l_irr was estimated from the 0.1–1000 keV flux (F_{irr,0.1-1000 keV} reported in Table E.1) of the extrapolated best-fit Comptonised incident spectrum and assuming the constraints on the distance from Atri et al. (2020). For the albedo we assumed a value of a = 0.5, which is an average of the values observed in the hard state of other BHXRBs (~ 0.3−0.7 Zdziarski & De Marco 2020). We accounted for the uncertainty on this parameter in the computations of the errors (see below). We kept the reflection fraction fixed at (as expected for an isotropic, static X-ray source, covering the inner disc), since time-averaged fits yield estimates very close to this value. In order to compute the 90 per cent confidence errors on the inferred values of R_in we randomised the measured irradiating fluxes from time-averaged spectra (according to a log-normal distribution, Uttley et al. 2005) and the best-fit values of kT_in,covar (according to a normal distribution) within their range of uncertainties. The error was then inferred from the resulting distribution of R_in values. In the computation of the errors we also included the uncertainty on the assumed values of f_col and a by considering the extreme values of these parameters within the ranges f_col ~ 1.3−2.2 (e.g. Salvesen & Miller 2021) and a ~ 0.3−0.7 (Zdziarski & De Marco 2020). Given that this computation yields a lower limit on R_in, we conservatively accounted for the uncertainties on the BH mass and source distance (Sect. 1) by considering their 1σ confidence values that minimise the inner radius (Sect. 1).

Figure 7 reports the results as a function of spectral hardness. These are consistent with a significant disc-truncation throughout most of the analysed observations, with a net decrease in the disc inner radius (by a factor of ~5–6 at most) between the first and the last observation.

Fig. 5 Unfolded best-fit TBabs[deltadisk + nthComp] model for the high-frequency (7− 30 Hz) covariance spectrum of combined observations O104-O105. The residuals observed at ~1–2 keV are likely due to NICER’s calibration systematics (NICER’s effective area shows features due to O, Si, and Au at 0.5, 1.8, and 2.2 keV).

Fig. 6 Best-fit values of spectral parameters from fits of high-frequency covariance spectra with the TBabs[deltadisk + nthComp] model. The parameters are plotted as a function of spectral hardness (i.e. the ratio between the 4–12 and the 2–4 keV total count rates). Left panel: spectral slope of the hard Comptonisation component, Γcovar. Right panel: inner disc temperature, kTin,covar.

Fig. 7 Lower limits on the inner disc radius as a function of spectral hardness, as obtained from constraints on the quasi-thermal component in high-frequency covariance spectra. The horizontal dotted lines mark the ISCO for a Schwarzschild and Kerr BH. Results are obtained assuming an isotropic static X-ray source irradiating the inner disc (reflection fraction ) and with a small scale height compared to the disc radius.

5 Discussion

Our analysis shows that the geometry of the innermost accretion flow of MAXI J1820+070 changes significantly and continuously throughout the hard and hard-intermediate states of its 2018 outburst. This result was obtained via two independent X-ray spectral-timing methods applied to NICER data, namely the analysis of soft X-ray reverberation lags, and of the soft quasi-thermal component in high-frequency covariance spectra.

5.1 The evolution of the inner flow geometry

During the rise, the plateau, and the bright decline (Fig. 1) the soft X-ray lags in MAXI J1820+070 show a steady increase (decrease) of characteristic frequency ν₀ (intrinsic amplitude) as the source softens (Fig. 4, left panel). This trend is in qualitative agreement with the evolution of the frequency of the type-C QPO in the PSDs of the source (Buisson et al. 2019). Soft X-ray lags in BHXRBs are ascribed to disc thermal reverberation of incident hard X-ray photons (e.g. Uttley et al. 2011; De Marco et al. 2015b; Kara et al. 2019). Therefore, the behaviour observed in MAXI J1820+070 can be interpreted in terms of a reduction of relative distance between the Comptonising region and the disc as the source softens (such behaviour is schematically illustrated in Fig. 8 in the framework of the evolving inner radius of a truncated disc; see discussion below). However, during the hard-soft transition this trend breaks (Fig. 4, right panel). In particular, during observations at spectral hardness ≲ 0.28, ν₀ suddenly decreases (by a factor of ~2–3), implying an increase in the reverberation lag intrinsic amplitude. We notice that this occurs simultaneously to an intense hard (15–50 keV)X-ray flare, as detected in Swift/BAT light curves (Wang et al. 2020b), suggesting a major change in the innermost emitting regions. Moreover, during O197, when the transition took place (as identified via X-ray and radio analysis; Shidatsu et al. 2019; Homan et al. 2020) there is no evidence of hard lags (the lack of a hard lag is observed at least down to ~ 0.05 Hz), and a long (~15–20 ms) reverberation lag is observed to dominate the entire frequency band (Fig. 3 bottom right panel). We suggest that such a behaviour might be the consequence of a different component dominating the X-ray emission at this time, possibly associated with the observed (Bright et al. 2020; Espinasse et al. 2020; Wood et al. 2021) ballistic jet ejections from the source (see Fig. 8 and discussion below). State-of-the-art X-ray spectral-timing models may not be able to reproduce the soft X-ray lags, as well as the (apparently total) suppression of hard X-ray lags in this phase of the outburst.

At the frequencies of the soft X-ray reverberation lag covariance spectra show excess soft X-ray emission, well-modelled by a quasi-thermal component due to the disc (Fig. 5, Sect. 4.2). This indicates that this component is produced via X-ray heating, and that the thermal response of the disc to variable hard X-ray irradiation preserves a certain level of coherence. We found that the inner temperature of the disc black body component in the high-frequency covariance spectra of MAXI J1820+070 steadily increases as the source softens (Fig. 6, right panel). Following the method outlined in Zdziarski & De Marco (2020) and assuming a constant and relatively small scale height for the X-ray source, we used these measurements to obtain constraints on the inner radius of the irradiated disc region responsible for the observed reverberation lags. Our estimates are consistent with an evolution of R_in (Fig. 7), causing changes in the irradiation of the disc, and consequently in its temperature (Fig. 6, right panel). In most of the hard state (i.e. during the rise, plateau, and bright decline) the observed increase (decrease) in the irradiated inner disc temperature (radius) as the source softens is consistent with the scaling characterising the frequency of soft X-ray reverberation lags (a change by a factor of ~ 3−4 is observed in all these cases, Figs. 4, 6 right panel, and 7). This behaviour is qualitatively in agreement with predictions of disc-truncated models (e.g. Esin et al. 1997; Done et al. 2007). However, the observed trend characterising the irradiated inner disc temperature and radius persists till transition, not showing any clear break as instead seen for the soft lags. In other words, these two independent diagnostics of geometry appear to stop probing the same regions when the source is close to transition.

Therefore, we interpret the behaviour near the transition to soft-intermediate states as due to a drastically different configuration of the innermost accretion flow than that characterising the preceding phases of the outburst. In particular, by tracing backthe motion of the discrete, relativistic ejecta from the central X-ray source detected in radio and X-ray observations of MAXI J1820+070, Bright et al. (2020) and Espinasse et al. (2020) estimate a launch date around MJD 58 306, while an additional earlier ejection (MJD 58 305.6) has been recently reported (Wood et al. 2021). The latter coincides with the last of our observations(O197). This suggests a link between these discrete ejections, the long reverberation lag and the suppression of hard X-ray lags.

The scenario we propose to explain the observed phenomenology is illustrated in the simplified sketch of Fig. 8. The inner radius of the disc decreases steadily throughout most of the hard state, as implied by the increase (decrease) of lag frequency (amplitude) and as traced by the inner temperature of the irradiated inner disc in high-frequency covariance spectra. Close to transition, while the inner radius is consistent with settling near (or at) the ISCO (as traced by the steady evolution of the quasi-thermal component in high-frequency covariance spectra), the structure of the X-ray source undergoes major changes (as traced by the break in the scaling of soft X-ray lags and the suppression of hard X-ray lags). Indeed, the intrinsic reverberation lag depends on the location of the region(s) of hard X-ray energy dissipation and the reprocessing region in the disc. If the emission due to internal shocks in a ballistic jet becomes significant close to the transition, then the dissipation of hard X-rays would predominantly occur in a larger or more distant region, which would cause the irradiation of a larger area of the disc. This would qualitatively explain the generally lower observed characteristic lag frequencies (Fig. 4), and in particular, the very long reverberation lag (Fig. 3 bottom rightmost panel) at state transition. Moreover, if the hard X-ray lags commonly observed during the outburst are produced in the inner accretion flow (e.g. as proposed in models of propagating accretion rate fluctuations; e.g. Lyubarskii 1997), a change in the dominant process of hard X-ray emission would also explain the observed suppression of these hard lags. In particular, this is expected if the dominant hard X-ray component is largely decoupled from fluctuations generated in the accretion flow (e.g. if this component originates in shocks along the jet, Vincentelli et al. 2018). We observe that the soft X-ray lags start deviating from the trend observed at higher spectral hardness about 4 days before the transition (Fig. 4, right panel), namely on MJD 58 302. This might indicate that the major changes in the inner accretion flow leading to the launch of relativistic ejecta start occurring a few days before the ejection. Interestingly, Buisson et al. (2021) reported the detection of a much longer (~20 s) soft X-ray lag ~0.4 day after the transition (i.e. in soft-intermediate states). The authors interpreted it as due to the disruption and refilling of the inner disc following the jet ejection.

We again emphasise that the estimates of the disc inner radius through Eq. (2) (Fig. 7) are strictly valid for a X-ray source with a height above the disc small compared to the disc radius (for example a corona above the inner disc or a lamppost located at a few R_g from the BH, Sect. 4.2). Moreover, our computations assume that the height of the source remains approximately constant. If this condition does not apply, then variations in the height of the X-ray source are expected to play a significant role too in the observed changes of disc illumination. Recently, Kara et al. (2019) and Buisson et al. (2019) proposed a scenario whereby the vertical extent of the X-ray source, rather than the inner edge of the accretion disc, decreases throughout the hard state, driving the observed changes in the thermal reverberation lags. In particular, this was proposed in order to justify the remarkable stability of the broad component of the Fe K line in the hard and hard-intermediate states of MAXI J1820+070. However, this scenario entails constant illumination of the inner parts of the disc, which ultimately implies that the geometry of the regions of the X-ray source closer to the disc (and responsible for the irradiation of the inner disc) does not change significantly. Our results appear to disfavour such a scenario because in this case we would expect to measure a constant inner disc temperature, while we find the temperature to increase steadily (Fig. 6, right panel). Therefore, our results appear not in agreement with the scenario proposed by Kara et al. (2019) and Buisson et al. (2019). On the other hand, given the above mentioned assumptions of our analysis, we cannot completely rule out the possibility that the X-ray source varies in its height (rather than in its vertical extent) during the outburst, thus contributing to the observed trend of inner disc temperature. However, this scenario seems hard to reconcile with the observed evolution of type-C QPOs in this source (Buisson et al. 2019). Therefore, we favour a solution whereby the hard state is dominated by a decrease in the inner radius. This is also supported by a recent re-analysis of four NuSTAR observations of MAXI J1820+070 at different hard state epochs (Zdziarski et al. 2021b), according towhich the lack of a spectral evolution of the broad component of the Fe K line is indicative of modest relativistic effects, as expected when the inner radius does not reach the ISCO. The apparent lack of a clear relativistic component in the high-frequency covariance spectra of the source (Fig. 5) further corroborates this interpretation.

Fig. 8 Sketch of the physical scenario proposed to explain the observed evolution of soft X-ray reverberation lags and of the quasi-thermal component from irradiation in high-frequency covariance spectra observed in MAXI J1820+070.

5.2 The truncation radius of the disc

Our analysis of high-frequency covariance spectra (Sect. 4.2) yields a lower limit on the disc inner radius for each of the analysed observations (Fig. 7). These estimates depend on the value of the reflection fraction parameter . In our treatment the parameter is used to account for a possible reduction of hard X-ray flux reaching the disc for causes other than changes of disc radius or X-ray source height (Zdziarski & De Marco 2020). Assuming isotropic illumination of the disc we found truncation at ≳30 R_g at the beginning of the monitoring, corresponding to bolometric luminosities of a few per cent of L_E (Zdziarski et al. 2021b). The disc turns out to be significantly truncated (≳ 10 R_g) throughout most of the bright hard and hard-intermediate states, and appears to reach close to the ISCO only at transition. Assuming a lower value of (e.g. as found in Zdziarski et al. 2021b) would reduce the inferred value for R_in by a factor , thus not affecting our inference of a disc significantly truncated throughout most of the analysed phases of the outburst. However, even accounting for the possible effects of a smaller , the data remain consistent with a disc at the ISCO of a low spinning BH during the transition (O197). This is in agreement with other independent results suggesting the BH in MAXI J1820+070 has spin a < 1 (Buisson et al. 2019; Atri et al. 2020).

As discussed in Sect. 4.1, the soft X-ray reverberation lags can also be used to infer an estimate of the disc truncation radius. However, this requires a more complicated modelling approach that relies on the simultaneous fit of variability, lag and time-averaged spectra (e.g. Mahmoud et al. 2019; Mastroserio et al. 2019), and depends on the assumptions on the geometry of the Comptonising region. This problem has been widely discussed in the literature (Uttley et al. 2011; De Marco et al. 2017; Wang et al. 2020a), and can be understood by considering the intrinsic soft lag amplitude implied by its characteristic frequency ν₀ (Fig. 4), and the distance that this amplitude maps. The corresponding Euclidean distance is , where θ is the angle between the line of sight to the X-ray source and the straight path linking the X-ray source to the disc reprocessing region. The case θ = π∕2 approximates a truncated disc illuminated by a central source (in this case the derived distance coincides with the inner disc radius). Assuming θ = π∕2 and i = 70°, our estimates of ν₀ (Sect. 4.1 and Table E.1) yield Euclidean distances between the X-ray source and the reprocessing region in the disc ranging between a few tens and a few hundreds of R_g. These distances are a factor of ~3−10 times larger than the disc inner radius inferred from covariance spectra. In fact the distance inferred from reverberation lags is a weighted average among all possible light paths between the X-ray source and the disc (i.e. corresponding to the mean lag of the impulse response function; e.g. Wilkins & Fabian 2013). Therefore, in case of complicated geometries (i.e. more intricate impulse response function), such as for an extended X-ray source, the link between the measured lag and the geometrical parameters of the disc and of the X-ray source may be quite complex. This has been also shown in Mahmoud et al. (2019), where the inner truncation radius resulted significantly smaller (by a factor of ~4) than the distance mapped by the observed reverberation lag. This would explain why the soft X-ray lags measured in MAXI J1820+070 apparently imply distances larger than the disc inner radius inferred from covariance spectra.

It is interesting to note that such large discrepancies between the distances implied by X-ray reverberation lags and thedisc inner radius inferred by other independent methods is generally not observed in active galactic nuclei (AGN; e.g. Fabian et al. 2009; De Marco et al. 2013; Marinucci et al. 2014; Kara et al. 2016; Alston et al. 2020). This might be due to a substantial difference between the geometry of the X-ray source in hard state BHXRBs ascompared to that of those AGN (all type-1 non-jetted Seyferts or narrow line Seyferts) for which X-ray reverberation has been measured. In particular, dissipation of hard X-rays in these AGN likely occurs in a more compact region (e.g. Fabian et al. 2015), thus reducing the number of possible light paths and thus the width of the impulse response function. Such inference is also in agreement with the fact that the characteristic timescales of variability appear systematically offset towards lower frequencies in hard state BHXRBs as compared to AGN (e.g. Körding et al. 2007).

Finally, it is important to point out that in the previous discussion we make the implicit assumption that the detected reverberation lags are dominated by the light travel distance between the X-ray source and the disc, while other timescales are negligible. In particular, the thermalisation timescale should also contribute to the observed delays. However, this quantity scales inversely with the density and directly with the temperature of the disc surface layers irradiated by the X-ray source (e.g. Collin et al. 2003). For certain values of these two parameters (e.g. density ≲ 10¹⁸cm⁻³ for temperatures of ~1 keV) the thermalisation timescale may be of the order of the observed lags (i.e. ~ 0.1−1 ms). However, while the density of the disc is currently a not well-constrained parameter, recent studies suggest it to be quite high in BHXRBs (≳ 10¹⁹⁻²⁰cm⁻³, Tomsick et al. 2018; Jiang et al. 2019; Zdziarski & De Marco 2020). If so, the corresponding reprocessing time is negligible compared to the measured lag. Moreover, contrary to what observed, if the lags were dominated by the thermalisation timescale, the observed process would be highly non-linear (incoherent; e.g. Vaughan & Nowak 1997).

6 Conclusions

We performed a systematic X-ray spectral-timing study of the first part of the 2018 outburst of the BHXRB system MAXI J1820+070. Our analysis covers the hard and hard-intermediate states as well as the transition to the soft-intermediate states. Our main results are as follows:

• The soft X-ray thermal reverberation lags show a steady evolution throughout most of the hard state (specifically during the rise, plateau, and bright decline). This evolution is indicative of a decrease in the relative distance between the X-ray source and the disc as the source softens.

• The quasi-thermal component associated with the soft X-ray reverberation lags and observed in high-frequency covariance spectra shows a steady increase in the inner temperature. This is a consequence of changes in the hard X-ray irradiation of the inner disc and is consistent with a decrease in the disc truncation radius, in agreement with the observed evolution of soft X-ray lags. The disc is consistent with being highly truncated at the beginning of the outburst (R_in ≳ 30 R_g) and with remaining significantly truncated throughout most of the hard and hard-intermediates states (R_in ≳ 10 R_g).

• Major changes occur before the transition to the soft-intermediate states. The thermal reverberation lags suddenly increase (decrease) in amplitude (frequency), breaking the trend observed during the preceding phases of the outburst. This suggests an increase in the relative distance between the X-ray source and the disc. At the same time, the hard X-ray lags intrinsic to the primary continuum get suppressed. On the other hand, the observed trend characterising the quasi-thermal component in high-frequency covariance spectra persists until transition, not showing any clear break. This trend is consistent with the inner edge of the disc increasingly approaching the ISCO, and possibly settling there at (or close to) transition.

In order to explain the observed X-ray spectral-timing phenomenology, we propose a scenario whereby the inner radius of the disc decreases steadily throughout the hard and hard-intermediate states. This scenario is an alternative to the contracting vertical extent of the corona proposed in Kara et al. (2019) and appears more in line with the detected steady increase in the irradiated disc temperature. While the inner disc maintains this trend close to the transition to soft-intermediate states, reaching (or nearing) the ISCO, the X-ray source undergoes major structural changes. In particular, we propose that in this phase the dissipation of hard X-rays predominantly occurs in the ballistic jet (possibly through internal shocks). This would qualitatively explain the behaviour of soft X-ray reverberation lags and the suppression of hard X-ray lags (provided the latter are intrinsically produced in the accretion flow). Interestingly, the soft X-ray lags start deviating from the trend observed during most of the hard state about 4 days before the appearance of a type-B QPO (which marks the hard-to-soft transition; Homan et al. 2020) and the ejection of relativistic plasma (as later detected at radio and X-ray wavelengths; Bright et al. 2020; Espinasse et al. 2020), possibly representing the precursor to such events. Future intensive monitoring campaigns of BHXRB systems with NICER will allow us to further test this scenario.

Acknowledgements

This work is part of the EU-funded “BHmapping” project. While this work was being reviewed, Wang et al. (2021) was published, also showing the behaviour of reverberation lags at transition. The authors thank G. Mastroserio for providing the version of RELTRANS used in this paper, M. Mendez for helpful discussions on fits of variability spectra, G.K. Jaisawal for help with NICER data reduction, and the anonymous referee for their useful comments. B.D.M. acknowledges support from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 798726 and Ramón y Cajal Fellowship RYC2018-025950-I. A.A.Z. and M.D. acknowledge support from the Polish National Science Centre under the grants 2015/18/A/ST9/00746 and 2019/35/B/ST9/03944. G.P. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 865637). T.M.B. acknowledges financial contribution from the agreement ASI-INAF n.2017-14-H.0 and PRIN-INAF 2019 n.15. We acknowledge support from the International Space Science Institute (Bern).

Appendix A Stationarity

The condition of stationarity was verified by measuring the PSD in the energy band 2–11 keV. We extracted PSDs from light curve segments of 40 s, and measured an average PSD for each observation, or part of observation characterised by a constant number of FPMs (the assumption here is that the source is stationary during the time spanned by each dataset). The light curves have a time bin of 0.4 ms (Sect. 2), and thus each PSD covers the frequency interval from 0.025 Hz up to a Nyquist frequency of 1250 Hz. This allowed us to obtain a good sampling of both the broadband X-ray variability components and the Poisson noise contribution. The latter was fit at high frequencies (ν > 300 Hz) with a constant model and subtracted out. The PSDs were normalised using fractional rms-squared units (Belloni & Hasinger 1990; Miyamoto et al. 1992). Only observations visually showing compatible PSD shape (i.e. characteristic frequencies and normalisation) were combined for our X-ray spectral-timing analysis (i.e. an average cross-spectrum was computed from which the lag and covariance spectrum was estimated, Sect. 4). Examplesof PSDs showing the non-stationary behaviour of the source are shown in Fig. A.1.

Fig. A.1 Selection of PSDs of MAXI J1820+070 during different phases (colour-coded according to Fig. 1 and reported in the labels) of the first part of the 2018 outburst (two examples are shown for each phase). The PSDs are extracted inthe energy band 2–11 keV and reveal a strong and complex evolution of the X-ray timing properties of the source.

Appendix B Observed soft X-ray lag amplitude and its evolution as a function of spectral hardness

The computations of frequency-dependent lags reported in Sect. 4.1 were double-checked using three different codes for spectral-timing analysis. Two of the codes were independently implemented by us in the past (e.g. Belloni et al. 2005; De Marco et al. 2015b, 2017; Belloni et al. 2020). In addition, we processed the data using the public software for spectral-timing analysis Stingray⁶ (Huppenkothen et al. 2019a,b). In all the cases the observed amplitude of the soft X-ray lags results a factor of ~ 3 longer than reported in Kara et al. (2019). The existence of such a discrepancy was verified using the same energy bands (i.e. 0.5–1 keV vs. 1–10 keV) as in Kara et al. (2019). As an example, in Fig. B.1 we contrast the lag-frequency spectrum of observation O106 obtained from the three different codes and the spectrum reported in Kara et al. (2019). We notice that, while the interpretation given in Kara et al. (2019) of a contracting corona in the hard state does not depend on the observed amplitude of the soft X-ray lags, the longer lags might not be easily reconciled with their conclusion of a disc extending close to the ISCO.

For further comparison with Kara et al. (2019), in Fig. B.2 we also report the evolution of the maximum observed soft X-ray lag amplitude (from the best-fit RELTRANS models) as a function of spectral hardness. This amplitude does not show the same steady trend as inferred from the lag characteristic frequency ν₀ (Fig. 4 and Sect. 4.1). The rise is characterised by an increase in the observed maximum amplitude of the lag. Then, during the plateau this amplitude decreases (as previously reported in Kara et al. 2019 for selected observations). During the last part of the plateau, the entire bright decline and the beginning of the hard-soft transition the observed lag amplitude does not have a well-defined trend, and it is consistent with being constant. Finally, it increases during the last observations of the hard-soft transition phase (Fig. B.2, right panel). We argue that this behaviour and its apparent inconsistency with the evolution of ν₀ is a consequenceof the concurrent complex evolution of the hard X-ray lags, which are likely to have a strong influence on the observed maximum amplitude of the soft lag (Fig. 2). It is worth noting that, contrary to what affirmed in Kara et al. (2019), the observed amplitude of the longest soft lags detected in MAXI J1820+070 (mostly at the beginning of the outburst, Fig. B.2), are consistent with that of the soft lag detected in GX 339-4 during its bright hard state (De Marco et al. 2015b, 2017).

Fig. B.1 Lag-frequency spectrum of observation O106 between energy bands 0.5–1 keV and 1–10 keV independently computed with the code used in this paper (implemented in IDL and used e.g. in De Marco et al. 2015b, 2017), with the GeneralHigh-Energy Aperiodic Timing Software (GHATS; used e.g. in Belloni et al. 2005), and with Stingray (developed by Huppenkothen et al. 2019a,b and publicly available). The scatter in the data is due to slight differences in the extraction and binning. These lags are compared to those reported in Kara et al. (2019) for the same observation.

Fig. B.2 Evolution of observed maximum soft lag amplitude (as inferred from the best-fit models) in lag-frequency spectra as a function of spectral hardness. For comparison with Fig. 4, the left panel reports results from the rise, plateau, andbright decline phases, while the right panel reports results from the hard-soft transition phase as compared to the other phases (grey symbols).

Appendix C Independence of ν₀ on the hard X-ray lags of the continuum

In Sect. 4.1 we used the first frequency at which the reverberation lag drops to zero, ν₀, as an unbiased diagnostic of the intrinsic reverberation lag amplitude. Here we show that the hard X-ray lags intrinsic to the X-ray continuum do not affect this frequency. For simplicity, we assume that one band, usually called channel-of-interest (CI, as in Uttley et al. 2014) contains both reprocessed and continuum photons, while the other band, usually called the reference band (REF), contains only continuum photons. This assumption is particularly appropriate for the energy bands used in our analysis, where the soft band is chosen to have contribution from the thermal component but includes contribution from the Comptonisation continuum as well, while the hard band is dominated by the primary Comptonised emission and has negligible contributionfrom the thermal/reflection component. It can be easily shown that in this case, the observed phase lag is: (C.1)

where ω = 2πν, Δτ_revω is the phase shift due to reverberation between reprocessed photons in the CI and continuum photons in the REF band; Δτ_contω is the phase shift of continuum photons in the CI with respect to continuum photons in the REF band; and f_R is the ratio between reprocessed and continuum photons in the CI. According to Eq. C.1, the observed phase lag is zero, ϕ_obs = 0, when (C.2)

The hard X-ray lags intrinsic to the continuum show a steadily decreasing trend with frequency (Fig. 2). This behaviour is commonly observed in BHXRBs (e.g. Nowak et al. 1999; Pottschmidt et al. 2000; Uttley et al. 2011; De Marco et al. 2015b), and predicted by models (e.g. Arévalo & Uttley 2006). Given the tendency of hard X-ray lags to get suppressed at high frequencies, it is reasonable to assume that their contribution is negligible (Δτ_cont ~ 0) around ν₀. In this caseEq. C.2 shows that the observed phase lag is zero when Δτ_revω = nπ, with n =1, 2.., so that the first zero-crossing frequency can be used to infer the intrinsic amplitude of the reverberation lag (see also Mizumoto et al. 2018). However, we verified that ν₀ is not significantly biased by the presence of hard lags even when the assumption that Δτ_cont ~ 0 around ν₀ does not strictly hold. To this aim we fixed the value of the reverberation lag at Δτ_rev = −0.01s (so that the first zero-crossing frequency in the absence of continuum hard lags is ν₀ = 50 Hz) and solved Eq. sectionlinking C.2 for three different values of Δτ_cont (namely Δτ_cont = 10⁻³s, 10⁻⁴s, and 10⁻⁵s; the highest value of 10⁻³s is chosen so as to resemble the value typically measured at the frequencies where the reverberation lag starts to dominate, Fig. 2). The results are shown in Fig. C.1 (where the zero-crossing points are solutions of Eq. C.2; here f_R was fixed at 0.5). It is clear that for Δτ_cont ≲ 10⁻³s the shift induced on ν₀ by the presence of hard X-ray lags is negligible. The maximum shift (about 8 Hz) of ν₀ compared to the Δτ_cont = 0 case is obtained for Δτ_cont = 10⁻³s; however, this shift is smaller than the estimated uncertainties on ν₀ (Table E.1).

Interestingly, observations of MAXI J1820+070 show that the amplitude of the hard lags at a given fixed (high) frequency increases as the source softens (Fig. 2). This behaviour is generally observed in BHXRBs (e.g. Kylafis & Reig 2018). However, according to Fig. C.1, the longer the hard lags, the larger would be the shift induced on ν₀ towards lower frequencies. This implies that, if driven by the observed changes in the hard lags, we would expect to observe a shift of ν₀ towards lower frequencies as the source softens. This is at odds with one of our main results (Fig. 4). Thus we conclude that the hard lags in MAXI J1820+070 do not affect our estimates of ν₀ and do not influence its observed increasing trend with decreasing spectral hardness.

Fig. C.1 Solutions of Eq. C.2 for three different values of Δτcont (and fixed Δτrev = −0.01s and fR = 0.5). The first high-frequency zero- crossing point is the solution that corresponds to ν0.

Appendix D Fits of lag-frequency spectra

Fits of lag-frequency spectra were performed using the RELTRANS spectral-timing model (Ingram et al. 2019) in the version presented inMastroserio et al. (2019; 2020), which also allows for the fit of hard lags. We let the X-ray source height, h, and the disc inner radius, R_in free to vary. We also let the parameters regulating the spectral pivoting of the inverse-Compton spectrum (Mastroserio et al. 2018)free to vary in order to account for contributions from the hard lags in the fitted frequency band. We let the BH mass vary within its estimated range of uncertainties (Sect. 1). We fixed the BH spin at its maximum value (in order to allow for the largest range of values for the disc inner radius, and thus the highest values of ν₀), the column density at N_H ~ 1.4 × 10²¹cm^-2 (Kajava et al. 2019; Dziełak et al. 2021), the boosting factor and the normalisation at 1. We also assumed Γ =1.6, log(ξ∕erg cm s^-1) ~ 3.3 and E_cut = 30 keV. The latter parameters were fixed in order to speed-up the fits. Nonetheless, letting them free to vary does not affect our estimates of ν₀ as this parameter, for a given value of BH mass, depends only on the height of the X-ray source and on the inner radius of the disc. The errors on ν₀ were obtained considering the 90 per cent confidence interval of the h and R_in parameters within Xspec.

Appendix E Estimating the total irradiating flux and the thermally re-emitted fraction

We used relxillCp (v1.4.3; Dauser et al. 2016 and references therein) to model the broadband continuum in time-averaged spectra. The spectrum is absorbed by a column of cold gas (modelled with TBabs) with column density fixed to N_H ~ 1.4 × 10²¹cm^-2 (Kajava et al. 2019; Dziełak et al. 2021). The emissivity index was fixed at 3, the spin parameter was fixed at its maximum value a = 0.998, and the high energy cut off, as done for the fit of covariance spectra, was fixed at kT_e = 30 keV (Zdziarski et al. 2021b). The spectral index of the Comptonised incident spectrum, the ionisation parameter, the reflection fraction, and the disc inner radius were left free to vary. Finally, the iron abundance and inclination angle were determined from the simultaneous fit of observations O103, O137, O170, and O197 (chosen so as to sample different phases of the outburst). We obtained a best-fit value of for the inclination and an upper limit of A_Fe < 0.51 for the iron abundance. Thus we fixed these parameters to i = 71° and A_Fe = 0.5 in the individual fits of the single observations. The fits of the single observations yield high values of reflection fraction consistent with . The resulting 0.1–1000 keV irradiating fluxes are reported in Table E.1.

While time-averaged spectra yield the total irradiating flux (which contributes to the observed inner disc temperature), high-frequency covariance spectra yield the true thermal response of the inner disc to hard X-ray irradiation, with negligible contribution from disc internal dissipation. This follows from the fact that while disc emission due to internal dissipation may be variable on the fastest timescales, this emission is more easily suppressed near the region where it is produced in a geometrically thin optically thick accretion flow (e.g. Churazov et al. 2001; Mushtukov et al. 2018), thus not contributing significantly to the overall high-frequency variability. Moreover, in order to contribute to covariance spectra, the variable intrinsic disc emission should be linearly correlated with the observed primary hard X-ray continuum, most likely through Compton up-scattering (e.g. Nowak & Vaughan 1996) or downward propagation of mass accretion rate perturbations (e.g. Lyubarskii 1997; Kotov et al. 2001; Bollimpalli et al. 2020). However, both processes would produce hard X-ray lags, whereas in our case the lag-frequency spectrum is dominated by soft X-ray lags at high frequencies (Fig. 3 and Sect. 4.1). Therefore, we expect any contribution from disc internal dissipation producing variable soft X-ray emission and a coherent hard X-ray response to be negligible in high-frequency covariance spectra.

In order to measure the variable thermally reprocessed flux emitted by the disc we can fit high-frequency covariance spectra with the model thcomp(deltadisk), where the use of the convolution model thcomp (Zdziarski et al. 2020) allows us to account for the contribution from both the up-scattered and direct disc photons. For example, applying this model to the high-frequency covariance spectrum of O103, we find a thermally reprocessed-to-total variable flux fraction of about 10 per cent. Such a fraction is significantly lower than expected if all of the observed flux in the variable primary continuum was reprocessed in the disc (either reflected or thermalised, and assuming an albedo of ~ 0.3–0.7). This is an interesting result that does not have an immediate explanation, thus deserving future investigations. Nonetheless,we argue that a more complex physical situation than that of a single power law emitting-region may determine the observed high-frequency covariance spectra. In this event, the fraction of Comptonised photons that effectively irradiates the disc should be lower, with negligible contribution from other (possibly inner) regions of the accretion flow (e.g. as in the scenarios proposed in Zdziarski et al. 2021a and Dziełak et al. 2021).

References

All Tables

All Figures

Fig. 1 HID of MAXI J1820+070 during its 2018 outburst (from MJD 58 189 to MJD 58 305). The hardness is computed as the ratio between the 4–12 keV and 2–4 keV NICER count rates. On the y axis we report the total 2–12 keV NICER count rate re-normalised for the number of used FPMs. Highlighted in different colours and symbols are the phases of the outburst that are the focus of our work. The grey empty circles show the subsequent spectral evolution of the source up until MJD 58 404. The epochs of the appearance of a type-B QPO (MJD ~ 58 305.7; Homan et al. 2020) and the launch of relativistic discrete radio ejecta (MJD 58 305.6 ± 0.04; Wood et al. 2021, and ~MJD 58 306, Bright et al. 2020; Espinasse et al. 2020) are marked, respectively, by the bigger star and orange arrows.

In the text

Fig. 2 0.5–1 keV vs. 2–5 keV lag-frequency spectrum of some of the analysed observations of MAXI J1820+070 during the rise, plateau, bright decline, and hard-soft transition phases (with colour codes and symbols as in Fig. 1). The positive lags dominating at low frequencies (≲2 Hz) are hard X-ray lags. The plot illustrates the hump-like structures characterising these low-frequency hard lags and their complex evolution throughout the first part of the outburst. At high frequencies (≳2 Hz) the (negative) soft X-ray lag becomes dominant (see inset, with the y axis in linear scale).

In the text

	Fig. 3 0.5–1 vs. 2–5 keV lag-frequency spectra of selected observations of MAXI J1820+070. The observations illustrate the diverse behaviour of soft (negative) lags during different phases of the outburst (as also reported in the insets). The solid red curve is the best-fit RELTRANS model. The blue arrow marks the frequency range used for the extraction of covariance spectra.
In the text

Fig. 4 Evolution of ν0 in lag-frequency spectra as a function of spectral hardness. Left panel: results from the rise, plateau, and bright decline phases, which are observed to follow the same trend. The dotted black line is the best-fit model. Right panel: results from the hard-soft transition phase as compared to the other phases (grey symbols), which is observed to break from the observed trend at higher spectral hardness.The arrows are 90 per cent lower and upper limits.

In the text

	Fig. 5 Unfolded best-fit TBabs[deltadisk + nthComp] model for the high-frequency (7− 30 Hz) covariance spectrum of combined observations O104-O105. The residuals observed at ~1–2 keV are likely due to NICER’s calibration systematics (NICER’s effective area shows features due to O, Si, and Au at 0.5, 1.8, and 2.2 keV).
In the text

	Fig. 6 Best-fit values of spectral parameters from fits of high-frequency covariance spectra with the TBabs[deltadisk + nthComp] model. The parameters are plotted as a function of spectral hardness (i.e. the ratio between the 4–12 and the 2–4 keV total count rates). Left panel: spectral slope of the hard Comptonisation component, Γcovar. Right panel: inner disc temperature, kTin,covar.
In the text

Fig. 7 Lower limits on the inner disc radius as a function of spectral hardness, as obtained from constraints on the quasi-thermal component in high-frequency covariance spectra. The horizontal dotted lines mark the ISCO for a Schwarzschild and Kerr BH. Results are obtained assuming an isotropic static X-ray source irradiating the inner disc (reflection fraction ) and with a small scale height compared to the disc radius.

In the text

	Fig. 8 Sketch of the physical scenario proposed to explain the observed evolution of soft X-ray reverberation lags and of the quasi-thermal component from irradiation in high-frequency covariance spectra observed in MAXI J1820+070.
In the text

	Fig. A.1 Selection of PSDs of MAXI J1820+070 during different phases (colour-coded according to Fig. 1 and reported in the labels) of the first part of the 2018 outburst (two examples are shown for each phase). The PSDs are extracted inthe energy band 2–11 keV and reveal a strong and complex evolution of the X-ray timing properties of the source.
In the text

Fig. B.1 Lag-frequency spectrum of observation O106 between energy bands 0.5–1 keV and 1–10 keV independently computed with the code used in this paper (implemented in IDL and used e.g. in De Marco et al. 2015b, 2017), with the GeneralHigh-Energy Aperiodic Timing Software (GHATS; used e.g. in Belloni et al. 2005), and with Stingray (developed by Huppenkothen et al. 2019a,b and publicly available). The scatter in the data is due to slight differences in the extraction and binning. These lags are compared to those reported in Kara et al. (2019) for the same observation.

In the text

	Fig. B.2 Evolution of observed maximum soft lag amplitude (as inferred from the best-fit models) in lag-frequency spectra as a function of spectral hardness. For comparison with Fig. 4, the left panel reports results from the rise, plateau, andbright decline phases, while the right panel reports results from the hard-soft transition phase as compared to the other phases (grey symbols).
In the text

	Fig. C.1 Solutions of Eq. C.2 for three different values of Δτcont (and fixed Δτrev = −0.01s and fR = 0.5). The first high-frequency zero- crossing point is the solution that corresponds to ν0.
In the text