A FOURIER-TRANSFORMED BREMSSTRAHLUNG FLASH MODEL FOR THE PRODUCTION OF X-RAY TIME LAGS IN ACCRETING BLACK HOLE SOURCES

The high-energy emission from many accretion-powered X-ray binaries is characterized by strong variability on timescales extending from milliseconds to weeks or longer. The softest X-ray emission (including a very soft thermal component) is observed during phases of rapid accretion, and a hard X-ray component (along with associated strong radio emission) is observed during periods of relatively slow accretion. The quiescent spectrum in the hard state is characterized by a power-law X-ray continuum, which is most likely due to the thermal Comptonization of soft X-rays in a hot (∼10^8–9 K) corona (e.g., Shapiroet al. 1976; Sunyaev & Titarchuk 1980).

The X-ray variability observed on millisecond timescales allows us to glimpse the inner workings of the accretion disk, but at these short timescales, the signal-to-noise ratio is too low to study the variability using snapshots of the spectrum. A useful alternative is provided by the study of X-ray time lags, which are computed by combining the Fourier-transformed data streams in two energy channels. Fourier-based analysis is extremely efficient because all of the data in the observational time window is utilized, significantly reducing the uncertainty in the derived parameters.

The physical origin of the observed X-ray time lags is not yet completely understood. The time lags may be the result of the upscattering of soft seed photons by hot electrons, known as "Compton reverberation" (Payne 1980), which naturally produces hard time lags because the higher-energy photons spend more time upscattering in the plasma before escaping. This concept was explored by van der Klis et al. (1987) and Miyamoto et al. (1988), who concluded that the observed decrease in the time lag with increasing Fourier frequency in sources such as Cyg X-1 was inconsistent with the Comptonization model. However, the idea was reexamined later by Nowak et al. (1999a) and Hua et al. (1999, hereafter HKC), who used more sophisticated simulations that produced better agreement with the data, but the size of the scattering cloud and the implied energy budget remained problematic.

In this Letter, we reconsider the possible role of thermal Comptonization in the formation of the X-ray time lags observed from black hole candidates. We develop a rigorous method for the computation of X-ray time lags based on the exact solution for the Fourier transform of the Green's function describing the radiation signal emerging from a spherical Compton scattering corona. By exploiting the linearity of the problem, we are able to explore a wide variety of injection scenarios, providing a useful alternative to Monte Carlo simulations, which are often either unreliable or noisy at high Fourier frequencies (and short timescales) due to sampling errors. We confirm that the Comptonization of monochromatic seed photons injected into a homogeneous scattering cloud is unable to account for the observed dependence of the X-ray time lags on the Fourier frequency. However, we find that the reprocessing of broadband (bremsstrahlung) injection reproduces the observed frequency dependence within the constraints of an acceptable cloud size and energy budget.

In the complex cross-spectrum (CCS) method first introduced by van der Klis et al. (1987), intensity variations in two energy channels _hard and _soft are used to extract phase shifts and associated time lags. The CCS method is implemented by Fourier transforming the light curves in the hard and soft energy channels, h(t) and s(t), respectively, to obtain the corresponding Fourier-transformed light curves H(ν_f) and S(ν_f), where ν_f is the Fourier frequency. The two transformed light curves are used to construct the complex cross spectrum, defined by C = S*(ν_f)H(ν_f), and the resulting time lag is given by

$$\begin{equation} \delta t=\frac{\arg (C)}{2\pi \nu _f}=\frac{\arg \left(S^{*}H\right)}{2 \pi \nu _f} , \end{equation} \tag{ 1 }$$

where $\arg (C)$ gives the phase lag (Nowak et al. 1999a). One can easily demonstrate that if a perfect time delay between the two channels is introduced, such that h(t) = s(t − Δt), where Δt > 0, then the time lag derived using the CCS method described above is δt = Δt, as expected. It is important to note that time lags are only created during transients, and are never produced during the steady-state reprocessing of radiation in a scattering cloud with constant properties.

Our focus here is on the determination of the effect of electron scattering on the signal measured by a distant observer during a transient in which a sudden flash of seed radiation is injected into a hot corona overlying an accretion disk. In this first exact treatment of the process, we restrict our attention to the case of a spherical, homogeneous, and isothermal cloud. Once the solution for the Fourier transform of the Green's function is in hand, we can compute the corresponding theoretical prediction for the X-ray time lags using Equation (1).

