Large and complex X-ray time lags from black hole accretion disks with compact inner coronae

Phil Uttley, Julien Malzac

► To cite this version:

Phil Uttley, Julien Malzac. Large and complex X-ray time lags from black hole accretion disks with compact inner coronae. Monthly Notices of the Royal Astronomical Society, 2025, 536 (4), pp.3284-3307. 10.1093/mnras/stae2514. hal-04356179

HAL Id: hal-04356179

https://hal.science/hal-04356179v1

Submitted on 18 Jul 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Large and complex X-ray time lags from black hole accretion discs with compact inner coronae

Phil Uttley and Julien Malzac<sup>2,3</sup>

Accepted 2024 October 14. Received 2024 October 13; in original form 2023 December 14

## ABSTRACT

Black hole X-ray binaries in their hard and hard-intermediate states display hard and soft time lags between broad-band noise variations (high-energy emission lagging low-energy and vice versa), which could be used to constrain the geometry of the disc and Comptonizing corona in these systems. Comptonization and reverberation lag models, which are based on light-travel delays, can imply coronae that are very large (hundreds to thousands of gravitational radii, $R_g$) and in conflict with constraints from X-ray spectral modelling and polarimetry. Here, we show that the observed large and complex X-ray time lags can be explained by a model where fluctuations are generated in and propagate through the blackbody-emitting disc to a relatively compact (~10 $R_g$) inner corona. The model naturally explains why the disc variations lead coronal variations with a Fourier-frequency dependent lag at frequencies < 1 Hz, since longer variability time-scales originate from larger disc radii. The propagating fluctuations also modulate successively the coronal seed photons from the disc, heating of the corona via viscous dissipation and the resulting reverberation signal. The interplay of these different effects leads to the observed complex pattern of lag behaviour between disc and power-law emission and different power-law energy bands, the energy-dependence of power-spectral shape, and a strong dependence of spectral-timing properties on coronal geometry. The observed spectral-timing complexity is thus a natural consequence of the response of the disc-corona system to mass-accretion fluctuations propagating through the disc.

**Key words:** accretion, accretion discs – X-rays: binaries – X-rays: individual: MAXI J1820 + 070.

## 1 INTRODUCTION

The X-ray spectra of hard state black hole X-ray binaries (BHXRBs) show two primary continuum components: low-temperature ($k_BT$ < 0.4 keV) blackbody emission from the accretion disc and power-law-like emission produced by inverse Compton scattering of ‘seed’ photons in hot ($k_BT_e$ ~ 100 keV) thermal plasma of optical depths $\tau \sim 1$ (e.g. Zdziarski et al. 1998; Makishima et al. 2008). For moderate- to high-luminosity hard states, the seed photons should primarily originate from the disc (Done, Gierliński & Kubota 2007). Observed X-ray reflection features further show that the disc reprocesses at least tens of per cent of the power-law luminosity, implying that the power-law emitting region has some vertical extent compared to the thin disc, although the exact geometry is still being debated (e.g. Buisson et al. 2019, Zdziarski et al. 2021). Thus, we can refer to this hot, Compton-scattering region as a ‘corona’ without making any further assumptions about its geometry.

Unlike the soft state which shows stable disc emission (fractional rms ~ 1 per cent or smaller, e.g. Belloni et al. 2005), both disc and coronal emission in the hard state are highly variable, showing correlated aperiodic noise variability over a broad range of timescales (Wilkinson & Uttley 2009). Significant constraints on the

relationship of the disc and corona as well as their size scale and relative geometry can be provided by X-ray spectral-timing, such as the measurement of Fourier time lags and power spectra for different energy bands. Fig. 1 shows example lags and power spectra for disc (0.5–0.9 keV) and power-law-dominated (2–3 keV, 5–10 keV) energy bands, for hard state data obtained for the black hole X-ray binary MAXI J1820 + 070¹ with the Neutron star Interior Composition Explorer (NICER; Gendreau et al. 2016).

The hard lags between power-law continuum photons (henceforth, PL-PL hard lags) with variations at higher photon energies lagging those at lower energies, were the first to be discovered and studied in some detail (in Cyg X-1, Miyamoto et al. 1988; Cui et al. 1997; Nowak et al. 1999). The lags are typically < 1 per cent of the variability time-scale but depend on Fourier frequency $v$, with time lags decreasing with increasing Fourier frequency roughly as $\tau \propto v^{-0.7}$, albeit with a ‘stepped’ appearance (e.g. see Nowak 2000; Misra et al. 2017). The lags also depend on energy $E$, with the time lag between energies $E_2$ and $E_1$, $\tau \propto \log(E_2/E_1)$ (Nowak et al. 1999; Kotov, Churazov & Gilfanov 2001). The evolution of PL-PL hard lags through the hard state for multiple BHXRBs shows a

© 2024 The Author(s).

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

<footnote n="1"><sup>1</sup>Anton Pannekoek Institute, University of Amsterdam, Science Park 904, NL-1098 XH Amsterdam, the Netherlands</footnote>

<footnote n="2"><sup>2</sup>Université de Toulouse, UPS-OMP, IRAP, Toulouse, France</footnote>

<footnote n="3"><sup>3</sup>CNRS, IRAP, 9 Av. colonel Roche, BP44346, F-31028 Toulouse cedex 4, France</footnote>

<footnote n="1">¹ObsID 1200120104, with 6 ks exposure obtained 2018 March 15 and processed using the approach outlined in Wang et al. (2022). See Uttley et al. (2014) for details on the analysis methods.</footnote>

<footnote>*E-mail: p.uttley@uva.nl</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** MAXI J1820 + 070 hard state frequency-dependent lags (top panel) and power-spectra (lower panel) measured from data obtained with NICER on 2018 March 15, plotted as a function of Fourier frequency for the 0.5–0.9, 2–3, and 5–10 keV bands. The X-ray spectrum for this data set is shown and discussed in Appendix D.</figcaption></figure>

general trend of average lag increasing with the coronal power-law continuum photon index Γ (Pottschmidt et al. 2003; Grinberg et al. 2014; Altamirano & Méndez 2015; Reig et al. 2018).

The power-law versus disc lags (henceforth PL-D lags) were first discovered in XMM–Newton data for GX 339-4 (Uttley et al. 2011) and show apparently ubiquitous behaviour across many sources, which is notably more complex than seen for the PL-PL lags (e.g. De Marco et al. 2015b, 2017; Wang et al. 2022). At low frequencies the lags are positive and substantially larger than the PL-PL lags seen for the same ratio of energies at the same frequencies. At high-frequencies however the lags turn over to become negative, i.e. soft lags, of order 1–10 ms, with the soft lag amplitude tracking the timing and spectral-state evolution closely from the hard state through to the hard-intermediate state (HIMS) to soft-intermediate state (SIMS) transition (Kara et al. 2019; De Marco et al. 2021; Wang et al. 2021, 2022).

Hard state power spectra also show a significant energy dependence. Besides the larger low-frequency variability amplitude at lower energies<fnref n="2" />, seen in Fig. 1, the power-spectrum shows a steeper high-frequency slope at lower energies (Nowak et al. 1999; Grinberg et al. 2014). This effect is particularly clear for the disc band (see also De Marco et al. 2015a).

To date there is no model which can explain the complex pattern of spectral-timing behaviour of both disc and power-law emission. For example, PL-PL hard lags and power-spectral energy-dependence have been attributed to the differential light-travel delays expected from inverse-Compton scattering of photons to different energies in a spherical corona (Kazanas, Hua & Titarchuk 1997; Hua, Kazanas & Cui 1999) or jet (Reig, Kylafis & Giannios 2003; Giannios, Kylafis & Psaltis 2004) of radially varying density. Since they link the lags to

the spectral-formation process, jet Comptonization models are able to explain correlations between the continuum power-law photon index Γ and PL-PL hard lag amplitude in terms of the jet geometric parameters (Kylafis et al. 2008). However, such models imply very large vertically (and horizontally) extended power-law emitting regions, which predict substantially flatter emissivity profiles than are inferred from reflection model fits to data (e.g. see Fig. 4 in Reig & Kylafis 2021 for the calculated flat illumination pattern). Furthermore, pure Comptonization models do not posit a physically-motivated origin of the variability, assuming simply that seed photons vary only where they are released at the base of the corona. This should lead to very short delays between the disc seed photons and lower-energy Comptonized photons, contrary to the large observed PL-D hard lags.

A solution to these difficulties can be found if we assume that variations are driven by mass-accretion fluctuations which propagate through the accretion flow, i.e. the PL-D lags are linked to the travel time of fluctuations between the disc and a central corona. Propagating accretion fluctuations are a natural feature of turbulent accretion flows (Lyubarskii 1997) and can explain a broad variety of aspects of XRB variability, such as their broadband power-spectral shapes (Churazov, Gilfanov & Revnivtsev 2001; Ingram & Done 2011), linear rms-flux relations and lognormal flux distributions (Uttley, McHardy & Vaughan 2005; Heil, Vaughan & Uttley 2012), the similarity between neutron star and black hole timing behaviour and its evolution (Sunyaev & Revnivtsev 2000; Gardenier & Uttley 2018) and XRB spectral-timing properties (Kotov et al. 2001; Arévalo & Uttley 2006).

To date, propagating fluctuation models have focussed on the PL-PL hard lags, by assuming that the corona corresponds to a hot inner accretion flow, i.e. within the disc truncation radius, with harder emission produced more centrally. Lags are produced by the propagation delays between radii with softer and harder emission, decreasing for smaller radii, which correspond to higher variability frequencies generated in the flow. Successively more complex versions of the model have been proposed to explain the energy and frequency dependence of the lag and the power spectrum by invoking radial discontinuities in: the emission spectrum (Rapisarda et al. 2016; Mahmoud & Done 2018a); input variability signal (Mahmoud & Done 2018b) and speed of propagation (Kawamura et al. 2022). However, these models assume ad hoc prescriptions for the radial dependence of the power-law continuum that best model the data, so they do not have explanatory power for, e.g. the Γ-lag correlation. They also do not explain the PL-D lags well (or they do not explain them at all), nor do they consider the role of disc seed photons in Comptonization.

Recent interest in the use of X-ray spectral-timing to constrain the disc-corona geometry has focussed on X-ray reverberation, where disc reprocessing of the illuminating coronal emission produces disc emission which lags the coronal variations by the light-travel delay from corona to disc. Reverberation thus provides a natural explanation of the transition from hard to soft PL-D lags seen at high frequencies (Uttley et al. 2011). The role of reverberation in producing the short-term soft lags was later confirmed with the NICER detection of delays in broadened iron Kα reflection (Kara et al. 2019). Current state-of-the-art reverberation models calculate light travel delays using relativistic ray-tracing and apply them to sophisticated reflection models to fit the lags over a broad energy range (Mastroserio, Ingram & van der Klis 2018; Mastroserio et al. 2021). For simplicity, a lamppost source geometry is commonly used, but more extended coronal geometries are being investigated (Lucchini et al. 2023). However, the large soft lags push current

<footnote n="2">²The low-frequency energy-dependence corresponds to the short-time-scale ‘softer-when-brighter’ behaviour of moderate–high luminosity hard-state BHXRBs (Skipper, McHardy & Maccarone 2013; Skipper & McHardy 2016).</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Summary of the propagating fluctuation lag model presented in this work. A fluctuation signal will produce a combination of the effects A, B, and C (shown *left*) depending on the radius of origin of the signal. Since the radius of origin determines the signal frequency, frequency-dependent lags are seen (*right*, plotted schematically with sinh⁻¹ scaling to show positive and negative lags on the same axis). For signals originating at large radii (*A*) the large propagation delays dominate the lags, leading to positive power-law versus disc lags as *A* is followed by *B* and *C*. However, for signals produced within the seed emitting region, the disc emission is dominated by the reverberation component (dot-dashed arrows) which is produced mostly in *C* and lags the soft power-law emission produced mostly at *B* by the propagation delay between *B* and *C*. Frequency-dependent hard lags between power-law bands are also produced due to the delay between *B* and *C*. For the highest frequencies, which are produced in the corona, propagation-related delays are negligible and light-travel delays may start to dominate (not shown in the calculated lags).</figcaption></figure>

reverberation models to infer large coronal heights, ranging from tens of $R_g$ in the hard state to hundreds of $R_g$ towards the end of the HIMS, prior to the state transition (Wang et al. 2021). While not as large as predicted by pure Comptonization models of PL-PL hard lags, these coronal heights remain difficult to reconcile with reflection spectral fits. Furthermore, the models do not include intrinsic disc variability and rely on fitted empirical prescriptions for power-law pivoting to account for the observed hard lags. Reconciling reverberation delays with propagation delays to model other spectral-timing behaviour remains challenging (see Mahmoud, Done & De Marco 2019 for a recent example).

In summary then, individual lag models fall well short of a complete explanation of the complex spectral-timing behaviour observed in black hole X-ray binaries in the hard state. To some extent, Comptonization, propagation, and reverberation must all play a role in solving the problem, but current implementations of these effects are limited because they neglect the other aspects, e.g. either the physical origin of the continuum spectrum (in the case of propagation or reverberation models) or the origin of variability itself (in the case of Comptonization models).

In this paper, we show how to resolve these difficulties by implementing a model which includes the effects of propagation and reverberation together with a self-consistent approximation for power-law continuum formation in the corona. We do so by making the simplifying assumption that, on the observed variability time-scales, the Comptonizing corona is in thermal equilibrium with respect to seed and coronal heating variations, with normalization and $\Gamma$ set by the seed photon flux and the ratio of seed photon to coronal heating luminosities (as expected from conservation of both photon number and total power in the corona).

Our model is shown schematically in Fig. 2. Following previous propagation models, mass accretion fluctuations are generated with radially dependent time-scales in both the disc and corona, which we treat as an inner hot flow but with a vertical extent described by the coronal geometry. Fluctuations propagating through the disc naturally produce the observed frequency-dependent PL-D hard lags

and suppression of soft band variability at high frequencies. Seed photons are produced in the disc and are modulated by accretion fluctuations before those same fluctuations reach the corona and modulate the coronal heating rate. This means that the coronal power-law spectrum pivots in response to any given fluctuation, with steepening in response to slower fluctuations from further out in the disc and hardening in response to variations propagating through the corona. This effect naturally produces hard PL-PL lags on long time-scales as well as flatter high-frequency power spectra at harder power-law continuum energies.

Although our model was originally conceived to simultaneously explain only the PL-D and PL-PL hard lags, the inclusion of reverberation of coronal luminosity on the disc naturally produces large soft PL-D lags at high Fourier frequencies, even for compact (~10 $R_g$) coronae and under the assumption of negligible light-travel delays. These lags arise from the delay between disc reverberation, which is driven by the total coronal power and the soft power-law emission, which is driven by seed photon variations and thus precedes the peak in coronal power due to the propagation delay of fluctuations between the inner disc and corona. The model can naturally explain changes in PL-PL hard lags and PL-D soft lags, in terms of changes in the visibility of disc seed photons from the corona that are linked to changes in coronal geometry.

Here in this work, we introduce the model in its most basic form, considering the lags between disc total luminosities and mono-energetic power-law bands, and treating the seed photons as being mono-energetic. We consider the dependence of lags on disc photon energy briefly in the Appendix, but will carry out a full energy-dependent treatment and fits to data in later papers. In Section 2, we use an impulse response function approach to show that hard power-law lags with a log-linear energy dependence are automatically produced by the corona when seed photons vary before coronal heating. In Section 3, we show how this effect is automatically realized by a variable accretion flow and we describe a numerical implementation of this model in Section 4. In Section 5, we show how the model reproduces most of the key spectral-timing results
discussed so far. We discuss the model implications and its limitations in Section 6, before setting out our conclusions in Section 7. In the Appendix, we demonstrate some key results and use numerical Monte Carlo simulations of variable accretion to demonstrate the validity of our impulse response approach in the non-linear variability regime, before also considering the energy-dependence of disc lags and the consistency of our model geometries with hard state spectral-energy distributions.

## 2 TIME RESPONSE OF CORONAL EMISSION

We consider Comptonizing coronae in thermal equilibrium, where electrons are heated via some unspecified heating mechanism with heating rate given by a luminosity $L_h$. This assumption is justified for coronae with light-crossing times < 100 $R_g/c$ and the expected moderate optical depths ($\tau \sim 1$), and is observationally supported by the high coherence of variations in different energy bands in the hard state, which imply that we see temporally connected emitting regions, with any delays described by simple linear transforms (Nowak et al. 1999). The electrons are cooled via inverse Compton scattering of seed photons with luminosity $L_s$. We then assume that the instantaneous photon index $\Gamma$ is set by the ratio of the time-dependent input seed luminosity to the coronal heating luminosity as follows:

$$ \Gamma(t) = \Gamma_0 \left( \frac{L_s(t)}{L_h(t)} \right)^\beta, \quad (1) $$

