Extending the theory of propagating fluctuations: the first fully relativistic treatment and analytical Fourier–Green’s functions

ABSTRACT

The aperiodic variability ubiquitously observed from accreting black hole X-ray binary systems is generally analysed within the framework of the so-called ‘theory of propagating fluctuations’. In this paper we derive the Fourier transforms of the Green’s function solutions of the thin disc equations. These solutions suffice to describe all possible solutions through standard convolution techniques. Solutions are found for both Newtonian discs and general relativistic solutions with a vanishing ISCO stress. We use this new relativistic theory to highlight the Kerr black hole spin dependence of a number of observable variability properties of black hole discs. The phase lags, coherence, and power density spectra of Kerr discs are shown to be strong functions of black hole spin. Observations of the aperiodic variability of black hole accretion sources may now, at least in principle, offer a new avenue to directly constrain black hole spins.

accretion, accretion discs, black hole physics, X-rays: binaries

1 INTRODUCTION

Accreting black hole systems generally exhibits pronounced temporal variability, a result of their fundamentally turbulent nature (Balbus & Hawley 1991). The power spectral densities of the luminosity emergent from accretion discs reveal fluctuations whose root mean square variation on short time-scales is found to be linearly proportional to its mean evolving over longer time-scales, a property of the log-normal distribution (Uttley & McHardy 2001; Uttley, McHardy & Vaughan 2005).

Different accreting sources display a remarkable similarity in their variability properties, despite the vast range of both length and time-scales involved, across a broad population of sources. In any individual source variability is observed over many temporal orders of magnitude: from time-scales as rapid as the local dynamical time-scale of the innermost disc edge up to many orders of magnitude longer than any physical process which is involved in the direct production of the disc luminosity. The length-scales spanned by different sources which show similar variability structure is vast: ranging from the compact discs in Galactic X-ray binaries (e.g. Gleissner et al. 2004), to the very large discs in active galactic nuclei (e.g. Vaughan et al. 2011).

This behaviour has a robust observational grounding, and has been confirmed with observations at various different frequencies in black hole accretion disc sources, including both AGN in X-ray (Gaskell 2004; Vaughan et al. 2011), AGN in optical (Lyutyj & Oknyanskij 1987), and X-ray binaries at both X-ray (Gleissner et al. 2004) and optical (Gandhi 2009) frequencies. We also note that non black hole disc sources, e.g. cataclysmic variables (Scaringi et al. 2012), and young stellar objects (Scaringi et al. 2015) also show the same variability structure.

The typical variability structure of an accreting system is the following. The observed power density spectrum has a broad (aperiodic) component which is generally well described by a twice-broken power-law. In the case of some black hole binaries, these power density spectra display a narrow feature peaking in the range ∼0.1–10 Hz’s, which are known as (type C) quasi-periodic oscillations (QPOs). We will discuss only the aperiodic variability of black hole sources in this paper, and not QPOs, which will have different physical origins. Observed light curves of the same sources taken in different energy bands are found to correlate with each other, but ‘harder’ (higher energy) X-ray variability usually lags ‘softer’ (lower energy) X-ray variability (Priedhorsky et al. 1979; Nolan et al. 1981). The magnitude of the hard-soft time-lag depends on Fourier frequency, but is typically found to be of the order of 1 per cent of the variability time-scale. Some sources on the other hand show negative time lags, with the soft-band variability lagging the hard-band variability (McHardy et al. 2007; Emmanoulopoulos, McHardy & Papadakis 2011; De Marco et al. 2013). This type of behaviour has been detected in both supermassive and stellar mass black holes. Finally, the variability in emission from different energy bands is found to be coherent at low frequencies, but becomes increasingly incoherent at higher frequencies (Nowak et al. 1999).

This observed aperiodic variability structure is rather naturally explained by the so-called ‘theory of propagating fluctuations’ (first put forward by Lyubarskii 1997), in which fluctuations in the disc’s alpha parameter (or equivalently surface density) are excited at all radii in the accretion flow. These fluctuations then evolve throughout the disc and are observed as stochastic variability in the source’s light curves. As fluctuations are excited at every disc radii different time-scales are injected into the accretion flow corresponding to the local evolution time-scales of the individual fluctuations at each distance from the central object (Lyubarskii 1997; Churazov, Gilfanov & Revnivtsev 2001; Ingram 2016). The propagating fluctuations model naturally explains, for example, the lagging of ‘hard’ (higher energy) X-ray variability behind ‘softer’ (lower energy) bands (Kotov, Churazov & Gilfanov 2001; Ingram & van der Klis 2013). Excitations sourced in the outer, cooler, disc regions first produce softer flares before subsequently propagating inwards, sourcing harder flares in the hotter innermost disc regions.

In this model, the power spectral shape of the broad-band aperiodic noise depends on both the noise generating process and the response of the accretion flow. The magnetorotational instability (MRI; Hawley & Balbus 1991; Balbus & Hawley 1998) is the underlying noise generator, which produces variability everywhere in the disc due to the interactions between magnetic fields and a differentially rotating gas. The response of an accretion flow to intrinsic fluctuations is governed by a diffusion equation (Lynden-Bell & Pringle 1974). Solving this diffusion equation for a δ-function perturbation (i.e. calculating the Green’s function) is an important step in calculating the resulting power spectrum of the mass accretion rate at radii which differ from where the noise originated.

On the technical level, simplified Green’s functions were used for modelling the propagating fluctuations by Lyubarskii (1997) and Kotov et al. (2001), where the fluctuations were assumed to evolve in an additive manner. The more physical model of multiplicative fluctuations in the accretion flow were first considered by Ingram & van der Klis (2013) and Ingram & Done (2011), but in each of these cases only inward propagation was considered. However, fluctuations in (for example) an accretion flow’s surface density do not simply evolve inwards towards the central object; they must conserve the flow’s total angular momentum. This means that some material must propagate outwards in the flow, soaking up the angular momentum of the material propagating inwards, and fluctuations in the inner disc must therefore effect the properties of the outer disc.

This important conceptual point was first analysed in detail by Mushtukov, Ingram & van der Klis (2018), who developed a framework for analysing the combined effects of inward and outward propagating fluctuations in an accretion flow. Mushtukov et al. (2018) used this new framework to show that this additional outward propagation has potentially important observational effects, including that propagating fluctuations can give rise not only to hard time lags as previously shown, but may also produce negative lags (softer bands lagging harder bands) at high frequencies, a routinely observed effect which had previously been attributed to reprocessing.

The Mushtukov et al. (2018) framework involves computing the Fourier transform of the Green’s function solutions of the classical thin disc equations (which we shall henceforth call the Fourier–Green’s functions). This Fourier–Green’s integral has a long history in the literature, having been first written down in the original Lyubarskii (1997) paper, it has since reappeared in various different forms in a number of subsequent works, but has only ever been solved numerically (see also Mushtukov et al. 2019).

In this work, we present two important advances in the theory of propagating fluctuations. The first is deriving the exact solution of the Fourier–Green’s integral, which turns out to be surprisingly simple in its final form. With the exact analytical solutions of the Fourier integral now at hand, various properties of these solutions may be derived. In Section 2, for example, we demonstrate that high frequency variability in the mass accretion rate is suppressed as exp (−Δf^1/2), where Δ(x, x′) is a function of the magnitude of the difference between the two disc locations x and x′. The high and low frequency asymptotic behaviour of the power spectrum resulting from mass accretion rate variability is also determined, and related to the intrinsic variability in the disc surface density/alpha parameter. On a practical level, knowledge of these exact solutions will rapidly speed up, and improve the accuracy of, the process of fitting analytical models of accretion variability to observational data; the numerical cost of Fourier transforming thin disc Green’s functions had previously been substantial.

The second, and potentially more important, development is in presenting the first analysis of the Fourier–Green’s solutions of the general relativistic thin disc equation. In this paper, we solve the Fourier–Green’s integral for a general Kerr metric, under the assumption that the dynamical disc stress vanishes at the ISCO. These solutions depend implicitly on the central black hole’s spin through their dependence on the space–time’s ISCO radius, and for the first time the effects of the black hole’s spin on the observed variability structure of an accretion flow can be examined. In the later sections of this paper we demonstrate that a number of observable variability properties of black hole discs are relatively strong functions of black hole spin, and that observations of the aperiodic variability of black hole accretion discs may now, at least in principle, offer a new avenue to directly determine black hole spins.

The layout of this paper is the following. In Section 2, we derive and analyse the formal solutions of the Fourier–Green’s integral. In Section 3, we specialize the analysis specifically to the solutions of the Newtonian disc equation, while in Section 4 we discuss the corresponding relativistic solutions. In Section 5, we introduce the Mushtukov et al. (2018) framework for relating these results to directly observable quantities, before showing in Section 6 that the black hole spin imprints strong signals on to the observed aperiodic variability structure of observed light curves. We conclude in Section 7, with some technical results presented in Appendices.

The reader interested in the application of this analysis for observational constraints on black hole spins may wish to skip directly to Section 5.

2 Thin disc Green’s functions in the Fourier domain: general solution and properties

2.1 The disc evolution equations

The evolution of the disc surface density is described by a diffusion-type equations which fundamentally arise from the turbulent transportation of angular momentum within the disc. The evolution equation for a Newtonian theory of gravity was first analysed in detail by Lynden-Bell & Pringle (1974), while the relativistic equation was derived in Balbus (2017). We will introduce the general relativistic form of the evolution equation at this point, as it is simpler to take the Newtonian limit (r ≫ GM/c²) than the alternative.

The coordinates used to describe the relativistic thin disc equation are the cylindrical Boyer–Lindquist representation of the Kerr metric: r (radial), ϕ (azimuthal), and z (vertical). The governing equation describes the evolution of the azimuthally averaged, height-integrated disc surface density Σ(r, t). The contravariant four velocity of the disc fluid is U^μ; the covariant counterpart is U_μ. The specific angular momentum corresponds to U_ϕ, a covariant quantity. There is an anomalous stress tensor present, |${W^r_\phi }$|⁠, due to low level disc turbulence, which is a measure of the correlation between the fluctuations in U^r and U_ϕ (Eardley & Lightman 1975; Balbus 2017). This is, as the notation suggests, a mixed tensor. |${W^r_\phi }$| serves both to transport angular momentum as well as to extract the free-energy of the disc shear, which is then thermalised and radiated from the disc surface, both assumed to be local processes.

Under these assumptions the governing disc equation can be expressed in the following compact form

$$\begin{eqnarray} {\partial \zeta \over \partial t} = { {W^r_\phi }\over (U^0)^2}{\partial \ \over \partial r} \left({U^0\over U^{\prime }_\phi } \left[ {\partial \zeta \over \partial r}\right] \right). \end{eqnarray}$$

(1)

here the primed notation ′ denotes an ordinary derivative with respect to r, and

$$\begin{eqnarray} \zeta \equiv {r \Sigma {W^r_\phi }\over U^0} . \end{eqnarray}$$

(2)

The Newtonian limit corresponds to |$U^0 = 1, U_\phi ^{\prime } = \sqrt{GM/4r}$|⁠. Naturally, one must specify a functional form of the disc’s turbulent stress tensor |${W^r_\phi }$| to derive solutions of this equation, for the remainder of this paper we shall consider stress parametrizations of the form

$$\begin{eqnarray} {W^r_\phi }= w \left({r\over r_0}\right)^\mu , \end{eqnarray}$$

(3)

which is a popular and analytically tractable choice.

2.2 Solution of the Fourier integral

The Green’s function solutions of both the Newtonian (Lynden-Bell & Pringle 1974) and general relativistic (Mummery 2023) thin disc equations are of the general form

$$\begin{eqnarray} G(x, x_0, t) = {q(x) \over t} \exp \left({- g(x)^2 - g(x_0)^2 \over 4 t}\right) I_{\nu } \left({g(x)g(x_0)\over 2 t} \right) , \end{eqnarray}$$

(4)

where G(x, x₀, t) describes the evolution of the variable ζ for an initial delta-function spike located at t = 0, x = x₀. In this expression x is a radial coordinate, which is typically normalized by some characteristic scale, either the ISCO (in the relativistic case), or the initial radius r₀ of the spike. In this expression I_ν is the modified Bessel function of the first kind, of order ν. The index ν is related to the stress parametrization through

$$\begin{eqnarray} \nu = 1/(3-2\mu). \end{eqnarray}$$

(5)

For the particular case of the Newtonian Green’s functions these functions are particularly simple

$$\begin{eqnarray} q(x) \propto x^{1/4}, \quad g(x) \propto x^{1/4\nu } . \end{eqnarray}$$

(6)

The functions g(x), g(x₀) as defined here have dimensions of |$\sqrt{\rm time}$|⁠. Physically, the amplitude of these functions correspond to the (square root of the) time-scale with which a perturbation at x propagates to the inner disc edge. This introduces a natural scale into these expressions, which will be discussed further in the following sections. It will be useful to note however that this solution is in fact a Laplace-mode superposition, of the form (Lynden-Bell & Pringle 1974; Gradshteyn et al. 2007; Balbus 2017; Mummery 2023)

$$\begin{eqnarray} G(x, x_0, t) = \int _0^\infty q(x) J_\nu (\sqrt{s} g(x)) J_\nu (\sqrt{s} g(x_0)) \exp ({-st}) \, {\rm d}s . \end{eqnarray}$$

(7)

The function denoted J_ν is an ordinary Bessel function of the first kind. The Fourier transform of G(x, x₀, t), denoted |$\widetilde{G}(x, x_0, f)$| is defined by the complex integral

$$\begin{eqnarray} \widetilde{G}(x, x_0, f) \equiv \int _0^\infty G(x, x_0, t) \exp (-2\pi i f t) \, {\rm d} t , \end{eqnarray}$$

(8)

where we have used the fact that G(x, x₀, t < 0) = 0. When written in terms of the Laplace mode superposition, this integral becomes

