Chaotic and stochastic processes in the accretion flows of the black hole Xray binaries revealed by recurrence analysis
Center for Theoretical Physics, Polish Academy of Sciences, Al. Lotnikow 32/46, 02668 Warsaw, Poland
email: psukova@cft.edu.pl
Received: 8 June 2015
Accepted: 2 December 2015
Aims. Both the well known microquasar GRS 1915+105, as well as its recently discovered analogue, IGR J170913624, exhibit variability that is characteristic of a deterministic chaotic system. Their specific kind of quasiperiodic flares that are observed in some states is intrinsically connected with the global structure of the accretion flow, which are governed by the nonlinear hydrodynamics. One plausible mechanism that is proposed to explain this kind of variability is the thermalviscous instability that operates in the accretion disk. The purely stochastic variability that occurs because of turbulent conditions in the plasma, is quantified by the power density spectra and appears in practically all types of sources and their spectral states.
Methods. We pose a question as to whether these two microquasars are one of a kind, or if the traces of deterministic chaos, and hence the accretion disk instability, may also be hidden in the observed variability of other sources. We focus on the black hole Xray binaries that accrete at a high rate and are, therefore, theoretically prone to the development of radiation pressureinduced instability. To study the nonlinear behaviour of the Xray sources and distinguish between the chaotic and stochastic nature of their emission, we propose a novel method, which is based on recurrence analysis. Widely known in other fields of physics, this powerful method is used here for the first time in an astrophysical context. We estimate the indications of deterministic chaos quantitatively, such as the Rényi’s entropy for the observed time series, and we compare them with surrogate data.
Results. Using the observational data collected by the RXTE satellite, we reveal the oscillations pattern and the observable properties of six black hole systems. For five of them, we confirm the signatures of deterministic chaos being the driver of their observed variability.
Conclusions. We test the method and confirm the deterministic nature of variability in the microquasars GRS 1915+105 and IGR J170913624. Furthermore, we also find significant traces of nonlinear dynamics in three other sources: GX 3394, XTE J1550564 and GRO J165540, particularly in the diskdominated soft state, as well as in the intermediate states at the rising and declining phase of the outburst. Only for the source XTE J1650500 no observation with such variability pattern was found. This is possibly due to the global accretion rate in this source being too small for the limitcycle instability to develop.
Key words: black hole physics / accretion, accretion disks / chaos / Xrays: binaries
© ESO, 2016
1. Introduction
The high energy radiation emitted by black hole Xray binaries originates in an accretion disk, where the viscous dissipation of the gravitational potential energy is responsible for both heating the plasma and transporting the angular momentum. Hardly any of the observed black hole systems exhibit the radiation flux being constant with time. Instead, most of the sources undergo fast and complicated variability patterns on different timescales. The variations that are purely stochastic in their nature are expected since the viscosity of the accretion disk is connected with its turbulent behaviour, which is induced by magnetic instabilities. Furthermore, when the variations occur close to some frequencies, the excess in the power spectrum known as the quasiperiodic oscillations (QPOs) appears. The mechanism that creates QPOs is generally not known and different scenarios for their origin have been proposed, including: i) oscillations of the accretion disk/torus (Titarchuk & Osherovich 2000; Chakrabarti & Manickam 2000); ii) shock fronts inside the accreting material (Das 2003; Suková & Janiuk 2015a); and iii) nonlinear resonances between the radial and orbital epicyclic motions (Abramowicz et al. 2006), which possibly also link the high frequency and low frequency QPOs (see the review by Belloni & Stella 2014). Among the nonresonant class of models, the LenseThirring precession is often invoked and used to explain the twin kHz QPOs as well as the low frequency QPOs, and was originally proposed for the neutron star Xray binaries in their Z and Atoll classes (Stella & Vietri 1998).
In the models where the type of global conditions occur whereby the system finds itself in an unstable configuration, large amplitude fluctuations around the fixed point solution are induced. The variability of the disk that reflects its global evolution and that is governed by the nonlinear differential equations of hydrodynamics, will not be purely stochastic. Instead, the observed behaviour of the disk is characterised by a process, the nature of which is deterministic chaos. This type of model, which explains the quasiperiodic flares with the global evolution of the accretion disk, rather than focusing on the flow oscillations in some particular radii, is not based on the effects of resonances.
The flarelike events of the microquasars, such as those observed in IGR 170913624 or GRS 1915+105 in their “heartbeat” states, are very strong, almost coherent, lowfrequency QPOs. In these sources, they are so orderly and quasicoherent that it is possible to see the oscillations, even directly in their light curves. Generally, it is not possible to identify the QPOs by simply looking at their light curves. Usually the Fourier transform is applied to the light curve to study the power spectrum and to identify the QPO as a small peak, or a feature, just up from the Poisson noise of the power spectrum. Here we apply a novel mathematical method to study the quasiperiodicity and to determine the nature of the variability of the Xray sources, i.e., if it is stochastic or deterministic. In the latter case, we propose that the underlying mechanism of the deterministic chaotic behaviour of the accretion flow is an underlying system of nonlinear equations that governs its global evolution.
1.1. Accretion disk instability due to the nonlinear hydrodynamics
The problem of chaotic variations in Xray binaries accretion disks is interesting and still unsolved from both a theoretical and observational point of view. First, the classic theory of accretion, as proposed by Shakura & Sunyaev (1973) and Lightman & Eardley (1974), predicts that the disk is unstable and undergoes limitcycle oscillations if the viscous stress tensor scales with the total (gas plus radiation) pressure, and the global accretion rate is large enough for the radiation pressure to dominate.
In general, the accretion disks may undergo limitcycle oscillations around a fixed point because of the two main types of thermalviscous instabilities. These are induced either by the domination of radiation pressure in the innermost parts of the disk, which occurs for high accretion rates, or by the partial ionization of hydrogen in the outer disk parts, at appropriate temperatures. Both these instability types have been known for over 40 years in theoretical astrophysics.
Hydrogen ionization instability operates in the outer regions of the accretion disk, where the temperatures are in the range of log T = 3.5−4 [K] and the hydrogen is partially ionized. Under such conditions the opacities in the plasma depend inversely on density and temperature, and nonlinear processes in the accretion flow are induced. This type of instability leads to outbursts of dwarf novae, and soft Xray transients, and generally its characteristic timescales are on the order of months in the case of stellarmass accreting objects (for a review see Lasota 2001).
Another type of nonlinear process that operates on much shorter timescales, is the instability that is induced by dominant radiation pressure. In the classical theory of Shakura & Sunyaev (1973), the accretion flow structure is based on α prescription for viscous energy dissipation. It assumes that the nonzero component T_{rφ} of the stress tensor is proportional to the total pressure. The latter includes the radiation pressure component, which scales with temperature as T^{4} and blows up in hot disks for large accretion rates. This in turn affects the heating and cooling balance, between the energy dissipation and radiative losses. If the accretion rate is small, then most of the disk is dominated by gas pressure and stable. For large enough accretion rates, a zone appears where some of its annuli are dominated by radiation pressure and, therefore, unstable. The larger the global accretion rate, the more annuli are affected by the instability, starting at the inner edge which is the hottest.
If there was no stabilizing mechanism, the disk that is dominated by radiation pressure will not survive. This is because, in this kind of solution, the decreasing density leads to runaway temperature growth. Consequently, the local accretion rate increases and more material is transported inwards. The disk annulus empties because of both the increasing accretion rate and the decreasing density, so there is no selfregulation of the disk structure.
However, the so called “slimdisk” solution, where advection of energy provides additional source of cooling, acts as a stabilising branch. In this way, the advection allows the disk to survive and oscillate between the hot and cold states. This kind of oscillating behaviour leads to periodic changes in the disk luminosity.
The theoretical calculations of the global, longterm disk evolution that are based on the vertically averaged structure, clearly show the oscillatory type of behaviour with characteristic limitcycles (Janiuk et al. 2002). The numerical computations performed in the 3D shearingbox simulations using various codes either deny (Krolik et al. 2007) or confirm (Jiang et al. 2013; Hirose et al. 2009) that the magnetohydrodynamic stress tensor scales with the total pressure. Therefore the issue of αviscosity parametrization and contribution of the radiation pressure to the energy dissipation is still an open problem. Moreover, the 3D simulations do not tackle the global structure of the flow, and we are therefore not able to verify them with the observations of the astrophysical sources directly.
1.2. Observational support for the accretion instability
For many years, the only source that undoubtedly showed oscillations characteristic of the accretion disk instability owing to the dominant radiation pressure, was the microquasar GRS 1915+105 in some of its observed states (Belloni et al. 2000; Janiuk et al. 2002). The source was thus thought to be one of a kind, until the very recent discovery of another microquasar, IGR J170913624, which was found to be an analogue of the previously known microquasar (Capitanio et al. 2012; Altamirano & Belloni 2012). The recent hydrodynamical simulations of the global accretion disk evolution confirmed that the quasiperiodic flarelike events observed in IGR J17091, in some of its variability classes, are also in good quantitative agreement with the radiation pressure instability model (Janiuk et al. 2015).
However, for sources other than the two microquasars mentioned above, no such detailed analysis was performed, neither observationally nor theoretically. However, as proposed by Janiuk & Czerny (2011), at least eight of the known BH Xray binaries should have large enough Eddington accretion rates for radiation pressure instability to develop. These other black hole binaries can also be promising objects for radiation pressure instability, because their Eddington ratio is above ~0.15. The particular value of this critical accretion rate depends on whether the viscous torque parametrization is adopted in a modified version, such as α, instead of αP_{tot}. Nevertheless, for the stellar mass black hole sources, the detailed analysis and verification of available observations of various sources should give constraints for the presence of nonlinear behaviour of their accretion flows.
Finally, because the instability timescale is scaling with the black hole mass, the intermittency in quasars on timescales of hundreds to thousands of years is likely to be of a similar origin as found in the microquasars. Here, some indirect arguments are also proposed for the intermittent type of activity in the supermassive black holes, e.g. from the studies of radio maps (KunertBajraszewska et al. 2010). Also, statistical studies were undertaken and provide useful evidence for episodic activity in AGN (Czerny et al. 2009).
1.3. Our current approach
In the current work, we aim to tackle the problem of the stochastic versus the deterministic nature of the black hole accretion disk variability from an analytical and observational point of view. We perform the recurrence analysis, which is a powerful tool, which is used to study time series and known to work in a broad range of applications, ranging from economy to geophysics and medicine (Eckmann et al. 1987; Marwan et al. 2007). The recurrence analysis method for processing the time series was also used recently to study chaos in motion of test particles (Kopáček et al. 2010; Suková 2011; Suková & Semerák 2012). Here, for the first time, we show that this method has great potential for studying time series observed from astrophysical sources.
Esssentially, the recurrence analysis method aims to reveal the dynamical parameters of the system from the observed time series. These parameters, such as the Rényi’s entropy (see Eq. (A.4) and the following text) or the maximal Lyapunov exponent, which characterises the rate of stretching or contracting of the attractor (Ott 2002), may give an indication that the underlying variability is deterministic in nature. This is possible if the structure evolution equations of the dynamical system do not contain time explicitly, but are nonlinear and have unstable steady state solutions. The socalled “recurrence plot” contains information about time correlation and its form has of an array of points in an NxN square for a time series u_{i}, where i = 1 ...N. The time series is used to construct the orbits x(i) in the system’s ddimensional phase space (where x may relate to the system’s physical state, such as its density or temperature, and d relates to the number of nonlinear equations that govern the evolution of its structure). The point in the recurrence plot is marked whenever the trajectory returns close to itself, so that x(j) is close to x(i), i.e. closer than some assumed radius ϵ of an ddimensional sphere. The plot is by definition symmetric with respect to the diagonal i = j, and the diagonal lines are not infinite, if only the variability is not completely periodic.
We study the occurrence of long diagonal lines in the recurrence plots of observed data series and compare it to the series of surrogate data. The latter are produced numerically to have the same mean and variance (shuffled surrogates, see below Sect. 3.1) or even the same power spectrum (IAAFT surrogates), albeit created by random processes as the permutation of the original time series. The surrogate data have already been used by Misra et al. (2006) as a countercheck with the correlation dimension technique and the singlevalue decomposition technique in the context of the nonlinear variability in the microquasar GRS 1915+105. Here too, our estimated dynamical invariants, e.g. the second order Rényi’s entropy, which is a measure of the positive Lyapunov exponent and an indication of deterministic chaos, are determined for both the observed series and the surrogates. The example of recurrence plots of several observations and their surrogate is given in Fig. A.3.
We analyse the sample of BHXBs to determine the nature of their Xray lightcurves. Our goal is to find hints of deterministic chaotic behaviour in the accreting black hole systems, as observed by the advanced Xray satellites that have good temporal resolution. Thus we verify the results of the nonlinear behaviour of the GRS 1915+105 (Misra et al. 2006) using the RXTE data, and we pursue our analysis with other sources. In particular, we study IGR J170913624 in four of its variability states. In addition, we also study the sample of several less certain candidates for the globally unstable accretion disk evolution, and ultimately we try to answer whether the two microquasars, mentioned above, are unique in their nature. The results of our analysis should allow for deeper insight into the accretion problem and to distinguish whether the variability is governed mainly by the set of nonlinear differential equations that describe the global disk structure and hence the fundamental physics, or possibly by local changes in the flow, and hence environmental conditions that are different for each source.
The article is organised as follows. In Sect. 2, we present our sample of black hole Xray binary sources and we briefly describe the RXTE dataextraction procedure. In Sect. 3, we briefly describe the recurrence analysis method and its general properties, while details of our method are given in Appendix A. We summarise the main results for the microquasar IGR J170913624 in Sect. 4 and in Appendix B we show in details how the method was applied to this source. In Appendix C, we give some more details about the significance of the method, which depends on the strength of the noise that is always present in the observed data sets. In Sect. 5 we continue to apply the recurrenceanalysis method to the rest of our sample of the black hole Xray binaries. Finally, in Sect. 6, we discuss our results and we present conclusions.
2. Observations
We analysed several tens of observations of six black hole Xray binaries, in which radiation pressure instability is considered to develop (Janiuk & Czerny 2011; Janiuk et al. 2015). We expect these lightcurves of these sources to exhibit signatures of nonlinear behaviour. Below we describe briefly the sources that belong to our sample and the method that we used for the extraction of the RXTE archival data.
2.1. Our sample
Observations of IGR J170913624.
The studied sources are the microquasars IGR J170913624, GRS 1915+105, GRO 165540, XTE J1650500, XTE J1550564, and GX 3394. Apart from the first one, which went into outburst toward the end of 2011, we selected them from the list given in (Janiuk & Czerny 2011) according to their presumably high accretion rate in the Eddington units and also their high count rate. We also included the source GX 3394 because a higher accretion rate was reported by Zdziarski et al. (2004) when combined with the results of Hynes et al. (2004) than was assumed by Janiuk & Czerny (2011), and because this source shows also other interesting features such as low frequency QPOs with changing frequency.
The first two systems are known to exhibit flarelike events in some states (see the references afterwards) and their similar nonlinear behaviour has been studied by different methods (Misra et al. 2004, 2006; Altamirano et al. 2011a; Pahari et al. 2013a,b). We decided to built our method on the sample of IGR J170913624, whose spectral states have been studied by Pahari et al. (2014). The light curve of IGR J170913624 during its heartbeat state has been modelled by Janiuk et al. (2015) with the model of the disk undergoing oscillations that are induced by radiation pressure instability. The wind ejected from the accretion disk regulates the amplitude of the oscillations and increase its frequency so that the simulations are in agreement with the observations. Moreover, the stronger winds, which are detected in some of the nonheartbeat states via the Xray spectroscopy, suppress the oscillations completely. Therefore, the idea of nonlinear instability producing the flares in IGR J170913624 is strongly supported by this modeling. Furthermore, we confirm the results of Misra et al. (2004) for some of the observations of GRS 1915+105 by our method.
Later we applied our method to several observations of the other four sources. We did not carry out an extended study of these sources, but aimed at finding several examples of observations that do or do not show traces of nonlinear behaviour. We took into account the observations with a high enough count rate and we focused on the rising and declining phases of the outbursts. At the time, the given source goes through the spectral state transition and this can be correlated with an unstable stage of the accretion disk caused by the internal instability.
Even though our sample of these four objects is limited and serves as a brief example rather than an indepth survey, we found indications of nonlinearity in several cases, apart from the two aforementioned microquasars.
Fig. 1 Lightcurve of the microquasar IGR J170913624, for the observation ID 96420010602, extracted in the energy range 2−10 keV with time bin dt = 0.5 s. We computed the significance of the nonlinear dynamics detected in this observation . 