The reprocessing of N₀ seed photons injected simultaneously at radius r₀ and time t₀ with energy ₀ is governed by the time-dependent transport equation (e.g., Becker 1992, 2003)

$$\begin{eqnarray} {\partial f_{{\rm G}}\over \partial t}\! &=& {n_e \sigma _{\rm T}c \over m_e c^2} {1 \over \epsilon ^2}{\partial \over \partial \epsilon }\left [\epsilon ^4 \left (f_{{\rm G}}+ kT_e{\partial f_{{\rm G}}\over \partial \epsilon }\right)\right ] + {1 \over r^2} {\partial \over \partial r}\left (\kappa _0 r^2 {\partial f_{{\rm G}}\over \partial r}\right) \nonumber \\ &&+{N_0 \, \delta (\epsilon -\epsilon _0)\delta (r-r_0)\delta (t-t_0) \over 4\pi r_0^2\epsilon _0^2} , \end{eqnarray} \tag{ 2 }$$

where f_G(, r) is the Green's function describing the distribution of photons with energy at radius r inside the cloud, and T_e and n_e denote the electron temperature and number density, respectively. The terms on the right-hand side of Equation (2) represent thermal Comptonization, spatial diffusion, and photon sources, respectively, and κ₀ = c/(3n_eσ_T) denotes the spatial diffusion coefficient. The Green's function is related to the total photon number density, n_r, via

$$\begin{equation} n_r(r) = \int _0^\infty \epsilon ^2 \, f_{{\rm G}}(\epsilon,r) \, d\epsilon . \end{equation} \tag{ 3 }$$

In the case of a homogeneous, isothermal scattering corona considered here, κ₀ and T_e are both constants, and it is convenient to work in terms of the dimensionless energy x ≡ /kT_e, the dimensionless temperature Θ ≡ kT_e/m_ec², the dimensionless time p ≡ t(c²/3κ₀), and the scattering optical depth τ ≡ n_eσ_Tr = (c/3κ₀)r, in which case the transport Equation (2) can be rewritten as

$$\begin{eqnarray} {\partial f_{{\rm G}}\over \partial p} &=& {1 \over 3\tau ^2} {\partial \over \partial \tau } \left (\tau ^2 {\partial f_{{\rm G}}\over \partial \tau }\right) + {\Theta \over x^2}{\partial \over \partial x} \left [x^4\left (f_{{\rm G}}+ {\partial f_{{\rm G}}\over \partial x}\right)\right ] \nonumber \\ &&+ {N_0\,\delta (x-x_0)\delta (\tau -\tau _0)\delta (p-p_0) \over 4\pi \tau _0^2 x_0^2 (m_e c^2)^3 \Theta ^3 \ell _0^3} , \end{eqnarray} \tag{ 4 }$$

where τ₀ is the injection optical depth, x₀ is the dimensionless injection energy, p₀ is the dimensionless injection time, and ℓ₀ = 3κ₀/c denotes the (constant) electron scattering mean free path in the corona.

We define the Fourier transform pair (f_G, F_G) using

$$\begin{eqnarray} F_{{\rm G}}(x,\tau,\omega) &\equiv & \int _{-\infty }^\infty e^{i \omega p} f_{{\rm G}}(x,\tau,p) \, dp , \nonumber \\ f_{{\rm G}}(x,\tau,p) &\equiv & {1 \over 2 \pi } \int _{-\infty }^\infty e^{-i \omega p} F_{{\rm G}}(x,\tau,\omega) \, d\omega , \end{eqnarray} \tag{ 5 }$$

where ω = ν_f(2πℓ₀/c) is the dimensionless Fourier frequency. We can operate on Equation (4) with $\int _{-\infty }^\infty e^{i\omega p}dp$ to obtain

$$\begin{eqnarray} {-} i \omega F_{{\rm G}}&=& {1 \over 3\tau ^2} {\partial \over \partial \tau } \left (\tau ^2 {\partial F_{{\rm G}}\over \partial \tau }\right) + {\Theta \over x^2}{\partial \over \partial x} \left [x^4\left (F_{{\rm G}}+ {\partial F_{{\rm G}}\over \partial x}\right)\right ] \nonumber \\ &&+ {N_0\,\delta (x-x_0)\delta (\tau -\tau _0)e^{i \omega p_0} \over 4\pi \tau _0^2 x_0^2 (m_e c^2)^3 \Theta ^3 \ell _0^3} , \end{eqnarray} \tag{ 6 }$$

