A Systematic Study of Variability in a Sample of Ultraluminous X-ray Sources

We present results from a study of short term variability in 19 archival observations by XMM-Newton of 16 Ultraluminous X-ray Sources (ULXs). Eight observations (six sources) showed intrinsic variability with power spectra in the form of either a power law or broken power law-like continuum and in some cases quasi-periodic oscillations (QPOs). The remaining observations were used to place upper limits on the strength of possible variability hidden within. Seven observations (seven sources) yielded upper limits comparable to, or higher than, the values measured from those observations with detectable variations. These represented the seven faintest sources all with f_x<3x10^-12 erg cm^-2 s^-1. In contrast there are four observations (three sources) that gave upper limits significantly lower than both the values measured from the ULX observations with detectable variations, and the values expected by comparison with luminous Galactic black hole X-ray binaries (BHBs) and Active Galactic Nuclei (AGN) in the observed frequency bandpass (10^-3 - 1 Hz). This is the case irrespective of whether one assumes characteristic frequencies appropriate for a stellar mass (10 M_sun) or an intermediate mass (1000 M_sun) black hole, and means that in some ULXs the variability is significantly suppressed compared to bright BHBs and AGN. We discuss ways to account for this unusual suppression in terms of both observational and intrinsic effects and whether these solutions are supported by our results.


INTRODUCTION
Ultraluminous X-ray sources (ULXs), point-like sources with luminosities greater than 10 39 erg s −1 , were first identified by the Einstein Observatory (Fabbiano 1989). Initial studies of short term variability indicated that they were likely to be accreting compact objects (Okada et al. 1998). Both spectral analysis and long and short term variability studies suggested similarities between ULXs and Galactic black hole X-ray binaries (BHBs) (comprising a stellar mass black hole accreting matter from a partner star (Miller & Colbert 2004;Fabbiano 2006;Roberts 2007)). The Xray spectra of ULXs can be fitted by models commonly used for BHBs (Colbert & Mushotzky 1999), and display changes in their energy spectra similar to the state changes observed in BHBs (La Parola et al. 2001). The high luminosities are greater than that predicted by applying the Eddington limit allowed by spherical accretion onto a typical stellar mass black hole ( ∼ < 20 M ⊙ ). ULXs are, however, mostly observed outside the nucleus of their host galaxies and thus are unlikely to be super massive black holes ( ∼ > 10 5 M ⊙ ), typically observed in the centre of galaxies as active galactic nuclei (AGN). Their high observed luminosities and the low appar-⋆ lmh38@star.le.ac.uk ent disc temperatures measured from the energy spectra have led to suggestions that these sources represent an until-recently unknown class of black holes with masses between those of stellarmass black holes and AGN around 100 − 10000 M ⊙ : intermediate mass black holes (IMBHs) (Colbert & Mushotzky 1999). In contrast a number of theories have been put forward suggesting that the high luminosities may be caused by a stellar mass black hole emitting or accreting in different ways. Emission processes such as anisotropic emission (King et al. 2001), relativistic beaming from jets (Georganopoulos et al. 2002), or accretion onto the BH over the expected limit (Begelman 2002) could explain some observed features. Constraining the masses of these objects would be an important step forward in understanding them.
A number of different relationships between BH mass, accretion rates and features of power spectral densities (PSDs) common to both BHBs and AGN have already been identified (McHardy et al. 2006;McClintock & Remillard 2006;van der Klis 2006;Done & Gierliński 2005;Markowitz et al. 2003;Uttley et al. 2002;Belloni & Hasinger 1990). If these key features are seen in the PSDs of ULXs then it may provide a way to estimate the black hole mass. BHBs are observed in a number of different "states" traditionally defined by the dominance of components within the energy spectra, which also show distinctly different PSDs. It is now recognised however that the accretion rate may in fact be a more effective parameter to discern them (McClintock & Remillard 2006). Hardness-Intensity diagrams following BHBs through outbursts have revealed a form of hysteresis path that these sources follow as they move between different states (Miyamoto et al. 1995).
A typical outburst cycle begins in a 'low/hard' state in which the X-ray spectrum is dominated by a hard power law, with a weak contribution from a thermal component originating in the accretion disc. The X-ray PSD shows band-limited noise that may be approximated by a doubly broken power law or multiple broad Lorentzians (Nowak 2000;McClintock & Remillard 2006). In this state the PSD often displays quasi-periodic oscillations (QPOs). Radio emission from a steady jet is frequently observed. At the peak of outburst the source may reach a 'very high state' in which the X-ray spectrum shows a steep power law and a strong thermal component. The PSD typically features a complex, broad continuum and strong QPOs. The third main state is the 'high/soft' state characterised by a dominant soft, thermal spectrum. In this state the PSD appears as a broken or bending power law with an index of −2 above the break and −1 below it, which may extend to very low frequencies (Reig et al. 2003). QPOs are not usually seen and the radio emission is greatly reduced. Towards the end of outburst the source returns to the 'low/hard' state and the luminosity decreases until the source reaches quiescence.
Observations of ULXs, particularly with XMM-Newton, have allowed for X-ray spectral and timing analysis which appears to show some of the features observed in BHBs. The PSDs from XMM-Newton are most sensitive to frequencies around which we would expect to see a high frequency break (present in both the high/soft and low/hard states) if the mass of the black hole is sufficiently large (i.e. the IMBH model predicts a break at ∼ 10 −4 − 10 −2 Hz).
Observations in all three states of BHBs, as well as all wellconstrained AGN data, show a high-frequency break above which the PSD decays as a steep power law (index ∼ < −2). This break frequency may be a key diagnostic feature in establishing the properties of the system. McHardy et al. (2006) showed that the break frequencies in AGN and BHBs display a correlation dependent on the mass of the black hole and its accretion rate T Break ∝ M 1.12 Edd . BHBs in the low/hard state have a second lower frequency break from a flat power law to an index of −1. Körding et al. (2007) have altered and extended the analysis from McHardy et al. (2006) to BHBs in the low/hard state through modelling the PSD as a series of Lorentzians (apparent when the PSD y-axis is plotted as P × f ) and then forming a relationship between accretion, BH mass and the Lorentzian relating to the high frequency break. QPOs have been detected in two ULXs, M82 X-1 and NGC 5408 X-1, indicating that the systems may have similar variability properties to BHBs (Strohmayer & Mushotzky 2003;Strohmayer et al. 2007). By contrast, other studies of ULX variability have indicated that these sources often do not show strong short term variability (Feng & Kaaret 2005;Roberts 2007). However a break in the PSDs of ULXs NGC 5408 X-1 and M82 X-1, from a flat power law to -1 dependence, similar to that seen in the low/hard state of BHBs, has been identified below the QPOs by Soria et al. (2004) and ; Mucciarelli et al. (2006) respectively. Identifying the higher characteristic break frequency in the PSDs of ULXs could allow for a more accurate determination of their masses, and a deeper understanding of the nature of these particular sources.
This paper describes a power spectrum analysis of the time series obtained from 19 XMM-Newton observations of ULXs. Sev-eral of these observations have been analysed and discussed separately elsewhere (see section 2), but this paper presents the first uniform analysis and presentation of the power spectra from a reasonable-sized sample of good quality ULX observations (i.e. observations with good exposure time of sources with high X-ray fluxes). As noted above, some ULX observations show no evidence for rapid X-ray variability, but in almost all published studies (e.g. Feng & Kaaret 2005) it is not clear whether the lack of detected variations is due to insufficient data or an intrinsically low variability amplitude (by comparison with well studied BHBs and AGN, and also the variable ULXs). Indeed, the careful analysis of Ho II X-1 by  showed the ULX to be intrinsically under-variable. This raises some interesting questions, but it is not clear whether the lack of strong variability is common in ULXs: in all other previous studies where no variability was detected, there was no attempt to constrain the amplitude of intrinsic variation allowed. In the present study we analysed the power spectra of all 19 observations, and where no variability was detected, we followed the approach used by  to constrain the power spectrum. This allowed us to compare the variability properties of all the ULXs, whether variability could be detected or not, in a way not possible using the results of previously published studies.
The rest of this paper is organised as follows. Section 2 describes our sample selection criteria and search process. In Section 3 we discuss the basic data reduction techniques used and methods of temporal analysis. Section 4 explains our results in some detail including investigations into the upper limits of variation hidden within the power spectra. Finally Section 5 gives a brief discussion and conclusion.

