A full spectral-timing model to map the accretion flow in black hole binaries: the low/hard state of MAXI J1820+070 Free

ABSTRACT

1 INTRODUCTION

Spectral studies of the low/hard state in black hole binaries lead to conflicting conclusions about the nature and geometry of the accretion flow. The hard Comptonised continuum spectra seen in this state cannot arise from an isotropically emitting corona over the inner accretion disc (Poutanen, Veledina & Zdziarski 2018 and references therein). A truncated disc/hot flow geometry is most often proposed to solve the energetic constraints, with the truncation radius decreasing for softer spectra as the source increases in brightness towards the transition to the high/soft (disc dominated) state (see e.g. Done, Gierliński & Kubota 2007). However, this is at odds with more detailed spectroscopy around the iron line, where the line profile appears so broad as to require an untruncated disc (e.g. García et al. 2015, 2018; Xu et al. 2018). The additional information encoded in the fast (0.01–|$100\, {\rm s}$|⁠) large amplitude variability which is characteristic of the low/hard state can potentially break the spectral controversy by giving an independent set of diagnostics into the nature and geometry of the state (e.g. Uttley et al. 2014).

One key advantage of the truncated disc geometry is that it also gives a framework to explain the overall variability properties. A decreasing truncation radius between the disc and hot flow as the spectrum softens in the transition from the low/hard (Compton dominated) to high/soft (disc dominated) spectra means a decreasing characteristic time-scale for all variability. This is seen in the obvious decrease in time-scale picked out by both the low-frequency break in the broad-band fast variability power spectrum and the prominent low-frequency (0.1–|$10\, {\rm Hz}$|⁠) quasi-periodic oscillation (QPO) seen in the Comptonised emission (Wijnands & van der Klis 1999). Both these can be quantitatively explained in a model where the inner radius of the disc is the outer radius of the hot flow, with the low-frequency break in the broad-band variability set by fluctuations in mass accretion rate stirred up on the viscous time-scale of the outer flow, while the QPO is Lense–Thirring (relativistic vertical) precession of the entire hot flow (Ingram, Done & Fragile 2009; Ingram & Done 2011, 2012).

There are other diagnostics of the fast time-scale variability which again can be most easily interpreted in the truncated disc/hot inner flow geometry. Combined spectral and timing constraints show that fluctuations in higher energy bands generally lag behind the same fluctuations seen in lower energy bands, and that the size of this lag decreases as the time-scale of the fluctuations decreases (Miyamoto & Kitamoto 1989; Revnivtsev, Gilfanov & Churazov 1999; Wilkinson & Uttley 2009; Grinberg et al. 2014; Uttley et al. 2014). Compton scattering should imprint a lag as a function of energy due to the additional light travel time required for each successive scattering order, but this is much shorter than observed lags and has no dependence on fluctuation time-scale (Miyamoto & Kitamoto 1989; Nowak, Wilms & Dove 1999). Instead, the observed lags are now generally interpreted as being related to the propagation time-scale of fluctuations through the accretion flow. Slow fluctuations are produced at the larger radii in the flow and have to propagate down through the entire flow to reach the innermost radii, so have the longest time lags. Conversely, faster fluctuations are produced at smaller radii, so have a shorter distance to travel to reach the innermost regions, and hence shorter time lag (Lyubarskii 1997; Kotov, Churazov & Gilfanov 2001; Arévalo & Uttley 2006). Propagation lags do not give any observable signature if the accretion flow is homogeneous, emitting the same spectrum at each radius. The additional key assumption required to explain the propagation lags in the data is that the flow is also radially stratified in its energy spectrum, with softer spectra emitted at larger radii, and harder spectra emitted closer to the black hole. This is a natural consequence of the truncated disc/hot inner flow geometry, as parts of the flow which are closest to the disc will see more seed photons for Compton cooling, so the spectrum will be softer than that close to the black hole, where fewer seed photons from the disc penetrate and where the energy released from gravity is highest (Axelsson & Done 2018).

Radial stratification of the Comptonised emission also gives a possible solution to the dilemma of the iron line profile. A sum of softer and harder Comptonisation components gives a spectrum which is subtly curved. Fitting a multitemperature Comptonisation continuum with a single temperature model leaves broad residuals. The iron line modelling assumes that all these residuals arise from reflection on the disc, but while the ionised reflection models do contain curvature, this is accompanied by strong atomic line features. Only by dramatically smearing the lines by extreme relativistic effects can these be smoothed into a pseudo-continuum to match the curvature in the data (Makishima et al. 2008, Kolehmainen, Done & Díaz Trigo 2014, Basak et al. 2017, Zdziarski et al. 2021a). The derived extreme gravity is then an artefact of continuum complexity.

However, the iron line also gives an additional way to determine the inner disc edge. It is produced by the (fluctuating) hard X-rays illuminating the disc, so its emission must be lagged by the light traveltime. This is difficult to see in the data as the iron line is a fairly small feature, and is produced at energies where standard CCD X-ray detectors are not so sensitive. However, X-ray illumination of the disc also produces a reflected continuum (Compton hump, peaking at 30 keV) and (quasi)thermalized emission from the X-ray heated disc photosphere, peaking below 1 keV. The latter can easily be studied with CCD instruments as long as the interstellar column to the source is low, and ‘soft lags’ – so called as the disc reprocessed emission is at soft energies and lags the hard X-ray continuum variability – are now clearly seen (Uttley et al. 2011). Their evolution to shorter time-scales as the source makes the transition from the low/hard to high/soft state supports the truncated disc models (De Marco et al. 2015, 2017, 2021).

Despite the overall successes of the truncated disc/hot inner flow geometry, developing fully quantitative models for all these spectral-timing features has proved difficult. On the spectral fitting alone, a complex continuum is not quite enough to make the truncated disc reflection models be the best statistical fit to the data (Buisson et al. 2019; Zdziarski et al. 2021b). Similarly, it is difficult to fully model the power spectra and the lags as a function of both frequency and energy (Rapisarda et al. 2016; Mahmoud & Done 2018a; Mahmoud, Done & De Marco 2019).

Here, we use NICER and NuSTAR data from the extremely bright recent outburst of MAXI J1820+070 (Kawamuro et al. 2018; Tucker et al. 2018) in its low/hard state. The inclination of the source has been measured as 66–81° from the orbit (Torres et al. 2020) and as 64 ± 5° from the jet (Wood et al. 2021). The source is relatively nearby (distance |$\sim 3.0\, {\rm kpc}$|⁠; Atri et al. 2020), and has only small interstellar absorption along its line of sight (galactic absorption |$\sim 10 ^{21}\, {\rm cm}^{-2}$|⁠; Uttley et al. 2018) so the spectrum can be studied down to the softest energies. NICER can handle the high count rates required for fast timing studies better than any other current detector, so the combination of object and instrument gives the best data so far on the disc in this controversial state.

We develop a spectral-timing model including both the propagating mass accretion rate fluctuations and reverberation through a full truncated disc/radially stratified hot inner flow geometry. Propagation and reverberation together play a key role in the fast variability, and various models combining them have been developed and applied to both AGN (Wilkins & Fabian 2013; Gardner & Done 2014; Wilkins et al. 2016) and black hole binaries (Rapisarda, Ingram & van der Klis 2017; Mahmoud et al. 2019). The most self-consistent way of spectral-timing modelling is to simultaneously fit both the spectral and Fourier spectral properties. However, this is extremely expensive in CPU time (Mahmoud et al. 2019). Instead we build on earlier techniques which fix the spectral components from detailed fits to the time-averaged energy spectra, and then couple these to a timing model of the variability which includes both propagation and reverberation (e.g. Gardner & Done 2014; Mahmoud & Done 2018a). This approach gives a quantitative model which by definition fits all the major features of the energy spectra. Here, we build a variability model using the energy spectral components and show that this can also broadly fit the power spectra as a function of energy and the lags seen between different energy bands as a function of frequency, including the switch between propagation (soft lead) and reverberation (soft lag).

