Stochastic Variability in X-Ray Emission from the Black Hole Binary GRS 1915+105

The Astronomical Journal, 143:148 (10pp), 2012 June  C 2012. doi:10.1088/0004-6256/143/6/148 The American Astronomical Society. All rights reserved. Printed in the U.S.A. STOCHASTIC VARIABILITY IN X-RAY EMISSION FROM THE BLACK HOLE BINARY GRS 1915+105 Yuriy S. Polyakov1 , Joseph Neilsen2,4 , and Serge F. Timashev3,5 2 1 USPolyResearch, Ashland, PA 17921, USA; ypolyakov@uspolyresearch.com MIT Kavli Institute for Astrophysics and Space Research, Cambridge, MA 02139, USA 3 Karpov Institute of Physical Chemistry, Moscow 103064, Russia Received 2012 January 25; accepted 2012 April 16; published 2012 May 14 ABSTRACT We examine stochastic variability in the dynamics of X-ray emission from the black hole system GRS 1915+105, a strongly variable microquasar commonly used for studying relativistic jets and the physics of black hole accretion. The analysis of sample observations for 13 different states in both soft (low) and hard (high) energy bands is performed by flicker-noise spectroscopy (FNS), a phenomenological time series analysis method operating on structure functions and power spectrum estimates. We find the values of FNS parameters, including the Hurst exponent, flicker-noise parameter, and characteristic timescales, for each observation based on multiple 2500 s continuous data segments. We identify four modes of stochastic variability driven by dissipative processes that may be related to viscosity fluctuations in the accretion disk around the black hole: random (RN), power-law (1F), one-scale (1S), and two-scale (2S). The variability modes are generally the same in soft and hard energy bands of the same observation. We discuss the potential for future FNS studies of accreting black holes. Key words: accretion, accretion disks – binaries: close – black hole physics – methods: data analysis – methods: statistical – X-rays: individual (GRS 1915+105) Online-only material: color figure in the power spectra are treated as superpositions of broken power laws and/or Lorentzians (Belloni et al. 2002). The second method is based on the analysis of phase lag spectra for “soft” (lower) and “hard” (higher) energy bands, which are believed to be related to two different sources of X-ray emission in the accretion flow; lags between energy bands can help constrain the emission geometry. The third method uses ratios between X-ray light curves in hard and soft energy bands as a proxy for studying the variability of the energy spectrum over time. One of the essential features observed in many X-ray-emitting astrophysical systems is strong apparent randomness on a range of timescales. While the methods above have proven successful in discerning and quantifying different variability modes, they could not elucidate the origin of the apparent randomness, which remains an open question. For instance, it may be related to variation of some external parameter, such as stochastic fluctuations of the viscosity or mass accretion rate in the accretion disk (Lyubarskii 1997; Titarchuk et al. 2007; Wilkinson & Uttley 2009; Uttley et al. 2011); some complex variability may be due to chaotic (deterministic) fluctuations produced by inner-disk instabilities (Misra et al. 2004, 2006). Other interpretations are possible, and a complete picture of energy-dependent variability has yet to emerge. The origin of apparent randomness in each particular case can be studied using the methods of complexity science, such as deterministic chaos, multifractal theory, flicker-noise spectroscopy (FNS), fractional calculus, cellular automata theory, and other approaches. In this paper, we apply a phenomenological complexity science method to analyze one of the most remarkable astrophysical systems, the Galactic microquasar GRS 1915+105, which consists of a black hole with mass 14 ± 4 M , where M is the mass of the Sun, in a 33.5 day orbit around a donor star with mass 1.2 ± 0.2 M (Greiner et al. 2001). GRS 1915+105 is famous for its relativistic jets (Mirabel & Rodrı́guez 1994) and strong X-ray variability in 14 different states (Belloni et al. 2000; Klein-Wolt et al. 2002; Hannikainen et al. 2005), which are characterized by highly 1. INTRODUCTION A microquasar is a binary system in which a normal star orbits a compact object, either a neutron star or a black hole with a typical mass of 10 solar masses. Under certain conditions, matter may be transferred from the normal star onto the compact object in the form of an accretion disk (for theory, see Shakura & Sunyaev 1973). The subsequent release of gravitational potential energy can lead to the formation of relativistic jets (Blandford & Znajek 1977), as well as strong, variable emission across a broad range of wavelengths, from radio waves to X-ray and γ -ray (Tananbaum et al. 1972; Mirabel & Rodrı́guez 1994). Due to their short characteristic timescales compared to supermassive black holes (i.e., black holes having a typical mass of several million solar masses), microquasars are regarded as key astrophysical objects for studying relativistic jets and the physics of the accretion process (Sams et al. 1996; Mirabel & Rodrı́guez 1999; McHardy et al. 2006). The diverse physical processes in microquasars may be studied indirectly by examining temporal variability in emission across the electromagnetic spectrum. Detailed analysis of temporal variations may be used to estimate the geometry and structure of a binary system and to test various astrophysical models for microquasars (Belloni 2010). This variability is commonly characterized (van der Klis 2006 and references therein) with estimates of the power spectrum of the light curve (Belloni et al. 2002, 2006; Qu et al. 2010; Maccarone et al. 2011; Neilsen et al. 2011), frequency-dependent time lags between different energy bands (Reig et al. 2000), or a combination of total count rate and color–color diagrams (colors are ratios of X-ray count rates in different energy bands; Belloni 2010). In the first method, any apparent noise components and quasi-periodic oscillations 4 Also at Harvard-Smithsonian Center for Astrophysics, Harvard University, Cambridge, MA 02138, USA. 5 Also at Institute of Laser and Information Technologies, Russian Academy of Sciences, Troitsk, Moscow Region 142190, Russia. 1 The Astronomical Journal, 143:148 (10pp), 2012 June Polyakov, Neilsen, & Timashev signals of complex systems characterized by nonlinear interactions, dissipation, and inertia (Timashev et al. 2010a; Timashev & Polyakov 2007). Conceptually, the method separates the analyzed signal into three components: low-frequency regular component corresponding to system-specific “resonances” and their interferential (nonlinear) contributions, stochastic “jump” (random walk) component at larger frequencies corresponding to a dissipative process of anomalous diffusion, and stochastic highest-frequency “spike” component corresponding to inertial (non-dissipative) effects. The idea of separating the stochastic dynamics of a complex nonlinear system into dissipative random-walk and non-dissipative inertial components stems from numerous simulations performed by various cellular automata models demonstrating that the stochastic dynamics of these processes is associated with intermittency, consecutive alternation of rapid (inertial) changes in the values of dynamic variables on small time intervals with small (dissipative) variations of the values on longer time intervals (Bak 1997; Wolfram 2002; Schuster 1995). These studies showed that the occurrence of dissipative and inertial changes in such dynamics can be responsible for the power-law distributions observed in nature, such as Guttenberg–Richter and Zipf–Parreto laws and flicker noise. The FNS method further assumes that correlations present in sequences of different irregularities, such as spikes and jumps, may be treated as main information carriers in the stochastic variability of a complex dissipative system (Timashev & Polyakov 2007). This assumption has much in common with how wavelet analysis evaluates the inverse transform only for those time intervals where the first derivative of the signal experiences steepest changes, i.e., intervals with the highest degree of irregular behavior (Walnut 2002). The spike-like irregularities are naturally associated in FNS with the inertial component, and jump-like irregularities with a dissipative process of anomalous diffusion (random walks), which is illustrated in Figure 2 of Timashev & Polyakov (2007). It should be pointed out that in practical signal analysis, Kolmogorov transient structure functions and power spectrum estimates contain only integral contributions from stochastic irregularities corresponding to a specific evolution hierarchy level and distributed with certain probability densities. If one analyzes the signal at a much higher sampling frequency, each of the irregularities observed at a smaller sampling frequency becomes a complex structure with its own spike- and jump-like irregularities. Further differences and similarities between FNS and other methods for the analysis of the dynamics of complex systems, such as nonlinear dynamics (deterministic chaos), wavelet analysis, and fractional calculus, are discussed in depth elsewhere (Timashev & Polyakov 2007). As FNS is a phenomenological method, its main goal is to estimate from analyzed series certain model-independent stochastic parameters corresponding to general dissipative and inertial processes. These parameters can then be used to develop specific models or refine existing ones. The validity of FNS assumptions for many applications was proven by effectively applying the method to the analysis of structure and dynamics for various electrochemical, geophysical, medical, and physical systems (Descherevsky et al. 2003; Telesca et al. 2004; Hayakawa & Timashev 2006; Timashev 2006; Ida et al. 2007; Timashev & Polyakov 2007; Timashev et al. 2009, 2010a, 2010b, 2012; Ryabinin et al. 2011; Mirsaidov et al. 2011). For instance, it was recently demonstrated by FNS that the dynamics of X-ray emission from GRS 1915+105 and Cygnus X-1, recorded in the timeframe from 1996 January 1 to 2005 December 31 structured low-frequency variability in the X-ray light curve on timescales ranging from seconds to hours, with amplitudes often well in excess of 100% (Belloni et al. 2000). Accordingly, it provides a useful case for studying the underlying processes behind the X-ray variability. Apparent random behavior in the X-ray light curves was examined for 12 of the states using nonlinear time series analysis techniques (Misra et al. 2004, 2006; Harikrishnan et al. 2011). It was suggested that three types of underlying mechanisms may take place: chaotic, stochastic (random), and “nonstochastic” (chaotic + colored noise), which implies that several of the 12 variability states could be explained in terms of a deterministic system of ordinary nonlinear differential equations. The analysis presented by Misra et al. (2004) was based on the estimation of correlation dimension, which is complicated in this case by the presence of high-frequency Poisson noise in all light curves. Moreover, because the correlation dimension can approach a constant value for both chaotic data and colored noise, it is possible that colored noise in the light curve data for GRS 1915+105 may be incorrectly classified as chaotic (Osborne & Provenzale 1989). A more advanced procedure including surrogate data analysis, singular value decomposition (SVD) technique, and correlation dimension was applied by Misra et al. (2006) in an attempt to obtain more conclusive results. The analysis of SVD plots for real and surrogate data revealed some apparent attractors in the real X-ray variability of GRS 1915+105. However, contamination by noise made it difficult to determine whether the attractors were chaotic or periodic (i.e., limit cycles). Correlation entropy and multifractal analysis were added to the nonlinear time series analysis of GRS 1915+105 data by Harikrishnan et al. (2011). This allowed the authors to identify colored noise and estimate the range of fractal dimensions for each observation. It should be noted that in addition to the challenges mentioned above, the light curve data are generally nonstationary and may exhibit a range of characteristic timescales, which further complicates the application of nonlinear time series analysis to these data and makes the results inconclusive, as pointed out by Harikrishnan et al. (2011). Some of these problems can be overcome by utilizing alternative methods of signal analysis for complex systems that are not based on chaos theory. In this study, we apply one such approach, flicker-noise spectroscopy (Timashev et al. 2010a, 2012; Timashev & Polyakov 2007; Timashev 2006, 2007), a phenomenological method operating on Kolmogorov transient structure functions and autocorrelation-based estimates of power spectrum, to examine the apparent randomness in the X-ray light curves of GRS 1915+105, and compare our findings with previously reported results obtained by nonlinear time series analysis (Harikrishnan et al. 2011). In contrast to chaos theory, the FNS method separates between regular and stochastic components, handles multiple characteristic timescales, and provides a mechanism for excluding the effect of Poisson noise at highest frequencies, thus addressing most of the challenges mentioned above. 2. FLICKER-NOISE SPECTROSCOPY Here, we will only deal with the basic FNS relations needed to understand the parameterization procedure. FNS is described in more detail elsewhere (Timashev 2006, 2007; Timashev & Polyakov 2007, 2008; Timashev et al. 2010a). It should be noted that the term “stochastic,” here and throughout in this paper, refers to random variability in the 2 The Astronomical Journal, 143:148 (10pp), 2012 June Polyakov, Neilsen, & Timashev interpolations for both types of stochastic irregularities have the following functional form (Timashev 2006): with a sampling interval of approximately 1.5 hr, manifested a one-scale process of anomalous diffusion (Timashev et al. 2010a). Therefore, it would be logical to apply the FNS method to the analysis of the GRS 1915+105 X-ray emission data at much smaller (subsecond) sampling intervals, which is the topic of this study. In FNS, all introduced parameters for signal V (t), where t is time, are related to the autocorrelation function ψ(τ ) = V (t)V (t + τ )T −τ , Ss (f ) ≈ Ss (0) , 1 + (2πf T0 )n (5) where Ss (0), T0 , and n and are phenomenological parameters, which have different meanings for spike- and jump-like irregularities. In other words, the resultant interpolation is a superposition of expressions (5) for these two components. When the contributions of spikes and jumps to the overall stochastic component are comparable in the highest-frequency range or the contribution of spikes is negligible, which is true for the data analyzed in this study, Equation (5) can be readily used as the resultant interpolation. In this case, Ss (0) is the phenomenological parameter corresponding to the initial (plateau) value of power spectrum estimate, T0 is the correlation time for jumpand spike-like irregularities after which the self-similarity observed in power spectrum estimate breaks down, and n is the flicker-noise parameter characterizing the rate of loss of correlations in the series of high-frequency irregularities in time intervals T0 . It should be noted that the spectral interpolations for spike- and jump-like irregularities were derived using the theory of generalized functions based on different assumptions (Timashev 2006), and generally correspond to different physical processes. The negligible contribution of the spike component in this particular case is attributed to the high level of Poisson noise at highest frequencies, as will be illustrated later in the paper. It should also be pointed out that mathematically Equation (5) represents a superposition of flat spectrum at lower frequencies with power law at higher frequencies (broken power law with a smooth transition), and T0 is related to the reciprocal of the characteristic frequency used in astrophysical spectral analysis (Belloni et al. 2002). The structure function interpolation for the case of one characteristic scale in the sequence of jump-like irregularities (random walks) has the following form (Timashev 2006): (1) where τ is the time lag parameter (0 < τ  TM ), TM is the upper bound for τ (TM  T /2), and T is the averaging window. The angular brackets in relation (1) stand for the averaging over time interval T − τ :  T −τ 1 (· · ·)T −τ = (· · ·) dt. (2) T −τ 0 The averaging over interval T − τ implies that all the characteristics that can be extracted by analyzing functions ψ(τ ) should be regarded as the average values on this interval. To extract the information contained in ψ(τ ) (V (t) = 0 is assumed), the following transforms, or “projections,” of this function are analyzed: cosine transforms (power spectrum estimates) S(f ), where f is the frequency,  TM V (t)V (t + t1 )T −τ cos(2πf t1 ) dt1 , (3) S(f ) = 2 0 and its difference moments (Kolmogorov transient structure functions) of the second order Φ(2) (τ ),   (4) Φ(2) (τ ) = [V (t) − V (t + τ )]2 T −τ . The information contents of S(f ) and Φ(2) (τ ) are generally different, and the parameters for both functions are needed to solve parameterization problems. By considering the intermittent character of signals under study, interpolation expressions for the stochastic components Ss (f ) and Φs (2) (τ ) of S(f ) and Φ(2) (τ ), respectively, were derived using the theory of generalized functions by Timashev (2006). It was shown that the stochastic components of structure functions Φ(2) (τ ) are formed only by jump-like (random walk) irregularities, and stochastic components of functions S(f ), which characterize the “energy side” of the process, are formed by spike-like (inertial) and jump-like irregularities. It should be noted that τ in Equations (1)–(4) is considered as a macroscopic parameter exceeding the sampling period by at least one order of magnitude. This constraint is required to derive the expressions and separate out contributions of dissipative jump-like and inertial (non-dissipative) spike-like components. The basic idea in parameterization is to use two interpolation expressions for the stochastic components. The first expression is used to determine the spectral contribution of the stochastic components of signal V (t) and exclude the contribution of the low-frequency regular component to the parameters related to jump- and spike-like stochastic irregularities. The second interpolation, which deals with the stochastic component in the structure function, is used to determine the parameters characterizing the series of jump-like irregularities, which correspond to anomalous diffusion (Timashev et al. 2010a). For the case of one characteristic scale in the sequences of jump-like and spike-like irregularities, the spectral 2 −1 Φ(2) (H ) · Γ (H, τ/T1 )]2 , s (τ ) ≈ 2σ [1 − Γ (6) ∞ where Γ(s, x) = x exp(−t) · t s−1 dt and Γ(s) = Γ(s, 0) are the complete and incomplete gamma functions, respectively (x  0 and s > 0); σ is the standard deviation of the measured dynamic variable with dimension [V]; H is the Hurst exponent, which describes the rate at which the dynamic variable “forgets” its values on the time intervals that are less than the correlation time T1 . At H = 0.5, interpolation (6) corresponds to Brownian motion (Fickian diffusion); at H < 0.5, to subdiffusion; and at H > 0.5, to superdiffusion (Lévy diffusion or Lévy flights; Timashev et al. 2010a). Subdiffusion and superdiffusion are nonlinear stochastic processes that are generally modeled by fractional Fokker–Planck equations (Tsallis & Lenzi 2002; Chen et al. 2010). Expression (6) at H = 0.5 was shown to correspond to the diffusion equation with integrodifferential boundary conditions (Timashev et al. 2010a). For the case of two characteristic scales, the resultant spectral interpolation can be written as Ss (f ) ≈ Ss1 (0) Ss2 (0) + , 1 + (2πf T01 )n1 1 + (2πf T02 )n2 (7) where indices 1 and 2 in subscripts denote respective scales (Parkhutik & Timashev 2000). 3 The Astronomical Journal, 143:148 (10pp), 2012 June Polyakov, Neilsen, & Timashev Figure 1. First 2500 s record of continuous activity for light curves of 15 sample observations (soft X-ray energy band). The structure function for two characteristic scales takes the following piecewise form: Φ(2) s (τ )   Γ(H1 , τ/T11 ) 2 ≈ 2σ1 1 − , τ < τ0 , Γ(H1 ) 2  Φ(2) s (τ ) ≈ 2σ1 + 2σ2 2 2 1− 0 Γ(H2 , τT−τ ) 12 Γ(H2 ) characterized by a standard deviation within a specific averaging window. In this case, the real signals at specific averaging intervals and sampling frequencies should be regarded as quasistationary with certain values of standard deviation and other phenomenological parameters. It should be noted that the values of phenomenological parameters may vary on different quasistationary intervals. The detailed FNS parameterization algorithms for one- and two-scale cases in discrete form, which were adapted for the X-ray emission data, are described in Appendices A and B, respectively. (8) 2 , τ  τ0 , (9) where indices 1 and 2 in subscripts denote respective scales, and τ0 is the value of τ corresponding to the transition to the second scale. Expressions (8)–(9) were derived based on the assumption that the characteristic scales are much apart from each other, which justifies the use of a linear superposition instead of the more complex expression employed by Parkhutik & Timashev (2000). In some applications, it is also useful to introduce parameters −1 −1 Ss (T0−1 ) (for one-scale case) and Ss (T01 ), Ss (T02 ) (for twoscale case) to account for the “intensity” of jump- and spike-like irregularities on characteristic frequency intervals (Timashev et al. 2012). FNS interpolations were derived for a wide-sense stationary signal in which the phenomenological parameters are the same at each level of the system evolution hierarchy. At the same time, they can also be applied to the analysis of real signals, which are generally nonstationary, but can be 3. RESULTS AND DISCUSSION In order to sample the stochastic variability in GRS 1915+105, we analyze 15 of the observations of this system (13 different X-ray states) made with the Rossi X-Ray Timing Explorer over the last 15 years, which are illustrated in Figure 1. Each observation consists of several ∼2500 s intervals separated by gaps. The specific observations are selected to study at least one long instance of each of GRS 1915+105’s known states (except for the ω state, for which we did not find a long enough continuous data set). From each observation, we extract X-ray light curves in two energy bands (soft: 2–6 keV, and hard: 15–60 keV) with sampling frequencies fd of 16 and 128 Hz. We also extract X-ray background light curves in both energy bands to subtract from each time bin, and we correct the light curve 4 The Astronomical Journal, 143:148 (10pp), 2012 June Polyakov, Neilsen, & Timashev Table 1 Parameters for the Light Curves of 12 GRS 1915+105 Observations with 1S, 1F, or RN Behavior (fd = 16 Hz) Observation ID Class Nci a Soft X-Rays b 90105-06-01-00 40703-01-24-00 60405-01-02-00 20402-01-54-00 20402-01-37-01 92092-03-01-00 80127-01-02-00 30183-02-02-02 20187-02-01-00 50125-01-04-00 50405-01-01-00 20402-01-13-00 4 3 3 3 2 4 4 3 2 3 2 3 γ κ ρ δ λ μ ξ φ α ν ν χ μV (s−1 ) Type PN%c 2790 1443 1548 2134 1409 3814 889 611 581 1102 848 387 1S 1S 1S 1F 1F 1F 1F 1F 1S 1F 1F RN 15 20 15 15 25 20 10 40 10 5 5 ... Hard X-Rays b n H T0 (s) T1 (s) μV (s−1 ) Type PN%c n H T0 (s) T1 (s) 1.73 2.55 2.15 1.2 1.91 1.7 1.47 1.36 1.53 1.29 1.37 ... 0.37 0.77 0.58 0.1 0.46 0.35 0.25 0.21 0.27 0.16 0.18 ... 13 5 6 ... ... ... ... ... 50 ... ... ... 17 5 7 ... ... ... ... ... 86 ... ... ... 122 89 100 160 182 148 44 49 129 107 92 119 1S 1S 1S 1F 1F 1F 1F 1F RN RN RN RN 30 35 35 50 30 45 65 70 ... ... ... ... 2.15 2.61 1.91 1.04 1.57 1.94 3.06 1.88 ... ... ... ... 0.59 0.86 0.47 0.06 0.29 0.48 0.92 0.3 ... ... ... ... 5 3 2 ... ... ... ... ... ... ... ... ... 6 3 3 ... ... ... ... ... ... ... ... ... Notes. a Number of intervals in observation. b Average count rate. c Percentage of Poisson noise in the highest-frequency range—determined for the plot in double logarithmic scale. Table 2 Parameters for the Light Curves of Three GRS 1915+105 Observations with 2S Behavior (fd = 128 Hz; Notes Same as in Table 1) Observation ID Class Nci Band μV (s−1 ) Type n1 H1 T01 (s) T11 (s) n2 H2 T02 (s) T12 (s) 92092-01-01-00 β 3 Soft Hard 1604 121 2S 2S 1.73 1.76 0.53 0.58 0.22 0.16 0.20 0.14 3.07 2.49 1.06 0.66 12 30 10 36 20186-03-02-01 θ 4 Soft Hard 1738 228 2S 2S 2.53 2.19 0.88 0.75 0.13 0.11 0.13 0.10 2.3 2.5 0.65 0.7 34 38 42 46 20402-01-45-02 θ 4 Soft Hard 1783 232 2S 2S 2.38 1.9 0.82 0.66 0.13 0.13 0.13 0.11 2.21 2.68 0.58 0.74 39 35 56 44 noise is significant in the frequency range close to 16 Hz. The other three observations, which have a relatively high signal-toPoisson-noise ratio in that frequency range, are parameterized for a higher sampling frequency fd = 128 Hz (Table 2). The pa−1 −1 rameters Ss (T0−1 ), Ss (T01 ), Ss (T02 ), σ , σ1 , and σ2 are not listed in the tables because their values cannot be reliably determined when a significant Poisson noise level is present. Our analysis shows that the noise behavior for the observations considered in this study exhibits four different modes of stochastic variability: “random” (RN), power-law (1F), onescale (1S), and two-scale (2S). The RN variability, which is illustrated in Figure 2(a), corresponds to values of n less than unity and the overall variation of spectral power in the considered frequency range by less than two orders of magnitude. The 1F variability is characterized by a power-law (1/f n ) dependence at low and intermediate frequencies and a flat Poisson noise level at highest frequencies, as can be seen in Figure 2(b). This implies that the stochastic process for this variability mode exhibits a form of self-similar behavior and is far from reaching the steady state (characteristic timescale) at the averaging interval and frequency range considered here. The 1S variability follows the interpolation (5) at low- and intermediate-frequency ranges and exhibits a Poisson noise at highest frequencies, as illustrated in Figure 2(c). The 2S variability is adequately described by interpolation (7) and contains two characteristic timescales, as can be seen in Figure 2(d). In Figure 3, we present the structure functions (panels (a) and (c)) and their stochastic components (panels (b) and (d)) for sample 1S and 2S variability modes, respectively. count rates for detector dead time. Taking into consideration that all analyzed time series exhibit Poisson noise in the highestfrequency range, we add a term to our interpolations ((5)–(9)) accounting for this noise. For the relatively low frequencies considered here, the expression for the dead-time-corrected Poisson noise level (Zhang et al. 1995; Morgan et al. 1997) can be approximated by a constant value (in general, it is frequency dependent; see Appendix C). Considering only continuous 2500 s intervals in two energy bands, we analyze 94 time series in total. We perform the analysis at TM = T /9 ≈ 280 s to make sure the autocorrelation-based power spectrum estimate is close to the actual power spectrum value. The spectral FNS parameters are also found for the power spectrum estimate calculated by the Welch method (Welch 1967; averaging eight windows of size T /9 with a 50% overlap). The Welch method is similar to methods for computing power spectra commonly used in X-ray astrophysics (Morgan et al. 1997). The difference between the parameter values found by the two methods rarely exceeded 10%, which implies that the autocorrelation-based power spectrum estimate is adequate for our analysis. For the 15 observations analyzed in this study, the variations of FNS parameters within a given observation were not significant in most cases (generally less than 10%), which suggests that the parameterization method is reliable and that specific variability modes last longer than the observations studied here. Table 1 lists the observation-averaged values of parameters for the 12 observations where the contribution of Poisson 5 The Astronomical Journal, 143:148 (10pp), 2012 June Polyakov, Neilsen, & Timashev Figure 2. Typical power spectrum plots for each mode of stochastic variability (thick black line is spectral FNS interpolation; thin gray line is experimental power spectrum estimate): (a) RN, class χ , observation 20402-01-13-00, 1st interval, soft X-ray energy band, n ≈ 0.52; (b) 1F, class μ, observation 920920-03-01-00, 3rd interval, soft X-ray energy band; (c) 1S, class γ , observation 90105-06-01-00, 3rd interval, soft X-ray energy band; (d) 2S, class θ , observation 20402-01-45-02, 4th interval, soft X-ray energy band. Figure 3. Sample structure function plots for variability modes 1S and 2S (1 (blue)—experimental curve; 2 (green)—interpolation; 3 (red)—resonant component of structure function): (a) structure function plot for 1S, class γ , observation 90105-06-01-00, 3rd interval, soft X-ray energy band and (b) its stochastic component of structure function; (c) structure function plot for 2S, class θ , observation 20402-01-45-02, 4th interval, soft X-ray energy band and (d) its stochastic component of structure function. (A color version of this figure is available in the online journal.) 6 The Astronomical Journal, 143:148 (10pp), 2012 June Polyakov, Neilsen, & Timashev The stochastic component is obtained by subtracting the resonant component from the overall structure function. For the 1S mode, the stochastic component is described by an anomalous diffusion process reaching a steady state (Figure 3(b); Timashev et al. 2010a), while for the 2S case a two-scale process is implied (Figure 3(d)). It can be seen that the parameters for two instances of the same class (ν in Table 1 and θ in Table 2) stay approximately same, which suggests that a particular stochastic variability (certain variability mode and ranges for parameter values) may characterize a class of observations. However, an exhaustive analysis of all X-ray emission observations for GRS 1915+105 is needed to draw definite conclusions in this context. For most of analyzed intervals corresponding to variability modes 1F, 1S, and 2S, we noticed that n ≈ 2H + 1 (n1 ≈ 2H1 + 1, n2 ≈ 2H2 + 1), as expected for fractional Brownian motions (Malamud & Turcotte 1999). We also observed that the characteristic timescales in power spectrum estimate and structure function were close to each other for most of the intervals characterized by variability modes 1S and 2S. Both of these facts imply that the random-walk component plays the dominant role in the stochastic variability of X-ray emission, and the contribution of highest-frequency inertial spike irregularities in these data is not significant (Timashev & Polyakov 2007). This can most likely be attributed to the “erasing” effect of Poisson noise at highest frequencies, which is illustrated in Figure 2 and by the values of parameter PN% in Table 1. Tables 1 and 2 as well as interval-specific analysis suggest that there is no statistically significant difference between the FNS parameter values for soft and hard X-ray energy bands in the same interval. The only differences are in the count rate magnitude (it is smaller for hard X-rays) and the frequency range dominated by the Poisson noise (PN% is higher for hard X-rays). Soft X-ray emission is nominally associated with the thermal spectrum of the accretion disk around the black hole, while hard X-ray emission is generally believed to be related to Compton scattering in a hot electron corona around the disk. The fact that the FNS parameters for these energy bands are close in value could suggest that at these frequencies, there is a common mechanism producing stochastic variability. Future analysis in more energy bands could help determine whether the FNS parameters are truly energy independent. Let us compare the classification in Tables 1 and 2 with the results obtained by the nonlinear time series analysis of total count rates for sample observations corresponding to 12 different states (Harikrishnan et al. 2011, Table 1). As total count rates are dominated by the soft X-ray energy band, we compare the results of Harikrishnan et al. (2011) with the soft X-ray classification presented in this paper. The noise variability for class χ is characterized as random by both analyses. Class φ is marked as random by Harikrishnan et al. (2011), but is assigned to 1F here. The difference, in our opinion, is due to the significant effect of the Poisson noise within half of the considered frequency range (high values of PN%), which was not separately taken into account in the nonlinear time series analysis. Class γ was marked as random by Harikrishnan et al. (2011), but is categorized as 1S by FNS. This discrepancy could be due to the plateau in the low-frequency range of the power spectrum (Figure 2 in our study; Harikrishnan et al. 2011, Figure 10): although it may mimic white noise in the nonlinear time series analysis, in FNS the plateau corresponds to the loss of correlations (at τ > T0 ) within a sequence of stochastic irregularities. Class δ is categorized as random by Harikrishnan et al. (2011), but is marked as 1F in FNS. Considering that the value of n is close to unity in this case (boundary range), additional sample observations for this class would be needed to determine which method is more accurate. Classes κ, λ, and μ were labeled by Harikrishnan et al. (2011) as “deterministic nonlinear + colored noise.” All three are marked as 1F or 1S in this study, which corresponds to a variant of “colored noise.” In classes where there are strong periodic or quasiregular (“resonant”) components (e.g., α, β, θ , ν, κ, λ, and ρ), it is difficult to compare our categorizations directly to the results of Misra et al. (2006) and Harikrishnan et al. (2011), because in contrast to those studies, we subtract the resonant components prior to performing our analysis. In other words, the regular (deterministic) behavior observed in SVD plots (Misra et al. 2006; Harikrishnan et al. 2011) could be dominated by the resonant periodicities, rather than the stochastic component studied here. Moreover, classes β and θ contain two separate stochastic scales in different frequency ranges, which would be difficult to capture in terms of chaos theory, originally developed for systems of nonlinear ordinary differential equations with one characteristic scale. The above analysis demonstrates that the random-walk component plays the dominant role in the stochastic variability of modes 1F, 1S, and 2S. The random-walk component can be interpreted in terms of a dissipative process of anomalous diffusion (Timashev et al. 2010a). It has been argued that stochastic fluctuations of the viscosity in the accretion flow may be responsible for X-ray variability on a wide range of timescales, as well as the presence of multiple characteristic timescales (Lyubarskii 1997; Titarchuk et al. 2007; Wilkinson & Uttley 2009; Uttley et al. 2011). This is because variations in the viscosity lead to changes in the accretion rate, which propagate toward the black hole, resulting in flickering X-ray emission. In principle, the radial dependence of the amplitude of the viscous fluctuations can determine the slope of the power spectrum, and different components in the accretion flow (e.g., corona, inner accretion disk, and outer accretion disk) can contribute at different frequencies. The stochastic modes described here may be related to these viscous fluctuations in the accretion disk, in which case a more comprehensive analysis with FNS could provide a useful characterization of X-ray states from the perspective of stochastic variability. 4. CONCLUDING REMARKS Our analysis shows that the consideration of different frequency ranges in FNS allows one to analyze the stochastic variability in complex signals of GRS 1915+105, despite the significant contribution of Poisson noise in the highest-frequency range, high low-frequency quasi-regular variability, and presence of multiple stochastic timescales. As a proof of concept for future study, the examination of sample observations for 13 different states demonstrates that there are at least four types of stochastic variability in the X-ray emission of GRS 1915+105: random (RN), power-scale (1F), one-scale (1S), and two-scale (2S). The last three are related to random walk processes interpretable in terms of a dissipative process of anomalous diffusion. In this case, the random walk processes could be related to viscosity fluctuations in the accretion disk. Other stochastic variability modes are also possible in the frequency ranges under consideration. For example, it is natural to expect a stochastic variability mode representing an intermediate category between 1S and 2S, i.e., a superposition of one full scale (1S) and self-similar (1F) behavior in different frequency 7 The Astronomical Journal, 143:148 (10pp), 2012 June Polyakov, Neilsen, & Timashev where q = 0. . .M. For q = 1. . .M − 1, Sd (q) should be multiplied by 2, which is the standard procedure for discrete Fourier transforms to take into account the spectral values in the second half of the frequency range. Here, relations q = 2ffd−1 M and Sd (q) = S(f ) × fd describe the equivalence between the discrete and continuous forms of power spectrum estimate. 5. Calculate Ssd (0) as the average value of the power spectrum for the points 1 and 2 (point 0, which corresponds to the zero frequency, is not used in calculating Ssd (0)): ranges. This is borne out by astrophysical observations, which exhibit a wide range of extremely complex power spectra, and are frequently modeled as the sum of a number of different components (Belloni et al. 2002; van der Klis 2006). An exhaustive analysis of all available X-ray emission observations for GRS 1915+105 is needed to get a complete list of variability modes, obtain more accurate values of FNS parameters for each known class of the binary system, and revisit the established classification from the stochastic perspective. A similar analysis of other accreting black holes could help to identify the effect of the extreme variability of GRS 1915+105 on the FNS parameters, and to determine the true relationship of these parameters to stochastic processes in the disk. It would also be interesting to study the variability at other wavelengths, e.g., in observations of the radio jet with the Ryle Telescope (Pooley & Fender 1997; Fender et al. 1997) or the EVLA, via FNS parameterization and cross-correlation with X-rays, in order to compare stochastic variability in different accretion processes. Ssd (0) = by the method of nonlinear least-squares fitting to determine the values of parameters n and T0d . The constant term PS is an estimate of the Poisson noise level calculated as the average value of Sd for 100 highest-frequency points. The fitting is done on the basis of a double logarithmic scale, dividing the entire series into a set of equal intervals (in the calculations presented in this study we took 200 equal intervals). We used the trust-region algorithm for nonlinear square fitting, which is built in MATLAB v.7 or higher (Branch et al. 1999). 7. Separate out the resonant component: APPENDIX A PARAMETERIZATION ALGORITHM FOR ONE-SCALE CASE IN DISCRETE FORM Srd (q) = Sd (q) − Ssd (q), q = 0. . .M. Consider a time series Vd (k). The subscript d here and below is used to denote the discrete form of expressions. Let Nt be the number of points corresponding to the selected averaging interval T, M be the number of points used in estimating the autocorrelation function. In this case, the parameterization procedure can be written as follows. 1. Calculate the arithmetic mean for the signal: 1 Nt ψrd (p) = (A1) k=1 2. Subtract the arithmetic mean from the series Vd (k): − Vd (k) = Vd (k) − μV . Nt −p − (A8) 9. Calculate the difference moment for the resonant component: − 1 Nt − p 1 {Srd (0) + Srd (M)(−1)p } 2M ⎧ ⎫ M−1 ⎨ 1 π pq ⎬ + Srd (q) cos 2 . 2M ⎩ q=1 M ⎭ (A2) 3. Calculate the autocorrelation function for the series Vd : ψd (p) = (A7) 8. Calculate the autocorrelation function for the resonant component as the inverse discrete cosine transform of Srd (q). When q = 1. . .M − 1, divide Srd (q) by 2 to take into account the spectral values in the second half of the frequency range. Then calculate the inverse cosine transform: Nt Vd (k). (A5) 6. Interpolate |Sd (q)| given by Equation (A4) using the expression Ssd (0) + PS (A6) Ssd (q) = 1 + (π Mq T0d )n J.N. was supported in part by the Harvard University Graduate School of Arts and Sciences, Chandra grant AR0-11004X, and the National Aeronautics and Space Administration through the Smithsonian Astrophysical Observatory contract SV3-73016 to MIT for support of the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under the contract NAS8-03060. μV = Sd (1) + Sd (2) . 2 Φ(2) rd (p) = 2[ψrd (0) − ψrd (p)], p = 0. . .M. − Vd (k) Vd (k + p), p = 0. . .M. k=1 (A9) (2) The continuous equivalent of Φ(2) rd (p) is Φr (τ ). 10. Calculate the difference moment for the experimental series: (A3) The autocorrelation interval M should not be higher than Nt /4 (higher values of M will result in the loss of statistical information in estimating the autocorrelation function). To go from discrete form to the continuous one, one can use the following expression: p = Nt τ/T . 4. Calculate the discrete cosine transform of the autocorrelation function: Φ(2) d (p ) 1 = Nt − p Nt −p  − Vd (k) − − Vd (k 2 + p) . (A10) k=1 11. Calculate the difference moment for the random component: (2) (2) (A11) Φ(2) esd (p) = Φd (p) − Φrd (p). M−1 πqp Sd (q) = ψd (0) + ψd (M)(−1)q + 2 , ψd (p) cos M p=1 (2) The continuous equivalent of Φ(2) esd (p) is Φes (τ ). (A4) 8 The Astronomical Journal, 143:148 (10pp), 2012 June Polyakov, Neilsen, & Timashev The parameter τ0d is estimated as the value of p after which 2 Φ(2) esd (p) − PΦ exceeds 2σ1 by at least 3%. Then we use interpolation   2 p−τ0d , ) Γ(H 2 T (2) 12d Φsd (p) = 2σ12 + 2σ22 + PΦ (B3) 1− Γ(H2 ) 12. Determine the parameters σ , H, T1d , PΦ by fitting Φ(2) esd (p) into Equation (A11) to the interpolation expression of the anomalous diffusion type (Timashev et al. 2010a): 2 −1 2 Φ(2) sd (p) = 2σ × [1 − Γ (H ) × Γ(H, p/T1d )] + PΦ , (A12) ∞ where Γ(s, x) = x exp(−t)t s−1 dt, Γ(s) = Γ(s, 0), using the same least-squares fitting method as in step 6. Here, the constant term PΦ is an estimate of the Poisson noise contribution to the structure function. −1 13. Calculate Ssd (T01d ) using Equation (A6). 14. After the values of all six FNS parameters—σ , T0d , T1d , −1 H, n, Ssd (T0d )—are determined, calculate the dimensional −1 values for T0d , T1d , Ssd (T0d ): T0 = T0d ×Δt, T1 = T1d ×Δt, −1 −1 Ss (T0 ) = Ssd (T0d ) × Δt, where Δt = fd−1 . 15. Calculate the relative error Φ in the interpolation of difference moment Φ(2) d (p): for p > τ0d to find the values of σ2 , H2 , and T12d . −1 −1 At step 13, calculate Ssd (T01d ) and Ssd (T02d ) using Equation (B1). At step 14, calculate the dimensional values for T01d , T02d , −1 −1 T11d , T12d , Ssd (T01d ), Ssd (T02d ): T01 = T01d ×Δt, T02 = T02d ×Δt, −1 −1 T11 = T11d × Δt, T12 = T12d × Δt, Ss (T01 ) = Ssd (T01d ) × Δt, −1 −1 −1 Ss (T02 ) = Ssd (T02d ) × Δt, where Δt = fd . Step 15 is the same as in the one-scale case (Appendix A). APPENDIX C APPROXIMATION FOR DEAD-TIME-CORRECTED POISSON NOISE LEVEL EXPRESSION  M    (2)  (2) (2) Φd (p) − Φrd (p) − Φsd (p) Φ = p=1 M  p=1 The formula for computing the Poisson power corrected for instrument dead time is given in Zhang et al. (1995) and Morgan et al. (1997). Adapting this for power spectra based on the autocorrelator, we have     τd SP (f ) = A × 2 1 − 2r0 τd 1 − 2τb   (Nt − 1) τd r0 τd −2 cos(2π tb f ) Nt tb  2 f ) sin(π τ vle 2 , (C1) + 2rvle r0 τvle π τvle f × 100%. (A13) Φ(2) d (p) Here, the error is determined as the ratio of the difference of areas between the experimental structure function and the total interpolation function to the area of the experimental structure function. The areas are calculated by numerical integration using the rectangle method because the original series have a rather large number of points. The parameterization is successful if Φ  10% (Timashev et al. 2010a). APPENDIX B PARAMETERIZATION ALGORITHM FOR TWO-SCALE CASE IN DISCRETE FORM where A is a scale factor between FNS units and standard astrophysical power spectra (where pure Poisson noise has SP (f ) = 2), r0 is the count rate per proportional counter unit (PCU), Nt is the number of frequency points, rvle is the rate of very large events (VLE; e.g., cosmic rays) per PCU, tb is the bin size, τd is the dead time per event, τvle is the dead-time window for each VLE. In our case, rvle = 200 counts−1 s−1 PCU−1 , tb = 0.0625 or 0.0078 s, τd = 10 μs, and τvle = 150 μs. It is more convenient for analysis to write Equation (C1) as Steps 1–5 are the same as in the one-scale case (Appendix A). In step 6, the following interpolation is used: Ssd (q) = Ss1d (0) Ssd (0) − Ss1d (0) + + PS . q n1 1 + (π M T01d ) 1 + (π Mq T02d )n2 (B1) Here, the method of nonlinear least-squares fitting is applied in log–log scale to determine the values of Ss1d (0), T01d , n1 , T02d , and n2 . Steps 7–11 are the same as in the one-scale case (Appendix A). In step 12, we first use the following interpolation to estimate the values of parameters for the first scale:  Γ(H1 , p/T11d ) (2) 2 Φsd (p) = 2σ1 × 1 − Γ(H1 )   2 σ2 Γ(H2 , p/T12d ) + + PΦ . (B2) 1− σ1 Γ(H2 ) SP (f ) = B + C cos[Df ] + E, where (C2)    τd , B = 2A 1 − 2r0 τd 1 − 2τb   τd (Nt − 1) , C = −2A r0 τd Nt tb D = 2π tb ,   sin(π τvle f ) 2 2 ≈ 2Arvle r0 τvle . π τvle f It can be easily shown for our data that B C + E (usually by at least two orders of magnitude). This implies that in our study we can use the constant-value approximation 2 E = 2Arvle r0 τvle The nonlinear least-squares fitting method in log-log scale is applied to find the values of parameters σ1 , H1 , T11d , σ2 , H2 , and T12d . The parameter PΦ is estimated as the value of Φ(2) esd (p) at p = 1. The set corresponding to the lower value of T11d and T12d is used as the values for σ1 , H1 , and T11d . SP (f ) ≈ B. 9 (C3) The Astronomical Journal, 143:148 (10pp), 2012 June Polyakov, Neilsen, & Timashev REFERENCES Pooley, G. G., & Fender, R. P. 1997, MNRAS, 292, 925 Qu, J. L., Lu, F. J., Lu, Y., et al. 2010, ApJ, 710, 836 Reig, P., Belloni, T., van der Klis, M., et al. 2000, ApJ, 541, 883 Ryabinin, G. V., Polyakov, Y. S., Gavrilov, V. A., & Timashev, S. F. 2011, Nat. Hazards Earth Syst. Sci., 11, 541 Sams, B. J., Eckart, A., & Sunyaev, R. 1996, Nature, 382, 47 Schuster, H. G. 1995, Deterministic Chaos: An Introduction (Weinheim: WileyVCH) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 Tananbaum, H., Gursky, H., Kellogg, E., et al. 1972, ApJ, 177, L5 Telesca, L., Lapenna, V., Timashev, S., Vstovsky, G., & Martinelli, G. 2004, Phys. Chem. Earth, 29, 389 Timashev, S. F. 2006, Russ. J. Electrochem., 42, 424 Timashev, S. F. 2007, Fliker-Shumovaya Spektroskopiya: Informatsiya v khaoticheskikh signalakh (Flicker-Noise Spectroscopy: Information in Chaotic Signals) (Moscow: Fizmatlit) Timashev, S. F., Panischev, O. Y., Polyakov, Y. S., Demin, S. A., & Kaplan, A. Y. 2012, Physica A, 391, 1179 Timashev, S. F., & Polyakov, Y. S. 2007, Fluctuation Noise Lett., 7, R15 Timashev, S. F., & Polyakov, Y. S. 2008, Int. J. Bifurcation Chaos, 18, 2793 Timashev, S. F., Polyakov, Y. S., Misurkin, P. I., et al. 2010a, Phys. Rev. E, 81, 041128 Timashev, S. F., Polyakov, Y. S., Yulmetyev, R. M., et al. 2009, Laser Phys., 19, 836 Timashev, S. F., Polyakov, Y. S., Yulmetyev, R. M., et al. 2010b, Laser Phys., 20, 604 Titarchuk, L., Shaposhnikov, N., & Arefiev, V. 2007, ApJ, 660, 556 Tsallis, C., & Lenzi, E. 2002, Chem. Phys., 284, 341 Uttley, P., Wilkinson, T., Cassatella, P., et al. 2011, MNRAS, 414, L60 van der Klis, M. 2006, in Rapid X-ray Variability, ed. W. H. G. Lewin & M. van der Klis (Cambridge: Cambridge University Press), 39 Walnut, D. F. 2002, An Introduction to Wavelet Analysis (Boston, MA: Birkhäuser) Welch, P. 1967, IEEE Trans. Audio Electroacoust., 15, 70 Wilkinson, T., & Uttley, P. 2009, MNRAS, 397, 666 Wolfram, S. 2002, A New Kind of Science (Champaign, IL: Wolfram Media) Zhang, W., Jahoda, K., Swank, J. H., et al. 1995, ApJ, 449, 930 Bak, P. 1997, How Nature Works: The Science of Self-organized Criticality (Oxford: Oxford Univ. Press) Belloni, T. (ed.) 2010, The Jet Paradigm: From Microquasars to Quasars (Berlin: Springer) Belloni, T., Klein-Wolt, M., Méndez, M., et al. 2000, A&A, 355, 271 Belloni, T., Psaltis, D., & van der Klis, M. 2002, ApJ, 572, 392 Belloni, T., Soleri, P., Casella, P., et al. 2006, MNRAS, 369, 305 Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433 Branch, M. A., Coleman, T., & Li, Y. 1999, SIAM J. Sci. Comput., 21, 1 Chen, W., Sun, H., Zhang, X., et al. 2010, Comput. Math. Appl., 59, 1754 Descherevsky, A. V., Lukk, A. A., Sidorin, A. Y., Vstovsky, G. V., & Timashev, S. F. 2003, Nat. Hazards Earth Syst. Sci., 3, 159 Fender, R. P., Pooley, G. G., Brocksopp, C., et al. 1997, MNRAS, 290, L65 Greiner, J., Cuby, J. G., & McCaughrean, M. J. 2001, Nature, 414, 522 Hannikainen, D., Rodriguez, J., Vilhu, O., et al. 2005, A&A, 435, 995 Harikrishnan, K. P., Misra, R., & Ambika, G. 2011, Res. Astron. Astrophys., 11, 71 Hayakawa, M., & Timashev, S. F. 2006, Nonlinear Process. Geophys., 13, 255 Ida, Y., Hayakawa, M., & Timashev, S. 2007, Nat. Hazards Earth Syst. Sci., 7, 479 Klein-Wolt, M., Fender, R. P., Pooley, G. G., et al. 2002, MNRAS, 331, 745 Lyubarskii, Y. E. 1997, MNRAS, 292, 679 Maccarone, T. J., Uttley, P., van der Klis, M., et al. 2011, MNRAS, 413, 1819 Malamud, B. D., & Turcotte, D. L. 1999, J. Stat. Plan. Inference, 80, 173 McHardy, I. M., Koerding, E., Knigge, C., et al. 2006, Nature, 444, 730 Mirabel, I. F., & Rodrı́guez, L. F. 1994, Nature, 371, 46 Mirabel, I. F., & Rodrı́guez, L. F. 1999, ARA&A, 37, 409 Mirsaidov, U., Timashev, S. F., Polyakov, Y. S., et al. 2011, Analyst, 136, 570 Misra, R., Harikrishnan, K. P., Ambika, G., et al. 2006, ApJ, 643, 1114 Misra, R., Harikrishnan, K. P., Mukhopadhyay, B., et al. 2004, ApJ, 609, 313 Morgan, E. H., Remillard, R. A., & Greiner, J. 1997, ApJ, 482, 993 Neilsen, J., Remillard, R. A., & Lee, J. C. 2011, ApJ, 737, 69 Osborne, A., & Provenzale, A. 1989, Physica D, 35, 357 Parkhutik, V., & Timashev, S. 2000, Russ. J. Electrochem., 36, 1221 10 View publication stats