where $\Gamma_0$ is the photon index for equal seed and heating luminosities and $\beta$ is an index which may be derived from the specific Comptonization model (e.g. see Pietrini & Krolik 1995 and Beloborodov 2001).

Since the model is linearized, the precise form of the relation between photon index and cooling/heating ratio does not matter as much as the direction of the relation (which sets the sign of the lags). A positive correlation between $\Gamma$ and the seed-to-heating luminosity ratio, arises naturally from conservation of both total luminosity and photon number in the corona. The assumption of thermal equilibrium is valid when the heating and cooling time-scales within the corona are small compared to the time-scales of variability of seed and heating luminosities. This should be the case for the relatively compact coronae considered here.

To calculate the variation of a quantity $f(t)$ in response to a perturbing signal $s(t)$, we require the impulse response function (henceforth 'impulse response') for that quantity $g(\tau)$ (where $\tau$ is the time delay), so that the perturbation in $f(t)$, $\delta f(t) = \int_{-\infty}^{\infty} s(t-\tau) * g(\tau)d\tau$. We will now derive the impulse responses for power-law continuum slope and flux variations, in terms of the coronal seed and heating variations.

For our simplified analytical calculation, we make the approximation that the seed photons are monochromatic with energy $E_s$, so that the total photon number is $L_s/E_s$. Assuming that total photon numbers are conserved during the Comptonization process and that the cut-off photon energy $E_{cut} \gg E_s$, the photon flux density variation at an energy $E$, $N(E, t)$ is given by

$$ N(E, t) = \frac{L_s(t)(\Gamma(t) - 1)}{E_s^2} \left(\frac{E}{E_s}\right)^{-\Gamma(t)}. \quad (2) $$

By linearizing equations (1) and (2), i.e. assuming small perturbations in the time-variable parameters we obtain the impulse responses for $\Gamma(t)$ and $N(E, t)$:

$$ g_{\Gamma}(\tau) = \langle \Gamma \rangle \beta \left( \frac{g_s(\tau)}{\langle L_s \rangle} - \frac{g_h(\tau)}{\langle L_h \rangle} \right), \quad (3) $$

$$ g_{\text{pl}}(E, \tau) = \langle N(E) \rangle \left( [1 - u(E)] \frac{g_s(\tau)}{\langle L_s \rangle} + u(E) \frac{g_h(\tau)}{\langle L_h \rangle} \right), \quad (4) $$

where $u(E) = \beta \langle \Gamma \rangle (\ln(E/E_s) - \frac{1}{\langle \Gamma \rangle - 1})$ and angle brackets denote the time-averaged quantities which are perturbed by changes in the seed and heating luminosities. Those changes are themselves quantified in terms of their own impulse responses, $g_s$ and $g_h$, which will depend on the physical system producing the seed and heating luminosities and their variations. We will determine these functions for the case of propagating accretion fluctuations in the next section.

The $u(E)$ term contains the energy-dependence of the impulse response shape, as it weights the heating and seed impulse responses differently according to the energy of power-law photons being considered. Effectively, $u(E)$ produces pivoting of the power-law spectrum about an energy $E_{piv}$, which corresponds to $g_{\text{pl}}(E_{\text{piv}}, \tau) = 0$ and is given by

$$ \ln \left( \frac{E_{\text{piv}}}{E_s} \right) = \left[ \beta \langle \Gamma \rangle \left( 1 - \frac{g_h(\tau)}{\langle L_h \rangle} / \frac{g_s(\tau)}{\langle L_s \rangle} \right) \right]^{-1} + \frac{1}{\langle \Gamma \rangle - 1}, \quad (5) $$

For $\beta = 1/6$ (Beloborodov 2001) and typical observed photon indices $1.4 < \langle \Gamma \rangle < 2.5$, if a change in seed luminosity dominates, the minimum pivot energy is $E_{\text{piv}} \approx 21-885 E_s$ (with the smallest pivot energy for the steepest spectra), which corresponds to $g_h(\tau) = 0$. As $g_h(\tau)$ increases relative to $g_s(\tau)$, the pivot energy increases rapidly until $g_h(\tau)/\langle L_h \rangle = g_s(\tau)/\langle L_s \rangle$, when $E_{\text{piv}} = \infty$ and only the power-law normalization changes, while photon index does not vary. If $g_h(\tau)/g_s(\tau)$ continues to increase, the pivot switches to very low energies, increasing towards a limiting maximum energy $E_{\text{piv}} \approx 2-12 E_s$ (with steeper spectra corresponding to lower maximum pivot energies), which corresponds to $g_s = 0$. Thus, differences between the time responses of the seed and heating luminosities lead to changes in the spectral pivot-point and significant time-delay dependent changes in where the largest and smallest power-law flux variations are seen. This mechanism naturally leads to energy-dependent lags in the power-law emission, as well as energy dependence of the power-law flux variability amplitude.

## 3 IMPULSE RESPONSES FROM A VARIABLE ACCRETION FLOW

We assume that the driving signals are mass accretion rate fluctuations propagating through the disc. For simplicity, we will assume simple non-diffusive propagation, so that the amplitudes and power-spectral shapes of mass-accretion fluctuations are preserved as they propagate inwards. This approximation can be applied to 'classical' accretion discs provided that the fluctuation time-scale is comparable to or greater than the viscous time-scale at the radius it is generated (Churazov et al. 2001; Mushtukov, Ingram & van der Klis 2018).

The response delay $\tau^3$ depends on the radial drift velocity of the accretion variations. Since the final model will combine signals from multiple radii, it is useful to set $\tau = 0$ to be the time when a propagating signal reaches the innermost radius of the accretion flow, $r_{\text{in}}$, which allows us to define a fixed, negative, propagation time-delay for each radius $r$ (where $r$ is expressed in units of the gravitational radius $R_g$):

$$ \tau(r) = \int_{r_{\text{in}}}^{r} \mathrm{d}\tau = \int_{r_{\text{in}}}^{r} \mathrm{d}r / v_r, \quad (6) $$

<footnote n="3">³Throughout this paper we use $\tau$ to denote a relative time delay, while $t$ is used to represent a time or time-scale.</footnote>
where the radial drift velocity $v_r = -\alpha(h/r)^2 r^{-1/2}$, for a disc with scale-height $h$ and dimensionless viscosity parameter $\alpha$. Note that $v_r$ is expressed as a fraction of the speed of light $c$, so the delays are naturally expressed in units of the light-crossing time for 1 $R_g$.

For simplicity, we assume that light-travel time lags, e.g from the disc to the corona, are negligible compared to the lags due to propagating fluctuations which we consider here. We have also assumed that $r_{in}$ is the inner zero-torque boundary of the whole accretion flow (including the coronal part), which we assume to correspond to the innermost stable circular orbit (ISCO). It's also important to note that the radial dependence of propagation delays (and time-scales of variability, discussed later) will be significantly affected by the radial scaling of the disc aspect ratio $h/r$.

We define all luminosities as fractions of the total dissipated power from the flow. The power dissipated by the accretion flow in an annulus $dr$ at radius $r$ is therefore defined, using the usual radial dependence (e.g. Frank, King & Raine 2002) as

$$f_{\text{diss}}(r)dr = \frac{r^{-2}[1 - (r_{\text{in}}/r)^{1/2}]dr}{\int_{r_{\text{in}}}^{r_{\text{out}}} r^{-2}[1 - (r_{\text{in}}/r)^{1/2}]dr}, \quad (7)$$

where $r_{out}$ is the outer radius of the flow. At this stage, we speak in terms of the accretion flow rather than the accretion disc, because we assume that the central compact corona is itself powered by accretion and may be thought of as a hot inner flow with radius $r_{cor}$, which corresponds to the truncation radius of the blackbody-emitting disc.

## 3.1 Coronal heating and seed luminosity

If we assume that all power dissipated within $r_{cor}$ goes into heating the corona, we can define the coronal heating impulse response function $g_h(\tau)$ as

$$g_h(\tau)d\tau = \begin{cases} f_{\text{diss}}(r)dr, & \text{if } r_{\text{in}} < r \le r_{\text{cor}} \\ 0, & \text{otherwise} \end{cases} \quad (8)$$

Seed photons are provided by the disc. If the fraction of disc emission from a radius $r$ which is intercepted by the corona is $f_{d\to c}(r)$, the seed photon impulse response due to disc photons produced by viscous dissipation is

$$g_{s,\text{diss}}(\tau)d\tau = \begin{cases} f_{d\to c}(r)f_{\text{diss}}(r)dr, & \text{if } r_{\text{cor}} < r \le r_{\text{fluc}}, \\ 0, & \text{otherwise} \end{cases} \quad (9)$$

where $r_{fluc}$ is the starting radius of the accretion fluctuation (which we assume to begin in the disc, i.e. outside the corona). Note that we include the additional subscript 'diss' to make a distinction between disc (and seed) emission which is produced directly by dissipation in the flow, and that produced by reprocessing of the coronal power-law emission by the disc. The calculation of $f_{d\to c}(r)$ depends on the assumed coronal geometry and is discussed in Appendix A. For the present work, we assume that the disc is the only source of seed photons to the corona. It would be simple however, to allow some fraction of the viscous dissipation powering the corona to be converted into internal seed photons, e.g. via synchrotron emission if the corona is magnetized (Veledina, Poutanen & Vurm 2013).

A fraction of the emitted total coronal luminosity (which is the sum of seed and heating luminosities) will be intercepted and reprocessed by the disc to produce an X-ray reverberation signal which has a short, light-travel time-delay with respect to the coronal X-ray variations. Some fraction of this reverberation component will return to the corona and can provide an extra source of seed photons. The fraction, $f_{\text{return}}$, of the total coronal luminosity which returns to the corona in

this way is

$$f_{\text{return}} = 2\pi \int_{r_{\text{cor}}}^{r_{\text{out}}} f_{d\to c}(r) f_{c\to d}(r) r dr, \quad (10)$$

where $f_{c\to d}(r)$ is the fraction per unit area of coronal photons intercepted by the disc, as a function of disc radius (see Appendix A for details of its calculation).

The returning photons will be inverse Compton-scattered in the corona and may undergo additional reprocessing in the disc to contribute higher orders of returning luminosity. Assuming negligible light-travel delays to maintain coronal energy balance, we obtain for the time-dependent total returning luminosity $L_{\text{return}}(t)$:

$$\begin{aligned} L_{\text{return}}(t) &= f_{\text{return}} (L_{\text{h}}(t) + L_{\text{s,diss}}(t) + L_{\text{return}}(t)) \\ &= f_{\text{return}} \left( \frac{L_{\text{h}}(t) + L_{\text{s,diss}}(t)}{1 - f_{\text{return}}} \right) \end{aligned} \quad (11)$$

A few tens of per cent of the returning luminosity is in the form of backscattered, predominantly harder X-rays (the so-called reflection hump). The remaining luminosity, besides a few per cent in isolated fluorescent emission lines such as Fe Kα, may be thermalized by the disc, or comprise of softer reprocessed emission, which for higher disc densities common in X-ray binaries can start to approximate the smooth blackbody continuum shape of the disc thermal emission (Ross & Fabian 2007; García et al. 2016; Mastroserio et al. 2021). For simplicity, we assume that this latter 'thermal reverberation' component provides returning seed photons to the corona at similar energies to the seed photons from viscous dissipation in the disc, with luminosity equal to a fraction $f_{\text{therm}}$ of the total returning luminosity. We implicitly assume that the remaining returning photons either do not affect the coronal spectral shape or may heat or cool the corona in a way which is subsumed by $f_{\text{therm}}$, which acts as a kind of efficiency factor for the contribution of returning luminosity towards coronal cooling. Using this factor and equation (11), we infer the impulse response of the seed luminosity associated with thermal reverberation:

$$g_{s,\text{rev}}(\tau)d\tau = \frac{f_{\text{therm}} f_{\text{return}}}{(1 - f_{\text{return}})} [g_{s,\text{diss}}(\tau) + g_{h}(\tau)] d\tau. \quad (12)$$

## 3.2 Direct disc emission and thermal reverberation

The direct disc emission is the non-reprocessed disc blackbody emission which is due to viscous dissipation, observed at infinity (averaging over all observer viewing angles).

$$g_{d,\text{diss}}(\tau)d\tau = \begin{cases} (1 - f_{d\to c}(r))f_{\text{diss}}(r)dr, & \text{if } r_{\text{cor}} < r \le r_{\text{fluc}} \\ 0, & \text{otherwise} \end{cases} \quad (13)$$

To conserve luminosity, the direct disc emission is reduced by the fraction of the disc emission which is intercepted by the corona, although we note that for a complete treatment the angular dependence of this component should be accounted for. We must also account for the part of the disc reverberation signal considered in the previous section which goes directly to the observer instead of returning to the corona. This component will track the total coronal luminosity variation with negligible delay in our propagating fluctuations model. The total fraction of coronal luminosity that is intercepted and reprocessed by the disc, $f_{\text{rev}}$, is equal to

$$f_{\text{rev}} = 2\pi \int_{r_{\text{cor}}}^{r_{\text{out}}} f_{c\to d}(r) r dr. \quad (14)$$

The reverberation signal is enhanced by the contribution of the returning luminosity to the coronal emission and must also be
corrected for it (since the returning luminosity is not observed directly as reverberation). We also only consider the (quasi-)thermal part of the reverberation signal, which contributes to the total impulse response of disc thermal emission. The resulting impulse response of disc thermal reverberation is

$$g_{\mathrm{d,rev}}(\tau)d\tau = \frac{f_{\mathrm{therm}}(f_{\mathrm{rev}} - f_{\mathrm{return}})}{(1 - f_{\mathrm{return}})} [g_{\mathrm{s,diss}}(\tau) + g_{\mathrm{h}}(\tau)] d\tau. \quad (15)$$

### 3.3 Total impulse responses, luminosities, and power-law flux impulse response

The final total disc and seed photon impulse responses are

$$g_{\mathrm{d}}(\tau)d\tau = [g_{\mathrm{d,diss}}(\tau) + g_{\mathrm{d,rev}}(\tau)]d\tau, \quad (16)$$

$$g_{\mathrm{s}}(\tau)d\tau = [g_{\mathrm{s,diss}}(\tau) + g_{\mathrm{s,rev}}(\tau)]d\tau. \quad (17)$$

The time-averaged luminosities for the heating, disc, and seed components can be determined by integrating the impulse responses over all radii (implicitly setting $r_{\mathrm{fluc}} = r_{\mathrm{out}}$ for the purposes of this calculation), e.g. $\langle L_{\mathrm{h}} \rangle = \int_{\tau(\tau_{\mathrm{in}})}^{\tau(\tau_{\mathrm{out}})} g_{\mathrm{h}}(\tau)d\tau$. Note that using the impulse response definitions described in Sections 3.1 and 3.2 the total luminosity reaching the observer is correctly conserved, i.e.:

$$\langle L_{\mathrm{total}} \rangle = \langle L_{\mathrm{d,diss}} \rangle + \langle L_{\mathrm{d,rev}} \rangle / f_{\mathrm{therm}} + (1 - f_{\mathrm{rev}}) (\langle L_{\mathrm{h}} \rangle + \langle L_{\mathrm{s,diss}} \rangle + \langle L_{\mathrm{s,rev}} \rangle / f_{\mathrm{therm}}) = 1, \quad (18)$$

where the factors including $f_{\mathrm{therm}}$ and $f_{\mathrm{rev}}$ correct, respectively, for the non-thermal part of the reflection spectrum and the fraction of coronal luminosity reaching the observer.

Finally, the impulse response of power-law emission can be calculated from the impulse responses and time-averaged luminosities for the coronal heating and total seed component derived above, using equation (4). The time-averaged photon index $\langle \Gamma \rangle$ is part of the pivoting term $u(E)$ and can be determined by using the seed and heating time-averaged luminosities in equation (1), which is appropriate for the linear approximation of power-law variability used in this work. Note that the time-averaged power-law flux can be calculated from the time-averaged quantities using equation (2), although since it is only a normalizing factor for the flux impulse response it does not have any influence on our results.

### 3.4 Calculation of spectral-timing products

The impulse responses derived above can be used to determine spectral-timing properties in response to a mass-accretion fluctuation signal which originates at a single radius $r_{\mathrm{fluc}}$. In reality, it is likely that propagating fluctuations originate over a range of radii, as surmised by the model of Lyubarskii (1997) and supported by the broadband nature of the observed power spectra and multi-time-scale rms-flux relation observations (e.g. see discussion in Uttley et al. 2005) as well as (general relativistic) magnetohydrodynamic simulations of accretion flows (Hogg & Reynolds 2016; Bollimpalli et al. 2020). Since the simple impulse response functions described above map delay times on to unique radii, it is simple to model multilocation signals by using the same impulse responses but with different starting delays. Given a starting radius $r_{\mathrm{fluc}}$ we can use equation (6) to define a starting delay time $\tau_{\mathrm{fluc}}$, which is negative since we define the innermost delay to be zero. The impulse response $g_{\mathrm{fluc}}$ 'seen' by that fluctuating signal is then related to the overall impulse response ($g(\tau)$) calculated for that component by

