Stochastic nonlinear self-oscillatory model of an accretion disk for the X-ray bursting of the microquasar GRS 1915+105

. The observed light curves of the microquasar GRS 1915+105 show three main characteristics, namely chaotic, stochastic and semi-stochastic (or nonstochastic). In order to probe the origin causing these diﬀerent features, we propose a stochastic nonlinear self-oscillatory model of an accretion disk, and use its oscillating luminosities with diﬀerent noise strengths and nonlinear factors to simulate the observational variability, and compare their correlation dimensions calculated by a nonlinear time series analysis. It is shown that when the stochastic noise strength is small and the nonlinear factor is large, the luminosity curves show chaotic behavior. On the contrary, it shows stochastic or semi-stochastic behavior when a stochastic factor dominates. The simulated data exhibits the same nonlinear and random characteristics as the observational sequences, and their corresponding correlation dimensions are also consistent. This indicates that the X-ray bursting of the microquasar should be related to a nonlinear self-oscillation of the accretion disk.


Introduction
It is well known that the microquasar GRS 1915+105 is a prominent black hole system exhibiting a large variable X-ray emission on different time scales and patterns ( Belloni et al. 1997;Chen et al.1997;Paul et al. 1997).Based on a large collection of multi-epoch RXTE observations, Belloni et al. (2000) classified its observed light curves into 12 temporal states, of which four classes show chaos, three are random, while five others show a deviation from being random (Mukhopadhyay 2004;Harikrishnan 2010).That is, the 12 temporal classes could be grouped into three main types, namely chaotic behavior, stochastic behavior and semi-stochastic (or nonstochastic) behavior (Misra et al. 2004(Misra et al. , 2006)), respectively.Many studies suggest that the complexities of light curves may be caused by the noise fluctuations in a chaotic system involving nonlinear processes, whose evolution can be described by several ordinary differential equations.For example, Massaro et al. proposed two nonlinear mathematical models, Fitzhugh-Nagumo (FhN) and Modified Hindmarsh-Rose (MHR) models, and used their solutions to reproduce several different classes of light curves of GRS 1915+105 (Massaro et al. 2014(Massaro et al. , 2020(Massaro et al. a, 2020(Massaro et al. b, 2020 c c).The MHR model described well the transition from stable to bursting states and the occurrence of low-frequency quasi periodic oscillations as a consequence of a transition from an unstable to a stable equilibrium.However, the physical explanation of these nonlinear differential equations is not completely understood.
As an important theoretical model, the oscillations of an accretion disk can interpret a large number of the observational features of black hole binary systems (Kato 2016).In the last decade, the stochastic oscillation of an accretion disk has attracted great attention in astrophysical studies, and has succeeded in explaining some observations.For example, a generalized Langevin equation with white or color noise has been used to describe the vertical oscillations of accretion disks, and the luminosity of stochastic oscillations could reproduce the observations of optical intra-day variability of BLLac objects ((Harkoet al. 2012((Harkoet al. , 2014;;Leung et al. 2011Leung et al. , 2014;;Long et al. 2016).Based on their model, we have discussed the stochastic resonance phenomenon in power spectral density curves, and considered that the resonance peak is an alternative explanation of the observed low-frequency quasi periodic oscillations (Wang et al. 2013(Wang et al. a, 2013(Wang et al. b, 2015(Wang et al. , 2016)).However, the light curves of this stochastic oscillating model only show stochastic behavior, not present periodicity and chaotic behavior.So Ou et al. (2014) developed a chaotic+stochastic oscillating accretion disk model, in which the chaotic factor was introduced into the stochastic oscillating disk.The stable response of this model is completely controlled by the external random and chaotic driving of the system, so it could simulate the light curves showing chaotic behavior of variability, and their correlated dimensions are consistent with that of the observational data.
The oscillations of an accretion disk surrounding a compact object are generally driven by magneto-hydrodynamic turbulence which is an intrinsically nonlinear process.The model for this process should be nonlinear and is expected to show qualitative changes in its behavior as a control parameter is varied.So in this paper, we improve the stochastic oscillation model mentioned above, and suggest a stochastic nonlinear self-oscillatory model of an accretion disk to simulate the light curves of GRS 1915+105.In contrast with linear stochastic oscillation model, self-oscillations of a constrained cycle can generate and maintain a regular mechanical periodicity without requiring a similar external periodicity to drive it (Jenkins 2013).Such constrained cycle oscillations usually occur in non-conservative dissipative systems and have been considered in the dynamics of Cepheids and RR Lyrae stars (Rudd, Rosenberg 1970;Buchler, Kovacs 1987;Das et al. 1996).Nakariakov et al. (2016) use a limit cycle solution corresponding to the self-oscillations to interpret the undamped quasi-harmonic kink oscillations of coronal loops, and successfully reproduce their observed properties.On some occasions, the complex variability of GRS 1915+105 is characterized by a long series of bursts expected by limit cycles (Castro-Tirado et al. 1992).In addition, nonlinear effects in the model may lead to overstable pseudo-periodic, or even chaotic motions.So we will calculate the dissipated energy of this selfoscillating system to simulate the X-ray variability classes of GRS 1915+105, and probe the intrinsic relation between the simulated data with observational variability by comparing their correlated dimensions.

Stochastic nonlinear self-oscillatory model of accretion disk
For a thin accretion disk around a black hole in contact with an isotropic and homogeneous external heat bath, the vertical oscillation of unit mass can be described by a generalized Langevin equation, given by Leung et al.( 2011) where z and M D are the displacement of oscillation and the mass of the accretion disk, respectively.µ expresses a viscosity coefficient and dF (t) is a stochastic excitation.The intrinsic angular frequency ω 2 0 can be written as (Titarchuk, Osherovich 2000) where, m = M 0 /M , M 0 and M are the masses of the central object and the sun.r in = R in /3R S is the innermost radius of the disk and R S = 2GM 0 /c 2 is the Schwarzschild radius.r adj = R adj /R in and r out = R out /R in are the adjustment radius and the outer radius of the disk, respectively.The index γ is either 3/5 or 3/4 (Shakura, Sunyaev 1973).
The accretion disk as a whole is subjected to viscous forces mainly from friction between the disk and the external heat bath.If the interaction of the disk with the surrounding moving medium is not stiff but slippery, this friction is a possible way to compensate for the dissipation losses (Nakariakov et al. 2016).Further, if we incorporate a cubic nonlinearity into the dissipation, then the damping term of the system may be expressed as −µ 0 dz dt + µ 0 az 2 dz dt (µ 0 and a are positive constants).The first term is the linear friction force, by which energy is transferred from the external heat bath to the disk, and the second one is the nonlinear dissipation force resulting in the energy loss of the system.
In addition, if we consider that the self-gravitational potential V z of the disk is 1 2 ω 2 0 z 2 + 1 4 ω 2 0 bz 4 , then the nonlinear gravitational term of a disk oscillation may be represented as ω 2 0 (z + bz 3 ).We also assume random excitation F (t) = DdW (t), where D refers to the noise strength and dW (t), defining the normalized δ-correlated Gaussian noise with the zero mean and unit variance, may be represented by the formal derivative of a Wiener process, W (t). Then we replace the equation (1) of the disk motion with a stochastic nonlinear differential equation, which can be written as (3) Self-oscillatory motions described by equation ( 3) are essentially dissipative, with the dissipative losses compensated by the continuous energy supply of the external heat bath and the radiation luminosity of the disk, representing the energy lost by the disk due to viscous dissipation and the presence of the random force, is given by: Obviously, because equation (3) contains nonlinear terms and a random term, the luminosity variables can show different behavior, which should be similar to different temporal states of GRS 1915+105.

Light curves of observation sequences and simulated data
We show in this section that the numerical solutions of the nonlinear oscillation system of equations ( 3) and ( 4) reproduce the main features of three main variability types of GRS 1915+105.Based on the PUBLIC in the RXTE TOO archive (at heasarc.gsfc.nasa.gov),here we choose three representative observation sequences whose observation IDs are 20402-01-02-02, 10408-01-40-00 and 10408-01-08-00, respectively.Their continuous light curves in the energy range 2-30keV binned at 1s intervals are plotted in Figures 1(a),(c) and(e), which show the stochastic behavior, semi-stochastic behavior and chaotic behavior, respectively.
In order to simulate the light curves of three observational sequences, we perform numerical computations of equations ( 3) and ( 4) by means of a Runge-Kutta fourth order integration routine, and obtain luminosities with different parameters.For GRS1915+105, we set the same disk parameters as Titarchuk and Osherovich (2000), r in = 1, r out = 10 4 and γ = 3/5, r in = 3, and adopted m = 14 (the mass range of this source is from 10M to 18M as estimated by Remillard and McClintock (2006)).Then the intrinsic oscillating frequency of the disk can be calculated, ω 2 0 ≈ 1.5.In the following we will calculate the self-oscillating luminosities of the accretion disk by adjusting the nonlinear term parameter b and the random strength D, with the other three kept fixed, µ 0 = 0.5, a = 1.0, ω 2 0 = 1.5.Figures 1(b), (d term makes the self-oscillating luminosity curves of the accretion disk appear chaotic, random and semi-random.The nonlinear factor makes the system practically appear as chaotic, but noise always suppresses it.Therefore, when the noise fluctuations are much smaller than the variability, the nonlinear factor is dominant, and the chaotic nature can be detectable.Similarly, the random and semi-random appearances are due to the dominance of the noise in the system.

Correlation dimension D 2 of light curves
In order to quantitatively compare the similarity between the simulated luminosities and the observed data, we use a nonlinear time series analysis to compute their correlation dimension D 2 , which is often used as a discriminating parameter for a hypothesis testing to detect nontrivial structures in the time series.

Correlation dimension D 2
The conventional method for the calculation of D 2 is the delay embedding method known as the Grassberger-Proccacia (GP) algorithm (Grassberger & Procaccia 1983;Harikrishnan et al.. 2006).In this technique, an embedding space of dimension M with delay vectors is constructed by splitting a discretely sampled scalar time series s(t i ) where τ is the delay time computed by means of the mutual information method (Kraskov et al. 2004), assuming that the vectors are not correlated.The correlation function for the embedded time series is the number of points that are within a distance R from the center, averaged over all the centers, and may be written as where N ν = N − (M − 1)τ is the total number of reconstructed vectors, N is the size of the data set and H is the Heaviside step function.The correlation dimension D 2 (M ) is defined as In order to test the computational schemes presented here, we first apply the D 2 analysis on the standard chaotic time series and the pure noise from the standard Lorenz attractor and stochastic white noise.Their correlation dimension curves calculated by this method are shown in Figure 2. We can see, as expected, that the D 2 (M ) values for the Lorenz system saturate at M = 2 to D sat 2 ≈ 2, which is close to the known value of 2.04.For the pure uncorrelated stochastic white noise, D 2 (M ) varies approximately linearly with M (D 2 ≈ M ) and no such saturation exists.This is consistent with the results of Misra et al. (2006), and confirms that the random data and the low-dimensional chaotic system can clearly be distinguished in this scheme.this ideal random curve, and there is the same saturation correlation dimension of D 2 = 3.693 for M = 12.
For observational sequences 20402-01-02-02 and 10408-01-40-00, their curves show a different deviation from the ideal random noise (see Figure 3(b)), and D 2 varies approximately linearly with M and no such saturation exists.It can be seen that for the observational sequence 20402-01-02-02, its D 2 (M ) curve is close to that of the ideal random noise, and is random dominated.While the other one deviates greatly from it and shows semi-stochastic behavior.As expected, the correlated dimensions of these observational variabilities are consistent with the results of the simulated luminosities for different model parameters D and b.

Summary and discussion
Many studies suggest that the variability of the micro quasar GRS 1915+105 may be chaotic in nature, this chaotic signature being suppressed in some cases due to the noise dominance when it appears as random or semi-random, but actually this may not be its fundamental characteristic.Therefore, in the present paper we have considered a stochastic nonlinear self-oscillatory equation to describe the vertical oscillations of an accretion disk around a black hole, and obtained the oscillating luminosities showing chaotic, random or semi-random behavior for different noise strengths and nonlinear parameters.Their properties are similar to those of the observed variability of GRS1915+105.
Our results show that the correlation dimensions between observational sequences and simulated data are identical approximately, which implies that the variability of the microquasar may be controlled by the nonlinear self-excited oscillation system; at the same time, it is also affected by the noise with varying intensity, resulting in a variety of variability types.So, as a possible astro-physical application, this model may turn out to be very useful to understand complicated temporal behavior of X-ray bursts.

Figure 2 .Figure 3 .
Figure 2. D2(M ) curves for the random noise (crosses) and for a Lorenz system (triangles), the straight line represents a D2 = M curve.