$$\begin{eqnarray} \widetilde{G}(x, x_0, f) &=& q(x)\int _0^\infty \Bigg [ \int _0^\infty J_\nu (\sqrt{s} g(x)) J_\nu (\sqrt{s} g(x_0)) \\ && \exp ({-st}) \, {\rm d}s \Bigg ] \exp (-2\pi i f t) \, {\rm d} t . \end{eqnarray}$$

(9)

As both integrals converge, we can swap the order of integration

$$\begin{eqnarray} \widetilde{G}(x, x_0, f) &=& q(x)\int _0^\infty \Bigg [ \int _0^\infty \exp (-st - 2\pi i ft) \, {\rm d}t \Bigg ] \\ && J_\nu (\sqrt{s} g(x)) J_\nu (\sqrt{s} g(x_0)) \, {\rm d}s \end{eqnarray}$$

(10)

which is more easily solved. Performing the t integral leaves

$$\begin{eqnarray} \widetilde{G}(x, x_0, f) = q(x) \int _0^\infty { J_\nu (\sqrt{s} g(x)) J_\nu (\sqrt{s} g(x_0)) \over s + 2\pi i f} \, {\rm d}s . \end{eqnarray}$$

(11)

By making the substitution |$u = \sqrt{s}$|⁠, this integral becomes

$$\begin{eqnarray} \widetilde{G}(x, x_0, f) = 2 q(x) \int _0^\infty { u J_\nu (u g(x)) J_\nu (u g(x_0)) \over u^2 + \beta ^2} \, {\rm d}u , \end{eqnarray}$$

(12)

where

$$\begin{eqnarray} \beta \equiv (1 + i) \sqrt{\pi f}. \end{eqnarray}$$

(13)

When written in this form the solution of the integral is a standard result, which can be found in the text of Gradshteyn et al. (2007)

$$\begin{eqnarray} \widetilde{G}(x, x_0, f) = 2 q(x) \left\lbrace \begin{array}{@{}l@{\quad }l@{}}& I_\nu (\beta g(x)) K_\nu (\beta g(x_0)) , \quad x\lt x_0, \\ \\ & I_\nu (\beta g(x_0)) K_\nu (\beta g(x)) , \quad x\gt x_0 . \end{array}\right. \end{eqnarray}$$

(14)

In this expression K_ν is the modified Bessel function of the second kind. The Green’s function solutions for the mass accretion rate, denoted |$G_{\dot{M}}$|⁠, are of interest in understanding the variability properties of black hole discs, as it is often assumed that variability in the mass accretion rate is directly communicated into variability in the locally emitted flux (e.g. Lyubarskii 1997; Ingram & van der Klis 2013; Ingram & Done 2012; Mushtukov et al. 2018). The mass accretion rate Green’s function has the following form (we again write this generally so as to consider both the Newtonian and relativistic solutions simultaneously)

$$\begin{eqnarray} G_{\dot{M}} (x, x_0, t) = p(x) {\partial \over \partial x} G(x, x_0, t) . \end{eqnarray}$$

(15)

In the Newtonian limit the function p(x) is simple

$$\begin{eqnarray} p(x) \propto x^{1/2}. \end{eqnarray}$$

(16)

In the Fourier domain

$$\begin{eqnarray} \widetilde{G}_{\dot{M}} (x, x_0, f) &\equiv &\int _0^\infty G_{\dot{M}}(x, x_0, t) \exp (-2\pi i f t) \, {\rm d} t, \\ &=& \int _0^\infty p(x) {\partial \over \partial x} \left[G(x, x_0, t)\right] \exp (-2\pi i f t) \, {\rm d} t, \\ &=& p(x) {\partial \over \partial x}\widetilde{G}(x, x_0, f), \end{eqnarray}$$

(17)

where in going to the final line we have used the fact that the x derivative and t integral commute. We therefore have the general solution

$$\begin{eqnarray} {1\over 2} \widetilde{G}_{\dot{M}}= p(x) \left\lbrace \begin{array}{@{}l@{\quad }l@{}}& K_\nu (\beta g(x_0)) \, {\partial _x} \left[ q(x) I_\nu (\beta g(x))\right] , \quad x\lt x_0, \\ \\ & I_\nu (\beta g(x_0)) \, {\partial _x} \left[ q(x) K_\nu (\beta g(x)) \right] , \quad x\gt x_0 , \end{array}\right. \end{eqnarray}$$

(18)

where we use the notation ∂_x ≡ ∂/∂x. As we shall demonstrate in Section 3, this equation further simplifies for the particular case of the Newtonian disc equations. Before we specialize to either the Newtonian or relativistic regimes, we analyse the asymptotic properties of these Fourier–Green’s function.

2.3 Asymptotic properties

The asymptotic (f → ∞ and f → 0) properties of |$\widetilde{G}_{\dot{M}}$| can be determined from this general formula. These asymptotic limits should be understood as the limiting behaviour at Fourier frequencies which are significantly larger (or smaller) than the characteristic accretion frequency associated with radius x. The characteristic accretion frequency for a Newtonian solution with perturbation at r₀ has the following form (Lynden-Bell & Pringle 1974):

$$\begin{eqnarray} f_{0} = {1 \over t_{\rm acc}(r_0)} = {(3 - 2\mu)^2 \over 2} \sqrt{w^2 \over GM_{\rm BH} r_0^3} . \end{eqnarray}$$

(19)

(We derive the relativistic analogue of this expression in Appendix A). The amplitude of the functions g(x) are |$1/\sqrt{f_0}$|⁠.

2.3.1 High frequency Fourier modes

The large frequency limit (f ≫ f₀) of the two Bessel functions are the following:

$$\begin{eqnarray} \lim _{\beta \rightarrow \infty } K_\nu (\beta g) = \sqrt{\pi \over 2\beta g} \exp (-\beta g) , \end{eqnarray}$$

(20)

and

$$\begin{eqnarray} \lim _{\beta \rightarrow \infty } I_\nu (\beta g) = \sqrt{1 \over 2\pi \beta g} \exp (+\beta g) . \end{eqnarray}$$

(21)

We therefore have the following simple result, valid in both of the x < x₀ and x > x₀ regimes:

$$\begin{eqnarray} &&\widetilde{G}_{\dot{M}}(x, x_0, f\rightarrow \infty) \sim {p(x) q(x) g^{\prime }(x) \over \sqrt{g(x)g(x_0)}} \\ &&\exp (-\sqrt{\pi f} (1 + i) \left| g(x) - g(x_0) \right|) , \end{eqnarray}$$

(22)

where |z| denotes the absolute value of z, and ′ denotes a derivative with respect to x. We see here that high frequency Fourier modes are exponential suppressed, with a suppression scale which depends on the inverse of the absolute magnitude of the difference between any two disc radii.

Physically, this corresponds to an exponential suppression of the propagation of modes with frequencies higher than the local accretion frequency. It is interesting to note that this expression is symmetric in the sign of x − x₀. This is an important result: at high Fourier frequencies (relative to the local accretion frequency), outward propagation of material is equally as important as the inward propagation of material. This is not true at all Fourier frequencies, as we demonstrate in future sections.

The fact that the suppression of high-frequency modes is proportional to exp (− Δf^1/2) is unsurprising, and is a result of the disc angular-momentum diffusion processes. To see this explicitly, consider the simple 1D diffusion equation

$$\begin{eqnarray} {\partial \psi \over \partial t} = D {\partial ^2 \psi \over \partial x^2} . \end{eqnarray}$$

(23)

In terms of the Fourier modes (⁠|$\psi = \int _{-\infty }^{+\infty } \widetilde{\psi }\exp (2\pi i f t) {\rm d} f$|⁠), this equation reads

$$\begin{eqnarray} 2\pi i f \widetilde{\psi }= D {\partial ^2 \widetilde{\psi }\over \partial x^2} , \end{eqnarray}$$

(24)

with a solution that is well behaved at large x of

$$\begin{eqnarray} \widetilde{\psi }(x, f) = A \exp \left(- \sqrt{\pi f \over D} (1 + i) \, x\right) . \end{eqnarray}$$

(25)

It is clear to see that a high frequency exponential exp (− f^1/2) suppression is a generic property of diffusive systems.

2.3.2 Low frequency Fourier modes

The small frequency limit (f ≪ f₀) of |$\widetilde{G}_{\dot{M}}$| can also be understood from the asymptotic properties of the Bessel functions K_ν and I_ν. In addition, we can understand a priori the asymptotic properties in the small frequency limit, as f → 0 corresponds physically to the integral

$$\begin{eqnarray} \widetilde{G}_{\dot{M}}(x, x_0, f\rightarrow 0) \rightarrow \int _0^\infty G_{\dot{M}}(x, x_0, t) \, {\rm d}t . \end{eqnarray}$$

(26)

The integral over all times of the mass accretion rate at radius x, initiated by a perturbation at radius x₀, must equal (as all of the disc material is eventually accreted) to

$$\begin{eqnarray} \widetilde{G}_{\dot{M}}(x, x_0, f\rightarrow 0) \rightarrow \left\lbrace \begin{array}{@{}l@{\quad }l@{}}-M_d, \, \, \, \quad x \lt x_0, \\ \\ 0, \, \, \, \, \, \, \, \quad \quad x \gt x_0. \end{array}\right. \end{eqnarray}$$

(27)

To prove this more rigorously, consider the small frequency limits of the two Bessel functions:

$$\begin{eqnarray} \lim _{\beta \rightarrow 0} K_\nu (\beta g) = {\Gamma (\nu) \over 2} \left({\beta g \over 2}\!\right)^{-\nu } + {\Gamma (-\nu) \over 2} \left(\!{\beta g \over 2}\!\right)^{+\nu } + {\cal O}\left((\beta g)^{2-\nu }\right) \\ \end{eqnarray}$$

(28)

and

$$\begin{eqnarray} \lim _{\beta \rightarrow 0} I_\nu (\beta g) = {1 \over \Gamma (\nu +1)} \left({\beta g \over 2}\right)^{+\nu } + {\cal O}\left((\beta g)^{2+\nu }\right) \end{eqnarray}$$

(29)

where Γ(n) is the usual gamma function. It will transpire that we are required to keep the first two terms in the expansion of K_ν to compute the correct f → 0 limit of these solutions. For x < x₀ (inward propagating modes) we have

$$\begin{eqnarray} \!\left| \widetilde{G}_{\dot{M}}(x \lt x_0, f\rightarrow 0)\! \right| \sim {p(x) \over \nu } \!\left({1 \over g(x_0)}\!\right)^\nu {\!\partial \over \partial x\!} \!\left[ q(x) (g(x))^\nu \!\right] , \end{eqnarray}$$

(30)

i.e. a constant which is independent of f. With the proper normalization of q(x) and g(x) this constant will be equal to the disc mass. For outward propagating modes (x > x₀) we have

$$\begin{eqnarray} &&\left| \widetilde{G}_{\dot{M}}(x\gt x_0, f\rightarrow 0) \right| \sim {p(x) \over \nu } \left({ g(x_0)}\right)^\nu {\partial \over \partial x} \left[ { q(x) \over (g(x))^{\nu }} \right] \\ && + p(x) {\Gamma (-\nu) \over \Gamma (\nu +1)} \left({ g(x_0) \over 4 }\right)^\nu |\beta |^{2\nu } {\partial \over \partial x} \left[ { q(x) (g(x))^{\nu }} \right] + {\cal O}\left(|\beta |^2 \right), \\ \end{eqnarray}$$

(31)

which means

$$\begin{eqnarray} \left| \widetilde{G}_{\dot{M}}(x \gt x_0, f\rightarrow 0) \right| \sim \left\lbrace \begin{array}{@{}l@{\quad }l@{}}f^\nu , \quad { q(x) \propto (g(x))^{\nu }}, \quad \nu \lt 1, \\ \\ f^1 , \quad { q(x) \propto (g(x))^{\nu }}, \quad \nu \gt 1, \\ \\ f^{0} , \quad {\rm otherwise} . \end{array}\right. \end{eqnarray}$$

(32)

Unsurprisingly, for the particular case of a Newtonian disc model q(x) ∝ x^1/4 and g(x) ∝ x^1/4ν, and the amplitude of the outward propagating low frequency modes goes to zero as a power-law of frequency in the limit of f → 0.

At low Fourier frequencies (relative to the local accretion frequency), only the inwardly propagating material significantly contributes to the local variability of the accretion rate. As conventional approaches typically neglect outward propagation, they are likely to be most accurate at low frequencies.

2.4 Discontinuity of the mass accretion rate Fourier–Green’s function at x = x₀

As can be seen from the above analysis, the properties of the mass accretion rate Fourier–Green’s functions must be discontinuous at the location of the perturbation x = x₀. This can be proved rather generally by using the Wronskian of the modified Bessel functions (Abramowitz & Stegun 1965):

$$\begin{eqnarray} {\cal W}\left[I_\nu , K_\nu \right] \equiv I_\nu (z) {{\rm d} \over {\rm d } z} K_\nu (z) - K_\nu (z) {{\rm d} \over {\rm d } z} I_\nu (z) = -{1 \over z} , \end{eqnarray}$$

(33)

which implies a magnitude of the discontinuity of

$$\begin{eqnarray} &&\widetilde{G}_{\dot{M}}(x \rightarrow x_0, x_0 \gt x, f) - \widetilde{G}_{\dot{M}}(x \rightarrow x_0, x_0 \lt x, f) \\ && = p(x) q(x) {{\rm d} \over {\rm d} x} \ln g(x) . \end{eqnarray}$$

(34)

Physically this discontinuity results from the eventual accretion of all of the disc material through the inner disc edge. The preferred mass flow direction (inwards) fundamentally breaks the x < x₀, x > x₀ symmetry of the mass accretion rate Fourier–Green’s functions.

2.5 The mass accretion rate power spectrum

As derived by Mushtukov et al. (2018) the power spectrum of the mass accretion rate at a particular radius, denoted |$S_{\dot{m}}(x, f)$|⁠, is given by

$$\begin{eqnarray} S_{\dot{m}}(x, f) = \int _{x_{\rm in}}^{x_{\rm out}} \left({1 \over x^{\prime }}\right)^2 l(x^{\prime }) \left| \widetilde{G}_{\dot{M}}(x, x^{\prime }, f) \right|^2 S_\Sigma (x^{\prime }, f) \, {\rm d}x^{\prime }, \end{eqnarray}$$