While the model still has some limitations, we show that a truncated disc at a few tens of gravitational radii reproduces the overall features of the power spectra and lag-frequency spectra. We also show how the propagated variability allows us to map out the accretion speed as a function of radius in the hot flow. This depends on the efficiency of the underlying angular momentum transport processes, which depends on the magnetic field geometry and strength. Our results favour models where the transport is dominated by a small-scale dynamo process (the magnetorotational instability: MRI; Balbus & Hawley 1991, 1998) giving Standard And Normal Evolution (SANE) rather than models with maximal magnetic flux (Magnetically Arrested Disc: MAD), where the dynamo is suppressed and accretion occurs via instabilities in the large-scale magnetic fields (Narayan et al. 2012). Thus spectral-timing models hold out the possibility of observationally testing models linking the jet to the hot flow.

2 OBSERVATION AND DATA REDUCTION

We study the bright low/hard state of the black hole binary MAXI J1820+070. We choose to use data from 2018-03-21, close to the peak in high energy flux, where there are simultaneous NICER and NuSTAR data (Buisson et al. 2019; Kara et al. 2019; Zdziarski et al. 2021b). These data are summarized in Table 1.

The NICER data were processed using the NICER data analysis software, NICERDAS 2019-06-19_V006a in HEASoft version 6.26, and CALDB 20200722. The unfiltered data were cleaned with the standard calibration and standard time screening through nicercal and nimaketime (Stevens et al. 2018). Then the data were barycentre-corrected with the tool barycorr. Following reported caveats by the NICER team, we excluded the occasionally noisy focal plane modules 14 and 34 from the created cleaned event file.

The energy spectrum and light curves were extracted with XSELECT. We excluded periods of high background by discarding episodes where the 13–|$15\, {\rm keV}$| light curve exceeds |$2\, {\rm counts\,s^{ -1}}$| (Ludlam et al. 2018). The power spectra, cross-spectra, and time lags as a function of frequency were calculated from light curves in segments of 65.536 s with 1 ms time bins (2¹⁶ points) and rebinned logarithmically with the constant factor of 1.2 by following the prescriptions described in Uttley et al. (2014) and Ingram (2019). Any segments with data gaps were discarded, which yields the smaller net exposure for timing analysis of |$2294\, {\rm s}$|⁠.

The NuSTAR data were reduced with the NuSTARDAS pipeline v.1.8.0 and CALDB 20200912. To filter passages through the South Atlantic Anomaly, we set saamode = strict and tentacle = yes. Following reported caveats by the NuSTAR team, we used statusexpr='STATUS= = b0000xxx00xxxx000’ to remove spurious triggers caused by bright sources. The source and background regions were extracted with the circle of 60 arcsec radius, in which the former was centred on the peak brightness and the latter was placed in the lowest apparent source contamination. The source is so bright in this observation that the background is negligible. However, the NuSTAR count rate is not high enough to give strong constraints on the high-frequency timing behaviour which we focus on here, so we use these data only to better define the spectral model. We use only data from a single detector module, FPMA, to illustrate the high-energy spectral shape and to avoid the cross-calibration between the FPMA and FPMB modules.

3 SPECTRAL-TIMING MODELS

A full spectral-timing model requires that we specify both the spectrum and its variability as a function of radius, with the propagation of fluctuations and reverberation through this radial structure.

Previous attempts at this have used a variety of different assumptions tailored to the observed bandpass. Ingram & Done (2012) modelled the timing properties assuming the spectrum was a homogeneous hot flow (single Comptonisation), with emissivity assumed to be a power law as a function of radius. The variability in this flow was assumed to be generated and propagated on the local viscous time-scale, which depends on radius. This picture can be linked to models of angular momentum transport in the flow. The MRI in particular causes fluctuations in all physical quantities, and fluctuations in density propagate inwards on the local viscous time-scale, producing fluctuations in mass accretion rate, which affects the observed luminosity (Lyubarskii 1997; Kotov et al. 2001; Arévalo & Uttley 2006). This model could fit to RXTE data, where the low-energy bandpass limit of >3 keV meant that the disc was not important spectrally. However, propagated disc variability is likely important in the hot flow even when not in the observed bandpass (Uttley et al. 2011). Not including a distinct low-frequency variability component from the disc in the hot flow variability means the model struggles to reproduce the breadth of the power spectrum, requiring uncomfortably large intrinsic variability power and spectral emissivity at very small radii (Ingram & Done 2011, 2012; Ingram & van der Klis 2013). This remains a problem even when the Comptonisation region is radially stratified, with softer Comptonisation continua produced at larger radii (Mahmoud & Done 2018a).

Instead, Rapisarda et al. (2016) show that including fluctuations from the disc can relieve these tensions, where the disc time-scale for fluctuations is discontinuously lower than the viscous time-scales in the hot flow at the disc truncation radius, consistent with physical expectations of the lower scale height disc structure compared to the hot flow. Importantly, this naturally gives the double-humped shape of the power spectra, which are otherwise very difficult to reproduce (Mahmoud & Done 2018b). Including lower frequency variability from the disc also means that the hot flow needs only to produce the higher frequency variability, which reduces the emissivity and variability requirements at small radii which are otherwise quite stringent (Mahmoud & Done 2018b).

However, the Rapisarda et al. (2016) models did not include full spectral information for the hot flow. Instead, they assumed the hot flow emissivity was a power law as a function of radius, and fit for the emissivity in each energy band. This allowed them to fit power spectra as a function of energy, and to calculate the lags between two energy bands, but does not use all the information contained in the energy spectrum.

Mahmoud & Done (2018a) used the energy spectrum to define the intrinsic spectral components (disc, soft, and hard Comptonisation), together with the reflected spectra. The intrinsic spectra were assumed to be strictly radially stratified in order to couple these with a variability model (generation and propagation of fluctuations) to make a full spectral-timing description of the data. They extended the formalism in Mahmoud & Done (2018b) and Mahmoud et al. (2019), but the bumpiness of their derived power spectra was predominantly due to the enhancement of variability/emissivity at specific radii.

Here, we revisit the approach of Mahmoud & Done (2018a), fitting for the spectral components and then using these to define a full spectral-timing model. But we couple this to the timing model incorporating the insight of Rapisarda et al. (2016) where the disc time-scale is discontinuously lower than that in the hot flow.

4 FITTING THE SPECTRAL COMPONENTS

The spectral fits were performed with the X-ray spectral fitting package XSPEC (Arnaud 1996). We use the energy ranges 0.5–|$10\, {\rm keV}$| and 3–|$78\, {\rm keV}$| for the NICER and NuSTAR data, respectively. Following the recommendation by the NICER team to handle calibration uncertainties,¹ we add 1 per cent systematic errors to the NICER data. In addition, we allow the relative normalization of NICER with respect to NuSTAR to be a free parameter, but assume the same spectral shape for both instruments.

Previous detailed spectral analyses of the NuSTAR data have shown that fits with a thermal disc and hot Comptonisation have some tensions, even when including two reflection components of different ionization and typical radii. Both Zdziarski et al. (2021b) and Buisson et al. (2019) show that this modelling leads to best-fitting solutions where one of the reflectors is strongly smeared, requiring R_in ∼ 2R_g, completely inconsistent with a truncated disc. Furthermore, there are other tensions in this fit, with the reflectors having surprisingly high Fe overabundance of Z_Fe ∼ 6 and surprisingly low inclination of ∼30°. Adding a second continuum component relieves these tensions, with iron overabundance reduced to Z_Fe ∼ 2 for the measured binary and jet inclination of ∼60–65°. Each continuum component is associated with its own reflection component (now with different shapes due to the differences in the spectral index as well as the different ionization). With this spectral complexity, the derived inner disc radius increases to ∼100 R_g, consistent with a truncated disc. None the less, the derived solution is still in the highly relativistic regime as the softer continuum lamppost source has a height of ∼2 R_g (Zdziarski et al. 2021b).

