Modelling the Energy Dependence of Black Hole Binary Flows
Abstract
We build a full spectraltiming model for the low/hard state of black hole binaries assuming that the spectrum of the Xray hot flow can be produced by two Comptonisation zones. Slow fluctuations generated at the largest radii/softest spectral region of the flow propagate down to modulate the faster fluctuations produced in the spectrally harder region close to the black hole. The observed spectrum and variability are produced by summing over all regions in the flow, including its emission reflected from the truncated disc. This produces energydependent Fourier lags qualitatively similar to those in the data. Given a viscous frequency prescription, the model predicts Fourier power spectral densities and lags for any energy bands. We apply this model to archival RXTE data from Cyg X1, using the timeaveraged energy spectrum together with an assumed emissivity to set the radial bounds of the soft and hard Comptonisation regions. We find that the power spectra cannot be described by any smooth model of generating fluctuations, instead requiring that there are specific radii in the flow where noise is preferentially produced. We also find fluctuation damping between spectrally distinct regions is required to prevent all the variability power generated at large radii being propagated into the inner regions. Even with these additions, we can either the power spectra at each energy, or the lags between energy bands, but not both. We conclude that either the spectra are more complex than two zone models, or that other processes are important in forming the variability.
keywords:
accretion, accretion discs – Xrays: binaries – Xrays: individual: Cygnus X11 Introduction
Black hole binaries (BHBs) show variability on a wide range of timescales. Over days, months and years, mass accretion rate changes drive changes in the energy spectrum. The most dramatic example of this behaviour is the spectral transition from the Comptonisationdominated (low/hard) spectra seen at low luminosities to the discdominated (high/soft) spectra at high luminosities (see e.g. Remillard & McClintock 2006). This has a very natural interpretation from the two stable solutions to the accretion flow equations: one which is hot, optically thin and geometrically thick (advection dominated accretion flow, ADAF; Narayan & Yi 1995) which can only exist at low mass accretion rates, and one which is cool, optically thick and geometrically thin (ShakuraSunyaev disc, SS; Shakura & Sunyaev 1973). The observed switch in spectral properties can therefore be explained by a switch between these two solutions at the maximum ADAF luminosity (Esin, McClintock & Narayan 1997; Done, Gierliński & Kubota 2007; hereafter DGK07).
There is also more subtle spectral evolution within the low/hard state. This can be explained by combining these two solutions into a composite structure, where the outer SS disc truncates at some radius to be replaced by a hot flow interior to this (truncated disc/hot flow models). In this geometry, decreasing the truncation radius as the mass accretion rate increases results in more disc seed photons incident upon the flow. This leads to more efficient Compton cooling of the flow, and naturally produces the softer Comptonised spectra with increasing luminosity, as observed (DGK07).
However, the broadband spectra show more complexity than the contributions from a simple truncated disc, a singletemperature Comptonisation region, and its reflection from that disc. This complexity can be fit by assuming that the Comptonisation is not at a singletemperature, but instead is radially stratified. Simple inhomogeneous flow models consisting of a truncated disc and two Comptonisation components can broadly fit the 0.2200 keV spectra seen in the low/hard states of Cyg X1 (Gierliński et al. 1997; Di Salvo et al. 2001; Makishima et al. 2008; Yamada et al. 2013; Basak et al. 2017).
Alternatively, the additional Xray component in low/hard state spectra can instead be modelled by a jet contribution (Markoff, Nowak & Wilms 2005; Nowak et al. 2011), or using the completely different geometry of an untruncated disc with highly relativistic reflection from a pointsource on the spin axis of the black hole (e.g. Rykoff et al. 2007; Reis, Fabian & Miller 2010; Fabian et al. 2014). Variations on the theme of the truncated disc/hot flow model where some of the hot flow electrons have a hybrid (thermal/nonthermal) electron distribution have also successfully fit the spectra (e.g. Poutanen & Coppi 1998; Ibragimov et al. 2005; Makishima et al. 2008; Poutanen & Vurm 2009; Nowak et al. 2011). However, the fast timing (0.01100 s) properties can break some of these spectral degeneracies by giving additional information on the source geometry. In particular, the evolution of the power spectral density (PSD) of the fast timing variability in the low/hard state strongly supports the truncated disc/hot flow geometry (Ingram & Done 2011, hereafter ID11). The PSD of the Comptondominated Xray emission shows bandlimited noise between low () and high () frequency breaks, often accompanied by a strong lowfrequency quasiperiodic oscillation (QPO) at . Both and increase together as the spectrum softens towards the transition (Wijnands & van der Klis 1999; KleinWolt & van der Klis 2008; Rapisarda, Ingram & van der Klis 2014), indicative of the decreasing characteristic radius predicted by the truncated disc models.
This geometry can be incorporated into a full timing model by assuming that density fluctuations are generated at all radii in the hot flow by the turbulent magnetic dynamo (magnetorotational instability: MRI, Balbus & Hawley 1998). These fluctuations propagate inwards on the viscous timescale, so that slow fluctuations stirred up at large radii modulate the faster fluctuations generated at smaller radii (Lyubarskii 1997). This process can reproduce the observed doublebroken power law shape of the low/hard state PSD, while Arévalo & Uttley (2006; hereafter AU06) also show that this behaviour is necessary and sufficient to produce the observed linear rmsflux relation. The correlated QPO can also be produced from the same geometry if the entire hot flow undergoes LenseThirring precession due to its misalignment with the black hole spin axis (Fragile & Meier 2009; Ingram et al. 2009; Liska et al. 2017). These propagating fluctuation/LenseThirring precession models have quantitatively fit the data from XTE J1550584 during its spectral transition, with the inner radius of the thin disc changing from (ID11; Ingram & Done 2012a, hereafter ID12a), while also correctly predicting the modulation of the iron line energy on the QPO period (Ingram & Done 2012b; Ingram et al. 2016).
Thus, the overall properties of the power spectra already strongly favour the truncated disc/hot flow geometry for the low/hard state, but they do not break the degeneracies between the different models for the Xray emission within this framework. However frequencydependent time lags are also observed between high and lowenergy Xray bands (Miyamoto & Kitamoto 1989; Nowak et al. 1999). These were first discovered in data above 2 keV, so they directly probe the structure of the hot flow rather than the disc emission. The lags show that flux variations are seen first in the softer Xrays, and later in the hard Xrays (hard lags), after a lag time which depends on the fluctuation frequency. This frequencydependence rules out a simple light travel time origin for the signal, such as the delay between successive Compton scattering orders, as this would produce a constant, very short hard lag (Miyamoto & Kitamoto 1989; Nowak et al. 1999). The light travel time between the source and disc in the reflection dominated spectra is also ruled out as the source of the hard lag, as this process results in only a constant, very short soft lag.
Instead, the observed frequencydependent hard lags can be qualitatively explained by the propagation of fluctuations through an inhomogeneous hot flow such as the one we have described, where the Compton spectrum is harder closer to the black hole (Kotov, Churazov & Gilfanov 2001). This results in a coupled spectraltiming model, where slower fluctuations are produced at larger radii, so have softer spectra. These fluctuations propagate down to smaller radii, modulating the harder spectra from these regions. The lag time for this propagation encodes the viscous timescale between the radii. Faster fluctuations are produced at smaller radii, so they have a shorter distance to propagate before they modulate the hardest spectra from the innermost region, giving the frequency dependent lag time (AU06). By contrast, in a model where the soft Xrays are from the jet, fluctuations would propagate down through the accretion flow which produces hard Comptonisation, and only then propagate up into the jet which produces soft Xray synchrotron. This would instead predict a soft lag, contrary to what is observed.
Hence the spectral lags impose additional constraints on the physical nature and geometry of the hot flow. Here we build a fully energydependent spectraltiming model of the simplest possible inhomogeneous hot flow interior to a truncated disc, where the flow is composed of only two Comptonisation regions of different temperature and optical depth. We quantitatively compare the predictions of this model to the best fast spectraltiming data currently available: the archival Rossi Xray Timing Explorer (RXTE) observations of Cygnus X1 in the low/hard state (see 2), as used by Nowak et al. (1999). The data span an energy range of 330 keV, so they are dominated by the emission from the hot flow, and exclude the truncated disc emission.
Model fitting to Xray lags has only recently become possible, and our approach is complementary to the few papers produced so far on this. Rapisarda et al. (2016, hereafter R16) use lower energy data (0.510 keV) from Swift to model the power spectra in a soft and hard band, and the lags between them for the BHB MAXI J1659152. This lower energy band means that they consider the intrinsic disc emission and its variability (Uttley et al. 2011) and how this propagates into the hot flow, which they assume is homogeneous. This model is adequate to describe that dataset as it does not extend above 10 keV. However, the same disc and homogeneous hot flow model fails to fit the RXTE data from XTE J1550564 (Rapisarda, Ingram & van der Klis 2017, hereafter R17), potentially because the higher energy range (230 keV) of these data mean that the crossspectral properties are sensitive to the structure within the hot flow, and such structure is not incorporated into their model.
We describe the data we compare to in 2, while 3 briefly details the single zone propagating fluctuation model. In 4  9 we systematically build our procedure to predict frequencydependent time lags, by applying different spectral components to different radial ranges in the propagating fluctuations model. Finally in 10, we discuss the successes and failures of our model prescription and directly tie these back to the nature and geometry of the Xray emission region close to the black hole.
2 Observations of Cygnus X1 in the Hard State
Cygnus X1 is typically the brightest low/hard state source, and so gives the best data for studies using high time resolution. The archival data from RXTE remains the best publically available data for studying the Comptonisation lags, due to its high effective area in the 330 keV bandpass. Many of the RXTE observations were taken in a mode with limited spectral resolution below 10 keV. However, there are 6 datasets taken in the ‘Generic Binned’ mode which has 15.6 millisecond time resolution with 64 energy bins across the entire RXTE PCA energy bandpass (standard channels 0249; B_16ms_64M_0_249 configuration) giving reasonable spectral resolution in the 310 keV band, which allows the broad iron line to be resolved (Revnivtsev, Gilfanov & Churazov 1999; Gilfanov, Churazov & Revnivtsev 2000).
We use three of these observations taken consecutively during 1996, with simultaneous data from the Proportional Counter Array (PCA) and the High Energy XRay Timing Experiment (HEXTE; ObsIDs: 10238010800, 102380107000, 10238010700, hereafter observations 13). We choose these as they have very similar time averaged spectra, with hardness ratios between the 610 and 36 keV bands of , and respectively. The remaining three observations in this mode are all somewhat softer, so we exclude them. All 5 Proportional Counter Units (PCUs) of the PCA were active during these epochs. Each observation is backgroundsubtracted (using background on 16 s time binning), Poisson noise is removed, and deadtime corrections are applied according to the standard procedure of Nowak et al. (1999).
Observations 13 also have statistically consistent power spectra at the level across the entire frequency range, so we coadd these observations to give 22.5ks of data for the timing analysis. However we use only Obs. 1 for spectral analysis, as the coaddition of spectral data with slightly different response matrices can lead to artefacts.
Even amongst observations restricted to the hard state, a range of ‘substates’ are seen in both the variability and the spectra (e.g. the hardintermediate state; DGK07). We would therefore like to place our observations in the wider context of states seen from Cygnus X1. Grinberg et al. (2014; hereafter G14) fit all the Cyg X1 data taken during the lifetime of RXTE with a phenomenological model of tbabs*(gaussian + highecut*bknpower), where the bknpower component approximates the Comptonised emission as a broken power law, parameterised by “soft" and “hard" photon indices, and respectively. Our data has a “soft" photon index of , which is the minimum found by G14, showing that this is one of the hardest states of Cyg X1 observed by RXTE. This extreme hard state is confirmed by the high fractional rootmeansquare variability (MuñozDarias, Motta & Belloni 2011; Heil, Vaughan & Uttley 2012) in the 215 keV band of .
For our analysis we use lightcurves in three energy bands: Low (3.134.98 keV), Mid (9.9420.09 keV) and High (20.0934.61 keV). We extract these using saextrct, ensemble averaging over 174 segments of 128 s length to derive power spectra and time lags which are far better constrained at high frequencies than previous modelcomparison studies (R16; R17).
3 The Propagating Fluctuations Model
The magnetorotational instability (MRI) threading the flow generates fluctuations in all quantities and on all timescales (Balbus & Hawley 1998). The stochastic variations in mass accretion rate propagate down through the Comptonising region, modulating all the faster fluctuations produced further in. We simulate a Comptonisation region extending from the thin disc truncation radius, , to the inner edge of the hot flow at , where all size scales are in units of . The flow is split into annuli, characterized by radius and width , logarithmically spaced such that is constant (ID11).
The largest amplitude fluctuations produced by any given radius have size , where is the thickness of the flow. For , this sets the local viscous time, , as the shortest timescale on which density fluctuations are generated at ; fluctuations on shorter timescales than this are damped by the response of the flow. This results in a break in the power spectrum of mass accretion rate fluctuations generated at of (Kotov, Churazov & Gilfanov 2001). The largest radius in the flow generates the slowest fluctuations, so the lowfrequency break in the observed PSD is .
However, translating this to an outer radius requires a functional form for the viscous timescale. This form is not yet clear. General Relativistic MagnetoHydrodynamical (GRMHD) simulations of the MRI currently predict that fluctuations can be generated on tentimes the Keplerian timescale, (Hogg & Reynolds 2017). However this predicts that the typical lowfrequency break seen in hardstate power spectra at Hz is produced by material at large distances, of order several hundred . This is in tension with results from spectral fitting to the iron line profile, which generally point to (Kolehmainen, Done & Diaz Trigo et al. 2014; Basak et al. 2017). This inconsistency is likely due to the limited physics currently incorporated into the GRMHD simulations. Typically these neglect radiative processes and the interface between the disc and hot flow (e.g. Liska et al. 2017). Until better simulations are available, we instead use a parameterised prescription where (ID11). For a thin SS disc, and (with ) while a selfsimilar ADAF adheres to the same scalings but with . More complex flows have , for example when including transonic effects in an ADAF (Narayan, Kato & Honma 1997) or in full MRI simulations (Fragile & Meier 2009).
We follow ID11, who determine and from fitting to the well known relationship between and (Wijnands & van der Klis 1999). This gives and , assuming that the QPO is indeed from LenseThirring precession of the entire hot flow and that . This is a simpler prescription than using a GRMHD surface density to derive as in ID12a and R16/17, and avoids the associated simulation uncertainties.
We assume that these stochastic mass accretion rate fluctuations are generated at each radius with random phase, but with a well defined power spectrum which is a zerocentered Lorentzian with a cutoff at ,
(1) 
where a tilde denotes the Fourier transform. The sinusoidal term on the right hand side describes the suppression of variability due to the time binning. The normalisation of our Lorentzian is selected such that all have a mean of and fractional variability , where and are the number of annuli and fractional variability generated per radial decade respectively.
Beginning at the outermost annulus, , we generate mass accretion rate fluctuations in the time domain, using the algorithm of Timmer and König (1996). For the outermost annulus we designate the accretion rate across the annulus as where is the mean mass accretion rate. These fluctuations propagate in to the next annulus, traveling a distance , which takes a time .
The response of the flow acts to smooth fluctuations on the lag timescale. We implement this via a moving average over the lag time across the light curve, such that the smoothed mass accretion rate is
(2) 
Taken together with time lags, the total propagated mass accretion rate function in the annulus is then
(3) 
until the annulus which is . These mass accretion rate functions are the fundamental quantity in several previous studies (e.g. AU06; ID11; ID12a) which accurately replicate the broken power law shape in BHB power spectra.
These works conventionally convert the mass accretion rate curves to light curves via where is the emissivity at annulus which can be parameterised in a number of ways depending on the assumptions made regarding energy dissipation. Instead, a key extension of our work is that the total energy dissipation is set by the gravitational energy release, with the photon energy dependence set by the different spectra generated at different radii.
4 Incorporating Energy Dependence
4.1 Spectral decomposition
The standard propagating fluctuations model assumes a constant spectral energy distribution (SED) across the entire hot flow. It is only the normalisation of this SED which varies in time according to the variability of the flow at each radius, while the shape is assumed to be invariant. However physically there is more energy from gravitational heating of the flow at smaller radii, and fewer seed photons cooling it, so we expect the inner regions to have higher temperatures and hence harder spectra (Poutanen & Veledina 2014). Since the viscous frequency is also a function of radius, this couples the spectral and timing properties so that the crossspectral statistics can be derived. This allows us to jointly compare the PSDs and time lags as a function of energy band as described by the propagating fluctuations model, simultaneously with the timeaveraged SED.
The simplest multicomponent flow is described by two main Comptonisation regions: one softer component close to the disc, and one harder close to the black hole (see Fig. 1). We therefore fit the time averaged SED, with 0.5% systematic errors, in xspec (version 12.9.1; Arnaud, Borkowski & Harrington 1996) with two Comptonisation components described by tbabs*(nthcomp+nthcomp) (Zdziarski, Johnson & Magdziarz 1996), and the combined reflection of these, tbabs*(kdblur*xilconv*twocomp). Here twocomp is a local model which adds the Comptonisation components together, so that reflection is explicitly calculated from the composite spectrum. Such a decomposition is motivated both by model simplicity, and by similar successful fits to Cyg X1 spectra (Gierliński et al. 1997; Di Salvo et al. 2001; Makishima et al. 2008; Basak et al. 2017). We also follow Makishima et al. (2008) and assume that both Compton components have the same electron temperature. The data and best fit model are shown in Fig. 2, with full parameters detailed in Table 1. The softer and harder Comptonisation components, (green) and (cyan), originate from the outer and inner regions of the flow respectively. Also included is the total reflection from the disc, , but we do not include the intrinsic or reprocessed disc emission as the energy of this is too low to make a significant contribution to the RXTE data above 3 keV.
Component  Parameter  Value 