SAMPLE SELECTION
We have performed a search for bright ULXs with a flux previously observed to be greater than 0.5×10 −12 erg s −1 . A sample was formed through use of sources previously identified in Stobbart et al. (2006), searching the ROSAT catalogues of Liu & Bregman (2005) and Liu & Mirabel (2005) and by carrying out an archival study of papers studying ULXs using the ADS database (ADS 1 ). A search was then made for long observations (> 25 ks) of regions including these sources available by February 2008 in the XMM − Newton public data archive (XSA 2 ) . The final sample consists of 19 observations from 16 sources with the potential to provide useful timing data (see Table 1).
As mentioned above, several of these observations have been discussed elsewhere in the literature, and here we briefly summarise these reports. The observations of NGC 1313 X-1 and X-2, NGC 5204 X-1, NGC 4945 X-2, NGC 4861 ULX and NGC 253 PSX-2 have no published discussion of their short-term variability.
The observations of NGC 55 ULX were known to display interesting variability (Stobbart et al. 2004), but this has not previously been analysed in detail. No variability was detected by Feng & Kaaret (2005) in the observations of NGC 4395 X-1, NGC 3628 X-1, M83 ULX and NGC 2403 X-1. Barnard et al. (2007) detected no variability in NGC 4559 X-1 (revising the previous report by Cropper et al. 2004). Only the remaining five observations have detailed power spectral analyses published to date. The PSDs from the two observations of M82 X-1 have been published in Strohmayer & Mushotzky (2003) and  respectively: these show a strong QPO and a break in the PSD. A similar QPO and spectral break has also been detected in the observation of NGC 5408 X-1 (Strohmayer et al. 2007).  claimed a modestly significant QPO in the power spectrum of the observation Ho IX X-1. Finally, as mentioned above,  presented a detailed analysis of the lack of variability present in the long observation of Ho II X-1.

