Research ArticleASTROPHYSICS

Accretion-induced variability links young stellar objects, white dwarfs, and black holes

See allHide authors and affiliations

Science Advances  09 Oct 2015:
Vol. 1, no. 9, e1500686
DOI: 10.1126/sciadv.1500686


The central engines of disc-accreting stellar-mass black holes appear to be scaled down versions of the supermassive black holes that power active galactic nuclei. However, if the physics of accretion is universal, it should also be possible to extend this scaling to other types of accreting systems, irrespective of accretor mass, size, or type. We examine new observations, obtained with Kepler/K2 and ULTRACAM, regarding accreting white dwarfs and young stellar objects. Every object in the sample displays the same linear correlation between the brightness of the source and its amplitude of variability (rms-flux relation) and obeys the same quantitative scaling relation as stellar-mass black holes and active galactic nuclei. We also show that the most important parameter in this scaling relation is the physical size of the accreting object. This establishes the universality of accretion physics from proto-stars still in the star-forming process to the supermassive black holes at the centers of galaxies.

  • accretion
  • black holes
  • white dwarfs
  • young-stellar objects
  • timing
  • variability
  • unification


Accretion is the process by which most objects in the Universe grow in mass: from young stellar objects (YSOs) still in the star-forming process, through accreting white dwarfs (WDs) and neutron stars (NSs), to stellar-mass black holes (BHs) and supermassive BHs at the center of galaxies [active galactic nuclei (AGN)]. Because the accreting material almost always carries excess angular momentum, the process usually leads to the formation of a rotating disc that slowly transports the material inward. However, we still do not fully understand the physical processes taking place in accretion discs. In particular, it is not clear whether these processes operate in the same way across different types of astrophysical systems. Nevertheless, several phenomenological similarities between accreting systems have been discovered over the past decade (112). This suggests that the physics of accretion may indeed be universal, irrespective of accretor mass, size, or type.

One powerful probe of the accretion process is provided by the brightness variability that is generated in the disc itself. These brightness variations produce observable light curves that are phenomenologically similar across AGN, stellar-mass BHs, and accreting WDs (1, 7, 8). For example, all of these systems—as well as accreting NSs—exhibit aperiodic red-noise variability that produces a linear correlation between the root mean square (rms) variability amplitude and the mean flux across a wide range of time scales (913). The variability properties are commonly analyzed in frequency space through the power spectral density (PSD). The power density, P(ν), gives the distribution of power (rms-squared variability amplitude) as a function of frequency, ν (or, equivalently, time scale t = 1/ν). In systems containing high mass flow rate accretion discs—such as luminous AGN, stellar-mass BHs (in luminous states), and accreting WDs of the nova-like class—the low-frequency (long time scale) PSDs are usually well described by a simple power law, P(ν) ∝ ν−α, with α ≈ 1. However, at frequencies above a characteristic break frequency, νb, the PSDs “bend” to a steeper slope (α ≥ 2) (1, 7, 8). The physical process that determines this characteristic frequency is not clear, but the variability probably originates from the inner accretion disc surrounding the accreting object (1416). In AGN and stellar-mass BHs, this region primarily emits in x-rays; in accreting WDs and YSO, it mainly emits in the optical/ultraviolet (UV).


Some types of YSOs (pre–main sequence stars accreting from a disc of gas and dust) also show strong signatures of variability in the optical and infrared due to accretion (17). To determine whether the variability properties of noncompact accretors are the same as those of compact ones, we analyzed the Kepler data of six YSOs (18) that were observed during campaign 2 of the repurposed Kepler/K2 mission and that exhibited aperiodic variability. Figure 1 shows the obtained light curve for one of these objects, the T Tauri star V866 Sco, whereas Fig. 2 shows the resulting PSD and rms-flux relation. The rms-flux relation is clearly linear, and the PSD has a broken power-law shape with an obvious bend at νb ≈ 1.56 × 10−6 Hz (tb = 1/νb ≈ 7.41 days). The light curves of the other five YSOs in our sample also produce linear rms-flux relations on a wide range of time scales as well as broken power-law PSDs. This strongly suggests that these variability properties are universal among disc-accreting astrophysical systems.

Fig. 1 Kepler/K2 light curve of V866 Sco.

Light curve of the YSO object V866 Sco (EPIC205249328) obtained with the Kepler/K2 mission during campaign 2 in long cadence mode (29.4 min). The light curve was extracted from the target pixel files provided by MAST by manually defining large target and background masks.

Fig. 2 Power spectrum and rms-flux relation of the accreting YSO V866 Sco.

