© ESO, 2016

1. Introduction

The high energy radiation emitted by black hole X-ray 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 quasi-periodic 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 non-resonant class of models, the Lense-Thirring 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 X-ray 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 quasi-periodic 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 flare-like events of the microquasars, such as those observed in IGR 17091-3624 or GRS 1915+105 in their “heartbeat” states, are very strong, almost coherent, low-frequency QPOs. In these sources, they are so orderly and quasi-coherent 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 quasi-periodicity and to determine the nature of the variability of the X-ray 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 X-ray 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 limit-cycle 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 limit-cycle oscillations around a fixed point because of the two main types of thermal-viscous 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 X-ray transients, and generally its characteristic timescales are on the order of months in the case of stellar-mass 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 non-zero 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⁴ 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 self-regulation of the disk structure.

However, the so called “slim-disk” 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, long-term disk evolution that are based on the vertically averaged structure, clearly show the oscillatory type of behaviour with characteristic limit-cycles (Janiuk et al. 2002). The numerical computations performed in the 3D shearing-box 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 J17091-3624, 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 quasi-periodic flare-like 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 X-ray 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 α$\sqrt{{\mathit{P}}_{\mathrm{gas}}{\mathit{P}}_{\mathrm{tot}}}$, 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 (Kunert-Bajraszewska 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 so-called “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 d-dimensional 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 d-dimensional 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 counter-check with the correlation dimension technique and the single-value 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 X-ray lightcurves. Our goal is to find hints of deterministic chaotic behaviour in the accreting black hole systems, as observed by the advanced X-ray 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 J17091-3624 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 X-ray binary sources and we briefly describe the RXTE data-extraction 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 J17091-3624 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 recurrence-analysis method to the rest of our sample of the black hole X-ray 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 X-ray 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

The studied sources are the microquasars IGR J17091-3624, GRS 1915+105, GRO 1655-40, XTE J1650-500, XTE J1550-564, and GX 339-4. 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 339-4 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 flare-like 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 J17091-3624, whose spectral states have been studied by Pahari et al. (2014). The light curve of IGR J17091-3624 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 non-heartbeat states via the X-ray spectroscopy, suppress the oscillations completely. Therefore, the idea of nonlinear instability producing the flares in IGR J17091-3624 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 in-depth survey, we found indications of nonlinearity in several cases, apart from the two aforementioned microquasars.

2.1.1. IGR J17091-3624

The microquasar IGR J17091-3624 is a moderately bright transient X-ray 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 quasi-periodic flare-like 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 J17091-3624 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 J17091-3624 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 X-ray binary, discovered in 1992 (Castro-Tirado et al. 1992). Its behaviour was interpreted using the time-dependent 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 so-called “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¹⁸−3 × 10¹⁹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 339-4

X-ray nova GX 339-4, discovered in 1972 (Lavrov et al. 1997), exhibits high level of X-ray 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)²(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.

2.1.4. GRO J1655-40

Low mass X-Ray binary GRO J1655-40 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¹⁷ 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 J1550-564

XTE J1550-564 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 J1650-500

XTE J1650-500 is the black hole X-Ray binary with black hole mass 4−7.3 M_⊙ (Orosz et al. 2004; Slaný & Stuchlík 2008). The object emits flares with non-thermal energy spectra that occur when the persistent luminosity is near 3 × 10³⁴ erg s^-1(d/4 kpc)² 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³⁷ 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 (standard-1 mode or generic single bit mode) we extract the data using saextrct command, with a proper bin size. The standard-1 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 time-scales.

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 J1655-40).

3. Analysis of the chaotic processes

Our aim is to reveal important information about the BHXBs on the basis of their X-ray light curves using the recurrence analysis, which is a well-established 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, X-ray 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 X-ray binaries.