ANALYSIS
Light-curves and energy spectra were extracted from the EPIC pn camera on XMM-Newton. In all observations the camera was operating in full-field mode excepting that of Ho IX X-1 (obs Id. 0200980101) which was taken in large window mode. The MOS cameras were not utilised due to their lower time resolution, making it harder to constrain the Poisson noise level during power spectral analysis. When the data were included they were found to add little to either the results or analysis. The data from each observation was extracted and processed using the XMM-Newton SAS version 7.1.0. Data was filtered leaving only events with PATTERN 4 and FLAG==0. Events were extracted either within the optimised radius of the source, as evaluated by the SAS, or a radius of our own definition if other factors (such as the source sitting close to a chip edge or near another X-ray point source) made this unsuitable. A background light curve was taken from a rectangular region on the same chip and as close to the source as possible.
Light-curves for the 0.3-10 keV band were extracted with a time resolution of 73.4 ms for all sources except for Ho IX X-1, the light curve of which had a time resolution of 96 ms accounting for the different camera operational mode. Most of the ULXs are near the aimpoint of the image and are clearly detected as separate point sources, with a few exceptions. Most sources were therefore extracted from a region centred on the co-ordinates given in Table 1. NGC 4395 X-1 is close to the edge of the chip; M82 X-1 and NGC 3628 X-1 are close to the centre of their host galaxies and other bright X-ray sources, causing difficulties in exactly identifying the source whilst excluding all others given the resolution of XMM-Newton. For M82 X-1 we followed Strohmayer & Mushotzky (2003) and extracted an 18" region centred on the co-ordinates for the point source identified clearly in Chandra observations. Analysis on power spectra taken from regions surrounding the bright ULX in the XMM-Newton observations by  has suggested that this is where the QPO originates. NGC 3628 X-1 was extracted in a 23" region centred on the co-ordinates given in Table 1. The second observation of NGC 55 ULX (obs. Id. 0028740101) is off-axis, causing a lower average count rate than observed in the first observation (obs. Id. 0028740201) and so the two observations were treated separately.
In order to remove background flares a light curve was extracted from the entire chip in the 10−15 keV band and then criteria applied to exclude regions where the count rate stayed continuously over three times the mean. We only considered flares longer than 30s. Any telemetry effects where the count rate dropped to zero for a period longer than 15 s were replaced with the mean count rate in the source file (this had only a small effect on the overall light curve and subsequent temporal analysis). The second half of the light curve from Ho IX X-1 (obs Id. 0200980101) was particularly affected by this. The continuous regions not affected by flaring were defined as good time intervals.
Luminosities were measured through fitting each of the en-ergy spectra with a multicoloured disc + power law model and are included in Table 1.

Timing Analysis
In order to estimate the PSD from an observation the available good time was split into at least five segments of equal length. The FFTs from each segment were averaged and the result rebinned in frequency giving a factor 1.1 increase in frequency per bin (i.e binning over f i → 1.1 f i ), but ensuring there were at least 20 points in each. This averaging process permits minimum-χ 2 fitting to be used as a maximum likelihood method (van der Klis 1989a; Papadakis & Lawrence 1993). The PSD was normalised to the fractional rms squared (van der Klis 1997). Visual examination revealed no QPOs in addition to those already reported in M82 X-1 (Strohmayer & Mushotzky 2003) and NGC 5408 X-1 (Strohmayer et al. 2007). Figure 1 shows a selection of eight of the sample of PSDs, these are all plotted with the same segment size for ease of comparison, it must be noted that some of the PSDs used for the analysis were produced with larger segments and this allows a few more points of data to be visible at low frequencies.
At high frequencies the PSD is expected to be dominated by Poisson noise, random fluctuations in the photon counts even in the absence of intrinsic variability in the source flux. This produces a constant (i.e. flat) additional component in the PSD (van der Klis 1989a) with an expected noise level (in the absence of detector nonlinearities) of approximately 2(S + B)/S 2 where S is the source count rate and B is the background count rate, calculated from the counts in the good time intervals (Uttley et al. 2002;Vaughan et al. 2003). In order to determine the Poisson noise level in each observation the high frequency region of the PSD ( 0.5 Hz) was fitted with a constant and the normalisation compared to the expected Poisson noise level. Generally these correspond within the ∼ 2% errors on the PSD estimate. However for the two brightest sources, Ho II X-1 and M82 X-1, the measured values were ∼ 3% below the prediction levels. These differences may be due to instrumental effects such as event pile-up or dead time (see also ).

RESULTS
The following subsections describe the results of fitting various simple models to the power spectral data.