$$g_{\mathrm{fluc}}(\tau) = \begin{cases} g(\tau) & \text{if } \tau \ge \tau_{\mathrm{fluc}} \\ 0, & \text{otherwise.} \end{cases} \quad (19)$$

If we then define $n_{\mathrm{fluc}}$ signals with Fourier power-spectra $P_{\mathrm{fluc},i}(\nu)$ (where $\nu$ is the Fourier frequency), originating at radii $r_{\mathrm{fluc},i}$ with corresponding time delays $\tau_{\mathrm{fluc},i}$, the resulting power-spectrum is given by

$$P(\nu) = \sum_{i=1}^{n_{\mathrm{fluc}}} P_{\mathrm{fluc},i}(\nu) |G_{\mathrm{fluc},i}(\nu)|^2, \quad (20)$$

where $G_{\mathrm{fluc},i}(\nu)$ is the Fourier transform of the impulse response for $\tau_{\mathrm{fluc},i}$ given by equation (19). The cross-spectrum between two components, labelled 1 and 2 (e.g. two energies in the power-law, or the disc versus a single power-law energy) is

$$C_{1,2}(\nu) = \sum_{i=1}^{n_{\mathrm{fluc}}} P_{\mathrm{fluc},i}(\nu) G_{\mathrm{fluc},i,1}(\nu) G_{\mathrm{fluc},i,2}^{*}(\nu) \quad (21)$$

with frequency-dependent time-lag given by the argument of the cross-spectrum (the phase lag) normalized by $2\pi\nu$.

For simplicity, we assume that the power spectra of the individual mass-accretion rate fluctuation signals $P_{\mathrm{fluc},i}(\nu)$ are described by Lorentzian functions parametrized by a peak frequency $v_{\mathrm{fluc},i}$, quality-factor $Q_i$ (approximately the frequency of the Lorentzian divided by its full-width at half maximum) and fractional rms rms$_i$. As discussed by Pottschmidt et al. (2003), the more common form of the Lorentzian function can then be calculated from these parameters, by determining the resonance frequency $v_{\mathrm{res},i} = v_{\mathrm{fluc},i} / \sqrt{1 + 1/(4Q_i^2)}$ and normalizing factor $R_i = \mathrm{rms}_i / \sqrt{0.5 - \tan^{-1}(-2Q_i)/\pi}$:

$$P_{\mathrm{fluc},i}(\nu) = \frac{2R_i^2 Q_i v_{\mathrm{res},i}}{\pi (v_{\mathrm{res},i}^2 + 4Q_i^2 (\nu - v_{\mathrm{res},i})^2)}. \quad (22)$$

The number and radial locations of Lorentzian signals and their parameters can be selected when setting up the model calculation and in principle can be arbitrary. However, the linear model described here does not account for viscous diffusion effects on mass accretion fluctuations as they propagate through the flow, which significantly suppresses variations on time-scales exceeding the local viscous propagation time-scale (e.g. see discussion in Churazov et al. 2001 and a more extensive treatment in Mushtukov et al. 2018). Therefore for simplicity, we will require that signal peak frequencies only correspond to time-scales which exceed the propagation delay through the flow from that radius.

The above equations and assumed input signal power spectra allow us to calculate the observed power spectra and lags expected from our model. It is important to note that this approach to combining signals from different radii assumes that the mass accretion rate variations combine additively instead of the multiplicative combination implicit in the propagating fluctuations model (Uttley et al. 2005). We will consider the effects of this simplification in Appendix B, by comparing our linear impulse response model predictions with predictions from numerical simulations of multiplicative propagating fluctuations.

## 4 NUMERICAL IMPLEMENTATION

To calculate the disc, seed and heating impulse responses, we define 200 contiguous, geometrically spaced radial bins of radii $r_i$ and widths $\Delta r_i$ between $r_{\mathrm{in}} = 2 R_g$ (corresponding to an ISO with dimensionless spin parameter $a = 0.94$) and $r_{\mathrm{fluc,max}}$, which corresponds to the outermost radius where mass-accretion fluctuations are generated in the disc. To account for the constant emission from the non-variable part of the disc, we further define a set of 200 contiguous, geometrically spaced disc radii from $r_{\mathrm{fluc,max}}$ to $r_{\mathrm{out}} = 1000 R_g$ which is used to calculate the constant parts of disc
and seed emission, which are included in the mean luminosities of those components, in addition to the contributions calculated from the impulse responses within $r_{fluc,max}$. For calculating the reverberation component of disc emission, we set $f_{therm} = 0.7$, corresponding approximately to 1 minus the disc albedo for X-ray reflection in the absence of soft X-ray reprocessing.

## 4.1 Mass accretion fluctuation time-scales and power-spectrum

We assume that mass accretion fluctuations are generated in both the disc and corona in each radial bin within $r_{fluc,max}$, i.e. there are $n_{fluc} = 200$ separate fluctuation signals. The fluctuation signals (see Section 3.4) have Lorentzian power spectra with peak frequencies $v_{fluc,i} = ((R_g/c)t_{fluc,i})^{-1}$ Hz and constant quality-factor $Q_i = 1$. The fluctuation time-scale $t_{fluc,i}$ in units of $R_g/c$ is given by

$$t_{fluc,i} = \begin{cases} s_i r_i^{(n_{t,d}+\frac{1}{2})}, & \text{if } r_{fluc,max} \ge r_i > r_{\text{cor}} \\ s_i f_{t,c} r_{\text{cor}}^{(n_{t,d}-n_{t,c})} r_i^{(n_{t,c}+\frac{1}{2})}, & \text{if } r_{\text{cor}} \ge r_i > r_{\text{in}} \end{cases} \quad (23)$$

where $s_i$ is a normalizing constant and $n_{t,d}$, $n_{t,c}$ are indices which correspond to the radial scaling of the fluctuation time-scale with respect to the local dynamical time-scale in the disc and corona respectively. The parameter $f_{t,c}$ is a scaling factor between time-scales of variability in the disc and corona. For an accretion flow described by viscosity parameter $\alpha$ (Shakura & Sunyaev 1973), $s_i r_i^{n_{t,d}} \propto \alpha^{-1}(r/h)^2$, so that (assuming constant $\alpha$) a radially constant aspect ratio $h/r$ corresponds to $n_t = 0$, while constant scale-height $h$ corresponds to $n_t = 2$. We choose this simple parametrization rather than make any physical assumptions about the radial dependence of $h$, because standard thin discs with $h \ll r$ are already known to be incompatible with propagation of mass-accretion fluctuations which could explain observed variability (Churazov et al. 2001). Therefore, we are implicitly assuming that hard-state inner, blackbody-emitting discs are not like geometrically thin, gas-pressure-dominated Shakura–Sunyaev discs and are somehow able to support rapid non-diffusive propagation of large variability signals. We discuss possible physical reasons for this difference in Section 6.2.

To obtain a reasonable match to observations, we fix $s_t = 1$ and $n_{t,d} = 2$, corresponding to a disc scale-height which is constant with radius. The main argument for a steep radial dependence of disc variability time-scale, i.e. $n_{t,d} = 2$, is that it better matches the observed energy-dependent low-frequency PL-D hard lags, which we demonstrate in Appendix C but do not consider further here, except to note that the spectral-timing properties discussed for $n_{t,d} = 2$ in Section 5 show similar properties for $n_{t,d} = 0$ or 1. We assume $f_{t,c} = 1$ which corresponds to the disc and corona fluctuation time-scales being equal at $r_{cor}$. We assume that the coronal time-scale scaling index is given by $n_{t,c} = \max\{0, n_{t,c}(t_{min})\}$, where $t_{min}$ is a time-scale which must equal or exceed the time-scale of variability fluctuations generated at $r_{in}$, from which $n_{t,c}(t_{min})$ can be calculated by rearranging the lower equation in equation (23). With this constraint, we limit the coronal scaling index $n_{t,c} \ge 0$ (i.e. the limiting case is a flow with constant aspect ratio), and the fastest variability time-scale produced in the flow to be $t_{fluc,1} \lesssim t_{min}$. Note that we also do not require the coronal variability time-scales to be coupled to coronal height or other aspects of the geometry apart from the radius, since the vertically extended coronal region which upscatters disc photons may be mostly distinct from the hot accretion flow which powers it.

The shapes of broadband BHXRB power spectra (e.g. Heil, Uttley & Klein-Wolt 2015) lead us to assume a natural form for the mass accretion fluctuations which is a doubly-bending power law with $P(\nu) \propto \nu^{-1}$ in the intermediate frequency range (i.e. flat in $\nu P(\nu)$). This shape can be approximately reproduced by weighting the squared fractional rms of each Lorentzian by the logarithmic change in peak frequency to the next radial bin<fnref n="4" />:

$$rms_i^2 = rms_{tot}^2 \frac{\log(v_{fluc,i}/v_{fluc,i+1})}{\log(v_{fluc,1}/v_{fluc,n_{fluc}})} \quad (24)$$

where rms$_{tot}$ is the total fractional rms of the mass accretion fluctuation signals, which we assume to be 40 per cent. The magnitude of rms is important for considering non-linear effects, both due to power-law photon index variations and the multiplicative nature of mass accretion variations implied by the rms-flux relation. However, in Appendix B we use numerical simulations to show that these effects remain negligible for accretion fluctuations with integrated rms of 40 per cent, consistent with the X-ray variability observed in hard state BHXRBs.

We determine the bend-frequencies of the mass-accretion fluctuation power-spectrum by setting the maximum and minimum fluctuation time-scales such that the frequency at the maximum fluctuation radius $r_{fluc,max}$ is $v_{fluc,max} = ((R_g/c)t_{fluc,i=n_{fluc}})^{-1}$ Hz and at the minimum fluctuation radius (terminating at $r_{in}$), $v_{fluc,min} = ((R_g/c)t_{fluc,i=1})^{-1}$ Hz. For the calculations shown here we choose $v_{fluc,max} = 0.03$ Hz (corresponding to $r_{fluc,max} = 46 R_g$ for $s_t = 1$, $n_{t,d} = 2$), and set $t_{min}$ to correspond to a frequency of 100 Hz, which sets the maximum fluctuation frequency to equal or exceed this value. Together, these choices will result in power-spectra similar to those seen in BHXRBs for at least part of the hard state.

## 4.2 Propagation delays

We assume that the propagation delay across a given radial bin $\Delta\tau_i$, scales with the local variability time-scale in that bin. This behaviour is expected if fluctuation time-scales scale with the local viscous time-scale, as expected in the propagating fluctuations scenario (Lyubarskii 1997; Churazov et al. 2001; Hogg & Reynolds 2016; Mushtukov et al. 2018). Based on equation (6) and incorporating additional scaling factors for propagation in the disc and corona, $f_{\Delta\tau,d}$ and $f_{\Delta\tau,c}$, we obtain the following relations for the propagation delay across a radial bin:

$$\Delta\tau_i = \begin{cases} s_t f_{\Delta\tau,d} r_i^{(n_{t,d}+\frac{1}{2})} \Delta r_i, & \text{if } r_{\text{fluc,max}} \ge r_i > r_{\text{cor}} \\ s_t f_{\Delta\tau,c} r_{\text{cor}}^{(n_{t,d}-n_{t,c})} r_i^{(n_{t,c}+\frac{1}{2})} \Delta r_i, & \text{if } r_{\text{cor}} \ge r_i > r_{\text{in}} \end{cases} \quad (25)$$

The resulting delays are used to calculate the impulse response functions following the approach outlined in Sections 3.1 and 3.2. Our default assumption is that $f_{\Delta\tau,d} = 1$ and $f_{\Delta\tau,c} = 0.1$. With these scalings and those of the disc and corona variability time-scales, we expect a smooth evolution of variability time-scale from disc to corona, but a transition to a fast accretion time-scale within the corona, commensurate with the expected geometrically thick and optically thin hot flow.

<footnote n="4">⁴For this purpose, we also include the frequency calculated using equation (23) for the non-fluctuating bin at $i = n_{fluc} + 1$.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 3. Top panel: Assumed radial dependence of fluctuation time-scale and total propagation delay of signals from that radius, for the fluctuation parameters described in Section 4. Lower panel: Power spectrum obtained by summing the Lorentzian power spectra of fluctuation signals from each radial bin, equivalent to the power spectrum of mass accretion fluctuations which would arrive at the innermost radius in the limit where those fluctuations are additive. In both cases, the split between the disc and coronal radii is indicated by a vertical grey dotted line.</figcaption></figure>

## 4.3 Summary of fluctuation timing properties

The radial dependencies of signal fluctuation time-scale and total propagation delay for the model assumptions described in this section are shown in the top panel of Fig. 3. Note that the propagation delay is always shorter than the variability time-scale, consistent with our assumption that viscous damping effects are small. The summed power spectrum of mass accretion fluctuations (obtained by summing the Lorentzian power spectra of fluctuations from each radius) is shown in the lower panel. Given our assumptions of a smooth join in variability time-scale from disc to corona, and constant fractional rms of accretion fluctuations produced across the variable region of both disc and corona, the power spectrum is also smooth with no feature linked to the transition from disc to corona. The observed energy-dependent steepening of BHXRB X-ray power spectra above 1 Hz is then caused by the radially extended nature of the disc and coronal emission (see Section 5 and also Arévalo & Uttley 2006).

## 4.4 Spectral-timing calculations

To enable rapid calculation of the Fourier transforms of impulse responses required for determining spectral-timing products, we approximate our impulse response functions to be sums of shifted top-hat functions, which correspond to each radial delay bin. That is, following equation (19), the continuous-time impulse response applied to fluctuations produced in the jth radial bin is

$$g_{\text{fluc},j}(\tau) = \sum_{k=1}^{j} g(\langle\tau\rangle_k)\text{rect}(\tau/\Delta\tau_k) \times \delta(\tau - \langle\tau\rangle_k), \quad (26)$$

where $g(\langle\tau\rangle_k)$ is the impulse response value for the bin with mean delay $\langle\tau\rangle_k = (\sum_{l=1}^k \Delta\tau_l) - \Delta\tau_k/2$, $\text{rect}(\tau/\Delta\tau_k)$ is the normalized

rectangular function centred on zero with total width $\Delta\tau_k$ and the product of these quantities is then time-shifted to $\tau = \langle\tau\rangle_k$ by convolution with a $\delta$-function. The Fourier transform of a time-shifted top-hat function is a phase-shifted sinc-function, so the Fourier transform of the impulse response for a signal originating in the jth radial bin can be calculated for any given frequency $\nu$, using

$$G_{\text{fluc},j}(\nu) = \sum_{k=1}^{j} g(\langle\tau\rangle_k) \Delta\tau_k \text{sinc}(\Delta\tau_k \nu) \exp(-i2\pi\nu\langle\tau\rangle_k), \quad (27)$$

where i here denotes the unit imaginary number. With this approach, we can use equations (20) and (21) to rapidly calculate the frequency-dependent spectral-timing properties for our chosen frequencies. We use a geometrically spaced grid of 500 frequencies over the range 0.003–200 Hz.

# 5 RESULTS

We calculate the power-law flux impulse responses for soft, medium, and hard photon energies $E_{\text{soft}} = 3E_s$, $E_{\text{med}} = 9E_s$, and $E_{\text{hard}} = 27E_s$. The power-law lags for energies defined in this way are independent of the choice of seed energy $E_s$. However, for a seed energy corresponding to the peak of a disc blackbody with maximum temperature $kT_{\text{max}}$, $E_s = 2.82kT_{\text{max}}$ so that the relevant energies for a hard state disc with $kT_{\text{max}} = 0.2$ keV are $E_s \simeq 0.56$ keV, $E_{\text{soft}} \simeq 1.7$ keV, $E_{\text{med}} \simeq 5.1$ keV, and $E_{\text{hard}} \simeq 15.2$ keV. To compare the power spectra, we normalize them by the summed mass accretion rate signal power spectrum, $P_{\sum m(r)} = \sum_{i=1}^{n_{\text{fluc}}} P_{\text{fluc},i}(\nu)$. By normalizing in this way, we remove some of the effects of the somewhat arbitrary choice of driving signals which dominate the power-spectral shape, while maintaining the energy-dependent effects which are a key outcome of the model.

## 5.1 Spectral-timing properties for a disc and spherical corona with $r_{\text{cor}} = 10 R_g$

To demonstrate and explain the X-ray spectral-timing properties predicted by our model, we first consider a spherical corona with radius $r_{\text{cor}} = 10 R_g$, with default fluctuation signal and propagation delay parameters given in Section 4. This coronal geometry predicts a mean photon index $\langle\Gamma\rangle = 1.47$, consistent with that seen in hard state black hole X-ray binaries. The frequency-dependent lags and power-spectra for different power-law energies (in terms of $E_s$) and the disc luminosity $L_d$ are shown by the solid lines in the upper-left panels of Fig. 4. The model shows a strong correspondence to a number of the key spectral-timing properties of hard state BHXRBs which we described in Section 1:

(i) At low frequencies, disc variations substantially lead power-law variations at all energies considered showing hard PL-D lags, but cross to large (~10 ms) soft PL-D lags above ~1 Hz, when measured with respect to the lower energy power-law photons. The high-frequency disc lags are much weaker relative to higher energy power-law photons, approaching zero then becoming positive (power-law lagging) as the power-law energy increases.

(ii) Power-law photons of different energies show frequency-dependent hard PL-PL lags, with a weaker and less complex frequency dependence than seen for PL-D lags.

(iii) Lags between power-law energies scale equal to or close to log-linearly with energy, over a broad frequency range (i.e. the same factor change in energy shows the same lag).
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Model spectral-timing properties and their relation to impulse response functions, for a spherical corona with radius $r_{cor} = 10 R_g$. Top left: Calculated frequency-dependent time-lag and power-spectral ratio to input fluctuation power spectrum, for the disc luminosity and several power-law energies. The solid lines show the full model calculations, while dotted lines show an approximation based on the impulse response centroid delays and emission fractions versus fluctuation signal frequency (see the text for details). The frequency dependences of these quantities are shown in the lower-left panel. The right panel shows the disc, seed and coronal heating impulse responses for this configuration, showing from top to bottom the part of the impulse responses (shaded) which responds to a signal at a given radius and corresponding frequency and delay (vertical dashed lines). The centroid delay values for each of the starting signal radii are shown as solid vertical lines in the colour corresponding to the component (disc, seed, and heating). The delays are plotted using symmetric logarithmic axes with delays > −100 $R_g/c$ plotted on a linear scale. The response units are luminosity per unit time delay, where luminosities are expressed as a fraction of the total dissipated accretion luminosity. The impulse responses are further multiplied by the absolute time delay, so that the delays which contribute the most to the impulse response appear as a peak.</figcaption></figure>

(iv) Variability at high-frequencies is suppressed and the amount of suppression is energy-dependent, with lower-energy power-law emission showing less high-frequency power than higher energy emission. Disc emission variability is strongly suppressed down to lower frequencies than the power-law emission.

To understand how these spectral-timing properties arise, consider the disc, seed and coronal heating impulse responses that are ‘seen’ by signals originating at different radii ($r_{fluc}$) in the accretion flow, with fluctuation frequency $v_{fluc}$ (see Section 3.4). These are illustrated in the right panel of Fig. 4 by shaded regions, which extend to the right of the delay time $\tau_{fluc}$ at the signal radius of origin, marked by vertical dashed lines. The impulse responses produced outside the coronal radius (illustrated by the vertical dot-dashed line) are generated primarily by viscous dissipation in the disc, with a weak reverberation component due to disc seed luminosity returning from the corona to

heat the disc. The impulse responses produced inside the coronal radius are powered by viscous dissipation in the corona, which produces coronal heating and in turn leads to a pure-reverberation component of disc and seed emission.

To better understand how these impulse responses lead to the complex spectral-timing behaviour predicted by the model, it is useful to make some simplifying approximations. For signals with variability time-scales which are long compared to the emission-weighted (centroid) delay $\langle \tau \rangle$, as is generally the case here, the time-lag between two emission components (denoted 1 and 2), $\tau_{1,2}(\nu) \simeq \langle \tau_2 \rangle - \langle \tau_1 \rangle$. The centroid delay is calculated using

$$ \langle \tau_* \rangle = \frac{\int_0^{\tau_{fluc}} \tau g_*(\tau) d\tau}{\int_0^{\tau_{fluc}} g_*(\tau) d\tau}, \qquad (28) $$

where * symbolizes the component, i.e. d(isk), s(ed), and h(eating). The seed and heating components are not directly observed but
contribute to the power-law emission, to give the power-law emission centroid delay:

$$ \langle \tau_{\text{pl}}(E) \rangle = \frac{[1 - u(E)] f_s \langle \tau_s \rangle + u(E) f_h \langle \tau_h \rangle}{[1 - u(E)] f_s + u(E) f_h}, \quad (29) $$

where $f_s$ and $f_h$ correspond to the 'emission fraction' of the total luminosity for that impulse response component, i.e. equivalent to the shaded areas in the right panels of Fig. 4:

$$ f_* = \frac{\int_0^{\tau_{\text{fluc}}} g_*(\tau) d\tau}{\langle L_* \rangle} \quad (30) $$

with power-law emission fractions given by $f_{\text{pl}}(E) = [1 - u(E)]f_s + u(E)f_h$. For time-scales which are long compared to the centroid delay, the observed variability power of a fluctuation signal scales with the intrinsic power multiplied by $f_*^2$ or $(f_{\text{pl}}(E))^2$.

The numerical values of the disc and power-law centroid delays and the squared emission fractions are shown as a function of fluctuation frequency in the lower-left panel of Fig. 4. To compare our approximate predictions for lag and power-spectral amplitude with the full model calculations, we plot the fluctuation frequency-dependent centroid delay differences and squared emission fractions as dotted lines in the upper-left panel. For our simplified calculations, we do not combine signals with overlapping power-spectra and simply assume that a signal dominates at its respective fluctuation frequency. Despite this approximation, the simplified calculation results are quite close to those obtained with the full Fourier treatment of combined signals.

For more extended impulse responses such as for the disc emission, the lags and relative power-spectral amplitude both decrease strongly with increasing fluctuation frequency due to the lag centroid (relative to the power-law emission) and emission fraction changes. The power-law emission on the other hand is a weighted combination of two more centrally concentrated impulse response components (seed and heating), leading to smaller lags and weaker evolution of lagas and power-spectra with fluctuation frequency, except at higher frequencies when the fluctuation signals are fully embedded in the regions producing seed and heating components. The energy-dependent weighting of the seed and heating components leads to more extended power-law impulse responses at lower energies where the seed emission contributes more strongly than heating, so that suppression of the high-frequency power-spectrum decreases with photon energy, as observed. For lower fluctuation frequencies the PL-PL lag amplitude scales linearly with the difference in $\log(E)$, also as observed (note that the lags are identical for the same ratio of energies considered). This effect arises as $f_s, f_h \to 1$ for larger fluctuation signal radii. In that limit, the lag for power-law energies $E_1$ and $E_2$ (calculated using the energy-dependent centroid delays obtained with equation 29) simplifies to

$$ \langle \tau(E_2) \rangle - \langle \tau(E_1) \rangle \approx \beta \langle \Gamma \rangle (\langle \tau_h \rangle - \langle \tau_s \rangle) \ln(E_2/E_1). \quad (31) $$

Therefore, the power-law emission lags scale linearly with the difference between heating and seed centroid delays (thus producing the observed positive, i.e. hard lags), the mean photon index, and the difference in $\log-E$. The latter dependence arises naturally from the power-law pivoting predicted by the model.

Our model also reproduces the negative PL-D lags that have previously been attributed to reverberation light-travel delays (Uttley et al. 2011). The PL-D negative lags arise primarily because the illumination of the disc by the corona produces a strong peak of disc reverberation emission at the smallest delays, which is not replicated for the seed emission, since only a small fraction of disc reverberation photons return to the corona. The combination of these

effects causes the disc centroid for small signal radii to become smaller (less negative) than the seed and soft power-law emission centroids, so that the power law versus disc lag becomes negative at higher Fourier frequencies, as observed. This relative shift in the disc and seed centroid delays can be seen in the top-right and middle-right panels of Fig. 4, as the fluctuation signal shifts from 0.3 Hz to 3 Hz.

These 'soft' lags are primarily due to the reverberation component of disc emission, but they are much larger than any expected light travel delay, as they are linked to the propagation of fluctuations through the inner disc radii which provide most of the seed photons. The disc lags disappear when comparing with higher energy power-law emission, also as observed, because of the increased effect of heating of the corona at those energies, which has an impulse response concentrated at small delays, similar to the reverberation component which is also driven by the coronal heating luminosity.

## 5.2 Effects of coronal geometry on impulse responses and spectral-timing properties

The coronal geometry has a significant impact on the predicted lags, for several reasons. First, the coronal visibility from the disc determines the amount of seed emission intercepted from different disc radii. A more vertically extended corona will intercept proportionately more photons from fluctuations further out in the disc, to produce a seed impulse response that extends to larger negative delays. Secondly, coronae that are more visible from the disc will illuminate the disc more strongly to produce stronger reverberation components in the disc and seed impulse response at small delays. The coronal radius also determines the distribution of accretion power between disc emission (and hence seed luminosity) and coronal heating, which affects the mean photon index as well as the relative strength of reverberation versus propagation components in the disc and seed impulse responses.

To examine the effect of coronal geometry on spectral-timing properties, we will consider two types of coronal geometry: spherical (with radius $r_{\text{cor}}$) and an inverted cone geometry. The inverted cone shape (strictly speaking an inverted conical frustum, i.e. a cone with the pointed end cut-off) is intended to represent a jet or wind-like corona or a cylindrically shaped 'hot flow'. It has a basal coronal radius in the disc plane $r_{\text{cor}}$ and coronal height above the disc plane $h_{\text{cor}}$. The radius at the top (widest) end of the cone is determined by the opening angle $\theta_{\text{cor}}$ which is the angle relative to the normal to the disc plane, such that $\theta_{\text{cor}} = 0^{\circ}$ gives a cylindrical corona. Unlike the spherical corona, the inverted cone geometry allows us to study the effect of changes in coronal visibility as seen from the disc, independent of changes in coronal radius.

We consider the effects of changing the coronal radius (for both the spherical and inverted-cone geometries), and the coronal height and opening angle (for the inverted-cone geometry). Fig. 5 shows the disc (solid lines), seed (dotted lines) and heating (dashed lines) impulse responses for these situations, with both the dissipated and reverberation contributions to the disc and seed components summed to give those impulse responses. To understand the differences in impulse responses between the different geometries, it will be useful to consider the visibility of the corona from the disc for the different geometries, which we can show by plotting the radially dependent fraction of disc photons intercepted by the corona, $f_{d \to c}(r)$, in Fig. 6.

Fig. 7 shows the resulting frequency-dependent time lags and power spectra for disc and power-law emission, for the different coronal geometries with impulse responses shown in the corresponding panels in Fig. 5. The lags are plotted for the soft power-law flux ($3E_s$) versus disc luminosity and hard ($27E_s$) versus medium ($9E_s$) and
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparison of impulse responses for different coronal geometries and parameter variations. **Top-left:** Changing coronal radius $r_{cor}$ for a spherical corona. Then clockwise from **top-right:** Inverted cone corona with changing coronal radius $r_{cor}$, changing opening angle $	heta_{cor}$, changing coronal height $h_{cor}$. For each changing parameter value the plots show the total disc (solid lines), seed (dotted lines), and heating (dashed lines) impulse responses. Vertical dashed lines show the radii and fluctuation signal frequencies associated with a given delay. For the top panels with varying $r_{cor}$, these radii correspond to different delays and are colour-coded accordingly. See Fig. 4 for further details.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Comparison of the radially dependent fraction of disc photons intercepting the corona $f_{d \to c}(r)$, for a spherical corona (solid lines) and an inverted cone corona (dashed, dot-dashed, and dotted lines) with different radii (left panel), heights (middle panel), and opening angles (right panel).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Comparison of spectral-timing properties, corresponding by position to the different coronal geometries with impulse responses shown in Fig. 5. For each plot, the top part of each panel shows time lags for power law at *E*<sub>soft</sub> versus disc (solid line), power-law at *E*<sub>med</sub> versus *E*<sub>soft</sub> (dashed), power-law at *E*<sub>hard</sub> versus *E*<sub>med</sub> (dotted). Positive lags indicate that the harder band lags the softer band (or the power-law lags the disc). The lower part of each panel shows the ratio of emission power spectrum to integrated *m* power spectrum for disc (solid line) and power-law flux at *E*<sub>med</sub> (dashed) and *E*<sub>hard</sub> (dotted). The signal peak frequencies at *r*<sub>cor</sub> are shown by the vertical dot-dashed lines.</figcaption></figure>

medium versus soft bands. For clarity, we only show power spectra for the disc and medium and hard power-law bands, since these latter bands can be observed free of disc emission in BHXRBs. We note for completeness that the soft power-law power spectrum continues the trend of greater high-frequency filtering at softer energies (see Fig. 4).

The largest changes in impulse response and corresponding lags (PL-D and different PL-PL combinations) are associated with changes in coronal radius for the inverted cone as well as spherical coronae. Increases in coronal radius naturally shift the inner edge of disc and seed impulse responses associated with disc propagation to

larger negative delays. The relative shift in disc and seed emission to larger radii produces greater suppression of high-frequency variability, which is seen in the power spectra. The resulting shift in the centroid delays increases the value of low-frequency PL-PL hard lags (as the seed impulse response is pushed to larger negative delays), but does not substantially change the PL-D hard lags at the lowest frequencies, since both disc and seed centroids are shifted in the same direction. However, the PL-D soft lags increase substantially in amplitude and shift to lower frequencies, thus decreasing the ‘crossing frequency’ from hard to soft lags (and corresponding ‘cut-off’ in the hard lags). This increase and shift in the soft lags is
due to the change in seed centroid and because for larger coronae, reverberation dominates the disc centroid at lower signal frequencies than for coronae with smaller radii.

The differences between spherical and inverted cone coronae can be understood in terms of their visibilities as seen from the disc, which are shown in the left panel of Fig. 6. The spherical coronae generally show smaller solid angles compared to the inverted cone coronae, which have a large surface 'facing' the disc. Thus, spherical coronae show systematically weaker seed emission than inverted cone coronae of the same radius, as well as relatively weaker 're-turning' seed reverberation components for delays produced inside the corona. The latter effect leads to more negative seed centroid values for large spherical coronae compared to large inverted cone coronae, which leads to larger PL-PL hard lags and PL-D soft lags (and lower crossing-frequencies) for large spherical coronae.

Comparing the effects of changing height and opening angle for inverted cone coronae with the same coronal radius is relatively simple and can be understood entirely in terms of the changing coronal visibility from the disc. Increases in coronal height lead to significant extensions of the propagation component of the seed impulse response to large delays, while the reverberation components do not change substantially. This pattern arises because the largest changes in visibility occur at large disc radii (see the middle panel of Fig. 6), which contribute the largest delays, while there is relatively minor change at small disc radii, which contribute most to reverberation. Thus, the seed centroid becomes substantially more negative with increasing coronal height, PL-PL hard lags and PL-D soft lags both increase, and the PL-D crossing-frequency decreases. At the same time, there is only minor change in the power-spectra, because the large change in seed centroid delay is mainly caused by a relatively small change in seed emission at very large propagation delays.

In contrast to height changes, coronal opening angle changes produce relatively little effect on the spectral-timing properties. A change in opening angle affects the coronal visibility from the disc in a similar way over a broad range of radii (see the right panel of Fig. 6). Thus, the seed impulse responses increase in amplitude with increasing opening angle, but the centroid delay increase is relatively small. The notable exception is in the high frequency power-law emission power spectra that show greater suppression for small opening angles. This effect is driven by the relatively large impact of coronal opening angle on the seed returning emission due to reverberation, which provides the only source of seed photons that can respond to high-frequency signals originating in the coronal part of the accretion flow.

## 5.3 Evolution of spectral-timing 'observables' with coronal geometry

Finally, it is useful to quantify and visualize the evolution of certain spectral-timing properties which approximate the properties which can be seen in observational spectral-timing and spectral data. These are $\langle\Gamma\rangle$, the time-averaged photon index; $\tau_{max,PP}$, the maximum (low-frequency) PL-PL hard lag (fluxes for $E_{hard}$ versus $E_{med}$, although the same lag is obtained for $E_{med}$ versus $E_{soft}$); $\tau_{min,PD}$, the minimum (or maximum negative, i.e. soft) PL-D lag (power-law flux at $E_{soft}$ versus disc luminosity); $\nu_{cr,PD}$, the zero-lag crossing frequency for PL-D lags; and finally $P_M/P_H(\nu = 5 \text{ Hz})$, the ratio of power at $E_{med}$ to that at $E_{hard}$, as measured at a frequency of 5 Hz.