where i² = −1.

When x ≠ x₀, Equation (6) is separable using F_λ = G(λ, τ)H(λ, x) and this yields the differential equations

$$\begin{equation} {1 \over \tau ^2} {d \over d\tau }\left (\tau ^2 {dG \over d\tau }\right) + \lambda G = 0 , \end{equation} \tag{ 7 }$$

$$\begin{equation} {1 \over x^2} {d \over dx}\left [x^4\left (H + {dH \over dx}\right)\right ] - {s \over 3 \Theta } H = 0 , \end{equation} \tag{ 8 }$$

where λ is the separation constant and s ≡ λ − 3iω. The solutions for G and H satisfying suitable boundary conditions in energy and radius are

$$\begin{equation} G(\lambda,\tau)={\sin (\tau \sqrt{\lambda }) \over \tau } , \end{equation} \tag{ 9 }$$

and

$$\begin{equation} H(\lambda,x)=(x_0 x)^{-2} e^{-(x+x_0)/2} \, M_{2,\mu }(x_{\rm min}) \, W_{2,\mu }(x_{\rm max}) , \end{equation} \tag{ 10 }$$

where M_{2, μ} and W_{2, μ} denote Whittaker functions, the constant μ is given by

$$\begin{equation} \mu \equiv \left ({9 \over 4} + {\lambda \over 3 \Theta } - {i \omega \over \Theta } \right)^{1/2} , \end{equation} \tag{ 11 }$$

and we have made the definitions x_min ≡ min (x, x₀), x_max ≡ max (x, x₀).

The spherical scattering corona has a finite size, with outer radius r = R. At the surface of the cloud, the distribution function f must satisfy the free-streaming boundary condition −κ₀(∂f/∂r) = cf. This expression yields a constraint on the spatial separation function G, which can be written in terms of the optical depth τ as

$$\begin{equation} {1 \over 3}{dG \over d\tau } + G \left |\right. _{\tau =\tau _*} = 0 , \end{equation} \tag{ 12 }$$

where τ_* ≡ n_eσ_TR = R/ℓ₀ is the scattering optical thickness from the center to the edge of the corona. The roots of Equation (12) are the eigenvalues λ_n. Since the eigenvalue equation is transcendental, the eigenvalues must be determined using a numerical root-finding procedure. The eigenvalues λ_n are all real and positive, and the corresponding values of μ are computed using Equation (11) by setting λ = λ_n. The associated eigenfunctions G_n and H_n are defined by

$$\begin{equation} G_n(\tau) \equiv G(\lambda _n,\tau) , \quad H_n(x) \equiv H(\lambda _n,x) . \end{equation} \tag{ 13 }$$

Application of the Sturm–Liouville theorem shows that the spatial eigenfunctions G_n form a complete orthogonal set.

The exact solution for the Fourier transform F_G of the Green's function describing the time-dependent photon distribution can now be expressed using the series

$$\begin{equation} F_{{\rm G}}(x,\tau,\omega)=\sum _{n=0}^\infty a_n G_n(\tau) H_n(x) , \end{equation} \tag{ 14 }$$

where the expansion coefficients a_n can be calculated by exploiting the orthogonality of the spatial eigenfunctions G_n, combined with the derivative jump condition

$$\begin{equation} \Delta \left [{\partial F_{{\rm G}}\over \partial x}\right ] = {-N_0 \, \delta (\tau -\tau _0) e^{i \omega p_0} \over 4 \pi \tau _0^2 x_0^4 (m_e c^2)^3 \ell _0^3 \Theta ^4} , \end{equation} \tag{ 15 }$$

obtained by integrating the transport Equation (6) over a small energy range around x = x₀. Solving for the expansion coefficients and combining the result with Equation (14), we find that the exact solution for the Fourier transform of the Green's function of the time-dependent spectrum is given by

$$\begin{equation} F_{{\rm G}} (x, \tau)= {N_0 \, e^{i \omega p_0} \, e^{x_0} \over 4 \pi \ell _0^3 (m_e c^2)^3} \sum _{n=0}^\infty {\Gamma (\mu \!-\!3/2) \, G_n(\tau _0) G_n(\tau) H_n(x) \over \Theta ^4 \, \Gamma (1+2\mu) \, \mathcal {I}_n} , \nonumber \end{equation} \tag{ 16 }$$