We show the obtained power spectrum of the accreting YSO star V866 Sco (AS 205A or EPIC205249328). The rms-normalized (42) power spectrum was estimated using Kepler data from campaign 2, comprising over 78 days of continuous observations with 29.4-min cadence. We mark the measured characteristic bend frequency with the vertical gray solid line, and 1σ limits with the gray-shaded region (19). The inset shows the linear rms-flux relation from the light curve in units of 105 electrons/s (black data points), computed on 2.5-hour time scales. The gray line is a linear fit to the data. We are able to measure characteristic bend frequencies for the other five YSOs in the Kepler sample (which include EPIC204181799, EPIC204395393, EPIC204830786, EPIC204908189, and EPIC205110000), as well as linear rms-flux relations, on a wide range of time scales.

If this is correct, the characteristic PSD break frequency, νb, may be expected to depend quantitatively—and exclusively—on the key physical parameters describing such systems. For black holes only, it has already been shown that νb can be predicted quantitatively across eight orders of magnitude in BH mass (1, 3, 4). More specifically, the BH scaling relation can be viewed as a plane in the parameter space given by the (logarithms of) BH mass, accretion rate, and break frequency. However, because the existing relation is derived solely from observations of accreting BHs, it is not clear whether and how it applies to other types of accreting systems. The same limitation has also prevented the unique identification of accretor mass as the underlying key scaling parameter, because, for black holes accreting at high rates, the characteristic inner disc radius and accretor mass linearly scale as RISCO ≈ 3GMBH/c2 (where G, c, MBH, and RISCO are the gravitational constant, speed of light, BH mass, and the innermost stable circular orbit of a BH with spin parameter a = 0.8, respectively; RISCO = 6GMBH/c2 and GMBH/c2 for nonrotating a = 0 and maximally spinning a = 1 BHs).

Motivated by the form of the scaling relation for BHs, which also allows degeneracy between mass and radius in BH samples, we make the ansatz that the break frequency should scale asEmbedded Image(1)

if accretion-induced variability is truly universal (where R, M, and Embedded Image are the characteristic radius, accretor mass, and accretion rate, respectively). To test for the existence of such a universal plane of accretion activity, we expand the existing BH sample by adding a set of eight accreting WDs. All of these belong to the nova-like class, in which the disc extends close to the WD surface. These systems are the key to breaking the BH mass-radius degeneracy because the radius of a WD is three orders of magnitude larger than RISCO for a comparably massive BH. Accreting NSs also display aperiodic accretion-induced variability, but we do not include them here because their PSDs display multiple breaks and quasi-periodicities (7, 13, 19), which make it difficult to uniquely identify the appropriate break frequency for them.

For our BH sample, we adopt accretor masses and break frequencies from the literature (1, 2023) and set the characteristic radius to RISCO with dimensionless spin parameter a = 0.8 (19). BH accretion rates are estimated from measured bolometric luminosities, Lbol, through the relation Embedded Image = Lbolc2, with η = 0.1 (24). Our WD sample includes two systems observed with the Kepler satellite (MV Lyrae and KIC 8751494), whose variability properties have already been studied previously (810). We add to these an additional six accreting WDs (BZ Cam, CM Del, KR Aur, RW Tri, UU Aqr, and V345 Pav) that have been observed with ULTRACAM (19, 25, 26). We use these data sets to estimate break frequencies with errors for all of our accreting WD systems (19). Where possible, we adopt published values for the WD masses (27) and use the theoretical WD mass-radius relationship to obtain the corresponding radii (28). If the WD mass is unavailable, we adopt a value of 0.75 solar masses (29). The mass accretion rate is taken to be 10−8 solar masses per year for all accreting WDs, which is a typical value for these systems (29). Our results are not sensitive to the exact parameter values we have adopted, and we also allow for conservative errors of 0.4 dex on Embedded Image and 0.2 dex on M and R across the whole sample (19).

We then perform a grid search in parameter space to determine the best-fit values of A, B, C, and D in Eq. 1 (19). The resulting fit is good [Figs. 3 and 4; see also Materials and Methods (19)] and consistent with a simple generalization of the previously established BH variability plane (1). More specifically, in our new generalized scaling relation, νb is most sensitive to the characteristic radius of the accreting object, rather than the accretor mass.

Fig. 3 Confidence contours for the mass, radius, and mass accretion rate indices.