We use the spectral model of Zdziarski et al. (2021b) (our data correspond to their epoch 2) with their two Comptonisation continua and reflection components as our baseline, i.e. we use the compps model continuum and corresponding reflkerrG ionised reflection models. These are required when fitting to the NuSTAR data as the standard nthcomp thermal Comptonisation models in XSPEC do not give the proper shape of the high-energy rollover (Niedźwiecki, Szanecki & Zdziarski 2019; Zdziarski et al. 2020). This impacts the reflected continuum as well, so the xillverCp reflection models also underestimate the high-energy continuum.

We fit first to the NuSTAR data alone and confirm that our fit matches that of Zdziarski et al. (2021b). We then extend the model down to the NICER energy range and see clear evidence for excess thermal emission from a disc component, as shown in Fig. 1 (a). We include this with a diskbb model, a simple sum of blackbody components, characterized by a peak temperature from the inner edge of the disc, kT_bb.

Figure 1.

Sequence of fits to NICER (magenta) and NuSTAR (orange) spectra. The upper and lower panels show the unabsorbed spectra and data-to-model ratios, respectively. The energy spectra are normalized to the NuSTAR effective area. (a) The NuSTAR spectrum alone is fitted with two Comptonisation components and reflection of these components from the truncated disc. The black solid line shows the total spectrum. The hard (soft) Comptonisation S_hc(E) (S_sc(E)) and its reflection S_hr(E) (S_sr(E)) are shown with blue (green) solid and dashed lines. (b) Both the NICER and NuSTAR spectra are fitted with the addition of a disc component S_d(E) (red). The unphysical low-energy cutoff of the Comptonisation and reflection components are corrected. A narrow Gaussian component centred at 6.4 keV (cyan) is included. In addition, three Gaussian components around 1.0, 1.9, and 2.4 keV (grey, dotted) are included to compensate for residuals. The derived parameters are summarized in Table. 2. The shaded energy bands, low (red), mid (green), and high (blue), are mainly used in the timing analysis.

Open in new tabDownload slide

However, the reflkerrG reflection model, which is used to define the direct Comptonisation spectrum as well as the reflected spectrum, fixes the seed photon of the direct Comptonisation spectrum at a blackbody temperature of |$1\, {\rm eV}$|⁠.² In addition, the reflection models do not include the seed photon temperature of the incident Comptonisation as a free parameter, and the reflection spectra are computed for a fixed seed photon temperature equivalent to a disc blackbody temperature of |$50\, {\rm eV}$|⁠.³ As a result, whereas the seed photons at |$kT_{\mathrm{bb}}\sim 0.22\, {\rm keV}$| are obtained from the spectral fitting, both the direct and reflected emission extend to lower energies than the seed photons in an unphysical way. This likely has an impact on attempts to measure the disc density from the low-energy extent of the reflected emission (e.g. Jiang et al. 2019).

In order to correct for the seed photon break, we make a local XSPEC model, xslowe, which multiplies both the reflected and direct Comptonisation spectra by the difference between the low-energy Comptonisation rollover for any given seed photon temperature and that assumed in the relxill and reflkerr models^{2, 3} (see also García et al. 2014, 2016). However, we stress that this is only an approximate correction factor. The reflection models of García et al. (2016) conserve energy, so the energy removed by suppressing the unphysical low-energy extent must emerge somewhere, most likely peaking around 3–4kT_bb like a blackbody or disc blackbody, i.e. around 1 keV. This lack of reprocessed emission in our best-fitting model becomes important when we consider the reverberation signal (see Section 7.3).

The residuals are dominated by a clear atomic feature at 1 keV from a combination of Ne X and Fe L shell transitions. An ∼5 per cent feature is seen at 1 keV in the best calculations of the disc photosphere at similar temperatures to those here of 0.22 keV (Davis & Hubeny 2006). Not all elements are ionized at these temperatures, so full radiative transfer through solar abundance material is required to calculate these features. However, these models are tabulated including relativistic corrections, assuming that the disc extends down to the innermost stable circular orbit. Thus they can be used to fit disc-dominated states (e.g. LMC X-3: Kubota et al. 2010), but the relativistic corrections are not appropriate for a truncated disc as assumed here. Hence, we cannot use these photosphere models to directly fit to our data, and therefore just add a Gaussian component at ∼1 keV.

The residuals then show a smaller set of features at 1.5–|$2.5\, {\rm keV}$|⁠, probably from highly ionized Si and S. These are not produced significantly in the photosphere models, so are more likely explained by the limitations in current reflection models discussed above. Again, we compensate for these using Gaussian components at 1.9 and |$2.4\, {\rm keV}$|⁠, respectively. It is important to model these features as otherwise the higher statistics of the NICER data pull the fit away from the NuSTAR.

There is then also a small residual around the iron line, so we additionally include a narrow Gaussian component to model this emission. This results in Δχ² = 152 for the addition of two free parameters, so is highly significant. The resulting line EW is ∼8 eV, similar to that seen in a similar bright low/hard state of GX339-4 (Done & Diaz Trigo 2010). We suggest that this is from the strong optical wind which is seen in this state (Muñoz-Darias et al. 2019; Sai et al. 2021), and stress that this should be included in all detailed models of the iron line profile (see discussion in Axelsson & Veledina 2021).

In summary, our employed model is written as tbabs*(diskbb+xslowe*reflkerrG(SC)+xslowe*reflkerrG(HC)+gauss+gauss+gauss + gauss). The unabsorbed data and model components, all normalized to the NuSTAR effective area, are shown in Fig. 1 (b). The galactic absorption is taken into account with tbabs. The blackbody component emitted from the truncated disc is calculated by diskbb. The Comptonisation components and relevant reflection components shown with solid and dashed lines in green (soft: SC) or blue (hard: HC) colour are calculated with xslowe*reflkerrG. The three Gaussian components to correct for residuals below 3 keV are shown in grey dotted lines, while the cyan solid line compensates for the additional narrow core iron line emission at |$6.4\, {\rm keV}$|⁠. The model parameters are tabulated in Table 2. The figure also shows the three energy bands which we will use for the power spectra, low (red, 0.5–|$1.5\, {\rm keV}$|⁠), mid (green, 1.5–|$5.0\, {\rm keV}$|⁠) and high (blue, 5.0–|$10.0\, {\rm keV}$|⁠).

Clearly (and model independently), the low-energy band has mainly disc emission with some Comptonisation, while the high-energy band has mainly Comptonisation and very little disc emission. Thus the low- and high-energy band power spectra should be distinctly different as they are dominated by quite different components. In our particular spectral model, the soft and hard Comptonisation components contribute about equally to the mid and high energy bands, which means their power spectra should be quite similar. We note that other spectral decompositions (e.g. a soft Comptonisation component which is more like a hotter blackbody component, as derived e.g. by Dziełak, De Marco & Zdziarski 2021 for these data or Di Salvo et al. 2001 for Cygnus X-1) would predict instead that the mid and high energy bands were dominated by different components so could have distinct variability.

We stress that none of the variability results presented in Section 5 depend on the detailed shape of the model components. It becomes important only in Section 7 where the detailed timing properties are computed using this spectral decomposition.

5 VARIABILITY DATA

5.1 Power spectra

Variability on a wide range of time-scales can be examined by power spectra. We calculate them from the NICER data, subtracting the Poissonian noise contribution (Uttley et al. 2014) and adopting the normalization such that the integral of the power spectra corresponds to the fractional variance (Belloni & Hasinger 1990; Miyamoto et al. 1991; van der Klis 1997; Vaughan et al. 2003). The power spectra in three energy bands, low (0.5–|$1.5\, {\rm keV}$|⁠, red), mid (1.5–|$5.0\, {\rm keV}$|⁠, green), and high (5.0–|$10.0\, {\rm keV}$|⁠, blue), are shown in the upper panel of Fig. 2. It is seen that they can be characterized by two broad-band variability components. We refer to the slower and faster components as P₁ and P₂, respectively.

