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

Astrophysical accretion is a universal process in objects from proto-stars to supermassive black holes.


INTRODUCTION
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 (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12). 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 (9)(10)(11)(12)(13). The variability properties are commonly analyzed in frequency space through the power spectral density (PSD). The power density, P(n), gives the distribution of power (rms-squared variability amplitude) as a function of frequency, n (or, equivalently, time scale t = 1/n). 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 lowfrequency (long time scale) PSDs are usually well described by a simple power law, P(n) º n −a , with a ≈ 1. However, at frequencies above a characteristic break frequency, n b , the PSDs "bend" to a steeper slope (a ≥ 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 (14)(15)(16). 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).

RESULTS
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 n b ≈ 1.56 × 10 −6 Hz (t b = 1/n 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.
If this is correct, the characteristic PSD break frequency, n 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 n 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 R ISCO ≈ 3GM BH /c 2 (where G, c, M BH , and R ISCO 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; R ISCO = 6GM BH /c 2 and GM BH /c 2 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 as if accretion-induced variability is truly universal (where R, M, and Ṁ are the characteristic radius, accretor mass, and accretion rate, respec-tively). 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 R ISCO 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,(20)(21)(22)(23) and set the characteristic radius to R ISCO with dimensionless spin parameter a = 0.8 (19). BH accretion rates are estimated from measured bolometric luminosities, L bol , through the relation Ṁ = L bol /hc 2 , with h = 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 (8)(9)(10). 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 Ṁ 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, n b is most sensitive to the characteristic radius of the accreting object, rather than the accretor mass.
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 rightand 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.

DISCUSSION
The leading theory for accretion-induced variability in compact objects is the "fluctuating accretion" model (14,15). In this picture, slow  (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 1s limits with the gray-shaded region (19). The inset shows the linear rms-flux relation from the light curve in units of 10 5 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.
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,(31)(32)(33). 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,(34)(35)(36) 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 R ISCO is valid.
Throughout this work, we assume R ISCO 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 R ISCO 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 . 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 (X 2 = 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 1s 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 n b º R a M b 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. 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 n 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. 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 not have 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 form where P(n) is variability power as a function of n; a and b are the power-law indices for the low-and high-frequency components, respectively, bending at frequency n 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, n b . We set the measured logn b as the mean of the logged distribution and 1s 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 n b for AGN and stellar-mass BHs estimated from x-ray data, whereas we estimate n b for accreting WDs and YSOs from optical data. Al-though 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 n 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 n b (and errors) and M for AGN and stellar-mass BHs (1, 20-23) and measure our own n b values with errors for accreting WDs and YSOs. In addition, we adopt bolometric luminosities, L bol , from the literature and translate these to mass accretion rates Ṁ = L bol /hc 2 , setting h = 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 Ṁ =10 −8 solar masses per year for all accreting WDs we set (29). The uncertainties in M, R, and Ṁ 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-s 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 Ṁ 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 Ṁ are driven by uncertainties in the distances to the sources (50), which require assuming that the novalike accreting WDs all have similar luminosities, resulting in a 0.4 dex error on Ṁ. 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 Ṁ (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 Χ 2 where E is the log of the predicted frequency, given by the model (Eq. 1).
Here, s n b , s R , s M , and s Ṁ are the errors on logn b , logR, logM, and log Ṁ, respectively. At least to within our adopted errors, the fit is good (X 2 = 38.41 for 37 df), and the coefficients of the logR, logM, and log Ṁ 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 Ṁ 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 Ṁ, 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 R ISCO estimates and the accretion efficiency h. The difference between these two extremes yields an about 0.8 dex difference on the BH predicted model frequencies (Fig. 4).