where the quadratic normalization integrals of the spatial eigenfunctions are defined by $\mathcal {I}_n\equiv \int _0^{\tau _*} \tau ^2 G^2_n(\tau)d\tau$, and the values of μ are computed by substituting the eigenvalues λ_n into Equation (11).

In the integrated model considered here, both the steady-state (quiescent) and the transient X-ray spectral components are produced via thermal Comptonization in the corona, but the source of seed photons is different in the two cases. We assume here that the steady-state quiescent spectrum is the result of the upscattering in the corona of soft seed photons continually injected into the corona from the underlying cool disk (note however that reprocessed bremsstrahlung emission produced in the corona itself may also contribute significantly to the high-energy X-ray spectrum). On the other hand, we assume that the transient component (responsible for the time lags) is the result of a sudden flash of seed photons injected into the corona at a particular radius, due to some instability, which may produce either monochromatic or broadband emission.

We can tie down the fundamental cloud parameters Θ and τ_* by comparing the theoretical quiescent spectrum with the observational X-ray data. The quiescent spectrum, $f_{{\rm G}}^{\rm S}(\epsilon,r)$, is computed by solving the steady-state transport equation

$$\begin{eqnarray} 0 &=& {n_e \sigma _{\rm T}c \over m_e c^2} {1 \over \epsilon ^2}{\partial \over \partial \epsilon }\left [\epsilon ^4 \left (f_{{\rm G}}^{\rm S}+ kT_e{\partial f_{{\rm G}}^{\rm S}\over \partial \epsilon }\right)\right ] \nonumber \\ &&+ {1 \over r^2} {\partial \over \partial r}\left (\kappa _0 r^2 {\partial f_{{\rm G}}^{\rm S}\over \partial r}\right) + {\dot{N}_0 \, \delta (\epsilon -\epsilon _0) \over (4/3) \pi R^3 \epsilon _0^2} , \end{eqnarray} \tag{ 17 }$$

where $\dot{N}_0$ denotes the rate at which photons with energy ₀ are injected uniformly throughout the cloud of radius R. Adopting the same dimensionless variables x and τ used in Equation (4) now yields

$$\begin{eqnarray} 0 &=& {1 \over 3\tau ^2} {\partial \over \partial \tau } \left (\tau ^2 {\partial f_{{\rm G}}^{\rm S}\over \partial \tau }\right) + {\Theta \over x^2}{\partial \over \partial x} \left [x^4\left (f_{{\rm G}}^{\rm S}+ {\partial f_{{\rm G}}^{\rm S}\over \partial x}\right)\right ] \nonumber \\ &&+ {\dot{N}_0\,\ell _0\,\delta (x-x_0) \over (4/3) \pi R^3 c (m_e c^2)^3 \Theta ^3 x_0^2} . \end{eqnarray} \tag{ 18 }$$

In analogy with Equation (6), we can separate Equation (18) for x ≠ x₀ in terms of the functions f_λ = G(λ, τ)K(λ, x), where G and K satisfy the differential equations

$$\begin{equation} {1 \over \tau ^2} {d \over d\tau }\left (\tau ^2 {dG \over d\tau }\right) + \lambda G = 0 , \end{equation} \tag{ 19 }$$

$$\begin{equation} {1 \over x^2} {d \over dx}\left [x^4\left (K + {dK \over dx}\right)\right ] - {\lambda \over 3 \Theta } K = 0 . \end{equation} \tag{ 20 }$$

The solutions for G and K are given by

$$\begin{equation} G(\lambda,\tau)={\sin (\tau \sqrt{\lambda }) \over \tau } , \end{equation} \tag{ 21 }$$

and

$$\begin{equation} K(\lambda,x)=(x_0 x)^{-2} e^{-(x+x_0)/2} \, M_{2,\sigma }(x_{\rm min}) \, W_{2,\sigma }(x_{\rm max}) , \end{equation} \tag{ 22 }$$

respectively, where

$$\begin{equation} \sigma \equiv \left ({9 \over 4} + {\lambda \over 3 \Theta }\right)^{1/2} . \end{equation} \tag{ 23 }$$

