A spectral-timing model for ULXs in the supercritical regime

Abstract

1 INTRODUCTION

Ultraluminous X-ray sources (ULXs) have been widely observed in the local Universe, with inferred isotropic luminosities above 10³⁹ erg s⁻¹ (Roberts 2007; Feng & Soria 2011). Those below ∼3 × 10³⁹ erg s⁻¹ can be readily associated with accretion on to stellar mass black holes (BHs) (∼10 M_⊙) accreting close to or at their Eddington limit (see Sutton, Roberts & Middleton 2013, and references therein). There is now strong evidence to support this assertion, with the discovery of extremely bright ballistic jets from a ULX in M31 (Middleton et al. 2013; Middleton, Miller-Jones & Fender 2014b), which unambiguously links the flow with Eddington rate accretion (Fender, Belloni & Gallo 2004), and the first dynamical mass measurement of the compact object in a ULX, from M101 ULX-1 (Liu et al. 2013). Observations of such ‘low-luminosity’ ULXs (Middleton et al. 2011b, 2012; Kaur et al. 2012; Soria et al. 2012) have revealed changes in the disc emission that may imply the creation of a radiation pressure supported, larger scaleheight flow in the inner regions (Middleton et al. 2012) or magnetic pressure support (Straub, Done & Middleton 2013). Although emission below ∼2 keV is generally heavily photoelectrically absorbed by material in the Galactic plane (e.g. Zimmermann et al. 2001), similar spectral behaviour may also be seen in a small number of Galactic BH X-ray binaries (BHBs) at high rates of accretion (e.g. Ueda, Yamaoka & Remillard 2009; Uttley & Klein-Wolt, in preparation). Such ‘extreme’ high state BHBs probably dominate the ULX population (Walton et al. 2011), yet a significant number of ULXs can still be found at higher luminosities. Those above 10⁴¹ erg s⁻¹ are dubbed hyperluminous X-ray sources (HLXs; Gao et al. 2003) and provide the best evidence (Farrell et al. 2009; Davis et al. 2011; Servillat et al. 2011; Webb et al. 2012) for a population of intermediate-mass BHs (IMBHs; Colbert & Mushotzky 1999). Such IMBHs (with masses above those expected from direct stellar collapse: >100s of M_⊙) could potentially be formed in globular clusters (Miller & Hamilton 2002, but see Maccarone et al. 2007), through capturing and tidally stripping a dwarf galaxy (King & Dehnen 2005) or mergers in young super star clusters (Portegies-Zwart, McMillan & Gerhard 2003; Portegies-Zwart et al. 2004).

ULXs that fall between the two categories, i.e. L_x = ∼3 × 10³⁹ to 1 × 10⁴¹ erg s⁻¹) remain contentious and have been proposed as possible locations for IMBHs accreting at low sub-Eddington rates (e.g. Miller, Fabian & Miller 2004; Strohmayer & Mushotzky 2009). Indeed, the brightest objects in this class, with L_{X, peak} > 5 × 10⁴⁰ erg s^{− 1}, have demonstrated observational properties consistent with IMBHs in the hard state (Sutton et al. 2012). However, for the less-luminous ULXs, several problems exist with this interpretation for the entire population (see King 2004 for a discussion of theoretical issues related to formation); namely the emission characteristics do not generally match the expectation for low rates of accretion where, in the case of a BH of mass 10²⁻⁵ M_⊙, the emission from the disc would still peak in the soft X-ray band and so the structure of the flow is not expected to deviate strongly from that in BHBs at similar rates (see Zdziarski et al. 1998; Remillard & McClintock 2006). As a result, we would expect such sources to display a hard spectrum up to >50 keV due to thermal Comptonization in an electron plasma arranged in some still-undetermined geometry. Instead, ULXs up to 2 × 10⁴⁰ erg s⁻¹ generally show spectra that cannot easily be reconciled with sub-Eddington accretion (Stobbart, Roberts & Wilms 2006; Gladstone, Roberts & Done 2009) showing a spectral break above ∼3 keV (recently confirmed by NuSTAR observations of a sample of luminous ULXs; Bachetti et al. 2013; Walton et al. 2013b, 2014; Rana et al. 2014) which, in at least one source, has been unambiguously associated with Eddington-rate accretion (Motch et al. 2014).

Obtaining a deeper and full understanding of the nature of these sources, requires consideration of both spectral and variability properties simultaneously, with the latter providing a complementary set of powerful diagnostics by which to make comparisons to better understood sources. This has proven to be valuable, e.g. whilst the presence of quasi-periodic oscillations (QPOs) in the light curves of five ULXs (to date) have been used as evidence in support of IMBHs (Strohmayer & Mushotzky 2003, 2009; Strohmayer et al. 2007; Rao, Feng & Kaaret 2010), the details of the variability do not generally appear to well-match this identification (Middleton et al. 2011a; Pasham & Strohmayer 2012, but see also Pasham, Strohmayer & Mushotsky 2014). Notably, the recent study of Sutton et al. (2013) has reinforced the idea of using variability properties (the fractional variability; Edelson et al. 2002) together with spectra to broadly characterize the properties of ULXs. This has demonstrated an apparent dependence of variability on spectral shape, distinctly unlike that expected from IMBH accretion. Instead, it has been argued that both the spectral (Feng & Kaaret 2009; Gladstone et al. 2009) and variability properties (Heil, Vaughan & Roberts 2009) of these contentious ULXs can be fully explained by a model for ‘supercritical’ accretion on to stellar mass BHs (or equally neutron stars; King 2009; Bachetti et al. 2014), where inclination and mass accretion rate are likely to be the key determining factors in appearance (Middleton et al. 2011a, 2014a; Sutton et al. 2013). In most supercritical accretion models, a large scaleheight, optically thick equatorial wind is predicted (and reproduced in magnetohydrodynamic, MHD, simulations – Ohsuga 2007; Ohsuga & Mineshige 2011) to be radiatively driven from the disc from within the ‘spherization radius’, R_sph which can be at large radii (depending on the mass transfer rate from the secondary; Shakura & Sunyaev 1973; King 2004; Poutanen et al. 2007 – hereafter P07). Unlike models for sub-Eddington accretion on to IMBHs, this model can reproduce the observed spectrum with the soft emission being associated with the wind (King & Pounds 2003; P07) and the hard emission originating in the innermost regions which have been stripped of material, revealing the hot disc underneath (potentially further distorted by turbulence, advection, spin and self-heating; Beloborodov 1998; Suleimanov et al. 2002; Kawaguchi 2003). As the wind is expected to be optically thick (τ ≥ 1: P07, at least near to the disc plane), the scattering probability is large such that, depending on inclination angle, we may expect geometrical beaming to amplify the hard emission (King 2009; P07) or scatter emission out of the line of sight. Should the wind be inhomogeneous/clumpy (Takeuchi, Ohsuga & Mineshige 2013, 2014), this scattering can theoretically imprint variability by extrinsic means (Middleton et al. 2011a) and produce the large fractional rms seen in many sources (Heil et al. 2009; Sutton et al. 2013). However, an explanation for how this variability mechanism operates and depends on key system parameters, such as inclination and mass accretion rate, has been lacking.

In this paper, we build on previous theory and recent spectral-timing analyses of ULXs (notably Sutton et al. 2013) and present a simple model for the variability originating in a radially propagating, inhomogeneous wind. This allows us to make a series of key predictions for the evolution with mass accretion rate (Sections 2 and 3) which we compare to observations (Sections 4, 5 and 6).

2 THE SUPERCRITICAL MODEL OF ACCRETION

For the benefit of the reader, we now summarize the key properties of the supercritical model for accretion, as discussed and developed by several key authors. Models describing the supercritical inflow (Shakura & Sunyaev 1973; King et al. 2001; King 2004, 2009, P07; Dotan & Shaviv 2011) differ in their precise details but the overall picture is one where the high mass transfer rate (⁠|$\dot{m}_{\rm 0}$| – in dimensionless units of Eddington accretion rate) from a close binary system results in the Eddington limit being reached at large radii. In the case where the mass transfer rate is much greater than the ‘critical’ value of (9/4) × Eddington (P07), an optically thick outflow (with the escape velocity, |$v_{\rm esc} > \sqrt{2GM/R}$|⁠, where M is the BH mass and R is the radial distance) can be launched from within R_sph ≈ |$\dot{m}_{\rm 0}R_{\rm in}$| (where R_in is the position of the inner edge in R_g = GM/c²) where the scaleheight of the disc (proportional to |$\dot{m}_{\rm 0}$|⁠; Shakura & Sunyaev 1973) exceeds unity. Mass is then lost through this outflow with the mass accretion rate decaying approximately linearly down to R_in, where the accretion rate is locally Eddington (in the limit of no advection; P07).

The geometry of the supercritical inflow can be broadly defined by three zones (depending on the optical depth of the outflow), as fully described in P07 (to which we direct the interested reader) and illustrated in Fig. 1.

• In the innermost regions (R < R_ph,in), the wind is essentially transparent (τ ≤ 1) and so the emission will appear as a distorted ‘hot disc’, peaking at a characteristic temperature, T_in (modified by some colour temperature correction, spin and overheating; see Beloborodov 1998; Suleimanov et al. 2002; Kawaguchi 2003).