Here, the underlying dynamical system being studied is the accreting gas governed by a global physical parameter, e.g. by the mass-accretion 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 mass-accretion 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 non-stochastic 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¹ (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₂ estimation and other post-processing 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². 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 ${\mathit{K}}_{\mathrm{2}}^{\mathrm{obs}}$ for the observation and $\mathrm{\{}{\mathit{K}}_{\mathrm{2}}^{\mathrm{surr}}{\mathrm{\}}}_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{100}}$ for the set of surrogates. The Rényi’s entropy of the order α is, in general, defined as follows: ${\mathit{K}}_{\mathit{\alpha }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{1}\mathrm{-}\mathit{\alpha }}\ln \sum _{\mathit{i}\mathrm{=}\mathrm{1}\mathit{,N}}{\mathit{p}}_{\mathit{i}}^{\mathit{\alpha }}\mathit{,}$(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₂, 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 $\mathrm{𝒮}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathrm{=}\frac{{\mathit{N}}_{\mathrm{sl}}}{{\mathit{N}}^{\mathrm{surr}}}{\mathrm{𝒮}}_{\mathrm{sl}}\mathrm{-}\mathrm{sign}\mathrm{\left(}{\mathit{Q}}^{\mathrm{obs}}\mathrm{\right(}\mathit{ϵ}\mathrm{)}\mathrm{-}\mathit{Q̅}\mathrm{surr}\mathrm{(}\mathit{ϵ}{\mathrm{\left)}}^{\mathrm{\right)}}\frac{{\mathit{N}}_{{\mathrm{𝒮}}_{\mathit{K}}}}{{\mathit{N}}^{\mathrm{surr}}}{\mathrm{𝒮}}_{{\mathit{K}}_{\mathrm{2}}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathit{,}$(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₂ 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 ${\mathrm{𝒮}}_{{\mathit{K}}_{\mathrm{2}}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathrm{=}\frac{\mathrm{\left|}{\mathit{Q}}^{\mathrm{obs}}\mathrm{\right(}\mathit{ϵ}\mathrm{)}\mathrm{-}\mathit{Q̅}\mathrm{surr}\mathrm{(}\mathit{ϵ}\mathrm{\left)}\mathrm{\right|}}{{\mathit{\sigma }}_{{\mathit{Q}}^{\mathrm{surr}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}}}\mathrm{·}$(3)The significance defined in this way expresses how much the value ${\mathit{K}}_{\mathrm{2}}^{\mathrm{obs}}$ differs from the mean value $\overline{{\mathit{K}}_{\mathrm{2}}}{}^{\mathrm{surr}}$ measured in the units of the standard deviation of the set $\mathrm{\{}{\mathit{K}}_{\mathrm{2}}^{\mathrm{surr}}{\mathrm{\}}}_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathrm{100}}$ in the logarithmic scale σ_Q^surr(ϵ). 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), non-stochastic, but linear behaviour (observations with a significant result with shuffled surrogates, but non-significant 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 $\overline{\mathrm{𝒮}\mathit{,}}$ and the number of thresholds for which the average is computed.

4. Case study of IGR J17091-3624

We used the sample of observation of IGR J17091-3624 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 03-01, classified as IVS, have low and almost constant values of mutual information (see Fig. A.1). Three observations (04-01, 04-03 (IVS) and 06-03 (ρ)) 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₂ (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.

The observations 06-03 and 05-01 do contain enough long lines, so that the significance can be computed. In both cases, non-significant 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 06-03 and 05-01, respectively. Thus, we classify these two observations as a non-stochastic 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 08-00 to 3.72 for ObsID 06-02. 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 J1655-40, XTE J1650-500, XTE J1550-564, and GX 339-4, 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 μ), non-stochastic 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₂ 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 J17091-3624, 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 J17091-3624 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.

For the black hole candidate, GX 339-4, 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 soft-intermediate 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 non-stochastic behaviour in all three cases. For the last observation of GX 339-4, a non-significant 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 J1550-564, as taken on September 08, 1998. The source has been previously studied by Sobczak et al. (2000), and the lightcurve of observation 30191-01-14-00 shown therein looks similar to the variable state of IGR J17091-3624, 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 J1550-564 show non-significant 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 non-stochastic behaviour of the source during these states.

Figure 5 shows an exemplary lightcurve of the source GRO J1655-40, as taken by RXTE on August 01, 1996 (10255-01-04-00B). 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 non-stationary 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 non-significant results. For the latter group, we computed also the shuffled significance. The obsID 90058-16-07-00A again shows a non-significant result, hence it is of stochastic origin. Three other observations (90058-16-03-00A, 90058-16-03-00B, and 90058-16-02-00) yield moderate results , indicating a weakly non-stochastic origin. The shuffled surrogates of the last three observations (10255-01-04-00C, 20402-02-14-00A, and 20402-02-14-00B) have no long lines, leading to , which is strong evidence of non-stochastic behaviour of the source.

Finally, within the studied sample of observations of XTE J1650-500, no significant traces of nonlinear behaviour were found. Figure 6 shows one of its exemplary lightcurves, from the observation ID 60113-01-08-02, 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 X-ray sources. They are the black hole candidates, in which their X-ray luminosity originates in an accretion disk. The latter may be responsible for a short-timescale, quasi-regular 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 thermal-viscous instabilities occur. These instabilities lead to global limit-cycle 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.

Two microquasars, GRS 1915+105, and the recently discovered IGR J17091-3624, are already firm candidates for the limit-cycle 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, disk-dominated, 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 J17091-3624, 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 X-ray binaries other than these two microquasars, the situation is not so obvious. However, Janiuk & Czerny (2011) discussed the theoretical aspects of the thermal-viscous 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 X-ray transient source XTE J1650-500 displayed its only outburst in the year 2001, and there were not many observations available for that source. Homan et al. (2003) detect the high-frequency variability in XTE J1650-500 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 disk-jet-corona 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 X-ray 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 X-ray binary as a prospective candidate for disk instability-induced oscillations, because of these timescales. However, these flares have non-thermal 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^-5L_Edd. Therefore they are not good candidates for radiation pressure instability. Tomsick et al. (2003) also report on aperiodic oscillation with a 14-day time scale. They relate that these variations in X-ray 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 limit-cycle 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 X-ray satellites.

In XTE J1550-564, 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 J1550-564 is classified as a microquasar on the basis of its large-scaled moving jets, detected at X-ray and radio (Corbel et al. 2002), similar to GX 339-4 (Corbel & Fender 2002). We expect these two objects to have similar disk variability characteristics and, along with the two well-studied microquasars, nonlinear dynamical processes should also occur in these sources. Our current analysis confirms these expectations.

Initially, the spectral state of GRO J1655-40 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 X-ray outburst from 1996−1997, estimated the unabsorbed X-ray 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 J1655-40 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 J1655-40 and showed that, during the high/soft state, its spectrum is dominated by the soft thermal emission from the accretion disk. Comparing the two above-mentioned 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 & Klein-Wolt (2015) report on an unusual soft state of GRO J1655-40, observed by the RXTE during its 2005 outburst. Chandra X-ray grating observations have revealed a high mass-outflow 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 J1655-40 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 X-ray binaries detected by RXTE satellite.

• We developed a method to distinguish between stochastic, non-stochastic 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 J17091-3624, 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 instability-induced 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 well-studied binary GRS 1915+105, we also found significant traces of nonlinear dynamics in three sources (GX 339-4, XTE J1550-564, and GRO J1655-40).

• 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 low-dimensional 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 thermal-viscous instability and is undergoing induced limit-cycle 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 X-ray 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 X-ray 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 DEC-2012/05/E/ST9/03914 from the Polish National Science Center.

References

Appendix A: Recurrence analysis of the observed time series of IGR J17091-3624

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 ${y}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathrm{\left\{}\mathit{x}\mathrm{\right(}\mathit{t}\mathrm{\left)}\mathit{,x}\mathrm{\right(}\mathit{t}\mathrm{+}\mathrm{\Delta }\mathit{t}\mathrm{\left)}\mathit{,x}\mathrm{\right(}\mathit{t}\mathrm{+}\mathrm{2}\mathrm{\Delta }\mathit{t}\mathrm{\left)}\mathit{,}\mathit{...}\mathit{,x}\mathrm{\right(}\mathit{t}\mathrm{+}\mathrm{(}\mathit{m}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{\Delta }\mathit{t}\mathrm{\left)}\mathrm{\right\}}\mathit{,}$(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 J17091-3624. 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 well-known 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.

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₂, 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: ${{R}}_{\mathit{i,j}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathrm{=}\mathrm{\Theta }\mathrm{(}\mathit{ϵ}\mathrm{-}\mathrm{\parallel }{{x}}_{\mathit{i}}\mathrm{-}{{x}}_{\mathit{j}}\mathrm{\parallel }\mathrm{)}\mathit{,} \mathit{i,j}\mathrm{=}\mathrm{1}\mathit{,}\mathit{...}\mathit{,N,}$(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, 02-00 (SIMS/HIMS – upper left corner on the left panel in first row) and 06-02 (ρ-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.

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, $\mathit{P}\mathrm{\left(}\mathit{ϵ,l}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i,j}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}^{\mathit{N}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{R}}_{\mathit{i}\mathrm{-}\mathrm{1}\mathit{,j}\mathrm{-}\mathrm{1}}\mathrm{(}\mathit{ϵ}\mathrm{\left)}\mathrm{\right)}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{R}}_{\mathit{i}\mathrm{+}\mathit{l,j}\mathrm{+}\mathit{l}}\mathrm{(}\mathit{ϵ}\mathrm{\left)}\mathrm{\right)}\underset{\mathit{k}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{0}}{\overset{\mathit{l}\mathrm{-}\mathrm{1}}{\mathrm{􏽙}}}{\mathit{R}}_{\mathit{i}\mathrm{+}\mathit{k,j}\mathrm{+}\mathit{k}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathit{.}$(A.3)Based on the histogram, more quantifiers can be defined (Marwan et al. 2007). Here we focus mainly on ${\mathit{L}}_{\max }\mathrm{=}\max \mathrm{\left(}\mathrm{\right\{}{\mathit{l}}_{\mathit{i}}{\mathrm{\}}}_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}\mathrm{)}$, which is the length of the longest diagonal line (except for the main diagonal), and on the estimate of the second-order Rényi’s entropy K₂. 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 ${\mathit{K}}_{\mathrm{2}}\mathrm{=}\mathrm{-}\underset{\mathrm{\Delta }\mathit{t}\hspace{0.17em}\mathrm{\to }\hspace{0.17em}\mathrm{0}}{\lim }\underset{\mathit{ϵ}\hspace{0.17em}\mathrm{\to }\hspace{0.17em}\mathrm{0}}{\lim }\underset{\mathrm{\Delta }\mathit{l}\hspace{0.17em}\mathrm{\to }\hspace{0.17em}\mathrm{\infty }}{\lim }\frac{\mathrm{1}}{\mathit{l}\mathrm{\Delta }\mathit{t}}\ln \sum _{{\mathit{i}}_{\mathrm{1}}\mathit{,}\mathit{...}\mathit{,}{\mathit{i}}_{\mathit{l}}}{\mathit{p}}_{{\mathit{i}}_{\mathrm{1}}\mathit{,}\mathit{...}\mathit{,}{\mathit{i}}_{\mathit{l}}}^{\mathrm{2}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathit{,}$(A.4)where the phase space is coarse-grained into m-dimensional hyper-cubes of size ϵ. Here, p_i1,...,il describes the probability that the first point of the trajectory x(t = Δt) is inside the box i₁, the second point x(t = 2Δt) is in the box i₂ 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 phase-space boxes).

More precisely, K₂ 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 ${\mathit{p}}_{\mathrm{c}}\mathrm{\left(}\mathit{ϵ,l}\mathrm{\right)}\mathrm{~}{\mathit{ϵ}}^{{\mathit{D}}_{\mathrm{2}}}{\mathit{e}}^{\mathrm{-}\mathit{l}\mathrm{\Delta }\mathit{t}{\mathit{K}}_{\mathrm{2}}}\mathit{,}$(A.5)hence we can estimate the value of K₂ 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₂ 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 ϵ₂ = ϵ₁ + Δϵ, we can also estimate the value of the correlation dimension D₂ by ${\mathit{D}}_{\mathrm{2}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathrm{=}\ln \left(\frac{{\mathit{p}}_{\mathit{c}}\mathrm{\left(}{\mathit{ϵ}}_{\mathrm{1}}\mathit{,l}\mathrm{\right)}}{{\mathit{p}}_{\mathit{c}}\mathrm{\left(}{\mathit{ϵ}}_{\mathrm{2}}\mathit{,l}\mathrm{\right)}}\right)/\left(\frac{{\mathit{ϵ}}_{\mathrm{1}}}{{\mathit{ϵ}}_{\mathrm{2}}}\right)\mathrm{·}$(A.6)However, to obtain reasonable results for D₂, much more data points are needed than for estimating K₂ (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₂ by means of the recurrence analysis.

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 “should-be the recurrence point” out from the ϵ-neighbourhood, and on the other hand, the edges on the lines appear when noise brings the “should-not-be 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₂ is overestimated. An example of the decreasing dependence of the K₂-estimate on the recurrence rate with Δt = 7 s, and three different embedding dimensions m₁ = 10, m₂ = 20, and m₃ = 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₂. 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 J17091-3624, only the observations 05-01 and 06-03 can be assigned to this subgroup, the example for 05-01 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, ${\mathit{L}}_{\max }^{\mathrm{obs}}$ is higher than ${\mathit{L}}_{\max }^{\mathrm{surr}}$ for all surrogates for a wide range of ϵ. This behaviour indicates nonlinear features of the background dynamics.

To quantify these results, we compute the estimate of K₂ for observations with enough long lines and we use this quantity as our discriminating quantity. The computation goes as follows.

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₂ 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 06-02 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 ${\mathit{K}}_{\mathrm{2}}^{\mathrm{surr}}$ 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 ${\mathit{K}}_{\mathrm{2}}^{\mathrm{obs}\mathit{,}}$ and the mean of estimates of Rényi’s entropy for ensembles of the surrogates $\mathit{K̅}\begin{array}{c}\mathrm{surr}\\ \mathrm{2}\end{array}$. However, as can be seen from relation A.5, the quantity K₂ is non-negative and behaves in an exponential manner (for 10-times longer diagonal lines, K₂ drops down by one order). Therefore we first introduce the quantity Q = ln(K₂), which scales linearly, and then define the significance stemming from the estimate of K₂ as ${\mathrm{𝒮}}_{{\mathit{K}}_{\mathrm{2}}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathrm{=}\frac{\mathrm{\left|}{\mathit{Q}}^{\mathrm{obs}}\mathrm{\right(}\mathit{ϵ}\mathrm{)}\mathrm{-}\mathit{Q̅}\mathrm{surr}\mathrm{(}\mathit{ϵ}\mathrm{\left)}\mathrm{\right|}}{{\mathit{\sigma }}_{{\mathit{Q}}^{\mathrm{surr}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}}}\mathit{,}$(B.1)where σ_Q^surr is the standard deviation of the set ${\mathit{Q}}_{\mathit{i}}^{\mathrm{surr}}$ (see Fig. B.5).

When the observed data yields ${\mathit{K}}_{\mathrm{2}}^{\mathrm{obs}}$, 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 ${\mathit{K}}_{\mathrm{2}}^{\mathrm{obs}}$ is lower than $\mathit{K̅}\begin{array}{c}\mathrm{surr}\\ \mathrm{2}\end{array}$. 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 .

On the other hand, if ${\mathit{K}}_{\mathrm{2}}^{\mathrm{obs}}$ is higher than $\mathit{K̅}\begin{array}{c}\mathrm{surr}\\ \mathrm{2}\end{array}$, 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 $\mathrm{𝒮}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathrm{=}\frac{{\mathit{N}}_{\mathrm{sl}}}{{\mathit{N}}^{\mathrm{surr}}}{\mathrm{𝒮}}_{\mathrm{sl}}\mathrm{-}\mathrm{sign}\mathrm{\left(}{\mathit{Q}}^{\mathrm{obs}}\mathrm{\right(}\mathit{ϵ}\mathrm{)}\mathrm{-}\mathit{Q̅}\mathrm{surr}\mathrm{(}\mathit{ϵ}\mathrm{\left)}\mathrm{\right)}\frac{{\mathit{N}}_{{\mathrm{𝒮}}_{\mathit{K}}}}{{\mathit{N}}^{\mathrm{surr}}}{\mathrm{𝒮}}_{{\mathit{K}}_{\mathrm{2}}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathit{,}$(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 ${\mathit{K}}_{\mathrm{2}}^{\mathrm{obs}}\mathit{>}\mathit{K̅}\begin{array}{c}\mathrm{surr}\\ \mathrm{2}\end{array}$ 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₂ 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 05-01 and 06-03 (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 non-significant 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 non-stochastic 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 X-ray 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₂ 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 non-significant 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