Figure 2.

Upper panel: power spectra in three energy bands calculated from the NICER data. The red smooth lines show the fit to the power with the sum (solid) of two broad-band variability components P₁ and P₂ (dashed). Lower panel: ratio of the power spectra of the mid (green) and high (blue) energy bands compared to the low-energy band. In both panels, the vertical error bars represent 1σ.

Open in new tabDownload slide

The lower panel of Fig. 2 highlights the much stronger high-frequency variability at higher energies by showing the ratio of the power spectra in the mid and high energy bands with respect to that in the low energy band. However, it also shows a very small increase in variability at low frequencies, around 0.03–0.04 Hz. This is primarily due to a weak QPO with fundamental around |$0.036\, {\rm Hz}$| seen at this time, together with its second harmonic at |$0.072\, {\rm Hz}$| (Buisson et al. 2019; Ma et al. 2021). The observed low-frequency variability in any energy band is then the sum of the broad-band noise P₁ and the QPOs. The QPOs are predominantly seen in the Comptonisation components so increase the ratio of the mid/high to low-energy power spectra at the QPO frequencies. Without the contribution from the QPO, we see that the underlying low-frequency component P₁ is similar in shape and normalization in all energy bands, peaking at |$0.08\, {\rm Hz}$|⁠.

By contrast, the energy-dependent behaviour of P₂ appears to be intrinsic. It has long been known that the higher energy bands have more variability in short time-scales, this is generally modelled in data above |$3\, {\rm keV}$| as having the same shape but higher normalization. The current data are the first to enable us to deeply explore the behaviour of P₂ at lower energy bands, where the thermal emission from the disc dominates the energy spectrum. Including the low-energy band makes it clear that the frequency at which the P₂ peak occurs is changing as well as its normalization. We study the energy dependence of P₂ in more detail below.

5.2 Modelling the power spectra to characterize the shift in peak frequency of P₂ with energy

We now use power spectra over narrower energy bands to look at the behaviour of P₂ in more detail. We fit these power spectra with the sum of the two broad-band variability components P₁, P₂. Note that neglecting the weak QPOs does not affect the estimates of the peak frequency of P₂.

To obtain the energy dependence of the peak frequency, we perform the fit to the power spectrum in each energy band with the sum of two variability components, P₁(f) + P₂(f). For P₁, we fix the peak frequency, f_c1, and width, Δf₁, as we find that these are almost constant for different energy bands. The free parameters left are thus A₁, A₂, Δf₂, and β₂. We show the energy dependence of the peak frequency as well as index β₂ governing the asymmetry of P₂ is shown in Fig. 3, where the errors are derived by a Monte Carlo technique (Press et al. 1992) from 1000 simulated power spectra.

Figure 3.

Energy dependence of the peak frequency (upper) and the index of the asymmetric Lorentzian β₂ (lower) for the higher frequency variability component P₂. The vertical error bars represent 68 per cent confidence intervals derived from the Monte Carlo technique. The data in 7.5–|$10.0\, {\rm keV}$| band are removed due to large uncertainties.

Open in new tabDownload slide

It is clear that the majority of the shift in peak frequency of P₂ with energy occurs below |$2\, {\rm keV}$|⁠. This corresponds to the energy at which the spectrum makes a transition from being dominated by the thermal disc to the Comptonisation (see Fig. 1). Thus this shift in peak frequency should be related to the transition between the disc and hot flow. By contrast, β₂ shows a more constant decline with energy, but crosses unity (standard zero centred Lorentian) at around 2 keV. At higher energies β₂ is less than unity, meaning there is more high-frequency power than in a standard zero-centred Lorentian. At lower energies, the index β₂ is greater than unity, so as well as the decrease in peak frequency of P₂, there is a sharper drop in power spectrum at frequencies beyond the peak (see low-energy power spectrum, red in Fig. 2). This is the first time that these properties have been seen in such detail. They show a shift in behaviour at the energy where the spectrum shows the transition between emission from the disc and hot flow, so they give us a new way to explore what is happening in this controversial region.

5.3 Lag-frequency spectrum

Another important aspect of the variability is time lags as a function of frequency as these show the causal connection between different spectral components. The lag-frequency spectrum between the soft band (0.5–|$1.5\, {\rm keV}$|⁠) and a total hard band (2.0–|$10.0\, {\rm keV}$|⁠, extended down to 2 keV to maximize counts) is shown in Fig. 4. The soft bands are primarily related to the disc while the hard is dominated by Comptonisation. Positive time lags mean variations in hard photons lag behind those in soft photons, whereas negative time lags mean variations in soft photons lag behind those in hard photons.

Figure 4.

Lag-frequency spectra between the soft band 0.5–|$1.5\, {\rm keV}$| and hard band 2.0–|$10.0\, {\rm keV}$|⁠. The positive time lag means fluctuations in the hard band lag behind those in the soft band and vice versa. The symmetric logarithmic scale is employed for the representation of the time lags, i.e. linear between |$-10^{-3}$| and |$10^{-3}\, {\rm s}$| and logarithmic otherwise. The vertical error bars represent 1σ. The grey dashed line shows the time lags of zero.

Open in new tabDownload slide

At longer time-scales, the hard band lags the soft band. On the other hand, the causal connection switches around |$\sim 3\, {\rm Hz}$|⁠, giving rise to the soft band lagging the hard band at shorter time-scales. Such low-frequency hard lags and high-frequency soft lags are typically observed in black hole binaries (Uttley et al. 2011; Kara et al. 2019). They are explained by propagating mass accretion rate fluctuations to give the hard lags (e.g. Kotov et al. 2001) and thermal reverberation of the Comptonised emission illuminating the disc making the soft lags. In the following sections, we quantitatively model the lag-frequency spectra as well as power spectra by combining propagation and reverberation through the radially stratified spectral components.

6 VARIABILITY MODEL FORMALISM

Our model builds on the analytical scheme proposed by Ingram & van der Klis 2013, in which the power spectrum and cross spectrum of the local mass accretion rate are calculated first and converted into those of the flux. Thus, we mainly focus on aspects of prescription that differ from those in Ingram & van der Klis 2013. The full technical method to calculate the power spectra and time lags of the flux is described in Appendix  A, but we outline the various steps below.

6.1 Accretion geometry and spectral stratification

We assume that there are three spectral components, a thermal component from the disc, S_d(E), together with soft, S_sc(E), and hard, S_hc(E), Comptonisation components from the hot flow. Modelling the hot flow with two Comptonisation components comes from the spectral decomposition performed in Section 4. It is also supported by the fast variability since the hard time lags typically observed above 2 keV are difficult to explain with a single Comptonisation component (e.g. Mahmoud & Done 2018a for Cygnus X-1). We set the time-averaged shapes of these from the results of the spectral fit (Fig. 1b, Table 2) All these S are photon flux spectra, since the important quantity for variability is the number of photons.

We also assume that these spectral components are radially stratified, where the variable disc region extends from r_out down to r_ds, the soft Comptonisation component from r_ds to r_sh, and the hard Comptonisation component from r_sh to r_in, as shown in Fig. 5. In this geometry, we refer to the truncation radius as the outer radius of the variable disc r_out. All radii expressed with the small letter, r, are given in units of the gravitational radius R_g = GM/c². We fix |$M=8\, \mathrm{ M}_{\odot }$| (Torres et al. 2019).

Figure 5.

Assumed geometry in the model. The emission regions of the three spectral components, that is the thermal component S_d(E) and two Comptonisation components S_sc(E), S_hc(E), are radially stratified. Each region is spectrally homogeneous. The variable flow is located between r_in and r_out, whereas the thin disc located further does not contribute to the variability. The local fluctuations propagate inwards throughout the variable flow with the speed of the radial velocity v_r(r), which is determined by the radius from the central object, r, and viscous frequency f_visc(r).

Open in new tabDownload slide

6.2 Mass accretion rate fluctuations

