Hot-spot model for accretion disc variability as random process
Key Words.:Accretion, accretion-discs – Black hole physics – Galaxies: active – X-rays: binaries
Aims:Theory of random processes provides an attractive mathematical tool to describe the fluctuating signal from accreting sources, such as active galactic nuclei and Galactic black holes observed in X-rays. These objects exhibit featureless variability on different timescales, probably originating from an accretion disc.
Methods:We study the basic features of the power spectra in terms of a general framework, which permits semi-analytical determination of the power spectral density (PSD) of the resulting light curve. We consider the expected signal generated by an ensemble of spots randomly created on the accretion disc surface. Spot generation is governed by Poisson or by Hawkes processes. The latter one represents an avalanche mechanism and seems to be suggested by the observed form of the power spectrum. We include general relativity effects shaping the signal on its propagation to a distant observer.
Results:We analyse the PSD of a spotted disc light curve and show the accuracy of our semi-analytical approach by comparing the obtained PSD with the results of Monte Carlo simulations. The asymptotic slopes of PSD are at low frequencies and they drop to at high frequencies, usually with a single frequency break. More complex two-peak solutions also occur. The amplitude of the peaks and their frequency difference depend on the inherent timescales of the model, i.e., the intrinsic lifetime of the spots and the typical duration of avalanches.
Conclusions:At intermediate frequencies, the intrinsic PSD is influenced by the individual light curve profile as well as by the type of the underlying process. However, even in cases when two Lorentzians seem to dominate the PSD, it does not necessarily imply that two single oscillation mechanisms operate simultaneously. Instead, it may well be the manifestation of the avalanche mechanism. The main advantage of our approach is an insight in the model functioning and the fast evaluation of the PSD.
It is widely accepted that massive black holes with accretion discs reside in cores of active galactic nuclei, where most activity originates and X-rays are produced (e.g., Blandford et al. 1990). The observed light curves, , show irregular, featureless fluctuations with a very complex behaviour, practically at every studied frequency (Gaskell et al. 2006). Variability has been traditionally analysed by the Fourier method (Feigelson & Babu 1992). Remarkably, a number of similarities appear between the properties of massive black holes in galactic nuclei and those in X-ray binaries, suggesting that some kind of universal rescaling operates according to central masses of these systems (Mirabel & Rodríguez 1998). This concerns also the X-ray power spectra (e.g., Markowitz et al. 2003; McHardy et al. 2006).
Light curves can be characterised by an appropriate estimator of the source variability which, in the mathematical sense, is a functional: . We accept the idea that the signal resulting from a spotted accretion disc is intrinsically stochastic, likely originating from turbulence. Hence, is a random value. It can be a number (for example, the ‘rms’ characteristic), function of a single variable (for example, the power-spectral density – PSD) or of many variables (e.g., poly-spectra, rms–flux relation, etc). The appropriate choice depends on the type of information we seek and the quality of data available. The PSD is a traditional and widely utilised method to examine variable signals, and the AGN light curves are no exception. A typical signal can be represented by a broad band PSD with the tendency towards flattening at low frequency (Lawrence et al. 1987; McHardy & Czerny 1987; Lawrence & Papadakis 1993; Mushotzky et al. 1993; Uttley et al. 2002).
There is an ongoing debate about the overall shape of the PSD and the occurrence of the break frequency or, possibly, two break frequencies at which the slope of the PSD can change (Nowak et al. 1999; Markowitz et al. 2003). In the case of a widely-studied Seyfert galaxy, MCG–6-30-15, McHardy et al. (2005) have closely examined the slope of PSD, namely its bending, with RXTE and XMM-Newton data. It is worth noticing that the accurate fits to the X-ray sources seem to exhibit a multi-Lorentzian structure rather than a simple power law. The same is true for the best studied example, the Seyfert 1 galaxy Akn 564 (McHardy et al. 2007).
It was proposed (Abramowicz et al. 1991; Zhang & Bao 1991; Wiita et al. 1992) that hot spots contribute to the variability of the AGN variability, or that they could even be the dominant process shaping the variability pattern. These spots should occur on the disc surface following its intermittent irradiation by localised coronal flares (Galeev et al. 1979; Merloni & Fabian 2001; Czerny et al. 2004). Here, the “spots” represent a somewhat generic designation for non-axisymmetric features evolving on the disc surface in connection with flares. They share the bulk orbital motion with the underlying disc. The observed signal is thus modulated by relativistic effects as photons propagate towards a distant observer.
Various schemes have been discussed in which the fluctuations of the disc emissivity at different points of space and time are mutually interconnected in some way. In particular, the avalanche model (Poutanen & Fabian 1999; Życki 2002; Życki & Niedźwiecki 2005) seems to be physically substantiated within the framework of magnetically-triggered flares and spots. It is also a promising model capable to reproduce, for example, the broken power-law PSD profiles. Notice, however, that other promising ideas were proposed (e.g., Mineshige et al. 1994; Lyubarskii 1997), provoking the question of whether a common mathematical basis could reflect the entire range of models and provide us with general constraints, independent of (largely unknown) model details.
We add to this model by applying the method of random point processes (Cox & Miller 1965). Interestingly enough, a rather formal approach can provide useful analytical formulae defining the basic form of the expected power spectrum. Apart from this practical aspect, we suggest that the concept discussed here offers much better insight into various influences that shape the expected form the power spectrum. These are very attractive features especially with respect to avalanche models, which may have different flavours, typically with a vast parameter space, thus proving very demanding to examine in a systematic manner.
Even more important is that the adopted formalism provides a very general tool and allows for a broader perspective on different mechanisms of variability (Pecháček & Karas 2007). We develop the idea in a systematic way and give the explicit correspondence between our approach and some of the above-mentioned and widely-known scenarios (Abramowicz et al. 1991; Poutanen & Fabian 1999). This description provides semi-analytical solutions, convenient to search through a broad parameter space. Our results can help to identify how the intrinsic properties of individual flares and the relativistic effects influence the overall PSD. In particular, we can identify those situations in which a doubly-broken power law occurs.
We consider stochastic processes (e.g., Feller 1971; Gardiner 1994) in the framework appropriate for modelling the accretion disc variability. In particular, in Sect. 2 we consider a simple model of a spotted accretion disc constrained by the following three assumptions about the creation and evolution of spots: (i) each spot is described by its time and place of birth (, and ) in the plain of the disc; (ii) every new occurrence starts instantaneously; afterwards, the emissivity decays gradually to zero (the total radiated energy is of course finite); and (iii) the intrinsic emissivity is fully determined by a finite set of parameters which form a vector, , defining the light curve profile. Later on, we will consider the modulation of the intrinsic emission by the orbital motion and relativistic lensing. The disc itself has a passive role in our considerations; we will treat it as a geometrically thin, optically thick layer lying in the equatorial plane.
Because of the apparently random nature of the variability, we adopt a stochastic model in which the creation of spots is governed by a random process. The assumption that spots are mutually statistically independent seems to be a reasonable (first) approximation, however, we find that we do need to introduce some kind of relationship between them. This connection is discussed in Sect. 3. The statistical dependence among spots can be introduced in several ways. In Sect. 4, we explore in detail the specific models of interrelated spots using the formalism of Hawkes-type processes. Conclusions are summarised in Sect. 5. Finally, in the Appendix we provide some mathematical prerequisites, which the reader may find useful to understand the general background of the paper, and we also summarise the mathematical notation.
2 Models of stochastic variability
2.1 Orbiting spots and relativistic effects
We will apply our investigations to models where the signal is produced by point-like orbiting spots (circular Keplerian motion along the azimuthal direction). The intrinsic emission, produced in the local co-orbiting frame of the spot, is influenced by the Doppler effect and gravitational lensing, which cannot be ignored at typical distances of several units or tens of gravitational radii. Photons emitted at different moments and positions experience different light-travel time on the way towards the observer, so the observed timing properties should reflect this specific modulation. We adopt the Schwarzschild metric for the gravitational field and employ the method of transfer function (Pecháček et al. 2005, 2006) to describe the light amplification (or dilution); is inclination angle of the observer ( deg corresponds the edge-on view of the disc plane). The periodical modulation of the observed signal is included in the transfer function . An implicit relation holds for the phase,
where the effect of time delays is taken into account. The modulation by is superimposed onto that due to intrinsic timescales present in the signal from spots. Then, for the final flux measured by a distant observer, we write .
In the case of an infinitesimal surface element with intrinsically constant and isotropic emissivity , orbiting with Keplerian orbital frequency , the flux measured by a distant observer varies just as changes along the orbit.
We remind the reader that the mass of central black holes in galactic nuclei is in the range of – Mass of the accretion disc is at least three orders of magnitude smaller than the black hole mass, so we neglect it in our calculations (the accretion disc self-gravity may be important for its intrinsic structure, but the direct gravitational effect on light is quite small; Karas et al. 1995). Hence, the gravitational field can indeed be described by a vacuum black-hole spacetime (Misner et al. 1973). We use geometrical units with . Transformation to physical scales can be achieved when the mass of the central black hole is specified because Keplerian frequency scales inversely with the gravitational radius. The gravitational radius of a massive black hole is pc, and the corresponding characteristic time-scale is sec, where the mass . All lengths and times can be made dimensionless by expressing them in units of , so they can be easily scaled to different masses. For example, for the Keplerian orbital period we obtain , where the radius is expressed in units of and is given in seconds.
Let us note that the intrinsic timescales of the spot evolution and of avalanches (both timescales will be discussed below) need not to be directly connected with the Keplerian orbital period. This internal freedom of the model can help to bring the predicted frequency of the breaks of the PSD profile in agreement with the data.
2.2 Model driven by a general point process
Now we will describe the process of the creation of spots from the statistics point of view. Let us consider a signal of the form
The underlying process consists of a sequence of events which, in general, can be either mutually independent, or there can be some statistical dependence among them. Naturally, the latter case will be more complicated and interesting. In Eq. (2), is the light curve profile of a single event; and are the time offsets (determined by the moment of the ignition of the spot) and the initial time delay (which is an arbitrary but fixed value); is the Heaviside function; and is a non-negative function of variables, and . Hereafter, we omit the explicit dependency on .
Quantities , , , and are random values. The vector determines the duration and shape of each event ( is time of ignition; parameter determines the initial phase of the periodical modulation of the –th event; and is the corresponding initial time-offset). These assumptions bring the formulation of the problem close to the framework studied by Brémaud et al. (2002; 2005). We will calculate the power spectrum of this process directly from Eq. (100) in the Appendix.
We remark that for the amplitudes of individual events we assume the identical values (at each given radius). This restriction is imposed only for the sake of definiteness of our examples; the formalism can deal with a distribution of amplitudes. Indeed, we do not impose any serious constraint on the model because the information about the level of the fluctuating signal can be adjusted by setting the frequency of the events (Lehto 1989). A simple demonstration of this concept is shown in Fig. 1. This plot illustrates how the model light curve arises from the elementary components. Naturally, we can approach such decomposition from another angle, by investigating how the total light curve can be expressed in terms of some basic profile. It is important to realise that, for the purposes of our present paper, light curves are of secondary importance. Instead, our calculations allow us to proceed from the distribution of the onset times and the characteristics of individual flares directly to the power spectral density, which stands as the primary characteristic of the source signal.
Equation (2) represents a very general class of random processes. By applying the Fourier transform, we find
where denotes the convolution operation. In Eq. (2), we sum together a set of all events (this infinite sum could in general pose problems with convergence, however, as we will see later, we can restrict the sum to a finite set of events without any loss of generality). The Fourier transform of a single event, , is then
Periodical function can be now expanded in a series, , and the expanded form substituted in the incomplete Fourier transform of .
Knowing the incomplete Fourier transform of , we can calculate its squared absolute value and perform the averaging over all realisations of the process. This can be simplified by assuming that every single event quickly decays. In principle, between and the process is influenced by all events ignited during the whole interval , however, because of the fast decay of the events, the interval can be restricted to , where is a sufficiently large positive constant. In other words, every realisation of the process on can be described by a set of points in -dimensional space with .
Values of the initial time delay and phase are functions of initial position of each spot, i.e.
Fourier transform of the resulting signal is then
where is the Fourier transform of the event light curve, corrected for the time delay. Every realisation of this process is completely determined by the set of points .
Defining the function
we can write
Due to Campbell’s theorem (105),
where is density of the second-order moment measure corresponding to the random point process of , is a Cartesian product of four sets representing the domains of definitions, i.e.
2.3 Model driven by the Poisson process
To start with a simple example, we assume mutually independent events with uniformly distributed ignition times. In other words, in this subsection we assume that there is no relationship among different spots – neither in their position nor in the time of birth (spots are statistically-independent). The intensity and the second-order measure are
where is the mean rate of events. Other parameters are treated as independent marks with common distribution . The second-order measure is
The result of integration (9) can be simplified for . The task reduces to evaluation of two limits (see Pecháček 2008, for details). After somewhat lengthy calculations we obtain a general formula for the power spectrum:
Here, we remark that the general relativity effects are included in the Fourier coefficients . Knowing them in advance (i.e., pre-calculating the sufficient number of the coefficients that are needed to achieve the desired accuracy) helps us to produce the PSD efficiently. But we will start by neglecting these relativistic effects, so that we can clearly reveal the impact of the intrinsic timescales of individual spots and the form of the avalanche process.
2.3.1 Example 1: Markov chain
Let us consider box-shaped events with exponentially-distributed life-times, i.e.,
where is the characteristic function ( for , otherwise); is the probability density; and is a constant (we will set for simplicity). This again represents a process of the type (2), but at the same time it can be considered as a continuous-time Markov chain with discrete states (Cox & Miller 1965). The process PSD is then given by Eq. (14) with
In this example, by setting the transfer function we completely “switch-off” the periodical modulation and the relativistic effects.
Coefficients are given by the relation
which for the constant function leads to , . Fourier transform of the profile function is
Substituting into the general formula (14) we find
Notice that the theory of Markov chains is usually formulated in another way, independent of our previous calculations (see Cox & Miller 1965). One can thus obtain the PSD of the Markov process by two different approaches and verify the results. Such comparison gives the same result, as expected.
3 Introducing a relationship among spots
The assumption that spots are statistically independent seems to be a reasonable (first) approximation. However, the actual ignition times and spot parameters should probably depend on the history of a real system. The statistical dependence among spots can be introduced in several ways. In this section we discuss different models where an existing spot gives, with a certain probability, birth to new spots. In this way a single spot at the beginning can trigger a whole avalanche of its offsprings. This avalanche can be in principle of arbitrary length, although, to obtain an infinitely long avalanche with non-diverging rate of new spots one would have to fine-tune the parameters. In order to avoid the unlikely fine-tuning and to obtain a stationary process, we will assume the occurrence of many spontaneous spots distributed by the Poisson process that keeps triggering new avalanches of finite duration.
The first example can be called “Chinese process”. By definition, an existing spot gives birth to exactly one new spot with probability . In other words, every event produces at most one offspring. The spot of the th generation is always ignited later than the spot of the th generation. As mentioned above, spontaneous spots arise randomly, according to Poissonian process. In the simplest version of this model, delays between the parent spot and its lineal descendant are random values obeying the probability density .
More generally, every spot can deliver new spots, where is a random value with Poisson distribution and the mean . We describe this situation in terms of (i) standard Hawkes’ process, and (ii) the pulse avalanche model. The temporal distribution of new spots is now governed either by the function of the Hawkes’ process (Hawkes 1971), or in the case of avalanches (Poutanen & Fabian 1999). Spots of different generations can appear at the same time.
The difference between the Chinese process and the latter two processes is schematically sketched in Fig. 2. Mathematically, all three examples belong to the class of cluster processes.
3.1 The cluster processes
Point processes are characterised by the generating functional, , which is defined by its action (Daley & Vere-Jones 2003)
where is a function satisfying the condition .
The functional satisfies various relations which can be derived in close analogy with the theory of generating functions of random variables. For our purposes it will be useful to expand into a series in terms of factorial measures,
The cluster processes consist of two point processes. The first one is connected with the counting measure and defines the central points of the clusters. The second process spreads new points around the central point according to the random measure . The complete counting measure is then .
Let be the generating functional of a cluster with the center at ,
If the cluster center process is Poissonian, the latter formula simplifies to
where is the intensity measure of the center process. One can expand the functional (24) in terms of the factorial measure of the cluster,
For a stationary process, the intensity measure stays constant, , and gives the density of the centers. All factorial moments are shift-invariant in the sense
where is density of . As a result of the shift invariance, density of the first-order moment must be also constant. Moreover, we can always choose in Eq. (28) equal to one of , so the function depends on only variables. We define the reduced factorial moments,
and from Eq. (27) it follows that
Stationarity of the process implies that the second-order measure density depends only on the difference of its arguments,
where is an even function.
The second-order measure of a marked cluster process is
almost identical with that of the Poissonian process. There is only one additional term associated with the function . The resulting spectrum is given by a somewhat lengthy, but still explicit formula. We find stationary cluster processes to be particularly promising as a general scheme, encompassing a broad range of models as special cases. The PSD is
The reduced quadratic factorial moment appears in the formulae for power spectra of cluster processes expressed in terms of . This is directly related to the two-dimensional Fourier transform of the quadratic factorial measure as
We will discuss the form of , the expression for , and the resulting PSD in different situations. But before that, we still need to show how the model properties can be detailed in term of marks.
3.2 Marks as a way to specify the model properties
Until now the variability patterns have been restricted only by very general properties of the assumed process. This means that the model is kept in a very general form. However, formulae (14) and (33) are too general for any practical use. Their main value is after defining special cases. Then these formulae can be readily applied to derive the analytical form of the PSD. Such special cases are conveniently defined by means of marks. We discuss possible choices of the mark distribution, .
We can simplify the situation by assuming axial symmetry. Therefore, all statistical properties should depend only on the radius (the azimuthal part of is constant). The distribution of marks has now the form
where is the reduced probability density of the remaining parameters. To illustrate this case more clearly, we will now examine the phenomenological model of Abramowicz et al. (1991) and Zhang & Bao (1991), which can be considered as representation of the spotted accretion disc.
3.3 Example 2: a spotted disc
Let be an average number of spots at radius , each of them shining with the average intensity for average duration . Let us further assume that these characteristics scale with the radius as power laws:
(here, s and s are constants). This setup falls within the category of models described by Eq. (2), with the profile function and the mark distribution . The first two relations (38) prescribe the conditional marginals of , i.e.,
where the integration limits can depend explicitly on . The third equation determines the amplitude of the profile function,
The dependency on is not determined uniquely. In order to obtain an explicit form of we have to go beyond the model of Abramowicz et al. (1991) by assuming the distribution of ,
The normalisation constant then equals
and for we find
is consistent with both equations. Constants , and are bound by the relation
Because and are positive, we can write
Therefore, by substituting back to Eq. (41), we obtain
where the definition domain is a set
We can set as a specific example. This value of the power-law index is special in the sense that neither short nor long timescales dominate the PSD, as we see from Eq. (42):
By substituting back to Eq. (41) we obtain
Knowing the concrete form of the distribution of marks, we perform the integration over . As does not explicitly depend on azimuthal angle, the integration is simplified. Denoting
we rewrite the PSD formula (14) in the final form
We find the coefficients by direct evaluation,
with . We note that the term (55) corresponds to the effect of delay amplification in terminology of Dovčiak et al. (2008). It influences the observed signal from a source moving (i.e., orbiting) near a black hole.
4 Results for the model PSD
4.1 Model driven by the Chinese process
Let us denote the probability that an existing spot generates a new one, and the probability that a family of spots consists of exactly members. The value obeys the geometrical distribution, .
We interpret probability density of the delay between successive spots as a mean number of first-generation spots that occur at the ignition time , where is a moment of ignition of the parent spot. Analogically, is the mean number of second-generation spots. For a sequence of spots, we obtain the intensity measure
where is the th convolutionary power of . We can write in the form
Convolution of two functions can be calculated via the Fourier image. Defining we get
From this we find
The meaning of integral (59) is the average number of spot offsprings in the whole avalanche. The meaning of the last integral is the average duration of the avalanche.
Calculation of the quadratic measure is a less intuitive procedure. We start from the generating functional (21) of the process,
Substituting in the expansion (22), we obtain the Fourier image
Equation (64) can be cast in the form
4.2 Model driven by the Hawkes process
Hawkes’ process (Hawkes 1971; Brémaud & Massoulié 2002) is more complicated than the previous example because it consists of two types of events. First, new spots are generated by the Poisson process operating with intensity (let us denote the moment of ignition). Second, an existing spot can give birth to new one at a later time, , according to the Poisson process with varying intensity .111Mathematically, a very similar model has been employed to describe the propagation of diseases through a population (Daley & Vere-Jones 2003). In this context, the function is called “infectivity”, for obvious reasons. We will adopt the same name for it, although the medical connotation is irrelevant here.
The mean number of events is
In analogy with Eq. (60) we define the characteristic time of avalanche . It can be proven that is related to the characteristic time of the infectivity :
For a stationary process, the first-order moment density is constant. Averaging the relation (66) we find
The meaning of is the mean number of the offsprings. Clearly, it satisfies the normalisation .
The generating functional of the cluster of the Hawkes process fulfils the integral equation,
Substituting and expanding both sides into the series (25) we find
To complete the calculation we solve the integral equation (70) for . Because this is a linear convolutional integral equation, it can be solved efficiently by using the Fourier transform:
For the Fourier transform of the quadratic factorial measure we find
Comparing this equation describing the Hawkes process PSD with the corresponding Eq. (65) for the case of the Chinese process, we reveal a subtle difference between the two mechanisms. It turns out that the high-frequency limit is identical for both of them, however, the difference grows as one proceeds towards the low-frequency end of the PSD domain.