Fig. 8 shows these quantities obtained from the frequency-dependent lags and power-spectra calculated for the spherical and inverted cone geometries. For the spherical corona, we calculate for 36 values of $r_{cor}$ spaced geometrically from 5 to 30 $R_g$ (corresponding

to edges of the disc radial bins). For the cone geometry, we calculate for each of three values of $r_{cor}$ (6, 10, and 20 $R_g$), 160 combinations of $h_{cor}$ and $\theta_{cor}$ evaluated along gridlines of fixed $h_{cor} = \{5, 10, 20, 30\} \, R_g$, and $\theta_{cor} = \{0, 15, 30, 45\}$ degrees. The coronal parameter values corresponding to set points or vertices on each grid line are labelled in the plot. For the inverted cone coronae, we label for clarity only those vertices that are clearly visible and for the $r_{cor} = 6 \, R_g$ and $r_{cor} = 20 \, R_g$ cases only: the corresponding grid points for $r_{cor} = 10 \, R_g$ can be inferred from those shown for the other coronal radii.

The time-averaged photon index $\langle\Gamma\rangle$ depends on the ratio of seed to coronal heating luminosity (equation 1) and thus decreases with increasing coronal radius $r_{cor}$, which changes the balance of viscous dissipation from disc emission and hence seed photons, towards coronal heating. Increases in coronal opening angle increase the seed photon flux substantially, while increases in coronal height only increase the seed flux marginally (from larger radii), hence $\theta_{cor}$ is strongly correlated with $\langle\Gamma\rangle$ while $h_{cor}$ is only weakly correlated with $\langle\Gamma\rangle$.

For the reasons described earlier in this subsection, the lags and crossing-frequency depend most strongly on coronal radius and height but not on opening angle. The strong dependencies on coronal radius and height result in strong (anti-)correlations between the minimum PL-D lag, crossing-frequency and maximum PL-PL lag, which can be seen in the lower two panels of Fig. 8. The energy-dependent power-spectral evolution is more complex however, showing a maximum effect for intermediate coronal radii (between 10 and 15 $R_g$). This maximum effect on energy dependence of the power spectra at 5 Hz seems to correspond to the point where this frequency corresponds to the signal frequency at $r_{cor}$. For signal frequencies produced within $r_{cor}$, coronal heating and seed impulse responses track each other closely, since the seed emission for delays produced within the corona is produced by reverberation driven by the coronal heating power. Therefore, the medium and hard band power spectra stop diverging in shape at signal frequencies produced within the corona, since they depend on an energy dependent combination of emission fractions $f_h$ and $f_s$.

## 5.4 Dependence on variability and propagation time-scale parameters

For the previous analyses, we have fixed the disc and coronal variability time-scale normalization and radial scaling index to $s_t = 1$ and $n_{t,d} = 2$, respectively (see equation 23) and the disc and coronal propagation time-scales to be equal to 100 per cent and 10 per cent of the variability time-scale, respectively (corresponding to $f_{\Delta\tau}$, $d = 1$ and $f_{\Delta\tau}$, $c = 0.1$, see equation 25). These values are chosen to replicate reasonably well the observed spectral-timing properties (see previous subsections and Appendix C). However, it is interesting to examine the effect of the time-scales for accretion flow variability and propagation on the observed spectral-timing properties.

The left panel of Fig. 9 shows the effects of changing the normalization of disc variability time-scale and propagation time-scale relative to the default situation (shown in blue), while maintaining fixed maximum fluctuation frequency (and power-spectral break), $\nu_{fluc,max} = 0.03 \text{ Hz}$. The main effect of reducing $s_t$ (orange curves) is to reduce (proportionate with $s_t$) the maximum PL-PL hard lags and PL-D soft lags, while increasing the zero-crossing frequency of PL-D lags and the frequencies of the roll-overs seen in the PL-PL hard lags and the power spectra. These changes occur because these spectral-timing features are produced over a range of radii which are fixed by the coronal geometry. The variability and propagation
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Dependence of spectral-timing observable quantities on coronal geometry for a spherical corona and inverted cone coronae (with $r_{\text{cor}} = \{6, 10, 20\} R_g$). The observable quantities are $\langle\Gamma\rangle$, the time-averaged photon index; $\tau_{\text{max,PP}}$, the maximum (low-frequency) PL-PL hard lag ($E_{\text{hard}}$ versus $E_{\text{med}}$); $\tau_{\text{min,PD}}$, the minimum (or maximum negative, i.e. soft) PL-D lag; $\nu_{\text{cr,PD}}$, the zero-lag crossing frequency for PL-D lags and $P_M/P_H(\nu = 5 \text{ Hz})$, the ratio of power at $E_{\text{med}}$ to that at $E_{\text{hard}}$ at $\nu = 5 \text{ Hz}$. Coronal parameter values are labelled for the spherical and largest/smallest inverted cone coronae. See the lower-left panel for the legend for all panels.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Left panel: Comparison of spectral-timing properties for the inverted cone corona with different scaling constants for accretion flow variability and propagation time-scales, assuming disc variability time-scale radial-scaling index $n_{t,d} = 2$ (see equations 23 and 25). The dot-dashed lines show the frequencies at $r_{\text{cor}}$ for the two cases of $s_t$. Right panel: Comparison of spectral-timing properties for different $n_{t,d}$. To ensure that fluctuation frequencies are identical at $r_{\text{cor}} = 10 R_g$ (dot-dashed lines), we set $s_t = \{100, 10, 1\}$ for $n_{t,d} = \{0, 1, 2\}$.</figcaption></figure>