We divide the variable flow between r_in and r_out into N_r rings logarithmically such that Δr_n/r_n is constant, where r_n and Δr_n (n = 1, ⋅⋅⋅, N_r from outer to inner rings) are the central radius of the nth ring measured from the central object and the distance between central radii of neighbouring rings. The number of rings N_r affects the resolution of the model.

6.3 Variability of the flux

7 MODELLING THE VARIABILITY

7.1 Matching the behaviour of P₂

We first explore whether the propagating fluctuations scenario can indeed match the observed shift in power spectral shape of P₂ with energy. We do not attempt to reproduce P₁ in this subsection.

Ingram & Done (2011) set the parameters B = 0.03 and m = 0.5 from the relation between the QPO frequency (assumed to be from Lense–Thirring precession) and low-frequency break in the broad-band variability (assuming this is set by fluctuations propagating from the outer edge of the hot flow). Here, instead we allow B and m in the hot flow to be free parameters as we do not model the low-frequency break (which is from P₁) in this subsection.

We assume r_in = 6 and set r_out = 45 for consistency with the next subsection where we discuss the origin of P₁ in more detail. The power spectra derived from a model with B = 45, m = 1.65 are shown in the upper panel of Fig. 6. Other model parameters are F_var = 0.55, N_r = 100, |$\Delta t=5\, {\rm ms}$|⁠, N = 2¹⁴, where Δt is the sampling interval and N the number of sampling points. The peak frequency and index β derived from the fit to the power spectra with the asymmetric Lorentzian (2) are shown in the lower panels of Fig. 6.

Figure 6.

Comparison of power spectra (a) and the energy dependence of the peak frequency (b) and asymmetry (c) between the model (solid) and observation (stepped). Energy bands in the power spectra are the same as in Fig. 2. The continuous viscous frequency throughout the variable flow is assumed in the model. Since the power spectra of low-energy bands calculated from the model are bumpy at high frequencies, making the derivation of the index β somewhat arbitrary, the indices at low-energy bands are omitted.

Open in new tabDownload slide

The energy-dependent power spectrum of P₂ as well as its peak frequency and index β are overall well reproduced. Both the shift in peak frequency and the change in the asymmetry with respect to energy are produced by the different proportion of the three spectral components, i.e. the emission of higher energies comes from the regions closer to the central object, which have variability in higher frequencies. The model particularly succeeds in reproducing the observed behaviour of P₂ above |$2\, {\rm keV}$|⁠, where the emission from the hot inner flow dominates the total energy spectrum. Thus, the propagating fluctuations scenario combined with the inhomogeneous hot flow explains well the behaviour of P₂ above |$2\, {\rm keV}$|⁠, though it has a sharper decrease below 2 keV than seen in the data.

However, the much larger discrepancy is that the current picture completely ignores the lower frequency variability which makes P₁. This is because the adjusted viscous frequency is fast compared to the prescription of Ingram & Done (2011) (B = 0.03, m = 0.5). If instead we change B and m to reproduce P₁, it is not at all easy to reproduce P₂ even with an extremely strongly centrally peaked emissivity and small inner radius (Ingram & Done 2011; Mahmoud & Done 2018a). Motivated by the work of Rapisarda et al. (2016), we then introduce a discontinuous jump in viscous frequency to reproduce both P₁ and P₂ simultaneously.

7.2 Discontinuity of the viscous time-scale between the disc and flow

The power spectra calculated with B_f = 6.0, m_f = 1.2, F_var = 0.8, and the same parameter values as in Fig. 6 are shown in Fig. 7. The double-humped shape of the power spectra matches the observation, where P₁ and P₂ originate from the disc and hot flow, respectively. The inwards propagation of fluctuations transmits P₁ into the mid (green) and high (blue) energy bands, where the emission from the hot flow dominates the energy spectrum.

Figure 7.

Power spectra obtained with a discontinuous viscous frequency prescription between the variable disc and hot flow. Lines and colours are the same as in Fig. 6 (a).

Open in new tabDownload slide

The clear defect of the model is that the low-energy band (red) lacks variability in the 0.3–|$3\, {\rm Hz}$| range. This is because it is dominated by the disc which only varies on the long time-scales of P₁. There is a very small contribution from the soft and hard Comptonisation regions in this energy band, so there is a small contribution to higher frequency power from P₂, but nowhere near enough to explain the data. We conclude that the power spectrum at energies where the disc emission is dominant cannot be explained by the propagating fluctuations process alone.

7.3 Including reflection and reprocessing

The hot flow illuminates the disc, so the fast variability of the Comptonisation components is reprocessed/reflected by the disc. The process is required not only by the lack of faster variability component P₂ in the low-energy band in the power spectrum, but by the detection of the soft lag, which shows that some part of the low-energy band emission lags behind the Comptonisation components, at high frequencies (see Fig. 4).

Using only the reflected emission from the spectral fits, S_sr(E) and S_hr(E), does not give sufficient signal for (negative – soft lags hard) reverberation lags to ever dominate over the (positive – hard lags soft) propagation lags. This shows that there must be a larger contribution of reprocessing than is included in the current spectral model. This is similar to the situation seen in the Narrow Line Seyfert 1 galaxy, PG 1211+026 (Gardner & Done 2014). Here, the solution was to realize that the spectral models of ionized reflection are not always appropriate, especially where the seed photons for Comptonisation are not close to the assumed (hardwired) temperature of |$50\, {\rm eV}$| (see Section 4).

There are four main free parameters which affect the power spectra and lag-frequency spectra, namely the average reverberation time lag, t₀, the fraction of reprocessed disc flux, f_rep, and the radial dependence of the propagation speed in the hot flow, set by B_f and m_f. The combination of these four parameters which best encompasses all of the data is shown in Fig. 8. How the model changes with each parameter is shown in Appendix  D.

Figure 8.

Power spectra (upper) and lag-frequency spectrum (lower) obtained with a discontinuous viscous frequency prescription and implementation of the reflected emission. Stepped and smooth lines show the observation and model, respectively. In the upper panel, colours are the same as in Fig. 6. In the lower panel, the energy bands are the same as in Fig. 4. For comparison, the lag-frequency spectra calculated with |$t_{0}=2.8\, {\rm ms}$| (magenta) and |$t_{0}=6.1\, {\rm ms}$| (orange) are shown. Larger time lags fit better to the maximum reverberation frequency, but exacerbate the discrepancy below 3 Hz.

Open in new tabDownload slide

This best overall model reproduces both the power spectra as a function of energy and the lag-frequency spectrum fairly well. The reprocessed emission not only produces more variability at high frequencies in the low-energy band as required, but also forms the soft lag at high frequencies. The fraction of the reprocessed emission is set to f_rep = 0.4, which yields r_sh = 17.8 and r_ds = 32.1 for r_in = 6 and r_out = 45. The impulse response is set by |$t_{0}=4.3\, {\rm ms}$|⁠, |$\Delta t_{0}=5.0\, {\rm ms}$|⁠, where the reprocessing starts at |$t_{0}-\Delta t_{0}/2 = 1.8\, {\rm ms}$| corresponding to the light crossing of 45R_g for M = 8 M_⊙. We also show the effect of changing t₀ on the lag-frequency spectra in Fig. 8, with |$t_{0}=2.8\, {\rm ms}$| (magenta) and |$6.1\, {\rm ms}$| (orange). Plainly the observed switch between propagation and reverberation at ∼3 Hz actually favours longer time lags rather than shorter ones, but these exacerbate the discrepancy below 3 Hz. All other parameter values are the same as in Fig. 7.

We check the consistency of the derived impulse response with the derived accretion flow geometry because the top-hat impulse response is one of the simplest assumptions of the reflection/reprocessing. In reality, the impulse response is affected by many factors, such as the geometry of the illuminating source to the disc and the inclination angle of the accretion disc to the observer. The derived light crossing of 45R_g does not mean the disc truncation of 45R_g, even if it is assumed that the reflection/reprocessing happens at the truncated disc.