Open with DEXTER 
2.1.1. IGR J170913624
The microquasar IGR J170913624 is a moderately bright transient Xray binary (peak flux level at ~20 mCrab in the range 20−100 keV) discovered in 2003 (Kuulkers et al. 2003). The 2011 outburst was the brightest one ever observed from IGR J17091, and the source flux increased up to 120 mCrab in the range 2−10 keV (Capitanio et al. 2012). The RXTE/PCA data showed quasiperiodic flarelike events occurring at a rate of between 25 and 30 mHz (Altamirano et al. 2011b), which resembled the “heartbeat” variation observed previously in the BH binary GRS 1915+105.
We tested our method on the sample of observations of IGR J170913624 that were selected from those analysed in Pahari et al. (2014) and which, according to their classification, belong to different spectral classes – hard intermediate state (HIMS), soft intermediate state (SIMS), intermediate variable state (IVS), and variable state (also denoted as ρ state). All the observations belong to proposal number 96 420, target 01 and were taken by the RTXE satellite between March and April 2011. The studied observations of IGR J170913624 are summarised in Table 1 and an exemplary lightcurve of the source is shown in Fig. 1.
2.1.2. GRS 1915+105
GRS 1915+105 is a very well known low mass Xray binary, discovered in 1992 (CastroTirado et al. 1992). Its behaviour was interpreted using the timedependent evolution of an accretion disk, which is thermally and viscously unstable (Taam et al. 1997; Janiuk et al. 2000). Belloni et al. (2000) classified 14 classes of variability for GRS 1915+105, out of which some are socalled “heartbeat states” and exhibit regular amplitude periodic oscillations, like the one shown in the lightcurve in Fig. 2. The flux emission during the heartbeat states of the source reaches the Eddington luminosity (Done et al. 2004). The average mass accretion rate of approximately 10^{18}−3 × 10^{19}g s^{1} (0.03−0.63 in Eddington units) in GRS 1915+105 was estimated by Neilsen et al. (2011).
2.1.3. GX 3394
Xray nova GX 3394, discovered in 1972 (Lavrov et al. 1997), exhibits high level of Xray activity. The distance was estimated as ≳7 kpc by Zdziarski et al. (2004), who also report that the luminosity varied during its outbursts (15 outbursts from 1987 to 2004) between 0.01−0.25L_{Edd}(d/ 8 kpc)^{2}(10 M_{⊙}/M). Similar values were also confirmed by Kolehmainen et al. (2011). However, Hynes et al. (2004) proposed the possibility of the location on the far side of the Galaxy with a distance d ≳ 15 kpc, which would imply that the source reaches Eddington luminosity in the peaks.
Fig. 2 Lightcurve of the microquasar GRS 1915+105, for the observation ID 20402010300, extracted in the full energy range 2−60 keV with time bin dt = 0.5s from Standard 1 data in the ρ state (classification according to Misra et al. 2004). 