Variable Sources
In the absence of any intrinsic variability from the source the power spectrum is expected to be constant with frequency resulting from Poisson noise in the photon counting detector. A reasonable test for the presence of intrinsic variability is therefore a χ 2 goodness-offit test comparing the data to a constant. The results are given in Table 2. Sources with a p−value 0.05 were considered to show significant intrinsic variability (with this threshold one might expect ∼ 1 false detection in a sample of 19 observations). In fact six sources (eight observations) showed significant variability: Ho IX X-1, NGC 5408 X-1, NGC 1313 X-1, NGC 1313 X-2, NGC 55 ULX and M82 X-1. Investigation of the PSD from NGC 4559 X-1 revealed it to be flat with no indication of the variability features analysed in Cropper et al. (2004), this is believed to be due to artificial variability caused by the analysis method as discussed in Barnard et al. (2007).   (2)  Source Name With the exception of the second observation of M82 X-1 and the observation of NGC 1313 X-1, which display variability in the power spectrum close to the Poisson noise level Mucciarelli et al. 2006), the variable objects all show a clear rise towards low frequencies compared to the observations with flat PSDs (e.g. Ho II X-1, NGC 5204 X-1). The continuing rise seen in XTE observations of M82 X-1 with a power law of -1 (Kaaret et al. 2006;) supported by evidence of breaks seen in NGC 5408 X-1 suggests that the PSDs are best fit with a broken power law model over a broad Lorentzian. For this reason power law models were fitted to the variable sources, although the limited number of points at low frequencies made it difficult to constrain the model parameters. The results are given in Table 3.
Ho IX X-1, NGC 1313 X-2, NGC 55 ULX and M82 X-1 (obs. Id. 0112290201) were fitted best with a power law (plus a constant for the Poisson noise). NGC 5408 X-1, M82 X-1 (obs Id 0206080101) and NGC 1313 X-1 were fitted better with a broken power law (plus a constant) as used in Strohmayer et al. (2007) to model the PSD from NGC 5408 X-1. This is the first report of a broken power law PSD in NGC 1313 X-1 and in order to test the plausibility of the detection we have attempted to fit the PSD with a power law (+ constant) model, and compared the χ 2 value from this fit to that of a broken-power law (+ constant) model which has two further free parameters. Analysing the difference between the two χ 2 values both through direct comparison and the F-test suggest that the spectral break at 0.09 Hz is significant and necessary to the fit with p null ≈ 10 −4 . Fig. 2 shows the standardised residuals after fitting these two models.
In the cases of M82 X-1 and NGC 5408 X-1 the obvious QPOs were modelled with narrow Lorentzians (see below). The power law-like regions can all be fitted with slopes between −2 and −0.5. For those sources with previous analyses of their variability (M82 X-1, NGC 5408 X-1 and Ho IX X-1) the power law slopes, break frequencies and Lorentzian parameters found agree with the previous results, with the exception of the QPO identified in Holmberg IX X-1 by  which was not reproduced with high significance although a continuum PSD is visible.