Fig. 9 shows the comparison of impulse responses (a), power spectra in the low-energy band (b), and lag frequency spectra (c) for different formulations of the impulse response. The black lines represent the case of the top-hat impulse response adopted above. The red and blue solid lines represent the more realistic treatments of the lamp-post and extended illuminating source geometries, respectively. In the lamp-post geometry, a source height of h = 10 in units of R_g is assumed, while in the extended geometry, we consider a cylinder with the inner radius of r_in = 6, the outer radius of r_out = 45, and the height of h = 36. The reflection/reprocessing is assumed to happen at the outer thin disc from r_out = 45 to r_disc = 400, and contributions from further radii to the flux are ignored due to the small illuminating flux from the primary source (Mahmoud et al. 2019; Chainakun et al. 2021). The magenta dashed line is calculated for the same lamp-post source geometry as in the red line (h = 6), but the thin disc is assumed to be untruncated (r_out = 6), where the reflection/reprocessing occurs from r_out = 6 to r_disc = 400. The inclination angle of 66° is used, as in the spectral fitting. Details of the calculation of the more realistic impulse responses are described in Appendix  E. Note that the model is not self-consistent for the red solid line and magenta dashed line as it assumes both the radially extended hot flow and the disc truncation of 45R_g in calculating the propagating fluctuations part.

Figure 9.

Comparison of impulse responses h(t) (a), power spectra in the low-energy band of 0.5–|$1.5\, {\rm keV}$| (b), and lag-frequency spectra between 0.5–|$1.5\, {\rm keV}$| and 2.0–|$10.0\, {\rm keV}$| (c) for different formulations of the impulse response for reflection and reprocessing. Parameter values other than t₀ and Δt₀ are common (see the bottom row of Table 3). The black solid lines represent the top-hat impulse response. The red and blue solid lines are obtained with lamp-post and extended illuminating source geometries, respectively. The outer thin disc is assumed to be truncated in common. The magenta dashed line is calculated with the lamp-post source geometry and untruncated thin disc. In (a), the impulse responses are normalized such that the integral over time in units of second becomes unity. In (b) and (c), the observation results are shown for comparison.

Open in new tabDownload slide

The top-hat impulse response and that for the extended source geometry give similar power spectra and lag-frequency spectra, which demonstrates that the simple top-hat impulse response characterized by |$(t_{0}, \Delta t_0)=(4.3\, {\rm ms}, 5.0\, {\rm ms})$| is a reasonable approximation to more realistic impulse responses for truncated radii of 45R_g. The mean reverberation delay time for the best fit top-hat is 4.3 ms, whereas that for the self-consistent extended source/truncated disc with r_out = 45 is very similar at 4.9 ms. This is predominantly set by r_out rather than the emissivity as a lamp-post at h = 10 gives a very similar mean delay of 4.7 ms for r_out = 45, whereas the delay reduces to 1.8 ms for a lamp-post at h = 10 irradiating a disc with r_out = 6. The untruncated disc then gives shorter soft lags, so these dominate over the propagation lag at higher frequencies than for the truncated disc models. Thus, it is more difficult to explain the observed soft lag with the untruncated disc geometry in our model.

However, there are still some limitations to our model. It underestimates the propagation lags in the 0.3–|$3\, {\rm Hz}$| frequency range. More fundamentally, the power spectra do not reproduce the shift in peak frequency of P₂ with energy since the reprocessed disc variability now matches the higher energy power spectra rather than being shaped by propagation. These issues are discussed further in Section 8.2.

8 DISCUSSION

Our model produces the double-humped power spectrum by assuming there are intrinsic low-frequency fluctuations on the inner edge of the truncated disc, with a distinct jump to higher frequency fluctuations in the hot flow. The turbulent disc region may be quite limited in radial extent, making only a single Lorentzian, P₁, whereas the hot flow spans a larger range in radii, so the intrinsic fluctuations produce Lorentzians of different frequencies. Propagation of the fluctuations means that the power spectrum at any radius carries the imprint of variability generated at all larger radii. Thus at the inner edge of the disc, the power spectrum is only P₁, whereas at any radius, r, in the hot flow it is the sum over all radii in the hot flow of the much higher frequency Lorentzians at each radius, plus P₁ from the truncation region.

We couple this radial stratification of variability with radial stratification of the spectral components: the disc is (a sum of) blackbodies, whereas the Comptonisation in the hot flow is softer closer to the disc and harder closer to the central object. The combination of the radial stratification predicts that the blackbody has power spectrum P₁, the soft Comptonisation has power spectrum which is P₁ plus the sum of Lorentzians from the outer radius of the hot flow, r_ds down to r_sh, and the hard Comptonisation power spectrum is this plus the additional noise generated from r_sh–r_in. This picture produces the observed asymmetric shape (broader than a single Lorentizan) of the second hump in the power spectrum, P₂, and its increasing high-frequency extent with increasing energy due to the increasing contribution from the hard Comptonisation. The hot flow illuminates the disc, producing a fast variable reprocessed thermal component which adds to the intrinsic, slowly variable disc emission.

The power spectrum in any energy band is then determined by the contribution of each spectral component in each energy band. Similarly, the lag between correlated variability at a given frequency in any two energy bands will also be determined by a combination of propagation and reverberation lags in each spectral component and the changing contribution of the different spectral components with energy. Our model has a light traveltime to the reverberating disc of ∼45R_g, as required in the truncated disc models and consistent with the inferred intrinsic thermal component. This combined spectral-timing model fits the time-averaged energy spectrum (by construction) but can also reproduce most of the features of the power spectra as a function of energy and the lags as a function of frequency.

Our timing model can be used to test hot flow models in terms of the radial velocity (accretion speed) v_r(r) or equivalently, the viscous time-scale (accretion time-scale) t_visc(r) (see equation 6). We compare our results to various hot flow models, but stress that our timing model probably does not reproduce the observed properties sufficiently well to reach a robust conclusion.

8.1 Linking the viscous time-scale to accretion flow models

The comparison of viscous frequency as a function of radius derived from the timing model and that of various hot flow models are shown in Fig. 10. Generically, these have radial velocity v_r = α(H/R)²v_ϕ(R), where α, R, H, and v_ϕ(R) are the viscosity parameter, radius, scale height, and rotation velocity of the flow, respectively (Kato, Fukue & Mineshige 2008). Thus the purely analytic, self-similar Advection Dominated Accretion Flow (ADAF) solutions of Narayan & Yi (1994), which have constant α and H/R, predict v_r ∝ v_ϕ ∝ r^−1/2 (orange dash–dotted line), whereas our derived model has a much steeper dependence. However, including the stress-free inner boundary condition in a global ADAF model (orange dashed line) recover a good match for α ∼ 0.01–0.03 (Narayan, Kato & Honma 1997; Popham & Gammie 1998) as the flow must accelerate inwards.

Figure 10.

Viscous frequency as a function of radius from the central object. Whereas the left axis shows it in units of Hz, the right axis in units of c/R_g. The black hole mass of M = 8 M_⊙ is assumed. The solid lines show the measure from our model (see Fig. 8), and red, green, and blue colours correspond, respectively, to the disc, soft Comptonisation, and hard Comptonisation regions. The orange lines show the numerical (dashed) and self similar (dash–dotted) solutions from the ADAF model with α = 0.01 (Narayan et al. 1997). The prediction from the SANE and MAD models (chunk S6 and M5 in Narayan et al. 2012, respectively) is shown with the cyan and brown dashed lines. The prediction from the JED-SAD model with the Mach number, defined as the accretion speed to sound speed ratio, being unity is shown with the magenta dashed line (Marcel et al. 2018a, b). The Keplerian frequency f_K(r) is shown with the grey dotted line for a reference.

Open in new tabDownload slide