Open with DEXTER 
2.1.4. GRO J165540
Low mass XRay binary GRO J165540 was discovered on July 27, 1994 with BATSE during the outburst and soon after its optical eclipsing counterpart was also found. Orosz & Bailyn (1997) found the mass of the black hole M = 7.02 ± 0.22 M_{⊙} and estimated the averaged mass transfer rate to be 3.4 × 10^{9} M_{⊙} yr^{1} = 2.16 × 10^{17} g s^{1}. A strong highly ionized accretion disk wind is reported by Miller et al. (2008) on April 1 during the 2005 outburst. The origin of the wind is likely not driven by the thermal process.
2.1.5. XTE J1550564
XTE J1550564 was first discovered by the ASM on board RXTE on 1998 September 7. During the 1998 outburst, the source exhibited a low and intermediate frequency QPO with evolving frequency in the range 0.08−13 Hz (Chakrabarti et al. 2009) and the accretion rate was a significant fraction of Eddington luminosity (Dunn et al. 2010; Heil et al. 2011).
2.1.6. XTE J1650500
XTE J1650500 is the black hole XRay binary with black hole mass 4−7.3 M_{⊙} (Orosz et al. 2004; Slaný & Stuchlík 2008). The object emits flares with nonthermal energy spectra that occur when the persistent luminosity is near 3 × 10^{34} erg s^{1}(d/4 kpc)^{2} in 2.8−20 keV range (Tomsick et al. 2003). Homan et al. (2006) estimate the distance as 2.6 ± 0.7 kpc based on the luminosity during the spectral state transition, assuming the black hole mass M = 4 M_{⊙}, for which the Eddington luminosity L_{Edd} ~ 1.1 × 10^{37} erg s^{1} in 0.5 keV−10 MeV.
2.2. Data analysis
We analysed the online data from the RXTE satellite. We extracted the time series using Heasoft 6.16 high energy astrophysics software package. Because we want to compare our results for IGR J17091 with the classification of Pahari et al. (2014), we use the filtering mentioned in that paper , taking data from all xenon layers from the PCU counter number 2 and we processed the events from channels 5 to 24 (2−10 keV approximately). Afterwards we selected several observations for the other five sources, which were observed by RXTE between 1996 and 2010. These observations were made in different observational modes with different data available. For each source we selected data in a lower energy band or in full energy range and normalized it to count/s/PCU. These sources and details about the data extraction are summarised in Table 3 in Sect. 5, together with the most important results of our theoretical analysis.
We adjust the proper binning time to minimize the error and at the same time to not lose the information about oscillations at the scale of several seconds. We use a minimal binning time of between 0.032 s and 0.5 s. The exact value depends on the flux and on the timescale of the variability in the data. Therefore, for the flux on the order of hundred counts per second we use the minimal binning size 0.5 s, if the average flux per PCU exceeds 500 cts/s. For 1000 cts/s, we can drop the value to 0.2 s, 0.1 s, or even 0.032 s. However, the variability is slower in some cases and we can increase the time bin above the minimal level.
While extracting the data, we use the xdf tool, filter the data with xtefilt, and we create the “good time interval”. For the science array format data files (standard1 mode or generic single bit mode) we extract the data using saextrct command, with a proper bin size. The standard1 mode data are collected in the full energy spectrum of PCA (2−60 keV), while the single bit mode contains one channel band in the lower energy range (it usually contains channels 0−35). For the science event format data files (generic event data mode), we extracted the lightcurves with the seextrct command with the energy range corresponding approximately to 2−10 keV. The resulting count rate is normalized to the number of PCUs that were used for the data extraction. The data mode used for the extraction is indicated in Table 3 in Sect. 5.
The background subtraction is ignored in our analysis. This is because it can be safely neglected for the bright sources. This is also because the background is not measured in the RXTE observations, but is simulated numerically. We do not, therefore, have a background measurement in the function of time, while the essential part of our analysis is to catch up the changes of about 0.5 s timescales.
We have not done the spectral analysis for the observations. The spectra were studied and in some cases classified by other authors. In Table 3, we provide their results or our guess for the spectral state that is based on the published plots (e.g. on the hardness ratio plot from Miller et al. 2008, for the 2005 outburst of GRO J165540).
3. Analysis of the chaotic processes
Our aim is to reveal important information about the BHXBs on the basis of their Xray light curves using the recurrence analysis, which is a wellestablished method for studying the properties of dynamical systems based on the behaviour of their phase space trajectories. More specifically, it looks at the “recurrences” in the phase space, which are the moments when the phase space trajectory returns close to itself, so that the points x(t_{i}) and x(t_{j}) are close. The long diagonal lines (of length l) are formed in the recurrence plot that represents the recurrence matrix (see relation A.2), when the pairs of successive points on the trajectory are close (∥ x(t_{i})−x(t_{j}) ∥ <ϵ, ∥ x(t_{i + 1})−x(t_{j + 1}) ∥ <ϵ,..., ∥ x(t_{i + l−1})−x(t_{j + l−1}) ∥ <ϵ). These lines correspond to the case, where the later part of the trajectory evolves inside a tube of ϵ diameter along the previous piece.
In the case of a periodic orbit, the recurrence plot would consist of infinite diagonals spaced by the period T. Stochastic processes yield randomly distributed points, which do not tend to form any lines or structures. A chaotic process produces lines with finite length, which is related to the maximal Lyapunov exponent, because for chaotic systems nearby trajectories eventually diverge outside the ϵtube. Therefore, one of the basic but important quantifiers of the recurrence analysis is the length of the longest diagonal line L_{max}. We quantify it simply as the number of points that make up the line.
Recurrence analysis can be supplied with the time delay reconstruction of the phase space trajectory (Takens 1981), so that the measured time series of one dynamical quantity could be used to study the underlying dynamical system. From this point of view, the method is suitable for the analysis of the astrophysical data obtained by, for example, Xray satellites, for which only the lightcurve is available. Therefore, we are also motivated to explore the abilities of recurrence analysis in the context of variability of Xray binaries.
Here, the underlying dynamical system being studied is the accreting gas governed by a global physical parameter, e.g. by the massaccretion rate. If the accretion is stationary, the produced flux has a constant mean and the variability is due to stochastic fluctuations. But if the global parameters are such that the intrinsic instability of an accretion disk develops, the massaccretion rate locally varies in time, which leads to deterministic changes of the flux, possibly also overlaid with some stochastic fluctuations.
A detailed description of the main quantities of the recurrence analysis used in this paper is given in Appendix A, where we show several examples of the technique. For more information, please refer to the extensive study of the method given in Marwan et al. (2007) and references therein. The details of the method, definitions of the terms, and quantities, together with a discussion about the possible drawbacks and difficulties arising from the presence of noise are given in Appendices B and C.
3.1. General approach
The recurrence analysis technique needs a careful setting of parameters. For noisy short data series, this kind of estimation is always complicated. Therefore instead of a direct estimate of these parameters, we have incorporated the method of surrogate data, as proposed by Theiler et al. (1992).
We first pose a null hypothesis about the measured time series, e.g. that the data are the product of temporally independent white noise or linearly autocorrelated Gaussian noise. Then we produce a set of surrogate series, which fulfill this hypothesis, but contain as much stochasticity as is possible (at given levels of the hypothesis), hence there is no further hidden nonlinear dynamics.
In the case of the temporally independent white noise, the surrogates are created as the random permutation of the values of the time series, thus such surrogates (we call them “shuffled surrogates”) have the same mean and variance as the original time series. In the case of the linearly autocorrelated Gaussian noise, the surrogate series are constructed this way so that they have the same distribution of values (flux per second per PCU for the observed lightcurves) and almost the same spectrum (Theiler et al. 1992). This can be achieved by an iterative algorithm called iterative amplitude adjusted Fourier transform algorithm (IAAFT) or its stochastic version SIAAFT, which is described by, e.g. Schreiber & Schmitz (2000); Venema et al. (2006). The algorithm iteratively replaces the amplitudes and magnitudes of the Fourier coefficients of the surrogate series (starting from the white noise) by corresponding values obtained from the original time series, so that at the end the surrogate series is a permutation of the original time series with almost the same spectrum. Throughout this paper, we refer to data series created by this procedure as simply “surrogates”, or IAAFT surrogates. For one experimental data series we typically construct one hundred surrogates. Afterwards, we apply the recurrence analysis both to the real experimental lightcurves and to the corresponding surrogates and we compare the obtained properties for the real and artificial data.
If the quantity measured for the real data differs significantly from the value obtained for its surrogate, we can reject our initial null hypothesis. This means that the nonlinear behaviour and perhaps chaotic nature of the system appears in the case of IAAFT surrogates, or nonstochastic but linear behaviour in the case of shuffled surrogates. Luckily, there is actually no need for the discriminating criterion to be a particular physical quantity, hence the majority of problems with accurate choice of the parameters for the analysis are overcome by our approach. The discriminating criteria are provided by the recurrence analysis of the time series.
3.2. Numerical method
For some tasks, including the creation of the IAAFT surrogates, we employ several procedures from the publicly available software package TISEAN^{1} (Hegger et al. 1999; Schreiber & Schmitz 2000). The shuffled surrogates are produced by our own code written in IDL. The recurrence analysis is performed on the basis of the software package described by Marwan et al. (2007); Marwan (2013), which also yields the cumulative histogram of diagonal lines. The linear regression for K_{2} estimation and other postprocessing is done by our IDL codes.
Below, we summarise the main steps to handle the data:

For every observation we extract the lightcurve (discussion aboutthe time and energy bin is given in Sect. 2.2 and inTable 3).
We rescale the extracted light curve to have zero mean and unit variance before producing the surrogates to simplify the comparison between different observations using the procedure rescale from TISEAN.
We use the procedure mutual from TISEAN to make an educated guess of the time delay Δt = kdt (see Appendix A and Fig. A.1).
We use the procedure false_nearest from TISEAN to make an educated guess of the embedding dimension m (see Appendix A and Fig. A.2).
We produce 100 surrogates (shuffled ones are produced by our own code, IAAFT ones are produced by the procedure surrogates from TISEAN) for each observation.
We find the recurrence threshold ϵ_{min} and ϵ_{max} such that the recurrence rate (RR) of the observation is ~1% and 25%, respectively^{2}. We create a sequence of ϵ in the range (ϵ_{min},ϵ_{max}) with constant difference. The number of used thresholds ranges from between 10 to 40, depending on the length of the data series (and thus on CPU time required for the analysis).
We construct the file with RQA quantifiers and the cumulative histograms of diagonal lines, using the program rp described by Marwan et al. (2007); Marwan (2013) with the found parameters k,m for all ϵ’s for the observation and the surrogates.
For each ϵ we compute the estimate of the Rényi’s entropy for the observation and for the set of surrogates. The Rényi’s entropy of the order α is, in general, defined as follows: (1)where p_{i} is the probability that the random variable has a value of i, and basically describes the randomness of the system. For a detailed explanation of the meaning of entropy K_{2}, see relations A.4 and A.5 in Appendix A.
We compute the average significance with respect to the shuffled surrogates and the averaged significance with respect to the IAAFT surrogates. The definition of significance of chaotic process is taken as (2)where N_{sl} is the number of surrogates, which have only short diagonal lines, and N^{surr} is the total number of surrogates, Q^{obs} and Q^{surr} are the natural logarithms of K_{2} entropy for the observed and surrogate data, respectively, and is the significance computed only from the surrogates, which have enough long lines according to the relation (3)The significance defined in this way expresses how much the value differs from the mean value measured in the units of the standard deviation of the set in the logarithmic scale σ_{Qsurr(ϵ)}. For further details, see Appendix B and the discussion before Eq. (B.2).
We distinguish between fully stochastic behaviour ( those observations that show no significant difference compared to the shuffled surrogates, or those that do not have enough long lines in the RP), nonstochastic, but linear behaviour (observations with a significant result with shuffled surrogates, but nonsignificant result with IAAFT surrogates, ) and nonlinear, possibly chaotic behaviour (observations with significant result with respect to the IAAFT surrogates, ).
In Tables 2 and 3, we provide the number of points of the observed time series, which is used for the analysis, the length of the longest line for several recurrence rates (the longer the diagonal lines are compared to the time delay k, the more regularly the system behaves), the averaged significance and/or and the number of thresholds for which the average is computed.
Results of our analysis for IGR J170913624.
4. Case study of IGR J170913624
We used the sample of observation of IGR J170913624 as described in Sect. 2.1.1 and in Table 1 to test the capabilities of recurrence analysis and to reveal important information about the dynamical properties of the source. Details of the procedure are given in Appendices A and B, main results are described here and summarised in Table 2.
In our study, we include several observations from the hard intermediate state (HIMS), soft intermediate state (SIMS), intermediate variable state (IVS), and variable state (also denoted as ρ state) of the source, see Pahari et al. (2014). All observations classified by Pahari et al. (2014) as SIMS and HIMS and the observation 0301, classified as IVS, have low and almost constant values of mutual information (see Fig. A.1). Three observations (0401, 0403 (IVS) and 0603 (ρ)) show a smooth decrease but no second maximum. All other ρ observations show oscillations with periods ranging from 10 s to 20 s.
We found that observations with no oscillations of mutual information show no lines long enough for computation of K_{2} (the criteria for the length and number of lines are described in Appendix B) and their L_{max} behaves similarly to those of both types of surrogates. Therefore, these observations are consistent with the hypothesis of temporally independent white noise (see Fig. A.6) and the source did not show any dynamical evolution during that time.
Table of RXTE observations of other sources.
The observations 0603 and 0501 do contain enough long lines, so that the significance can be computed. In both cases, nonsignificant results are obtained with respect to the IAAFT surrogates. However, we obtain significance with respect to shuffled surrogates for both, computing the average using N_{ϵ} = 13 and N_{ϵ} = 15 different thresholds for 0603 and 0501, respectively. Thus, we classify these two observations as a nonstochastic linear process. As mentioned in Appendix C, this can also mean, that these observations refer to the period, when the source was following a regular or “sticky trajectory” of the generally nonlinear dynamical system. This kind of orbit evolves for very long time very close to a regular island, sharing its main dynamical properties, but eventually departs to another part of the phase space, differing significantly from the regular motion (for more details, see e.g. Semerák & Suková 2012). Pahari et al. (2014) classified these observations as the variable ρ state.
The rest of the ρ state observations show significant results with . The values range between 1.63 for ObsID 0800 to 3.72 for ObsID 0602. We claim for these observations that the source was showing nonlinear behaviour.
By way of illustration, we also provide information about the length of the longest diagonal line for different values of RR in Table 2. Higher ratio L_{max}/k means that the system evolves similarly for a longer time compared to the period of the oscillations (i.e. that the phase trajectory in two different times evolves within ϵtube). Further details are given in Appendix B.
5. Using the method on other sources
The last step of our study is to apply our method to other sources listed in Sect. 2. We have chosen several observations of five objects, GRS 1915+105, GRO J165540, XTE J1650500, XTE J1550564, and GX 3394, which are summarised in Table 3.
At first, we examined observations of GRS 1915+105, a well studied source. Several papers about its lightcurves with respect to the nonlinear nature of the source have been published. We picked the observations studied by Misra et al. (2004), who classified them into three different groups according to their correlation dimension, namely those that show the chaotic behaviour (spectral classes β, λ, κ, and μ), nonstochastic behaviour (spectral classes θ, ρ, ν, α, and δ), and stochastic behaviour (spectral classes Φ, γ and χ). An exemplary lightcurve taken in the class ρ is shown in Fig. 2.
In agreement with their results, we found for all observations belonging in the first and second group and we did not obtain diagonal lines long enough for the K_{2} estimation for observation from class Φ. The averaged significance ranges between 1.6 and 5.1 for different observations. The variability of GRS 1915+105 is slower than in the case of IGR J170913624, so that the time delay Δt = kdt is longer and ranges between 10 s to as much as 125 s. Therefore, the ratio of the number of points of the data set and the embedding delay N/k is lower than for IGR J170913624 for a similar length of observation (a lower number of “cycles” is seen), hence lower significance could be expected for the same type of behaviour. The relatively high values of averaged significance thus provides quite strong evidence for the nonlinear behaviour of this source.
Fig. 3 Lightcurve of the source GX 3394. The observation ID is 95409011605 and the source was in its softintermediate state. Data were extracted in the energy range 2−10 keV (channels 5−25). We computed the significance of the nonlinear dynamics detected in this observation . 

Open with DEXTER 
Fig. 4 Lightcurve of the source XTE J1550564. The observation ID is 30188060300 and the data were extracted in the full energy range from standard1 mode data. We computed the significance of the nonlinear dynamics detected in this observation . 

Open with DEXTER 
For the black hole candidate, GX 3394, we analysed four observations from four different spectral classes classified according to Nandi et al. (2012). In Fig. 3, we show the exemplary lightcurve, as taken on September 24, 2010, when the source was in its softintermediate state. Here we found the averaged significance equal to 2.90. For the observation taken on March 26, 2010, the averaged significance is equal to 2.19. The third observation yields , which is close to our chosen threshold for the nonlinear behaviour. We have also computed the significance with respect to the shuffled surrogates , which show nonstochastic behaviour in all three cases. For the last observation of GX 3394, a nonsignificant result is obtained with respect to both types of surrogates, hence the variability of the source in this state is of stochastic origin.
Figure 4 shows the lightcurve of the source XTE J1550564, as taken on September 08, 1998. The source has been previously studied by Sobczak et al. (2000), and the lightcurve of observation 30191011400 shown therein looks similar to the variable state of IGR J170913624, but on much smaller timescales. Thus, we needed to extract this lightcurve with a very small time bin (dt = 0.032 s), which increases the influence of noise. However, this data set covers a much higher number of cycles (it has the highest N/k ratio in our sample). Our analysis of this observation yields the averaged significance , which is the highest among the sample. By way of comparison, we also provide the significance with respect to the shuffled surrogates, which is . In addition, three other observations yield quite high significance around 2.7. For the lightcurve presented in Fig. 4, the computed significance of the nonlinear dynamics is .
The other two observations of XTE J1550564 show nonsignificant results with IAAFT surrogates ( and ). However, their RPs contain quite long lines compared to their shuffled surrogates, which do not show any lines that are longer than Δt with the same recurrence parameters ( for all ϵ), so that in both cases. This indicates the nonstochastic behaviour of the source during these states.
Figure 5 shows an exemplary lightcurve of the source GRO J165540, as taken by RXTE on August 01, 1996 (10255010400B). For this data set, the significance of the nonlinear dynamical process that governs its variability is above . However, in this case the high value of significance could partially be caused by the fact that, during this observation, several episodes with different behaviours occur, so that the evolution is not stationary. In the case of nonstationary data, the theoretical problem is much more complicated and, thus, using the IAAFT surrogates that are created with the same spectrum may not be appropriate (Theiler et al. 1992).
In summary, for this source we found six observations with and seven others with nonsignificant results. For the latter group, we computed also the shuffled significance. The obsID 90058160700A again shows a nonsignificant result, hence it is of stochastic origin. Three other observations (90058160300A, 90058160300B, and 90058160200) yield moderate results , indicating a weakly nonstochastic origin. The shuffled surrogates of the last three observations (10255010400C, 20402021400A, and 20402021400B) have no long lines, leading to , which is strong evidence of nonstochastic behaviour of the source.
Fig. 5 Lightcurve of the source GRO J165540. The observation ID is 10255010400 and the data were extracted in the energy range 2−13 keV (channels 5−35). We computed the significance of the nonlinear dynamics detected in this observation . 

Open with DEXTER 
Finally, within the studied sample of observations of XTE J1650500, no significant traces of nonlinear behaviour were found. Figure 6 shows one of its exemplary lightcurves, from the observation ID 60113010802, taken on September 24, 2001, for which the computed significance of chaotic process was the highest. In any case, it was always below . For this source the significance with respect to the shuffled surrogates is also always below . Therefore we did not find any observations of this source, which departs from stochastic behaviour.
6. Discussion and conclusions
6.1. Discussion
In this article, we have used the recurrence analysis method to study the nonlinear behaviour of several Xray sources. They are the black hole candidates, in which their Xray luminosity originates in an accretion disk. The latter may be responsible for a shorttimescale, quasiregular oscillations of the emitted energy flux, through some physical process related to its hydrodynamical evolution. The possible trigger for nonlinear behaviour could be the moment when the disk enters the parameter region, where the thermalviscous instabilities occur. These instabilities lead to global limitcycle oscillations of the disk configuration (including density, mass accretion rate, temperature, etc.). This type of unstable evolution of the disk translates into nonlinear patterns in the light curve.
Fig. 6 Lightcurve of the source XTE J1650500. The observation ID is 60113011802 and the data were extracted in the energy range 2−10 keV (channels 5−24). We computed the significance of the nonlinear dynamics detected in this observation . 

