The rapid Xray variability of NGC 4051
Abstract
We present an analysis of the high frequency Xray variability of NGC 4051 ( ) based on a series of XMMNewton observations taken in 2009 with a total exposure of ks (EPIC pn). These data reveal the form of the power spectrum over frequencies from Hz, below the previously detected power spectral break, to Hz, above the frequency of the innermost stable circular orbit (ISCO) around the black hole ( Hz, depending on the black hole spin parameter ). This is equivalent to probing frequencies kHz in a stellar mass ( ) black hole binary system. The power spectrum is a featureless power law over the region of the expected ISCO frequency, suggesting no strong enhancement or change in the variability at the fastest orbital period in the system. Despite the huge amplitude of the flux variations between the observations (peaktopeak factor of ) the power spectrum appears to be stationary in shape, and varies in amplitude at all observed frequencies following the previously established linear rmsflux relation. The rmsflux relation is offset in flux by a small and energydependent amount. The simplest interpretation of the offset is in terms of a very soft spectral component that is practically constant (compared to the primary source of variability). One possible origin for this emission is a circumnuclear shock energised by a radiatively driven outflow from the central regions, and emitting via inverseCompton scattering of the central engine’s opticalUV continuum (as inferred from a separate analysis of the energy spectrum). A comparison with the power spectrum of a long XMMNewton observation taken in 2001 gives only weak evidence for nonstationarity in power spectral shape or amplitude. Despite being among the most precisely estimated power spectra for any active galaxy, we find no strong evidence for quasiperiodic oscillations (QPOs) and determine an upper limit on the strength of a plausible QPO of per cent rms in the Hz range, and per cent in the Hz range. We compare these results to the known properties of accreting stellar mass black holes in Xray binaries, with the aim further of developing a ‘black hole unification’ scheme.
keywords:
galaxies: active – galaxies: individual: NGC 4051 – galaxies: Seyfert – Xrays: galaxies1 Introduction
If strong gravity dominates the dynamics of the inner accretion flows around black holes then an elementary consequence is scale invariance: many important aspects of accretion onto supermassive black holes ( ) in Active Galactic Nuclei (AGN) should be fundamentally the same as for stellar mass black holes ( ) in Galactic Black Hole systems (GBHs). (See e.g. Shakura & Sunyaev 1976 and Mushotzky et al. 1993.) Over the past few years several similarities have been observed in the Xray variability of nearby AGN and GBHs, supporting the idea of black hole unification (Fender et al., 2006). In particular, the power spectrum (often power spectral density, PSD), and interband Xray time lags appear to have similar frequency dependence and relative amplitudes in AGN and GBHs but with characteristic frequencies scaled (to first order) inversely with mass, (see e.g. Lawrence et al., 1987; Hayashida et al., 1998; Edelson & Nandra, 1999; Uttley et al., 2002; Vaughan et al., 2003; Markowitz et al., 2003; MHardy et al., 2004; Done & Gierliński, 2005; MHardy et al., 2006; MHardy et al., 2007; Gierliński et al., 2008; Arévalo et al., 2008). Another property of the Xray variability common to both types of system is a linear rmsflux relation, a scaling of the rms amplitude of variability with average source flux that appears to hold over a very wide range of timescales (Uttley & MHardy, 2001; Vaughan et al., 2003; Gaskell, 2004; Gleissner et al., 2004; Uttley et al., 2005).
Arguably among the most important recent breakthroughs in the study of stellar mass black holes was the discovery of quasiperiodic oscillations (QPOs) at high frequencies ( Hz) in GBHs (see Remillard et al., 1999; Strohmayer, 2001; Remillard et al., 2002, 2003; van der Klis, 2006; Remillard & McClintock, 2006). These are the fastest variations observed from GBHs and their high frequencies, close to the Keplerian frequency of the innermost stable circular orbit (ISCO; ), indicate an origin in the central regions of the accretion flow. High frequency QPOs (HF QPOs) therefore carry information on the strongly curved spacetime close to the event horizon and should, in principle, constrain the fundamental parameters of black holes: mass and spin. As a product of strong gravity, HF QPOs should also be present in AGN, but until recently there were no robust detections, due mainly to the insufficient length of previous XMMNewton observations (see Vaughan & Uttley, 2005, 2006). Recently a QPO was detected in the Seyfert galaxy RE J1034+396 (Gierliński et al., 2008). This is a major discovery that adds further support to the general idea that AGN are scaledup GBHs, but this one detection provides limited diagnostic power because it is so far unique and the other properties of the source are not well understood (Middleton et al., 2009; Middleton & Done, 2010)
In parallel with these Xray advances, studies of transient radio jets from GBHs have led to the development of a scheme that unifies jet production with accretion state (Fender et al., 2004). Relativistic jets must be launched from the region dominated by strongfield gravity, and so scale invariance implies this accretion state/jet unification scheme should extend to AGN (Heinz & Sunyaev, 2003), a hypothesis supported by observations (Merloni et al., 2003; Falcke et al., 2004; Körding et al., 2006). From these results a new paradigm is emerging in which accretion mode, Xray spectrum, highfrequency timing properties and jet production for both GBH “states” and AGN “types” may be unified into a single framework for the activity cycles of accreting black holes.
Here we discuss the results obtained from a series of XMMNewton observations of the bright, highlyvariable Seyfert 1 galaxy NGC 4051. These observations were designed to provide some of the best constraints to date on the highfrequency Xray variability properties of any AGN and improve our understanding of black hole unification generally. This particular target was chosen because of a combination of factors. NGC 4051 (redshift , distance Mpc; Russell, 2002) has a relatively low and welldetermined black hole mass, (Denney et al., 2009, 2010), displays strong, persistent Xray flux and spectral variability (see e.g. Lawrence et al., 1987; Papadakis & Lawrence, 1995; MHardy et al., 2004; Uttley et al., 2004; Ponti et al., 2006; Terashima et al., 2009), and a soft Xray spectrum rich in emission and absorption features (Collinge et al., 2001; Pounds et al., 2004; Ogle et al., 2004; Steenbrugge et al., 2009).
2 Observations and data analysis
NGC 4051 was observed by XMMNewton times over days during MayJune 2009, and previously during 2001 and 2002 (see table 1). The total duration of the useful data from the 2009 campaign is ks, giving EPIC pn source counts (using PATTERN ) or EPIC source counts (using PATTERN pn and MOS combined).
The raw data were processed from Observation Data Files (ODFs) following standard procedures using the XMMNewton Science Analysis System (SAS v10.0.2). The EPIC data were processed using the standard SAS processing chains to produce calibrated event lists. For each observation, source events were extracted from these using a arcsec circular region centred on the the target, and background events were extracted from a rectangular region on the same chip but not overlapping with the source region. (In order to obtain a precise background estimate the background region was times larger than the source region.) Examination of the background time series showed the background to be relatively low and stable throughout most of the observations, with the exception of the final few ks of each observation where the background rate increased as the spacecraft approached the radiation belts at perigee (each of the observations occurred towards the end of a spacecraft revolution). The EPIC data were taken using the small window mode of each camera and the medium blocking filter to reduce pileup and optical loading effects, respectively. Nevertheless, during bright intervals NGC 4051 is sufficiently bright to cause moderate pileup, especially in the MOS, which can affect power spectrum estimation (Tomsick et al., 2004). In order to mitigate any distorting effects all results are based on pn data PATTERNS only, but where possible checked for consistency with the MOS data. Response matrices were generated using RMFGEN v1.55.2 and ancillary responses were generated with ARFGEN v.1.77.4.
Obs ID.  Obs. date  Duration  