We also look at global MRI simulations of hot flows, where the angular momentum transport is self-consistently calculated by the magneto-rotational instability (MRI), so α is a function of radius. The Standard And Normal Evolution’ (SANE) model of Narayan et al. (2012) gives viscous time-scales (cyan dashed line in Fig. 10) which match well to our derived viscous time-scales. Other SANE models (Bollimpalli et al. 2020) also show similar viscous time-scales to those derived from Narayan et al. (2012) and thus match our spectral-timing model.

However, recent simulations have focused instead on Magnetically Arrested Discs (MAD; Narayan, Igumenshchev & Abramowicz 2003; Narayan et al. 2012). These have maximal magnetic flux pinned on to the black hole horizon, such that the magnetic field is strong enough that its pressure exceeds the ram pressure of the infalling material. In 1D calculations, this halts the accretion, but full 3D calculations show that accretion takes place instead via the magnetic interchange instability (White, Stone & Quataert 2019; Xie & Zdziarski 2019). The motivation for such states is that they power strong jets. MAD models have faster radial velocity than the SANE models as the large-scale magnetic fields give more efficient angular momentum transport than the small-scale dynamo (Narayan et al. 2012). However, this means that they predict an accretion speed which is too fast to match our derived constraints (brown dashed lines in Fig. 10).

The jet can also be incorporated in a rather different way, by assuming a large scale height magnetic field exists rather than calculating it ab initio (the Jet Emitting Disc models (JED); Marcel et al. 2018a). If this field is responsible for transporting angular momentum in the hot inner flow then the effective α ∼ 3. However, this can still give a good match to our velocity profile (magenta dashed line in Fig. 10) as these models have a rather small-scale height such that H/R ∼ 0.1 contrasting with the larger scale height global ADAF, SANE and MAD flows.

However, while the SANE and JED models both give a good match to our derived radial velocity, there is still a potential observable difference between them in their predicted sound speed, c_s ∼ (H/R)v_ϕ(R) (Kato et al. 2008). The SANE models have α < H/R, and thus c_s > v_r, so the sound waves called the bending waves can cross the flow before they are damped out by viscosity. The flow can then globally precess if it is misaligned with the black hole spin (the wave-like regime; e.g. Ingram & Motta 2019). This is the most compelling current model for the origin of the low-frequency QPO (Ingram et al. 2009). By contrast, the JED/MAD models have α > H/R, and thus c_s < v_r, so viscosity damps out any misalignment torques into a stable warp, which does not precess (the diffusive regime), ruling out Lense–Thirring precession as the origin of the QPO in the JED models (Marcel & Neilsen 2021). If global precession really switches off sharply at α = H/R then a Lense–Thirring origin of the QPO favours the SANE models. If instead there is a more gradual transition at α = H/R (so that precession is damped but still possible) then the JED models could give an origin for both the QPO and the jet.

8.2 Limitations of the variability model

The timing model predicts time lags significantly lower than that observed for 0.3–|$3\, {\rm Hz}$|⁠, and overestimates the frequency at which the soft (negative) time lag reaches its maximum (Fig. 8, lower panel). It is difficult to solve these discrepancies within the current spectral decomposition (fixed to that in Fig. 1). There could be a better match to the timing properties from a different set of spectral components. We fixed the emissivity of the whole variable flow, which remains unclear, to that of the standard disc throughout the paper. Different emissivity yields different timing properties, although it could be uncomfortably fine-tuned (Ingram & Done 2011, 2012). However, we note that the different emissivity does not significantly change the truncation radius in our model since it is derived from the low-frequency break in the power spectra, which simply marks the viscous frequency at the truncation radius (Ingram & Done 2011). The extensive exploration of the spectral-timing analysis is beyond the scope of the paper.

Our final model also does not reproduce the observed shift in peak frequency of P₂ in power spectra, as shown in Fig. 8 (upper panel). This is a qualitative discrepancy between the model and observation, and not solved as long as we simply attribute P₁ and P₂ to the variable disc and hot flow, respectively. In this regard, we conclude that our current model lacks some fundamental factor(s) that impacts the fast variability. The discrepancies are seen particularly at low energies where the disc emission significantly contributes to the spectrum, and mid frequencies between P₁ and P₂. This suggests a complicated interplay between the truncated disc and hot flow at their transition region. The complete discontinuity of viscous time-scale at the transition radius r_ds may be an oversimplification. Due to the variability of the seed photons supplied by the inner edge of the truncated disc, the spectral shape of Comptonisation components is likely to vary in time (Kotov et al. 2001; Veledina 2016, 2018; Mastroserio, Ingram & van der Klis 2018), which is not taken into account in the current model.

We show here a proof of concept that the variability data allow us to probe the dynamics of the X-ray emitting region in such a way as to enable us to distinguish the nature and geometry of the hot flow, and how it connects to the jet. We plan to improve the timing model and perform a more comprehensive spectral-timing analysis in future works.

9 CONCLUSIONS

We have performed a full spectral-timing analysis for MAXI J1820+070 in the low/hard state. This is based on the spectral components identified in the very detailed NuSTAR study of Zdziarski et al. (2021b), though we show the data favour an additional narrow iron line which is likely contributed from the strong wind seen in optical spectra in the state. We extended the model down to the NICER energy bandpass, including thermal emission from the disc. We see evidence for atomic features in this component also, as expected from the low-temperature disc photosphere at ∼0.22 keV.

Analysing the fast variability, we find the standard double-humped power spectrum seen in this state, indicating two variability components related to different time-scales, P₁ and P₂, which behave differently with respect to energy: whereas the component P₁ maintains both its shape and amplitude for different energies, the peak frequency shifts and the high-frequency variability increases with energy in the component P₂. The NICER data allow us for the first time to investigate this behaviour with good statistics below |$2\, {\rm keV}$|⁠.

To explore the observed timing properties, we have developed a timing model based on the propagating fluctuations of the mass accretion rate. The spectral decomposition derived from the fit is used to define the energy spectra emitted at each radius in the flow. We have developed a full spectral-timing model which incorporates discontinuity of the viscous frequency at the transition radius between the truncated disc and hot flow, as well as a spectrally inhomogeneous hot flow. The discontinuity of the viscous frequency naturally reproduces the double-humped shape corresponding to P₁ and P₂ in the power spectra, as shown by Rapisarda et al. (2016), but our full spectral-timing model also allows us to explore the observed energy dependence of P₂. We find that this can be explained by propagating fluctuations through the radially stratified (soft and hard Comptonisation components) hot flow.

We also implement reflection/reprocessing of the Compton components on the disc. Using the disc truncation at a few tens of gravitational radii, we show that this can give a good overall match to the lag-frequency spectrum, with the switch between soft leads (propagation) at low frequencies and soft lags (reverberation) at high frequencies, as well as fitting to the shape of the energy-dependent power spectra. Our results on the reverberation size scale are consistent with our assumed truncated disc/hot inner flow geometry. While our timing model does not fit all aspects of the data, we show that reverberation from a much smaller inner disc radius (as derived from some iron line fits: Buisson et al. 2019) makes these discrepancies worse.

The viscous time-scale constrains the combination of effective viscosity and scale height of the flow. The SANE models (Narayan et al. 2012), where the MRI dynamo transports the angular momentum, have fairly low effective viscosity and large-scale height, and these match well to our constraints. By contrast, the MAD models (Narayan et al. 2012) have similar scale heights, but higher effective viscosity from the large-scale magnetic fields, so predict viscous time-scales which are substantially faster than those observed. None the less, the JED models (Marcel et al. 2018a) with large-scale magnetic fields linked to the jet can match our derived viscous time-scales, where a smaller scale height flow offset the higher efficiency of the large-scale magnetic torques. However, these have a sound speed which is slower than the viscous speed. This is important as Lense–Thirring precession, currently the best model for the low-frequency QPO (Ingram et al. 2009), depends on misalignment torques transmitted by bending waves. These travel at ≈ the sound speed, but are damped out if the viscous speed is faster than the sound speed (Ingram & Motta 2019). If this really is able to completely suppress the global precession (Marcel & Neilsen 2021) then we favour models where the angular momentum transport in the hot flow is not caused by a large scale, ordered magnetic field linked to the jet, but instead is powered by the small-scale MRI dynamo. If instead there is more complex behaviour around the transition from the wave-like to diffusive regime, then models which also tie the jet into the hot flow are clearly attractive.