Establishing limits on the variability
The other 11 observations show no strong evidence for intrinsic variability, despite a range of observation lengths and source fluxes, and therefore sensitivity to variability. In these cases the variability amplitude was constrained, following the procedure of , by assuming the intrinsic power spectrum has one of two generic forms common to the PSDs of BHBs and AGN.
The first model is a broken power law model. It is well established that BHBs and AGN show PSDs that are usually dominated by a broad-band noise component with a form that resembles to first-order a broken power law. Most sources show a frequency range in which the spectral index (PSD slope) is approximately −1 -noise with this PSD are known as "flicker noise" (Press 1978). The integral of a f −1 spectrum diverges (if rather slowly) to low and high frequencies, and so must break to flatter and steeper slopes respectively. Therefore, we may use a broken power law as an approximate model of the generic broad-band noise PSD.
In many AGN and "high/soft" state BHBs the low frequency break (to an index flatter than −1) is not observed and must occur at very low frequencies (Reig et al. 2003;Markowitz et al. 2003;Uttley et al. 2002;Reig et al. 2002;Churazov et al. 2001). This means a good approximation to the observed PSDs can be made by considering a singly broken power law (BPL) with indicies of −1 and −2 below and above some break frequency, respectively.
where A is a normalisation, f b is the (high frequency) break frequency, and C P is a constant to account for the Poisson noise.
If we represent this spectrum in f × P( f ) form the f −1 part of the intrinsic (i.e. after subtracting the contribution due to Poisson noise) spectrum ( f < f b ) becomes flat, with a constant level It has been suggested that the value of C 1/ f is fairly constant across BHBs and AGN, with values in the range 0.005 − 0.03 Papadakis 2004). Of course the break frequencies f b vary widely from source to source -from 10 −7 Hz for AGN with BH masses ∼ 10 8 M ⊙ , to ∼ 10 Hz in stellar mass BHBs -and on long timescales the break frequency may vary by as much as a decade in frequency for a single source (Belloni & Hasinger 1990;Uttley et al. 2002;Belloni et al. 2002).
For the non-varying ULX observations there was generally no observational constraint on f b and so C 1/ f was estimated (or constrained) as a function of break frequency over the range 10 −4 − 1 Hz. Figure 3 (left panel) shows a contour plot of the upper limit on C 1/ f as a function of f b for NGC 5204 X-1 (obs. Id. 0405690201). The contours show the value of C 1/ f at which the χ 2 fit statistic increased by 1.0, 2.71, and 6.63 over the minimum. These correspond to the 68.3, 90 and 99 per cent quantiles of the χ 2 distribution with ν = 1 degrees of freedom, and are conventionally used to delimit confidence regions on one parameter. Here they are used to define approximate limits on the magnitude of C 1/ f as a function of f b . Table 4 gives 90% limits on the value of C 1/ f for each of the sources in the sample calculated assuming two different values of f b , namely 10 −3 Hz and 1 Hz. For the variable sources for which an unbroken power law was a good fit we estimate the value of C 1/ f by fitting a BPL model. Where the observed power law index was ≈ −2 (i.e. in NGC 55 ULX and NGC 1313 X-2) it was assumed that the high frequency break, f b , must lie below the observed frequency range ( ∼ < f ew × 10 −3 Hz), and so we assumed f b = 10 −4 Hz. If, on the other hand, the observed power law index was close to −1 (i.e. Holmberg IX X-1, and the first observation of M82 X-1), it was assumed that the f b must lie close to, or higher than, the highest frequency at which ULX variability can be detected above the Poisson noise. In these cases we assumed f b = 1 Hz.
In the cases of NGC 1313 X-1, NGC 5408 X-1 and the second observation of M82 X-1, where a break in the power law is detected, the observed indices at low frequencies are rather low (closer to 0 than −1) indicating that the observed PSD covers the low-intermediate part of the model. Here it is flat, breaking from an index of 0 to −1 the low frequency flattening discussed earlier.
In these cases the intermediate-high part of the PSD was again assumed to break from an index of −1 to −2 at f b = 1 Hz, while the low frequency index of 0 was fixed, the location of the lower frequency break was left as a free parameter.
The second model represents the band-limited noise (BLN) commonly observed in BHB low/hard states (Nowak 2000;Belloni et al. 2002). In this model the variability is mostly limited to ∼ 1 − 2 decades of temporal frequency (van der Klis 2006). A simple model for this type of noise is a broad (incoherent) Lorentzian: where R 2 is a normalisation term (approximately integrated power, Nowak 2000), f 0 is the resonant frequency, and Q is the quality factor. Again, C P is a constant to account for the Poisson noise. The Q-factor defines the frequency width of the Lorentzian, Q ≈ f 0 /HWHM where f 0 is the central Lorentzian frequency and HWHM its half width at half maximum. Typically a feature with Q > 2 is classified as a QPO and Q < 2 as band-limited noise (van der Klis 1989aKlis , 2006. In order to model band-limited noise Source     a value of Q = 0.5 was assumed. The same model form but with Q = 5.0 was used to model QPO features. Of course, given the quality of the data the resonant frequency f 0 of the BLN could not be constrained, and so the strength of any BLN was estimated (or constrained) as a function of f 0 spanning 10 −4 − 1 Hz, in the same manner as with the BPL model. In those observations that showed no obvious QPOs the same procedure was used to put limits on the strength of any QPOs as a function of QPO frequency. The centre and right panels of Fig. 3 show the upper limits on amplitudes of the BLN and QPO models for NGC 5204 X-1 (obs Id. 040569201). These were found through minimising the fit to the models over 1000 different values of either the central Lorentzian frequency or the break frequency over the XMM-Newton frequency range. Table 4 gives the approximate limits on the amplitudes of the two models (BLN and BPL) assuming two representative values of the characteristic frequency in each case (1 mHz and 1 Hz). The limits were calculated using a ∆χ 2 = 2.71 criterion (which corresponds to the 90 per cent limit of a χ 2 distribution with ν = 1). The fainter sources (such as NGC 4945 X-2, NGC 4395 X-1, NGC 4861 ULX, NGC 3628 X-1 and M83 ULX) have upper limits such that variability from these sources would have to be particularly strong to have been detectable. For example, in NGC 4395 X-1 a BLN spectrum could have R ≈ 20% fractional rms and this would barely reach the 90% confidence limit even assuming a characteristic frequency best-suited to detection ( f 0 ∼ 10 −3 Hz). For those sources in which variability was detected the measured values of the amplitude are given for each model. The observation of NGC 5408 X-1 gives the amplitude of a BLN spectrum (assuming f 0 ≈ 10 −3 Hz) as R ≈ 18% whereas the upper limit on the amplitude of the same model in the case of NGC 5204 X-1 (obs Id. 0405690201) gives R ∼ < 7% even though the two observations show comparable Poisson noise levels.

Limits on QPOs
The data were examined for possible QPOs by re-binning the FFT data in both frequency and time using various different bin widths. As discussed above a narrow (Q = 5) Lorentzian was used to model a QPO. For each observation the change in χ 2 was calculated as a function of the two free parameters R 2 and f 0 . An improvement of ∆χ 2 = 11.3 was considered indicative of a QPO. This process clearly recovered the known QPOs in NGC 5408 X-1 and M82 X-1 but did not find any further candidates. The known QPOs show typically R 2 ∼ 10% and f 0 ∼ 2 − 5 × 10 −2 Hz. By comparison the limits on the QPO strength derived from observations of several other bright ULXs, such as NGC 5204 X-1 and Ho II X-1, show that such QPOs would have been detected if they were present (see Fig. 3). The QPO claimed in Ho IX X-1 by  was not reproduced with high significance. In these sources QPOs similar to those in NGC 5408 X-1 and M82 X-1 sources would have to be fairly muted if present.

Source
Broken Power Law (BPL) Band limited noise (BLN) Name C 1/ f (10 −4 ) C 1/ f (10 −4 ) R 2 (10 −2 ) R 2 (10 −2 ) R 2 m (10 −2 ) 10 −3 Hz 1Hz 10 −3 Hz 1Hz (1) (2)  Table 4. Upper limits on the fitted model amplitudes using two different power spectral models. The estimate of C 1/ f from the BPL model at f b = 10 −3 Hz for Ho II X-1 differs from the value found in  differs due to the use of different binning for the fitting process. Columns: (1) Source Name; (2) C 1/ f from the broken power law model with f b = 10 −3 Hz; (3) For non-varying sources: C 1/ f from the broken power model with the f b = 1Hz; For varying sources: C 1/ f from fit to model according to shape of PSD: See text; (4) Upper limit on fractional rms within a Lorentzian of coherence Q=0.5 fit at 10 −3 Hz; (5) Upper limit on fractional rms within a Lorentzian of coherence Q=0.5 fit at 1 Hz; (6) Fit normalisation values to Band limited noise when present. All limits correspond to ∆χ 2 = 2.71 (which is approximately 90% confidence).
• Of these six, two show strong QPO features (NGC 5408 X-1 and M82 X-1, as previously identified) and all show a continuum spectrum consistent with a power law or broken power law model.
• Of the remaining 10 sources, the intrinsic variability amplitudes were constrained by comparing the data to simple PSD models based on BHB observations.
• The limits obtained for several of the observations (e.g. NGC 4559 X-1, NGC 5204 X-1, Ho II X-1) are such that the strength of the intrinsic variability, in the frequency range observed, must be substantially lower than the observed power seen in the variable sources.

Break Frequencies and BH mass
It is in principle possible to draw inferences about the nature of the putative accreting black hole in the ULXs if analogies can be formed between their observed properties and those of BHBs and AGN. Perhaps the simplest and best studied relation linking BHBs and AGN is the correlation between the PSD break frequencies and BH mass (e.g. Markowitz  The 3 PSDs of NGC 5408 X-1, M82 X-1 (obs. ID. 0206080101) and NGC 1313 X-1 show a break from an index of ∼ −1 above the break to ∼ 0 below, which resembles the 'low frequency break' observed in low/hard BHB observations (Belloni & Hasinger 1990;Nowak 2000;McClintock & Remillard 2006). Soria et al. (2004) have used the break frequency in NGC 5408 X-1 to derive an expected BH mass of ∼ 100 M ⊙ . For M82 X-1  measured the PSD break frequency in relation to the QPO frequency and drew a comparison with BHB features in order to estimate the BH mass at 25−520 M ⊙ . Casella et al. (2008) applied the relationship in Körding et al. (2007), which extends the analysis from McHardy et al. (2006) into systems in a low/hard state, to estimate a BH mass between 115 − 1305 M ⊙ for NGC 5408 X-1 and 95 − 1260 M ⊙ for M82 X-1.
By contrast the break in NGC 1313 X-1 at 0.09 Hz is comparable to those observed in BHBs, for example the observed break of Cyg X-1 in the low/hard state is 0.08 Hz. Cyg X-1 is thought to have a BH mass of ∼ 10M ⊙ . However the observed luminosity of Cyg X-1 in this state is only ∼ 4×10 37 erg/s (Nowak et al. 1999) a factor of 100 lower. Either NGC 1313 X-1 hosts a commensurately larger BH mass than Cyg X-1, or it does not but is accreting at a much higher relative rate (L/L Edd > 1), or at least appears to be due to anisotropic emission (or some combination of these). In either case the flat PSD below 0.09 Hz is difficult to account for. The flat slope is most commonly observed in 'low/hard' state BHBs, but this state is thought to occur when L/L Edd ∼ < 0.05 (McClintock & Remillard 2006). Assuming NGC 1313 X-1 to be in such a state would therefore require M BH ∼ > 600M ⊙ , yet the frequency of the break matches well the expectation for a M BH ∼ 10M ⊙ black hole. This supports the suggestion that ULXs may represent a different, 'ultraluminous' state (Roberts 2007).
NGC 1313 X-2, NGC 55 ULX and Ho IX X-1 show excess power rising to low frequencies below ∼ 5 × 10 −3 Hz. As can be seen in Figure 1, the PSD from the first observation of NGC 55 ULX (obs. ID. 0028740201) shows a steep red noise region at low frequencies, also visible in the PSD from the second observation (obs. ID. 0028740101). When the two observations are fitted with a power law + constant model, spectral indices of 1.75±0.2 and 1.96±0.4 are found respectively, with no sign of a low frequency break. This would suggest, assuming the system is in the high/soft state, that the high frequency break would have to be below 10 −4 Hz, indicative of a BH with a mass in excess of 10 6 M ⊙ (Roberts 2007) more typical of AGN. However it must be pointed out that dips in the light-curve have been identified (Stobbart et al. 2004) and that these features affect the PSD on timescales commensurate with the strong variability power seen in its PSD. The power laws seen in the PSDs from Ho IX X-1 and NGC 1313 X-2 have spectral indices of ∼0.5 and ∼1.91 respectively. The few rising points visible in the PSD of NGC 1313 X-2 do not allow for accurate fitting and so little can be derived from this value. For Ho IX X-1 however, the slope is much clearer and we have carried out a basic analysis. We have been able to constrain the likelihood of a brake to a slope of -2 to be above 0.1 Hz. This would indicate a BH with a mass below 1000 M ⊙ (Roberts 2007).
The possibility of estimating the masses of ULXs by comparing characteristic (e.g. break or QPO) frequencies with those observed in BHBs and AGN, is dependent on there being a suitably strong analogy mapping PSD characteristics between the source, together with a simple mass-frequency scaling relation (McHardy et al. 2006). The recent discovery of an apparent high frequency QPO within a AGN (Gierliński et al. 2008) indicates that these features are common across the range of observed sources. If however ULXs do not display variability in forms typical to either BHBs or AGN and instead exist in a new 'Ultraluminous state' (Roberts 2007) with distinct and different characteristics to those states observed in AGN and BHBs, then these mass estimates will no longer be reliable.

Variability vs Flux and Luminosity
The two main properties of an observation that determine the sensitivity to intrinsic spectral power at a given frequency are the length of the observation T (the number of 'cycles' of that frequency that were observed) and the count rate of the source (which determines the Poisson noise in the observation). Denoting the mean count rate as I we may approximate the significance (in units of σ) of an excess in power over the Poisson noise level by n σ van der Klis 1989b) where P S ( f ) is the power density due to the intrinsic variations in the source flux, and ∆ f ( f ) gives the frequency binning. Thus, with similar observation lengths T and frequency re-binning ∆ f ( f ) the sensitivity should scale as n σ ∝ IP S ( f ) and so if a detection is made when n σ exceeds some number the limit on the power that can be detected scales as P S ( f ) ∝ I −1 . Fig. 4 (left panel) shows the estimated variability amplitudes (detections and upper limits) against source flux. Also on these plots are typical values for C 1/ f observed in a BHB Cyg X-1 and the AGN Mrk 335, illustrating the consistency in C 1/ f between BHBs and AGN. The upper limits clearly improve (i.e. give tighter limits) at higher fluxes, as expected. In particular, they indicate it is in general not possible to detect even strong variability (C 1/ f > 0.05) in sources fainter than ∼ 10 −12 erg s −1 cm −1 . Fig. 4 (right panel) shows the same variability amplitudes against observed X-ray luminosity, with no strong trend apparent. Above 10 39 erg s −1 (i.e. definite ULX luminosities) there is a range of ∼ 30 in the variability amplitude, confounding any efforts to identify source states on the basis of luminosity-variability correlations. Considering only those source for which variability was detected (filled circles) there appears to be a correlation suggesting a decrease in variability with increasing luminosity. However there are two reasons why this could be misleading. The first is that the two highest estimated values of C 1/ f are the observations of NGC 55 ULX. This source shows dipping episodes in its variability (Stobbart et al. 2004), not seen in other ULXs, and this anomalous behaviour enhances the variability amplitude. Removing these points substantially decreases the apparent significance of any luminosity-amplitude anti-correlation. Secondly, including the limits from observations that did not have variability detections shows that some lower luminosity sources (∼ 10 39 erg s −1 ) must have very low amplitude variability, opposite to the prediction of a general luminosity-amplitude anti-correlation. When these points are considered there appears to be no obvious correlations within the plot.

Are some ULXs significantly less variable?
The upper limit contours from four bright sources (Ho II X-1, NGC 2403 X-1, NGC 5204 X-1 and NGC 4559 X-1) show that these sources are significantly less variable than others (such as NGC 55 ULX, NGC 5408 X-1, etc.). Fig. 5 compares the data from NGC 5204 X-1 with the model obtained by fitting the data from NGC 5408 X-1 (these two have similar luminosities and count rates, and hence Poisson noise levels in their PSDs). The figure clearly shows the variability seen in NGC 5408 X-1 to be absent or greatly suppressed in NGC 5204 X-1.
The difference between the strength of variability observed in many ULXs and that expected by analogy with BHBs was confirmed using simulated data. Two PSDs of Cyg X-1 -one showing BLN ("low/hard state") and one showing BPL ("high/soft state") 3 -were used as examples of BHB variability. Random time series were generated from these PSDs (using the method of Timmer & König 1995) and given the same mean count rate and exposure time as the XMM-Newton observation of Holmberg II X-1. Poisson noise was included at the appropriate level, and the simulated data were analysed in the same manner as the real ULX data. Variability was clearly detected (i.e. a constant was rejected by a χ 2 goodness-of-fit test) in the simulated data for both PSD shapes. The same was true when the PSD frequencies were reduced (by factors of up to 10 5 ) prior to generating the time series -i.e. simulating the effects of larger black hole masses by scaling the variability timescales.
The results confirm that if Holmberg II X-1 was variable with a PSD shape and amplitude similar to that seen in Cyg X-1 its variability should have been detected (at high significance) in the XMM-Newton observations. This is true whether the PSD shape is BPL or BLN, and holds when the characteristic variability frequencies are reduced as might be expected for IMBHs. The results of the simulations confirm the the PSD fitting results (discussed above), namely that the apparent lack of variability in the observations of  sources such as Holmberg II X-1 and NGC 5204 X-1 must be intrinsic to the sources, not a result of insufficient observations.

Interpreting the lack of variability
There are two fundamentally different ways to account for the lack of variability in the ULX observations: as an observational effect caused by the limited energy or frequency bandpass of the observations, or as a result of the sources being intrinsically under-variable. We consider some specific possibilities below. The energy spectra of BHBs in their higher accretion rate states are often dominated by a (quasi-)thermal component associated with the accretion disc. By comparison with the nonthermal (power law-like) hard X-ray continuum, this spectral component shows little short term variability (Churazov et al. 2001). LMC X-3 is an example of a BHB that often shows very low amplitude short term variability in the "high/soft" state (dominated by thermal disc emission; Nowak et al. 2001). See also a discus-sion of this point in the context of BHB-AGN comparisons in Done & Gierliński (2005). It could be that a similar component dominates the spectra of the under-variable ULXs within the XMM-Newton bandpass. If this was the case we would expect the thermal emission to peak around ∼ 1 − 2 keV, higher than expected for IMBH models, which predict a much cooler temperature ∼ 0.2 keV (Miller et al. 2003Stobbart et al. 2006), i.e. this solution requires BH masses similar to Galactic BHBs. One problem with this hypothesis is that the energy spectra of ULXs often appear more complex than the thermal-dominated spectra of "high/soft" BHBs when viewed in the highest quality XMM-Newton observations (Stobbart et al. (2006); Gladstone et al. (2009)), although even then there are some ULXs that do appear to have modified disc-like spectra (e.g. Vierdayanti et al. (2006); Hui & Krolik (2008)). A prediction of this hypothesis is that the sub-variable ULXs should show stronger variability at higher energies, where the spectrum is dominated by variable non-thermal emission.  suggested that the lack of variability seen in Holmberg II X-1 may be related to the so-called 'χ-state' observed in GRS 1915+105, where the variability power is concentrated in a BLN spectrum at high frequencies (> 0.1 Hz) and there is often little variability on longer timescales (Belloni et al. 2000). The integrated power over the 10 −4 − 0.1 Hz range is typically only ∼ 10 −3 (fractional variance; see e.g. Zdziarski et al. 2005) which is close to the detection limits of the ULX observations. Again, this explanation requires BH masses comparable to that of GRS 1915+105 (∼ 14 M ⊙ ; Greiner et al. 2001) as a significantly larger BH mass would shift the BLN variability into the frequency bandpass of the ULX observations (10 −4 − 0.1 Hz) where it would be detectable. A prediction of this hypothesis is that ULXs should show significant variability power at higher frequencies (e.g. 0.2 − 20 Hz), in the form of a BLN component and possibly strong QPOs, as does GRS 1915+105.
The above explanations hypothesise that the variability in ULXs is similar to that of BHBs (and AGN) but is often not observed due to the limited energy and/or frequency bandpass of the available data. Another possibility is that the X-ray emission from ULXs is intrinsically less variable (on frequencies > 10 −3 Hz) than expected by analogy with other accreting systems.
The timescales involved are sufficiently large to rule out at least one model for the suppression of variability within the ULX system: scattering in an extended region. Recent studies of the X-ray spectra of ULXs (e.g. Done & Kubota 2006;Stobbart et al. 2006;) have put forward the idea of an optically thick corona of hot electrons that modifies the emerging energy spectrum. The timescale for the diffusion of photons through an electron scattering medium is of order ∆t ∼ Rτ es /c, where R is the size of the scattering region and τ es ( 1) is the optical depth to electron scattering (see e.g. Nowak et al. 1999;Miyamoto et al. 1988, and references therein). On timescales shorter than ∆t the variability will be strongly attenuated. The ULX observations show weak or absent variability on timescales at least as large as ∼ 10 2 s, which would require Rτ es ∼ 100 lt-s. This is implausibly large; even with M BH ∼ 100 M ⊙ an optical depth of τ es ∼ 10 − 20 (see Gladstone et al. 2009) would need a scattering region as large as R ∼ 2 × 10 9 m (∼ 7 × 10 3 gravitational radii), comparable to the binary separation of Cyg X-1. This analysis suggests electron scattering is not a viable mechanism for suppressing variability in ULXs on such long timescales, although this is still a promising model for the energy spectra which require a much smaller emission region than needed here. Of course, it remains plausible that ULXs are examples of an accretion state that is simply not found (or very rare) among the commonly observed BHBs and AGN, and which displays very stable X-ray emission. Nevertheless, any mechanism for explaining the lack of variability in some ULXs faces a challenge: the X-ray spectra of the variable (e.g. NGC 5408 X-1) and nonvariable (NGC 5204 X-1) ULXs are quite similar Roberts et al. 2006;Stobbart et al. 2006) while their variability properties appear very different (see Fig. 5). This apparent inconsistency requires further study to understand fully, multiple observations of a few ULXs could allow for further clarification of this problem.
Longer observations with XMM-Newton would allow the frequency range to be extended to lower frequencies and could reveal in more detail the PSDs of ULXs such as NGC 55 ULX and NGC 1313 X-2, allowing better comparisons with BHBs and AGN. Observations of ULXs with more powerful X-ray missions (such as Astro-H and the proposed International X-ray Observatory) are needed before ULX variability can be studied at the higher frequencies and in the higher energy bands needed to test the hypotheses discussed above.