ID  rev.  (YYYYMMDD)  (ks) 
0109141401  0263  20010516  121958 
0157560101  0541  20021122  51866 
0606320101  1721  20090503  45717 
0606320201  1722  20090505  45645 
0606320301  1724  20090509  45548 
0606320401  1725  20090511  45447 
0606321301  1727  20090515  32644 
0606321401  1728  20090517  42433 
0606321501  1729  20090519  41813 
0606321601  1730  20090521  41936 
0606321701  1733  20090527  44919 
0606321801  1734  20090529  43726 
0606321901  1736  20090602  44946 
0606322001  1737  20090604  39756 
0606322101  1739  20090608  43545 
0606322201  1740  20090610  44453 
0606322301  1743  20090616  42717 
Time series were extracted from the filtered events for the source and background regions using a range of time bin sizes and energy ranges. These were all corrected for exposure losses, as catalogued in the Good Time Information (GTI) header contained in each event file, except where less than per cent of a bin was exposed, in which case the data were calculated by linear interpolation from the nearest ‘good’ data either side and appropriate Poisson noise was added (this was needed for only per cent of the total exposure time from the pn, for s time bins). The outputs of this process were regularly sampled and uninterrupted, exposurecorrected time series for source and background. The background time series was then subtracted from the corresponding source time series (after scaling by the ratio of the ‘good’ extraction region areas). Figure 1 shows the resulting time series from the separate 2009 observations.
3 Power spectrum
The power spectrum of the 2009 data was estimated using standard methods (van der Klis, 1989). In particular, time series were extracted with binning s and divided into uninterrupted ks intervals. For each interval a periodogram was computed in units of relative power (Miyamoto et al., 1991; van der Klis, 1997; Vaughan et al., 2003), the Poisson noise level^{1}^{1}1The expected power density due to measurement uncertainties in this case is , where is the sampling rate and is the mean square error. See Vaughan et al. (2003). was subtracted, and then the periodgrams were averaged^{2}^{2}2 We have confirmed that the power spectral estimation and modelling was not significantly biased by sampling distortions – commonly known as ‘aliasing’ and ‘red noise leakage,’ see van der Klis (1989); Uttley et al. (2002); Vaughan et al. (2003). Aliasing is minimal for contiguously binned data (van der Klis, 1989), like those used here, and does not occur for the Poisson noise that dominates the highest frequencies probed. Leakage could in principle have been a problem, resulting from the finite duration of each observation, but in practice was not significant because the power spectrum was sufficiently flat () on timescales comparable to the observation length. This was confirmed using an analysis discussed in section 4, footnote 4.. Figure 2 shows the resulting power spectrum estimate after rebinning in logarithmic frequency intervals. The model fitting was performed using XSPEC v12.6.0 (Arnaud, 1996) using as the fit statistic, all data were binned to ensure that Gaussnormal statistics applied, and errors quoted correspond to per cent confidence limits (e.g. using a criterion for one parameter) unless otherwise stated.
The previous analysis of the power spectrum of NGC 4051 by MHardy et al. (2004) found a smoothly bending power law continuum provided an adequate description of the data
(1) 
with indices of at low frequencies and above Hz. This model provided an excellent fit to the 2009 power spectrum, as shown in Figure 2, with the lower frequency index fixed at (since it is much better constrained by the RXTE data; MHardy et al., 2004). The bestfitting values for the free parameters of the model were and Hz. The fit statistic was with degrees of freedom (dof), giving a value of . Allowing to vary provided no improvement in the fit () and the index itself was poorly constrained (with a per cent confidence interval of ), and restricting to vary only within the range (the confidence interval determined by MHardy et al. (2004) from the RXTE data) provided very little change in the confidence interval for compared to the fixed model. There is no obvious indication, in either the overall fit statistic or the residual plot, of additional power spectral components. We discuss this point in more detail in Section 5.
A sharply broken power law model also provided a good fit, with for dof (), slightly worse than the smoothly bending model. The index below the break was again fixed to , and the bestfitting break frequency was Hz, and a high frequency index of . These parameters, and the overall fit quality, are very similar to those of the bending power law model, indicating that the data are rather insensitive to the detailed shape of the change from flat to steep index. By contrast, the power spectrum is less well explained in terms of an exponentially cutoff power law, which gives for dof (), but only by allowing which is inconsistent with the longterm RXTE results (MHardy et al., 2004). Assuming gave an unacceptable fit ( for dof; ).
The model fitting was repeated using a Bayesian scheme by treating the likelihood function as and assigning simple prior distributions to the parameters. See Vaughan (2010), and references therein, for a brief introduction to these ideas. Uniform priors were assigned to and the overall normalisation, and to (equivalent to a Jeffreys’ prior), although in practice using uniform or Jeffreys’ prior made no significant difference to the results. As expected, the parameter values at the posterior mode were almost identical to those giving the bestfit using the classical method. But having the problem restated in Bayesian terms allowed for an exploration of the posterior distribution of the parameters, with the parameters specifying the continuum and the observed data, using a Markov Chain Monte Carlo (MCMC) scheme to randomly draw sets of parameters from the posterior^{3}^{3}3Specifically, five chains of length were generated (after removing an equal length “burnin” period), using a MetropolisHastings algorithm with a multivariate Cauchy proposal distribution. Convergence of the chains was confirmed using the statistic (Gelman & Rubin, 1992; Gelman et al., 2004). The results from each chain were then combined to produce a large sample of posterior draws.. The marginal density for the bend frequency, , as computed using the MCMC data, is shown in Figure 3; the per cent credible interval on was Hz.
The limited durations of the individual 2009 observations makes it difficult to directly probe frequencies lower than Hz. But for completeness the power spectrum was also estimated using longer uninterrupted intervals ( ks duration, of which there were in 2009, one for each observation). However, the estimates of power density at the lowest frequencies were not normally distributed in this case, due to the small number of estimates contributing to the average (Papadakis & Lawrence, 1993). These data were not included in the power spectral fitting, but the lowest frequency points are shown in Figure 2 and are clearly consistent with an extrapolation of the bestfitting model.
As a test for a possible high frequency cutoff in the power spectrum the bending power law model was modified by including an exponential cutoff (the highecut model in XSPEC). This provided only a very small improvement in the overall fit statistic ( for additional free parameters) with only poorly constrained parameters (cutoff frequency and folding frequency ). The power spectrum estimated from the high flux data (see next section) has slightly higher signaltonoise above mHz, but again shows no significant difference in the quality of the fit between the bending power law with and without an exponential cutoff at high frequencies.
3.1 Energy dependence of the power spectrum
The dependence of the power spectrum on photon energy was investigated by dividing the data into three broad energy bands – –, – and – – and estimating and fitting the power spectrum for each of these. A simultaneous fit to all these spectra using the same continuum model (equation 1) gave a rather poor fit ( with dof, ) with the parameters tied between the three spectra. But allowing to be different for each energy band, while keeping tied between them, gave a much better fit ( with dof, ); the improvement is significant in an test, with , indicating the high frequency power spectrum is energy dependent. The slopes estimated for the low, medium and high energy bands were , and , respectively. The flattening of the power spectrum at higher energies is very similar to that detected by MHardy et al. (2004) from the 2001 data, and previously observed in several other Seyfert galaxies (e.g. Nandra & Papadakis, 2001; Vaughan & Fabian, 2003). Allowing the bend frequency to be different for each energy band, but with tied between them, provided a considerably worse fit ( with dof) than the energydependent model. Furthermore, allowing both and to differ between energy bands did not improve significantly on the variable fit ( with dof). In this case the bestfitting values for the low and medium energy bands were very similar while the hard band gave a somewhat lower value, but with a much larger confidence interval. We conclude there is strong evidence for a change in the high frequency slope with energy, but the data do not resolve any difference in the location of the bend.
4 rmsflux relation
The observations shown in Figure 1 clearly show very different mean levels and amplitudes of variability, with noticeably stronger variations occurring preferentially during brighter intervals (e.g. rev1730), and the lowest flux intervals being extremely stable (e.g. rev1739). This behaviour is exactly as expected for sources that obey the rmsflux relation (Uttley & MHardy, 2001; Uttley et al., 2005).
In order to test this explicitly and in detail the power spectrum was estimated in different flux bins. As before, the time series were divided into uninterrupted ks intervals, the periodograms computed and the contribution of the Poisson noise subtracted, but using absolute units of power (i.e. without normalising to the squared mean rate, see Vaughan et al., 2003). The intervals were sorted according to their mean flux, and the periodograms were then averaged in flux bins to provide estimates of the power spectrum as a function of mean flux. Figure 4 shows the power spectral estimates when the data were divided into three flux bins. There is clearly a large difference in the overall amplitude (in absolute units) but little change in shape, as revealed by the ratio of bright and faint power spectra. Direct fitting of the high and low flux power spectra gives the same result: the bestfitting values for and estimated for the high flux data fall within the per cent confidence intervals derived from the low flux data. The same effect was seen in Cygnus X1 by Uttley & MHardy (2001) (see their figure ), and is consistent with there being a single rmsflux relation independent of the frequencies used to compute the rms.
The detailed rmsflux relation was then examined explicitly by dividing the data into ks uninterrupted intervals, computing a periodogram for each (again, in absolute units), subtracting the Poisson noise from each, and then averaging the results in several flux bins^{4}^{4}4 We have checked that the power spectra estimated using short time series intervals (e.g. ks) are not significantly biased by leakage (see Uttley et al., 2002, and references therein) in two different ways. First we compared the average power spectrum estimated from the real data using ks, ks and ks intervals and confirming they are all consistent. Secondly, by simulating data (using the bestfitting bending power law model) with the same sampling pattern as the real data, and estimating the bias between the power spectral estimate and the model. (). The mHz rms, was then estimated for each flux bin by integrating the corresponding power spectrum estimate over this frequency range. The resulting rmsflux data are shown in Figure 5. (The uncertainties on the rms estimates were calculated based on the variance of the sum of the averaged periodograms, which are themselves calculated based on the chisquare distribution of periodogram ordinates. The method is explained more fully in Heil, Vaughan & Uttley in prep.)
The dependence of rms on flux is clearly very close to linear over the full flux range, as previously noted in this source and accreting black holes generally (e.g. Uttley & MHardy, 2001; Uttley et al., 2005, Heil, Vaughan & Uttley in prep.). The bestfitting linear function – , where is the gradient and is the offset on the flux axis – gave a rather high fit statistic ( for dof) but the residuals show little structure. The bestfitting linear function has a negative offset in the rms axis, equivalent to a positive offset on the flux axis. The energy spectrum of this offset was computed by dividing the pn data into different energy bands ( logarithmically spaced) and for each one calculating the rmsflux data and fitting with a linear model. The spectrum of the flux offset is shown in Figure 6, along with some representative pn spectra from bright (rev1730), intermediate (rev1743) and faint (rev1739) flux observations. The spectra were “fluxed” for display purposes, that is, they have been normalised to flux density units by dividing out the effective area curve (see Appendix A for details). The spectrum of the offsets from the rmsflux relation is clearly soft: it is significantly positive in the lower energy bands, decreasing until consistent with zero intercept above keV. Most interestingly, it has a very similar shape and strength to the faint flux spectrum below keV.
5 A search for quasiperiodic oscillations
The power spectral fitting discussed above provided no clear evidence for narrow features, i.e. QPOs. Such features are usually described in terms of Lorentzian “lines” seen in the power spectrum in addition to the continuum. In order to assess the evidence for a QPO more rigorously we defined two competing hypotheses: the null hypothesis of a continuum power spectrum, and the alternative hypothesis which includes an additional Lorentzian component.
(2) 
(3) 
where the continuum model is defined above (equation 1) and the QPO will be described by a Lorentzian profile
(4) 
The parameters of the continuum are and the QPO parameters are , where is the centroid frequency, is the width parameter (, where is the Full Width at Half Maximum), and is the amplitude parameter (equal to the total rms in the limit of high ). See Pottschmidt et al. (2003) for details of this parameterisation. In the present case we have assumed two different values, namely and , which are fairly typical values for high frequency QPOs in GBHs (e.g. Remillard et al., 2002, 2003). Using these two representative values rather than allowing to be a free parameter made the QPO search and upper limit calculations (discussed later) much more efficient.
One of the simplest ways to select between these two competing models is to assess the improvement in the fit statistic upon the addition of the QPO, i.e. the difference between the minimum for and . For a given dataset a large reduction (improvement) in the fit statistic between and is taken to be evidence in favour of a QPO. For the 2009 data, and the full keV energy band, the reduction in statistic was only (assuming , and a very similar value assuming ), with bestfitting parameters mHz and per cent. The modest improvement in the fit falls short of our threshold for a significant detection, as discussed below. The QPO search was repeated for the three energy bands discussed above, with very similar results. The – keV band gave the largest improvement in the fit upon the addition of a QPO ( assuming ), but this still falls short of a detection in the test.
The lack of detection may itself be an interesting result if the observation was sufficiently sensitive, and to assess this requires a calculation of the upper limit on the strength of a QPO. The first step is to define a precise detection procedure and then calibrate its statistical ‘power’ as a function of QPO strength. That would allow us to determine the weakest QPOs that would be detected with reasonable probability.
5.1 Detection procedure
The search for QPOs in binned power spectral data is no different from the search for “lines” in any other one dimensional spectrum, and suffers from the same statistical and computational challenges. See Freeman et al. (1999), Protassov et al. (2002), Park et al. (2008) for elaboration of the issues.
The difference in the minimum chisquare statistic for each model was used as a test statistic. For a dataset we compute
(5) 
The nuisance parameters – specifying the continuum and the QPO – are eliminated by minimisation, i.e. . This test statistic is closely related to the Likelihood Ratio Test (LRT), as has been discussed by Protassov et al. (2002). A large value of the test statistic (i.e. a large improvement in the fit upon adding a QPO to the model) is taken as evidence for a QPO.
The critical value of for detection, , was calculated by a Monte Carlo method as discussed in Appendix B. In the present case we found for (this is for , using gave which is practically the same). In other words, assuming the continuum model to be true (with parameters represented by the posterior distribution derived above from the real data), and no QPO, a reduction in the fit statistic of will occur with probability (and have a posterior predictive value ).
5.2 Upper limit procedure
In order to define an upper limit on the strength of possible QPOs we use the definition of an upper limit in terms of the power^{5}^{5}5Here we follow the notation used by Kashyap et al. (2010) and use to denote the power of the test, i.e. the probability of meeting the detection criterion assuming the effect is real (e.g. there is a QPO present). of a detection procedure, as discussed very clearly in the recent paper by Kashyap et al. (2010). A QPO of a given frequency and strength has a probability of detection using the above method (for fixed ).
We used a Monte Carlo method, as discussed in Appendix B, to estimate the detection probability for a range of different QPO locations and strengths ( and values). The estimate is simply the fraction of simulations (for each and ) that gave a significant detection defined by . Contours of constant on the  plane define the upper limit on the strength of a QPO as a function of frequency . Figure 7 shows the contours for , and assuming (solid curves) and (dotted curves). The middle of each set of curves represents the weakest QPO that could be present (at each frequency) and have a probability of being detected. This limit extends below per cent for Hz, and below per cent for Hz. The shape and amplitude of this curve is similar to that obtained using the simpler confidence contour method used previously by Vaughan & Uttley (2005).
6 Long term stationarity
The power spectral estimation and fitting described in section 3 was repeated for the 2001 data. Figure 8 shows a comparison of the 2001 and 2009 power spectra which are in broad agreement in shape and overall (fractional) amplitude. The bestfitting bending power law model ( for dof, ) had a slightly higher bend frequency ( Hz with a per cent confidence interval of Hz, and a credible interval of Hz) and highfrequency index () compared to the 2009 data. The marginal distribution of from both datasets are shown in Fig 3. These appear to suggest there was a slight change in the power spectral shape between the two sets of observations. It is for this reason that we have not combined the 2001 and 2009 data to produce a composite power spectrum.
7 Discussion and conclusions
We have presented an analysis of the rapid Xray flux variability of the lowmass Seyfert 1 galaxy NGC 4051 covering the range to Hz, based on a series of XMMNewton observations taken during 2009. This is among the best Xray power spectra for any active galaxy, and the large range of fluxes sampled during the observations has allowed for an examination of the fluxdependence of the variability. Below we discuss the results and, where possible, relate them to the known behaviour of GBHs (for which we use a fiducial mass of ).
7.1 Power spectral shape in relation to other sources
The power spectrum of the variations is well described by a bending power law model that bends from an index of to around Hz. This is fairly consistent with the MHardy et al. (2006) relation
(6) 
where is the break timescale^{6}^{6}6 Note that MHardy et al. (2006) defined as the break frequency estimated by fitting power spectra with a sharply breaking power law model. This difference is of little consequence to the present discussion as we obtained very similar estimates for the characteristic timescale using smoothly bending and sharply broken models (section 3). in days ( Hz), , and erg s is the bolometric luminosity. The average keV Xray flux during the 2009 observations was erg s cm which gives erg s at an assumed distance of Mpc (Russell, 2002), and (assuming a ; Elvis et al., 1994). The black hole mass is thought to be (Denney et al., 2009), therefore . Together these give a predicted bend frequency of Hz ( day), within a factor of the observed value. In fact the estimated from the 2009 data is closer to the value expected from the MHardy et al. (2006) relation than that estimated by MHardy et al. (2004), when calculated using the revised black hole mass estimate (Denney et al., 2009, 2010).
The high frequency power spectrum above the bend rarely provides strong, unambiguous signatures of the accretion state of GBH systems, which tend to differ more at lower frequencies where features such as low frequency QPOs, flattopped noise and broad Lorentzians occur (see e.g. van der Klis, 2006; McClintock & Remillard, 2006; Belloni, 2010). MHardy et al. (2004) suggested NGC 4051 may be analogous to a soft state (particularly like that seen in Cygnus X1) on the basis of the low frequency noise (a spectrum extending over two decades or more), the high frequency slope () and the value of (which tends to be higher in the soft state, making the ratio smaller).
For the 2009 data we found a very similar same high frequency shape and therefore is consistent with this identification. Ignoring any dependence and simply scaling frequencies inversely with mass (), the bend frequency for NGC 4051 correponds to Hz in a typical GBH. This is close to the Hz estimated from a soft state observation of Cygnus X1 by MHardy et al. (2004).
The state analogy is supported by radio observations. At the highest available resolution the source is resolved into a core and two oppositely directed sources separated from the core by arcsec (Ulvestad & Wilson, 1984; Kukula et al., 1995; Ho & Ulvestad, 2001; Giroletti & Panessa, 2009). This structure is similar to the core and hotspots commonly seen in FanaroffRiley II radio galaxies, although no radio jet has been detected in NGC 4051. Jones et al. (2010) recently used radio monitoring observations and found no evidence for strong variability of the core flux. Despite the presence of these radio sources, Jones et al. (2010) showed that NGC 4051 falls below the ‘fundamental plane’ connecting the radio and Xray luminosities to black hole mass in radio loud AGN and hardstate GBHs (Merloni et al., 2003; Falcke et al., 2004; Körding et al., 2006), i.e. it is less radioloud than would be expected for a hardstate or radioloud system (see Fig. 14 of Jones et al., 2010). This is consistent with the ‘soft state’ interpretation of NGC 4051 mentioned above.
The XMMNewton data reveal the power spectrum of NGC 4051 up to frequencies as high as Hz (see Fig 2 and 4), a factor above the bend frequency . Assuming a black hole mass as above, the gravitational radius is km, and the period of the innermost stable circular orbit (ISCO) should be s, depending on the dimensionless spin parameter, (see e.g. van der Klis, 2006, and references therein). The power spectrum of NGC 4051 above shows no indication of deviating from a power law up to the highest frequencies (but see section 7.3), and in particular shows no signature of the ISCO period, which might be expected to produce an excess of variability power due to Doppler boosting of “hot spots” (Revnivtsev et al., 2000; Syunyaev, 1973). There is no sign of the additional power spectral break seen at very high frequencies in the hard state of Cygnus X1 by Revnivtsev et al. (2000)^{7}^{7}7Above the break the power spectrum of NGC 4051 is sufficiently steep, , that the integrated power in the power spectra of and both converge at high frequencies. The rate of change of the luminosity may be limited by the radiative efficiency of accretion (Fabian, 1979). The power spectrum is therefore not required to become steeper at even higher frequencies.. Indeed, assuming linear scaling of timescales with black hole mass, we detect power at frequencies corresponding to kHz for a GBH, far above the fastest variations yet detected from these sources (Revnivtsev et al., 2000; Strohmayer, 2001).
7.2 Stationarity and the existence of “states”
The location of the power spectrum bend frequency is lower than that reported by MHardy et al. (2004) from their analysis of the 2001 XMMNewton data. The per cent confidence intervals from these two analyses – mHz and Hz – do not overlap, although some of this may be due to the interval reported by MHardy et al. (2004) being an underestimate (see Mueller & Madejski, 2009). The difference between the two results may in part be due to differences in the data reduction and power spectrum estimation techniques; the reanalysis of the 2001 data presented here (section 6) gives a lower estimate and shows the difference between the two observations to be only moderately statistically significant (compare also the Bayesian posterior distributions of Fig 3). If this is a real effect, it is different from the expected scaling of with accretion rate in the central region implied by the MHardy et al. (2006) relation, since the Xray luminosity in the month leading up to the 2009 observation (based on the RXTE monitoring) was per cent higher than in the month leading up to the 2001 observations, but the bend frequency was per cent lower.
The separation of years between the two observations corresponds to a separation of only minutes for a GBH; on these timescales GBH power spectra are approximately stationary except during very rapid state transitions. By the same linear scaling each of the ks XMMNewton observations is equivalent to s, and the day span of the observations covers the equivalent of only s for a typical GBH. Unless NGC 4051 is drastically more nonstationary that most GBHs we might therefore expect all the XMMNewton observations to be representative of almost exactly the same state and variability properties.
During both the 2001 and 2009 observations the variability does, to a very good approximation, follow a single rmsflux relation over the full flux range (see Figure 5) and the shape of the power spectrum remains roughly the same at different flux levels (see Figure 4). The lack of variability during low flux intervals (e.g. rev1739) and the extreme variability during high flux intervals (rev1730) are opposite extremes of this relation. These apparently different behaviours are therefore not the result of distinct variability states^{8}^{8}8 Not in the sense normally applied to GBHs where different states are often characterised by distinct power spectral shapes and amplitudes as well as different energy spectra. See van der Klis (2006); Remillard & McClintock (2006); McClintock & Remillard (2006)., despite their very different energy spectra (see Figure 6). (See also Uttley et al. 2003.) The variations in energy spectral shape are also continuous, as discussed by Vaughan et al. (2010, in prep.).
7.3 Absence of QPOs
We found no strong evidence (using a test) for QPOs (additional, narrow components) in the 2009 power spectrum, and we place limits of per cent rms on the strength of a QPO in the Hz range, and per cent over the Hz range. The Hz bandpass over which these data are sensitive to weak QPOs corresponds to Hz for a GBH, where HF QPOs are often found, and the upper limit of per cent on the QPO strength is comparable to the strength of typical HF QPOs (Remillard et al., 2003). Such HF QPOs been detected in GBHs only in intermediate states (particularly the soft intermediate state; Belloni, 2010), but are not always present (at least not at detectable levels), and when detected may be strongly energydependent. The lack of detection in NGC 4051 therefore does not constitute strong evidence against an analogous state in this AGN.
Despite not being formally significant in an test, the bestfitting narrow Lorentzian QPO component has a frequency of Hz and rms per cent. This frequency is very close to that found previously by Vaughan & Uttley (2005) from fitting the power spectrum of the 2001 data. This may be a coincidence in the random sampling fluctuations around a smooth continuum power spectrum. But it is around this frequency that the signaltonoise in the data is highest, and so model fitting will be most sensitive to small ( per cent) deviations from the smooth power law model. The high frequency power law may be only an approximation to a more complex continuum power spectrum, the structure of which is on the edge of detectability in these data.
It is well known that HF QPOs in GBHs can be strongly energy dependent, and frequently stronger at higher energies, and in fact we did find the higher energy band (– keV) gave slightly stronger evidence for a possible QPO. But this could also indicate that this energy band contains a stronger contribution from a second, variable emission component with a different power spectrum from that which dominates the softer energy bands.
7.4 A quasiconstant emission component
The flux offsets in the energyresolved rmsflux relations, , reveals the spectrum of a quasiconstant component (QC) of the emission. A simple model comprising a single power law () modified by Galactic absorption (Elvis et al., 1989) gave a very good fit, with for dof. A Comptonising plasma model (compTT) with optical depth and seed photon energy eV, which has an approximately power law shape, also gave a good fit ( for dof). The simplest interpretation is that this represents an emission component that does not vary at all during the 2009 observations; emission that remains if the highly variable emission component(s) drop to zero flux. We note that the Chandra data taken during an extended low flux period in 2001 and the XMMNewton observation taken during a low flux period in 2002 were previously interpreted in terms of an unresolved ( pc), nonvarying soft spectral component and a strongly variable, harder spectral component (Uttley et al., 2003, 2004). However, the QC emission is not required to be absolutely constant over the observations, it is only required to have a far lower (absolute) amplitude (at mHz) and be weakly correlated (or uncorrelated) with the highly variable component that dominates brighter flux spectra (e.g. rev1730) and drives the rmsflux relation.
The spectral shape and absolute flux level of this component are similar to the lowest flux spectrum taken during 2009 (rev1739) below keV, as shown in Figure 6. This strongly suggests that most of the soft Xray continuum emission observed during low flux periods is due to the soft QC component. This may be the tail of the spectrum of inverseCompton scattering of opticalUV continuum photons in shockheated gas where the ionised outflow collides with the interstellar medium of the host galaxy, as recently predicted to occur in AGN by King (2010). Indeed, the velocityionisation structure of the outflow as resolved by the RGS data provide further evidence to support this scenario (Pounds & Vaughan, 2011).
Acknowledgements
We thank the referee, Iossif Papadakis, for a thorough and constructive referee’s report. PU is supported by an STFC Advanced Fellowship, and funding from the European Community’s Seventh Framework Programme (FP7/20072013) under grant agreement number ITN 215215 ”Black Hole Universe”. This research has made use of NASA’s Astrophysics Data System Bibliographic Services, and the NASA/IPAC Extragalactic Data base (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This paper is based on observations obtained with XMMNewton, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA).
References
 Arévalo et al. (2008) Arévalo P., MHardy I. M., Summons D. P., 2008, MNRAS, 388, 211
 Arnaud (1996) Arnaud K. A., 1996, in G. H. Jacoby & J. Barnes ed., Astronomical Data Analysis Software and Systems V Vol. 101 of Astronomical Society of the Pacific Conference Series, XSPEC: The First Ten Years. p. 17
 Belloni (2010) Belloni T. M., 2010, in T. Belloni ed., Lecture Notes in Physics, Berlin Springer Verlag Vol. 794 of Lecture Notes in Physics, Berlin Springer Verlag, States and Transitions in Black Hole Binaries. p. 53
 Collinge et al. (2001) Collinge M. J., Brandt W. N., Kaspi S., Crenshaw D. M., Elvis M., Kraemer S. B., Reynolds C. S., Sambruna R. M., Wills B. J., 2001, ApJ, 557, 2
 Denney et al. (2010) Denney K. D., Peterson B. M., Pogge R. W., Adair A., Atlee D. W., AuYong K., Bentz M. C., Bird J. C., et al. 2010, ApJ, 721, 715
 Denney et al. (2009) Denney K. D., Watson L. C., Peterson B. M., Pogge R. W., Atlee D. W., Bentz M. C., Bird J. C., Brokofsky D. J., et al. 2009, ApJ, 702, 1353
 Done & Gierliński (2005) Done C., Gierliński M., 2005, MNRAS, 364, 208
 Edelson & Nandra (1999) Edelson R., Nandra K., 1999, ApJ, 514, 682
 Elvis et al. (1989) Elvis M., Wilkes B. J., Lockman F. J., 1989, AJ, 97, 777
 Elvis et al. (1994) Elvis M., Wilkes B. J., McDowell J. C., Green R. F., Bechtold J., Willner S. P., Oey M. S., Polomski E., Cutri R., 1994, ApJS, 95, 1
 Fabian (1979) Fabian A. C., 1979, Royal Society of London Proceedings Series A, 366, 449
 Falcke et al. (2004) Falcke H., Körding E., Markoff S., 2004, A&A, 414, 895
 Fender et al. (2004) Fender R. P., Belloni T. M., Gallo E., 2004, MNRAS, 355, 1105
 Fender et al. (2006) Fender R. P., Körding E., Belloni T., Uttley P., MHardy I., Tzioumis T., 2006 Vi microquasar workshop: Microquasars and beyond. Proceedings of Science: SISSA, Trieste (Italy), (arXiv:0706.3838)
 Freeman et al. (1999) Freeman P. E., Graziani C., Lamb D. Q., Loredo T. J., Fenimore E. E., Murakami T., Yoshida A., 1999, ApJ, 524, 753
 Gaskell (2004) Gaskell C. M., 2004, ApJL, 612, L21
 Gelman et al. (2004) Gelman A., Carlin J. B., Stern H. S., B. R. D., 2004, Bayesian Data Analysis (2nd ed). Chapman & Hall, London
 Gelman et al. (1996) Gelman A., Meng X.L., Stern H. S., 1996, Statistica Sinica, 6, 733
 Gelman & Rubin (1992) Gelman A., Rubin D. B., 1992, Statistical Science, 7, 457
 Gierliński et al. (2008) Gierliński M., Middleton M., Ward M., Done C., 2008, Nature, 455, 369
 Gierliński et al. (2008) Gierliński M., Nikołajuk M., Czerny B., 2008, MNRAS, 383, 741
 Giroletti & Panessa (2009) Giroletti M., Panessa F., 2009, ApJL, 706, L260
 Gleissner et al. (2004) Gleissner T., Wilms J., Pottschmidt K., Uttley P., Nowak M. A., Staubert R., 2004, A&A, 414, 1091
 Hayashida et al. (1998) Hayashida K., Miyamoto S., Kitamoto S., Negoro H., Inoue H., 1998, ApJ, 500, 642
 Heinz & Sunyaev (2003) Heinz S., Sunyaev R. A., 2003, MNRAS, 343, L59
 Ho & Ulvestad (2001) Ho L. C., Ulvestad J. S., 2001, ApJS, 133, 77
 Hurkett et al. (2008) Hurkett C. P., Vaughan S., Osborne J. P., O’Brien P. T., Page K. L., Beardmore A., Godet O., Burrows D. N., Capalbi M., Evans P., Gehrels N., Goad M. R., Hill J. E., Kennea J., Mineo T., Perri M., Starling R., 2008, ApJ, 679, 587
 Jones et al. (2010) Jones S., McHardy I., Moss D., Seymour N., Breedt E., Uttley P., Kording E., Tudose V., 2010, ArXiv eprints
 Kashyap et al. (2010) Kashyap V. L., van Dyk D. A., Connors A., Freeman P. E., Siemiginowska A., Xu J., Zezas A., 2010, ApJ, 719, 900
 King (2010) King A. R., 2010, MNRAS, 402, 1516
 Körding et al. (2006) Körding E. G., Jester S., Fender R., 2006, MNRAS, 372, 1366
 Kukula et al. (1995) Kukula M. J., Pedlar A., Baum S. A., O’Dea C. P., 1995, MNRAS, 276, 1262
 Lawrence et al. (1987) Lawrence A., Watson M. G., Pounds K. A., Elvis M., 1987, Nature, 325, 694
 Markowitz et al. (2003) Markowitz A., Edelson R., Vaughan S., Uttley P., George I. M., Griffiths R. E., Kaspi S., Lawrence A., MHardy I., Nandra K., Pounds K., Reeves J., Schurch N., Warwick R., 2003, ApJ, 593, 96
 McClintock & Remillard (2006) McClintock J. E., Remillard R. A., 2006, in Lewin, W. H. G. & van der Klis, M. ed., Compact stellar Xray sources Black hole binaries. pp 157–213
 MHardy et al. (2007) MHardy I. M., Arévalo P., Uttley P., Papadakis I. E., Summons D. P., Brinkmann W., Page M. J., 2007, MNRAS, 382, 985
 MHardy et al. (2006) MHardy I. M., Koerding E., Knigge C., Uttley P., Fender R. P., 2006, Nature, 444, 730
 MHardy et al. (2004) MHardy I. M., Papadakis I. E., Uttley P., Page M. J., Mason K. O., 2004, MNRAS, 348, 783
 Meng (1994) Meng X.L., 1994, Annals of Statistics, 22, 1142
 Merloni et al. (2003) Merloni A., Heinz S., di Matteo T., 2003, MNRAS, 345, 1057
 Middleton & Done (2010) Middleton M., Done C., 2010, MNRAS, 403, 9
 Middleton et al. (2009) Middleton M., Done C., Ward M., Gierliński M., Schurch N., 2009, MNRAS, 394, 250
 Miyamoto et al. (1991) Miyamoto S., Kimura K., Kitamoto S., Dotani T., Ebisawa K., 1991, ApJ, 383, 784
 Mueller & Madejski (2009) Mueller M., Madejski G., 2009, ApJ, 700, 243
 Mushotzky et al. (1993) Mushotzky R. F., Done C., Pounds K. A., 1993, Ann. Rev. A&A, 31, 717
 Nandra & Papadakis (2001) Nandra K., Papadakis I. E., 2001, ApJ, 554, 710
 Nowak (2005) Nowak M., 2005, Ap&SS, 300, 159
 Ogle et al. (2004) Ogle P. M., Mason K. O., Page M. J., Salvi N. J., Cordova F. A., MHardy I. M., Priedhorsky W. C., 2004, ApJ, 606, 151
 Papadakis & Lawrence (1993) Papadakis I. E., Lawrence A., 1993, MNRAS, 261, 612
 Papadakis & Lawrence (1995) Papadakis I. E., Lawrence A., 1995, MNRAS, 272, 161
 Park et al. (2008) Park T., van Dyk D. A., Siemiginowska A., 2008, ApJ, 688, 807
 Ponti et al. (2006) Ponti G., Miniutti G., Cappi M., Maraschi L., Fabian A. C., Iwasawa K., 2006, MNRAS, 368, 903
 Pottschmidt et al. (2003) Pottschmidt K., Wilms J., Nowak M. A., Pooley G. G., Gleissner T., Heindl W. A., Smith D. M., Remillard R., Staubert R., 2003, A&A, 407, 1039
 Pounds et al. (2004) Pounds K. A., Reeves J. N., King A. R., Page K. L., 2004, MNRAS, 350, 10
 Pounds & Vaughan (2011) Pounds K. A., Vaughan S., 2011, MNRAS, in press, arXiv:1012.0998
 Protassov et al. (2002) Protassov R., van Dyk D. A., Connors A., Kashyap V. L., Siemiginowska A., 2002, ApJ, 571, 545
 Remillard et al. (2003) Remillard R., Muno M., McClintock J. E., Orosz J., 2003, in P. Durouchoux, Y. Fuchs, & J. Rodriguez ed., New Views on Microquasars Xray QPOs in blackhole binary systems. p. 57
 Remillard & McClintock (2006) Remillard R. A., McClintock J. E., 2006, Ann. Rev. A&A, 44, 49
 Remillard et al. (1999) Remillard R. A., Morgan E. H., McClintock J. E., Bailyn C. D., Orosz J. A., 1999, ApJ, 522, 397
 Remillard et al. (2002) Remillard R. A., Muno M. P., McClintock J. E., Orosz J. A., 2002, ApJ, 580, 1030
 Revnivtsev et al. (2000) Revnivtsev M., Gilfanov M., Churazov E., 2000, A&A, 363, 1013
 Rubin (1984) Rubin D. B., 1984, Annals of Statistics, 12, 1151
 Russell (2002) Russell D. G., 2002, ApJ, 565, 681
 Shakura & Sunyaev (1976) Shakura N. I., Sunyaev R. A., 1976, MNRAS, 175, 613
 Steenbrugge et al. (2009) Steenbrugge K. C., Fenovčík M., Kaastra J. S., Costantini E., Verbunt F., 2009, A&A, 496, 107
 Strohmayer (2001) Strohmayer T. E., 2001, ApJL, 552, L49
 Syunyaev (1973) Syunyaev R. A., 1973, SvA J, 16, 941
 Terashima et al. (2009) Terashima Y., Gallo L. C., Inoue H., Markowitz A. G., Reeves J. N., Anabuki N., Fabian A. C., Griffiths R. E., Hayashida K., Itoh T., Kokubun N., Kubota A., Miniutti G., Takahashi T., Yamauchi M., Yonetoku D., 2009, PASJ, 61, 299
 Tomsick et al. (2004) Tomsick J. A., Kalemci E., Kaaret P., 2004, ApJ, 601, 439
 Ulvestad & Wilson (1984) Ulvestad J. S., Wilson A. S., 1984, ApJ, 285, 439
 Uttley et al. (2003) Uttley P., Fruscione A., MHardy I., Lamer G., 2003, ApJ, 595, 656
 Uttley & MHardy (2001) Uttley P., MHardy I. M., 2001, MNRAS, 323, L26
 Uttley et al. (2002) Uttley P., MHardy I. M., Papadakis I. E., 2002, MNRAS, 332, 231
 Uttley et al. (2005) Uttley P., MHardy I. M., Vaughan S., 2005, MNRAS, 359, 345
 Uttley et al. (2004) Uttley P., Taylor R. D., MHardy I. M., Page M. J., Mason K. O., Lamer G., Fruscione A., 2004, MNRAS, 347, 1345
 van der Klis (1989) van der Klis M., 1989, in Ögelman H., van den Heuvel E. P. J., eds, Timing Neutron Stars . p. 27
 van der Klis (1997) van der Klis M., 1997, in G. J. Babu & E. D. Feigelson ed., Statistical Challenges in Modern Astronomy II Quantifying Rapid Variability in Accreting Compact Objects. p. 321
 van der Klis (2006) van der Klis M., 2006, in Lewin, W. H. G. & van der Klis, M. ed., Compact stellar Xray sources Rapid Xray Variability. pp 39–112
 Vaughan (2010) Vaughan S., 2010, MNRAS, 402, 307
 Vaughan et al. (2003) Vaughan S., Edelson R., Warwick R. S., Uttley P., 2003, MNRAS, 345, 1271
 Vaughan & Fabian (2003) Vaughan S., Fabian A. C., 2003, MNRAS, 341, 496
 Vaughan et al. (2003) Vaughan S., Fabian A. C., Nandra K., 2003, MNRAS, 339, 1237
 Vaughan & Uttley (2005) Vaughan S., Uttley P., 2005, MNRAS, 362, 235
 Vaughan & Uttley (2006) Vaughan S., Uttley P., 2006, Advances in Space Research, 38, 1405