(35)

where l(x′) is the radial scale over which the initial perturbations can be considered as coherent with one another, and |$S_\Sigma (x, f)$| is the power spectrum of the initial surface density perturbations at radius x. In general it is possible that fluctuations in certain regions of the disc might be coherent over a larger radial scale than in other regions, and so in general l(x′) is a function of radius. In this work, however, we shall assume that all disc radii have the same coherence length. Mushtukov showed numerically that at high and low frequencies the power |$P(f) = f S_{\dot{m}}(f)$| behaved like a power-law in frequency, assuming that the input surface density power spectrum |$S_\Sigma$| was Lorentzian:

$$\begin{eqnarray} S_\Sigma (x, f) = {2 p \over \pi } {f_{\rm br} \over f^2 + f_{\rm br}^2} , \end{eqnarray}$$

(36)

where p and f_br are the total power and break frequency (simply parameters of the model) at radius x, respectively.

As the focus of this work is on single epoch variability the parameter p, the amplitude of the locally driven power in the surface density perturbations, is treated as a simple constant. However, in a real physical system p has an important property: it scales quadratically with the local average (background) mass accretion rate |$\dot{M}^2$|⁠. This is a property of the propagating fluctuations model, which was first shown by Lyubarskii (1997), and was expanded upon by Mushtukov et al. (2018). This behaviour of p is important for understanding the variability properties of accreting sources at different stages of their evolution, as it ultimately results in the well known linear flux-rms relation (Uttley & McHardy 2001; Uttley et al. 2005; see Mushtukov et al. 2018 for a proof of this statement).

Finally, in addition to the leading order term presented here, there is formally an additional non-linear term present in equation (35), which arise from local fluctuations of the surface density superimposing on top of existing variability in the accretion flow. Mushtukov et al. (2018) found that the inclusion of these non-linear terms alters the quantitative properties of the local power spectrum at the |$1{{\ \rm per\ cent}}$| level, and do not alter the qualitative properties of the theory. The amplitude of deviation does however grow with the injected variability amplitude, and so this simplification should be kept in mind. These non-linear effects are beyond the scope of the present work, and will not be considered further.

In the following sections, we rigorously prove various numerical results presented in Mushtukov et al. (2018).

2.5.1 High frequency power spectrum

For f → ∞ the inward and outward propagating modes both have transfer functions which behave like

$$\begin{eqnarray} &&\left| \widetilde{G}_{\dot{M}}(x, x_0, f\rightarrow \infty) \right| \sim {p(x) q(x) g^{\prime }(x) \over \sqrt{g(x)g(x_0)}} \\ && \exp (-\sqrt{\pi f} \left| g(x) - g(x_0) \right|) , \end{eqnarray}$$

(37)

thus

$$\begin{eqnarray} && S_{\dot{m}}(x, f\rightarrow \infty) \sim \left({p(x) q(x) g^{\prime }(x) \over \sqrt{g(x)}}\right)^2 \int _{x_{\rm in}}^{x_{\rm out}} \left({1 \over x^{\prime }}\right)^2 {l(x^{\prime }) \over g(x^{\prime })} \\ && \exp (-2\sqrt{\pi f} \left| g(x) - g(x^{\prime }) \right|) \, S_\Sigma (x^{\prime }, f) \, {\rm d}x^{\prime }, \end{eqnarray}$$

(38)

which, provided that the high frequency fall-off of |$S_\Sigma (f)$| is weaker than exp (−Af^1/2), is an integral which can be performed by Laplace’s method. The leading order behaviour is simply

$$\begin{eqnarray} S_{\dot{m}}(x, f\rightarrow \infty) \sim f^{-1/2} S_\Sigma (x, f\rightarrow \infty) . \end{eqnarray}$$

(39)

This result highlights that at the highest frequencies, the observed variability is dominated by locally driven variability, with exponentially small contributions from distant disc regions. For the particular case of a Lorentzian power spectrum for Σ, we have

$$\begin{eqnarray} P(f\rightarrow \infty) \equiv f S_{\dot{m}}(x, f\rightarrow \infty) \sim f^{-3/2}. \end{eqnarray}$$

(40)

This exact behaviour was discovered numerically by Mushtukov et al. (2018).

2.5.2 Low frequency power spectrum

For f → 0, the inward propagating modes have transfer functions which become independent of frequency

$$\begin{eqnarray} \left| \widetilde{G}_{\dot{M}}(x\lt x_0, x_0, f\rightarrow 0) \right| \sim f^0 , \end{eqnarray}$$

(41)

thus

$$\begin{eqnarray} S_{\dot{m}}(x, f\rightarrow 0)& \sim &\int _{x}^{x_{\rm out}}F_1(x^{\prime }, x) \, S_\Sigma (x^{\prime }, f) \, {\rm d}x^{\prime } \\ && + \, f^\alpha \int _{x_{\rm in}}^{x} F_2(x^{\prime }, x) \, S_\Sigma (x^{\prime }, f) \, {\rm d}x^{\prime } , \end{eqnarray}$$

(42)

where F₁(x, x′) and F₂(x, x′) are independent of frequency, and α ≥ 0. The leading order behaviour of the mass accretion rate power spectrum is therefore given, in the low frequency limit, simply by that of the initial surface density perturbations at radii larger than x (the first integral)

$$\begin{eqnarray} S_{\dot{m}}(x, f\rightarrow 0) \sim \int _{x}^{x_{\rm out}}F_1(x^{\prime }, x) \, S_\Sigma (x^{\prime }, f\rightarrow 0) \, {\rm d}x^{\prime } . \end{eqnarray}$$

(43)

As expected from the earlier analysis (Section 2.3), only those perturbations initialized at larger radii x₀ > x which then propagate inwards contribute to the power spectrum in the low frequency limit.

For the particular case of a Lorentzian power spectrum for a, we have

$$\begin{eqnarray} P(f\rightarrow 0) \equiv f S_{\dot{m}}(f\rightarrow 0) \sim f^{1} , \end{eqnarray}$$

(44)

as found numerically by Mushtukov et al. (2018).

3 Newtonian Fourier–Green’s functions

The particular Green’s function solutions for a Newtonian theory of gravity have

$$\begin{eqnarray} q(x) \propto x^{1/4},\quad p(x) \propto x^{1/2}, \quad g(x) \propto x^{1/4\nu }. \end{eqnarray}$$

(45)

It transpires that the derivative in equation (18) simplifies greatly, a result of the identities

$$\begin{eqnarray} {{\rm d} \over {\rm d}z} I_l(z) = I_{l-1}(z) - {l \over z} I_l(z), \end{eqnarray}$$

(46)

and

$$\begin{eqnarray} {{\rm d} \over {\rm d}z} K_l(z) = - K_{l-1}(z) - {l \over z} K_l(z). \end{eqnarray}$$

(47)

We therefore have the remarkably simple result

$$\begin{eqnarray} \widetilde{G}_{\dot{M}} = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}+ A \beta x^{(1-\nu)/4\nu } K_\nu \left(\beta \epsilon x_0^{1/4\nu }\right) I_{\nu - 1}\left(\beta \epsilon x^{1/4\nu }\right)\!, \,\, x \lt x_0, \\ \\ - A \beta x^{(1-\nu)/4\nu } I_\nu \left(\beta \epsilon x_0^{1/4\nu }\right) K_{\nu - 1}\left(\beta \epsilon x^{1/4\nu }\right) \!, \,\, x \gt x_0 , \end{array}\right. \end{eqnarray}$$

(48)

where (units where GM = c = w = 1; see also equation 19)

$$\begin{eqnarray} A = {x_0^{1/4} \epsilon M_d }, \quad \epsilon = 2\nu \sqrt{2x_0^\mu } , \quad \mu = {3-1/\nu \over 2} , \end{eqnarray}$$

(49)

and we remind the reader

$$\begin{eqnarray} \beta \equiv (1 + i) \sqrt{\pi f}. \end{eqnarray}$$

(50)

Note that Newtonian gravity is entirely scale free, and x here should be thought of as the disc radius suitably normalized by some arbitrary radial scale in the problem of interest. For relativistic systems the solution is explicitly a function of r/r_g, where r_g = GM/c² is the gravitational radii of the black hole.

In Fig. 1, we display the amplitude of the Newtonian Fourier–Green’s function, as a function of Fourier frequency, for a number of different disc radii x (listed in caption), and x₀ = 10. It is clear to see that the properties of the Fourier–Green’s functions are as predicted by the asymptotic analysis of the previous section. At high Fourier frequencies there is an exponential suppression of the Fourier mode amplitude, with modes at disc radii x closer to x₀ having substantially higher amplitude at high frequencies when compared to disc radii further from x₀ (contrast the solid curves with the dot–dashed curves). At low Fourier frequencies the inward propagating modes (x < x₀) the Fourier–Green’s functions approach 1 (when normalized by the disc mass), while outward propagating mode (x > x₀) approach zero as a power-law in frequency.

The amplitude of the Newtonian Fourier–Green’s function of the mass accretion rate. Red curves are for inward propagating modes x < x0, and blue curves show outward propagating modes x > x0. For this figure, we take x0 = 10, and the inward propagating modes have x = 1 (dot–dashed), 2 (dotted), 4 (dashed), and 8 (solid), while the outward propagating modes have x = 11 (solid), 18 (dashed), 25 (dotted), and 30 (dot–dashed). The Fourier–Green’s functions are normalized so that the total accreted mass is equal to 1.

Figure 1.

The amplitude of the Newtonian Fourier–Green’s function of the mass accretion rate. Red curves are for inward propagating modes x < x₀, and blue curves show outward propagating modes x > x₀. For this figure, we take x₀ = 10, and the inward propagating modes have x = 1 (dot–dashed), 2 (dotted), 4 (dashed), and 8 (solid), while the outward propagating modes have x = 11 (solid), 18 (dashed), 25 (dotted), and 30 (dot–dashed). The Fourier–Green’s functions are normalized so that the total accreted mass is equal to 1.

Open in new tab Download slide

This low-frequency asymptotic behaviour is further demonstrated in Fig. 2, where we plot the amplitude of the Fourier–Green’s functions for different indices ν, for x = 20, x₀ = 10. At low frequencies each amplitude approaches a power-law in frequency, with power-law index given by equation (32; black dashed curves). In Fig. 3, we examine the effects of varying the stress index on the inward propagating Fourier–Green’s functions. While less severe than for the outward propagating modes, the stress parametrization does quantitatively effect the properties of the Fourier–Green solutions of inward propagating modes.

The amplitude of the Newtonian Fourier–Green’s function of the mass accretion rate, produced with x0 = 10 and x = 20, for a variety of different indices ν.

Figure 2.

The amplitude of the Newtonian Fourier–Green’s function of the mass accretion rate, produced with x₀ = 10 and x = 20, for a variety of different indices ν.

Open in new tab Download slide

The amplitude of the Newtonian Fourier–Green’s function of the mass accretion rate, produced with x0 = 10 and x = 9, for a variety of different indices ν. The stress index modifies the properties of the mass accretion rate Fourier–Green’s functions at intermediate frequencies.

Figure 3.

The amplitude of the Newtonian Fourier–Green’s function of the mass accretion rate, produced with x₀ = 10 and x = 9, for a variety of different indices ν. The stress index modifies the properties of the mass accretion rate Fourier–Green’s functions at intermediate frequencies.

Open in new tab Download slide

The angle of the Newtonian Fourier–Green’s function is defined as

$$\begin{eqnarray} \Phi _{\widetilde{G}_{\dot{M}}} = {\rm arg}\left(\widetilde{G}_{\dot{M}}\right) \equiv \tan ^{-1} \left({{\cal I} \left[ \widetilde{G}_{\dot{M}} \right] \over {\cal R} \left[ \widetilde{G}_{\dot{M}} \right]} \right) , \end{eqnarray}$$

(51)

where |${\cal I}[z]$| and |${\cal R}[z]$| denote the imaginary and real parts of the complex variable z, respectively, and the appropriate branches of the tan ⁻¹ function are chosen so that −π < arg(z) < π. The angle of these Fourier–Green’s functions is related to the time lag, a more useful observable quantity, via

$$\begin{eqnarray} t_{\rm lag} = { \Phi _{\widetilde{G}_{\dot{M}}} \over 2\pi f} . \end{eqnarray}$$

(52)

The angle of the Newtonian Fourier–Green’s functions, as a function of disc radius x, behave qualitatively as follows. At low Fourier frequencies, the angle of the Fourier–Green’s functions are approximately constant, whereas at high Fourier frequencies the angle of the Fourier–Green’s functions cycle through π → −π (or −π → π for x < x₀) as a function of increasing radius. The angle of the Fourier–Green’s functions are discontinuous at x = x₀.

The cycling of the Fourier–Green’s function angle at high frequencies is simple to understand analytically, using the results of the proceeding section:

$$\begin{eqnarray} \widetilde{G}_{\dot{M}} (x, x_0, f\rightarrow \infty) \sim C \exp (- i \sqrt{\pi f} \left| g(x) - g(x_0) \right|) , \end{eqnarray}$$

(53)

where C is purely real, meaning

$$\begin{eqnarray} \tan \left(\Phi _{\widetilde{G}_{\dot{M}}} \right) \rightarrow \tan \left[- \sqrt{\pi f} \left(g(x) - g(x_0) \right) \right], \end{eqnarray}$$

(54)

resulting in a cyclic behaviour.

The mass accretion rate power density spectrum |$S_{\dot{m}}(x, f)$|⁠, defined in Section 2.5, is shown in Fig. 4. In Fig. 4, we take an accretion rate fluctuation coherence length l(x′) = 1, and in the Laplacian input for |$S_\Sigma$| we take p = 1, and |$f_{\rm br} = \sqrt{GM/r^3}$|⁠, the Keplerian frequency at disc radius r. Also displayed as black dashed curves are the high and low frequency asymptotic results derived in Section 2.5. The amplitude of the accretion rate power density spectrum increases with decreasing disc radius, a result which is particularly true at high Fourier frequencies.