• At larger radii, the outflow – seen effectively as an extension to the large scaleheight inflow – is expected to be optically thick for mass accretion rates above a few times Eddington. As a result of advection in this optically thick material (as the photon diffusion time-scale is longer than the viscous time-scale), the radial temperature profile (as seen at the last scattering surface at τ = 1) is broadened from R^−3/4 to R^−1/2 (see Abramowicz et al. 1988), and the result is a smeared blackbody extending from T_ph,in to T_sph. As the material in the wind is outflowing and we have assumed this is an extension of the inflow, the viscously dissipated energy emerges at larger radii (∼ twice the radius at which it is generated; P07). However, this does not affect the temperature profile (and nor in principal does having the wind disconnected from the inflow due to stratification).

• At radii greater than R_sph, the optical depth of the wind falls as 1/r (P07) such that the underlying emission begins to emerge at approximately the radius at which it is generated, i.e. the emission is a quasi-thermal blackbody with a peak temperature of T_ph, which we shall refer to hereafter as the ‘outer photosphere’ (and for high mass transfer rates will peak in the UV rather than X-ray band – equation 38 of P07).

Figure 1.

Schematic representation of the supercritical model for bright (>3 × 10³⁹ erg s⁻¹) ULXs (P07) where an optically thick wind is launched from the surface of a large scaleheight flow. This flow is driven from the disc from R_sph (≈|$\dot{m}_{\rm 0}R_{\rm in}$|⁠; Shakura & Sunyaev 1973) down to R_ph,in within which the outflow is effectively transparent and the underlying disc emission can emerge (which may be further distorted by advection, spin and overheating; P07, Beloborodov 1998; Suleimanov et al. 2002; Kawaguchi 2003). Variability due to obscuration/scattering can be introduced when the wind is inhomogeneous – expected on large scales via radial propagation of flux (Lyubarskii 1997) and seen on small scales in 2D MHD simulations (Takeuchi et al. 2013, 2014). Depending on the orientation of the observer (given by θ), the emission can be stochastically boosted by scattering into the line of sight (position 1) or reduced by scattering out of the line of sight (position 2). The spectrum will also depend heavily on θ, with the emission from the hot inner region becoming progressively geometrically beamed with smaller θ and Compton down-scattered at larger θ.

Open in new tabDownload slide

2.1 Effect of a static wind

The detailed physical properties of the wind require the global geometry and radiative transfer (e.g. Abolmasov, Karpov & Kotani 2009) to be fully modelled. Whilst beyond the scope of this work, we can already make some simple deductions. The large scaleheight flow and wind will naturally subtend a large solid angle to the hottest inner regions (Shakura & Sunyaev 1973; P07) and so a correspondingly large fraction of the high-energy flux will illuminate the material. The latest 2D MHD simulations (Takeuchi et al. 2013, 2014) show that the optically thick wind will have inhomogeneities due to the Rayleigh–Taylor (or other hydrodynamic instability) with a size of the order of 10 R_s (R_s = 2GM/c²). Given typical particle densities expected for these clumps (up to ∼10¹⁷ cm⁻³; Middleton et al. 2014a), the wind will have a highly photoionized skin (logξ > 6 out to hundreds of R_g for a 10 M_⊙ BH) down to an optical depth of unity (where most of the scatterings occur). As a result, we expect the surface material of the wind (clumpy or otherwise) in the inner regions to be approximated to an electron plasma with a large optical depth (τ ≥ 1 – see P07) such that large amounts of incoming flux from the inner regions will be scattered into the cone of the wind. This naturally results in geometrical beaming (proportional to |$\dot{m}_{\rm 0}^2$|⁠; King & Puchnarewicz 2002; King 2009) where, from simple inspection of the assumed geometry (a hemispherical wind of opening angle θ_w, see Fig. 1), we expect a scattered fraction of S_f ≈ cos θ_w/(1 − cos θ_w). As an example, for a wind launching angle of 45° (Ohsuga & Mineshige 2011), we would then expect a factor of ∼2.4 more flux to be scattered to an observer viewing down the cone of the wind.

Scattering will not only amplify the emission from the innermost regions when looking into the cone of the wind but also change the energy of the scattered photons: as the wind plasma temperature is less than that of the hot disc, i.e. hν > kT_e, the highly energetic incident photons will lose energy to the plasma in the inflow via Compton down-scattering (Sunyaev & Titarchuk 1980) and to the outflowing wind via both recoil and bulk effects (e.g. Titarchuk & Shrader 2005; Laurent & Titarchuk 2007). At larger inclinations (i.e. sight lines intercepting the wind), the effect of beaming diminishes and the number of scatterings towards the observer – and therefore energy loss – will increase (as the optical depth to the observer will be larger), although this will be complicated by the density profile of the wind which will not be constant as radiation momentum transfer (plus the lack of hydrostatic equilibrium) will cause the initial clumps to expand and density drop as they move away from the launching point (e.g. Takeuchi et al. 2013).

2.2 On the origin of variability in ULXs

The toy model of Middleton et al. (2011a) ascribed the observed large amplitude, short time-scale (tens to thousands of seconds) variability seen in NGC 5408 X-1 to obscuration by individual clumps, generated through radiative-hydrodynamic instabilities in the wind (now reproduced in simulations; Takeuchi et al. 2013, 2014). In the specific case of NGC 5408 X-1, this requires our line of sight to the inner regions to intercept the wind (see Middleton et al. 2014a). In general, however, we may expect the variability in ULXs to originate through two mechanisms and below we provide the framework for these.

2.2.1 Method 1: clumps

Should |$\dot{m}_{\rm 0}$| increase, we should expect some change in the filling factor of the clumps. An increase may tend to lower the amount of imprinted variability, but as the size and launching of individual clumps (on which the variability relies) is a function of the radiative-hydrodynamic instabilities (see Takeuchi et al. 2014), an accurate understanding is beyond the scope of the work here.

2.2.2 Method 2: propagating flux

The above description assumes a wind that is ‘steady state’ (i.e. matter is being launched at a constant rate at a given radius). However, in sub-Eddington inflows, the observed (often considerable) variability (e.g. Remillard & McClintock 2006; Belloni 2010; Muñoz-Darias, Motta & Belloni 2011; Heil, Uttley & Klein-Wolt 2014a,b) originates as a result of inwardly propagating variations in mass accretion rate/surface density through the viscous inflow (Lyubarskii 1997; Ingram & Done 2012) which leads to the universally observed rms–flux relation in all accreting sources (Uttley, McHardy & Vaughan 2005, Arévalo & Uttley 2006; Heil & Vaughan 2010; Heil, Vaughan & Uttley 2012; Scaringi et al. 2012). This can then lead to the observed power density spectrum (PDS) shape common to BHBs formed from a series of convolved Lorentzians, damped above the local viscous time-scale (Ingram & Done 2012).

If we assume that propagation of fluctuations still occurs in the larger scaleheight supercritical flow, then this underlying and inherent stochasticity will drive radial variations in the mass-loss in the form of the clumpy wind (see equation 26 of P07; Fig. 2). Should the variability originate via obscuration/scattering by clumps in the launched material, the propagating fluctuations will imprint additional variability on to the emission by changing the radially dependent ‘global’ amount of mass-loss and scattering. Additionally, should the inner disc not be geometrically thin, then the variability propagated down should also emerge directly. Such an origin would seem to argue against an IMBH interpretation where such discs should be geometrically thin and intrinsically stable (Churazov, Gilfanov & Revnivtsev 2001).

Figure 2.

Schematic slice of the inflow and wind. At large radii close to R_sph, the wind loss is affected by propagated variability (Lyubarskii 1997) through the sub-critical disc into the large scaleheight inflow. As the mass-loss does not conserve the variability but instead removes it from the flow, the next inner radius has less with which to convolve (vertical arrows). Therefore, within R_sph (≈ |$\dot{m}_{\rm 0}R_{\rm in}$|⁠), the propagated variability (which can emerge from both the hot inner disc if not geometrically thin and via scattering of this radiation) is damped. Scatterings by individual clumps (Method 1) – which can imprint an extra source of variability on relatively short physical time-scales – will average out for a large number of elements, i.e. at small θ, but at larger θ an imprint is expected to remain.

Open in new tabDownload slide

To determine the likely time-scales imprinted by this process, we can consider the viscous time-scale at R_sph for representative values, e.g. H/R (disc scaleheight) ≈1, α (viscosity parameter) = 0.01, |$\dot{m}_{\rm 0}$| = 1000–10 000 (e.g. as expected for an SS433 type system; Fuchs, Koch-Miramond & Ábrahám 2003) and for a 10–100 M_⊙ BH, ν_sph ≈ 0.1–0.001 mHz. This is at the edge or beyond the typically observable bandpass (due to current X-ray count rate limitations) and so we should not expect large amounts of ‘rapid’ variability (e.g. tens to hundreds of seconds) to be imprinted by this process as a result of dampening above ν_sph. However, it is clear that even where R_sph is large, the time-scales of both methods could, in principle, overlap.

Importantly, as opposed to the contribution from single wind elements (Method 1), at small inclinations the imprinted non-Gaussian variability via propagated fluctuations will not average out. The total variability for a single ULX will then be a combination of the two methods – although at small inclinations we should only expect variability by this second process – and a PDS resulting from the combination. Although it may be tempting to produce test PDS from the combination of the two methods (specifically for those ULXs at moderate inclinations), we caution that the relative normalizations are unknown and doing so would lead to misleading results. In spite of the unknowns, we can make predictions resulting from variations in |$\dot{m}_{\rm 0}$| within a combined model of the spectrum and variability as presented in the following subsection and Section 3.