nthComp  
(keV)  
norm  
nthComp  
(keV)  
norm  
xilconv  relative refl norm  
log() 
4.2 Spectraltiming model
In all simulations we assume that Cyg X1 has a black hole of mass, , and a dimensionless spin parameter of 0.85 (Kawano et al. 2017). The inner radius is set to the approximate ISCO size implied from the spin of Cyg X1, so that .
The timeaveraged spectrum, , emitted from each radius is given by the expression
(4) 
here is the transition radius between the soft and hard Comptonisation regions. This is analytically derived from an assumed emissivity, , where is an inner boundary condition, such that the luminosity ratio between the two components matches that observed, such that
(5) 
The light curves produced by the standard propagating fluctuations model at each radius are then made energydependent and renormalised such that their timeaverage is the flux for that energy bin and radius, yielding
(6) 
The summation limits implied by ‘region’ are to for and to for . This normalisation guarantees that if Eq. 6 is time averaged and summed over all radii, the observed energy spectrum is reproduced.
We match our spectra to the data as closely as possible by converting these fluxes to count rates using the detector effective area and galactic absorption . The count rate is then expressed
(7) 
where is the Thompson crosssection.
In practice, Eq. 7 describes a threedimensional matrix, which can be operated on in different ways to obtain a variety of statistics. For instance, the total count rate in each energy band can be obtained by summing the matrix in Eq. 7 over all radii, and over the energy band of interest, yielding
(8) 
We use this quantity to produce the model power spectral and crossspectral statistics, which we then fit to their analogues from the data.
We first use the viscous model of ID11, i.e. a frequency prescription with and as discussed in 3. Tying the viscous frequency at the outer radius to the lowfrequency break in the data so that Hz , we obtain an outer radius of , which is consistent with the range of disc truncation radii found from spectral fitting of 1320 (Basak et al. 2017). The fiducial model of ID11 had an emissivity described by and a stressed inner boundary condition, . Coupling this with our decomposition of the time averaged spectrum through Eq. 5 gives .
We calculate the light curves on a time binning of ms (matched to the timing mode resolution of RXTE) and simulate s for each realisation, ensemble averaging over realisations. All simulations use radial bins, and we require in order to match the slope and amplitude of the lowfrequency break. A summary of all parameter values used in the simulations in this paper can be found in Table 2.
Fig. 3 shows the model PSD from this simulation, where it is clear that this a poor match to the data. Overall, all energy bands show far too little highfrequency power. The model also predicts that the PSDs of all energy bands are similar, while the data shows that the Low band dominates at all frequencies below Hz. We analytically explore the factors determining the shapes of these PSDs below.
5 Analytic power spectral models
The pioneering work of Ingram & van der Klis (2013, hereafter IK13) show how the PSD can be analytically calculated by considering how propagated PSDs are constructed in Fourierspace, and how they are weighted by the emissivity in calculating the final count spectra. We will now adapt their procedure to reflect our simulations, including the light curve weightings according to the energy spectrum.
In the following we denote the PSDs generated in annulus , as while those which are propagated from all outer annuli down to are denoted . Propagation causes the noise in (closer to the black hole) to be modulated by the noise in , lagged by the viscous timescale. Since this lag time is small compared to the generator timescale in , then it is almost perfectly coherent between and so that the power is additive, and . Generalising this, the propagated PSDs are described by
(9) 
where the assumed selfsimilar nature of the fluctuations means that all the individual have the same amplitude. Fig. 4a (dashed lines) shows the generator PSDs of the individual annuli, with the soft region in green and the hard in cyan. The solid green and cyan lines show the propagated PSDs of the total soft and hard regions respectively. This shows the clear difference in highfrequency extent of the PSDs of the two regions, with the hard region producing substantial additional power above Hz.
Our mass accretion rates are converted to counts in a given band using the emissivity prescription and SED decomposition described in 4.1. This effectively weights the propagated mass accretion rate from each annulus by a factor, , given by
(10) 
The count spectrum for that band can then be written
(11) 
Since the mean count rate of is normalised to , the mean count rate in a given energy band is then
(12) 
We now drop the superscript on for notational convenience, and take the rmsnormalised PSD:
(13)  
Including decoherence due to the propagation lag results in the full PSD expression of
(14)  
where the weights now have a spectral dependence in our case. The second term in Eq. 14 arises since the propagated noise at interferes with the propagated power spectra found at all outer radii. If this noise were not lagged between radii, this interference would be purely constructive and the cosine term would reduce to unity for all frequencies. However the time lag causes a component of the PSDs to interfere destructively, and this suppression is expressed by the lagdependent cosine factor. Here, is the total time lag between two annuli so that
(15) 
where the second equality comes from the fact that the radial bin size is logarithmic, and is therefore a constant.
Eq. 14 shows that each bandspecific PSD is a weighted combination of the propagated PSDs in the outer and inner regions, where the weighting factor depends on the fraction of hard and soft spectral components in that band. Fig. 2 shows that the hard spectrum contributes more to the higher energy bands but is always a fairly small fraction of the total emission. Hence in Fig. 4b (dashed lines), the analytic PSD of each band is dominated by the soft region. This results in banddependent PSDs which are highly similar, with very little highfrequency power. We also show the full simulation output as the solid lines for comparison. Smoothing has been included in the simulations but is neglected in the analytic form, but the good agreement below Hz indicates that this effect is negligible.
It is clear that this combination of spectral decomposition, emissivity and viscosity will be unable to produce PSDs in each band which are close to the observations. Since the SED decomposition is fit prior to the spectraltiming model, we now explore how a better match to the power spectra can be achieved by varying the emissivity and viscous frequency prescription. We will start by reproducing the lowfrequency Lorentzian hump at 0.2 Hz.
6 Varying the emissivity and viscous frequency prescription
6.1 A physically motivated emissivity
Gravity gives an expected emissivity with , while the innermost stable circular orbit motivates a stressfree (SF) inner boundary condition which we approximate as . Together with the standard , viscosity prescription, we will hereafter refer to this parameter set as the fiducial model. In Figs. 4c and 4d we show the effect on the PSDs of applying this new emissivity, with a new set by the integrated SED components via Eq. 5. It is clear when comparing Figs. 4a and c, that power at high frequencies has been suppressed in both the soft (green) and hard (cyan) bands. This is because the new emissivity weights the emission to larger radii, so the highfrequency contribution to the variability from the smallest radii is decreased. We also see that the simulated soft region power is now even lower than that from the hard region, even at frequencies below Hz.
In Fig. 5a, we fit the model with this new emissivity to the data, using to match to the lowfrequency break amplitude, but keeping all other parameters the same. This matches very well to the lowfrequency PSD hump in the Mid and High bands, although it does not match the significantly higher amplitude of the Low band since this remains a total propagation model. However, the rest of the power spectrum is completely unmatched as the new, less centrally peaked emissivity means that more of the emission arises from larger radii, so the PSD is weighted more to lower frequencies. In effect, the emissivity defines an envelope which suppresses all power above 0.51 Hz (see the Appendix of AU06 for an energyindependent treatment).
6.2 Different viscous frequencies, same radial range
The viscous parameterisation of ID11 assumed so far gives a maximum possible peak frequency of Hz, although the finite width of the Lorentzians means that there is some power of even higher frequency generated near . However, this highfrequency variability is suppressed, as the emissivity profile prevents these radii from producing a significant proportion of the total luminosity.
Increasing the maximum viscous frequency associated with the flow from would instead allow the PSD to extend to higher frequencies, while maintaining a gravitational emissivity. By diverging from , , we will break the  relation, but we nevertheless explore this in order to better understand the effects of varying the viscous frequency prescription.
We first maintain the size scale of the region (, ) and emissivity ( with the SF boundary condition) so stays constant at , but we now fit and such that the PSD amplitudes at Hz are approximated. This yields and . This keeps tied to , but now gives Hz. We see in Fig. 5b that, although it cannot match the peak structure seen in the data, this viscosity prescription can produce the observed highfrequency power. However, it has no physical motivation.
6.3 ADAF viscous frequencies, large radial range
The transonic ADAF models do make physical predictions about , predicting and for (Narayan, Kato & Honma 1997). The very high ion temperature of the ADAF means that the sound speed and hence the radial velocity is high, so a much larger radial scale is required to produce the lowfrequency break observed. We find a best fit of and , which gives . The PSDs produced by this very different parameter set (Fig. 5c) are indeed equivalent in all essential features to those of the standard size scale assumed in Fig. 5b, due to the similar viscous frequency ranges spanned by the models.
This is a key degeneracy. Without any external information to set the viscous frequency prescription (such as assuming that the QPO is set by LenseThirring precession) then the size scale of the region cannot be determined from the PSD. The data do show time lags between bands, however, so we now explore whether those time lags can break this degeneracy.
7 Time Lags
So far we have only investigated which elements of the observed PSDs can be replicated by this energy dependent model. However crossspectral statistics including time lags can also be extracted from our simulations. These give additional information to that contained in the power spectrum; a good match to the PSDs does not necessarily imply a good fit to the crossspectral lag (and vice versa), so any complete energydependent model must match both the power spectra and energy spectrum simultaneously with the time lags.
Figs. 6ac show the time lags for each of the three power spectra shown in Figs. 5ac, all of which used the physically motivated emissivity with and the SF inner boundary condition. It is striking that the fiducial prescription (a: , , , ), which has an excellent match to the lowfrequency Mid and Highband power spectra (Fig. 5a), also has an excellent match to the lowfrequency lags (Fig. 6a). This frequency prescription came from fitting the relation in ID11, so the good match to the lowfrequency lag amplitude therefore provides additional support for the assumed LenseThirring origin of the QPO. However, this model completely fails to match the lags above 2 Hz, because frequencies with Hz are produced only in the hard region. For frequencies above 2 Hz, the variability contribution in both bands therefore comes entirely from the hard region, so there are no spectral lags even though the fluctuations themselves are lagged.
Conversely, the viscosity prescription which gives higher frequencies over the same size scale can mostly match the PSD (Fig. 5b) but underpredicts the lags at all frequencies (Fig. 6b). This underprediction in the lag occurs because the viscous speed in the flow goes as . Compared to the fiducial prescription, is now higher (so is faster) for all , resulting in shorter lags. This prescription does produce significant lags up to a higher frequency however. This is because we now have Hz, so fluctuations slower than this are found in both the soft and hard regions, giving measurable lags at these frequencies.
The larger size scale ADAF model (Fig. 6c) produces very similar lags to those of Fig. 6b, highlighting the fact that degeneracies on size scale can remain even when incorporating time lags. We illustrate here why this occurs using a physically intuitive derivation of the maximum lag between the High and Low bands, but in Appendix A we extend this to all frequencies using the formalism of IK13.
The radial velocity is not constant, so the raw time lag, . The lowest frequency component, , propagates down through the entire flow. Light curves are calculated by weightsumming over the flow, and all fluctuations therefore appear to initiate at some radius as seen in the Low band, and arrive some time later at some radius as seen in the High band. and are the emissivityweighted averages of all radii in the soft and hard regions respectively, so that
(16) 
The maximum raw lag is then the propagation time between these radii
(17)  
This lag is then diluted by the soft SED contribution in the High band and the hard SED contribution in the Low band (Uttley et al. 2014), so we obtain
(18) 
where is the ratio of integrated soft flux to integrated hard flux in the High band and is the ratio of integrated hard flux to integrated soft flux in the Low band. The predicted values for the maximum lags in Figs. 6a, b and c (indicated by the solid magenta lines) are 0.09, 0.009 and 0.009 s respectively, which are generally consistent with the simulation results. Higher frequencies demand a full crossspectral treatment since they are more prone to interference (see Appendix A and IK13), however this simple result confirms that two very different size scales/viscosity prescriptions can predict indistinguishable lags.
Fundamentally, the larger size scale and higher velocity of the ADAF prescription (Fig. 6c) is degenerate with the smaller size scale and lower velocity of the Fig. 6b prescription, as both are tuned to match the breaks in the PSD. This is a direct consequence of assuming that the fluctuations are generated and propagated on the timescale. However, there are some models which do not require that the propagation and generation timescales are the same. We explore these below.
7.1 Decoupling the timescales of propagation and generation
So far, we have assumed that fluctuations are produced on the same timescale at which they propagate. We now decouple the generating timescale from the propagation timescale in order to determine the effect this has on the time lags. Section 3 argued that the largestscale coherent fluctuations are generated over distances , so that the maximum coherent timescale was . However considering the flow over all azimuths means that a better estimate for the generating timescale would be the timescale of fluctuations with are coherent around the entire annulus, i.e. , rather than as assumed thus far. Eq. 1 then becomes
(19) 
while the propagation timescale remains .
We assume the fiducial source size with , and , and attempt to find a viscosity prescription which recovers a PSD similar to that in Fig. 5b (which has and ). As the generating frequency is now , instead of , we start by dividing by and leaving unchanged, but the altered effect of smoothing means that there is a better match for slightly different parameters, with and , and . Fig. 7 shows the PSD and lags from this model. Interestingly, this has an even worse match to the observed lags than Fig. 6b, because the lags are even shorter relative to the generating timescale than before.
An alternative approach one might consider would be to generate fluctuations on the 10 timescale suggested by GRMHD simulations (Fragile & Meier 2009; Hogg & Reynolds 2017), while still propagating on set by the relation. However, generating on a 10 timescale would require an inner flow radius much smaller than the ISCO size of 2.5 in order to produce the required highfrequency power, so we do not explore this here.
In summary, the only case which approaches both the lowfrequency PSD and lowfrequency lags is the fiducial model in Figs. 5a and 6a, whereby the generating and propagation timescales are set equal to , derived from assuming that the QPO originates as LenseThirring precession of the hot flow. However, this model fails to explain the observed power at higher frequencies. This higher frequency power is also concentrated in a distinct ‘hump’ around 2 Hz, unlike the smooth PSD produced in the propagation models with faster viscous timescales. Pure propagation models with self similar fluctuations cannot produce such humps, and neither can they explain how the Lowenergy band can have more power than the Mid and High bands at most frequencies. We now explore a new family of models which allow the fractional variability to vary with radius, to see whether these can reproduce those essential features of the data.
8 Variability as a function of radius
The power spectra of the data are inherently ‘bumpy’, and this is a generic feature in both Cyg X1 (Churazov, Gilfanov & Revnivtsev 2001; Axelsson et al. 2008; Torii et al. 2011; G14) and other sources, such as GX 3394 (Nowak 2000). Veledina (2016) proposes an idealised model to explain the bumpy PSDs from the interference of two radially separated, lagged Compton continua. However our results so far have shown that it becomes much more difficult for interference to produce the observed peaks if we consider the extended nature of the source and the generation of fluctuations at all radii. Alternatively, R16 suggest that the hump structure can be produced by considering fluctuations in the truncated thin disc at . However the frequencies of these humps at or Hz are not easily consistent with any expected thin disc timescale. Instead the hump frequencies are more compatible with the viscous timescale within the flow itself. We therefore adapt our radially stratified model to allow enhanced variability at specific radii in the flow, in order to reproduce the observed PSD structure.
We keep the fiducial prescriptions for the viscous frequency (, , , ), and emissivity ( and the SF boundary condition) as these match well to the observed lowfrequency hump. is then allowed to vary with radius, so that the additional variability can be incorporated.
We first assume that there is enhanced turbulence at a specific radius , and derive this radius from the viscous frequency of the second peak in the PSD, i.e. Hz. From this we obtain , placing it at the inner edge of the soft region since . All annuli in the flow apart from the one containing have . The one which contains requires in order to match the amplitude of the 2 Hz peak in the Mid band PSD.
Fig. 8a shows that this additional power produces a divergence of the PSDs in different energy bands, reaching a maximum amplitude difference at Hz. However it is worth noting that since this is a loglog plot, the Mid and High bands are actually much less distinct than the Low. The divergence arises because is situated close to , so only a small fraction of the soft region is affected by this additional variability, whereas it all propagates through the hard region. The Mid and High bands sample mostly from the hard region, while the Low band samples mostly from the soft region, resulting in this deficiency in power at high frequencies in the Low band.
Adding variability at also does not reproduce the desired ‘hump’ structure at 2 Hz. Instead the model PSDs are smooth from to Hz. This is due to propagation, since the noise generated in the soft region propagates coherently to , and so adds constructively to the additional noise. To obtain the observed decrease in the PSD from Hz requires that fluctuations are damped as they propagate, even more strongly than the smoothing in Eq. 2 (see also R17).
The enhanced fluctuation power at also underpredicts the highfrequency power above Hz, which rises to a third hump at Hz. This third Lorentzian peak is commonly seen in the power spectra of Cyg X1 (Pottschmidt et al. 2003; Axelsson et al. 2008; Axelsson & Done, in preparation), and potentially in other sources (e.g. GX 3394; Nowak 2000) indicating that the process driving this additional noise may be a fundamental physical mechanism in the Comptonising region.
We therefore add a second enhanced variability component at with amplitude to match the third Lorenzian peak at 8 Hz. However, cannot now simply be derived from Hz due to the increased effects of interference. Instead we fit for this, and find a best fit to the Midband PSD for these parameters of and . The resulting PSDs are shown in Fig. 9. The Mid band power spectrum is now fairly well matched (apart from the dip between 0.31 Hz), as are the lags, but the Low and Highenergy PSDs are far from the observed data.
These results collectively support a model which relies on additional turbulence at characteristic positions in the flow to produce the highfrequency observed power. However certain key features have yet to be reproduced. The models so far have assumed that all variability from the outer regions is propagated, uninterrupted, into the inner regions such that Eq. 3 is applicable throughout the flow. However, the observed drop in power in the 0.31 Hz requires damping of the fluctuations. We now explore whether this can also finally reproduce the dominance of the Lowenergy power spectrum over the High and Mid bands, whilst maintaining the lags of Fig. 9.
9 Damped Soft Variability
The observed dominance of the Lowenergy power has been seen in previous studies of Cyg X1 (Grinberg et al. 2014), as well as in other BHBs including SWIFT J1753.50127 and GX 3394 (Wilkinson & Uttley 2009), so it is a generic feature of the data.
Total propagation has so far prevented the Lowenergy band from having more variability power at any frequency than the High band, since the lowfrequency variability power dominating the Low band modulates the highfrequency power which dominates the High band. To suppress the transfer of lowfrequency power to the High band then requires that some fraction of the variability power in the soft region fails to propagate into the hard region. However, we also require that the coherence between the fluctuations in each region is maintained. The soft variations must therefore map onto the innerregion variability after propagation, although with smaller amplitude and a time delay.
Physically, unpropagated noise could arise if part of the variability comes from disc seedphoton fluctuations. If the soft Comptonisation region becomes optically thick then it would shield the hard Comptonisation region from this variability component. Alternatively, part of this seed photon variability could be produced by a turbulent, clumpy transition between the truncated disc and hot flow, perhaps induced by instabilities in a shearing layer between the Keplerian disc and subKeplerian flow. These clumps might then evaporate, or be shredded by the MRI turbulence as they propagate inwards.
PSDs  Lags  