Open with DEXTER 
Two microquasars, GRS 1915+105, and the recently discovered IGR J170913624, are already firm candidates for the limitcycle oscillations that are also called also the “heartbeat” state (see e.g. Neilsen et al. 2011). Our present analysis confirms that the variability in these sources is significantly governed by the nonlinear dynamics of accretion process. Other sources of that kind, which accrete presumably at high accretion rates and exhibit soft, diskdominated, states, have not yet been studied extensively enough to tackle this problem quantitatively. The question however remains, whether GRS 1915+105, and later its analogue, IGR J170913624, are the only two sources of this kind, or whether this kind of behaviour is common to all accreting black holes, if the physical conditions that trigger the nonlinearity are also found there. For Xray binaries other than these two microquasars, the situation is not so obvious. However, Janiuk & Czerny (2011) discussed the theoretical aspects of the thermalviscous instabilities in the accretion disks and present observational constraints for a sample of Galactic black hole systems. In particular, for the radiation pressure instability, a few other sources were hinted at, based on their estimated Eddington ratios and observed properties. It should be emphasised that the accretion rate reaching exactly the Eddington rate is not required by the theory for the radiation pressure instability to develop. Therefore, we aim to search for the nonlinear hydrodynamics of the accretion disk also hidden in the data for these other sources.
The black hole Xray transient source XTE J1650500 displayed its only outburst in the year 2001, and there were not many observations available for that source. Homan et al. (2003) detect the highfrequency variability in XTE J1650500 and interpret it as the orbital frequency at the innermost stable orbit around a Schwarzschild black hole, while Kalemci et al. (2003) suggest that the accretion flow geometry is that of the diskjetcorona system. No firm estimation exists in the literature about the value of the accretion rate in this source. Tomsick et al. (2003) detect in this source the Xray flares that have durations of between 62 and 215 s and peak fluxes that are 5−24 times higher than the persistent flux. Consequently, Janiuk & Czerny (2011) tentatively hints at this Xray binary as a prospective candidate for disk instabilityinduced oscillations, because of these timescales. However, these flares have nonthermal energy spectra, and possibly do not originate in the accretion disk/corona at all. They also occur at luminosities about a thousand times lower than the peak outburst luminosity, probably at the level about 2 × 10^{5}L_{Edd}. Therefore they are not good candidates for radiation pressure instability. Tomsick et al. (2003) also report on aperiodic oscillation with a 14day time scale. They relate that these variations in Xray light curve with the ionization instability. However, it is not possible to study such a long variability using our present method because no continuous observation of the required length (at least one order of magnitude longer than the period) is available. Our current analysis is based on earlier observations, up to October 2001, with a higher count rate. This showed that, in this source, the significance of the underlying nonlinear dynamics is very low. We would therefore rather reject it as a candidate for the limitcycle type of oscillations. This source requires further monitoring to be able to make more firm conclusions about the nature of its variability.
For the other sources in our study, the significance of the underlying deterministic chaotic process, being the intrinsic reason for the accretion disk’s observed variability, is much higher. This result is in full agreement with the observational facts from 20 years of data collection by Xray satellites.
In XTE J1550564, at the peak of the outburst, the luminosity is close to L_{Edd} (Heil et al. 2011), and the QPO that were observed in this source during its 1998 outburst were found to be linked to the disk count rate (Sobczak et al. 2000). XTE J1550564 is classified as a microquasar on the basis of its largescaled moving jets, detected at Xray and radio (Corbel et al. 2002), similar to GX 3394 (Corbel & Fender 2002). We expect these two objects to have similar disk variability characteristics and, along with the two wellstudied microquasars, nonlinear dynamical processes should also occur in these sources. Our current analysis confirms these expectations.
Initially, the spectral state of GRO J165540 was not determined (i.e. the low, high, or very high state) for its first outburst in 1994, based on the BATSE observation (Crary et al. 1996). However, Remillard et al. (1999) after analysing data spanning the Xray outburst from 1996−1997, estimated the unabsorbed Xray luminosity to above ~0.2L_{Edd} during the peak. This luminosity and accretion rate is quite sufficient for the instability of the accretion disk, which develops within the innermost 60−80 Schwarzschild radii from the black hole and leads to moderate luminosity oscillations (Janiuk & Czerny 2011). Méndez et al. (1998) studied the evolution of GRO J165540 through the high, intermediate, and low states. This source at the beginning of its decay might have even shown signatures of a very high state, just like other black hole candidates. Sobczak et al. (1999a) presented the full spectral analysis of the RXTE data for the 1996−1997 outburst of GRO J165540 and showed that, during the high/soft state, its spectrum is dominated by the soft thermal emission from the accretion disk. Comparing the two abovementioned black hole binaries, Sobczak et al. (2000) find that in both sources the QPO frequency correlates with the disk flux, and hence with the rate of mass accretion through the disk, thereby confirming that the intrinsic mechanism of this variability is related to the accretion process. Most recently, Uttley & KleinWolt (2015) report on an unusual soft state of GRO J165540, observed by the RXTE during its 2005 outburst. Chandra Xray grating observations have revealed a high massoutflow accretion disk wind in this state, which can be related to the accretion disk (in)stability evolution (Janiuk et al. 2015). Indeed, Miller et al. (2008) show that during its 2005 outburst, GRO J165540 ejected massive winds, possibly driven by a magnetic process in the accretion disk.
The fact that we have found nonlinear dynamics hidden in the lightcurves of almost all the studied sources means that the evolution of the accretion disk, and possibly the corona, is important for outgoing radiation and that the accretion flows are influenced by a physical instability, at least at some specific states.
In Appendix C we summarise the results presented in (Suková & Janiuk 2015b), where we test the capabilities of our method on two numerical trajectories. We show how the method works with two simulated trajectories of a complicated nonlinear system (motion of the test particle in the field of a black hole surrounded by thin massive disk given by exact solution of Einstein equations). One trajectory is a regular orbit, while the other one is chaotic. The chaotic orbit shows higher significance of nonlinear dynamics. We studied the influence of noise on the output of the method. While the significance for regular orbit quickly drops below unity, the significance of chaotic motion remains high up to comparable variance of the noise and the data.
6.2. Conclusions

We applied recurrence analysis to the observations of six blackhole Xray binaries detected by RXTE satellite.

We developed a method to distinguish between stochastic, nonstochastic linear, and nonlinear processes, using the comparison of the quantification of the recurrence matrices for the real and surrogate data.

We tested our method on the sample of observations of the microquasar IGR J170913624, whose spectral states were provided by Pahari et al. (2014). Results with for the “heartbeat” variable ρ state were obtained, which indicates nonlinear dynamics. Hence, our results also corroborate the validity of simulations of radiation pressure instabilityinduced oscillation of the accretion disk in this source, as reported by Janiuk et al. (2015).

We examined several observations of five other microquasars. Besides the wellstudied binary GRS 1915+105, we also found significant traces of nonlinear dynamics in three sources (GX 3394, XTE J1550564, and GRO J165540).

The nonlinear behaviour of the lightcurve during some observations provides evidence that the temporal evolution of the accretion flow in the binaries is governed by a lowdimensional system of nonlinear equations and the dynamics are chaotic. Hence the variability in the lightcurve is driven by a dynamical system of accreting gas that is described by a few parameters (e.g. the accretion rate) rather than by stochastic variations coming from random processes (flares) in the disk. A possible explanation for this is that the accretion disk is in the state prone to the thermalviscous instability and is undergoing induced limitcycle oscillations.