We show the confidence contours for the dependency of the characteristic bend frequency on the accretor mass, radius, and mass accretion rate indices. We have fitted the model, log νb = A log R + B log M + C log Embedded Image + D, where νb represents the measured characteristic bend frequencies, and R, M, and Embedded Image are the radii, masses, and mass accretion rates, respectively (all using cgs units), to determine the best-fit values for A, B, C, and D. The obtained fit is good (X2 = 38.41 for 37 df) and consistent with the previously obtained fit (1), resulting in A = −2.07 ± 0.11, B = 0.043 ± 0.17, C = 0.95 ± 0.22, and D = −3.07 ± 2.61. The dark gray lines show the previously obtained fit using a BH-only sample (1), with 1σ errors marked with light gray regions. Given that the previous fit was solely based on BHs, a degeneracy between mass and radius existed, such that it was essentially fitting for νbRaMb and the finding that a + b ≈ −2. The contours refer to the fit including accreting WDs, AGN, and soft-state Galactic BHs and allowing the radius to be an additional free parameter. We thus show that the dominant parameter in the scaling law is a characteristic radius in the inner disc, and not the mass of the accretor as previously thought.

Fig. 4 Edge-on projection of our sample on the accretion variability plane.

We show that the predicted characteristic bend frequencies, derived by inserting the observed masses, radii, and mass accretion rates into the best-fit relationship to the combined supermassive BHs, stellar-mass BHs, and accreting WDs, agree very well with the observed characteristic bend frequencies. If the predicted and observed bend frequencies are identical, then an object will exactly lie on the black line. We display supermassive BHs as filled green squares, stellar-mass BHs in their soft state as filled red squares, and accreting WDs as filled blue squares. We additionally show the position of the YSO V866 Sco with the filled magenta circle using the best estimate for its magnetospheric radius at 6.4 solar radii. The error on the predicted νb assumes strict upper and lower limits for the disc truncation radius, from 0.1 to 2 AU solar radii (30). We thus demonstrate that the variability plane of accreting systems extends from supermassive BHs all the way to YSOs.

To further test the universal nature of the scaling relation, we also show in Fig. 4 the one YSO in our sample—V866 Sco—for which published estimates for the mass, accretion rate, and inner disc radius are available (30). The inner discs of some YSOs are thought to be truncated by the magnetosphere of the accreting proto-star. In V866 Sco, the truncation radius has been estimated by fitting the spectral energy distribution (SED) of the system and found to be about 1 to 10 times larger than the radius of the proto-star itself (30). We therefore adopt log R* = 0.5 ± 0.5 (where R* is the stellar radius) as the characteristic radius for this system. Such truncation radius estimates are model-dependent and uncertain, but V866 Sco nevertheless lies close to the universal variability plane. In Fig. 4, the point corresponding to V866 Sco shifts to the right—and hence closer to the predicted location—with decreasing inner disc radius. This may suggest that the inner disc radii of YSOs lie at the lower end of the ranges inferred from the SED fits.


The leading theory for accretion-induced variability in compact objects is the “fluctuating accretion” model (14, 15). In this picture, slow accretion rate fluctuations generated in the outer disc propagate inward, where they combine with the progressively faster fluctuations generated at smaller radii. This results in a multiplicative coupling of the variability on all time scales and naturally produces both linear rms-flux relations and broken power-law PSDs (14, 15, 3133). Our discovery of the same variability characteristics in YSOs establishes the universality of these features across all types of accreting systems, compact or otherwise. This is qualitatively consistent with the fluctuating accretion model because the process it envisages is scale-free (in the sense that it does not rely on a particular system mass or size).

The universal plane of accretion-induced variability is a further significant step forward in understanding the physics of accretion because it quantitatively links all types of disc-accreting astrophysical systems. For example, the identification of a characteristic inner-disc radius as the dominant parameter—as opposed to accretor mass or accretion rate—has obvious implications for the physical interpretation of accretion-induced variability, in general, and the fluctuating accretion model, specifically. Our work thus provides a solid basis for the development of a universal physical model of astrophysical accretion on all scales.


Sample selection

We attempt to include all relevant sources in our sample. Our sample includes two stellar-mass BHs observed during their soft state (Cyg X-1 and GRS 1915+105) (3, 7, 22, 3436) with strong enough variability to measure break frequencies. We also include all unabsorbed, bright AGN with good XMM-Newton observations (having more than 40-ks exposures) (21). In addition, we include both nova-like accreting WDs observed during the Kepler mission and all six nova-like accreting WDs observed with ULTRACAM (see Table 1 for a full list) in the sample of Barros (25). All systems used are thought to be in an analogous state to the BH (high)soft states, where the accretion discs are thought to extend very close to the accreting objects. Therefore, our assumption that the inner-disc edge of both stellar-mass BHs and AGN extends close to RISCO is valid.