Figure 4.

The power density spectrum of the mass accretion rate multiplied by the Fourier frequency, at a number of disc radii x (denoted in legend), as a function of Fourier frequency. Displayed as black dashed curves are the high and low frequency asymptotic results derived in Section 2.5. The amplitude of the accretion rate power density spectrum increases with decreasing disc radius, a result which is particularly true at high Fourier frequencies.

Open in new tab Download slide

The complex cross-spectrum, which measures the correlation between variability at disc radii x₁ and x₂, is given by (Mushtukov et al. 2018)

$$\begin{eqnarray} C_{\dot{m}}(x_1, x_2, f) &=& \int _{x_{\rm in}}^{x_{\rm out}} \widetilde{G}_{\dot{M}}(x_1, x^{\prime }, f) \widetilde{G}^\dagger _{\dot{M}}(x_2, x^{\prime }, f) \\ && \times \, \left({1 \over x^{\prime }}\right)^2 l(x^{\prime }) S_\Sigma (x^{\prime }, f) \, {\rm d}x^{\prime }, \end{eqnarray}$$

(55)

where z^† denotes the complex conjugate of z. We define the coherence function

$$\begin{eqnarray} {\rm Coh}_{\dot{m}}(x_1, x_2, f) \equiv {\left| C_{\dot{m}}(x_1, x_2, f) \right|^2 \over S_{\dot{m}}(x_1, f) S_{\dot{m}}(x_2, f)} \end{eqnarray}$$

(56)

which is limited to the range |$0 \lt {\rm Coh}_{\dot{m}}(x_1, x_2, f) \lt 1$|⁠. As its name suggests, the coherence function measures the level of coherence between two variable quantities, in this case the variability at radii x₁ and x₂ in the disc. Note that Coh = 1 corresponds to a fully coherent fluctuations at both radii, while Coh = 0 represents incoherent fluctuations.

The coherence function is plotted as a function of Fourier frequency in Fig. 5, for a number of disc radii x₁, and x₂ = 50. At high enough Fourier frequencies the coherence between any two disc radii decays exponentially. At low Fourier frequencies the coherence between disc radii tends to unity. Clearly the coherence function is a complicated function of Fourier frequency at intermediate frequency scales.

The coherence function of the mass accretion rate, between disc radii x1 (displayed on legend), and x2 = 50. At high enough Fourier frequencies the coherence between any two disc radii decays exponentially. At low Fourier frequencies the coherence between disc radii tends to 1.

Figure 5.

The coherence function of the mass accretion rate, between disc radii x₁ (displayed on legend), and x₂ = 50. At high enough Fourier frequencies the coherence between any two disc radii decays exponentially. At low Fourier frequencies the coherence between disc radii tends to 1.

Open in new tab Download slide

In Fig. 6, we plot the angle of the mass accretion rate cross spectrum between disc radii x₁ (displayed on legend), and x₂ = 50. At low Fourier frequencies disc variability at radii x₁ leads variability at disc radii x₂ > x₁, whereas disc variability at x₁ > x₂ lags variability at radii x₂. At high Fourier frequencies variability at small radii (e.g. the blue curve) can lag, rather than lead, variability at larger radii. At the highest Fourier frequencies the cross spectrum cycles from −π → π, as a function of frequency.

The angle of the mass accretion rate cross spectrum between disc radii x1 (displayed on legend), and x2 = 50. At low Fourier frequencies disc variability at radii x1 leads variability at disc radii x2 > x1, whereas disc variability at x1 > x2 lags variability at radii x2. At high Fourier frequencies variability at small radii (e.g. the blue curve) can lag, rather than lead, variability at larger radii. At the highest Fourier frequencies, the cross spectrum cycles from −π → π, as a function of increasing frequency.

Figure 6.

The angle of the mass accretion rate cross spectrum between disc radii x₁ (displayed on legend), and x₂ = 50. At low Fourier frequencies disc variability at radii x₁ leads variability at disc radii x₂ > x₁, whereas disc variability at x₁ > x₂ lags variability at radii x₂. At high Fourier frequencies variability at small radii (e.g. the blue curve) can lag, rather than lead, variability at larger radii. At the highest Fourier frequencies, the cross spectrum cycles from −π → π, as a function of increasing frequency.

Open in new tab Download slide

In this section, we have calculated a number of the properties of the Newtonian Fourier–Green’s functions, showing that their properties are exactly as predicted by the analytical analysis of the previous section. For the remainder of this paper we shall focus on the properties of relativistic discs.

4 Relativistic Fourier–Green’s functions

The Fourier–Green’s functions of the previous section are solutions of the Newtonian disc equation. A large number of observed variable disc systems, however, are those discs evolving around black holes, and must therefore be described by a relativistic theory. Analytical solutions of the relativistic disc equations have recently been derived by Mummery (2023), and are discussed below.

4.1 Relativistic Green’s functions in the time domain

Balbus (2017) demonstrated that in full general relativity the governing disc equation can be expressed in the following compact form

(57)

Here the primed notation ′ denotes an ordinary derivative with respect to r, and

$$\begin{eqnarray} \zeta \equiv {r \Sigma {W^r_\phi }\over U^0} . \end{eqnarray}$$

(58)

Only two of the orbital components of the disc’s flow appear in this equation. These are the time dilation factor

$$\begin{eqnarray} U^0 = \frac{1+a\sqrt{{r_g}/{r^3}}}{\left({1 - {3r_g}/{r} + 2a\sqrt{{r_g}/{r^3} } }\right)^{1/2}} , \end{eqnarray}$$

(59)

and the circular orbit angular momentum gradient

$$\begin{eqnarray} U_\phi ^{\prime } = \frac{\sqrt{GM} \left(a\sqrt{r_g} + r^{{3}/{2}} \right) \left(r^2 - 6r_g r - 3a^2 + 8a\sqrt{r_g r}\right)}{2r^4 \left(1 - {3r_g}/{r} + 2a\sqrt{{r_g}/{r^3}} \right)^{{3}/{2}}}. \end{eqnarray}$$

(60)

In these expressions, a is the black hole’s angular momentum parameter (having dimensions of length), M is the black hole’s mass, r_g = GM/c² the gravitational radius, and G and c are Newton’s constant and the speed of light, respectively.

Clearly, the relativistic disc equation (57) is extremely algebraically complex, and in fact the general relativistic thin disc equation does not have exact solutions that can be written in terms of elementary functions. However, in a recent work Mummery (2023) derived the leading order general relativistic Green’s function solution, valid for the case of a vanishing stress at the innermost stable circular orbit (ISCO). The Mummery (2023) Green’s function solution has the same functional form as the Newtonian solutions discussed above, but with (note that these solutions are in units where G = M = w = 1, see Appendix A for the relevant solutions presented in physical units)

$$\begin{eqnarray} g(x) &=& {x^\alpha \over 2 \alpha } \sqrt{1 - {2\over x}}\left[1 - {x^{ - 1} \over { (\alpha - 1)}} {}_2F_1\left(1, {3\over 2}-\alpha ; 2-\alpha ; {2\over x}\right) \right] \\ && + {2^{\alpha - 2} \over \alpha (\alpha - 1)}\sqrt{\pi } {\Gamma (2-\alpha) \over \Gamma ({3/ 2} - \alpha)} , \end{eqnarray}$$

(61)

$$\begin{eqnarray} q(x) = x^{1/4} \sqrt{x^{-\alpha } g(x)} \exp \left({1 \over 2 x}\right) \left[1 - {2\over x}\right]^{5/4 - 3/8\alpha } , \end{eqnarray}$$

(62)

and

$$\begin{eqnarray} p(x) = {x^{1/2} \exp \left(-{1/ x}\right) \over 1 - {2/ x}} , \end{eqnarray}$$

(63)

where

$$\begin{eqnarray} x \equiv {2 r \over r_I}, \quad \alpha = {1\over 4\nu } , \end{eqnarray}$$

(64)

and ₂F₁(a, b; c; z) is the hypergeometric function. For a full description of the approximations employed in deriving this solution see Mummery (2023). Note that in the x → ∞ limit, the above solutions revert to their Newtonian form.

While the above solutions are not formally exact, in Fig. 7 we plot the numerically (blue dashed curve) and analytically (green solid curve) computed |$r\Sigma {W^r_\phi }$| profiles, assuming an initial radius of r₀ = 50r_g and Kerr angular momentum parameter a = 0. The curves are plotted at dimensionless times t/t_visc = 0.003, 0.015, 0.045, 0.09, 0.15, 0.225, 0.45, and 4.5, the curves at later times are identifiable through their decreasing peak amplitude. It is remarkable how accurately the analytical Green’s function solutions described above reproduces the properties of the full numerical solutions.

$Upper: The Green’s function solution of the variable $y \equiv r \Sigma {W^r_\phi }$, for a Kerr black hole with spin a = 0. The blue dashed curves are the numerical solutions of the full general relativistic disc equations, while the green solid curves are the analytical solution of Mummery (2023). The initial radius was r0 = 50rg and the curves are plotted at dimensionless times t/tvisc = 0.003, 0.015, 0.045, 0.09, 0.15, 0.225, 0.45, and 4.5. Lower: The absolute value of the difference between the numerical and analytical Green’s function solutions. To allow a proper comparison at each time both the numerical and analytical Green’s functions are renormalized to have a peak amplitude of 1. This figure is reproduced from Mummery (2023).$

Figure 7.

Upper: The Green’s function solution of the variable |$y \equiv r \Sigma {W^r_\phi }$|⁠, for a Kerr black hole with spin a = 0. The blue dashed curves are the numerical solutions of the full general relativistic disc equations, while the green solid curves are the analytical solution of Mummery (2023). The initial radius was r₀ = 50r_g and the curves are plotted at dimensionless times t/t_visc = 0.003, 0.015, 0.045, 0.09, 0.15, 0.225, 0.45, and 4.5. Lower: The absolute value of the difference between the numerical and analytical Green’s function solutions. To allow a proper comparison at each time both the numerical and analytical Green’s functions are renormalized to have a peak amplitude of 1. This figure is reproduced from Mummery (2023).

Open in new tab Download slide

In Fig. 8, we show the Green’s function of the mass accretion rate, plotted as a function of radius, for the spin parameter a = 0, at a number of different dimensionless times denoted on each plot. The initial radius in both cases was taken to be r₀ = 25r_g. A value of |$\dot{M}(r, t) \lt 0$| denotes mass inflow (towards the ISCO), while |$\dot{M}(r, t) \gt 0$| denotes mass outflow (some mass must move outwards within the disc so as to conserve the total angular momentum of the flow). The normalization of the accretion rate was chosen so that the time-integrated ISCO accretion rate was equal to 1.

The Green’s function solution of the mass accretion rate for a Schwarzschild black hole (a = 0). The initial radius was r0 = 25rg and the curves are plotted at the dimensionless times denoted in the legend. The normalization of the accretion rate was chosen so that the time-integrated ISCO accretion rate was equal to 1. This figure is reproduced from Mummery (2023).

Figure 8.

The Green’s function solution of the mass accretion rate for a Schwarzschild black hole (a = 0). The initial radius was r₀ = 25r_g and the curves are plotted at the dimensionless times denoted in the legend. The normalization of the accretion rate was chosen so that the time-integrated ISCO accretion rate was equal to 1. This figure is reproduced from Mummery (2023).

Open in new tab Download slide

4.2 Relativistic Green’s functions in the frequency domain

With this solution in hand, we may write down the Fourier–Green’s function solutions of the relativistic disc equation by substituting the above definitions into equation (18). The full expressions for the relativistic Fourier–Greens functions are rather complex, and are presented in Appendix A. These solutions have the same gross properties in the Fourier domain as their Newtonian analogues, but they differ in the details, as we now discuss.

In Figs 9 and 10, we plot the amplitude of the general relativistic Fourier–Green’s functions of the mass accretion rate, for inward (Fig. 9) and outward (Fig. 10) propagating modes. In black, we plot the corresponding Newtonian Fourier–Green’s solutions. We note that for inward propagating modes the relativistic Fourier–Green’s modes are more strongly suppressed at high frequencies, but are larger at intermediate frequencies, than their Newtonian analogues. In contrast, for outward propagating modes the relativistic Fourier–Green’s modes are more strongly suppressed at all frequencies.

The amplitude of the general relativistic Fourier–Green’s functions of the mass accretion rate (red), for inward propagating modes x < x0. This calculation was for a Kerr black hole with a = 0.5. Newtonian modes are displayed in black. For this figure, we take r0 = 10rg, and the inward propagating modes have r/rg = 4 (dotted–dashed), 5 (dotted), 7 (dashed), and 9 (solid). The Fourier–Green’s functions are normalized so that the total accreted mass is equal to 1. We note that the relativistic Fourier–Green’s modes are more strongly suppressed at high frequencies, but are larger at intermediate frequencies.

Figure 9.

The amplitude of the general relativistic Fourier–Green’s functions of the mass accretion rate (red), for inward propagating modes x < x₀. This calculation was for a Kerr black hole with a = 0.5. Newtonian modes are displayed in black. For this figure, we take r₀ = 10r_g, and the inward propagating modes have r/r_g = 4 (dotted–dashed), 5 (dotted), 7 (dashed), and 9 (solid). The Fourier–Green’s functions are normalized so that the total accreted mass is equal to 1. We note that the relativistic Fourier–Green’s modes are more strongly suppressed at high frequencies, but are larger at intermediate frequencies.

Open in new tab Download slide

The amplitude of the general relativistic Fourier–Green’s functions of the mass accretion rate (blue), for outward propagating modes x > x0. This calculation was for a Kerr black hole a = 0.5. Newtonian modes are displayed in black. For this figure, we take r0 = 10rg, and the outward propagating modes have r/rg = 11 (dot–dashed), 13 (dotted), 15 (dashed), and 17 (solid). The Fourier–Green’s functions are normalized so that the total accreted mass is equal to 1. We note that the outward propagating relativistic Fourier–Green’s modes are more strongly suppressed at all frequencies, when compared to Newtonian modes.