time-scales associated with these fixed radii simply scale with $s_t$. In contrast the PL-D lags at low frequencies are produced at radii which increase as $s_t$ is reduced (to maintain the fixed $v_{\text{fluc,max}}$, so they remain relatively independent of $s_t$.

In contrast, reducing the propagation time-scales relative to $s_t$ (green curves), by reducing $f_\lambda \tau$, $d$ and $f_\Delta \tau$, $c$, leads to a proportionate reduction in the amplitudes of all lags while largely maintaining the crossing-frequency of PL-D lags and the roll-over frequencies of PL-PL lags and all power spectra. Due to our requirement that propagation time-scales in the disc should be similar to the variability time-scale, it is unlikely that significant changes of propagation time-scale relative to variability time-scale will occur in real accretion flows. However, significant changes of variability time-scale could be associated with changes in $h/r$ (or possibly viscosity parameter $\alpha$) and could therefore produce evolution of spectral-timing behaviour similar to that expected from changing $r_{\text{cor}}$ (e.g. see Fig. 7, upper panels), but without the concomitant changes in $\Gamma$.

The right panel of Fig. 9 shows the spectral-timing properties for an inverted cone corona with three different values of the radial scaling index of variability time-scale $n_{t,d}$, assuming the same fluctuation frequencies at $r_{\text{cor}}$. For standard accretion disc variability time-scales, we recall here that $n_{t,d} = 0$ corresponds to a constant aspect ratio $h/r$ of the accretion flow with radius, while $n_{t,d} = 2$ corresponds to constant $h$. The overall spectral-timing properties depend only weakly on $n_{t,d}$, when compared to changes in coronal radius and height (e.g. see Fig. 7). This means that the model results shown here are qualitatively robust to changes in $n_{t,d}$, but also implies that it may be difficult to determine the correct value of $n_{t,d}$ from observations. However, the range of disc variable disc radii changes dramatically with $n_{t,d}$, with maximum fluctuation radius $r_{\text{fluc,max}} = \{355, 86, 46\} R_g$ for $n_{t,d} = \{0, 1, 2\}$, assuming a 0.03 Hz minimum fluctuation frequency for this example. These differences have strong implications for the energy dependence of the disc impulse responses and corresponding spectral-timing properties, which we discuss in Appendix C.

# 6 DISCUSSION

We have shown that a model where mass accretion fluctuations propagate through the accretion disc to a central and relatively compact corona can explain a wide range of the observed hard state spectral-timing properties from accreting black holes. To enable this match to observations, the effects on the coronal power-law spectrum of variable disc seed photons and coronal heating should be included, as well as reprocessing of coronal emission by the disc. The variability is produced by a plausible mechanism, namely turbulence in the accretion flow which leads to propagating mass-accretion fluctuations, with longer time-scale variations being generated at larger radii. The latter process is strongly implied by other observable properties such as the rms-flux relation (Uttley et al. 2005) and the broad-band nature of the aperiodic variability (Lyubarskii 1997; Ingram & Done 2011), and appears naturally in (GR)MHD simulations of accretion flows (Hogg & Reynolds 2016; Bollimpalli et al. 2020).

Beyond these main assumptions, no special features are required for the model to explain many observed spectral-timing properties of BHXRBs along with changes in those properties, e.g. due to changing coronal geometry. We consider it a major strength of the model that it is both physically simple and consistent with the current disc plus inner corona picture of the innermost regions of accreting black holes. To demonstrate this for the MAXI J1820 + 070 hard state data set shown in Fig. 1, we show in Appendix D that the time-averaged decomposition of disc and coronal luminosity in that observation are broadly consistent with the disc truncation radii and coronal geometries considered by our lags model. If hard state spectral energy distributions (SEDs) are consistent with disc truncation radii of ~20 $R_g$ or less, it is likely that the HIMS and SIMS are consistent with successively smaller truncation radii which overlap the $r_{\text{cor}}$ values assumed in Section 5.

The disc-corona picture considered by our model is supported by a wide range of observational and theoretical results, even if the nature
of the corona is still to be determined (Done et al. 2007; Ingram & Motta 2019; Kalemci, Kara & Tomsick 2022). We further suggest that to prevent the significant spectral-timing effects shown here, we would need to discard one or more of these assumptions which underpin our understanding of other key aspects of BHXRB spectra and variability in luminous hard states. The spectral-timing features we explain here can thus be seen as a natural consequence of the response of the disc-corona system to mass-accretion fluctuations propagating through the disc.

We now consider some of the key implications of our results for the lags predicted by Comptonization and reverberation models, the nature of the variable accretion disc and the evolution of coronal geometry through the different spectral states. We will end with some discussion of possible improvements and extensions to the model.

## 6.1 Light-travel delays due to Comptonization and reverberation

In Section 2, we showed how a time-delay between correlated variations in coronal seed photon and heating luminosities leads naturally to energy-dependent time lags between flux variations at different power-law energies. The lag scales with the difference in centroid delays of the heating and seed impulse response (equation 31), so any physical scenario where there is a delay between these components should lead to a large lag in addition to that caused by other effects such as Compton scattering light-travel delays. An obvious corollary of this effect is that models which seek to explain power-law lags exclusively with Compton scattering delays (e.g. Kazanas et al. 1997; Reig et al. 2003; Karpouzas et al. 2021; Bellavita et al. 2022) must include a plausible mechanism to explain why coronal heating and seed variations should be coincident in time.<fnref n="5" />

Our model shows that Comptonization delays from large coronae are not needed to explain the observed PL-PL hard lags, which can be explained by much smaller coronae, and in fact seem to require relatively small coronae, or the expected lags would be much larger than observed. Such small coronae should produce Comptonization delays < 1 ms, which could become significant at high Fourier frequencies where propagation-induced delays become small.

Reverberation light-travel delays are not included in our model, where we assume them to be negligible compared to the propagation-related delays. This assumption is consistent with the compact coronae that are required to explain the data in our model, but such delays must be present, even if they are dominated by the expected propagation-related hard and soft lags over most of the observable frequency range. Fig. 7 shows that the propagation-related soft lags go to zero at frequencies of tens of Hz or higher. This arises because the seed and reverberation components are directly powered by coronal heating for impulse response delays produced inside the corona, resulting in zero lag between the power-law and reverberation emission. This result suggests that ‘clean’ light-travel reverberation signals are probably best-observed at high frequencies, ~ 100 Hz, where any seed photon variations should be driven primarily by reverberation. Because of signal-to-noise limitations at high frequencies, such light-travel reverberation signals may be inaccessible to current X-ray spectral-timing instruments such as

NICER but they will be detectable with proposed large-area missions such as the enhanced X-ray Timing and Polarimetry (eXTP) mission (De Rosa et al. 2019; Zhang et al. 2019) and the Spectroscopic Time-Resolving Observatory for Broadband Energy X-rays (STROBE-X; Ray et al. 2018).

## 6.2 The nature of the accretion disc

The spectral-timing properties of hard state black hole X-ray binaries can be explained if the variability over a broad range of Fourier frequencies is caused by propagating mass-accretion fluctuations generated in the blackbody-emitting disc. Our results thus support with a quantitative model the scenario proposed by Wilkinson & Uttley (2009) and Uttley et al. (2011). The speed of the propagating fluctuations is much faster than expected for a classical thin ($h \ll r$) disc however. For disc accretion, the scaling of radial infall time in our model parametrization gives $s_t f_{\Delta t,d} r^{n_{t,d}} \simeq \alpha^{-1}(r/h)^2$, so that for our assumed {$s_t, f_{\Delta t,d}$} = 1, $n_{t,d}$ = 2 and $\alpha$ = 0.1, we infer $h \simeq 3$ across the range of variable disc radii (e.g. 46 $\gtrsim$ r $\gtrsim$ 10).

It is important to consider whether such a fast accretion disc remains significantly optically thick to Compton scattering in order to produce the observed blackbody disc spectrum. To estimate the vertically integrated optical depth of the disc ($\tau_z$), we assume constant mass accretion rate $\dot{M}$ and use the continuity equation for $\dot{M}$ in terms of disc surface density $\Sigma$, radius $R = r R_g$, and radial velocity $v_r$ (e.g. see Frank et al. 2002), i.e. $\dot{M} = 2\pi R \Sigma v_r$. From equation (25), we obtain $v_r = (s_t f_{\Delta t,d})^{-1} r_i^{-(n_{t,d}+\frac{1}{2})} c$. The vertically integrated optical depth for electron scattering of photons is $\tau_z = \Sigma \kappa_{es}$ where $\kappa_{es}$ is the electron-scattering opacity of the gas. We can therefore rewrite the mass accretion rate in terms of $\tau_z, r, black hole mass M, and various constants as

$$ \dot{M} = \left( \frac{2\pi G M}{\kappa_{es} c} \right) \tau_z (s_t f_{\Delta t,d})^{-1} r^{(\frac{1}{2}-n_{t,d})}. \quad (32) $$

We can also write $\dot{M}$ in terms of the radiated bolometric luminosity $L$ expressed as a fraction of the Eddington luminosity $L_{Edd} = 4\pi G Mc/\kappa_{es}$, and the radiative efficiency of the accretion flow $\eta$:

$$ \dot{M} = \left( \frac{4\pi G M}{\kappa_{es} c} \right) \eta^{-1} \frac{L}{L_{Edd}}. \quad (33) $$

Substituting and rearranging the above equations yields the following expression for $\tau_z$:

$$ \tau_z = 2\eta^{-1} \left( \frac{L}{L_{Edd}} \right) s_t f_{\Delta t,d} r^{(n_{t,d}-\frac{1}{2})}. \quad (34) $$

If we assume $\eta = 0.25$ (consistent with $r_{in} = 2 R_g$) then for {$s_t, f_{\Delta t,d}$} = 1, $n_{t,d}$ = 2, and $L/L_{Edd} = 0.1$, we obtain $\tau_z \simeq 25$ at $r = 10 R_g$, with $\tau_z \propto r^{3/2}$. Note that this calculation is independent of the vertical scale-height of the disc and confirms robustly that radial infall times that are fast enough to explain the observed lags, also remain consistent with the observed blackbody emission. Therefore, the accretion disc required to explain propagation delays is both fast and optically thick.

The inferred radial scaling of disc time-scales $n_{t,d} = 2$, means that the disc variability time-scale is proportional to the Keplerian time-scale multiplied by $r^2$, which implies constant $h$ versus $r$ if the variability time-scale behaves as a standard accretion disc viscous time-scale. Such discs have not been studied using GRMHD simulations, which tend to assume constant $h/r$ by construction, since they focus on the hot inner accretion flows commonly associated with the hard state or with radiatively inefficient accretion. Besides being strongly inconsistent with the energy-dependent PL-D hard lags

<footnote n="5">⁵Note that in their Comptonization model for quasi-periodic oscillation spectral-timing properties, Karpouzas et al. (2021) and Bellavita et al. (2022) do include separate time-variable coronal heating and seed photon oscillations, but these are set to be in phase with one another, so that Comptonization light-travel delays dominate the lags.</footnote>
(see Appendix C), constant $h/r$ ($n_{t,d} = 0$) in the blackbody-emitting optically thick disc would also imply that disc vertical optical depth decreases with radius (see equation 34), which is inconsistent with the observed disc blackbody emission. If constant $h$ versus $r$ is a better fit to the data, it would be useful to investigate whether such discs can be produced in GRMHD simulations of inner accretion flows. It is worth noting that constant $h$ corresponds to the radiation-pressure dominated inner disc of Shakura & Sunyaev (1973), although such discs are generally associated with higher luminosities than observed in the hard state. Radiation-pressure dominated discs are thermally and viscously unstable (Lightman & Eardley 1974; Shakura & Sunyaev 1976), but a strong vertical magnetic field could stabilize the disc (Mishra et al. 2022), and might be expected in the hard state.

In our model, we have not explicitly considered the possibility that mass accretion fluctuations originate at discrete radii in the accretion flow, which has been suggested as a way to explain the step-like features observed in some PL-PL lag-frequency spectra and the distinct broad bumps seen in the power spectra (Arévalo & Uttley 2006; Mahmoud & Done 2018a, b). Nor have we considered a discontinuity in the fluctuation time-scale, e.g. at the disc truncation radius, as an origin for these features (Rapisarda et al. 2016). In fact, the smoothness of the PL-D lag-frequency dependence seems to disfavour models where fluctuations originate from just a few distinct radii or where there is a sharp change in variability time-scale. For example, the predicted bumps in the PL-D lag-frequency spectrum, and their inconsistency with the data, can be seen in examples of such models shown in Kawamura et al. (2022). One possibility is that PL-PL lags on long time-scales are associated with a different phenomenon than propagation, e.g. short-term variability of the coronal geometry (e.g. see Bollemeijer et al. 2024). Thus, on certain time-scales, the lags may be enhanced and we see the effect as steps in lag versus frequency.

## 6.3 Evolution of coronal geometry

In our model, the power-law versus disc (PL-D) soft lags and the PL-PL lags between power-law bands encode information about the temporal response of seed photons to mass accretion fluctuations, which is strongly dependent on coronal geometry. Observationally, these lags are seen to evolve through the hard state in ways which are strongly correlated with the X-ray spectral shape (Altamirano & Méndez 2015; Reig et al. 2018) and timing properties (Pottschmidt et al. 2003; Wang et al. 2022). We can use the lag model presented here, in particular the results in Fig. 8, to interpret these changes in terms of the changing coronal geometry, subject to the caveat that the model presented so far remains incomplete (see Section 6.4). We also implicitly assume that spectral-timing changes are linked to coronal geometry rather than changes in accretion flow variability properties such as time-scales or maximum fluctuation radius. As shown in Section 5.4, changes in variability time-scales can mimic the effects of changes in $r_{cor}$, but our model does not predict the evolution of $\Gamma$ in this case. Changes in maximum fluctuation radius may impact the maximum PL-PL hard lag, but should not significantly affect the high-frequency spectral-timing behaviour such as PL-D soft lags and crossing-frequency.

First, it is notable that changes in coronal radius result in PL-PL and PL-D lag changes that are anticorrelated with the power-law photon index. This result seems to be contrary to the observed trend of increasing lags with increasing photon index (Pottschmidt et al. 2003; Reig et al. 2018). Such evolution is better explained with a simultaneous increase in both coronal height and opening angle, since neither coronal parameter on its own can cause a strong increase

in both lag and photon index at the same time (see the top-left panel of Fig. 8). One should be cautious, however, since the observed lags are usually calculated from the cross-spectrum averaged over a broad frequency range for the purposes of correlating them with changes in photon index. Therefore, changes in the lag-frequency shape may also be important, but these effects are not captured in the simple analysis of the maximum lag amplitude used to make Fig. 8. For example, above 10 Hz, the PL-PL hard lags appear to be strongly anticorrelated with coronal radius for the inverted cone geometry (see Fig. 7). Changes in coronal radius would lead to a strong positive correlation of PL-PL hard lags with photon index at those frequencies.

Wang et al. (2022) showed for a sample of transient BHXRBs, the strong connection of PL-D soft lags with the broadband power-spectral shape parametrized by the power-spectral hue (Heil et al. 2015). Our model does not predict the power-spectral properties, which depend on the assumed radial dependence of variability time-scales and the maximum radius of fluctuations (which sets the break-frequency of the power spectrum). Without speculating about the physical connection between timing and coronal geometry, we note that the timing properties also link to the spectral shape of coronal emission, with photon index increasing from the hard state to the SIMS (Reig et al. 2018).

In the analysis of Wang et al. (2022), PL-D soft lags are seen to first increase in amplitude from < 1 ms to 2 ms during the hardest part of the hard state, then decrease back to < 1 ms before increasing sharply to > 10 ms as the source evolves from the HIMS to SIMS, where broadband-noise variability is quenched (Belloni et al. 2005; Heil et al. 2015). Studying the top-right panel of Fig. 8, we see that the decrease in PL-D soft lag during the hard state could possibly be linked to a small drop in disc truncation radius if it is associated with an increase in photon index, while the following large increase in lag during the HIMS may be explained with a combination of increasing coronal height and opening angle, which could also explain the PL-PL hard lag evolution.

Thus, the overall picture proposed by Wang et al. (2021, 2022), of increasing coronal height during the HIMS, remains applicable when considering our model. However the required increase in height is much smaller than that inferred from light-travel time lags. PL-D soft lags exceeding 10 ms can be produced by inverted cone coronae with height ~ 10 $R_g$, which appear to be consistent with the relatively small coronal heights inferred from modelling of reflection features. Such coronae would also not appear as being jet-like, so they are likely to be consistent with observed X-ray polarization signatures which imply a Comptonizing region which is extended orthogonally to the radio jet axis (Krawczynski et al. 2022; Veledina et al. 2023; Ingram et al. 2024).

Our model provides a natural explanation for the behaviour of BHXRBs in the HIMS which show only large PL-D soft lags across the full observed frequency range (Wang et al. 2022). In Fig. 8, the lower-left panel shows that the PL-D soft lag and lag crossing-frequency (from hard to soft lags) are tightly correlated. Therefore, if the crossing-frequency is smaller than the minimum fluctuation frequency produced at the largest variable radius in the disc, i.e. the power-spectral break frequency, hard lags are never seen. For soft lags > 10 ms this will occur for power spectral breaks ~ 1 Hz or smaller. We demonstrate this effect in Fig. 10. An interesting observational example is shown by Wang et al. (2022) for GX 339-4 (see their fig. 3), which shows in its intermediate state large (> 30 ms) soft lags extending well below the break frequency at ~ 0.5 Hz. Kawamura, Done & Takahashi (2023) have interpreted the large negative lags extending to low frequencies as a consequence
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Comparison of spectral-timing properties for inverted cone coronae with radius $r_{cor} = 10 R_g$ and maximum variable disc radii of $46 R_g$ and $23 R_g$ producing fluctuation frequencies $v_{fluc,max}$ of 0.03 and 0.3 Hz, respectively. For coronal height $h_{cor} = 20 R_g$ the positive PL-D hard lags seen for $v_{fluc,max} = 0.03$ Hz disappear when the maximum variable radius decreases, with only PL-D soft lags seen for $v_{fluc,max} = 0.3$ Hz. PL-D hard lags return if the corona then shrinks to $h_{cor} = 10 R_g$.</figcaption></figure>

of seed photons with energies significantly lower than those from the inner disc, e.g. due to cyclotron emission in the hot coronal flow. However, in our model interpretation such seed photons are not required. Instead, this situation would correspond to the largest variable disc radius lying well inside the seed-producing region of the disc. This scenario might occur if the variable region of the disc shrinks, or the corona increases in height, or through some combination of these changes.

## 6.4 Extensions and improvements to the model

An important caveat to our model is that it treats the seed photons as if they are mono-energetic, and furthermore does not include the inevitable spectral overlap of power-law and disc emission components which should be a consequence of the disc providing seed photons to the corona. Both these aspects can be dealt with using a full energy-dependent spectral-timing treatment to predict lag versus energy spectra among other quantities. Such a treatment would also account for the energy-dependence of the disc lags, which we touch on briefly in Appendix C. The radial temperature-dependence of the disc leads to a strong energy dependence of the PL-D lags which can be used to constrain the radial scaling of disc variability time-scales. For our current model it is encouraging to note that for the time-scale radial scaling index $n_{t,d} = 2$ assumed here, the lags of disc bolometric luminosity are similar to the lags at $E = 2.6k_B T_{max}$ which is close to the disc blackbody flux peak (see the lower panel of Fig. C1). Thus the PL-D lag results inferred in this work for total disc luminosity are unlikely to deviate much from those obtained at the disc peak, when a full spectral treatment is used.

Our treatment of coronal geometry is currently very simple, since we consider the corona to be a solid body for the purposes of calculating the solid angle of the corona as seen from the disc and the radial dependence of coronal flux illuminating the disc. In reality, the

corona is unlikely to have a sharp 'edge' and one should also include optical depth effects on the fraction of seed photons intercepted by the corona, as well as obscuration of part of the direct inner disc emission from the observer. For calculating coronal illumination of the disc as well as observed power-law emission from a given direction, one should account for angular dependence of coronal emission, which will be especially relevant for non-spherical coronae such as the inverted cone considered here. A detailed treatment would require input from fully-relativistic Comptonization models, but it is unlikely that the basic results shown here would change significantly, only the more detailed dependence of spectral-timing properties on coronal parameters. This is because the key requirement for PL-PL hard lags and PL-D soft lags to be produced by propagation is that power-law spectral pivoting takes place in response to seed and coronal heating variations. Provided the corona is in thermal equilibrium in response to variations, spectral pivoting will occur due to photon number and luminosity conservation, independent of the details of Comptonization models.

So far we have only considered seed photons which originate from the disc, but it is also possible that seed photons might be internally generated in the corona, e.g. via internal synchrotron emission if the corona is magnetized (Veledina et al. 2013). It is easy to include such emission in our model by converting some fraction of accretion power dissipated into the corona into seed luminosity. This would have the effect of pushing the seed impulse response centroid closer to the heating centroid value, which would reduce the PL-PL hard lags and PL-D soft lags. Therefore, in the context of our model, the existence of significant internal seed photons is constrained by the observed large hard and soft lags. However, a quantitative treatment would require the spectral shape of internal seed photons to be accounted for. An interesting comparison comes from accreting neutron stars, which, assuming they possess a corona similar to black holes, are expected to produce a significant source of X-ray seed photons from the neutron star boundary layer with the accretion flow. This means that neutron stars in equivalent hard states should not show large PL-PL hard lags or PL-D soft lags, as a result of the central source of seed photons which pushes the seed impulse centroid to zero delay (Basak & Uttley in preparation).

Finally, we note that a complete model for lags should account for the effects of strong gravity and relativistic velocities close to the black hole. Relativistic beaming due to orbital motion will significantly affect the direct disc emission seen by the observer (especially when considering energy-dependent disc lags), but not the seed impulse response seen from the central corona. However, for small disc radii, light-bending could play a significant role in enhancing the seed flux from the inner disc into the corona and at the same time amplify the coronal illumination of the inner disc while reducing illumination of larger disc radii, which would also increase the seed flux due to reverberation. Both effects would shift the seed centroid to smaller delays and thus reduce the amplitudes of both the PL-PL hard lags and the PL-D soft lags, while also reducing the PL-D hard lags due to the enhanced disc reverberation component. A natural step for the model would be to incorporate the relativistic effects on seed photons and direct disc emission with the expected reverberation light-travel delays from a fully ray-traced system, e.g. following a similar approach to the RELTRANS model of Ingram et al. (2019) and Mastroserio et al. (2021), albeit with an extended coronal geometry rather than a single lamppost (e.g. see Lucchini et al. 2023). For full consistency, it would be necessary to calculate reverberation delays, coronal disc illumination, and seed photon illumination of the corona using the same coronal geometry throughout.
# 7 CONCLUSIONS

We have shown that the complex and evolving pattern of Fourier-frequency and energy-dependent X-ray time lags seen in hard and HIMS BHXRBs can be explained by propagation of mass accretion fluctuations through the disc to the inner corona. Specifically,

(i) Lower-frequency fluctuations are generated at larger disc radii, propagating through the disc to the corona to produce the observed large low-frequency delays of coronal power-law emission with respect to disc blackbody emission, which decrease with increasing frequency.

(ii) The disc emission provides coronal seed photons, which vary in response to propagating fluctuations before the corona is heated by the same fluctuations. This leads to softening and then hardening of the power-law spectrum, producing the observed hard lags between power-law bands which increase log-linearly with energy separation and are smaller at high Fourier-frequencies.

(iii) Some fraction of coronal luminosity intercepts and heats the disc to produce blackbody reverberation emission. Since variations in internal coronal heating dominate the variability of bolometric coronal luminosity while the low-energy coronal power-law variability is dominated by seed photon variations, the disc reverberation signal shows large delays compared to the soft power-law emission, which depend on the propagation time from the seed-emitting part of the disc to the corona. At high frequencies, these delays dominate the power-law versus disc lags to produce the observed 'soft' lags.

(iv) The smaller fraction of disc emission originating from small radii suppresses the disc and seed photon variability at high frequencies, which leads naturally to energy-dependent suppression of high-frequency variability, with the disc and softer power-law energies showing less high-frequency variability, consistent with observations.

(v) Because of the role of seed photons in these effects, the spectral-timing properties strongly depend on coronal geometry, with coronae that are more extended (radially or horizontally) showing larger power-law versus power-law hard lags and power-law versus disc soft lags. However, even soft lags > 10 ms can be explained with relatively compact coronae (heights and radii ~ 10 $R_g$) because they are associated with propagation delays, not light-travel times.

(vi) The observed relatively weak dependence of power-law versus disc hard lags on disc photon energy, with lags increasing down to the power-spectral low-frequency break, imply that disc variability is generated over a relatively small range of radii (a few tens of $R_g$). This could be explained if the disc scale height *h* is constant with radius.

The relatively compact coronae which we infer from the lags could help to reconcile the measurements of large hard and soft lags with results obtained from spectral fitting, which imply relatively small disc truncation radii in the hard state, and X-ray polarimetry, which further implies that the X-ray emitting region is not significantly vertically extended. Because the model links both the power law versus power law and power law versus disc lags with the coronal geometry, it can be used to provide a powerful independent check on spectral and polarimetric fitting results. However, a more exciting prospect would be to simultaneously fit the spectral-timing properties, X-ray spectrum and polarimetry from a BHXRB to obtain the strongest possible constraints on the disc-corona geometry and accretion flow properties. Additional strong constraints could be provided by incorporating the expected reverberation light-travel lags which should appear at frequencies ~ 100 Hz where propagation effects are washed out, and which will be revealed by proposed large-

area X-ray spectral-timing missions such as eXTP (Zhang et al. 2019) and STROBE-X (Ray et al. 2018).

# ACKNOWLEDGEMENTS

We would like to thank Niek Bollemeijer for providing the MAXI J1820 + 070 spectral-timing and spectral data shown in this paper, and the referee Chris Done for many valuable comments and suggestions. This work was supported by the Action Thématique Hautes Energies of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES.

# DATA AVAILABILITY

The code and Jupyter notebooks used to calculate the model spectral-timing properties and produce the plots shown in this paper are available at the following DOI: https://doi.org/10.5281/zenodo.14038609. This research has made use of NICER data obtained through the High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center. The data used in this article are available in HEASARC at https://heasarc.gsfc.nasa.gov/docs/archive.html.

# REFERENCES

* Altamirano D., Méndez M., 2015, MNRAS, 449, 4027

* Arévalo P., Uttley P., 2006, MNRAS, 367, 801

* Atri P. et al., 2020, MNRAS, 493, L81

* Bellavita C., García F., Méndez M., Karpouzas K., 2022, MNRAS, 515, 2099

* Belloni T., Homan J., Casella P., van der Klis M., Nespoli E., Lewin W. H. G., Miller J. M., Méndez M., 2005, A&A, 440, 207

* Beloborodov A. M., 2001, Adv. Space Res., 28, 411

* Bollemeijer N. et al., 2024, MNRAS, 528, 55

* Bollimpalli D. A., Mahmoud R., Done C., Fragile P. C., Kluźniak W., Narayan R., White C. J., 2020, MNRAS, 496, 3808

* Buisson D. J. K. et al., 2019, MNRAS, 490, 1350

* Churazov E., Gilfanov M., Revnivtsev M., 2001, MNRAS, 321, 759

* Cui W., Zhang S. N., Focke W., Swank J. H., 1997, ApJ, 484, 383

* Dauser T., Garcia J., Parker M. L., Fabian A. C., Wilms J., 2014, MNRAS, 444, L100

* De Marco B. et al., 2017, MNRAS, 471, 1475

* De Marco B., Ponti G., Muñoz-Darias T., Nandra K., 2015a, MNRAS, 454, 2360

* De Marco B., Ponti G., Muñoz-Darias T., Nandra K., 2015b, ApJ, 814, 50

* De Marco B., Zdziarski A. A., Ponti G., Migliori G., Belloni T. M., Segovia Otero A., Dzielak M. A., Lai E. V., 2021, A&A, 654, A14

* De Rosa A. et al., 2019, Sci. China Phys. Mech. Astron., 62, 29504

* Done C., Gierliński M., Kubota A., 2007, A&AR, 15, 1

* Frank J., King A., Raine D. J., 2002, Accretion Power in Astrophysics, 3rd edn. Cambridge Univ. Press, Cambridge

* García J. A., Fabian A. C., Kallman T. R., Dauser T., Parker M. L., McClintock J. E., Steiner J. F., Wilms J., 2016, MNRAS, 462, 751

* García J. et al., 2014, ApJ, 782, 76

* Gardenier D. W., Uttley P., 2018, MNRAS, 481, 3761

* Gendreau K. C. et al., 2016, in den Herder J.-W. A., Takahashi T., Bautz M., eds, Proc. SPIE Conf. Ser. Vol. 9905, Space Telescopes and Instrumentation 2016: Ultraviolet to Gamma Ray. SPIE, Bellingham, p. 99051H

* Giannios D., Kylafis N. D., Psaltis D., 2004, A&A, 425, 163

* Grinberg V. et al., 2014, A&A, 565, A1

* Heil L. M., Uttley P., Klein-Wolt M., 2015, MNRAS, 448, 3339

* Heil L. M., Vaughan S., Uttley P., 2012, MNRAS, 422, 2620

* Hogg J. D., Reynolds C. S., 2016, ApJ, 826, 40

* Hua X.-M., Kazanas D., Cui W., 1999, ApJ, 512, 793

* Ingram A. R., Motta S. E., 2019, New A Rev., 85, 101524
Ingram A., Bollemeijer N., Veledina A., Dovčiak M., Poutanen J., Egron E., Russell T. D., et al., 2024, ApJ, 968, 76

Ingram A., Done C., 2011, MNRAS, 415, 2323

Ingram A., Mastroserio G., Dauser T., Hovenkamp P., van der Klis M., García J. A., 2019, MNRAS, 488, 324

Kalemci E., Kara E., Tomsick J. A., 2022, in, Handbook of X-ray and Gamma-ray Astrophysics, Springer, Singapore, 9

Kara E. et al., 2019, Nature, 565, 198

Karpouzas K., Méndez M., García F., Zhang L., Altamirano D., Belloni T., Zhang Y., 2021, MNRAS, 503, 5522

Kawamura T., Axelsson M., Done C., Takahashi T., 2022, MNRAS, 511, 536

Kawamura T., Done C., Takahashi T., 2023, MNRAS, 525, 1280

Kazanas D., Hua X.-M., Titarchuk L., 1997, ApJ, 480, 735

Kotov O., Churazov E., Gilfanov M., 2001, MNRAS, 327, 799

Krawczynski H. et al., 2022, Science, 378, 650

Kylafis N. D., Papadakis I. E., Reig P., Giannios D., Pooley G. G., 2008, A&A, 489, 481

Lightman A. P., Eardley D. M., 1974, ApJ, 187, L1

Lucchini M. et al., 2023, ApJ, 951, 19

Lyubarskii Y. E., 1997, MNRAS, 292, 679

Mahmoud R. D., Done C., 2018a, MNRAS, 473, 2084

Mahmoud R. D., Done C., 2018b, MNRAS, 480, 4040

Mahmoud R. D., Done C., De Marco B., 2019, MNRAS, 486, 2137

Makishima K. et al., 2008, PASJ, 60, 585

Makishima K., Maejima Y., Mitsuda K., Bradt H. V., Remillard R. A., Tuohy I. R., Hoshi R., Nakagawa M., 1986, ApJ, 308, 635

Mastroserio G. et al., 2021, MNRAS, 507, 55

Mastroserio G., Ingram A., van der Klis M., 2018, MNRAS, 475, 4027

Mishra B., Fragile P. C., Anderson J., Blankenship A., Li H., Nalewajko K., 2022, ApJ, 939, 31

Misra R. et al., 2017, ApJ, 835, 195

Mitsuda K. et al., 1984, PASJ, 36, 741

Miyamoto S., Kitamoto S., Mitsuda K., Dotani T., 1988, Nature, 336, 450

Mushtukov A. A., Ingram A., van der Klis M., 2018, MNRAS, 474, 2259

Nowak M. A., 2000, MNRAS, 318, 361

Nowak M. A., Vaughan B. A., Wilms J., Dove J. B., Begelman M. C., 1999, ApJ, 510, 874

Pietrini P., Krolik J. H., 1995, ApJ, 447, 526

Pottschmidt K. et al., 2003, A&A, 407, 1039

Rapisarda S., Ingram A., Kalamkar M., van der Klis M., 2016, MNRAS, 462, 4078

Ray P. S. et al., 2018, in den Herder J.-W. A., Nikzad S., Nakazawa K., eds, Proc. SPIE Conf. Ser. Vol. 10699, Space Telescopes and Instrumentation 2018: Ultraviolet to Gamma Ray. SPIE, Bellingham, p. 1069919

Reig P., Kylafis N. D., 2021, A&A, 646, A112

Reig P., Kylafis N. D., Giannios D., 2003, A&A, 403, L15

Reig P., Kylafis N. D., Papadakis I. E., Costado M. T., 2018, MNRAS, 473, 4644

Ross R. R., Fabian A. C., 2007, MNRAS, 381, 1697

Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337

Shakura N. I., Sunyaev R. A., 1976, MNRAS, 175, 613

Skipper C. J., McHardy I. M., 2016, MNRAS, 458, 1696

Skipper C. J., McHardy I. M., Maccarone T. J., 2013, MNRAS, 434, 574

Sunyaev R., Revnivtsev M., 2000, A&A, 358, 617

Timmer J., Koenig M., 1995, A&A, 300, 707

Torres M. A. P., Casares J., Jiménez-Ibarra F., Álvarez-Hernández A., Muñoz-Darias T., Armas Padilla M., Jonker P. G., Heida M., 2020, ApJ, 893, L37

Uttley P., Cackett E. M., Fabian A. C., Kara E., Wilkins D. R., 2014, A&AR, 22, 72

Uttley P., McHardy I. M., Vaughan S., 2005, MNRAS, 359, 345

Uttley P., Wilkinson T., Cassatella P., Wilms J., Pottschmidt K., Hanke M., Böck M., 2011, MNRAS, 414, L60

van Paradijs J., Lewin H. G., 1986, A&A, 157, L10

Veledina A., et al., 2023, ApJL, 958, L16

Veledina A., Poutanen J., Vurm I., 2013, MNRAS, 430, 3196

Wang J. et al., 2021, ApJ, 910, L3

Wang J. et al., 2022, ApJ, 930, 18

Wilkinson T., Uttley P., 2009, MNRAS, 397, 666

Wilms J., Allen A., McCray R., 2000, ApJ, 542, 914

Zdziarski A. A., Dzielak M. A., De Marco B., Szanecki M., Niedźwiecki A., 2021, ApJ, 909, L9

Zdziarski A. A., Johnson W. N., Magdziarz P., 1996, MNRAS, 283, 193

Zdziarski A. A., Poutanen J., Mikolajewska J., Gierlinski M., Ebisawa K., Johnson W. N., 1998, MNRAS, 301, 435

Zhang S. et al., 2019, Sci. China Phys. Mech. Astron., 62, 29502

Zycki P. T., Done C., Smith D. A., 1999, MNRAS, 309, 561

## APPENDIX A: CORONAL GEOMETRY CALCULATIONS

For simplicity, we calculate the fractions of photons intercepted by the corona from the disc, and vice versa, by assuming that the corona interacts with light as a solid body. Neglecting gravitational lightbending and Doppler beaming, the fraction of disc photons from a radius $r$ that are intercepted by the corona, $f_{d→c}(r)$, depends on the solid angle subtended by the corona as seen from $r$, as well as a weighting by the cosine of the angle between the direction to the corona and the normal to the disc plane (Lambert's cosine law):

$$ 
\begin{aligned}
f_{d \to c}(r) &= \int_{\Omega_{\text{cor}}} \hat{\mathbf{n}}_{\text{d}} \cdot \mathbf{d}\Omega_{\text{cor}} / \int_{2\pi} \hat{\mathbf{n}}_{\text{d}} \cdot \mathbf{d}\Omega_{\text{sph}} \\
&= \frac{1}{\pi} \int_{\Omega_{\text{cor}}} \hat{\mathbf{n}}_{\text{d}} \cdot \mathbf{d}\Omega_{\text{cor}}.
\end{aligned} 
\quad \text{(A1)}
$$

Here, $\Omega_{\text{cor}}$ and $\Omega_{\text{sph}}$ are the solid angles of the visible surface of the corona seen from the location on the disc at $r$ and a sphere of visibility centred on the same location, $\hat{\mathbf{n}}_{\text{d}}$ is the unit vector normal to the disc plane, and $\mathbf{d}\Omega = \hat{\mathbf{p}} \cdot \mathbf{d}\Omega$ where $\mathbf{d}\Omega$ is the solid angle of a surface element of the sphere of visibility for a visible (i.e. not self-obscured) surface element of the corona as seen from the disc location at $r$ and $\hat{\mathbf{p}}$ is the unit vector pointing towards that disc location from the surface element. By assuming cylindrical symmetry, the calculated fraction applies to emission from the entire disc radius $r$ and not just the single disc location calculated for.

We illustrate the situation for calculating solid angles of surface elements in Fig. A1. The solid angle of a coronal surface element with area $dA_e$ as seen from the disc at a given location, is given by $d\Omega = \hat{n}_e \cdot \hat{p} \cdot dA_e / p^2$, where $\hat{n}_e$ is the unit vector normal to the surface element and $p$ is the magnitude of the position vector $p$ from the disc location to the surface element, i.e. $p = p \cdot \hat{p}$. Although some special cases can be solved analytically to obtain $f_{d→c}(r)$ (such as the spherical corona, see below), for most cases we must solve numerically by summing over all the coronal surface elements (denoted by the index $i$) which are visible from a location on the disc at a given radius:

$$ f_{d \to c}(r) = \frac{1}{\pi} \sum_{\hat{n}_{e,i} \cdot \hat{p}_i > 0} \hat{n}_{d,i} \cdot \hat{p}_i \hat{n}_{e,i} \cdot \hat{p}_i \frac{dA_{e,i}}{p_i^2}. \quad \text{(A2)} $$

The requirement $\hat{n}_{e,i} \cdot \hat{p}_i > 0$ ensures that the element's outward face is visible from the disc location, provided the geometry is such that no part of the corona blocks the view of the outward face. The dependence on $r$ is not explicitly shown but arises in $p_i$ and $\hat{p}_i$.

To determine the fraction of coronal emission intercepted by a unit surface area of the disc at radius $r$, $f_{d→c}(r)$, we assume that the corona has a uniform intrinsic surface brightness, with luminosity per unit surface area given by the total coronal luminosity divided by the total emitting surface area of the corona $A_{\text{cor}}$. With this simplifying assumption, the fraction of total coronal luminosity, $d f_e$ reaching a disc element of area $dA_d$ can be calculated in a similar way to the contribution of a coronal element to the fraction of disc luminosity
<figure><img src="image_13.png" /><figcaption>**Figure A1.** Setup for calculating the solid angle of a coronal surface element of area $dA_\epsilon$ as seen from a disc surface element of area $dA_d$, and vice versa. See Appendix A for details of the calculation.</figcaption></figure>

intercepted by the corona (see Fig. A1):

$$d f_{\epsilon} = \frac{1}{\pi} \hat{n}_{\epsilon} \cdot \hat{p} \hat{n}_{d} \cdot \hat{p} \frac{dA_{d}}{p^{2}} \frac{dA_{\epsilon}}{A_{\text{cor}}}. \qquad (A3)$$

Summing over all visible coronal elements, and multiplying by $1/dA_d$ to obtain the fraction of coronal emission intercepted by unit surface area of the disc, we obtain

$$f_{c \to d}(r) = \frac{1}{\pi} \sum_{\hat{n}_{\epsilon,i} \cdot \hat{p}_i > 0} \hat{n}_{\epsilon,i} \cdot \hat{p}_i \hat{n}_{d,i} \cdot \hat{p}_i \frac{dA_{\epsilon,i}}{A_{\text{cor}} p_i^2}. \qquad (A4)$$

Comparing this equation with equation (A2), we find the simple result:

$$f_{c \to d}(r) = \frac{f_{d \to c}(r)}{A_{\text{cor}}}. \qquad (A5)$$

For a spherical corona of radius $r_{cor}$, we can solve equation (A1) analytically to obtain:

$$f_{d \to c}(r) = \frac{1}{\pi} \left[ \arcsin\left(\frac{r_{\text{cor}}}{r}\right) - \sqrt{\left(\frac{r_{\text{cor}}}{r}\right)^2 - \left(\frac{r_{\text{cor}}}{r}\right)^4} \right] \qquad (A6)$$

with $f_{c \to d}(r)$ obtained from equation (A5) assuming $A_{cor} = 2\pi r_{cor}^2$ (since we only consider the coronal hemisphere on the observable side of the disc).

For the inverted cone corona modelled in this work, we solve equation (A2) numerically by splitting the corona into surface elements defined by a grid of 1000 × 1000 bins in azimuth (from 0 to $2\pi$) and coronal height (from 0 to $h_{cor}$). To calculate $f_{c \to d}(r)$ we use $A_{cor} = \pi (r_{top}^2 + (r_{cor} + r_{top})\sqrt{(r_{top} - r_{cor})^2 + h_{cor}^2})$, where $r_{top}$ is the radius of the top surface of the corona ($r_{top} = r_{cor} + h \tan \theta_{cor}$).

Fig. 6 shows $f_{d \to c}(r)$ calculated for the different geometries and geometric parameters, to demonstrate the effects of changing coronal radius, height and opening angle.

## APPENDIX B: TEST OF LINEARITY WITH TIME-SERIES SIMULATIONS

We can test the validity of our use of linear impulse response functions to calculate the spectral-timing properties of our model using stochastic time-series simulations of a simplified variable accretion flow. Here, we follow the approach of Arévalo and Uttley (2006) and use the method of Timmer and Koenig (1995) to simulate a time series of fractional accretion rate modulation, $\dot{m}_i(t)$ (normalized to a mean $\langle \dot{m}_i(t) \rangle = 1$) for each annulus of the disc corresponding to radius $r_i$, using the Lorentzian parameters associated with each radius to determine the $\dot{m}_i(t)$ power spectrum. The physical accretion rate time-series $\dot{M}_i(t)$ for a given radius is produced by multiplying

<figure><img src="image_14.png" /><figcaption>**Figure B1.** Comparison of frequency-dependent spectral-timing properties for numerical simulations of our lags model using a simplified accretion flow with 40 per cent total rms variability of accretion rate, versus the results obtained from linear impulse response calculations, assuming a spherical corona with $r_{cor} = 10 R_g$. Top panel: Time lags for power law at $E_{soft}$ versus disc (solid line), power law at $E_{med}$ versus $E_{soft}$ (dashed), power-law at $E_{hard}$ versus $E_{med}$ (dotted). Positive lags indicate that the harder band lags the softer band (or the power-law lags the disc). Lower panel: power spectra for disc (solid line) and power-law flux at $E_{med}$ (dashed) and $E_{hard}$ (dotted). The signal peak frequency at $r_{cor}$ is shown by the vertical dot-dashed grey line.</figcaption></figure>

together the $\dot{m}(t)$ generated at that radius with the physical time-series generated for the preceding (outer) radius and shifting the resulting signal by the propagation delay across the radial bin, $\Delta\tau_i$:

$$\dot{M}_i(t + \Delta\tau_i) = \dot{m}_i(t) \dot{M}_{i+1}(t) \qquad (B1)$$

By repeating this calculation for successive inwards radii, we can generate an accretion fluctuation time series for each radius, from which we can determine the dissipation power time series for each radius, by multiplying by $f_{diss}(r_i)$, calculated from equation (7). This can in turn be used to generate the disc dissipation and coronal heating time series associated with each radius, the seed contribution from disc dissipation (using $f_{d \to c}(r_i)$ and from summing their totals over all radii we can determine the total heating, disc, and seed light curves, including reverberation contributions to the disc and seed contributions (by applying $f_{c \to d}(r_i)$). Finally, we can use the generated time-series directly in equation (1) to determine a time-series for photon index $\Gamma$, and use equation (2) to determine the power-law flux at a given energy.

This approach provides a numerical non-linear time-series method for calculating spectral-timing properties using our lags model. The method is non-linear for two reasons. Firstly the generated accretion rate time-series are formally non-linear, incorporating the rms-flux relation and log-normal flux distributions that are observed in the X-ray variability of accreting black holes (Uttley et al. 2005). Secondly, the power-law flux variations are generated from the full, non-linear versions of equations (1) and (2) rather than the linearized versions that are used to determine the impulse responses of the power-law flux.

Fig. B1 shows the spectral-timing results of our numerical simulation assuming a spherical coronal geometry ($r_{cor} = 10 r_g$) and with
disc and coronal variability and propagation time-scales as described in Section 4, except that we assumed a fluctuation frequency at the maximum variable radius $v_{fluc,max} = 0.1$ Hz, to reduce the required light curve lengths and thus save computation time. The rms amplitude of accretion variability is set to be 40 per cent, consistent with variations often observed in the hard state. For the simulation, 350 sets of light curves, each consisting of $2^{22}$ 14.1 µs time bins ($\approx 59.3$ s per light curve) were simulated and used to calculate the spectral timing products, which are further binned geometrically in frequency, with geometric binning factor 1.05. For comparison we show the corresponding results for the impulse response approach that is used in this work. The only significant deviations that can be seen at low frequencies are due to the intrinsic noise variability of the simulated time-series, which lead to fluctuations in the time-lags and the signal power (and hence the power spectra of emission). At higher frequencies this noise in the simulated data is smoothed over, but small shifts in the power spectra can be seen with power enhanced by up to ~20% at high frequencies in the simulated data. Only very small deviations are seen in the lags however. The increase in the high-frequency power from the simulated data is a known effect of the multiplicative variability process used to generate the light curves (Uttley et al. 2005). However, like the lags, power-spectral ratios between different bands are not changed significantly. We conclude that even for relatively large amplitudes of variability, our linear impulse response model is an excellent approximation to the full non-linear model.

## APPENDIX C: DISC ENERGY-DEPENDENT LAGS AS A CONSTRAINT ON THE RADIAL SCALING OF DISC VARIABILITY TIME-SCALES

In this work we calculate spectral-timing properties using the bolometric disc luminosity, mainly for simplicity and generalizability. A complete energy-dependent treatment of the lags involving the disc emission also requires consideration of the overlap of the disc and coronal emission in energy, including a more realistic seed photon spectrum. We leave this detailed treatment to a future work. However it is useful to note an important constraint that the energy-dependence of the disc impulse response places on the radial scaling index of the disc variability time-scale, $n_{t,d}$ [see equations (23) and (25) and Section 5.4].

Accretion disc emission is expected to take the form of a multitemperature blackbody, with radial temperature due to viscous dissipation alone scaling as $k_B T \propto r^{-3/4}$ for $r \gg r_{in}$ (Makishima et al. 1986). Since the disc radius for a fixed fluctuation signal frequency $v_{fluc}$ scales as $v_{fluc}^{-1/(n_{t,d}+3/2)}$, a given range of fluctuation frequencies should cover a greater range of disc temperature for smaller $n_{t,d}$ and thus show a stronger dependence of the disc spectral-timing properties on observed energy. In particular we expect to see a maximum PL-D hard lag which depends on photon energy and $n_{t,d}$, as the exponential cut-off of the blackbody spectrum limits the emission at a given energy to radii above a certain temperature, and hence limits the extent of the impulse response towards negative delays.

To determine how the energy-dependence of disc lags depends on $n_{t,d}$, we calculate the energy-dependent disc impulse responses using the following approach. First, we use the radially dependent disc dissipation and reverberation quantities defined in Section 3.2 to define a time-averaged total intensity at the radius $r_j$:

$$I(r_j) = \frac{f_{\text{diss}}(r_j)}{2\pi r_j} + \frac{f_{\text{therm}} f_{c \to d}(r_j) (\langle L_{\text{s,diss}} \rangle + \langle L_h \rangle)}{1 - f_{\text{return}}} \quad (C1)$$

We then assume that the total intensity is proportional to the blackbody intensity at that radius, which further assumes that coronal emission reprocessed by the disc contributes to and is reprocessed as a blackbody at the local disc temperature. This latter assumption may be a reasonable approximation to the reprocessed spectrum at high disc densities (e.g. see figures and discussion in Ross and Fabian 2007; García et al. 2016). Thus, we can define the radial disc temperature $k_B T(r_j)$ relative to the maximum radial disc temperature $k_B T_{max}$ (which is observable) as

$$k_B T(r_j) = k_B T_{\max} \left( \frac{I(r_j)}{\max I(r_j)} \right)^{\frac{1}{4}} \quad (C2)$$

For a given blackbody temperature $k_B T$ we can calculate what fraction of the blackbody photons are emitted within an energy range $E_1-E_2$, which for the constant component of disc emission is

$$f_{\text{bb,const}} = \frac{1}{2\zeta(3)(k_B T)^3} \int_{E_1}^{E_2} \frac{E^2}{(\exp(E/k_B T) - 1)} dE, \quad (C3)$$

where the term inside the integral is the blackbody photon rate spectrum $N_{bb}(E)$ and $\zeta(3) = 1.2020569...$ is the Riemann zeta function for argument 3, also known as Apéry's constant. Note that we calculate the fraction of total photon rate due to photons between $E_1, E_2$ rather than the fraction of luminosity in this range, because this is closer to the observed quantity which is usually a photon count rate.

For the disc impulse response used to calculate the time-variation of the emission, we cannot use the time-averaged blackbody spectrum, because although the blackbody bolometric intensity scales linearly with e.g. mass accretion rate or illuminating coronal flux, the blackbody spectrum changes shape in response to the corresponding temperature change. We must therefore obtain the linear change in blackbody photon flux spectrum $\delta N_{bb}(E) \approx \frac{dN_{bb}(E)}{d(k_B T)} \delta(k_B T)$ in response to a small temperature change $\delta(k_B T) = \frac{1}{4} k_B T \frac{\delta I}{I}$, where $\frac{\delta I}{I}$ is the relative change in bolometric intensity.<fnref n="6" /> The fraction of photons from this variable component contained within $E_1-E_2$ is then

$$f_{\text{bb, var}} = \frac{1}{6\zeta(3)(k_B T)^4} \int_{E_1}^{E_2} \frac{E^3 \exp(E/k_B T)}{(\exp(E/k_B T) - 1)^2} dE. \quad (C4)$$

With the above equation we can write the impulse response for the disc photon rate at radius $r_j$, with total time-delay $\tau_j$ and propagation delay across the radial bin $d\tau_j$, as

$$g_{\text{d,diss}}(\tau_j, E_1, E_2) d\tau_j = \frac{C_{\text{bb, var}} f_{\text{bb, var}}(r_j, E_1, E_2)}{k_B T(r_j)} g_{\text{d,diss}}(\tau_j) d\tau_j, \quad (C5)$$

where $C_{\text{bb, var}} = \frac{45\zeta(3)}{2\pi^4}$ and $C_{\text{bb, var}}/k_B T(r_j)$ corrects the bolometric luminosity of the impulse response to a scaled total photon rate, which is then converted by $f_{\text{bb, var}}$ to the rate between $E_1, E_2$.

Using a similar approach we also obtain the reverberation impulse response, which must be weighted and summed over all disc radii for each impulse response bin of the illuminating coronal emission:

$$g_{\text{d,rev}}(\tau_j, E_1, E_2) d\tau_j = \frac{C_{\text{bb, var}} f_{\text{therm}} (g_{\text{s,diss}}(\tau_j) + g_{\text{h}}(\tau_j))}{1 - f_{\text{return}}} \times \sum_i \frac{f_{\text{bb, var}}(r_i, E_1, E_2) A_i f_{c \to d}(r_i) (1 - f_{d \to c}(r_i))}{k_B T(r_i)} \quad (C6)$$

<footnote n="6">⁶See van Paradijs & Lewin (1986) for a related calculation of the change in blackbody flux spectrum seen in a thermonuclear X-ray burst.</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure C1.** Energy-dependence of power-law versus disc lags for an inverted cone corona with $n_{t,d} = 0$ and 2 (upper and lower panels) and for comparison (centre panel) the MAXI J1820 + 070 observed lags for three different energy bands covering the accretion disc emission, obtained from NICER ObsID 1200120104 (see Section 1 for further details).</figcaption></figure>

where $A_i$ is the area of the disc annulus for radial bin $r_i$, which is required because $f_{c→d}$ is the fraction of coronal emission intercepted by the disc *per unit disc area*. By adding together the above impulse responses we can obtain the impulse response of disc photon rate between $E_1, E_2$.

In Fig. C1, we show in the upper and lower panels the calculated PL-D lags for disc energies (relative to the maximum temperature) $2.6k_B T_{max}$, $3.5k_B T_{max}$, and $4.7k_B T_{max}$, as well as the bolometric luminosity $L_d$ used elsewhere in this paper. Lags are calculated

relative to a power-law energy $4.9E_s$. These energies all correspond to the median detected photon energies in each energy band used for the observational PL-D lag measurements shown in the centre panel, assuming a peak disc temperature $k_B T_{max} = 0.17$ keV for this observation (Bollemeijer private communication) and mono-energetic seed photons with energy $2.82k_B T_{max}$. The lags are calculated for two values of disc time-scale radial scaling index, $n_{t,d} = 0$ and 2 (corresponding respectively to constant aspect ratio $h/r$ and constant scale height $h$), assuming an inverted cone corona with $r_{cor} = 10 R_g$, $h_{cor} = 10 R_g$, $\theta_{cor} = 30^{\circ}$, and parameters as given in Section 4 but with $s_t = 100$ and $s_t = 1$ respectively, such that the disc variability time-scales are identical at $r_{cor} = 10 R_g$.

The $n_{t,d} = 0$ case shows a distinct energy-dependent low-frequency flattening which occurs at substantially higher frequencies than the minimum signal frequency produced in the disc. This is because, for this radial dependence of time-scales, at these energies only the inner variable part of the disc is emitting, and the extent of the variable emitting region (and hence the impulse response) is also energy-dependent. For $n_{t,d} = 2$, however, the lags continue to rise down to the lowest variable signal frequency in the disc, since the varying region is more compact for this radial scaling.<fnref n="7" /> Although some differences remain, the MAXI J1820 + 070 data appear to be much more consistent with the lag behaviour for $n_{t,d} = 2$, even though the model has not been fitted to match the data in any way.

## APPENDIX D: CONSISTENCY WITH THE HARD STATE X-RAY SPECTRAL ENERGY DISTRIBUTION

The model we present here does not consider the detailed energy dependence of spectral-timing properties and consistency with the time-averaged X-ray spectrum. However, it is still possible to compare model parameters with our expectations about the X-ray SEDs of BHXRBs in the hard state. In Section 1 (Fig. 1) and Appendix C (Fig. C1), we show NICER spectral-timing data from an observation (ObsID 1200120104) of MAXI J1820 + 070 in the hard state, which was obtained early in the 2018 outburst. Fig. D1 shows the corresponding unfolded NICER X-ray spectrum for this observation, obtained by fitting the spectrum using XSPEC 12.14 with a model consisting of a disc blackbody (DISKBB, Mitsuda et al. 1984) plus Comptonized emission (NTHCOMP; Zdziarski, Johnson and Magdziarz 1996; Życki, Done and Smith 1999) and two relativistically broadened reflection components (RELCONV\*XILLVERCP, Dauser et al. 2014; García et al. 2014) which represent a combination of close-in (radial emissivity index $\eta = 4$) and more distant ($\eta = 1$) reflection, both with inner disc radius and black hole spin fixed at $10 R_g$, $a = 0.94$ respectively, to match the geometry assumed in our lags model (Appendix C). The entire source spectrum is absorbed by neutral Galactic absorption (TBABS; Wilms, Allen and McCray 2000) with column density fixed $N_H = 1.4 \times 10^{21} \text{ cm}^{-2}$. The Comptonized seed photon spectrum was assumed to be a single-temperature blackbody with temperature $kT_{BB}$ set to equal that of the fitted disc blackbody inner disc temperature. A single-temperature seed is more appropriate for the compact coronae considered here, since the majority of seed photons originate from the inner disc. The Comptonization cut-off temperature $kT_e$ is fixed at 60 keV, to produce a high-energy spectrum consistent with that reported by Zdziarski et al. 2021 (see their Fig. 2) using INTEGRAL

<footnote n="7">⁷The fluctuation frequency at the largest variable radius, $v_{fluc,max} = 0.03$ Hz used here, corresponds to $r_{fluc,max} = 46 R_g$ for $n_{t,d} = 2$, versus $r_{fluc,max} = 355 R_g$ for $n_{t,d} = 0$.</footnote>
<figure><img src="image_16.png" /></figure>

**Figure D1.** NICER X-ray spectrum for MAXI J1820 + 070 observed on 2018 March 15 while the source was in the hard state (see Figs 1 and C1 for corresponding spectral-timing data). The spectrum has been fitted with a blackbody plus Comptonized blackbody spectrum and two reflection components, all absorbed by a neutral absorbing column (see text for details) and unfolded to remove the effects of the instrument response (care should therefore be taken when interpreting the sharpness of features in the spectrum). The solid lines show the fitted data (in grey) and (absorbed) model components in the range covered by the data, while the dashed lines show the extrapolation of these (unabsorbed) components over a broad range of the SED produced by the accretion flow.

high-energy X-ray data which overlaps in time with the March 18 NICER data. The reflection spectra have continuum shapes with photon index Γ and $kT_e$ tied to the values for the Comptonizing component (with Γ then left free to find the best-fitting value), while reflection normalization is left as a free parameter.

The goal of our model fit is to investigate the likely disc contribution to the luminosity rather than provide a serious model fit; however, for completeness, we note that with the default 1.5 per cent systematic errors included, the model provides an excessively good fit to the data ($\chi^2 = 735$ for 943 degrees of freedom). The fitted spectrum shows a disc blackbody with $kT_{BB} = 0.18$ keV, and Comptonization Γ = 1.64, consistent with the moderately luminous hard state shown here. It is important to note that although the current version of XILLVERCP can in principle model the reflection spectra for discs of hydrogen density $N > 10^{18}$ cm$^{-3}$ (García et al. 2016), consistent with BHXRB discs, these high densities produce significant emission features and sharp edges below 1 keV, which are strongly disfavoured in the fit, due to the smoothness of the soft X-ray excess which is better fitted by the disc blackbody. This is not surprising since the current XILLVERCP models assume illumination by a Comptonized continuum with a disc blackbody seed spectrum with $kT = 0.01$ keV$^8$, which extends to much lower energies than the spectra considered here. Furthermore, these models do not yet account for the significant upward vertical blackbody flux expected in the inner disc of hard state BHXRBs with the truncation radii and luminosities considered here. Since disc blackbody and Comptonized emission alone provide a reasonable fit to the low energy part of the spectrum, we obtain good fits by assuming reflection model density $N = 10^{15}$ cm$^{-3}$ and low ionization parameter $\log(\xi) \simeq 0$ for both reflection components.

<figure><img src="image_17.png" /></figure>

**Figure D2.** Top panel: Dependence on coronal radius of the disc luminosity (total including reverberation, solid), dissipated disc luminosity only, dashed) expressed as a fraction of coronal luminosity (sum of seed and heating luminosities including returning seed luminosity) and calculated for the spherical and inverted cone ($h_{cor} = 10 R_g$, $\theta_{cor} = 30^\circ$) geometries considered here. The broad bands show the range of inferred values of total disc luminosity fraction after accounting for disc inclinations ranging from 60 to 80 deg (see the text for details). The dotted line shows the value estimated for MAXI J1820 + 070 on 2018 March 15, based on the spectral fit shown in Fig. D1. Lower panel: Dependence of time-averaged photon index $\langle\Gamma\rangle$ on $r_{cor}$, with the fitted MAXI J1820 + 070 result shown with a dotted line.

Fig. D2 shows the fraction of coronal luminosity contributed by the disc blackbody emission due to viscous dissipation and reverberation, plotted as a function of disc truncation radius $r_{cor}$, for the spherical corona and the inverted cone corona with $h_{cor} = 10 R_g$ and $\theta_{cor} = 30^\circ$, together with the corresponding mean photon index $\langle\Gamma\rangle$. Solid lines show the result of combining disc dissipation and reverberation luminosity, while dashed lines show the result only for disc dissipation. Considering the integrated 0.01–1000 keV flux of the fitted disc blackbody and Comptonization components (after removing absorption), we obtain $L_d/L_{cor} \simeq 0.09$ for the MAXI J1820 + 070 data shown in Fig. D1. This is considerably lower than expected from the isotropic luminosities calculated from our model. However, MAXI J1820 + 070 is considered to be a high-inclination system, with the inclination $i$ estimated by jet-orientation (Atri et al. 2020) and binary orbit (Torres et al. 2020) lying conservatively in the range $i \simeq 60–80^\circ$. The inclination should further suppress the disc emission relative to the corona by a factor up to $\cos(i) \simeq 0.17–0.5$, since even for the opaque coronae assumed here, the apparent coronal flux will not be significantly affected by the high inclination due to its vertical extent and greater spherical symmetry compared to the disc. The inclination effect on the inferred total disc fraction is shown by the coloured bands in Fig. D2. Finally, the disc fraction may be further overestimated by screening of the far side of the disc by the corona (an effect which will be enhanced at high inclinations), so that the data may be consistent with radii $r_{cor} \simeq 10–20 R_g$, as inferred from our lag model for the hard state. $\langle\Gamma\rangle$ is very strongly dependent on the coronal geometry, so there is scope for a wide range of coronal radii which could be consistent with the observed disc luminosity fraction and $\Gamma$.

<footnote n="8">8 http://www.sternwarte.uni-erlangen.de/~dauser/research/relxill/relxill_docu_v1.0.pdf</footnote>