The timing model still has room for improvement, particularly for the low energy bands, where the disc emission dominates the energy spectrum. None the less, this study demonstrates the potential for spectral-timing models to address the fundamental nature and geometry of the accretion flow close to the black hole, and its link to the jet.

ACKNOWLEDGEMENTS

We thank G. Marcel and P. O. Petrucci for illuminating discussions about the JED models. We acknowledge valuable comments from the referee that helped to improve the manuscript. This work was supported by Japan Society for the Promotion of Science (JSPS) KAKENHI Grant Numbers 18H05463 and 20H00153, and World Premier International Research Center Initiative (WPI), Ministry of Education, Culture, Sports, Science and Technology (MEXT). TK acknowledges support from JSPS Overseas Challenge Program for Young Researchers, and Japan Science and Technology Agency (JST), Support for Pioneering Research Initiated by the Next Generation (SPRING), Grant Number JPMJST2108. CD acknowledges support from Science and Technology Facilities Council (STFC) through grant ST/P000541/1, and Kavli Institute for the Physics and Mathematics of the Universe (IPMU) funding from the National Science Foundation (No. NSF PHY17-48958).

DATA AVAILABILITY

The data underlying this article are available in HEASARC, at https://heasarc.gsfc.nasa.gov/docs/archive.html.

Footnotes

REFERENCES

APPENDIX A: MODEL CALCULATION OF POWER SPECTRA AND TIME LAGS

We summarize the practical method to calculate the power spectra and frequency-dependent time lags of the flux in our model. The method is mainly based on Ingram & van der Klis (2013). First, we take into account the inwards propagation of mass accretion rate fluctuations alone. The case, where the reflection by the disc is included, is described in the following appendices. Since the observation data are discrete time series, we should take discrete times t_j and frequencies f_m with integer subscripts j, m, rather than continuous t, f, in the formulation. However, in some contexts, we use the notation as the continuous time and frequency in order to simplify expressions.

APPENDIX B: MODEL INCLUDING REFLECTION AND REPROCESSING

First of all, taking the reprocessed emission into account changes the transition radii r_ds, r_sh because the total flux of the intrinsic disc emission decreases to (1 − f_rep)S_d(E) and so in the relation used to derive the transition radii (3), S_d(E) is replaced by (1 − f_rep)S_d(E). We summarize alterations in the variability of the flux below.

The former component x_dir(E, t) is written as the form of equation (7) with the slight modification of the weight (8), that is S_d(E) is replaced by (1 − f_rep)S_d(E). We define the modified weight as w_dir(r_n, E).

The formulation described above is extended to the case, where a wide energy band E_min < E ≤ E_max is taken. In the case, the weight for the direct component w_dir(E, t) is modified into w_dir((E_min, E_max), t), which is obtained from equation (9) with the substitution of (1 − f_rep)S_d(E) for S_d(E). The calculation of the normalization of the impulse response is described in Appendix  C. It is noted that S(E) in equation (8) should be replaced by |$S(E) A_{\mathrm{eff}} (E) \mathrm{e}^{-N_{\mathrm{H}}(E)\sigma _{\mathrm{T}}}$| to take the different statistics for different energy bins in observations into account.

APPENDIX C: AMPLITUDE OF THE IMPULSE RESPONSE OF THE REPROCESSING

APPENDIX D: COMPARISON OF TIMING PROPERTIES WITH DIFFERENT PARAMETER VALUES

We show how model parameters affect the timing properties, picking the important parameters: the delay time of the reflection t₀, the fraction of the reprocessed emission f_rep, and the power index of the viscous frequency in the hot flow m_f. As in Section 7.3, the model including both the discontinuity of the viscous frequency and the reflection is used. The power spectra at the low energy band of 0.5–|$1.5\, {\rm keV}$| and lag-frequency spectra between 0.5–|$1.5\, {\rm keV}$| and 2.0–|$10.0\, {\rm keV}$| calculated with various parameter values are shown in Fig. D1. We explain the dependence of timing properties on each model parameter in order below.

Figure D1.

Power spectra at the low energy band of 0.5–|$1.5\, {\rm keV}$| (upper) and lag-frequency spectra between 0.5–|$1.5\, {\rm keV}$| and 2.0–|$10.0\, {\rm keV}$| (lower) calculated from various values of model parameters, t₀ (left), f_rep (middle), and (B_f, m_f) (right). The black lines are calculated with the parameter values tabulated in the bottom row of Table 3. Other lines are calculated by varying one or two parameter(s) with respect to the black lines. Note that, although the time delay of the reflection t₀ and the fraction of the reprocessed emission f_rep should be related to the radii setting spectral regions in the accretion flow, r_sh, r_ds, and r_out (e.g. Appendix  B), these radii are fixed for simplicity. In the right column, the normalization parameter of the viscous frequency in the hot flow, B_f, is calculated such that the viscous frequency at the outer edge of the hot flow r_ds corresponds to the same value for different power indices m_f. For comparison, the observation results are shown with the stepped lines.

Open in new tabDownload slide

From the left-top panel, the power spectrum is insensitive to the delay time t₀ because the power spectrum does not contain the information on the phases of the Fourier components, which is affected by t₀. On the other hand, the lag-frequency spectrum is affected by t₀ at high frequencies, where the reflection contributes to the variability, as shown in the left-bottom panel. The longer time delay increases the amplitude of the soft lag and makes the frequency at which the switch from the hard to soft lag occurs lower.

From the middle column, both the power spectrum and lag-frequency spectrum are generally insensitive to the fraction of the reprocessed emission f_rep if we use the same radii setting the spectral regions in the accretion flow. The increase in f_rep slightly increases the variability of P₂ in the low energy band, which is generated as a result of the reflection. The increase in f_rep also slightly contributes to the shift of the time lags to the negative direction by decreasing the intrinsically variable disc emission resulting from the propagating mass accretion rate fluctuations.

The right column shows the dependence of timing properties on the power-index parameter of the viscous frequency in the hot flow m_f. The normalization B_f is calculated such that the viscous frequency at the outer edge of the hot flow, |$\lim _{r \nearrow r_{\mathrm{ds}}}f_{\mathrm{visc}}(r)=B_{\mathrm{f}}r_{\mathrm{ds}}^{-m_{\mathrm{f}}} f_{\mathrm{K}}(r_{\mathrm{ds}})$| (see equation 13), corresponds to the same value for different m_f. The bigger power-law index m_f means the larger variability at high frequencies, as shown in the right-top panel. The lag-frequency spectra in the right-bottom panel show complicated dependence on m_f. For the larger m_f, the transition from the hard to soft lag happens on a wider range of frequencies, although the frequency at which the time lag takes zero seems irrelevant to m_f.

Fig. D1 also points to the difficulty in reproducing the timing properties related to low energy bands, where the disc contributes to the energy spectrum. We plan to modify the timing model to explain the timing properties at these energy bands.

APPENDIX E: CALCULATION OF REALISTIC IMPULSE RESPONSES

We describe the calculation of impulse response of the outer disc for the X-ray illumination from the primary source with the lamp-post and extended geometries. For simplicity, we assume that the outer disc is flat with infinitesimal thickness and optically thick and that the extended source is optically thin and radiation from everywhere in the source reaches the outer disc. We formulate the case of the extended illuminating source first. The impulse response for the lamp-post geometry can be calculated in the same way and more easily, as stated later.

For the lamp-post source geometry, the calculation of impulse response becomes simpler. The integration for the source vanishes from equation (E1), and the consideration of emissivity of the source in equation (E2) is no longer necessary. Regardless of formulation, the amplitude of the impulse response h(t) is calculated from the time-averaged spectrum of the reflected/reprocessed emission (see Appendix  C for the case of top-hat impulse response).