Figure 10.

The amplitude of the general relativistic Fourier–Green’s functions of the mass accretion rate (blue), for outward propagating modes x > x₀. This calculation was for a Kerr black hole a = 0.5. Newtonian modes are displayed in black. For this figure, we take r₀ = 10r_g, and the outward propagating modes have r/r_g = 11 (dot–dashed), 13 (dotted), 15 (dashed), and 17 (solid). The Fourier–Green’s functions are normalized so that the total accreted mass is equal to 1. We note that the outward propagating relativistic Fourier–Green’s modes are more strongly suppressed at all frequencies, when compared to Newtonian modes.

Open in new tab Download slide

It is important to note that the outward propagating relativistic Fourier–Green’s modes must be treated with some care. Mathematically this results from the fact that the solutions described in the proceeding section are not exact, but are asymptotic ‘leading order’ solutions [see Mummery (2023) for a detailed discussion]. Physically care is required because at very late times the exact numerical and analytical solution begin to deviate, and an extremely small fraction of the initial disc mass is not accreted in these solutions. As a result

$$\begin{eqnarray} q(x) \ne A (g(x))^{\nu } \end{eqnarray}$$

(65)

and therefore (following the reasoning outlined in Section 2)

$$\begin{eqnarray} \lim _{f \rightarrow 0} \widetilde{G}_{\dot{M}}(x \gt x_0, f) \rightarrow \delta (x, x_0) M_d \ne 0 . \end{eqnarray}$$

(66)

The discrepancy is small, as a result of the high accuracy of these analytical solutions (Fig. 7), and typically only effects extremely small frequencies f/f₀ ≪ 10⁻⁶ to a small degree

$$\begin{eqnarray} \delta \ll 10^{-4}. \end{eqnarray}$$

(67)

As such, the conclusions derived in the remainder of this paper are not noticeably quantitatively or qualitatively effected by this slight discrepancy. The function δ(x, x₀) can be written in a closed form (Appendix A), and then simply subtracted off the low frequency Fourier–Green’s modes.

4.3 Differences between the Newtonian and relativistic Fourier–Green’s functions

One of the main differences between the Newtonian and relativistic Fourier–Green’s functions is highlighted in Figs 9 and 10: relativistic Fourier–Green solutions are more strongly suppressed at high Fourier frequencies than their Newtonian analogues. This can be understood by plotting the high-frequency mode suppression function Δ(r, r₀) ≡ |g(r) − g(r₀)| (Fig. 11). High frequency modes are suppressed as as exp (−Δf^1/2) (see Section 2). In Fig. 11, we plot Δ(r, r₀) as a function of black hole spin for r₀ = 20r_g. We see that relativistic modes are more strongly suppressed at high Fourier frequencies than Newtonian modes, but that higher black hole spins reduce this suppression.

The high-frequency mode suppression function Δ(r, r0) ≡ |g(r) − g(r0)|, which suppresses modes as exp (−Δf1/2) at high Fourier frequencies, as a function of black hole spin for r0 = 20rg. We see that relativistic modes are more strongly suppressed at high Fourier frequencies than Newtonian modes, but that higher black hole spins reduces this suppression.

Figure 11.

The high-frequency mode suppression function Δ(r, r₀) ≡ |g(r) − g(r₀)|, which suppresses modes as exp (−Δf^1/2) at high Fourier frequencies, as a function of black hole spin for r₀ = 20r_g. We see that relativistic modes are more strongly suppressed at high Fourier frequencies than Newtonian modes, but that higher black hole spins reduces this suppression.

Open in new tab Download slide

One interesting result that is highlighted in Fig. 9 is that the inwards propagating relativistic and Newtonian Fourier–Green’s functions differ at intermediate frequencies, with Newtonian modes having a more pronounced ‘kink’ at f ∼ 10⁻²f₀. This difference propagates through into different properties of the local accretion rate power spectra of relativistic and Newtonian discs (Fig. 12). In Fig. 12, we compute the local (at radius r) mass accretion rate power spectrum using equation (35), for a range of different black hole spins and also for Newtonian discs with different ‘inner radii’. To mimic an effective inner disc radii for the Newtonian solutions we set the input power of the surface density variability (see equation 36) to be zero for r ≤ r_in. It is clear to see in Fig. 12 that the relativistic and Newtonian Fourier–Green’s functions result in significantly different local accretion rate power spectra. While this is particularly true for the case of a Newtonian disc with small inner radius r_in = 0.01 (black dashed curve), this is still true for Newtonian discs with more realistic inner radii (red dashed curve).

Figure 12.

The power density spectrum of the mass accretion rate multiplied by the Fourier frequency, at a disc radius r = 10r_g for a number of black hole spins (denoted in legend), as a function of Fourier frequency. Displayed as black and red dashed curves are the equivalent Newtonian solutions with different ‘inner radii’ (the radius at which the input power is set to zero, see text). Local relativistic and Newtonian accretion rate power density spectra are noticeably distinct.

Open in new tab Download slide

5 OBSERVABLE QUANTITIES

The previous sections have focused on the properties of the Green’s function solutions of the Newtonian and relativistic disc equations in the Fourier domain. The Fourier properties of the local mass accretion rate are of course not by themselves directly observable, and in this section we focus our discussion onto the properties of the emitted disc flux, which is a more readily observable quantity.

The problem now becomes one of determining how mass accretion rate fluctuations are observed in the light curves of accreting sources. In a purely steady state flow the local energy available to be radiated is directly proportional to the local mass accretion rate, and in the literature it is generally assumed that the local variability in the accretion rate is directly proportional to the variability in the local emission (e.g. Lyubarskii 1997; Ingram & van der Klis 2013; Ingram & Done 2012; Mushtukov et al, 2018). This is unlikely to hold much beyond small fluctuations in |$\dot{M}$|⁠, as the locally radiated energy in a time dependent flow is not directly proportional to the accretion rate, but instead to the energy liberated by the local disc shear. In a time dependent accretion flow the local mass accretion rate can in fact be negative (see Fig. 8), while the locally emitted flux will of course remain positive.

In this work, in common with the literature, we will also make the assumption that the variability in the local accretion rate is directly proportional to the variability in the locally emitted photon flux. We stress however that it is important to bear in mind that this assumption can only hold for small variations in both the photon flux and mass accretion rate.

The emission in a certain band (which we shall call ‘hard’ h, or ‘soft’ s) is therefore assumed to be correlated with the local accretion rate, with some ‘emissivity profile’ s(r) and h(r) (e.g. Lyubarskii 1997; Ingram & van der Klis 2013; Ingram & Done 2012; Mushtukov et al, 2018) which specifies the constant of proportionality between the two quantities. Much of the variable emission in a typical observation of an X-ray binary system is sourced from a hot ‘corona’. The corona is a population of hot electrons, located close to the black hole, whose existence is inferred from the observation of a power law spectral component resulting physically from photons being Compton up-scattered as they pass through the electrons. The geometry of this corona is currently poorly understood, and could in principal be located above the disc, or as part of a hot ‘inner flow’. The advantage of working with emissivity profiles, as opposed to any particular physical model for the emission, is that under the assumption that coronal emission also scales locally with the mass accretion rate (perhaps as a result of variability in the seed photon field), variability in both thermal and non-thermal emission components may be modelled with additional degrees of freedom. We discuss potential extensions to this analysis in Section 7.

It is common in the literature to assume that these emissivity profiles are given by power laws of disc radius r (although there is no compelling justification for this), and in this work we shall parameterise our emissivity profiles with the following functional form

$$\begin{eqnarray} s(r) &= s_0 \left({r \over r_I}\right)^{-\gamma _s} \left({1 - \sqrt{r_I \over r} } \right), \end{eqnarray}$$

(68)

$$\begin{eqnarray} h(r) &= h_0 \left({r \over r_I}\right)^{-\gamma _h} \left({1 - \sqrt{r_I \over r} } \right) . \end{eqnarray}$$

(69)

Here r_I is the black hole’s ISCO radius, γ_{h, s} are phenomenological emissivity indices, and the final term |$1-\sqrt{r_I/r}$| enforces the vanishing ISCO stress condition used in deriving the relativistic Green’s functions. The constants s₀ and h₀ are arbitrary scaling factors included for dimensional reasons which we shall simply set equal to 1 for the remainder of this paper. We require γ_h > γ_s ≥ 2, so that the hard emission is produced at radii interior to the soft emission (γ_h > γ_s) and that the flux observed in either band falls off with radius at least as quickly as the total liberated photon flux (γ_{h, s} ≥ 2). The first requirement, that γ_h > γ_s, is in effect simply the definition of the ‘hard’ and ‘soft’ bands: the flux from the inner regions is emitted at higher photon energies, and will thus contribute more to harder (higher energy) bands. As the disc cools at larger radii, the relative contribution of that disc region to harder bands will be suppressed more than its contribution to softer bands.

In the following three subsections we shall present formal expressions for three readily observable quantities: the power density spectrum of a band, and the cross spectrum and coherence between different bands. It will then be demonstrated that each of these quantities is a relatively strong function of black hole spin, and could in principle be used to constrain the black hole spins of astrophysical sources.

5.1 The power density spectrum of a band

The first observable quantity we consider is the power density spectrum of the flux variability in a given observing band.

For simplicity we shall quote the results derived by Mushtukov et al. (2018), before discussing the assumptions inherent to the modelling. The following expression for the power density spectrum of a band (we here denote the band by h, for ‘hard band’ emission) assumes that the total luminosity available to be radiated in a given region of the accretion flow is proportional to the local mass accretion rate |$\dot{M}(r, t)$|⁠. Further assuming that the local variability of the mass accretion rate (⁠|$\dot{m}$|⁠) is small in comparison with the average mass accretion rate, the fluctuations of the flux in some energy band will also be proportional to |$\dot{m} (r, t)$|⁠. Under the above assumptions Mushtukov et al. (2018) derived the following expression for the power density spectrum of band h, denoted S_h(f)

$$\begin{eqnarray} S_h(f) = \int _{\cal D} \int _{\cal D} h(r_1) h(r_2) C_{\dot{m}} (r_1, r_2, f) \, {\rm d}r_1 \, {\rm d}r_2 , \end{eqnarray}$$

(70)

where we use the shorthand

$$\begin{eqnarray} \int _{\cal D} f(r) \, {\rm d}r\equiv \int _{r_{\rm in}}^{r_{\rm out}} f(r) \, {\rm d}r , \end{eqnarray}$$

(71)

to indicate an integral over the entire disc surface |${\cal D}$|⁠. We remind the reader that

$$\begin{eqnarray} C_{\dot{m}}(r_1, r_2, f) &=& \int _{\cal D} \widetilde{G}_{\dot{M}}(r_1, r^{\prime }, f) \, \widetilde{G}^\dagger _{\dot{M}}(r_2, r^{\prime }, f) \\ && \times \,\left({1 \over r^{\prime }}\right)^2 l(r^{\prime }) S_\Sigma (r^{\prime }, f) \, {\rm d}r^{\prime }, \end{eqnarray}$$

(72)

and so we have in reality a triple integral to compute

$$\begin{eqnarray} S_h(f) &=& \int _{\cal D} \int _{\cal D} \int _{\cal D} h(r_1) h(r_2) \widetilde{G}_{\dot{M}}(r_1, r^{\prime }, f) \, \widetilde{G}^\dagger _{\dot{M}}(r_2, r^{\prime }, f) \\ && \times \, \left({1 \over r^{\prime }}\right)^2 l(r^{\prime }) S_\Sigma (r^{\prime }, f) \, {\rm d}r^{\prime } \, {\rm d}r_1 \, {\rm d}r_2 . \end{eqnarray}$$

(73)

Note that, as argued in Mushtukov et al. (2018), as the amplitude of the local surface density variability scales with |$\dot{M}^2$|⁠, integrating S_h(f) over all frequencies, and then square rooting, one recovers the linear rms-flux relationship.

The physical reason behind the power density spectrum being related to a triple integral is the following. The flux in a given energy band is, per our assumptions, given by an integral of the local mass accretion rate, with a weighting function h(r), over the entire disc |$F_h(t) \propto \int _{\cal D} h(r^{\prime }) \dot{M}(r^{\prime }, t) \, {\rm d}r^{\prime }$|⁠. The power density spectrum corresponds to the square of the Fourier-transformed flux |$S_h(f) \equiv \widetilde{F}_h(f) \widetilde{F}^\dagger _h(f)$|⁠, which introduces the second integration. However, variability in the accretion rate at a given radius r is caused by the integrated contributions of surface density fluctuations at all disc radii, with differing levels of correlation encapsulated by |$C_{\dot{m}}(r_1, r_2, f)$|⁠. Summing each of these individual contributions introduces the third integral over the disc surface.

5.2 The cross spectrum between bands

Observations of an accreting system may be contemporaneously taken at numerous different photon energies (or ‘bands’). A natural observational question is then how strongly is the variability in the emission observed across a ‘soft’ band correlated with the variability in the emission observed across a ‘hard’ band. This correlation is quantified by the hard-soft cross spectrum, a quantity given explicitly by (Mushtukov et al. 2018)

$$\begin{eqnarray} C_{h, s}(f) = \int _{\cal D} \int _{\cal D} h(r_1) s(r_2) C_{\dot{m}} (r_1, r_2, f)\, {\rm d}r_1 \, {\rm d}r_2 . \end{eqnarray}$$

(74)

The physical cause of this triple disc integral is identical to that of the power density spectrum. The hard-soft cross spectrum is a complex quantity, and its phase encapsulates an important physical quantity, namely the phase-lag in fluctuations in the hard sate emission with respect to the soft state emission. This phase lag is explicitly equal to

$$\begin{eqnarray} \tan \Phi _{h,s}(f) = {{\cal I}[C_{h, s}(f) ] \over {\cal R}[C_{h, s}(f) ] } , \end{eqnarray}$$

(75)

where |${\cal I}[z]$| and |${\cal R}[z]$| represent the real and imaginary parts of z respectively. This quantity can also be equivalently expressed as a time lag between hard and soft fluctuations:

$$\begin{eqnarray} t_{{\rm lag}\, h,s }(f) = {\Phi _{h,s}(f) \over 2\pi f} . \end{eqnarray}$$

(76)

5.3 The coherence between bands

The final observable quantity we shall discuss in this paper is the coherence of the fluctuations in hard and soft observing bands, a quantity explicitly given by

$$\begin{eqnarray} {\rm Coh}_{h, s}(f) = {|C_{h, s}(f)|^2 \over S_h(f) S_s(f)} . \end{eqnarray}$$

(77)

The coherence function satisfies 0 ≤ Coh_{h, s} ≤ 1, where Coh = 1 corresponds to a fully coherent fluctuations in both bands, while Coh = 0 represents incoherent fluctuations.

6 THE BACK HOLE SPIN DEPENDENCE OF OBSERVABLE QUANTITIES

The relativistic Fourier–Green’s functions derived in this paper are functions of the black holes spin through their explicit dependence on the ISCO radius of the black hole (equation 64). This explicit spin dependence can be seen in Figs 13, 14, and 15, where we plot the amplitude of the relativistic Fourier–Green’s functions, the phase and coherence of mass accretion rate fluctuations for a number of different black hole spins. An interesting result seen in Figs 13 and 14 is that the higher the black hole spin is made the closer to the Newtonian Fourier–Green’s function (black dashed curves) the solution becomes. In Fig. 15, we highlight that the black hole spin quantitatively effects the properties of the local mass accretion rate.

The amplitude of the general relativistic Fourier–Green’s functions of the mass accretion rate, for inward propagating modes r < r0, for black holes of different Kerr spin parameters. Newtonian modes are displayed in black. For this figure, we take r0 = 20rg, and the inward propagating modes have r/rg = 15.

Figure 13.

The amplitude of the general relativistic Fourier–Green’s functions of the mass accretion rate, for inward propagating modes r < r₀, for black holes of different Kerr spin parameters. Newtonian modes are displayed in black. For this figure, we take r₀ = 20r_g, and the inward propagating modes have r/r_g = 15.

Open in new tab Download slide

The amplitude of the general relativistic Fourier–Green’s functions of the mass accretion rate, for outward propagating modes r > r0, for black holes of different Kerr spin parameters. Newtonian modes are displayed in black. For this figure, we take r0 = 15rg, and the inward propagating modes have r/rg = 20.

Figure 14.

The amplitude of the general relativistic Fourier–Green’s functions of the mass accretion rate, for outward propagating modes r > r₀, for black holes of different Kerr spin parameters. Newtonian modes are displayed in black. For this figure, we take r₀ = 15r_g, and the inward propagating modes have r/r_g = 20.

Open in new tab Download slide

Upper: The coherence function of the mass accretion rate at r1 = 10rg and r2 = 15rg for Kerr metrics with different spin parameters (denoted on legend). Lower: The phase of the mass accretion rate cross spectrum between disc radii r1 = 10rg and r2 = 15rg. At higher frequencies than displayed in this lower plot the phase wraps between −π and π. Both the coherence and phase of the mass accretion rate are relatively sensitive to the Kerr metric spin parameter.

Figure 15.

Upper: The coherence function of the mass accretion rate at r₁ = 10r_g and r₂ = 15r_g for Kerr metrics with different spin parameters (denoted on legend). Lower: The phase of the mass accretion rate cross spectrum between disc radii r₁ = 10r_g and r₂ = 15r_g. At higher frequencies than displayed in this lower plot the phase wraps between −π and π. Both the coherence and phase of the mass accretion rate are relatively sensitive to the Kerr metric spin parameter.

Open in new tab Download slide

In this section, we demonstrate how this spin dependence is also present in the integrated observable properties of the photon flux emitted from the disc surface. As an explicit example of this result, see Fig. 16, where we plot the power density spectrum of the hard band (upper left), the hard–soft coherence (upper right), the hard–soft phase lag (lower left) and time delay (lower right) for four different black hole spins (a = 0, blue; a = 0.5, black; a = +0.99, orange; and a = −0.99, green). For the time delay plot, we denote by solid dots positive lags (hard lags soft), and by crosses negative lags (hard leads soft). The units of the power spectrum are arbitrary and would be set in a physical system by the input power in the surface density variability, the units of frequency and time delay are scaled to a system with black hole mass M_BH = 10 M_⊙, and disc parameters α = H/R = 0.1. To scale these results to any other system the frequency axis should be scaled by a factor

$$\begin{eqnarray} N=\left({\alpha \over 0.1}\right) \left({H/R \over 0.1}\right)^2 \left({10M_\odot \over M_{\rm BH}}\right) , \end{eqnarray}$$

(78)

while the time-lags should be scaled by 1/N. This particular plot has emissivity parameters γ_s = 3, γ_h = 5. For the input surface density fluctuations we take a correlation length l(r) = 1r_g, an integrated power p = 1, and a break frequency equal to one percent of the local Keplerian frequency |$f_K = \sqrt{GM/r^{3}}$| (see equation 36).

Figure 16.

The power density spectrum of the hard band (upper left), the hard–soft coherence (upper right), the hard–soft phase lag (lower left) and time delay (lower right) for emissivity parameters γ_s = 3, γ_h = 5, and four different black hole spins (a = 0, blue; a = 0.5, black; a = +0.99, orange; and a = −0.99, green). For the time delay plot, we denote by solid dots positive lags (hard lags soft), and by crosses negative lags (hard leads soft). The units of the power spectrum are arbitrary and would be set in a physical system by the input power in the surface density variability, the units of frequency and time delay are scaled to a system with black hole mass M_BH = 10 M_⊙, and disc parameters α = H/R = 0.1.

Open in new tab Download slide

It is clear from Fig. 16 that the black hole spin has a substantial effect on the observed variability properties of a black hole accretion flow. This is a result of fundamental theoretical importance, and the key observational result of this paper.

As can be seen in the upper left panel of Fig. 16, the power density spectra of the accretion variability predicted by this theory are qualitatively simple. At high and low Fourier frequencies the observed power density spectra are given by simple power laws of frequency. The precise power-law indices of the high and low frequency slopes are determined by the input surface density variability profile, as discussed and derived in Section 2. It is interesting to note that the power density spectrum of a given band peaks at higher Fourier frequencies for larger black hole spins, but generally with a smaller magnitude.

The upper-right panel of Fig. 16 displays the coherence of the hard–soft variability. At low Fourier frequencies, the coherence function tends to a near-unity constant for each black hole spin, but at larger Fourier frequencies the coherence is in general an extremely non-trivial function frequency. It is interesting to note that roughly the frequency at which the phase lags turn from positive to negative (lower left-hand panel) the hard–soft variability becomes increasingly incoherent. We note two properties of the coherence as a function of black hole spin: the hard–soft variability is both inherently more coherent (larger Coh_h,s), and is inherently smoother as a function of Fourier frequency, for larger black hole spins.

A simple interpretation of this spin dependence of the coherence is that the emission from Kerr discs with higher spins primarily originates from a region of the disc which has a smaller radial extent when compared to lower spins. This is because the emission from highly spinning Kerr systems is dominated by the hottest and very innermost regions, which are physically close together. The coherence between two radii is a strong function of their separation, with larger separations having significantly lower coherence, and smaller separations a correspondingly larger coherence.

In the lower two panels, we plot the phase and time lags associated with the hard and soft bands. As has been discussed by many authors, it is natural that the variability in ‘hard’ bands lag the variability in ‘soft’ bands (solid dots, Fig. 16), a result of the fluctuations excited in the outer disc regions propagating inwards towards the hotter inner regions where the hard flux is generated. For these Fourier frequencies the time lags are roughly given by the typical time of propagation of fluctuations from outer radii, where the soft flux is predominantly produced, to inner radii, where the ‘hard’ flux is predominantly produced. However, negative phase lags (when the soft energy band variability lags hard energy variability) are possible at high frequencies (see lower left-hand and right-hand panels of Fig. 16). The negative lags are a result of the fact that the variability of the mass accretion rate at the inner radii can affect the variability at the outer radii through outward propagating modes (as argued first by Mushtukov et al. 2018).

These negative phase lags are only present at high Fourier frequencies, a result which can be understood with respect to the analytical analysis of Section 2. At high Fourier frequencies the Fourier–Green’s functions of inward and outward propagation are symmetric and therefore equally important. At low frequencies, however, outward propagation is suppressed (as a power-law in frequency), and inward propagation dominates resulting in positive lags.

We note that the predicted negative phase lags are comparable to the observed negative lags in stellar mass black hole systems (Uttley et al. 2011; De Marco et al. 2015) and AGN (Zoghbi et al. 2010; Walton et al. 2013; Alston, Done & Vaughan 2014). It is interesting that the fundamental diffusive accretion process in discs can in principle play a key role in the generation of negative time lags.

The precise quantitative results of Fig. 16 of course depend on the precise assumptions inherent to the modelling, namely the emissivity indices γ_s and γ_h. While varying these parameters changes the numerical values of each of the observed quantities (power spectrum, coherence, and phase lags), the qualitative spin dependence of the different quantities remains unchanged. This is demonstrated in Figs 17 and 18 where we again plot the power density spectrum of the hard band (upper left), the hard–soft coherence (upper right), the hard–soft phase lag (lower left) and time delay (lower right) for four different black hole spins (a = 0, blue; a = 0.5, black; a = +0.99, orange; and a = −0.99, green). In Fig. 17, we take γ_s = 2, γ_h = 3.5, while in Fig. 18 we take γ_s = 3, γ_h = 8.

Figure 17.

The power density spectrum of the hard band (upper left), the hard–soft coherence (upper right), the hard–soft phase lag (lower left) and time delay (lower right) for emissivity parameters γ_s = 2, γ_h = 3.5 and four different black hole spins (a = 0, blue; a = 0.5, black; a = +0.99, orange; and a = −0.99, green). For the time delay plot, we denote by solid dots positive lags (hard lags soft), and by crosses negative lags (hard leads soft). The units of the power spectrum are arbitrary and would be set in a physical system by the input power in the surface density variability, the units of frequency and time delay are scaled to a system with black hole mass M_BH = 10 M_⊙, and disc parameters α = H/R = 0.1.

Open in new tab Download slide

Figure 18.

The power density spectrum of the hard band (upper left), the hard–soft coherence (upper right), the hard–soft phase lag (lower left) and time delay (lower right) for emissivity parameters γ_s = 3, γ_h = 8 and four different black hole spins (a = 0, blue; a = 0.5, black; a = +0.99, orange; and a = −0.99, green). For the time delay plot, we denote by solid dots positive lags (hard lags soft), and by crosses negative lags (hard leads soft). The units of the power spectrum are arbitrary and would be set in a physical system by the input power in the surface density variability, the units of frequency and time delay are scaled to a system with black hole mass M_BH = 10 M_⊙, and disc parameters α = H/R = 0.1.

Open in new tab Download slide

The following ‘rules of thumb’ appear to describe well the spin-dependence of the observed flux variability, independent of the choice of precise emissivity profile

The magnitude of the maximum hard–soft phase lag decreases with increasing black hole spin.
The magnitude of the maximum negative hard–soft phase lag increases with increasing black hole spin.
The frequency at which the maximum negative hard–soft phase lag occurs increases with increasing black hole spin.
The power density spectrum of a given band peaks at higher Fourier frequencies for larger black hole spins.
The power density spectrum of a given band peaks with smaller magnitude for larger black hole spins.
The hard–soft variability is inherently more coherent (larger Coh_h,s) for larger black hole spins.
The hard–soft coherence is inherently smoother as a function of Fourier frequency for larger black hole spins.

The reason that the frequencies at which various key observational properties occur increase with black hole spin is simply as a result of the reducing ISCO radius of the more rapidly rotating Kerr space–time. These smaller radii have associated with them larger orbital and accretion frequencies, which are then observable in the variability signatures of these systems.

These rules of thumb are robust, and are not even dependent on the chosen functional form of the emissivity profiles. In Fig. 19, we once again plot the power density spectrum of the hard band (upper left), the hard–soft coherence (upper right), the hard–soft phase lag (lower left), and time delay (lower right) for four different black hole spins (a = 0, blue; a = 0.5, black; a = +0.99, orange; and a = −0.99, green). In this figure, however, we take emissivity profiles given by

$$\begin{eqnarray} s(r) &= s_0 \exp \left(-{r \over r_I}\right) \left({1 - \sqrt{r_I \over r} } \right) , \end{eqnarray}$$

(79)

$$\begin{eqnarray} h(r) &= h_0 \exp \left(- 3 {r \over r_I}\right) \left({1 - \sqrt{r_I \over r} } \right) . \end{eqnarray}$$

(80)

Figure 19.

The power density spectrum of the hard band (upper left), the hard–soft coherence (upper right), the hard–soft phase lag (lower left) and time delay (lower right) for exponential emissivity profiles (see text), and four different black hole spins (a = 0, blue; a = 0.5, black; a = +0.99, orange; and a = −0.99, green). For the time delay plot, we denote by solid dots positive lags (hard lags soft), and by crosses negative lags (hard leads soft). The units of the power spectrum are arbitrary and would be set in a physical system by the input power in the surface density variability, the units of frequency and time delay are scaled to a system with black hole mass M_BH = 10 M_⊙, and disc parameters α = H/R = 0.1.

Open in new tab Download slide

This could in fact be a more physically reasonable profile, as it may more accurately describe the suppression of hard and soft emission by the Wien-tail of the local blackbody disc emission function F(E) ∝ exp (− E_h,s/kT), E_h > E_s. In Fig. 19, we once again see that the precise numerical values of each of the observable quantities depends on the precise emissivity profile chosen, but that the gross spin-dependence of the variability is qualitatively unchanged.

7 FUTURE EXTENSIONS TO THE MODEL

The relativistic variability model presented in the previous section employed a number of simplifications. In this section, we recap and discuss these simplifications, their physical basis, and how they may be improved upon in future work.