Since Equations (7) and (19) are identical, it follows that the eigenvalues λ_n are exactly the same ones obtained in the solution for the Fourier transform F_G in Section 4. The corresponding eigenfunctions G_n and K_n are given by

$$\begin{equation} G_n(\tau) \equiv G(\lambda _n,\tau), \, \ \ \ K_n(x) \equiv K(\lambda _n,x) . \end{equation} \tag{ 24 }$$

The exact solution for the steady-state photon distribution function $f_{{\rm G}}^{\rm S}$ can now be expressed using the series

$$\begin{equation} f_{{\rm G}}^{\rm S}(x,\tau)=\sum _{n=0}^\infty b_n G_n(\tau) K_n(x) , \end{equation} \tag{ 25 }$$

where the expansion coefficients b_n are calculated using the same approach used to obtain the coefficients a_n in the expansion for F_G. In this case, integration of Equation (18) yields the derivative jump condition

$$\begin{equation} \Delta \left [{\partial f_{{\rm G}}^{\rm S}\over \partial x}\right ] = {- \dot{N}_0 \, \ell _0 \over (4/3) \pi R^3 c (m_e c^2)^3 \Theta ^4 x_0^4} . \end{equation} \tag{ 26 }$$

After some algebra, we find that the exact solution for the steady-state (quiescent) Comptonized spectrum is given by

$$\begin{equation} f_{{\rm G}}^{\rm S}(x,\tau) = {9 \dot{N}_0 \, e^{x_0} \over 4 \pi R^2} \sum _{n=0}^\infty {\Gamma (\sigma \!-\!3/2) \, \sin (\tau _*\sqrt{\lambda _n}) G_n(\tau) K_n(x) \over \Theta ^4 \, c \, (m_e c^2)^3 \lambda _n \ \Gamma (1+2\sigma) \, \mathcal {I}_n} , \end{equation} \tag{ 27 }$$

where the values of σ are computed by substituting the eigenvalues λ_n into Equation (23). The photon number flux measured at the detector, ${\cal F}_\epsilon$, can be computed from $f_{{\rm G}}^{\rm S}(x,\tau)$ using

$$\begin{equation} {\cal F}_\epsilon (\epsilon) = \left(R \over D\right)^2 c \ \epsilon ^2 f_{{\rm G}}^{\rm S}\left({\epsilon \over k T_e},\tau _*\right) , \end{equation} \tag{ 28 }$$

where D is the distance to the source. The emergent spectrum will be compared with the quiescent data for Cyg X-1 in Section 6.

The integrated model we have developed can be applied to any active galactic nucleus or accreting galactic black hole candidate using a two-step approach. First, the cloud temperature Θ and optical thickness τ_* are determined by comparing the photon flux computed using Equation (28) with the observed quiescent X-ray spectrum. Second, the Fourier transform of the emergent signal is computed, and the results are substituted into the CCS method using Equation (1) to calculate the associated time lags.

As an example, we use the integrated model to interpret the X-ray spectral and timing observations of Cyg X-1. In Figure 1(a), the steady-state spectrum (28) is compared with the quiescent X-ray spectrum observed in the low hard state of Cyg X-1 by Cadolle Bel et al. (2006) assuming D = 2.4 kpc. The solid curve in Figure 1(a) represents the steady-state spectrum resulting from the thermal Comptonization of monochromatic seed photons continually injected throughout the entire corona with energy ₀ = 0.1 keV, which approximates a T = 10⁶ K blackbody source. The temperature and optical thickness of the corona are Θ = 0.12 and τ_* = 2.5, respectively.

Zoom In Zoom Out Reset image size

In computing the associated time lags, we consider both monochromatic injection and bremsstrahlung injection. The Fourier transform in the bremsstrahlung case is computed using the integral convolution

$$\begin{equation} F_{\rm brem}(x) = \int _{x_{\rm abs}}^\infty F_{{\rm G}}Q N_0^{-1} dx_0 , \end{equation} \tag{ 29 }$$

where the bremsstrahlung source term is given by $Q \propto e^{-x_0}/x_0$ and x_abs = _abs/(kT_e) is the low-energy cutoff due to self-absorption in the source (Rybicki & Lightman 1979). The integral convolution in Equation (29) can be carried out analytically but we do not include the result here for brevity. The Fourier transform of the signal is computed using either Equation (16) or (29), for monochromatic or bremsstrahlung injection, respectively, and then the time lags are calculated using Equation (1).