Fig. 3    4.5  1  0.03  0.5  14.  2.5  3.1      0.45  1  1   
Fig. 5a  Fig. 6a  3.  SF  0.03  0.5  14.  2.5  5.4      0.59  1  1   
Fig. 5b  Fig. 6b  3.  SF  250.00  3.95  14.  2.5  5.4      0.74  1  1   
Fig. 5c  Fig. 6c  3.  SF  94.87  1.21  140.  6.  16.0      0.49  1  1   
Fig. 7a  Fig. 7b  3.  SF  3.35  14.  2.5  5.4      0.69  1  1    
Fig. 8a  Fig. 8b  3.  SF  0.03  0.5  14.  2.5  5.4  5.5    0.52  14  1   
Fig. 9a  Fig. 9b  3.  SF  0.03  0.5  14.  2.5  5.4  5.5  2.7  0.52  9  166   
Fig. 10a  Fig. 10b  3.  SF  0.03  0.5  14.  2.5  5.4  5.5  2.6  60 
This particular case decouples the fluctuationgenerator and propagation timescales so that we still have but now Eq. 19 describes the generator Lorentzians.
Instead this is as described in the text of 9.
We model these effects generically by suppressing the amplitude of the propagated fluctuations by a factor, , at , and assume that the generated fluctuations within the hard region are also smaller by this factor than those generated in the soft region. For annuli without enhanced variability (, we therefore have .
We use the fiducial frequency prescription (, ) and size scale (, ) as this was the case which best approximated the observed lags. Fig. 10 shows the best fit found when the model is extended to include the free parameter, . An optimal fractional variability of is found on a simulation which also has additional variability at the two characteristic radii of and where and . The optimal suppression factor is found to be . In Fig. 11, we display the best fitting profiles of Sections 8 and 9 for ease of comparison. Using this parameterisation, we find the best approximation yet for the relative amplitudes and shapes of the PSDs in each energy band, although the difference in lowfrequency power between the Low and Mid/High bands is still slightly underestimated.
However, Fig. 10b shows that the simulated lags are now a poor match to the data, severely underestimating those which are observed, particularly at low frequencies. This is because the low frequency fluctuations are now highly damped, so they do not propagate sufficiently into the hard region for the crossspectral lags to be significant. However, this prescription does reproduce the shape of the 2 Hz ‘step’ in the lags, further suggesting the presence of specific radii in the flow which produce enhanced variability.
The other key shortfall of this model lies in the magnitudes of and . Large magnitude values such as these cause the generated light curves from these regions to go negative, which is clearly unphysical. Instead any future model demands smaller values in the regions of enhanced turbulence, so the emissivity in these regions must also be enhanced to transmit this smaller variability into the simulated light curves. This feature will be applied in future work (Mahmoud & Done, in preparation). However the results we have shown here stand as a proof of concept that a nonuniform radialvariability profile is a key element in the timing behaviour of BHBs.
10 Conclusions
We build a full spectraltiming model of fluctuations propagating down through a two component Comptonisation region in the BHB low/hard state. We systematically explore the effects of changing the model parameters on the energy dependent PSD and lags, and compare these to some of the best available data from Cyg X1. We have fit to data only above 3 keV so that it is dominated by the flow, not by the intrinsic disc emission. The main results of this study can be summarised as follows:

The viscous frequency parameterisation is degenerate with the radial size scale of the Comptonising region. Time lags do not break this degeneracy without some external constraints from estimates of the truncated disc radius e.g. from spectral fitting of the broad iron line, a LenseThirring origin of the QPO, and/or light travel time lags. All of these support the ID11 prescription with and , and so require that the low/hard state modelled here has an inner disc truncation radius of .

Coupling this to a standard emissivity with and a stressfree inner boundary condition alone cannot produce the observed PSD using these parameters from a selfsimilar propagation model. This emissivity weights the observed power strongly to larger radii and hence lower frequencies, such that the significant variability observed above 0.5 Hz cannot be produced by this viscosity prescription alone.

Additional highfrequency power can only be produced in these models by assuming that there is enhanced turbulence within the flow, varying as a function of radius. This is also likely required to replicate the ‘steps’ in the lagfrequency spectrum.

The PSD shape at all energies is emphatically nonmonotonic, with a distinct dip in variability power between a lowfrequency peak at 0.2 Hz and one at 2 Hz in these Cyg X1 data. This distinct dip cannot be produced in any pure propagation model, and requires that variability from the outer flow is damped at some characteristic radius (or radii; see also R17).

The commonly observed low/hard state feature of Lowenergy band dominance of the PSD at low frequencies requires that damping is included in the physical model. Some of the turbulence generated in the outer regions of the flow is not propagated down into the inner regions of the flow. However this damping also suppresses the lags.
Our work adds to a growing understanding that the Comptonising region found in hardstate BHBs  far from being spectrallyhomogeneous and smoothly variable  is almost certainly stratified in boths its spectrum and its variability (Wilkinson & Uttley 2009; Veledina 2016; R16; R17; Basak et al. 2017). Clearly there are specific radii in the flow at which the fluctuations are enhanced and/or damped. These could be physically associated with the bending wave radius from a misaligned spinning black hole (Fragile & Meier 2009; Ingram et al. 2009), and/or the radius at which the jet is launched. A better understanding of the PSD and lags mean that we should be able to observationally trace the radii at which this physics operates.
However, even with these additional model features, the energydependent PSDs and lags cannot be fit simultaneously with a twoCompton component spectral decomposition. Nonetheless, what we develop here is a flexible framework in which to construct a full spectraltiming model for the data. In future work we will modify this model to include more detailed spectral decompositions with three Compton components, which have been suggested by the most sophisticated spectral fits (Yamada et al. 2013). We will also explore the effect of introducing a distinct fluctuation timescale at the discflow interface, to better approximate the variability in the potentially unstable discflow transition layer.
Acknowledgements
RDM acknowledges the support of an STFC studentship. RDM and CD thank Magnus Axelsson for helpful discussions on the spectral decompositions and doublepeaked PSD. We also thank the anonymous referee for helpful comments which helped to improve the manuscript. This research has made use of data obtained through the High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center.
References
 Arévalo & Uttley (2006) Arévalo, P. & Uttley, P., 2006, MNRAS, 367, 801
 Arnaud, Borkowski & Harrington (1996) Arnaud K., Borkowski K.J. & Harrington J.P., 1996, ApJ, 462, L75
 Axelsson et al. (2008) Axelsson M., Hjalmarsdotter L., Borgonovo L. & Larsson S., 2008, A&A, 490, 253
 Balbus & Hawley (1998) Balbus S.A. & Hawley J.F., 1998, RvMP, 70, 1
 Basak & Zdiarski (2016) Basak R. & Zdiarski A.A., 2016, MNRAS, 458, 2199
 Basak et al. (2017) Basak R., Zdziarski A.A., Parker M., Islam N., 2017, https://arxiv.org/abs/1705.06638
 Churazov, Gilfanov & Revnivtsev (2001) Churazov E., Gilfanov M. & Revnivtsev M., 2001, MNRAS, 321, 759
 Di Salvo et al. (2001) Di Salvo T., Done C., ycki P.T., Burderi L. & Robba N.R., 2001, ApJ, 547, 1024
 Done, Gierliński & Kubota (2007) Done C., Gierliński M. & Kubota A., 2007, A&ARv, 15, 1
 Esin, McClintock & Narayan (1997) Esin A.A., McClintock J.E. & Narayan R., 1997, ApJ, 489 (2), 865
 Fabian et al. (2014) Fabian A.C., Parker M.L., Wilkins D.R., et al., 2014, MNRAS, 439, 2307
 Fragile & Meier (2009) Fragile P.C. & Meier D.L., 2009, ApJ, 693, 771
 Gierliński et al. (1997) Gierliński M., Zdziarski A.A., Done C., et al., 1997, MNRAS, 288, 958
 Gilfanov, Churazov & Revnivtsev (2000) Gilfanov M., Churazov E. & Revnivtsev M., 2000, MNRAS, 316, 923
 Grinberg et al. (2014) Grinberg V., Pottshmidt K., Böck M., et al., 2014, A&A, 565, A1
 Heil, Vaughan & Uttley (2012) Heil L.M., Vaughan S. & Uttley P., 2012, MNRAS, 422, 3620
 Hogg & Reynolds (2017) Hogg J.D. & Reynolds C.S., 2017, ApJ, 834, 80
 Ibragimov et al. (2005) Ibragimov A., Poutanen J., Gilfanov M., Zdziarski A.A. & Shrader C.R., 2005, MNRAS, 362, 1435
 Ingram et al. (2009) Ingram A., Done C. & Fragile P.C., 2009, MNRAS, 397 (1), L101
 Ingram & Done (2011) Ingram A. & Done C., 2011, MNRAS, 415 (3), 2323
 Ingram & Done (2012a) Ingram A. & Done C., 2012, MNRAS, 419, 2369
 Ingram & Done (2012b) Ingram A. & Done C., 2012, MNRAS, 427, 934
 Ingram & van der Klis (2013) Ingram A. & van der Klis M., 2013, MNRAS, 434, 1476
 Ingram et al. (2016) Ingram A., van der Klis M., Middleton M., et al., 2016, MNRAS, 461 (2), 1967
 Kawano et al. (2017) Kawano T., Done C., Yamada S., Takahashi H., Axelsson M. & Fukuzawa Y., 2017, PASJ, 69(2), 36
 KleinWolt & van der Klis (2008) KleinWolt M. & van der Klis M., 2008, ApJ, 675, 1407
 Kolehmainen, Done & Diaz Trigo et al. (2014) Kolehmainen M., Done C. & Diaz Trigo M., 2014, MNRAS, 437, 613
 Kotov, Churazov & Gilfanov (2001) Kotov O., Churazov E. & Gilfanov M., 2001, MNRAS, 327, 799
 Liska et al. (2017) Liska M., Hesp C., Tchekhovskoy A., Ingram A., van der Klis M., Markoff S., 2017, https://arxiv.org/abs/1707.06619
 Lyubarskii (1997) Lyubarskii Y.E., 1997, MNRAS, 292, 679
 Makishima et al. (2008) Makishima K., Takahashi H., Yamada S., et al., 2008, PASJ, 60, 585
 Markoff, Nowak & Wilms (2005) Markoff S., Nowak M.A. & Wilms J., 2005, ApJ, 635, 1203
 Miyamoto & Kitamoto (1989) Miyamoto A. & Kitamoto S., 1989, Nature, 342, 773
 MuñozDarias, Motta & Belloni (2011) MuñozDarias T., Motta S. & Belloni T.M., 2011, MNRAS, 410, 679
 Narayan, Kato & Honma (1997) Narayan R., Kato S. & Honma F., 1997, ApJ, 476, 49
 Narayan & Yi (1995) Narayan R. & Yi I., 1995, ApJ, 452, 710
 Nowak et al. (1999) Nowak M.A., Vaughan B.A., Wilms J., Dove J.B. & Begelman M.C., 1999, ApJ, 510, 874
 Nowak (2000) Nowak M.A., 2000, MNRAS, 318, 361
 Nowak et al. (2011) Nowak M.A., Hanke M., Trowbridge S.N., et al., 2011, ApJ, 728, 13
 Pottschmidt et al. (2003) Pottschmidt K., Wilms J., Nowak M.A., et al., 2003, A&A, 407, 1039
 Poutanen & Coppi (1998) Poutanen J. & Coppi P.S., 1998, Phys. Scripta, T77, 57
 Poutanen & Vurm (2009) Poutanen J. & Vurm I., 2009, ApJ, 690, L97
 Poutanen & Veledina (2014) Poutanen J. & Veledina A., 2014, Space Sci. Rev., 183, 61
 Rapisarda, Ingram & van der Klis (2014) Rapisarda S., Ingram A. & van der Klis M., 2014, MNRAS, 440, 2882
 Rapisarda et al. (2016) Rapisarda S., Ingram A., Kalamkar M. & van der Klis M., 2016, MNRAS, 462, 4078
 Rapisarda, Ingram & van der Klis (2017) Rapisarda S., Ingram A. & van der Klis M., 2017, MNRAS, 469 (2), 2017
 Reis, Fabian & Miller (2010) Reis R.C, Fabian A.C & Miller J.M., 2010, MNRAS, 402, 836
 Remillard & McClintock (2006) Remillard & McClintock, 2006, in Compact Stellar XRay Sources, Ch. 4, Cambridge University Press, ed. Lewin, W.H.G. & van der Klis, M.
 Revnivtsev, Gilfanov & Churazov (1999) Revnivtsev M., Gilfanov M. & Churazov E., 1999, A&A, 347, L23
 Rykoff et al. (2007) Rykoff, E.S., Miller J.M, Steeghs D., Torres M.A.P., 2007, ApJ, 666, 1129
 Shakura & Sunyaev (1973) Shakura N.I. & Sunyaev R.A., 1973, A&A, 24, 337
 Timmer & König (1996) Timmer J. & König M., 1996, A&A, 300, 707
 Torii et al. (2011) Torii S., Yamada S., Makishima K, et al., 2011, PASJ, 63, S771
 Uttley et al. (2011) Uttley P., Wilkinson T., Cassatella P., Wilms E., Pottschmidt K., Hanke M. & Böck M., 2011, MNRAS, 414, L60
 Uttley et al. (2014) Uttley P., Cackett E.M., Fabian A.C., Kara E. & Wilkins D.R., 2014, Astron. Astrophys. Rev., 22, 72
 Veledina (2016) Veledina A., 2016, ApJ, 832, 2
 Wijnands & van der Klis (1999) Wijnands R., van der Klis M., 1999, ApJ, 522 (2), 965
 Wilkinson & Uttley (2009) Wilkinson T. & Uttley P., 2009, MNRAS, 397, 666
 Yamada et al. (2013) Yamada S., Makishima K., Done C., Torii S., Noda H. & Sakurai S., 2013, PASJ, 65, 80
 Zdziarski, Johnson & Magdziarz (1996) Zdziarski A.A., Johnson W.N. & Magdziarz P., 1996, MNRAS, 283, 193
Appendix A Generalised Lags
IK13 show that the PSD form in Eq. 14 can be adapted to yield an analytic form for the cross spectrum between a low and high band, . For with a mean of and rmsnormalisation, ignoring smoothing, this form becomes
(20)  
where
(21) 
The real and imaginary parts respectively are then
(22)  
(23)  
From these components the time lag is computed in generality as
(24) 
We show an example of this analytic lag in Fig. 12. Inconsistencies with the simulation output arise from the finite number of simulation realisations and the spatial resolution of the simulations.
Appendix B Analytic Case With Suppression of Soft Variability
In the case of damping of the amplitude of the variability on propagation from the soft to the hard region, Eq. 14 modifies to
(25)  
where
(26) 
(27)  
and
(28)  
respectively.