The principal simplification employed in this analysis is in relating the variability in the mass accretion rate to the variability in the observed photon field. As we have discussed, in this work we employ phenomenological emissivity profiles. The advantage of working with emissivity profiles, as opposed to any particular physical model for the emission is that, provided that both thermal and coronal emission scale locally with the mass accretion rate, variability in both thermal and non-thermal emission components may be modelled with the additional degrees of freedom provided by the emissivity profiles.

However, a more detailed treatment of the local emission is of interest, and we briefly discuss a possible modelling approach below. Assuming that the disc emits thermally with an effective temperature profile T_eff(r), the flux at an observed energy E is proportional to

$$\begin{eqnarray} F \propto \int _{\cal D} {2 \pi r \over \exp \left(E / k T_{\rm eff}(r)\right) -1 } \, {\rm d} r, \end{eqnarray}$$

(81)

where this expression neglects the effects of gravitational lensing, gravitational redshifts, and the Doppler boosting of radiation. The variable flux δF is given, assuming a small fluctuation in |$\dot{M}$| and to linear order, by

$$\begin{eqnarray} \delta F \propto \int _{\cal D} {\delta \dot{M} \over \dot{M}} \left({E \over kT_{\rm eff}(r) }\right) {2 \pi r \exp \left(E / k T_{\rm eff}(r)\right) \over \left(\exp \left(E / k T_{\rm eff}(r)\right) -1 \right)^2} \, {\rm d} r, \end{eqnarray}$$

(82)

where we have assumed that |$\delta T/T = \delta \dot{M} / 4\dot{M}$|⁠. Therefore, to linear order, the thermal flux from an accretion flow can be treated by the same methods developed here, with an ‘emissivity profile’ given by equation (82). The treatment of non-thermal radiation would be more complex, but a variable input thermal flux of the form given above could in principle be propagated through a Comptonizing electron population.

At this level of detail however additional relativistic effects will likely become important. This would include energy-dependent delays in the propagation of photons through the space–time of Kerr black holes, with photons emitted deeper in the gravitational well of the black hole following paths to the observer more severely warped by gravity. In addition, the effects of gravitational and Doppler energy shifting of photons will likely modify the results presented here, at the quantitative level. These higher order effects will all themselves be sensitive to the black hole’s spin parameter, and a more detailed treatment of photon propagation and its effects on observed variability are of real interest.

A further improvement of the analysis presented here would be in improving the treating of the input surface density fluctuations. In this work, we have assumed that the input surface density perturbations are described by a Lorentzian profile with phenomenological parameters p and f_br (equation 36). While physically motivated, it would be of interest in future works to calibrate this model input with numerical analyses of the fundamental disc equations (e.g. Hogg & Reynolds 2016; Turner & Reynolds 2021).

8 CONCLUSIONS

In this paper, we have presented two important advances in the theoretical framework for describing aperiodic variability from accreting sources (the so-called theory of propagating fluctuations). First, we present the exact analytical solutions of the Fourier integral of the Green’s functions of the classical thin disc equations. With analytical solutions now at hand, various asymptotic properties of these solutions may be derived. In Section 2, we demonstrated that high frequency variability in the mass accretion rate is suppressed as exp (−Δf^1/2), where Δ(x, x′) is a function of the magnitude of the difference between the two disc locations x and x′, and corresponds physically to the (square root of) the accretion propagation time between x and x′. The high and low frequency asymptotic behaviour of the power spectrum of variability are also determined, and related to the intrinsic variability in the disc surface density/alpha parameter.

We have demonstrated that the power spectrum of the local mass accretion rate spectrum is, at high Fourier frequencies, dominated by locally driven variability, with exponentially small contributions from distant disc regions. At high Fourier frequencies the inward and outward propagation of material are equally important, with the Fourier–Green’s functions symmetric in x − x′ in this limit. At low Fourier frequencies; however, the variability is dominated by perturbations sourced at radii which are more distant from the central object, which then propagate inwards. Outward propagation is suppressed at low frequencies, as a power law in frequency.

In addition, these exact solutions will rapidly speed up the process of fitting analytical models of accretion variability to observational data; the numerical cost of Fourier transforming thin disc Green’s functions had previously been substantial.

The second key development is in presenting the first analysis of the Fourier–Green’s function solutions of the general relativistic thin disc equation. In this paper, we have presented the Fourier–Green’s function solutions valid for in a relativistic theory of gravity, under the assumption that the dynamical disc stress vanishes at the ISCO. These solutions depend implicitly on the central black hole’s spin through their dependence on the space–time’s ISCO radius.

We use this new theoretical development to highlight the Kerr black hole spin dependence of a number of observable variability properties of black hole discs. The power density spectrum of the hard band (upper left), the hard–soft coherence (upper right), the hard–soft phase lag (lower left), and time delay (lower right) are displayed in Figs 16, 17, 18, and 19 for four different black hole spins (a = 0, blue; a = 0.5, black; a = +0.99, orange; and a = −0.99, green), and a number of different parameterisations of the disc emissivity. Clearly the black hole spin imparts a strong signal on to the observable variability properties of black hole disc systems.

While the precise choice of emissivity profile of the hard and soft bands quantitatively effect the system’s observed variability properties, the following ‘rules of thumb’ appear to describe well the spin-dependence of the observed flux variability, independent of the choice of precise emissivity profile. These rules of thumb may be of use even in systems where a detailed analysis of the variability is not performed.

The magnitude of the maximum hard–soft phase lag decreases with increasing black hole spin.
The magnitude of the maximum negative hard–soft phase lag increases with increasing black hole spin.
The frequency at which the maximum negative hard–soft phase lag occurs increases with increasing black hole spin.
The power density spectrum of a given band peaks at higher Fourier frequencies for larger black hole spins.
The power density spectrum of a given band peaks with smaller magnitude for larger black hole spins.
The hard–soft variability is inherently more coherent (larger Coh_h,s) for larger black hole spins.
The hard–soft coherence is inherently smoother as a function of Fourier frequency for larger black hole spins.

The results presented in this paper therefore open up the possibility of using the aperiodic variability observed from black hole accretion systems to constrain the central black hole’s spin, a parameter of fundamental observational and theoretical interest.

DATA ACCESSIBILITY STATEMENT

No observational data was used in producing this manuscript. Python scripts which compute the Newtonian and relativistic Fourier–Green functions are available at https://github.com/andymummeryastro/GR_prop_fluc.

ACKNOWLEDGEMENTS

I would like to thank Alexander Mushtukov for interesting discussions which initiated this work. I am particularly grateful to Adam Ingram for extremely illuminating discussions regarding the propagating fluctuation model. I am grateful to the reviewer, whose detailed report strengthened the analysis in a number of places. This work was supported by a Leverhulme Trust International Professorship grant [number LIP- 202-014]. For the purpose of Open Access, I have applied a CC BY public copyright licence to any Author Accepted Manuscript version arising from this submission.

References

Abramowitz

Stegun

I. A.

1965

Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables

Dover

New York

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Alston

W. N.

Done

Vaughan

2014

MNRAS

439

1548

10.1093/mnras/stu005

Crossref

Search ADS

Balbus

S. A.

2017

MNRAS

471

4832

10.1093/mnras/stx1955

Crossref

Search ADS

Balbus

S. A.

Hawley

J. F.

1998

RvMP

10.1103/RevModPhys.70.1

Crossref

Search ADS

Balbus

S. A.

Hawley

J. F.

1991

ApJ

376

214

10.1086/170270

Crossref

Search ADS

Churazov

Gilfanov

Revnivtsev

2001

MNRAS

321

759

10.1046/j.1365-8711.2001.04056.x

Crossref

Search ADS

De Marco

Ponti

Cappi

Dadina

Uttley

Cackett

E. M.

Fabian

A. C.

Miniutti

2013

MNRAS

431

2441

10.1093/mnras/stt339

Crossref

Search ADS

De Marco

Ponti

Muñoz-Darias

Nandra

2015

ApJ

814

10.1088/0004-637X/814/1/50

Crossref

Search ADS

Eardley

D. M.

Lightman

A. P.

1975

ApJ

200

187

10.1086/153777

Crossref

Search ADS

Emmanoulopoulos

McHardy

I. M.

Papadakis

I. E.

2011

MNRAS

416

L94

10.1111/j.1745-3933.2011.01106.x

Crossref

Search ADS

Gandhi

2009

ApJL

697

L167

10.1088/0004-637X/697/2/L167

Crossref

Search ADS

Gaskell

C. M.

2004

ApJL

612

L21

10.1086/424565

Crossref

Search ADS

Gleissner

Wilms

Pottschmidt

Uttley

Nowak

M. A.

Staubert

2004

A&A

414

1091

10.1051/0004-6361:20031684

Crossref

Search ADS

Gradshteyn

I. S.

Ryzhik

I. M.

Jeffrey

Zwillinger

2007

Table of Integrals, Series, and Products

Elsevier Academic Press

New York

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Hawley

J. F.

Balbus

S. A.

1991

ApJ

376

223

10.1086/170271

Crossref

Search ADS

Hogg

J. D.

Reynolds

C. S.

2016

ApJ

826

10.3847/0004-637X/826/1/40

Crossref

Search ADS

Ingram

Done

2012

MNRAS

419

2369

10.1111/j.1365-2966.2011.19885.x

Crossref

Search ADS

Ingram

2016

337

385

10.1002/asna.201612318

Ingram

Done

2011

MNRAS

415

2323

10.1111/j.1365-2966.2011.18860.x

Crossref

Search ADS

Ingram

van der Klis

2013

MNRAS

434

1476

10.1093/mnras/stt1107

Crossref

Search ADS

Kotov

Churazov

Gilfanov

2001

MNRAS

327

799

10.1046/j.1365-8711.2001.04769.x

Crossref

Search ADS

Lynden-Bell

Pringle

J. E.

1974

MNRAS

168

603

10.1093/mnras/168.3.603

Crossref

Search ADS

Lyubarskii

Y. E.

1997

MNRAS

292

679

10.1093/mnras/292.3.679

Crossref

Search ADS

Lyutyj

V. M.

Oknyanskij

V. L.

1987

Astron. Z.

465

McHardy

I. M.

Arévalo

Uttley

Papadakis

I. E.

Summons

D. P.

Brinkmann

Page

M. J.

2007

MNRAS

382

985

10.1111/j.1365-2966.2007.12411.x

Crossref

Search ADS

Mummery

2023

MNRAS

518

1905

10.1093/mnras/stac2846

Crossref

Search ADS

Mushtukov

A. A.

Ingram

van der Klis

2018

MNRAS

474

2259

10.1093/mnras/stx2872

Crossref

Search ADS

Mushtukov

A. A.

Lipunova

G. V.

Ingram

Tsygankov

S. S.

Mönkkönen

van der Klis

2019

MNRAS

486

4061

10.1093/mnras/stz948

Crossref

Search ADS

Nolan

P. L.

et al. ,

1981

ApJ

246

494

10.1086/158949

Crossref

Search ADS

Nowak

M. A.

Vaughan

B. A.

Wilms

Dove

J. B.

Begelman

M. C.

1999

ApJ

510

874

10.1086/306610

Crossref

Search ADS

Priedhorsky

Garmire

G. P.

Rothschild

Boldt

Serlemitsos

Holt

1979

ApJ

233

350

10.1086/157396

Crossref

Search ADS

Scaringi

Körding

Uttley

Knigge

Groot

P. J.

Still

2012

MNRAS

421

2854

10.1111/j.1365-2966.2012.20512.x

Crossref

Search ADS

Scaringi

Maccarone

T. J.

Kording

Knigge

Vaughan

Marsh

T. R.

Aranzana

et al. ,

2015

Sci. Adv.

e1500686

Crossref

Search ADS

PubMed

Turner

S. G. D.

Reynolds

C. S.

2021

MNRAS

504

469

10.1093/mnras/stab875

Crossref

Search ADS

Uttley

Wilkinson

Cassatella

Wilms

Pottschmidt

Hanke

Böck

2011

MNRAS

414

L60

10.1111/j.1745-3933.2011.01056.x

Crossref

Search ADS

Uttley

McHardy

I. M.

2001

MNRAS

323

L26

10.1046/j.1365-8711.2001.04496.x

Crossref

Search ADS

Uttley

McHardy

I. M.

Vaughan

2005

MNRAS

359

345

10.1111/j.1365-2966.2005.08886.x

Crossref

Search ADS

Vaughan

Uttley

Pounds

K. A.

Nandra

Strohmayer

T. E.

2011

MNRAS

413

2489

10.1111/j.1365-2966.2011.18319.x

Crossref

Search ADS

Walton

D. J.

et al. ,

2013

ApJ

777

L23

10.1088/2041-8205/777/2/L23

Crossref

Search ADS

Zoghbi

Fabian

A. C.

Uttley

Miniutti

Gallo

L. C.

Reynolds

C. S.

Miller

J. M.

Ponti

2010

MNRAS

401

2419

10.1111/j.1365-2966.2009.15816.x

Crossref

Search ADS

APPENDIX A: THE EXACT FORM OF THE RELATIVISTIC FOURIER–GREEN’s FUNCTIONS

The solution of the mass accretion rate Fourier–Green’s integral has the following general form

(A1)

where

$$\begin{eqnarray} \beta \equiv (1 + i) \sqrt{\pi f}. \end{eqnarray}$$

(A2)

The Mummery (2023) Green’s function solution is fully described by (in units where G = c = w = 1)

(A3)

$$\begin{eqnarray} q(x) = x^{1/4} \sqrt{x^{-\alpha } g(x)} \exp \left({1 \over 2 x}\right) \left[1 - {2\over x}\right]^{5/4 - 3/8\alpha } , \end{eqnarray}$$

(A4)

and

$$\begin{eqnarray} p(x) = {x^{1/2} \exp \left(-{1/ x}\right) \over 1 - {2/ x}} , \end{eqnarray}$$

(A5)

where

$$\begin{eqnarray} x \equiv {2 r \over r_I}, \quad \alpha = {1\over 4\nu } , \end{eqnarray}$$

(A6)

and ₂F₁(a, b; c; z) is the hypergeometric function. For the purposes of taking the derivative in the above expression it will be helpful to note that (Mummery 2023)