Appendix A Construction of ‘fluxed’ spectra
Xray spectra are conventionally given in terms of count rates per detector channel. These “raw” spectra are of limited use for data visualisation because of the distorting effect of (i) the energydependent efficiency of the telescope and detector (described by the effective area function), and (ii) blurring by the line spread function of the instrument (modelled by the redistribution matrix). The true source spectrum (in units of photon s cm keV) is related to the observed count rate in channel by a Fredholm integral equation:
(7) 
where is the instrument effective area (cm) and gives the probability that a photon of energy will be recorded in channel (i.e. the redistribution function). Here the observed spectrum is assumed to be backgroundsubtracted, and the instrumental response defined by is assumed to be linear, that is, independent of the value of which is a good approximation so long as pileup is negligible. This integral equation can be approximated using a discrete formulation in which the energy range is divided into bins :
(8) 
where is the average and and are the definite integrals of their continuous counterparts, over the energy range of each bin .
In general, methods that attempt to invert the above equations to recover from an observed spectrum are unreliable and one must turn to the ‘forward fitting’ method to build a model of the spectrum (see Arnaud, 1996). However, it is often useful to view even an approximate ‘fluxed’ spectrum, largely free from the distorting effects of the instrument efficiency, if only to give a better visual representation of the underlying source spectrum. Nowak (2005) discussed a simple method to accomplish this in which the observed count spectrum is normalised by an effective area curve that has been blurred using the redistribution function.
(9) 
This method matches the blurring of the effective area to the blurring of the observed spectrum to give an estimate of the spectrum in true flux units (e.g. photon s cm) that is, to a large extent, free from the distorting effects of the energydependent effective area. This process is essentially equivalent to calibrating the sensitivity using an observation of an ideal standard star (with a perfectly flat spectrum)^{9}^{9}9 This can be done trivially within XSPEC by defining a constant model spectrum (e.g. a power law with index and normalisation ) and plotting the “unfolded” spectrum. For the special case of a constant model this is equivalent to equation 9. The ‘unfolded’ plot shows the model in flux units () multiplied by the ratio of the observed spectrum to the ‘folded’ model, i.e. . the SAS task efluxer can also perform this transformation.. The resulting spectrum can easily be converted to conventional flux density units (e.g. erg s cm keV).
It is crucial that the effective area function is blurred by the redistribution matrix. The reason is that the effective area function may contain abrupt changes in efficiency (e.g. near photoelectric edges in the detectormirror system, such as the Kedges of O and Si, or the Medge of Au) that can be far sharper than the resolving power of the detector, which means that normalising the observed (blurred) spectrum by the exact (not blurred) effective area introduces very strong spurious features into the spectrum near the edges. Blurring the effective area curve matches the resolution to that of the data and suppresses this effect. The result is a spectrum with the energydependence of the detector/mirror efficiency removed (to a reasonable approximation), but which remains at the detector spectral resolution. We emphasise that this is not, even approximately, a deconvolved spectrum. It has the advantage over the popular ‘unfolding’ method (in XSPEC) of being independent of the spectral model, and therefore providing a more objective visualisation ‘fluxed’ spectrum.
Appendix B QPO detection and upper limits
This appendix presents first an outline of the method used to search for QPOs in the power spectral data, and second how upper limits were derived based on the detection method.
b.1 Calibrating the QPO detection procedure
An threshold detection procedure uses a test statistic to quantify evidence for a detection, and is calibrated such that a detection occurs when for some threshold value chosen to have a small probability of occurring by chance when is true. Specifically, the detection threshold is chosen such that , where is the chance of a ‘false detection’ (a ‘type I error’ when is true). Commonly used values are and . But in order to calculate we need to know the reference distribution of the test statistic on the assumption that the null hypothesis is true, . In general this may be estimated by calculating the distribution of the test statistic from a large number of Monte Carlo simulations of data generated assuming and distributed as expected for realistic data.
The above would be valid if the null hypothesis was simple, that is, did not contain any unknown (free) parameters. In the present case the null hypothesis does involve parameters estimated using the observed data, . We account for the uncertainties in their values by using random parameters drawn from the posterior distribution to generate the simulated data, using the MCMC discussed in section 3. This gives the distribution of the test statistic known as the posterior predictive distribution (see Rubin, 1984; Meng, 1994; Gelman et al., 1996, 2004; Protassov et al., 2002; Vaughan, 2010).
(10) 
In order to determine this distribution was mapped out using Monte Carlo simulations as follows:

repeat for simulated datasets:

Draw a set of parameter values from the posterior

simulate data from with parameters :

fit to get

fit to get

compute


use the distribution of to define for which
We used simulations to calibrate and found the critical value for an test was . (This was calculated using a QPO, but using gave a very similar result.) for the data and models discussed in this paper. In order to ensure the global minimum was found in each case (this can be tricky for which may produce multiple minima in ) we used the robust automated fitting scheme discussed by Hurkett et al. (2008, section 3.2.2)
b.2 Upper limit procedure
We wish to determine the upper limit on the QPO strength , which is the smallest value of such the probability of detecting the QPO (assuming it is present, i.e. ) is for some specific . See Kashyap et al. (2010) for a very clear introduction to the meaning and calculation of upper limits in astronomy.
In the present case we are searching for a QPO of unknown location (frequency) but the upper limit may vary as a function of this location – because the “background” continuum level and the uncertainty on the spectrum vary as a function of frequency – so we must calculate the limit for a set of different QPO locations . The probability of a detection, , which is (probability of ‘type II error’), is:
(11) 
and the upper limit is the smallest value of which satisfies
(12) 
for some chosen minimum power .
Given the defined above we can compute the upper limit using Monte Carlo simulations of realistic data under . These can be used to map out the distribution for different and values, and hence find that satisfies the above equations. In these terms the upper limit is the strength of the weakest QPO that would have a reasonable probability , for some specific such as or .
This method will, in general, be sensitive to the parameters of the continuum model, , which of course we do not know exactly. But, as before, we can average over the posterior distribution of these as derived from the data .
(13) 
Again, this can be done be tacking values from the posterior using the MCMC discussed in section 3.
The calculation of for a given pair of QPO parameters ( and ) was performed using a Monte Carlo scheme as follows.

loop over simulated datasets

Draw a set of parameter values from the posterior

simulate data from with continuum parameters and QPO parameters ():

fit to get

fit to get

compute


use distribution of to find
This process is repeated over a grid of QPO parameters, and , and at each frequency the smallest for which is the upper limit . For the present analysis we used logarithmically spaced values of QPO location , logarithmically spaced values of QPO strength , and at each pair generated simulations. This is a sufficient number of simulations to accurately define contours of e.g. on the  plane.
The detection procedure was calibrated using simulated data generated assuming . But to evaluate , and hence the upper limit, the data were simulated assuming . The former is concerned with the chance of spurious detection assuming no QPO present, the latter is concerned with the chance of detection when a QPO is present.