[For clarity, we reiterate that, as opposed to the origin of variability via Method 1, variability via propagating flux variations through the supercritical inflow can imprint on to the emission via the hot disc directly only if the scaleheight is not small (i.e. not a thin disc) and/or via scattering events from the global changes in radial mass-loss in the wind (which does not rely upon the inner disc being intrinsically variable).]

2.2.3 Effect of increasing |$\dot{m}_{\rm 0}$| on Method 2

Note in the above equation that the index would usually be −γ but as the PDS is moving to lower frequencies with increasing |$\dot{m}_{\rm 0}$|⁠, across a fixed bandpass, the ratio of observed powers is inverted. To help illustrate this point, we have shown the shift in PDS with increasing R_sph in Fig. 3.

Figure 3.

Illustration of the shift of the PDS to lower frequencies associated with an increase in |$\dot{m}_{\rm 0}$| and R_sph. Due to the time-scales at R_sph, we expect to see the diluted power above the break (at R_sph) where P ∝ ν^−γ. As the power shifts out of the observable bandpass (shown as the shaded area), the drop in power is equivalent to the increase in power that would have occurred had the bandpass not been fixed. As a result of the shift, we then find ΔlogP = γΔlogν or (P₂/P₁) = (ν₂/ν₁)^γ in the formulae below, where we use the ratio of the viscous frequencies at each R_sph in the derivations for ease (and the ready connection to the change in |$\dot{m}_{\rm 0}$|⁠).

Open in new tabDownload slide

We note that, in the above case, γ is intuitively also a function of |$\dot{m}_{\rm 0}$| as the effect of dampening is expected to scale with |$\dot{m}_{\rm w}$| (see equation 26 of P07). Whilst there are several assumptions, it is clear that the general trend between spectral hardness and power will be an anticorrelation when considering only Method 2, i.e. specifically for those ULXs viewed at smaller inclinations.

As a final caveat, we note that we have only considered the variability to be a modulation of the high-energy emission from the hot inner disc (via scattering and/or directly); however, stochastic variability via Method 2 may also lead to some variability of the intrinsic wind emission (the likely advection-dominated flow from R_sph to R_ph,in) by varying the column density of material and thereby changing the colour temperature correction (f_col). Although this too will tend to average out at small θ for individual elements (Method 1), the longer time-scale trends introduced by the radial inhomogeneities will leave a global imprint which, once again, should not average out. However, assuming that the soft emission has the same physical origin in all ULXs, then the lack of variability observed in NGC 5408 X-1 at these energies (Middleton et al. 2011a), implies that such variability is likely to be relatively weak.

3 A COMBINED SPECTRAL-TIMING MODEL

In order to make a set of predictions for how the properties of ULXs should evolve, we must consider both of the methods discussed above and the spectral-variability patterns we should expect to result from changes in |$\dot{m}_{\rm 0}$|⁠. The predictions which follow are, by necessity, only semiquantitative; a full and accurate quantitative picture can only be obtained from full radiative simulations which include the nature and impact of instabilities on the short time-scale variability (and are beyond the scope of this work).

For the following predictions, we assume that the inflow is supercritical (or is responding in the manner of being such; Bisnovatyi-Kogan & Blinnikov 1977) with the wind launched from the disc such that increasing |$\dot{m}_{\rm 0}$| will still increase the scaleheight of the inflow further. We also assume that emission from the subcritical disc beyond the wind is negligible (i.e. |$\dot{m}_{\rm 0}$| is large enough that emission beyond R_sph is out of the X-ray bandpass).

3.1 Source/population evolution: spectral hardness versus variability

The inclination of the observer's line of sight and how this intercepts the wind (θ in Fig. 1) is pivotal to the observed spectrum and power. We therefore dispense with past descriptive terms for ULXs based on the spectral shape alone and instead present an inclination-dependent description based on the line of sight indicated in Fig. 1.

Sources at small θ.

At small inclinations, where the wind does not enter the line of sight (position 1 in Fig. 1), the observer sees the maximum unobscured emission from the hot inner ‘disc’ as well as the scattered flux from the cone of the wind. We therefore expect a spectrum with a strong hard component with a beaming factor scaling as |$\dot{m}_{\rm 0}^2$| (King & Puchnarewicz 2002; King 2009). Such a spectrum would correspond to the ‘hard ultraluminous’ class in Sutton et al. (2013).

At low inclinations, we expect variability only through Method 2 with the variability emerging directly in the disc (if not geometrically thin) and/or via scattering of hot photons from the inner regions into the line of sight. Given the nature of the plasma, some proportion of the variability should be shifted to energies below the high-energy peak via down-scattering in the surface plasma of the wind (Titarchuk & Shrader 2005). We may also see variability on similar time-scales in the intrinsic wind emission possibly due to the changing radial density profile and therefore f_col, although as scattering events are predicted to dominate the emission (via beaming), we might not expect this component to dominate the fractional variability (although we may still detect its presence in an absolute variability spectrum).

At these inclinations, an increase in |$\dot{m}_{\rm 0}$|⁠, which leads to a larger R_sph and also smaller wind cone opening angle (θ_w in Fig 1, as scaleheight scales with |$\dot{m}_{\rm 0}$|⁠; Shakura & Sunyaev 1973), can lead to one of two possibilities (assuming no system precession).

• The wind remains out of the line of sight such that the spectrum will get increasingly beamed (i.e. brighter and harder) and the variability will drop as discussed in Section 2.2.3. The intrinsic shape of the hot component should vary little (P07).

• The wind enters our line of sight such that the spectrum becomes softer (as hard emission is both scattered away and down-scattered through the wind).

In both cases, as R_sph increases with |$\dot{m}_{\rm 0}$|⁠, we might expect a corresponding decrease in T_sph (equation 37 of P07), leading to a predicted anticorrelation with luminosity (see P07; King 2009). However, contributions to the soft emission by down-scattering and advection will complicate the evolution and this remains an important, unresolved issue (see Miller et al. 2013, 2014).

Sources at moderate θ.

At larger θ, the wind will start to enter the line of sight (position 2 in Fig. 1). As a result, a significant fraction of the hard X-ray flux will be scattered out of the line of sight such that the spectrum has a relatively smaller contribution from the hot inner regions, and that which does arrive is expected to be down-scattered such that the peak temperature of the hot component is cooler than when seen more directly, i.e. at small θ. Such a spectrum would correspond to the ‘soft ultraluminous’ class in Sutton et al. (2013).

As opposed to ULXs at small θ, the variability can originate by both methods presented in Section 2 (including possible changes in f). As the variability we can readily observe should have a large contribution from obscuration events (Method 1), it will likely peak at the energy of the (down-scattered) hot inner disc. As the fraction of flux towards the observer is lower, the integrated fractional variability (and therefore power) should be larger than for smaller θ.

We note that whilst it may appear incongruous to have both obscuration and down-scattering, given the covering fraction of the wind to the hot inner regions, it seems inescapable that down-scattering will occur in some less optically thick phase of this wind (we speculate that this could occur in the expanded clumps further from the launching point).

Unlike those ULXs viewed at small inclinations, we expect only one outcome from increasing |$\dot{m}_{\rm 0}$| (once again assuming no system precession).

• As we expect the wind to increasingly dominate our sight lines, more of the hard emission is beamed out of the line of sight whilst that in the observer's direction is increasingly Compton down-scattered. As a result, emission from the soft component should increase by a larger amount than the hard component. Variability via Method 2 will be increasingly damped, whilst the effect on Method 1 (clumps) is somewhat unclear (and will be a function of the generation of the instabilities); however, we predict that the filling factor will likely increase with increasing |$\dot{m}_{\rm 0}$| and so we might expect a corresponding decrease in this contribution as the variability will increasingly tend to average out. Irrespective of the method, any variable, scattered emission from the hot disc must also pass through the larger optical depth inflow/outflow where reprocessing will redistribute it to lower energies (as with the inner disc emission) and, as we assume the soft emission from the wind to be relatively stable, dilute the fractional imprint (and therefore power). The most extreme result for this evolution is the spectrum and variability tending towards those predicted for largest θ; as a corollary, it is a distinct possibility that the luminosity could drop below that considered to indicate a ‘bright’ ULX (>3 × 10³⁹ erg s⁻¹).

Sources at large θ.

In the case, where the observer is viewing at largest θ (position 3 in Fig. 1), the outer photosphere of the wind at R > R_sph, is predicted to dominate the observed emission (at T_ph; P07), which may be out of the X-ray bandpass and emit in the UV. Observing at such inclinations may then explain the extremely bright UV emission seen in ultraluminous UV sources (e.g. Kaaret et al. 2010), the supersoft ULX, M101 (Kong & Di Stefano 2005; Shen et al. 2015), and the UV excess of SS433 (Dolan et al. 1997). As the X-ray luminosity is expected to be low, the source may not fall into the empirical class of bright or even faint ULXs (see as possible examples Soria et al. 2010, 2014). Whilst we may not expect these to feature in ULX population studies, they represent an important set of ‘hidden’ ULXs, should provide interesting diagnostics for the winds (Middleton & Maccarone, in preparation) and are a necessary component of any model which hopes to describe the entire population. Although these sources are unlikely to feature in X-ray studies (although over time we may be able to see some sources evolve to being UV bright), it is useful to make predictions about the variability properties and effect of changing |$\dot{m}_{\rm 0}$|⁠. Due to the large optical depth in this direction (P07), we should not expect contributions to the variability from that emerging directly from the inner regions or via scattering and propagating through the inflow (Method 2). However, we should expect variability on the long viscous time-scales at R_sph. For an observer already at large θ, increasing |$\dot{m}_{\rm 0}$| will increase the soft emission only.