$$\begin{eqnarray} g(x) = {1\over 2} \int _2^{x} x^{\prime \alpha - 1} \sqrt{1 - {2\over x^{\prime }}} \, {\rm d}x^{\prime } , \end{eqnarray}$$

(A7)

and so

$$\begin{eqnarray} {\partial g \over \partial x} = {1\over 2} x^{\alpha - 1} \sqrt{1 - {2\over x}} . \end{eqnarray}$$

(A8)

Expanding in full we have

$$\begin{eqnarray} {1\over 2} \widetilde{G}_{\dot{M}} = p(x) {\cal A}_\nu {\cal B}_\nu {\partial q \over \partial x} + p(x) q(x) {\cal A}_\nu {\partial B_\nu \over \partial x}, \end{eqnarray}$$

(A9)

where for notational ease we define

$$\begin{eqnarray} {\cal A}_\nu \equiv \left\lbrace \begin{array}{@{}l@{\quad }l@{}}& K_\nu , \quad x\lt x_0, \\ \\ & I_\nu , \quad x\gt x_0 , \end{array}\right. \end{eqnarray}$$

(A10)

and

$$\begin{eqnarray} {\cal B}_\nu \equiv \left\lbrace \begin{array}{@{}l@{\quad }l@{}}& I_\nu , \quad x\lt x_0, \\ \\ & K_\nu , \quad x\gt x_0 . \end{array}\right. \end{eqnarray}$$

(A11)

Therefore

$$\begin{eqnarray} {1\over 2 p(x) q(x) } \widetilde{G}_{\dot{M}} = {\cal A}_\nu {\cal B}_\nu {\partial \ln q \over \partial x} + {\cal A}_\nu {\partial B_\nu \over \partial x} . \end{eqnarray}$$

(A12)

The Bessel derivative is simplified by noting

$$\begin{eqnarray} {\partial B_\nu \over \partial x} = \beta {\partial g \over \partial x} {{\rm d} {\cal B}_\nu (z) \over {\rm d} z}, \quad z \equiv \beta g(x) . \end{eqnarray}$$

(A13)

The following identities

$$\begin{eqnarray} {{\rm d} \over {\rm d}z} I_l(z) = I_{l-1}(z) - {l \over z} I_l(z), \end{eqnarray}$$

(A14)

and

$$\begin{eqnarray} {{\rm d} \over {\rm d}z} K_l(z) = - K_{l-1}(z) - {l \over z} K_l(z) , \end{eqnarray}$$

(A15)

suffice to give the quantity |${{\rm d} {\cal B}_\nu / {\rm d} z}$|⁠. Combining

$$\begin{eqnarray} {\partial g \over \partial x} &= {1\over 2} x^{\alpha - 1} \sqrt{1 - {2\over x}}, \end{eqnarray}$$

(A16)

$$\begin{eqnarray} {\partial \ln q \over \partial x} &= {1 - 2\alpha \over 4x} - {1\over 2 x^2} + {10\alpha - 3 \over 4 \alpha x^2 } {1 \over 1 - 2/x} + {1 \over 2g} {\partial g\over \partial x}, \end{eqnarray}$$

(A17)

with

$$\begin{eqnarray} {{\rm d} {\cal B}_\nu \over {\rm d} z} = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}& + I_{\nu -1} - {\nu \over \beta g} I_\nu , \quad x\lt x_0, \\ \\ & -K_{\nu -1} - {\nu \over \beta g} K_\nu , \quad x\gt x_0 , \end{array}\right. \end{eqnarray}$$

(A18)

gives

$$\begin{eqnarray} \widetilde{G}_{\dot{M}} = 2 p(x) q(x) \left[ {\cal C}_\nu {\partial \ln q \over \partial x} + \beta {\partial g\over \partial x} {\cal C}_\nu ^{\prime } \right], \end{eqnarray}$$

(A19)

where

$$\begin{eqnarray} {\cal C}_\nu \equiv {\cal A}_\nu {\cal B}_\nu , \quad {\cal C}_\nu ^{\prime } \equiv {\cal A}_\nu {{\rm d} {\cal B}_\nu \over {\rm d} z} . \end{eqnarray}$$

(A20)

As a final step, we must now re-insert the dimensionful quantities M_d, M_BH, G, c, etc., so that these results may be fit to observational data. The physical constraint that

$$\begin{eqnarray} \lim _{f \rightarrow 0} \widetilde{G}_{\dot{M}} (x \lt x_0, f) \rightarrow -M_d , \end{eqnarray}$$

(A21)

where M_d is the mass content of the perturbation, simplifies this procedure greatly. Taking this limit, and using the fact that

$$\begin{eqnarray} \lim _{\beta \rightarrow 0} I_\nu (\beta g) = {1 \over \Gamma (\nu +1)} \left({\beta g \over 2}\right)^{+\nu } + {\cal O}\left((\beta g)^{2+\nu }\right) , \end{eqnarray}$$

(A22)

and

$$\begin{eqnarray} \lim _{\beta \rightarrow 0} K_\nu (\beta g) \!=\!{\Gamma (\nu) \over 2} \left({\beta g \over 2}\right)^{-\nu } \!+\! {\Gamma (-\nu) \over 2} \left({\beta g \over 2}\right)^{+\nu } + {\cal O}\left((\beta g)^{2-\nu }\right) , \\ \end{eqnarray}$$

(A23)

gives the normalization constant

$$\begin{eqnarray} N_1 = -{M_d \over 4 p(x)} \left[\nu \left({g(x_0) \over g(x)}\right)^\nu {1 \over \partial _x q(x) + \nu q(x)/g(x) } x^{\alpha - 1} \sqrt{1 - {2\over x}} \right] , \\ \end{eqnarray}$$

(A24)

which must be multiplied to the above Green’s function. Finally, the combination βg(x) within the Bessel functions must be dimensionless, and thus the substitution

$$\begin{eqnarray} \beta \rightarrow N_2 \beta , \end{eqnarray}$$

(A25)

must be made. The constant N₂ is given by

$$\begin{eqnarray} N_2 = {1 \over \sqrt{f_0}} = {2\sqrt{2} \nu } \left({GM_{\rm BH} r_I^3 \over w^2}\right)^{1/4} \sqrt{1-{2 \over x_0}} x_0^{-\mu /2} \exp \left(-{1\over x_0}\right), \\ \end{eqnarray}$$

(A26)

which defines a fiducial viscous frequency associated with the Green’s function, f₀.

A1 The low-frequency correction factor δ(x, x₀)

As discussed in Section 4 of this paper, at very low frequencies f/f₀ ≪ 10⁻⁶ the analytical solutions of this paper differ from the numerical Fourier-Green’s solutions of the relativistic disc equations for outward propagating modes. Mathematically this results from the fact that the solutions derived in Mummery (2023) are not exact, but are asymptotic ‘leading order’ solutions [see Mummery (2023) for a detailed discussion]. Physically care is required because at very late times the exact numerical and analytical solution begin to deviate, and an extremely small fraction of the initial disc mass is not accreted in these solutions. As a result

$$\begin{eqnarray} \lim _{f \rightarrow 0} \widetilde{G}_{\dot{M}}(x \gt x_0, f) \rightarrow \delta (x, x_0) M_d \ne 0 . \end{eqnarray}$$

(A27)

The discrepancy is small, as a result of the high accuracy of these analytical solutions (Fig. 7), and typically only effects extremely small frequencies f/f₀ ≪ 10⁻⁶ to a small degree

$$\begin{eqnarray} \delta \ll 10^{-4}. \end{eqnarray}$$

(A28)

The function δ(x, x₀) can be written in closed form, and is equal to

$$\begin{eqnarray} \delta (x, x_0) = \left({g(x_0) \over g(x)}\right)^{2\nu } \left[{2g(x) \partial _x \ln q(x) - \nu x^{\alpha -1}\sqrt{1-2/x} \over 2g(x) \partial _x \ln q(x) + \nu x^{\alpha -1}\sqrt{1-2/x} } .\right] \end{eqnarray}$$

(A29)

The discrepancy can then be removed by taking, for example

$$\begin{eqnarray} \widetilde{G}_{\dot{M}}(x \gt x_0, f) \rightarrow \widetilde{G}_{\dot{M}}(x \gt x_0, f) - \delta (x, x_0) \exp \left(-100 \left({f \over f_0}\right) \right), \\ \end{eqnarray}$$

(A30)

which removes any discrepancies for f ≪ f₀ but does not effect the moderate-to-high frequency behaviour of the solutions.

APPENDIX B: STRESSED NEWTONIAN FOURIER–GREEN’s FUNCTIONS

In a Newtonian disc system with a central stress (torque at r = 0), the Green’s function solution of the disc equations is

$$\begin{eqnarray} G(x, x_0, t) = {q(x) \over t} \exp \left({- g(x)^2 - g(x_0)^2 \over 4 t}\right) I_{-\nu } \left({g(x)g(x_0)\over 2 t} \right) , \\ \end{eqnarray}$$

(B1)

where it is important to note that the index on the Bessel function is now negative. This solution is in once again more usefully written as a Laplace-mode superposition, of the form (Gradshteyn et al. 2007)

$$\begin{eqnarray} G(x, x_0, t) = \int _0^\infty q(x) J_{-\nu }(\sqrt{s} g(x)) J_{-\nu }(\sqrt{s} g(x_0)) \exp ({-st}) \, {\rm d}s , \\ \end{eqnarray}$$

(B2)

and therefore every step of the derivation performed in Section 2 still holds. The Fourier–Green’s function for a stressed disc is then

$$\begin{eqnarray} \widetilde{G}(x, x_0, f) = 2 q(x) \left\lbrace \begin{array}{@{}l@{\quad }l@{}}& I_{-\nu }(\beta g(x)) K_{-\nu }(\beta g(x_0)) , \quad x\lt x_0, \\ \\ & I_{-\nu }(\beta g(x_0)) K_{-\nu }(\beta g(x)) , \quad x\gt x_0 . \end{array}\right. \end{eqnarray}$$

(B3)

Finally, the derivative with respect to x of this function must be taken, and then the result multiplied by p(x) ∝ x^1/2. We again find a remarkably simple result

$$\begin{eqnarray} \widetilde{G}_{\dot{M}} = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}+ A \beta x^{(1-\nu)/4\nu } K_{-\nu }\left(\beta \epsilon x_0^{1/4\nu }\right) I_{1-\nu }\left(\beta \epsilon x^{1/4\nu }\right)\!, \,\, x \lt x_0, \\ \\ - A \beta x^{(1-\nu)/4\nu } I_{-\nu }\left(\beta \epsilon x_0^{1/4\nu }\right) K_{1-\nu }\left(\beta \epsilon x^{1/4\nu }\right)\!, \,\, x \gt x_0 , \end{array}\right. \end{eqnarray}$$

(B4)

where (units where G = c = w = 1)

$$\begin{eqnarray} A = {x_0^{1/4} \epsilon M_d }, \quad \epsilon = 2\nu \sqrt{2x_0^\mu } , \quad \mu = {3-1/\nu \over 2} , \end{eqnarray}$$

(B5)

and we remind the reader

$$\begin{eqnarray} \beta \equiv (1 + i) \sqrt{\pi f}. \end{eqnarray}$$

(B6)

Note that for the case of an inner disc stress the low-frequency limits of the mass accretion rate Fourier–Green’s functions are effectively reversed from the non-stressed case, and we find

$$\begin{eqnarray} \lim _{f \rightarrow 0} \widetilde{G}_{\dot{M}}(x, x_0, f) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}f^A, \quad x \lt x_0, \quad A \gt 0, \\ {\rm const}, \quad x \gt x_0. \end{array}\right. \end{eqnarray}$$

(B7)

In other words, all of the disc material is eventually expelled to infinity, with none of the matter crossing the boundary at x = 0.

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.

Download all slides

Month:	Total Views:
May 2023	11
June 2023	52
July 2023	41
August 2023	34
September 2023	14
October 2023	18
November 2023	12
December 2023	20
January 2024	22
February 2024	18
March 2024	22
April 2024	15
May 2024	24
June 2024	25
July 2024	21
August 2024	19
September 2024	25
October 2024	16
November 2024	8
December 2024	15
January 2025	10
February 2025	7
March 2025	14
April 2025	16
May 2025	3

Article Contents

Extending the theory of propagating fluctuations: the first fully relativistic treatment and analytical Fourier–Green’s functions

ABSTRACT

1 INTRODUCTION

2 Thin disc Green’s functions in the Fourier domain: general solution and properties

2.1 The disc evolution equations

2.2 Solution of the Fourier integral

2.3 Asymptotic properties

2.3.1 High frequency Fourier modes

2.3.2 Low frequency Fourier modes

2.4 Discontinuity of the mass accretion rate Fourier–Green’s function at x = x0

2.5 The mass accretion rate power spectrum

2.5.1 High frequency power spectrum

2.5.2 Low frequency power spectrum

3 Newtonian Fourier–Green’s functions

4 Relativistic Fourier–Green’s functions

4.1 Relativistic Green’s functions in the time domain

4.2 Relativistic Green’s functions in the frequency domain

4.3 Differences between the Newtonian and relativistic Fourier–Green’s functions

5 OBSERVABLE QUANTITIES

5.1 The power density spectrum of a band

5.2 The cross spectrum between bands

5.3 The coherence between bands

6 THE BACK HOLE SPIN DEPENDENCE OF OBSERVABLE QUANTITIES

7 FUTURE EXTENSIONS TO THE MODEL

8 CONCLUSIONS

DATA ACCESSIBILITY STATEMENT

ACKNOWLEDGEMENTS

References

APPENDIX A: THE EXACT FORM OF THE RELATIVISTIC FOURIER–GREEN’s FUNCTIONS

A1 The low-frequency correction factor δ(x, x0)

APPENDIX B: STRESSED NEWTONIAN FOURIER–GREEN’s FUNCTIONS

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only

2.4 Discontinuity of the mass accretion rate Fourier–Green’s function at x = x₀

A1 The low-frequency correction factor δ(x, x₀)