As pointed out in Misra et al. (2004), knowledge of the chaotic nature of the system can enlighten the fundamentals of the physical processes that yield the observed radiation. The rough estimate of the Rényi’s entropy, which is related to the maximal Lyapunov exponent, can be compared with estimates for magnetohydrodynamic simulations of the accretion flow to further confirm or refute the accretion models.
Our analysis of course does not answer the question of how the appearance or absence of nonlinear chaotic variability is related to the evolution of a given Xray binary during its outburst and the spectral state transitions. This level of analysis is extremely complex and currently beyond the scope of this article, but definitely worth further study. For now, our tentative picture would be that the purely stochastic variability occurs when the source is in its hard state, at the very beginning of the outburst. In the disk dominated soft state, as well as in the intermediate states, nonlinear variability, which is due to the chaotic process in the accretion disk, may appear. However, even in this state, the intrinsic coronal emission, which should be quite stochastic, as well as possibly reprocession of the hard Xray flux in the disk, may significantly affect the observed flux, and hence our results. On the other hand, if the corona above the disk forms through the evaporation of the upper layers of an accretion disk, then its variability will be related to the disk variability, with a possible short delay time (Janiuk & Czerny 2005).
Acknowledgments
We thank Ranjeev Misra, Piotr Zycki, Fiamma Capitanio, and Bozena Czerny for helpful discussions. We also thank the anonymous referee, whose comments helped us to make the paper more organised and clearer for the readers. This work was supported in part by the grant DEC2012/05/E/ST9/03914 from the Polish National Science Center.
References
 Abramowicz, M., Blaes, O., Horák, J., Kluźniak, W., & Rebusco, P. 2006, Class. Quant. Grav., 23, 1689 [NASA ADS] [CrossRef] [Google Scholar]
 Altamirano, D., & Belloni, T. 2012, ApJ, 747, L4 [NASA ADS] [CrossRef] [Google Scholar]
 Altamirano, D., Belloni, T., Linares, M., et al. 2011a, ApJ, 742, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Altamirano, D., Linares, M., van der Klis, M., et al. 2011b, ATel, 3225, 1 [NASA ADS] [Google Scholar]
 Belloni, T. M., & Stella, L. 2014, Space Sci. Rev., 183, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Belloni, T., KleinWolt, M., Méndez, M., van der Klis, M., & van Paradijs, J. 2000, A&A, 355, 271 [NASA ADS] [Google Scholar]
 Cao, L. 1997, Physica D: Nonlinear Phenomena, 110, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Capitanio, F., Del Santo, M., Bozzo, E., et al. 2012, MNRAS, 422, 3130 [NASA ADS] [CrossRef] [Google Scholar]
 CastroTirado, A. J., Brandt, S., & Lund, N. 1992, IAU Circ., 5590, 2 [NASA ADS] [Google Scholar]
 Chakrabarti, S. K., & Manickam, S. G. 2000, ApJ, 531, L41 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Chakrabarti, S. K., Dutta, B. G., & Pal, P. S. 2009, MNRAS, 394, 1463 [NASA ADS] [CrossRef] [Google Scholar]
 Corbel, S., & Fender, R. P. 2002, ApJ, 573, L35 [NASA ADS] [CrossRef] [Google Scholar]
 Corbel, S., Fender, R. P., Tzioumis, A. K., et al. 2002, Science, 298, 196 [NASA ADS] [CrossRef] [Google Scholar]
 Crary, D. J., Kouveliotou, C., van Paradijs, J., et al. 1996, ApJ, 463, L79 [NASA ADS] [CrossRef] [Google Scholar]
 Czerny, B., Siemiginowska, A., Janiuk, A., NikielWroczyński, B., & Stawarz, Ł. 2009, ApJ, 698, 840 [NASA ADS] [CrossRef] [Google Scholar]
 Das, T. K. 2003, ApJ, 588, L89 [NASA ADS] [CrossRef] [Google Scholar]
 Done, C., Wardziński, G., & Gierliński, M. 2004, MNRAS, 349, 393 [NASA ADS] [CrossRef] [Google Scholar]
 Dunn, R. J. H., Fender, R. P., Körding, E. G., Belloni, T., & Cabanac, C. 2010, MNRAS, 403, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Eckmann, J.P., Kamphorst, S. O., & Ruelle, D. 1987, Europhys. Lett., 4, 973 [NASA ADS] [CrossRef] [Google Scholar]
 Grassberger, P. 1983, Phys. Lett. A, 97, 227 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Hegger, R., Kantz, H., & Schreiber, T. 1999, Chaos: An Interdisciplinary Journal of Nonlinear Science, 9, 413 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Heil, L. M., Vaughan, S., & Uttley, P. 2011, MNRAS, 411, L66 [NASA ADS] [CrossRef] [Google Scholar]
 Hirose, S., Blaes, O., & Krolik, J. H. 2009, ApJ, 704, 781 [NASA ADS] [CrossRef] [Google Scholar]
 Homan, J., KleinWolt, M., Rossi, S., et al. 2003, ApJ, 586, 1262 [NASA ADS] [CrossRef] [Google Scholar]
 Homan, J., Wijnands, R., Kong, A., et al. 2006, MNRAS, 366, 235 [NASA ADS] [Google Scholar]
 Hynes, R. I., Steeghs, D., Casares, J., Charles, P. A., & O’Brien, K. 2004, ApJ, 609, 317 [NASA ADS] [CrossRef] [Google Scholar]
 Janiuk, A., & Czerny, B. 2005, MNRAS, 356, 205 [NASA ADS] [CrossRef] [Google Scholar]
 Janiuk, A., & Czerny, B. 2011, MNRAS, 414, 2186 [NASA ADS] [CrossRef] [Google Scholar]
 Janiuk, A., Czerny, B., & Siemiginowska, A. 2000, ApJ, 542, L33 [NASA ADS] [CrossRef] [Google Scholar]
 Janiuk, A., Czerny, B., & Siemiginowska, A. 2002, ApJ, 576, 908 [NASA ADS] [CrossRef] [Google Scholar]
 Janiuk, A., Grzedzielski, M., Capitanio, F., & Bianchi, S. 2015, A&A, 574, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jiang, Y.F., Stone, J. M., & Davis, S. W. 2013, ApJ, 767, 148 [NASA ADS] [CrossRef] [Google Scholar]
 Kalemci, E., Tomsick, J. A., Rothschild, R. E., et al. 2003, ApJ, 586, 419 [NASA ADS] [CrossRef] [Google Scholar]
 Kennel, M. B., Brown, R., & Abarbanel, H. D. I. 1992, Phys. Rev. A, 45, 3403 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Kolehmainen, M., Done, C., & Díaz Trigo, M. 2011, MNRAS, 416, 311 [NASA ADS] [Google Scholar]
 Kopáček, O., Karas, V., Kovář, J., & Stuchlík, Z. 2010, ApJ, 722, 1240 [NASA ADS] [CrossRef] [Google Scholar]
 Krolik, J. H., Hirose, S., & Blaes, O. 2007, ApJ, 664, 1045 [NASA ADS] [CrossRef] [Google Scholar]
 KunertBajraszewska, M., Janiuk, A., Gawroński, M. P., & Siemiginowska, A. 2010, ApJ, 718, 1345 [NASA ADS] [CrossRef] [Google Scholar]
 Kuulkers, E., Lutovinov, A., Parmar, A., et al. 2003, ATel, 149, 1 [NASA ADS] [Google Scholar]
 Lasota, J.P. 2001, New A Rev., 45, 449 [NASA ADS] [CrossRef] [Google Scholar]
 Lavrov, S. V., Borozdin, K. N., Aleksandrovich, N. L., et al. 1997, Astron. Lett., 23, 433 [NASA ADS] [Google Scholar]
 Lightman, A. P., & Eardley, D. M. 1974, ApJ, 187, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Marwan, N. 2013, Commandline Recurrence Plots, http://tocsy.pikpotsdam.de/commandlinerp.php, accessed: 20150721 [Google Scholar]
 Marwan, N., Carmen Romano, M., Thiel, M., & Kurths, J. 2007, Phys. Rep., 438, 237 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Méndez, M., Belloni, T., & van der Klis, M. 1998, ApJ, 499, L187 [NASA ADS] [CrossRef] [Google Scholar]
 Miller, J. M., Raymond, J., Reynolds, C. S., et al. 2008, ApJ, 680, 1359 [NASA ADS] [CrossRef] [Google Scholar]
 Misra, R., Harikrishnan, K., Mukhopadhyay, B., Ambika, G., & Kembhavi, A. 2004, ApJ, 609, 313 [NASA ADS] [CrossRef] [Google Scholar]
 Misra, R., Harikrishnan, K., Ambika, G., & Kembhavi, A. 2006, ApJ, 643, 1114 [NASA ADS] [CrossRef] [Google Scholar]
 Nandi, A., Debnath, D., Mandal, S., & Chakrabarti, S. K. 2012, A&A, 542, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Neilsen, J., Remillard, R. A., & Lee, J. C. 2011, ApJ, 737, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Orosz, J. A., & Bailyn, C. D. 1997, ApJ, 477, 876 [NASA ADS] [CrossRef] [Google Scholar]
 Orosz, J. A., McClintock, J. E., Remillard, R. A., & Corbel, S. 2004, ApJ, 616, 376 [NASA ADS] [CrossRef] [Google Scholar]
 Ott, E. 2002, Chaos in Dynamical Systems (Cambridge University Press) [Google Scholar]
 Pahari, M., Misra, R., Mukherjee, A., Yadav, J. S., & Pandey, S. K. 2013a, MNRAS, 436, 2334 [NASA ADS] [CrossRef] [Google Scholar]
 Pahari, M., Yadav, J. S., Rodriguez, J., et al. 2013b, ApJ, 778, 46 [NASA ADS] [CrossRef] [Google Scholar]
 Pahari, M., Yadav, J. S., & Bhattacharyya, S. 2014, ApJ, 783, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Remillard, R. A., Morgan, E. H., McClintock, J. E., Bailyn, C. D., & Orosz, J. A. 1999, ApJ, 522, 397 [NASA ADS] [CrossRef] [Google Scholar]
 Schreiber, T., & Schmitz, A. 2000, Physica D: Nonlinear Phenomena, 142, 346 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Semerák, O., & Suková, P. 2012, MNRAS, 425, 2455 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Slaný, P., & Stuchlík, Z. 2008, A&A, 492, 319 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sobczak, G. J., McClintock, J. E., Remillard, R. A., Bailyn, C. D., & Orosz, J. A. 1999a, ApJ, 520, 776 [NASA ADS] [CrossRef] [Google Scholar]
 Sobczak, G. J., McClintock, J. E., Remillard, R. A., et al. 1999b, ApJ, 517, L121 [NASA ADS] [CrossRef] [Google Scholar]
 Sobczak, G. J., McClintock, J. E., Remillard, R. A., et al. 2000, ApJ, 531, 537 [NASA ADS] [CrossRef] [Google Scholar]
 Stella, L., & Vietri, M. 1998, ApJ, 492, L59 [NASA ADS] [CrossRef] [Google Scholar]
 Suková, P. 2011, J. Phys. Conf. Ser., 314, 012087 [NASA ADS] [CrossRef] [Google Scholar]
 Suková, P., & Janiuk, A. 2015a, MNRAS, 447, 1565 [NASA ADS] [CrossRef] [Google Scholar]
 Suková, P., & Janiuk, A. 2015b, in Workshop Series of the Argentinian Astronomical Society, submitted [Google Scholar]
 Suková, P., & Semerák, O. 2012, in AIP Conf. Ser. 1458, eds. J. Beltrán Jiménez, J. A. Ruiz Cembranos, A. Dobado, A. López Maroto, & A. De la Cruz Dombriz, 523 [Google Scholar]
 Taam, R. E., Chen, X., & Swank, J. H. 1997, ApJ, 485, L83 [NASA ADS] [CrossRef] [Google Scholar]
 Takens, F. 1981, Detecting strange attractors in turbulence (Springer) [Google Scholar]
 Theiler, J., Eubank, S., Longtin, A., Galdrikian, B., & Farmer, J. D. 1992, Phys. D: Nonlinear Phenomena, 58, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Thiel, M., Romano, M. C., Kurths, J., et al. 2002, Phys. D: Nonlinear Phenomena, 171, 138 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Thiel, M., Romano, M. C., Read, P. L., & Kurths, J. 2004, Chaos: An Interdisciplinary Journal of Nonlinear Science, 14, 234 [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
 Titarchuk, L., & Osherovich, V. 2000, ApJ, 542, L111 [NASA ADS] [CrossRef] [Google Scholar]
 Tomsick, J. A., Kalemci, E., Corbel, S., & Kaaret, P. 2003, ApJ, 592, 1100 [NASA ADS] [CrossRef] [Google Scholar]
 Uttley, P., & KleinWolt, M. 2015, MNRAS, 451, 475 [NASA ADS] [CrossRef] [Google Scholar]
 Venema, V., Ament, F., & Simmer, C. 2006, Nonlinear Processes in Geophysics, 13, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Zdziarski, A. A., Gierlinski, M., Mikolajewska, J., et al. 2004, MNRAS, 351, 791 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Recurrence analysis of the observed time series of IGR J170913624
Fig. A.1 The dependence of mutual information on the time delay Δt for the set of observations of IGR J170913624 listed in Table 1. Observations with low and constant mutual information are plotted with black dashed lines, observations with oscillating mutual information are plotted with grey solid lines and those with intermediate behaviour are plotted by yellow dotted lines. The vertical line indicates the chosen time delay for our analysis. 

Open with DEXTER 
The recurrence analysis is based on the study of the times, when the trajectory returns close to itself in the phase space (closer than a certain recurrence threshold ϵ). However, the trajectory has to be reconstructed from observed time series with the time/delay technique. Hence the resulting phase space vector is given as (A.1)where x(t) is the time series, Δt is the embedding delay and m is the embedding dimension. When working with measured time series, the embedding delay Δt has to be an integer multiple of the binning time of the data dt, hence Δt = k dt. The determination of these parameters is a subtle issue (Kennel et al. 1992; Cao 1997; Thiel et al. 2004; Marwan et al. 2007). The usual practise is to determine the time delay as the first distinct minimum of the time delayed mutual information, for which we use the procedure mutual from the TISEAN package. In Fig. A.1 the dependence of the mutual information on the time delay Δt is given for the set of observations of IGR J170913624. The observations show two kinds of behaviour. Mutual information of some lightcurves is very low even for small time delays and is roughly constant for all values of Δt (black dotted lines), whereas the course of mutual information of others is smooth and shows oscillations with period ranging from 10 s to 20 s (grey solid lines). Three observations show intermediate behaviour, there is a smooth decrease of the mutual information, but no clear second maximum, the values then seem to stabilise around a constant value of the first minimum (yellow dotted line). Because we want to compare results from the analysis with the same parameters for all observations, we choose the time delay Δt = 7 s (for dt = 0.5 s is k = 14) as a reasonable value for each lightcurve.
To find the appropriate embedding dimension, the wellknown method is to determine the number of false nearest neighbours (Kennel et al. 1992). This procedure is based on the fact that, when the embedding dimension is too low, it can happen that the points, which are not close on the attractor, are projected close together and become false neighbours. When the dimension is increased sufficiently, the ambiguity of their mutual position is removed and their distance grows substantially. Hence, we set a tolerance threshold R_{tol} and denote the point as a false neighbour, if the distance in m and m + 1 dimension increase by a factor bigger than R_{tol}. In this way we investigate the change of distance for each point of the trajectory and its nearest neighbour in the reconstructed phase space. If the dimension is big enough to cover the whole attractor, the fraction of false neighbours goes to zero. We use the procedure false_nearest from the TISEAN package. However, because our data are both short and noisy, we have to skip the second criterion, which omits points being farther than the standard deviation of the data divided by the threshold. Otherwise we would have very small number of points to consider and the results would be meaningless. Therefore we are forced to also include points that are not very close neighbours (even when they are the nearest ones) because there is not enough really close points, owing to the noise and the limited amount of data. We also have to be careful about the size of the threshold R_{tol}, which cannot be very high, otherwise almost no nearest neighbour distance can grow by this factor due to the finite size of the attractor. This can affect the results slightly.
Fig. A.2 The dependence of the ratio of false nearest neighbours on the embedding dimension m for the set of observations of IGR J170913624 for the choice R_{tol} = 10. 

Open with DEXTER 
Fig. A.3 RP for several observations of IGR J170913624 from Table 1. 1. Row: left − SIMS/HIMS state 0200 (left upper corner in blue) and ρ state 0602 (right lower corner in red). Long diagonal lines appear in red RP. Right − the zoom of preceding plot. The diagonal lines are affected by the noise in such a way that the apparently long lines are composed of shorter squareshaped pieces. Next rows: in the red right lower halves of each plot RP are plotted for observations 0401 (2. row left, IVS state), 0402 (2. row right, ρstate), 0502 (3. row left, SIMS state) and 0701 (3. row right, ρstate). In the left upper half of each plot one of the corresponding surrogates is shown in green. 

Open with DEXTER 
We choose the threshold R_{tol} = 10. The corresponding false nearest neighbour ratio for the set of IGR’s observations is given in Fig. A.2. We set the embedding dimension m = 10, in which for almost all observations the false nearest neighbour ratio drops to zero. However, even m = 6 could be a good choice, because for this value there is only 1% or less of false neighbours.
These two methods are used for guessing some appropriate values of the two parameters, but the results of the recurrence analysis do not strongly depend on it. This is supported by Thiel et al. (2004), who claim that some of the dynamical invariants estimated by recurrence analysis, including second order Rényi entropy K_{2}, do not depend on those parameters.
Having the values for embedding delay and dimension determined, we can now proceed with the recurrence analysis. The basic object of the analysis is the recurrence matrix, defined as follows: (A.2)where x_{i} = x(t_{i}) are (N) points of the reconstructed phase trajectory and Θ is the Heaviside step function. The matrix thus contains only 1s and 0s and can easily be visualised by plotting a dot at the coordinates i, j whenever R_{i,j}(ϵ) = 1, which is called the recurrence plot (RP). From the definition (A.2) it is obvious, that the RP is symmetrical with respect to the main diagonal, which is formed by dots (for i = j) and which is omitted in further computations.
The comparison of the RPs for two observations, 0200 (SIMS/HIMS – upper left corner on the left panel in first row) and 0602 (ρclass – lower right corner of the same panel) is given in Fig. A.3. Because the plot is symmetrical with respect to the main diagonal, we show only one half of the RP for each observation in the upper/bottom trilateral half of the plot. Different geometrical structures contained in the RP could be seen in these two examples. The structures can be quantified based on the statistical properties of the recurrence matrix.
Fig. A.4 The cumulative histogram p_{c}(ϵ,l) for the observation 0602 computed for ϵ = 2.8. The length of the lines l is given in number of points of the lines (thus the total time of the line is given as t_{l} = lΔt). The estimate of the Rényi’s entropy K_{2} = 0.0699 ± 0.0003 is obtained as the linear regression of the central part of the cumulative histogram (indicated by a solid red line, the dashed red lines indicate the used part of the histogram). 

Open with DEXTER 
The most important entity for the recurrence analysis is the existence of long diagonal lines, which indicate periodic (and regular) behaviour. Such structures are much more prominent in the RPs of ρ observations. The quantification of the visual information is contained in the histogram of diagonal lines of a certain prescribed length l, (A.3)Based on the histogram, more quantifiers can be defined (Marwan et al. 2007). Here we focus mainly on , which is the length of the longest diagonal line (except for the main diagonal), and on the estimate of the secondorder Rényi’s entropy K_{2}. This entropy, also called correlation entropy or correlation exponent, collision entropy or just Rényi’s entropy (Grassberger 1983), is defined for the phase trajectory x(t) as (A.4)where the phase space is coarsegrained into mdimensional hypercubes of size ϵ. Here, p_{i1,...,il} describes the probability that the first point of the trajectory x(t = Δt) is inside the box i_{1}, the second point x(t = 2Δt) is in the box i_{2} and so forth to the lth point being in the box i_{l}. It can be shown (Marwan et al. 2007) that the square of this probability is connected with the existence of a diagonal line in the RP (the diagonal line means that the trajectory goes twice through the same sequence of phasespace boxes).
Fig. A.5 Rényi’s entropy K_{2} versus recurrence rate for the observation 0602 for three different embedding dimensions indicated in the figure. 

Open with DEXTER 
More precisely, K_{2} is related to the cumulative histogram of diagonal lines p_{c}(ϵ,l), describing the probability of finding a line of minimal length l in the RP, by the relation (A.5)hence we can estimate the value of K_{2} as the slope of the logarithm of the cumulative histogram versus l for constant ϵ (see an example in Fig. A.4). This estimation holds for the limit l → ∞, particularly for l> Δt. The important property of K_{2} entropy is that it is the lower estimate of the sum of positive Lyapunov exponents of the system, hence it is positive for chaotic systems (Marwan et al. 2007).
Considering formula (A.5) for two different values of threshold ϵ_{2} = ϵ_{1} + Δϵ, we can also estimate the value of the correlation dimension D_{2} by (A.6)However, to obtain reasonable results for D_{2}, much more data points are needed than for estimating K_{2} (Marwan et al. 2007). Because we are limited by the length of individual observations, which are not longer than several thousand seconds, which yields twice the data points for our chosen binning time dt = 0.5s, we were unable to determine the value of D_{2} by means of the recurrence analysis.
Fig. A.6 Dependence of the longest diagonal line L_{max} given in points on the recurrence threshold ϵ for the observation 0502 (thick red line). The ensemble of 100 surrogates are shown in grey and the ensemble of shuffled surrogates are shown in cyan. 

Open with DEXTER 
Besides the short duration of the observation, another difficulty is the presence of noise, which affects the structures contained in RP. This is documented in the right panel in the first row of Fig. A.3, which shows the zoom of the left panel for t ∈ (50,200) s. Here the apparently long lines in fact consist of shorter and broader pieces separated by gaps. The gaps appear when the noise takes the “shouldbe the recurrence point” out from the ϵneighbourhood, and on the other hand, the edges on the lines appear when noise brings the “shouldnotbe the recurrence point” inside. Hence, there are more shorter diagonal lines and less longer ones, thus the slope of the histogram is steeper.
The usual way to deal with the noise in data is to take a big enough radius of the neighbourhood, so that it covers the noisy perturbations of the position, which according to Thiel et al. (2002) and Marwan et al. (2007) is ϵ = 5σ, where σ is the standard deviation of the noise. However, a part of the noise is caused by the measurement process and another part stems from the stochastic processes behind the emission of the radiation. Moreover, the amount of noise and its standard deviation is unknown, but it is quite large comparing to the standard deviation of the whole data set (see Fig. 1). Hence, this choice of threshold would not be appropriate because it would be huge and almost every point would become the recurrence point.
We consider these thresholds, which yield the recurrence rate (RR – ratio of the recurrence points to all points of the recurrence matrix) within 1% to 25% and we have to bear in mind, that because of the noise, the slope of the cumulative histogram and K_{2} is overestimated. An example of the decreasing dependence of the K_{2}estimate on the recurrence rate with Δt = 7 s, and three different embedding dimensions m_{1} = 10, m_{2} = 20, and m_{3} = 30, is given in Fig. A.5.
Appendix B: Comparison between the surrogate data and observations of IGR J17091
We first study the appearance of long diagonal lines in RP. For the real data set and its surrogates, we compare the dependence of L_{max} on chosen ϵ. We obtain three different types of output.
For observations from states other than ρ and for the used range of thresholds, we do not see any lines longer than Δt, hence it is impossible to estimate K_{2}. Moreover, in Fig. A.6, we can see that the length of the longest line for the data (thick red line) and its surrogates (grey lines) do not differ. We conclude that the origin of these observations is sufficiently well described by the assumption of linearly autocorrelated Gaussian noise. We can also try to soften our null hypothesis and use the shuffled surrogates. It turns out, that all the observations with short lines only are also perfectly consistent with this null hypothesis, thus there is no evidence of any dynamics at all (in Fig. A.6 the shuffled surrogates are plotted by cyan color).
The next group of observations are those that do contain longer diagonal lines (usually for higher ϵ), but the values of L_{max} for data and IAAFT surrogates do not differ notably. Within our sample of observations of IGR J170913624, only the observations 0501 and 0603 can be assigned to this subgroup, the example for 0501 can be seen in Fig. B.1. However, comparing the data with the shuffled surrogates (cyan lines in Fig. B.1), we find a striking difference, meaning that these observations cannot be described as temporally independent, identically distributed, random data. Hence, in this case, our analysis suggests that the data comes from a linear dynamical process (linearly autocorrelated Gaussian noise). The rest of ρ observations show stronger difference between the data and its surrogates, is higher than for all surrogates for a wide range of ϵ. This behaviour indicates nonlinear features of the background dynamics.
Fig. B.1 The same as in Fig. A.6 for the observation 0501. The condition l> Δt on lines’ length used for K_{2}estimation is indicated by a black horizontal line. 

Open with DEXTER 
Fig. B.2 The same as Fig. A.6, but for the observation 0602. 

Open with DEXTER 
To quantify these results, we compute the estimate of K_{2} for observations with enough long lines and we use this quantity as our discriminating quantity. The computation goes as follows.
Fig. B.3 Comparison of recurrence rate for observation 0602 (value is indicated as the horizontal thick red line) and its hundred surrogates (number of surrogate is on the xaxis) computed for the same recurrence parameters (ϵ = 2.8,m = 10,Δt = 14dt). 

Open with DEXTER 
For a chosen value of ϵ and the chosen embedding dimension and delay (m = 10 and Δt = 7 s), we compute the recurrence matrix and the cumulative histogram for the data and for its surrogates. Usually, the recurrence rate for the surrogates is lower, which means that the shuffled order of points yields more distant points after the reconstruction of the phase space trajectory (see Fig. B.3); however, for higher recurrence thresholds this difference is smeared. Afterwards, we estimate K_{2} from the slope of the cumulative histogram for the central part of the histogram which satisfies l> Δt, p_{c}(ϵ,l) >N_{min} and l< 0.1Ndt. We also require that there is at least five points of the cumulative histogram, which satisfy our criteria and could be used for the linear regression. The value of N_{min} has to be chosen with respect to the length of the time series and the total number of lines. We usually choose N_{min} = 100. The example of cumulative histogram computed for the observation 0602 with ϵ = 2.8 is given in Fig. A.4, where the used region of the histogram is marked by two vertical dashed lines.
We execute the same analysis with the same parameters on the set of surrogates and we obtain the set of values (see Fig. B.4). According to Theiler et al. (1992), we define a significance of the obtained result as a weighted difference between the estimate of Rényi’s entropy for the observed data and the mean of estimates of Rényi’s entropy for ensembles of the surrogates . However, as can be seen from relation A.5, the quantity K_{2} is nonnegative and behaves in an exponential manner (for 10times longer diagonal lines, K_{2} drops down by one order). Therefore we first introduce the quantity Q = ln(K_{2}), which scales linearly, and then define the significance stemming from the estimate of K_{2} as (B.1)where σ_{Qsurr} is the standard deviation of the set (see Fig. B.5).
Fig. B.4 Comparison of estimates of K2 for observation 0602 (indicated by horizontal thick red line) and its hundred surrogates (number of surrogate is on the xaxis) computed for the same recurrence parameters (ϵ = 3.25,m = 10,Δt = 14dt). 

Open with DEXTER 
When the observed data yields , but the RP of its surrogates does not contain enough long lines, it is impossible to compute Q^{surr}. Obviously this should contribute to the significance in the case, that is lower than . However, the values of Q^{surr} are not available. We set the significance of the case, when all N^{surr} surrogates have only short lines, to be .
Fig. B.5 Comparison of Q for observation 0602 (indicated by the horizontal thick red line) and its hundred surrogates (number of surrogate is on the xaxis) computed for the same recurrence parameters (ϵ = 3.25,m = 10,Δt = 14dt). The mean value and the standard deviation of the ensemble of surrogates is indicated by horizontal lines. 

Open with DEXTER 
On the other hand, if is higher than , then some of the surrogates have lines longer than the observed data. Here, the coincident existence of surrogates, which do not have long lines at all, lessens the significance. This fact is incorporated into our procedure in the way that we add or subtract the two significances and according to the sign of the difference .
Finally, we obtain the total significance as the weighted difference of significances of these two cases, namely (B.2)where N_{sl} is the number of surrogates, which have only short lines, and is the total number of surrogates. Our definition of significance yields negative values for some cases, e.g. when and N_{sl} = 0. Since we expect the real nonlinear dynamical system to contain more long lines than the artificial surrogates, only positive values of significance serve as the indication of nonlinear dynamics. In fact, we use the value as the limit, above which we consider the observation to be produced by a nonlinear dynamical system.
The last ambiguous step in the calculation of the significance is the choice of the recurrence threshold ϵ. There are two possible ways how to treat this issue. One way is to choose some recurrence rate for the observation (e.g. RR ~ 10%), choose the threshold appropriately and compute the significance. This approach is suitable for quick analysis, however it is not very accurate. In this paper we adopted the other way, which lies in finding the range of ϵ, for which RR ∈ (1%,25%), computing the significance for a linearly spaced set of ϵ in this interval (typically, we used N_{ϵ} ~ 10−40 values of ϵ dependent on the number of points of the data set) and establishing the final decisive value as the average of the obtained significances . If for some ϵ the observed data set does not have enough long lines for the K_{2} estimation, this value is not taken into account in the evaluation.
In Table 2 we summarised the length of the longest diagonal line (expressed as the number of points) for three different recurrence rates RR ~ 5%, 10%, and 20%, the resulting averaged significance and the number of values of ϵ used for the averaging. All ρ observations except of 0501 and 0603 (recognised earlier as a special group with respect to L_{max} behaviour) have .
In the same way we can also compute the significance for the other null hypothesis, using the set of shuffled surrogates. This quantity serves as the indicator of departure from stochastic behaviour. Typically, we compute this quantity only for those observations, which have enough long lines but show nonsignificant results with the IAAFT surrogates. This is because the data sets which have significant results for the nonlinear dynamics should automatically have significant results for nonstochastic nature. However, to confirm this assumption, we also provide several examples of this quantity for observations displaying the nonlinear features in Sect. 5.
Appendix C: Testing the method with simulated time series
In general, our method can be applied to different kinds of time series, which are the product of any dynamical system, not only to Xray lightcurves. For testing, it is desirable to apply the method on time series whose nature is known and which has been studied by other means. Because the dynamical system behind the observed lightcurves is a priori not known, we prefer to use numerical data. We chose the numerical time series, which describe the motion of geodesic test particles in the field of a static black hole surrounded by a massive thin disk.The details about the testing procedure will be given in separate paper (Suková & Janiuk 2015b), here we only summarise the main results.
The analysis we performed on the numerical time series revealed what the recurrence analysis provides for regular and chaotic trajectories of the complicated nonlinear system, whose
phase space consists of regular KAM tori nested in the chaotic layers. We saw that, even though governed by the same global set of equations, the evolution of the orbits in the regular and chaotic parts of the phase space can differ significantly, which is also reflected by the patterns in the recurrence plots. In agreement with the fact that regular orbits in conservative systems have zero Lyapunov exponents, the estimate of K_{2} is much lower for the regular trajectory than for its chaotic counterpart.
Despite the expectation, that regular motion should not exhibit nonlinear behaviour in our analysis, the significance for the regular orbit is quite high, however it is still almost twice as low as for the chaotic one. We argue that the reason for this is the way the surrogates data are constructed, which means that they have exactly the same value distribution but they reproduce the spectrum only approximately, depending also on the available length of the data set. In the case of the regular motion, there are very narrow peaks in the spectrum and the error in reproducing such spectrum causes the very long diagonal lines to be broken.
However, when the strength of the added noise is increased, the significance of the regular orbit drops down very quickly, while the chaotic orbit shows solid signs of nonlinear behaviour, even in the case of when the added noise has the same variance as the original data. Yet, even for high noise levels, L_{max} for the regular orbit is higher than or similar to the chaotic orbit. Hence, observations which show nonsignificant results with the IAAFT surrogates, but very significant results for shuffled surrogates, and contain long diagonal lines compared to the time delay Δt, can possibly be explained as regular motion. Another possibility is that such observations cover the part of the chaotic trajectory that is called ‘sticky orbit’. Because the observations only cover a limited time interval, it is not possible to make a distinction between these two cases.
All Tables
All Figures
Fig. 1 Lightcurve of the microquasar IGR J170913624, for the observation ID 96420010602, extracted in the energy range 2−10 keV with time bin dt = 0.5 s. We computed the significance of the nonlinear dynamics detected in this observation . 

Open with DEXTER  
In the text 
Fig. 2 Lightcurve of the microquasar GRS 1915+105, for the observation ID 20402010300, extracted in the full energy range 2−60 keV with time bin dt = 0.5s from Standard 1 data in the ρ state (classification according to Misra et al. 2004). 

Open with DEXTER  
In the text 
Fig. 3 Lightcurve of the source GX 3394. The observation ID is 95409011605 and the source was in its softintermediate state. Data were extracted in the energy range 2−10 keV (channels 5−25). We computed the significance of the nonlinear dynamics detected in this observation . 

Open with DEXTER  
In the text 
Fig. 4 Lightcurve of the source XTE J1550564. The observation ID is 30188060300 and the data were extracted in the full energy range from standard1 mode data. We computed the significance of the nonlinear dynamics detected in this observation . 

Open with DEXTER  
In the text 
Fig. 5 Lightcurve of the source GRO J165540. The observation ID is 10255010400 and the data were extracted in the energy range 2−13 keV (channels 5−35). We computed the significance of the nonlinear dynamics detected in this observation . 

Open with DEXTER  
In the text 
Fig. 6 Lightcurve of the source XTE J1650500. The observation ID is 60113011802 and the data were extracted in the energy range 2−10 keV (channels 5−24). We computed the significance of the nonlinear dynamics detected in this observation . 

Open with DEXTER  
In the text 
Fig. A.1 The dependence of mutual information on the time delay Δt for the set of observations of IGR J170913624 listed in Table 1. Observations with low and constant mutual information are plotted with black dashed lines, observations with oscillating mutual information are plotted with grey solid lines and those with intermediate behaviour are plotted by yellow dotted lines. The vertical line indicates the chosen time delay for our analysis. 

Open with DEXTER  
In the text 
Fig. A.2 The dependence of the ratio of false nearest neighbours on the embedding dimension m for the set of observations of IGR J170913624 for the choice R_{tol} = 10. 

Open with DEXTER  
In the text 
Fig. A.3 RP for several observations of IGR J170913624 from Table 1. 1. Row: left − SIMS/HIMS state 0200 (left upper corner in blue) and ρ state 0602 (right lower corner in red). Long diagonal lines appear in red RP. Right − the zoom of preceding plot. The diagonal lines are affected by the noise in such a way that the apparently long lines are composed of shorter squareshaped pieces. Next rows: in the red right lower halves of each plot RP are plotted for observations 0401 (2. row left, IVS state), 0402 (2. row right, ρstate), 0502 (3. row left, SIMS state) and 0701 (3. row right, ρstate). In the left upper half of each plot one of the corresponding surrogates is shown in green. 

Open with DEXTER  
In the text 
Fig. A.4 The cumulative histogram p_{c}(ϵ,l) for the observation 0602 computed for ϵ = 2.8. The length of the lines l is given in number of points of the lines (thus the total time of the line is given as t_{l} = lΔt). The estimate of the Rényi’s entropy K_{2} = 0.0699 ± 0.0003 is obtained as the linear regression of the central part of the cumulative histogram (indicated by a solid red line, the dashed red lines indicate the used part of the histogram). 

Open with DEXTER  
In the text 
Fig. A.5 Rényi’s entropy K_{2} versus recurrence rate for the observation 0602 for three different embedding dimensions indicated in the figure. 

Open with DEXTER  
In the text 
Fig. A.6 Dependence of the longest diagonal line L_{max} given in points on the recurrence threshold ϵ for the observation 0502 (thick red line). The ensemble of 100 surrogates are shown in grey and the ensemble of shuffled surrogates are shown in cyan. 

Open with DEXTER  
In the text 
Fig. B.1 The same as in Fig. A.6 for the observation 0501. The condition l> Δt on lines’ length used for K_{2}estimation is indicated by a black horizontal line. 

Open with DEXTER  
In the text 
Fig. B.2 The same as Fig. A.6, but for the observation 0602. 

Open with DEXTER  
In the text 
Fig. B.3 Comparison of recurrence rate for observation 0602 (value is indicated as the horizontal thick red line) and its hundred surrogates (number of surrogate is on the xaxis) computed for the same recurrence parameters (ϵ = 2.8,m = 10,Δt = 14dt). 

Open with DEXTER  
In the text 
Fig. B.4 Comparison of estimates of K2 for observation 0602 (indicated by horizontal thick red line) and its hundred surrogates (number of surrogate is on the xaxis) computed for the same recurrence parameters (ϵ = 3.25,m = 10,Δt = 14dt). 

Open with DEXTER  
In the text 
Fig. B.5 Comparison of Q for observation 0602 (indicated by the horizontal thick red line) and its hundred surrogates (number of surrogate is on the xaxis) computed for the same recurrence parameters (ϵ = 3.25,m = 10,Δt = 14dt). The mean value and the standard deviation of the ensemble of surrogates is indicated by horizontal lines. 

Open with DEXTER  
In the text 