Summary.

Here, we summarize the predicted spectrum and variability properties for a source viewed at a given inclination to the wind.

• Small θ: the spectrum is expected to be hard with variability via Method 2 only. If |$\dot{m}_{\rm 0}$| increases, then the spectrum should get harder with a corresponding drop in variability. Should the wind enter the line o sight, then the spectrum may soften and, assuming |$\dot{m}_{\rm 0}$| is not too high to make the wind homogeneous, the variability may increase – analogous to viewing at moderate θ.

• Moderate θ: the spectrum is expected to be soft with variability via Method 1 and 2 but likely dominated by Method 1 (although we re-stress that the relative normalizations are unknown at this time). If |$\dot{m}_{\rm 0}$| increases, then the spectrum should get softer and the variability should decrease, both due to the wind tending towards homogeneity and the variability being suppressed by the stable soft component.

• Large θ: the spectrum should be extremely soft and possibly out of the X-ray bandpass altogether. Long time-scale variability only should be present.

A highly simplified version of the predicted evolution for the population is presented in Fig. 4 with spectral hardness (in an observing band which includes both the soft and hard components) plotted against fractional rms squared (in log space and under the reasonable assumption that the two methods of imprinting variability overlap in time-scales: see Section 2). The overall shape is an inflection (positive to negative gradient with spectral hardness), peaking where the sources are soft and highly variable, corresponding to those seen at moderate inclinations. From this inflection point, the variability drops as the spectrum hardens (dominated by those face-on sources: see equation 14 for a derivation of this in a single source) and drops towards softer spectral colours (as the wind increasingly dominates the emission). The grey points track the possible evolution of the spectrum softening out of a hard state whilst the variability does not substantially increase (or even drops), as a result of the wind tending towards homogeneity as |$\dot{m}_{\rm 0}$| increases. In practice, this will lead to increased scatter (depending on the likelihood of the situation) around the inflected path. We note that the predicted shape is distinctly unlike that expected from scaling the hardness-variability evolution of BHBs (Belloni 2010; Muñoz-Darias et al. 2011) to IMBH masses where M_BH > 100 s of M_⊙ would imply sub-Eddington states and increasing rms with spectral hardness.

Figure 4.

Simple schematics showing the possible evolutionary paths of a ULX in spectral hardness (in an observing band which includes both the soft and hard components) and power (fractional rms squared), which depend on the inclination the wind makes to the observer and the change in |$\dot{m}_{\rm 0}$| as described in section 3. Left: the evolution without precession is shown. When viewed at low inclinations to the wind (i.e. when the wind is not directly in the line of sight), the source is predicted to get harder with increasing |$\dot{m}_{\rm 0}$|⁠, with a drop in variability due to larger R_sph and increased dampening in a fixed frequency bandpass (giving a gradient in log space of ≈−3γ/4 – equation 14). Alternatively, should the wind start to intercept our line of sight, the spectrum will become softer due to beaming of hard emission away from, and down-scattering towards, the observer, with variability increasing via Method 1 (where N elements contributing to such events is lower than at small inclinations). Eventually, the wind will dominate, leading to substantial Compton down-scattering and a softer, less variable spectrum. The grey points indicate the evolution where an increase in mass accretion rate results in the wind tending towards homogeneity and suppressing variability via Method 1. Centre: the evolution with increasing |$\dot{m}_{\rm 0}$| including the effect of precession which can further change our view of the ULX. Right: the evolution including precession with decreasing |$\dot{m}_{\rm 0}$|⁠. Across a population, we would expect the combination of the possible tracks, weighted by likelihood. In effect, we expect this to be an inflected shape (peaking with those soft ULXs seen at moderate inclinations with large amounts of variability) with substantial scatter due to the effect of precession and wind tending towards homogeneity at larger |$\dot{m}_{\rm 0}$| (suppressing variability in Method 1). Importantly, above the softest, most variable sources we would expect a negative gradient, distinctly unlike the positive gradient expected from BHBs scaled to IMBH masses (see Belloni 2010; Muñoz-Darias et al. 2011).

Open in new tabDownload slide

3.2 Source/population evolution: effect of precession

In the previous subsection, we made predictions for how a source's spectral and variability properties should appear for a given inclination and how they might change in response to increasing |$\dot{m}_{\rm 0}$|without precession. However, precession – as seen in SS433 (see Fabrika 2004 for a review) – due to a ‘slaved-disc’, i.e. one where the donor star's rotational axis is inclined to its orbital axis (Shakura 1972; Roberts 1974; van den Heuvel, Ostriker & Petterson 1980; Whitmire & Matese 1980), must lead to changes in the spectral and variability properties. In the simplest case, where |$\dot{m}_{\rm 0}$| remains unchanged as the disc precesses, we should see evolution between each of the inclinations minus the effect of increased accretion rate on homogeneity of wind and position of R_sph. Changes in |$\dot{m}_{\rm 0}$|as well as precession would lead to greater complexity but can be deduced from the discussions above.

• Low→moderate θ with increasing |$\dot{m}_{\rm 0}$|⁠: as discussed in the subsection above, if the new |$\dot{m}_{\rm 0}$| is high enough to suppress Method 1 then the variability could drop whilst the spectrum softens. Otherwise, the variability will increase as the spectrum softens. The reverse evolution will occur for moderate→low θ with decreasing |$\dot{m}_{\rm 0}$|⁠.

• Low→moderate θ with decreasing |$\dot{m}_{\rm 0}$|⁠: variability increases as the spectrum softens. The reverse evolution will occur for moderate→low θ with increasing |$\dot{m}_{\rm 0}$|⁠.

• Moderate→high θ with increasing |$\dot{m}_{\rm 0}$|⁠: the spectrum will soften and the variability will decrease (and start to shift to longer time-scales).

• Moderate→high θ with decreasing |$\dot{m}_{\rm 0}$|⁠: the spectrum will soften and variability is also likely to decrease (or decrease with a flatter gradient).

We plot the possibilities above for both increasing and decreasing |$\dot{m}_{\rm 0}$| in Fig. 4. There are clearly several paths a given ULX may take depending on the presence of precession, |$\dot{m}_{\rm 0}$| and starting inclination, with the overall trend we should expect to observe across the population resulting from the combination of the possibilities weighted against their likelihood. Whilst precession will no doubt distort observational trends, by studying a large number of sources we can expect to average across the effect and still see an inflected shape across the population.

4 OBSERVATIONAL COMPARISON

The work of Sutton et al. (2013), goes some way to establishing a trend of variability with spectral hardness which (encouragingly) appears similar to our predictions. However, to rigourously test our model predictions and identify weaknesses, it is important to examine the spectral and variability evolution of a high-quality sample of ULXs – excluding those which can now be confidently associated with ‘normal’ modes of BHB evolution (Middleton et al. 2013) – in greater detail.

We select nine of the brightest ULXs (in flux), known to span a range of spectral shapes (Gladstone et al. 2009), that have been well observed by XMM–Newton to ensure the availability of high-quality data sets for spectral and timing analysis. Importantly, we exclude the ‘broadened disc’ class of sources from Sutton et al. (2013) which are considered to be BHBs experiencing normal mass transfer rates but with Eddington inflows in their most inner regions (Middleton et al. 2011b, 2012, 2013). The observational details of our sample are given in Table 1.