In Figure 1(b), the theoretical time lags are plotted as a function of the Fourier frequency ν_f and compared with the RXTE data obtained in the hard state of Cyg X-1 by Nowak et al. (1999a). The time lags are computed using the same values Θ = 0.12 and τ_* = 2.5 used to calculate the quiescent spectrum in Figure 1(a). The red dashed curve corresponds to the injection of monochromatic seed photons with energy ₀ = 0.1 keV at optical depth τ₀ = 2.5 into a cloud with radius R = 10, 000 km. The blue dot-dashed curve denotes the lags resulting from broadband (bremsstrahlung) injection at optical depth τ₀ = 2.5 with _abs = 1.6 keV into a cloud with radius R = 40, 000 km. The black solid curve represents the lags resulting from bremsstrahlung injection at optical depth τ₀ = 2.5 with _abs = 1.6 keV into a cloud with radius R = 30, 000 km. The bremsstrahlung time lag results fit the observational data fairly well, which would imply that the X-ray time lags are produced at much larger radii than the broad iron lines.

Based on the results plotted in Figure 1(b), we conclude that thermal Comptonization of monochromatic radiation is unable to reproduce the observed time lag profiles, in agreement with the conclusions reached by Miyamoto et al. (1988), Nowak et al. (1999b), and HKC. On the other hand, the two curves in Figure 1(b) corresponding to the reprocessing of bremsstrahlung seed radiation are a near match to the observed profiles at high frequencies, although the time lags plateau at low frequencies due to the homogeneous density distribution assumed here. HKC were able to remedy the problem at low frequencies by focusing on an inhomogeneous corona with electron density n_e∝1/r, although this approach leads to other problems as discussed below.

Compton upscattering in a hot corona produces time lags that increase logarithmically with increasing hard channel energy, and the approximately logarithmic dependence seen in the data provided part of the observational motivation for the simulations carried out by HKC. Our results confirm the logarithmic behavior, as expected, but there are some interesting differences between our results and those obtained by HKC. Our model requires a relatively high electron temperature, T_e ∼ 10⁸ K, at large radii, r ∼ 10³ GM/c², in order to explain the observed time lags, in agreement with other Compton scattering models (HKC; Poutanen 2001). This temperature is higher than some disk models predict at that distance, although we note that You et al. (2012) have developed a fully relativistic two-temperature disk-corona model that predicts electron temperatures T_e ∼ 10⁸ K on size scales comparable to the radius of the corona in our model. The heating problem is more severe in the model of HKC, which requires a hot corona extending out to ∼10^4–5GM/c², and it is unclear whether the dissipation mechanism proposed by You et al. (2012) can provide the required energy at that distance.

The qualitatively different behavior observed in the bremsstrahlung injection scenario stems from the fact that in the bremsstrahlung case, photons with energies inside both energy channel windows already exist in the injected spectrum. Hence, when the injection occurs close to the edge of the cloud, "prompt" photons inside the energy channel windows are able to escape immediately. The prompt escape phase represents the fastest timescale during the transient, and consequently this is the process that contributes to the high-frequency part of the Fourier transform. The resulting time lag diminishes with increasing Fourier frequency because the high- and low-energy prompt photons escape at the same rate, with no upscattering required, so there is no relative delay between the two channel energies at the highest Fourier frequencies.

Conversely, at low frequencies, the process becomes dominated by the longest timescale phenomenon, which is the exponential decay of the photon number density due to photons that remain in the plasma for a long time. This phase is dominated by the upscattering of very soft bremsstrahlung seed photons, which take longer to upscatter to the higher-energy channel. Hence the low-frequency behavior of the bremsstrahlung time lag profile is similar to what is observed in the case of monochromatic injection.

We have developed a rigorous analytical model with very few free parameters that is capable of simultaneously describing the quiescent X-ray spectrum from Cyg X-1, as well as the transient time lags, based on thermal Comptonization in a homogeneous corona. In future work, we plan to extend the analytical model to treat an inhomogeneous corona, with electron density n_e∝1/r. The homogeneous case treated here underestimates the time lags at low frequencies, and we expect the inhomogeneous density profile to result in a better fit at low fourier frequencies because the distribution of scattering times will be more accurately computed.

The authors are grateful to the anonymous referee for a number of insightful comments and suggestions.