Table 1 Samples and parameters used in this work.

The table lists the objects used in this work. We include the adopted masses, radii, mass accretion rates, and break frequencies, νb. Where these values have been taken from the literature, we provide the relevant reference. Where these values have been estimated in this work, we refer to Materials and Methods (19). All systems, except for the YSO V866 Sco, have been used for the fit shown in Figs. 3 and 4.

View this table:

Throughout this work, we assume RISCO for BHs with dimensionless spin parameter a = 0.8. Estimates of black hole spins suggest that there is likely to be a range of spin parameter values, but different methods seem to yield different estimates for the spin parameter, at least in some cases (37, 38). Deviations from nonspinning RISCO estimates will introduce scatter in the scaling relation, but we account for this by introducing systematic uncertainties to our sample. The estimated spins of the two stellar-mass BHs in the sample are both nearly maximal, and the estimated values of the AGN spins that exist for the AGN in the sample are mostly about a = 0.8, but most of the AGN in the sample do not have well-established spin estimates (39). If we set a = 0 (no spin), we recover slightly different scaling parameters, although they remain consistent with the previously derived mass-scaling relation (1).

Stellar-mass BHs in their (low)hard state, as well as accreting NSs, have already been shown to seemingly follow the original variability plane (3). However, all hard-state systems seem to be systematically offset from the general soft-state relation. Thus, we only consider soft-state analogs in this work to obtain a refined scaling relation that includes the radius dependence. Soft-state NSs often have very low amplitudes of variability, and much of the variability is coming from the very low frequency noise component, which is often suggested to be due to fusion processes on the NS itself (40). We will address the variability of hard-state systems, including NSs (in all states), in a later work.

We inspected all YSOs observed during campaign 2 of the Kepler/K2 mission under the GO program GO2056. These targets have been selected as displaying clear evidence of accretion-induced emission lines (18). The sample contained 71 systems observed in long cadence mode. We extracted light curves for 17 systems from this sample that did nothave neighbors close by in the target pixel files provided by the Mikulski Archive for Space Telescopes (MAST). This allowed us to create light curves using relatively large target and background masks, which mitigated the effect of spacecraft jitter in the resulting light curves (41). We visually identified 6 of these 17 light curves as displaying clear aperiodic variability. The light curve obtained for V866 Sco is shown in Fig. 1.

PSD fitting

The PSDs used in this study (both from Kepler and ULTRACAM data) were all estimated from evenly sampled sections of data using standard methods (42). Specifically, we computed the rms-normalized periodogram from each continuous section of data, merged the periodograms where appropriate, and averaged the geometrically spaced frequency bins.

We fit each individual PSD through weighted least squares with a bending power-law of the formEmbedded Image(2)where P(ν) is variability power as a function of ν; α and β are the power-law indices for the low- and high-frequency components, respectively, bending at frequency νb; R is a normalization constant; and N is a constant noise component that takes into account high-frequency power not intrinsic to the sources. Confidence intervals were computed by “bootstrap” resampling, which was performed 10,000 times to obtain a distribution for the characteristic bend frequencies, νb. We set the measured logνb as the mean of the logged distribution and 1σ error from its SD.

Scaling relation fitting

The characteristic break frequencies are thought to originate from the inner parts of the accretion disc. These emit most of their radiation in x-rays for AGN and stellar-mass BHs, whereas they emit mostly optical/UV light for accreting WDs and YSOs. Thus, we adopt values of νb for AGN and stellar-mass BHs estimated from x-ray data, whereas we estimate νb for accreting WDs and YSOs from optical data. Although this choice is physically motivated, we note that, for accreting WDs and stellar-mass BHs, optical and x-ray PSDs have already been shown to yield consistent νb for the same object (13, 43). Furthermore, although no such comparison has yet been possible for AGN, the high-frequency power-law slope has been shown to be similar between different AGN at both x-ray and optical wavelengths (44, 45).