For all observations, we re-process the data using sas v.12.0.1 and up-to-date calibration files. We apply standard data filters and flags for bad pixels and patterns (see the XMM–Newton user's handbook¹), and remove periods of high-energy background (based on the 10–15 keV count rate from the full field of view) associated with soft proton flaring.

We use xselect to extract EPIC-PN spectra and light curves from source regions not smaller than 30 arcsec radius (with the exception of the first six observations of NGC 6946 X-1 as the source is close to a chip-gap) and background regions chosen to avoid the pixel read-out direction and other sources in the field. We do not use the EPIC-MOS in this analysis as the fractional increase in data is small, and the instrumental responses differ at soft energies.

4.1 Spectral analysis

The 0.3–10 keV spectra of ULXs above ∼3 × 10³⁹ erg s⁻¹ generally favour a two-component model description (Sutton et al. 2013) with a high-energy (>3 keV) break, distinctly unlike the spectra observed from sub-Eddington accretion (Zdziarski et al. 1998; Stobbart et al. 2006; Gladstone et al. 2009; Bachetti et al. 2013). In this section, we model the ULX spectra assuming the soft emission originates in the wind and hard emission from the inner, distorted disc with down-scattering likely broadening the emission towards energies below the peak.

Although advection is likely to be important for the wind at large |$\dot{m}_{\rm 0}$|⁠, we note that even the highest quality data sets do not presently favour describing the soft emission with a model for an advection-dominated disc over a simple thin disc (though see Section 5). As a result, when fitting the data in xspec v 12.8 (Arnaud 1996), we use a quasi-thermal, multicolour disc blackbody (DISKBB; Mitsuda et al. 1984) to account for emission from the wind (i.e. the emission from R_ph,in to R_sph; e.g. Kajava & Poutanen 2009) and a broad model to describe the emission from the inner disc and its down-scattered component. NTHCOMP is an appropriate model for the latter as it can provide a variety of complex continuum shapes (Życki, Done & Smith 1999) whilst providing some useful and readily identifiable ‘physical’ parameters (e.g. the photon index of the spectrum, Γ, and high-energy rollover: 2–3 kT_e) although we have to be careful in how we interpret these in the sense of the physical model. We note that Walton et al. (2014) show that the high-energy tail (in at least Ho IX X-1 and Ho II X-1) following the rollover is not well described by a Wien tail. As this is out of our bandpass, it will not affect our spectral modelling (as our model hinges on the characteristic temperatures) although we speculate that this broadening may result from the heavily distorted inner disc emission (see also Tao & Blaes 2013).

We combine these emission components with a model for neutral absorption (TBABS) using appropriate abundance tables (Wilms et al. 2000) with the lower limit set at the Galactic column in the line of sight to each host galaxy (Dickey & Lockmann 1990). We fit this simple convolved model (TBABS*(DISKBB+NTHCOMP)) to all of the ULX data sets, keeping the high-energy component within a region of parameter space such that the high-energy rollover (which here refers to the inner disc peak temperature, i.e. 2–3kT_e ≈ kT_in) is in the observable bandpass (as seen in Stobbart et al. 2006; Gladstone et al. 2009). This is preferable to allowing the model parameters to roam freely which can often lead to the high-energy rollover being out of the bandpass, which in light of our model would be incorrect (i.e. we do not expect a further break in the spectrum, as indicated by NuSTAR observations; Bachetti et al. 2013; Walton et al. 2014). As a result, the index (Γ) is associated with the rise to the peak of the high-energy emission and allows for a contribution from down-scattering.

In using the NTHCOMP component to model emission from the hot inner disc and down-scattered component, we fix the seed photon temperature to be that of the wind to provide a realistic (albeit conservative) constraint: by energy balance, we expect the down-scattering to lead to temperatures ≳T_sph. For the highest quality (or softest) data sets, we note that the models favour the wind emission component dominating at soft energies (≲1 keV) over the hard component. However, due to poor data quality in several data sets or where the data does not obviously show a two-component structure, it becomes necessary to force the presence of a soft component dominating below ≲1 keV rather than allowing a fit with only a single (broad) component. As a consequence, the errors on the parameters in these models (denoted by a * or ** in Tables 2 and 3) are strictly only an approximation as they are found by fixing the temperatures (kT_d and kT_e) and/or the slope (Γ) for each component in turn (and rarely by fixing the normalization of the DISKBB component).

Across all observations, we generally obtain a fairly hard (Γ < 2.4) photon index and a rollover in the spectrum >3 keV (as reported in Gladstone et al. 2009; Bachetti et al. 2013). The best-fitting parameters of interest from model fitting are presented in Tables 2 and 3 (with the total absorbed luminosities also provided in Table 1) and spectral plots in Figs 5, 6 and 7.

Figure 5.

Best-fitting continuum models (for model parameters, see Tables 2 and 3) for the ULXs as labelled. For those where we have significant variability (≥3σ), the spectral models are shown in green. For the brightest and dimmest spectra (in unabsorbed 0.3–10 keV flux, indicated by thicker solid lines), we also plot the model components, DISKBB: red solid/dashed (bright/faint) lines and NTHCOMP: blue solid/dashed (bright/faint) lines. This allows for crude inspection of how the ULX has evolved from dimmest to brightest (although care must be taken as this may be misleading in certain cases: see NGC 1313 X-1 for an example in the discussion).

Open in new tabDownload slide

Figure 6.

As for Fig. 5.

Open in new tabDownload slide

Figure 7.

As for Fig. 5.

Open in new tabDownload slide

In order to make comparisons to our predictions for the spectral-timing evolution of ULXs as a population, we obtain the de-absorbed fluxes in two energy bands: 0.3–1 and 1–10 keV, respectively, by including a CFLUX component in the model fitting (TBABS*CFLUX*(DISKBB+NTHCOMP)); these values are given alongside the model parameters in Tables 2 and 3. We then determine the ‘colour’ from the ratio of hard to soft fluxes and determine the error on the colour from propagation of errors. Together with the fractional rms (see next section), we can subsequently assess the evolutionary pattern of ULXs in spectral hardness and variability.

4.2 Timing analysis

We obtain the PDS (normalized to be in [σ/mean]² units) by fast-Fourier transforming segments of background-subtracted, 0.3–10 keV, time series of length 1200 s and taking the average of the resulting periodograms (see van der Klis 1989 for a review). Short periods (<10 s) of instrumental dropouts are then corrected by linearly interpolating between the points either side of the gap. By integrating the (Poisson noise subtracted) power over a given frequency range and taking the square root, we are naturally left with the fractional rms (F_var; see Edelson et al. 2002). This is a more robust method of measuring F_var than by obtaining the variance of the light curve in the time domain as multiple measurements of the variance are taken and the true distribution can be measured (see Vaughan et al. 2003). We integrate the PDS over a frequency range 3 to 200 mHz – as this range encapsulates the broad-band noise variability seen in most ULXs (Heil et al. 2009) and ensures adequate statistics in the majority of cases – with the resulting values presented in Tables 2 and 3. We note that, although QPOs are seen in two of the sources presented here (NGC 5408 X-1 and NGC 6946 X-1) and lie within the frequency range over which our rms is calculated, in neither case will the variability of the QPO dominate over the broad-band noise (see e.g. Strohmayer & Mushotzky 2009) and so our values fairly reflect the underlying variability.

Constraining a value for F_var depends on both the level of Poisson noise (set to be 2/mean count rate in our normalization) and the length of an observation; for those observations where the former is high or the latter short, we may not be able to directly constrain the presence of variability. For observations with unconstrained fractional rms, upper limits were found by performing light curve simulations with an underlying PSD shape equivalent to that of a highly variable observation of NGC 6946 X-1 (Obs.ID: 0691570101), determined by fitting Lorentzians to the PDS. Simulated light curves were generated with the same count rates and duration as the real observations. In order to place upper limits on the power (with a similar shape to that of NGC 6946 X-1), hidden within each PDS, the input model was rescaled at 1 per cent intervals from 20–100 per cent of the normalization of the original PDS. For each rescaled input model, 200 light curves were generated and the average rms and its error calculated; the upper limit was then taken to be the point where the rms was just significant at a 3σ level and is given in Tables 2 and 3.

We plot the power (i.e. fractional rms squared) for each observation (where constrained) against the spectral hardness in Fig. 8. In plotting these together, we have assumed that the variability time-scales between sources are the same (reasonable in light of expected viscous time-scales: see Section 2) and that we have a homogeneous population of compact objects. Although differences in the mass transfer rate will lead to differences in the dampening of variability (as γ in Method 2 is expected to be a function of |$\dot{m}_{\rm 0}$|⁠), by including a large enough sample of sources (and observations) we can start to average over this difference.

Figure 8.

Spectral hardness (from the ratio of unabsorbed 1–10/0.3–1 keV fluxes) versus power (fractional rms squared) in log space (with 1σ errors). The colour scheme corresponding to the different sources is as follows: black: NGC 5408 X-1, red: NGC 6946 X-1, purple: NGC 1313 X-2, orange: NGC 55 ULX-1, green: NGC 1313 X-1, cyan: IC 342 X-1, blue: Ho IX X-1, yellow: Ho II X-1. This appears similar in overall shape to our predicted evolution shown in Fig. 4.

Open in new tabDownload slide

Although there is considerable scatter, the general trend appears to be broadly consistent with the predicted shape of spectral hardness versus power shown in Fig. 4. The amount of scatter is unsurprising given the range of possible evolutionary tracks (driven by the effect of precession and the wind tending towards homogeneity at large |$\dot{m}_{\rm 0}$|⁠: Section 3) and that there will be a distribution of masses. We note that should this association be correct, the positive gradient branch would appear to be underpopulated; this is likely due to selection effects as these will be underluminous and may not qualify as bright ‘ULXs’ in the traditional sense (see Section 3).

To estimate the significance of a correlation on the negative slope, we exclude NGC 55 ULX-1 – as the spectrum is significantly softer than any other ULX (see next section for more details) – and perform a Kendall's rank coefficient test on the data binned into five even logarithmic spectral hardness bins (which does not account for the errors on the data); however, this is only marginally (2σ) significant. We can however, attempt to rule out a positive slope, expected should we be observing a population of IMBHs (which we assume behave as BHBs at similar Eddington ratios; Belloni 2010; Muñoz-Darias et al. 2011). We determine the gradient of the slope through a least-squares fit to the un-binned data points (using average symmetrical errors), constraining the gradient to be negative at > 3σ (<−0.14 at 3σ). This would seem to rule out the presence of a positive slope although we caution that the fit is heavily influenced by the five points at high hardness and further data would be useful to confirm this.

Whilst the hardness-variability trend is consistent with the expectations of our model, we proceed to investigate individual source behaviour to further test whether the observations match (or can be broadly explained by) the detailed predictions discussed in Section 3.

5 ENERGY-DEPENDENT VARIABILITY

A clear prediction of our spectral-variability model is that at small θ, variability should be present mostly on long time-scales (see Section 2) and should be a complex combination of the variability directly from the disc (assuming that the disc is not thin; Churazov et al. 2001), Compton down-scattered emission from the scattering surface of the wind and changing f_col. At intermediate θ, we should see a slightly different energy dependence of the variability (as obscuration events are expected to dominate where m₀ is not too large) with the same spectrum as the inner disc emission (which arrives to the observer down-scattered).

The energy dependence of the variability in NGC 5408 X-1 has already been shown in Middleton et al. (2011a) with the fractional rms increasing with energy. However, it is clear that those spectrally hard sources where the variability is predicted to be created by scattering into the line of sight (i.e. those at small θ) all have low rms and so we cannot investigate the shape of the rms spectrum. We can improve on this situation by extracting the cross spectrum (Nowak et al. 1999) and selecting the linearly correlated/coherent variability relative to a reference band with high signal-to-noise rms. Normalizing by the excess variance in the reference band and plotting this against energy gives the covariance spectrum (Wilkinson & Uttley 2009) where the removal of uncorrelated Poisson noise significantly reduces the sizes of the errors on each data point (of variance) in each energy bin. We plot the covariance (extracted in the time domain) relative to the 1.5–3 keV band, for the best-constrained single observation of each source (with the exception of IC 342 X-1, NGC 1313 X-2 and NGC 5204 X-1 due to insufficient data quality) in Figs 9, 10 and 11 over two time-scales: long (0.9–3 mHz) and short (3–200 mHz). The short time-scale variability corresponds to the fractional rms values reported in Tables 2 and 3 and we expect to include the contributions from both methods of generating variability discussed in Section 2. The covariance spectra interrogated on long time-scales allow for a comparison (where variability on both time-scales can be constrained) but also allows us to investigate the correlated spectral variability for sources with low rms on the shorter time-scales. We plot these unfolded through, and plotted alongside, the de-absorbed, best-fitting time-averaged model for each observation to allow ease of comparison.

Figure 9.

Upper panels: covariance spectra relative to the 1.5–3 keV band, over two time-scales: long (red: 0.9–3 mHz) and short (black: 3–200 mHz) for NGC 5408 X-1 (Obs.ID: 0500750101) and NGC 6946 X-1 (Obs.ID: 0691570101) plotted with (and unfolded though) their best-fitting, de-absorbed, time-averaged spectral model. Lower panels: ratio plots from fitting the covariance spectra with the time-averaged, de-absorbed model with the relative amount of DISKBB to NTHCOMP fixed (black), free (green) and with the NTHCOMP component only, with free model parameters (blue; note that kT_e is fixed to its best-fitting value in the time-averaged spectral fitting). When the normalizations are free, the DISKBB component is not required; clearly, the variability is therefore associated with the hard component. The ratio plots also show that, whilst the variability spectrum is a good match to the high-energy component above ≈1 keV, at softer energies the model overpredicts the data as expected from our model assumptions (Section 3).

Open in new tabDownload slide

Figure 10.

As for Fig. 9, covariance spectra for Ho II X-1 (Obs.ID: 0561580401), which can only be extracted on long time-scales, and NGC 1313 X-1 (Obs.ID: 0405090101).

Open in new tabDownload slide

Figure 11.

As for Fig. 9, covariance spectrum of Ho IX X-1 (Obs.ID: 0657801801).

Open in new tabDownload slide

In the case of the softer sources (NGC 5408 X-1, NGC 6946 X-1 and to a lesser extent, Ho II X-1: Figs 9 and 10), we see a shape consistent with the variability originating in the hard component only, consistent with a model where the variability on short and long time-scales at such inclinations is dominated by obscuration of the high energy emission (Method 1). To demonstrate this more clearly, we fit the short-time-scale covariance data (long time-scales in the case of Ho II X-1) with the de-absorbed, time-averaged model where the component normalizations relative to one another are initially fixed (so that there is the same proportion of DISKBB to NTHCOMP as in the time-averaged, best-fitting model); this leads to large residuals (plotted as a ratio in the lower panels of Figs 9 and 10). We then proceed to free the normalizations, finding that the residuals improve and the DISKBB normalization tends to zero. It is apparent that, whilst the variability is generally well matched by our choice of hard component model above 1 keV, the spectral modelling may deviate somewhat at the very softest energies (most readily seen for NGC 6946 X-1). This is not unexpected; we forced the lowest temperature of the NTHCOMP component to be that of the wind via energy balance, however, this assumes the most extreme scenario where the down-scattering occurs by the coolest wind material at the largest radii. Instead, for sources at more moderate or low inclinations, the temperature of the wind plasma intercepting the hard photons is expected to be somewhat higher, leading to the low-energy peak (of the high-energy component) moving to higher energies (or in other words kT_in > kT_d in our model). As we are using thermal components in our spectral fitting rather than simple power laws, the effect on the spectral hardness will only be small (as the spectrum is down-turning to meet the peak of the seed photon distribution: Figs 5, 6 and 7) and so our spectral-variability trend (Fig. 8) remains robust. As an interesting aside, as the hard component rolls over at lower energies/faster than we have accounted for, it seems likely that the soft component could be of a different shape to that determined from the time-averaged spectra (which we will investigate elsewhere). In order to obtain the best description of the covariance spectra, we allow the parameters of the NTHCOMP component to be free with the exception of the peak temperature, kT_e, set to the best-fitting value in the time-averaged spectrum (Tables 2 and 3) to keep the rollover within the bandpass. We also restrict the lower limit of kT_d, to that of its time-averaged, best-fitting value as it is clear that the true temperature likely resides at or above this value. Fitting with the NTHCOMP component only, provides the greatest improvement as shown in the lower panels of Figs 9 and 10. The best-fitting parameters and their errors are given in Table 4, although we caution that these values should only be seen a representative of the spectral shape rather than a physical description. It is interesting to note that, in the case of NGC 6946 X-1, even this component alone does not appear to account for the shape of the covariance at the softest energies.

The covariance spectra of Ho IX X-1 and NGC 1313 X-1 (Figs 10 and 11) highlight the variable component for the hardest ULXs which we have not been able to previously study. Once again we show the residuals (as a ratio) from fitting the covariance (short time-scales in the case of NGC 1313 X-1 and long time-scales for Ho IX X-1) with the de-absorbed, time-averaged spectral model with a fixed proportion of normalizations, free normalizations and with the NTHCOMP component alone (with free model parameters, with the exception of kT_e). In the case of the highest data quality, i.e. NGC 1313 X-1, we see that the best description of the covariance spectrum is once again the hard component on its own (i.e. the DISKBB normalization tends to zero) with free model parameters (see Table 4). Based on our model, we may have expected the variability spectrum to have a contribution due to down-scattering below that of the hard component with the addition of variability in the soft component due to changing f_col. Although the covariance spectrum may appear slightly flatter than the de-absorbed, time-averaged model, consistent with expectation (appearing as an excess at ∼1–2 keV above the NTHCOMP component in the residuals for Ho IX X-1 in Fig. 11), the data quality is not yet sufficient for any unambiguous claim of such broadening.

Of particular interest for our model is the source NGC 55 ULX-1 as, due to its extremely soft spectrum, this may sit on the positive branch of the hardness-variability plot (Figs 4 and 8). The source has been studied by Stobbart, Roberts & Warwick (2004) who suggested that the ‘dips’ in the X-ray light curve could be analogous to those seen in the class of dipping neutron star binaries (i.e. obscuration by optically thick material in the outer disc), implying a large (yet small enough that the source is X-ray bright) inclination to the source. Notably, the dips get stronger with increasing energy and so also fit into our model as a ULX seen at moderate-to-high inclinations. Although the short time-scale variability is unconstrained in the first observation, on similarly long time-scales (0.03–3 and 0.02–2 mHz, respectively) both observations have sufficient variability for the covariance spectrum to be extracted and compared, and are shown in Fig. 12. The second observation, which is significantly softer than the first observation, is dimmer as we might expect if the wind has become stronger or the system has precessed. The covariance spectrum of the first observation is dominated by a hard component (as can be seen from the spectral residuals), with a peak temperature consistent with that of the time-averaged component, and appearing similar to those of other ULXs we have identified as moderately inclined (e.g. NGC 5408 X-1, NGC 6946 X-1). The second observation shows that the variability spectrum is considerably flatter and the peak has shifted to lower energies between observations. By fitting a simple blackbody (BBODY in xspec) to the covariance spectra, we find that the peak has indeed shifted from kT_obs = 0.71|$^{+0.06}_{-0.05}$| to 0.28|$^{+0.14}_{-0.07}$| keV and is no longer a reasonable match to the temperature of the hard component within errors (Table 3). Under the reasonable assumption that the variability component is of the same physical origin, such a change is not expected in any model where the variable emission arrives to the observer directly (e.g. a disc–corona model as seen in BHBs). The nature of the evolution implies that the variable component of the emission has reached us after being down-scattered by cooler/more optically thick material than intercepts the hard emission. How this occurs is presently unclear but further observations of NGC 55 ULX-1 (and similar sources) may help place constraints on the required geometry.

Figure 12.

Covariance spectra of NGC 55 ULX-1 on long time-scales (Obs.1: 0028740201, black: 0.03–3 mHz; Obs.2: 0655050101, purple: 0.02–2 mHz) which shows the decrease in peak temperature as the spectrum softens and dims, implying that we do not see the variable component of the emission directly. The lower panels show the residuals (as per the previous figures) for each observation in turn.

Open in new tabDownload slide

6 DISCUSSION

One of the leading models to explain ULX luminosities above ∼3 × 10³⁹ erg s⁻¹ is the formation of a powerful wind at the spherization radius as a result of a high mass transfer rate from the donor star (as seen in SS433, see Fabrika 2004, for a review). Radiative hydrodynamic instabilities are expected to be present leading to a clumpy structure (Takeuchi et al. 2013, 2014). When combined with inflowing propagating mass accretion rate fluctuations – now seen as the likely origin for the broad-band variability characteristics of accreting sources (Lyubarskii 1997; Uttley et al. 2005; Ingram & Done 2012; Scaringi et al. 2012) – this leads to the prediction of variability emerging directly (should the inner disc not be geometrically thin; Churazov et al. 2001) and/or via scattering events.

In Section 2, we presented two methods of imprinting variability which vary as a function of inclination angle and mass accretion rate. Combining these with the spectral predictions from P07, leads to a series of observational spectral/timing predictions, as outlined in Section 3. These can be complicated by disc precession but the overall trend we should expect across the population is one where the harder sources, seen at small inclinations are generally less variable than the soft sources seen at moderate inclinations. We also expect the variability to drop as the wind becomes increasingly dominant and the source becomes softer. The predicted, inflected path (shown in Fig. 4), differs from that expected should the variability originate in a high-energy Compton-scattering plasma, as seen in BHBs, where, crucially, the variability is greatest for the hardest source states (Belloni 2010; Muñoz-Darias et al. 2011).

We combine spectral hardness, obtained from the best-fitting models to the spectra of a sample of ULXs, with rms values obtained from the PDS and plot these in Fig. 8. There is considerable scatter in the plot – expected in light of different evolutionary paths and a distribution of masses. However, we find a significantly negative slope (when ignoring the softest data point – which we have assumed sits on the positive slope of the predicted trend) which is consistent with expectation and would disagree with the expectation of IMBH accretion based on scaling BHBs to IMBH masses (Belloni 2010; Muñoz-Darias et al. 2011).

We also consider whether the variability could be caused by a variable, hot inner disc and stable wind, i.e. no contribution via Method 1. We note that, in the case of BHBs a thin disc can appear variable but only if illuminated by a variable photon flux (Wilkinson & Uttley 2009), otherwise (under the assumption that variability becomes damped by the density of the disc) the disc would need to not be geometrically thin (Churazov et al. 2001), pointing to a situation unlike that expected for IMBHs. In Method 2, we consider that the propagating flux can emerge via scattering and/or from the inner disc directly. Assuming no contribution via Method 1 (i.e. obscuration by clumps), the implication would be that the soft sources are at a lower mass accretion rate so that the dampening due to mass-loss is less severe (in effect more of the PDS is in our bandpass: see Fig. 3 and Section 2.2.3). Whilst we cannot rule this out, we still invoke a wind in this situation and so it seems highly likely that radiative hydrodynamic instabilities will occur leading to clumps and increased variability via Method 1 for these sources.

New models incorporating inhomogeneous accretion disc emission (e.g. Dexter & Quataert 2012) based on the photon bubble instability (Gammie 1998; Begelman 2001) may also describe the X-ray emission (Miller et al. 2014) and may offer a viable alternative to distortion by radiatively driven winds. However, a key test of any such model will be their ability to explain and incorporate the changing variability properties which, as we have discussed, can naturally be explained by powerful radiatively driven winds.

In terms of individual source behaviour, which may differ from the overall trend due to individual characteristics e.g. precession, we can consider whether the spectral evolution with variability is consistent with our model and we do so on a source-by-source basis.

NGC 5408 X-1: the XMM–Newton observations of this source cover a luminosity range of only 7–9 × 10³⁹ erg s⁻¹ and show large and well-constrained fractional variability in each (Table 2). The X-ray spectrum changes very little, with model parameters remaining consistent throughout and similar fractional increases in soft and hard flux (from faintest to brightest: Fig. 5, Table 2). Although it would appear that the fractional variability drops with increasing flux in the hard component, possibly in keeping with our expectations of the wind tending towards homogeneity (see also Caballero-García, Belloni & Wolter 2013), this is not well constrained. However, the spectral and variability characteristics are fully consistent with a source that has been viewed at moderate inclinations such that the wind is in our line of sight with large amounts of variability (dominated by Method 1).

NGC 6946 X-1: there are clear spectral similarities between NGC 6946 X-1 and NGC 5408 X-1 (see Fig. 5) with the source showing only small total variations in flux (Table 2). We also detect large amounts of variability; indeed NGC 6946 X-1 appears to demonstrate even larger amounts of fractional variability than NGC 5408 X-1. It is therefore a reasonable assertion that we have viewed the ULX at similar inclinations. In addition to our analysis, NGC 6946 X-1 has been reported to have extreme (ultraluminous) UV emission (∼4 × 10³⁹ erg s⁻¹; Kaaret et al. 2010). Given the brightness this could be associated with down-scattering in the wind and/or emission from the outer photosphere, and would imply moderate to high inclinations (depending on accretion rate – P07 – and any precession at the time).

NGC 5204 X-1: the source appears to become considerably softer with increasing luminosity (Fig. 5) and the best-fitting model to the brightest observation would imply that the hard emission becomes much cooler (Table 2), although this is not well constrained. The evolution therefore appears consistent with a line of sight changing such that the observer views increasingly into the wind. As a result, the hard emission is progressively beamed out of the line of sight with the remaining emission in the direction of the observer being increasingly Compton down-scattered. In our model, we associate this behaviour with an inclination to the wind changing due to either precession or an increase in |$\dot{m}_{\rm 0}$|⁠.

NGC 1313 X-1: the evolution of this source (and X-2) has been the subject of many papers (e.g. Pintore & Zampieri 2012) and appears remarkably similar to that of NGC 5204 X-1 (Fig. 5) with a larger fractional increase in flux at soft energies than hard, although it traverses a larger range in luminosity (6–12 × 10³⁹ erg s⁻¹). Once again, we see a spectrum (Obs.ID 0205230401) that appears to have a considerably cooler and weaker hard component (Table 2). Although this is not the brightest observation, should we be viewing the source so that more of the wind enters the line of sight, then a large portion of the intrinsic flux may be beamed out of the line of sight and this may instead be at higher |$\dot{m}_{\rm 0}$| or precessed.

Unlike NGC 5204 X-1, we can constrain the presence of variability in two observations (Table 2). Both observations are at the lowest flux levels (excluding Obs.ID 0205230401), consistent with our model where the variability drops at higher luminosities (i.e. higher |$\dot{m}_{\rm 0}$|⁠) in either method.

Based on our model, we predict that our view of the ULX is evolving from low inclinations, where the wind does not substantially enter our line of sight so that the source appears hard with variability via Method 2, to one where the wind now enters our line of sight, the spectrum gets softer and variability becomes suppressed by increasing |$\dot{m}_{\rm 0}$|⁠.

NGC 1313 X-2: at the lowest fluxes, the source appears very similar to NGC 1313 X-1 and NGC 5204 X-1 (Figs 5 and 6, Table 3) implying once again that we are viewing at moderate inclinations to the wind. However, unlike the former sources, the spectrum becomes increasingly hard with luminosity (Fig. 6), suggesting that the inner regions are not significantly obscured by the wind and are instead geometrically beamed with increasing mass accretion rate, i.e. we are looking at small inclinations.

At lower fluxes, we detect variability (Table 3) which, given the hard spectrum, we attribute to Method 2. We note that, once again the variability does not appear in observations at higher luminosities even when the data quality is extremely high, e.g. Obs.ID 0405090101.

We also note the more recent observations presented by Bachetti et al. (2013) with two XMM–Newton observations of the source at low–moderate luminosities. Neither show constrained variability above 3 mHz (our band of interest) but the brighter observation may show variability on longer time-scales which could feasibly be explained by the source transiting along the softer branch of Fig. 4 (or a combination of precession with decreasing |$\dot{m}_{\rm 0}$|⁠).

HoII X-1: the spectra generally appear very similar to those of NGC 5408 X-1 and NGC 6946 X-1 (Figs 5 and 6) although span a much greater range in luminosity (3–10 × 10³⁹ erg s⁻¹); we therefore also predict that we have observed the source at moderate inclinations.

We note the substantial difference in levels of variability compared to NGC 5408 X-1 and NGC 6946 X-1 and, whilst we detect low fractional rms in a single observation of HoII X-1 (Obs.ID 0112520601), this is poorly constrained. Certainly, we do not detect variability in the highest quality data set (Obs.ID 0200470101) at the highest observed luminosity, which might imply that the wind has already tended towards homogeneity, i.e. a higher |$\dot{m}_{\rm 0}$|⁠. It is important to note that, in the case of HoII X-1, the detection of radio lobes (Cseh et al. 2014) suggest that the source is not likely to have been viewed at the smallest inclinations (the limit of not being strongly Doppler boosted gives a lower limit of 10°).

HoIX X-1: the evolution of HoIX X-1 has been studied by Vierdayanti et al. (2010) and appears similar to NGC 1313 X-2 (Fig. 6, Table 3). As the source is notable as one of the spectrally hardest ULXs, and appears to get harder with increasing luminosity, we predict that we are viewing at small inclinations to the wind so that it remains out of the line of sight throughout the existing observations. We attribute the general lack of variability (even given the high data quality) to dampening due to high |$\dot{m}_{\rm 0}$| (Method 2).

It is worth noting that Walton et al. (2014) observe spectral evolution of Ho IX X-1 across a broader energy bandpass using NuSTAR, in concert with XMM–Newton and Suzaku. The evolution may be consistent with a model of increasing hardness with luminosity (in the XMM–Newton bandpass) with the brightest epoch showing a cooler peak of the inner disc component, possibly due to increased Compton down-scattering in the (now more narrow) wind cone. However, given present model degeneracies, the exact nature of the evolution (in terms of spectral components) is not yet unambiguous.

IC 342 X-1: the spectra of IC 342 X-1 and their evolution with luminosity appear similar to those of Ho IX X-1 (Fig. 6), implying inclinations that do not intercept the wind. In addition, we find that only the dimmest observation with adequate data quality (Obs.ID: 0206890201) shows constrained variability (Table 3), which we associate with Method 2, implying that, as with HoIX X-1, the increase in brightness due to an increase in |$\dot{m}_{\rm 0}$| damps the variability.

NGC 55 ULX-1: the spectra of NGC 55 ULX-1 are the softest of our sample implying moderate to high inclinations to the wind. As variability is present, albeit at low levels, this implies a |$\dot{m}_{\rm 0}$| low enough so that variability has not yet been fully damped.

Whilst the first-order spectral-variability behaviour in individual sources as well as across the sample can be explained by our model, we can gain further insights by studying the nature of the energy dependence of the variability. Whilst past studies have used the fractional rms-spectrum (e.g. Middleton et al. 2011a), here we extracted the covariance spectra (Wilkinson & Uttley 2009) which have considerably smaller errors and allow a view of the correlated variability (relative to a reference band). In all cases, where statistics allow, we should expect the hard component to be variable due to obscurations (and potentially intrinsic disc emission should it not be thin) at larger inclinations, and an intrinsically variable hot inner disc (again, should it not be thin) and scatterings into the line of sight of the observer at smaller inclinations. The major difference we predict is that there could be some contribution to the shape of the covariance spectrum from down-scattering and changing f_col at small inclinations (although changes in f_col may also add variability to sources at larger inclinations, we expect the variability to be dominated by obscurations). In all cases, we do indeed see that the hard component rather than the soft is varying with the disc normalization tending to zero. This is potentially important as the covariance spectra of GX 339-4 and SWIFT J1753.5-0127 (Wilkinson & Uttley 2009) when in the hard state (analogous to many ULXs should they contain more massive IMBHs), shows the disc to be variable with a similar proportion of disc to power law as in the time-averaged spectrum (with the proportion of disc increasing to longer time-scales). This is clearly not the case in the sources here (see the black residuals in Figs 9–11) as the disc contributes negligible variability, and implies that we are not seeing the expected behaviour for IMBHs (>100 s of M_⊙). Unfortunately, the data quality is not presently high enough to confirm possible hints of down-scattering in the covariance spectra of those sources we have identified as viewed at small inclinations to the wind. However, based on the softening spectra, we predict that NGC 55 ULX-1 has been observed with the wind increasingly entering our line of sight by precession or increasing |$\dot{m}_{\rm 0}$|⁠. The covariance spectra also show constrained evolution and imply that the variable component in the second observation has intercepted cooler/optically thicker material to that seen by the rest of the hard emission, inconsistent with models where the hard component is seen directly (e.g. disc–corona models of BHBs).

The covariance spectra have proven valuable not only in allowing us to test consistency with predictions but also for highlighting the issues associated with spectral fitting which can often prove degenerate. In the case of our ULX sample, we note that the time-averaged, best-fitting, high-energy spectral component (NTHCOMP) overestimates the true component at soft energies. This results from tying the two spectral components together to ensure physical energy balance (Section 3) but does not account for the changing temperature profile of the wind plasma. As a result, the soft component is not modelled accurately (by DISKBB in our fitting) and could conceivably allow for contribution from advection (P07). Whilst this does not cause problems for the inferred spectral evolution with variability (nor the apparent lack of a disc component in the covariance spectrum), it highlights the need for caution when obtaining characteristic temperatures and spectral parameters from fitting to broad continua. However, the variability spectrum should provide a means to model the soft component more accurately in future, thereby allowing us to probe its nature more thoroughly. We note that, whilst the covariance spectra point towards issues in characterizing the soft component, the rollover at high energies should be unaffected. Whilst the nature of this component is still unclear, if we identify it with the hot inner disc (see e.g. Walton et al. 2014), then the change in position of the rollover across the population should result from the changing amount of wind down-scattering the peak to lower energies. Whilst this seems broadly consistent with those sources where we believe the wind to have entered our line of sight (e.g. NGC 5204 X-1, NGC 1313 X-1, NGC 55 ULX-1: see Figs 5 and 7), quantifying the change in temperature is a relatively complex issue but one we will address in a future work.

7 CONCLUSIONS

In this paper, we have presented a model that results from considering the structure of the supercritical accretion flow (Shakura & Sunyaev 1973; P07; King 2009), clumps formed through radiative-hydrodynamic instabilities (Takeuchi et al. 2014) and propagating mass accretion rate fluctuations (Lyubarskii 1997; Arevalo & Uttley 2006; Ingram & Done 2012). This invokes two methods of generating variability; via clumps themselves on short time-scales and via longer time-scale trends imprinted via propagating fluctuations. Combining the expected changes in the spectrum with the variability leads to a set of spectral-timing predictions which depend on inclination and accretion rate (Section 3 and Fig. 4). In order to test our predictions, we have analysed a sample of nine ULXs which are bright in flux (thereby providing the highest quality available statistics) and avoid those ULXs which are likely to be BHBs with less extreme mass transfer rates. Although there is considerable scatter, we find that the trend of spectral hardness with variability power – taking the sample as a whole – is consistent with the predicted evolution (Figs 4 and 8) and a significantly negative slope which is inconsistent with accretion on to IMBHs (if we assume these would follow the ubiquitous path in hardness-variability of BHBs; Belloni 2010; Muñoz-Darias et al. 2011). Such a finding is therefore consistent with the spectral properties being associated with supercritical accretion (as suggested by the mass measurement of a bright ULX; Motch et al. 2014).

We present the covariance spectra of those sources where data allows for the first time. We note that these are broadly consistent with the expectation of obscuration/scattering/intrinsic variability of the inner disc (i.e. only the hard component varying), although the quality is not sufficient for detecting the subtle deviations we predict should be present. Whilst the shapes of the covariance spectra are also consistent with a variable high-energy component/corona as seen in BHBs (although we strongly rule out a power-law shape for this component in the time-averaged spectra), we note that there appears to be very little/no contribution in any source by the soft component, inconsistent with what is seen in BHBs in the low-hard state (Wilkinson & Uttley 2009), where disc reprocessing leads to a variable component. Whilst this appears to indicate that the emission is not analogous to that seen in the low-state of BHBs, we also reiterate that, more importantly, the global trend in hardness-variability does not appear to match predictions for such a model-invoking analogous accretion on to IMBHs.

Whilst the covariance spectra in general do not show evidence that unambiguously links the sources to the model we have presented, in the specific case of NGC 55 ULX-1 we observe evidence for down-scattering of the radiation, which argues against the variable hard component being seen directly. Finally, the covariance spectra highlight the present degeneracies inherent in modelling ULX spectra; notably, the contribution of the soft component was generally underestimated, although future studies making use of such a technique should be able to better determine the shape of the soft component and more accurately constrain its evolution with luminosity.

Although other models may be able to explain some of the behaviours of ULXs, the one we present here benefits from also being able to explain the lack of narrow atomic features in the X-ray spectra of the hardest ULXs (e.g. Walton et al. 2012, 2013b) – explained by the inner, scattering surface of the wind cone being highly ionized by the strongly beamed flux – and possible presence of strong broad absorption features in the softest ULXs (Middleton et al. 2014a). We briefly note that variability caused by scattering from clumps in the manner described in our model will lead to variable illumination on to the optically thin expanded wind at larger radial distances (directly and indirectly via reprocessing). Although we recognize that there are multiple explanations, should the ionization state of this material change as a result, then we should expect changes in the (possibly broad) atomic features (Middleton et al. 2014a) which could then explain the ‘soft lag’ discovered in the cross spectrum of NGC 5408 X-1 (Heil & Vaughan 2010, De Marco et al. 2013). Another possible explanation is a change in the soft continuum in response to the variable hard emission (with the lag corresponding to the extra light travel time); the covariance spectra of NGC 5408 X-1 would imply that the soft component as a whole does not contribute significantly to the variability but we cannot rule out a ring of material with a small area being illuminated at the inner edge of the wind. Such behaviour may place independent constraints on the geometry and we will investigate this in a future, dedicated paper.

There is clearly still a great deal of work left to be done in exploring the more complex predictions of such a model and any outliers that could be associated with unusual and exciting objects (e.g. HLX-1: Farrell et al. 2009, M82 X-1: Pasham et al. 2014, M82 X-2: Bachetti et al. 2014). Key to making progress will be to utilize the new powerful spectral-timing techniques, investigate the nature of the broad-band spectra (as NuSTAR is demonstrating to great effect) and obtain high-quality spectra from multiple epochs, which next-generation missions such as ATHENA will regularly produce.

The authors thank the anonymous referee for helpful suggestions, and Tom Maccarone, Phil Uttley, Adam Ingram and Andrew King for useful discussion. MJM acknowledges support via a Marie Curie FP7 Postdoctoral scholarship. DJW is supported by ORAU under the NASA Postdoctoral Programme. TPR was funded as part of the STFC consolidated grant ST/K000861/1. This work is based on observations obtained with XMM–Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA.

REFERENCES