https://orcid.org/0000-0003-2658-6559
Abstract
Accreting compact objects show variations in source flux over a broad range of time-scales and in all wavebands. The light curves typically show a lognormal distribution of flux and a linear relation between flux and rms. It has been demonstrated that an exponential transform of an underlying (and unobserved) Gaussian stochastic process provides a very good description of the light curves with these observed properties. Recently, a non-stationary power spectrum was observed on fast time-scales (∼days) in the active galaxy, IRAS 13224–3809, as well as a non-lognormal flux distribution and non-linear rms–flux relation. Here, we investigate the effects of piecewise non-stationary power spectra on the resultant flux distribution and rms–flux relation. We demonstrate that the simple ‘exponentiation’ model successfully reproduces the observed quantities, even when the light curves are non-stationary. We also demonstrate how non-lognormal flux distributions and rms–flux relations inconsistent with a linear model can be erroneously produced from poorly sampled PSDs. This is of particular importance for AGN surveys where very long baselines are required to sample the PSD down to low enough frequencies.
1 INTRODUCTION
One of the defining characteristics of accreting compact objects is their variation in source flux, observed in many wavebands and over a broad range of time-scales. The variability is a stochastic process which can be described using a power spectral density (PSD), typically showing broad-band noise power over several decades in frequency. The origin of the variability and origin of emission mechanisms are not necessarily the same thing. It is therefore crucial to understand the properties of the variability process if we want to (i) understand the underlying accretion physics, and (ii) use the variability to perform spectral timing studies to infer important properties of the accretion flow.
Another apparent universal feature of accreting objects is the linear relationship between the short-term rms amplitude and mean flux variations on longer time-scales. Known as the linear rms–flux relation (Uttley & McHardy 2001), this has been observed in active galactic nuclei (AGN; e.g. Vaughan et al. 2011), neutron star (NS) and black hole (BH) X-ray binaries (XRBs; e.g. Gleissner et al. 2004; Heil, Vaughan & Uttley 2012), ultraluminous X-ray sources (ULXs; Heil & Vaughan 2010; Hernández-García et al. 2015), as well as the fast optical variability from XRBs (Gandhi 2009), blazars (Edelson et al. 2013), accreting white dwarfs (Scaringi et al. 2012; Dobrotka & Ness 2015; Van de Sande, Scaringi & Knigge 2015) and young stellar objects (Scaringi et al. 2015). Crucially, a given source exhibits a linear relation over a broad range of time-scales
Models for the variability process often only explain the shape of the PSD. However, this is not the best tool to distinguish among models, as many types of models can produce identical PSD shapes. The lognormal flux distribution, rms–flux relation and the resulting non-linear variability appear to be more fundamental features of the underlying process. These characteristics rule out many types of models to describe the process, including additive shot-noise, self-organized criticality, and any model which consist of multiple additive independent varying regions. The process is required to be multiplicative in order to produce the observed variability characteristics.
Uttley, McHardy & Vaughan (2005; hereafter UMV05) proposed the simple ‘exponentiation’ model to explain the key observed light-curve properties in a single energy band. Here, an underlying Gaussian distributed light curve is exponentiated to produce a non-linear light curve with the same properties as the data: a lognormal distribution in fluxes as well as a linear rms–flux relation. Importantly, the model produces the desired linear rms–flux relation over all time-scales
It is important to distinguish here models for the variability process with models for the underlying physics producing the observed variability. The leading physical model is based on the inward propagation of random accretion rate fluctuations in the accretion flow (e.g. Lyubarskii 1997; Kotov, Churazov & Gilfanov 2001; King et al. 2004; Zdziarski 2005; Arévalo & Uttley 2006; Ingram & Done 2010; Kelly, Sobolewska & Siemiginowska 2011, Ingram & van der Klis 2013; Scaringi 2014). In this model, slow fluctuations in the local mass accretion rate at large radii propagate through the accretion flow and modulate faster fluctuations at smaller radii. This multiplicative model inherently produces variations coupled together over a broad range in time-scales. Numerical accretion disc simulations (e.g. Cowperthwaite & Reynolds 2014) including the effects of magnetohydrodyamics (MHD; Hogg & Reynolds 2016) are starting to reproduce the observed flux distribution and rms–flux relation.
Another important property of a stochastic process is the concept of stationarity. A process is said to be strongly stationary when its mean and variance tend to some well-defined value on long time-scales. In the case of accreting objects, this ‘long’ time-scale should at least correspond to the time-scale required to measure the low-frequency rollover in the PSD. For AGN, observing constraints mean only the red noise part of the PSD at high frequencies is typically observed. They are therefore considered only weakly stationary, as the low-frequency break in the PSD is typically not observed: this roll-over to zero amplitude provides a well-defined mean at some finite duration. The full range of light-curve behaviour is therefore not observed in a typical observation. A process is considered non-stationary when one or more of its moments do not satisfy these conditions. The statistics literature has a strict definition of stationarity, whereas time series can be non-stationary in an infinite number of ways.
In Alston et al. (2019; hereafter A19) we revealed non-stationary variability on very fast time-scales (≲days) in the X-ray light curves of the AGN, IRAS 13224–3809. The non-stationarity in IRAS 13224–3809 manifested itself as changes in the PSD shape and normalization with time, in addition to the changes in mean flux over any given segment the PSD is measured. With a central black hole mass in IRAS 13224–3809 of |$\hbox{$M_{\rm {BH}}$}\sim 10^{6} \hbox{$\, \mathrm{M}_{\odot }$}$|, this is equivalent to the broad-band noise PSD changing after just ∼10 s in a |$10 \hbox{$\, \mathrm{M}_{\odot }$}$| BH XRB.
In addition to the strong non-stationarity, IRAS 13224–3809 displays a non-lognormal distribution of source flux, as well as a non-linear rms–flux relation on all observed time-scales. The relation was well modelled as rms ∝ flux2/3, meaning the source is fractionally more variable at lower source flux, as was evident in the PSD. Given the short duration the deviations from the expected observables occur on, an important question to ask is whether the simple exponentiation model of UMV05 still work in the case of non-stationary light curves, or is some other more complex process model required to produce the observed quantities?
In this paper we explore the effects of non-stationarity light curves with the exponentiation model and use this to describe the observed variability properties of IRAS 13224–3809. We will also demonstrate the effects of observational constraints on the true PSD for the measured flux distributions and rms–flux relations. This has big implications when using the variability properties for understanding the accretion process, in particular from AGN time domain surveys.
2 RESULTS
2.1 The exponentiation model
We start by demonstrating the exponentiation model using some simple broad-band PSDs and inspecting the resulting flux distributions and rms–flux relations. Light curves are simulated following the approach described in detail in UMV05 and Uttley, McHardy & Vaughan (2017; hereafter UMV17). The model produces light curves x(t) with a linear rms–flux relation on all time-scales by taking the exponential of a linear and stationary aperiodic light curve. We used the Timmer & König (1995) method to produce a zero-mean Gaussian distributed stationary light curve, g(t). The model light curve is then generated by x(t) = exp[g(t)], which is non-linear by design. We simulate (Poisson noise free) light curves of length T = 1 × 108 s and time resolution dt = 1 s. The non-linear transformation increases the total variance of the light curves, so these are then scaled to the variance of the input PSDs (UMV05). The light curves were then rebinned to dt = 50 s to alleviate aliasing effects (e.g. Uttley, McHardy & Papadakis 2002; Alston, Vaughan & Uttley 2013).
![Demonstration of the exponentiation model for producing stationary light curves with the properties seen in accreting compact objects. Panel (a) shows model PSDs formed of identical broad (q = νc/σL = 0.2) Lorentzians with normalizations differing by a factor of 3. The resulting exponentiated light-curve flux distributions are shown in panel (b) and the log transformed fluxes in panel (c). The 0.2 − 0.1 mHz rms–flux relations for each model PSD are shown in panel (d). Panel (e) shows the zero mean Gaussian distributed stationary light curve g(t). Panel (f) shows the light curve x(t) = exp[g(t)] after being scaled to a mean of unity and variance of the model PSD.](https://oup.silverchair-cdn.com/oup/backfile/Content_public/Journal/mnras/485/1/10.1093_mnras_stz423/1/m_stz423fig1.jpeg?Expires=1750196190&Signature=1bHgumVK72RQ7qHkzSRQXCQd-TtoSE5Ub34zi1jgZ-SRVV0y~LL0F3funiWKtFhZzi8XYCNIt5ulcTCZEXV28YjO1OaaWT27E4dgEp800irsp24lY724QSxKHLhXqpc7rrgZxcq1IXXwstON6RsrwlpJb5bkxbTvnSBAlTXHqoKok9J8ULGQdCJmAbPvGvntlV1MFt3cbzNX83IcSwvKCgoJbb3vpYyM8PPbCtaTkaaEqBIEkKMljc1INDp9yqchAI0Vxks0qcMQefy8mVZVlepHXvyqQ~ATZOFzcn9Y7M~RTyQ17iX8RkapouJ4Cp8KHFI2nAnR7Pbwm55Gls2AIw__&Key-Pair-Id=APKAIE5G5CRDK6RD3PGA)
Demonstration of the exponentiation model for producing stationary light curves with the properties seen in accreting compact objects. Panel (a) shows model PSDs formed of identical broad (q = νc/σL = 0.2) Lorentzians with normalizations differing by a factor of 3. The resulting exponentiated light-curve flux distributions are shown in panel (b) and the log transformed fluxes in panel (c). The 0.2 − 0.1 mHz rms–flux relations for each model PSD are shown in panel (d). Panel (e) shows the zero mean Gaussian distributed stationary light curve g(t). Panel (f) shows the light curve x(t) = exp[g(t)] after being scaled to a mean of unity and variance of the model PSD.
We determine the rms–flux relation of the simulated light curves following the procedure described in Vaughan et al. (2011), Heil et al. (2012), and A19. The simulated light curves are split into segments with length Tseg = 105 s, resulting in an rms estimate for each segment. The data points are then sorted on mean flux and binned to have >50 estimates per bin. The 0.2 − 0.1 mHz rms–flux relations are shown in Fig. 1 panel (d), which as expected are well described by a linear model. As well as having a larger spread in flux the model light curve with higher variance also has a larger spread in rms values. As the mean rate of both model light curves is identical, the model with larger variability power is fractionally more variable, resulting in an rms–flux relation with a steeper gradient.
2.2 A non-stationary model for IRAS 13224–3809
We now extend this simple exponentiation model to describe the non-stationary light curves in IRAS 13224–3809 A19. This manifested itself as changes in the PSD shape and normalization with time, in addition to the changes in mean flux over any given segment the PSD is measured. For any given time series, one of the difficulties in measuring any changes in the PSD with time reliably, is the trade-off between variances of the estimates (i.e. error bars) and the frequency resolution of the data (i.e. the lowest frequency the PSD can be measured to). In IRAS 13224–3809, several neighbouring XMM–Newtonobservations were required to provide enough signal to noise to show a statistically significant change in the PSD. The light curves in IRAS 13224–3809 can be described as piecewise non-stationary; within the epoch used to measure a PSD it can be treated as stationary as the changes in the PSD are too small to be significantly measured. For simplicity we use only the 0.3–1.2 keV (soft) band due to the strong energy dependence of the variability. The 1.2–10.0 keV (hard) band also shows the same non-stationary behaviour (A19).
We use the PSD fits from A19 (fig. 10) where the data are sorted into three flux bins of equal duration; the PSDs within a flux bin are statistically equivalent. The PSDs consist of two broad Lorentzian components; one peaking at ∼10−5 Hz and one at ∼10−3 Hz. One of the important findings of A19 was that the low-frequency component rolls over to a slope α < 1 down to low frequencies; we are seeing the bulk of the variability power, which is all occurring on time-scales measurable with current observations. As the source flux decreases the low-frequency component moves to higher frequency (∼10−4 Hz). The centroid frequency, νc of the high-frequency component remains approximately constant.
The model PSDs with rms normalization are shown in Fig. 2 (a), where it can be seen how the data are fractionally more variable at lower source flux (A19). In a similar manner to the previous section, we simulate (Poisson noise free) light curves of length T = 1 × 107 s and time resolution dt = 1 s from each of the 3 PSDs. Following the exponentiation the light curves are rebinned to a time resolution of dt = 50 s. The model light curves were then scaled to the observed count rate of 1.0, 1.6, and 3.1 |${\rm cts\, s^{-1}}$| for the low, medium, and high-flux bins, respectively. The light curves are then scaled to the variance of the respective input PSDs. The total model light curve is then the concatenation of the piecewise stationary light curves, producing a piecewise non-stationary light curve.
The exponentiation model applied to the piecewise non-stationary PSDs observed in IRAS 13224–3809. Panel (a) shows the rms normalized model PSDs for the low (grey), medium (blue), and high (yellow) flux bins. Panel (b) shows the flux distributions of the model light curves, shown as a probability density (normalized to unity). The black curve is the total model light curve and the red curve is the IRAS 13224–3809 data from A19. Panel (c) shows the log transformed flux distributions. Panel (d) is the unbinned 0.2–0.1 mHz rms–flux relation for each model light curve. Panel (e) shows the same data after rebinning to 100 estimates per flux bin. The dotted lines are the best-fitting linear models. Panel (f) shows the rms–flux relation for the total model light curve. The red curve is the best-fitting model from A19 and the black dotted line is the best-fitting linear model.
Fig. 2(b) shows the flux distributions for the three model light curves as well as their sum, representing what is captured in the ‘real’ observations. As expected from the results of the previous section, the individual flux distributions are lognormal (or Gaussian when log transformed). When the sum of all three pieces is captured, this produces a clear deviation from a lognormal relation. This is due to the mean values of each distribution differing as well as the variances. Also shown in Fig. 2(b) is the observed flux distribution from A19, with error bars proportional to |$\sqrt{N}$| data points per histogram bin. The model light curves are a remarkably good fit to the data, with the skew to low fluxes being captured by the model. There is slight disagreement with the data, but this is to be expected as the exponentiation is a very simplified model (UMV05, UMV17) and the simulation is one realization of a stochastic process.
The rms–flux relation is determined in a similar manner to before. In Fig. 2 panel (d) we show the unbinned rms–flux relation for each piecewise model simulation. The spread in values for each piecewise model can be seen, as well as the overlap in estimates from the individual model; the flux binning to form the relation is performed on each piecewise model light curve only. Panel (e) shows the binned (>50 estimates) rms–flux estimates for each piecewise model light curve. As individual models are stationary, they produce a linear rms–flux relation following the exponentiation, with the gradient of each relation related to the variance of the model light curve.
If we treat the three model light curves as one large model, we can perform the flux sorting and binning on the model as a whole, representing what is typically done with real data. Fig. 2 panel (f) shows the rms–flux relation for the total model light curves with >100 estimates per bin. The resulting relation shows a clear deviation to a linear model, with a linear model ruled out with p-value <10−5. Overplotted in panel (f) is the best-fitting model with rms ∝ flux2/3 to the 0.2–0.1 mHz data from A19. We note here that the observed non-linear rms–flux relations at other frequency bands are also reproduced by the model. As before, the model is not a perfect fit to the data, but captures the important aspects of being fractionally more variable at lower flux and flattening at higher fluxes.
2.3 Weakly stationary power spectra
In this section we look at the expected flux distributions and rms–flux relations from poorly sampled PSDs. In previous sections we have shown that in addition to the PSD, any valid variability model must be non-linear, and explain the flux distribution and rms–flux relation of the observations (see also UMV05 and UMV17). In some situations, observational constraints mean that the high frequency part of PSD is only observed, with the high-frequency break often not detected. This is particularly true for AGN as the characteristic variability time-scales scale inversely with the mass of the accreting object, so are shifted to lower Fourier frequencies. The light curves can then be described as being only weakly stationary. This may then lead to not sampling the full range of the flux distribution, and falsely detecting inconsistencies with the expectations of the exponentiation model.
To represent ‘real’ observations we take (non-overlapping) subsections from x0(t) with durations T1 = 108 s: x1(t), and T2 = 107: x2(t). These represent observations which measure the PSD down to just below and at νb, respectively. Fig. 3(a) shows the model PSD and effective sampling for each of the model light curves. Panels (b, c and d) show the flux distributions for each model light curve, along with the best-fitting three-parameter lognormal distribution. The full light curve x0(t) is well described by a lognormal. The ‘real’ light curves x1(t) and x2(t) clearly show a poor fit to this model, displaying a multimodal distribution which cannot be captured by a simple lognormal. Also apparent in these plots is how the means of the two sampled light curves are different to the original.
The exponentiation model applied to weakly stationary PSDs, which do not break to zero power at low frequencies (see text for details). Panel (a) shows the rms normalized model PSD using a bending power-law function νb = 1 × 10−7 Hz. The full PSD provides the model light curve x0(t) with duration T0 = 1010 s (blue). From this we take (non-overlapping) subsections to represent a ‘real’ observations, x1(t) with T1 = 108 s (red) and x2(t) with T2 = 107 s (green). Panels (b, c, and d) show the flux distributions of the model light curves and best-fitting lognormal distribution. Panels (e, f, and g) show the rms–flux relation for each model light curve, calculated at frequencies above νb. The solid lines are the best-fitting linear relations to each model light curve.
Fig. 3 panels (e, f, and g) show the rms–flux relation for each model light curve, calculated at frequencies above νb. The relation is calculated as before with ≥100 estimates per flux bin. As the rms–flux relation requires many estimates to average we use a different frequency range – hence different segment length – for each model light curve. The rms is calculated in the range [νb, νNyq] Hz, [10−5, νNyq] Hz, and [10−4, νNyq] Hz, where νNyq = 1/2dt = 5 × 10−3 Hz is the Nyquist frequency, for x0(t), x1(t), and x2(t), respectively. The rms–flux relation for x0(t) shows excellent agreement with a linear model. There is some scatter, but this is expected for a stochastic process and using the simple exponentiation model (Vaughan et al. 2003). Despite being sampled from a light curve with an intrinsic linear rms–flux relation, x1(t) and x2(t) show a large amount of scatter and poor fit to a simple linear model.
We note here that these subsections are random selections of the original realization of the long light curve. Different subsections show a variety of flux distributions with positive or negative skewness and bimodality. A variety of rms–flux relations are also apparent, with the deviation from linear on average increasing as the light curve duration decreases. The simulations presented here do not include additive Poisson noise (caused by measurement noise in real observations). The effect of Poisson noise is to increase the width of the flux distribution and increase the scatter in the rms–flux relation. However, this effect is marginal when the power on the observed time-scales is much greater than the Poisson noise level, as is typically the case with real observations.
3 DISCUSSION
Any model for the variability process in accreting systems needs to be able to produce observed flux distribution and rms–flux relation. The simple exponentiation model is successful at describing the stationary light curves observed in accreting compact objects (UMV05, UMV17). For the case of the AGN, IRAS 13224–3809, we have shown that the observed flux distribution and rms–flux relation can be reproduced using the exponential model. The light curves are made up of piecewise stationary epochs which themselves have the typical properties. The result of the sum of several non-stationary pieces is to distort the resulting observables.
This is important for understanding the source behaviour, as it means the underlying process is similar to that of other accreting sources, i.e. it has its origin in the accretion flow. Some property of the standard process is changing with time to produce the observed light curve properties. The question still remains of what is the origin of the non-stationarity in IRAS 13224–3809. This is something which can be understood in MHD accretion disc simulations (e.g. Hogg & Reynolds 2016), with the non-stationarity providing an extra handle for physical changes in the system. To first order, the resulting flux distribution and rms–flux relation are not the result of some extrinsic modification of the standard variability process, for example variable absorption, as this is very different or non-existent in many other types of accreting compact objects.
The changes in IRAS 13224–3809 correspond to changes in the PSD on ∼10 s in a |$10 \hbox{$\, \mathrm{M}_{\odot }$}$| XRB, with the frequency range explored here equating to ∼1–100 Hz. The shape of the broad-band PSD means that IRAS 13224–3809 is a likely analogue of the transition state known as the very-high/intermediate state in BH XRBs (A19 and references therein). Changes in the PSD and light curve behaviour on short time-scales are seen in this transitional state (e.g. Belloni et al. 2000; Nespoli et al. 2003). In particular, quasi periodic oscillations are observed to appear and change frequency on time-scales of ∼5 s (Nespoli et al. 2003). These fast changes have been linked to radiation pressure driven instabilities in the accretion disc structure, expected when the accretion rate is a large fraction of the Eddington rate (e.g. Szuszkiewicz & Miller 1998; Nayakshin, Rappaport & Melia 2000). Understanding the observational properties of these non-stationary PSDs in XRBs may be possible with new instruments, such as NICER.
The model provides a nice description of the observed data in IRAS 13224–3809, but there are small differences. These are expected given the exponentiation approach is very simplified, resulting in a trivial frequency coupling and a time-symmetric model light curve. An improved approach is to include higher order measures such as the bicoherence (UMV05; Maccarone & Coppi 2002). In addition, a slight energy dependence to the PSD was found in A19, which is only crudely replicated in the model PSDs and mean count rates here: the non-stationary behaviour may be different for neighbouring energy bands which is only approximated here using a broad energy band. This is the focus of a further paper in a series on IRAS 13224–3809.
The non-linear rms–flux relation observed in IRAS 13224–3809 is an example of Simpson’s paradox (Blyth 1972; Wagner 1982). This is the phenomenon where a particular trend appears in separate groups, but disappears or even reverses when the groups are combined. In the case of IRAS 13224–3809, the linear rms–flux trend is present over a short enough interval where the PSD is statistically stationary. When estimates from a piecewise non-stationary light curve are combined, the individual linear trends are lost in the flux sorting, resulting in a distorted overall relation.
Emmanoulopoulos, McHardy & Papadakis (2013) presented an algorithm for producing simulated light curves with the same properties as the light curves for many types of astrophysical sources, including accreting compact objects. In the special case where the desired output distribution of the simulation is lognormal, the simulated light curves also have a linear rms–flux relation, without the need for the exponentiation step. This arises from the fact the Timmer & König (1995) method only produces Gaussian distributed light curves, so the lognormality needs to be imprinted. However, like this method, the Emmanoulopoulos et al. (2013) algorithm also only generates stationary light curves for a given input PSD. In order to produce the same results for IRAS 13224–3809 in Section 2.2 a piecewise non-stationary PSD approach will also have to be used.
Smith et al. (2018) recently reported non-lognormal flux distributions in a sample of ∼20 Kepler optical AGN, with a variety of skewness and (bi)modality reported. This data also showed rms–flux relations inconsistent with a simple linear model. The authors inferred that some different underlying variability process or accretion mode was at play in these sources. However, the high-frequency PSD bend was poorly detected in these data, or at the lower end of the frequency window. The results in Section 2.3 show how the measured flux distributions and rms–flux relations cannot look as expected, even when those light curves are taken from a large light curve which is lognormally distributed and have a linear rms–flux relation by design. Caution must therefore be taken when making inferences about the underlying process when the PSD is only poorly sampled (i.e. is only weakly stationary). This is something which can be trivially assessed using simulations similar to those presented here, but specific to the PSD and count rate of the data in question. We note here that the optical emission from AGN can arise from multiple emission mechanisms and multiple emission regions, potentially providing a source of non-stationarity in the light curves.
4 CONCLUSIONS
We have demonstrated how the simple exponentiation model can successfully produce observed flux distribution and rms–flux relation. The atypical results observed in the AGN, IRAS 13224–3809, are well described by the model when the PSD is piecewise non-stationary. Given the similarities in the variability process observed across accreting compact objects, it is likely that similar effects will be observed in further sources. Understanding the cause of the non-stationarity will then provide further insights into the accretion process. There will indeed be some deviations from the model due to e.g. changes in source state, long term accretion rates, absorption, and obscuration events. These can of course all be accounted for in the simulation process.
With longer baseline AGN variability studies becoming available it will be tempting to use these in a similar manner to what has classically been achieved with X-ray timing studies (see e.g. Vaughan et al. 2016). Before any claims of disagreements with simple model can be made, it is important to check your results with simulated light curves with the correct intrinsic properties to see what the flux distribution and rms–flux relation should look like for the observed PSD.
ACKNOWLEDGEMENTS
WNA acknowledges discussions on accretion disc simulations with Chris Reynolds. WNA acknowledges support from the European Research Council through Advanced Grant 340442, ‘Feedback’. This paper is based on observations obtained with XMM–Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA).