We fit the accretion variability plane with a function of the form taken from Eq. 1 using centimeter-gram-second (cgs) units. We use published values of νb (and errors) and M for AGN and stellar-mass BHs (1, 2023) and measure our own νb values with errors for accreting WDs and YSOs. In addition, we adopt bolometric luminosities, Lbol, from the literature and translate these to mass accretion rates Embedded Image = Lbolc2, setting η = 0.1 (24) for the BH systems. Where possible, we adopt literature values for the masses of accreting WDs (27) and set M = 0.75 solar masses otherwise. We also set Embedded Image =10−8 solar masses per year for all accreting WDs we set (29). The uncertainties in M, R, and Embedded Image are not typically laid out well in the literature, but a general understanding of the scatter in the estimates of these parameters does exist. For the BHs, the scatter in M and R is mostly due to scatter in the M-σ relation (46) and uncertainty in the spin distribution, respectively, yielding errors of about 0.2 dex on both parameters. For the WDs, the systems have been taken to have typical values, and the uncertainties represent the breadth of the distribution of these values. Although the spread in M found in accreting WDs is about 0.14 dex (47), we adopt a slightly more conservative value of 0.2 dex accounting for the fact that the population used to obtain the measured spread on M contains very few accreting WDs of the nova-like class, which we use exclusively. The uncertainties in Embedded Image for BHs are mostly driven by uncertainties in the source radiative efficiencies (spin-dependent) and bolometric corrections (48, 49), which we find to be around 0.4 dex. For WDs, the uncertainties in Embedded Image are driven by uncertainties in the distances to the sources (50), which require assuming that the nova-like accreting WDs all have similar luminosities, resulting in a 0.4 dex error on Embedded Image. We therefore simply adopt a constant 0.2 dex for the uncertainties in M and R and a constant 0.4 dex for the uncertainty in Embedded Image (because this is the hardest to estimate accurately) for the whole sample. We note that modifying the errors on the sample does not affect the obtained fitted parameters for the scaling relation, but rather change the size of the obtained contour levels (Fig. 3).

To determine the best-fit values for the coefficients A, B, C, and D, we performed a total least squares parameter grid search (using errors on all variables) and determined the minimum value of Χ2Embedded Image(3)where E is the log of the predicted frequency, given by the model (Eq. 1). Here, Embedded Image, σR, σM, and σEmbedded Image are the errors on logνb, logR, logM, and logEmbedded Image, respectively. At least to within our adopted errors, the fit is good (X2 = 38.41 for 37 df), and the coefficients of the logR, logM, and log Embedded Image terms are consistent with the previously obtained fit for the BH-only sample (1), resulting in A = −2.07 ± 0.11, B = 0.043 ± 0.17, C = 0.95 ± 0.22, and D = −3.07 ± 2.61. By setting the errors on M, R, and Embedded Image to larger or smaller values, we recover consistent fit parameters as with our original choice of 0.2 dex on M and R and 0.4 dex on Embedded Image, albeit with larger or smaller error contours. Thus, our analysis is robust against the precise choice of adopted errors. We also recall that our best fit assumes BH spins of a = 0.8 for the whole sample. Deviations from this assumption using either a = 0 or a = 1 would affect both RISCO estimates and the accretion efficiency η. The difference between these two extremes yields an about 0.8 dex difference on the BH predicted model frequencies (Fig. 4).

This is an open-access article distributed under the terms of the Creative Commons Attribution license, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


Acknowledgments: We wish to thank G. Ponti, R. Baptista, J. Green, and M. Middleton, as well as the anonymous referee, for useful discussions relating to this project. Funding: S.S. acknowledges funding from the Alexander von Humboldt Foundation. S.V., C.K., V.S.D., and T.R.M. acknowledge funding from the Science & Technology Facilities Council. E.K. and E.A. acknowledge funding from the Dutch research council under grant no. 2013/15390/EW. S.C.C.B. acknowledges support from Fundacao para a ciencia e tecnologia through the Investigador FCT Contract No. IF/01312/2014. Some of the data presented in this paper were obtained from MAST. STScI is operated by the Association of Universities for Research in Astronomy Inc., under NASA contract NAS5-26555. Support for MAST for non-HST (Hubble Space Telescope) data is provided by the NASA Office of Space Science through grant NNX09AF08G and by other grants and contracts. This paper includes data collected by the Kepler mission. Funding for the Kepler mission is provided by the NASA Science Mission directorate. Author contributions: S.S. initiated the project; collected, reduced, and analyzed the data; interpreted the results; and prepared the manuscript. C.K., E.K., T.J.M., and S.V. equally contributed to project design and interpretation of the results and commented on, discussed, and edited the manuscript. T.R.M. and V.S.D. also commented on and discussed the manuscript. E.A., T.R.M., V.S.D., and S.C.C.B. observed and reduced the ULTRACAM data. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All Kepler and K2 data can be found through MAST. The raw ULTRACAM data are available on request from V.S.D. ( or T.R.M. (
View Abstract

Stay Connected to Science Advances

Navigate This Article