Quantifying the memory and dynamical stability of magnetar bursts

Magnetars are neutron stars with extremely strong magnetic fields exceeding $ \sim 10^{13} $ G [1−3]. They have long rotational periods that typically last for several seconds and gradually spin down due to electromagnetic radiation. Magnetars are observationally recognized as soft gamma repeaters (SGRs), which persistently emit hard X-rays and soft gamma-rays [4−6]. The magnetar bursts usually have duration of $ \sim 0.1-1 \; \rm{s} $ and peak luminosity in the range of $ \sim 10^{39}-10^{41} \; \rm{erg s}^{-1} $. Although it is widely accepted that the bursts are powered by the strong magnetic fields, the triggering mechanism remains unclear. Many theoretical models have been proposed to explain the triggering mechanism of SGRs, such as the crustquakes of neutron star [7] and the magnetic reconnection [8].

There are already several studies on the statistical properties of SGRs, particularly focusing on the distributions of energy and waiting time [9−15]. The cumulative energy distribution of 111 bursts from SGR 1806−20 was found to be well described by a power-law function with an index $ \gamma = 1.66 $, which is very close to the index $ \gamma \approx 1.6 $ of the earthquake Gutenberg–Richter power law [9]. Chang et al. [12] investigated the cumulative distributions of SGR J1550−5418 and found that the distributions of fluence, peak flux and duration are well fitted by a bent power law, while the distribution of waiting time follows a simple power law. In addition to the power-law distribution of magnetars bursts, the fluctuations of bursts show scale-invariant properties [12, 14−16]. The probability density functions of fluctuations are well described by the Tsallis q-Gaussian function for fluence, peak flux, and duration of 384 bursts from SGR J1550−5418, and the q values are consistent across different scale intervals, which indicates that the Tsallis q-Gaussian distribution of the fluctuations is scale-invariant [12].

The statistical similarity between SGRs and fast radio bursts (FRBs) has also been discussed [14, 15], which is observationally motivated by the association between the Galactic FRB 200428 [17, 18] and the hard X-ray burst from SGR 1935+2154 [19−23]. The repeating FRBs show properties similar to those of SGRs, including a power-law distribution for both energy and waiting time [24−29], as well as a scale-invariant Tsallis q-Gaussian distribution of their fluctuations [14, 27, 28, 30, 31]. The power-law distribution of energy and the scale-invariant fluctuations are predicted by the self-organized criticality (SOC) systems [32, 33], which provides a potential explanation for the burst properties observed in both SGRs and repeating FRBs.

The statistical properties of the time series of repeating FRBs and SGRs have also been investigated recently [30, 34, 35]. Particularly, Zhang et al. [34] employed the Pincus index and Lyapunov exponent to quantify the randomness and chaos of repeating FRB 20121102A and FRB 20190520B, and compared them with other natural phenomena such as pulsars, earthquakes, solar flares and Brownian motion. The repeating FRBs were found to exhibit high randomness and low chaos, mimicking the behavior of Brownian motion in the randomness-chaos phase space. The same dynamical stability analysis was also applied to magnetar bursts from SGR J1550-5418 and SGR J1935+2154. It was found that these magnetar bursts have a distinct separation from FRBs in time domain, but has no significant difference compared to FRBs in energy domain [35]. Another interesting phenomenon identified through the analysis of time series is the long-term memory, which has been observed in repeating FRBs [30, 36, 37], but has not yet been studied in magnetar bursts.

In this paper, we investigate the statistical properties of magnetar burst using four active SGRs, which consists of approximately 2000 bursts in total. We particularly focus on the long-term memory and dynamical stability of SGRs. The structure of this paper is arranged as follows: The dataset used in this paper is presented in Section 2. The long-term memory is studied in Section 3. The randomness and chaos are analyzed in Section 4. Finally we give the discussion and conclusions in Section 5.

In the analysis, we use the burst data from four SGR sources: SGR 1806−20, SGR 1900+14, SGR J1935+2154, and SGR J1550−5418. The first sample consists of 924 bursts from SGR 1806−20 detected by the Rossi X-ray Timing Explorer (RXTE) between 1996 and 2011. The second sample consists of 432 bursts from SGR 1900+14 detected by RXTE between 1998 and 2006. Both samples are available at the online database ¹ constructed by the high-energy astrophysics group at Sabancı University. The burst time and total count of each bursts can be found in the database. The third sample is 217 bursts from SGR J1935+2154 observed by the NICER telescope during the 1120s burst storm period on 2020 April 28 [38]. The database provides the burst start time, duration and flux of the bursts. We can also calculate the fluence by the product of duration and the time-averaged flux. The fourth sample is 384 bursts from SGR J1550-5418 observed in the three active episodes in 2008-2009 by the Gamma-ray Burst Monitor onboard the Fermi Gamma-ray Space Telescope (Fermi/GBM) [39]. The $ T_{\rm 90} $ start time and fluence are given in the database. All the four samples are shortly summarized in Table 1. The photon counts or fluence as a function of arrival time are depicted in Figure 1.

Figure 1.  (color online) Photon counts or fluence vs. the arrival time for all the four SGR samples. The arrival time of the first burst is set to zero.

In this paper, we investigate the long-term memory and dynamical stability of magnetar bursts through analyzing the time sequence of waiting time and energy. The waiting time is defined as the time difference between two successive bursts, $ \Delta T = T_{i+1} - T_i $. We use the burst time of SGR 1806−20 and SGR 1900+14, the burst start time of SGR J1935+2154, and the $ T_{\rm 90} $ start time of SGR J1550−5418 as the arrival time of the i-th burst $ T_i $. Note that each sample consists of different observing sessions. The duration of an observing session is expected to be no more than a half of the orbital periods of satellites. The orbital periods of satellites around the earth are 93, 93 and 95 minutes for RXTE, NICER and Fermi satellites, respectively. Therefore, we use a uniform criterion to discard the waiting time larger than 1 hour to avoid the long observing gaps. The burst energy is not a directly observable quantity. For a specific SGR source, the burst energy is proportional to the fluence or the photon counts. Rescaling the burst energy by a universal constant does not affect the results we discussed bellow. Hence we directly use the fluence or the photon counts as the representation of burst energy. Specifically, we use the time sequence of total count in SGR 1806−20 and SGR 1900+14 samples, and the time sequence of fluence (flux multiplied by duration) in SGR J1935+2154 and SGR J1550−5418 samples as the energy sequence in the following analysis.

The long-term memory observed in the time series data is a fascinating phenomenon which has already been identified in several natural events, such as earthquakes [40], solar flares [41] and repeating FRBs [30, 36, 37]. In this section, we present the measurement of long-term memory in magnetar bursts for the first time. We employ the rescaled range analysis (R/S analysis) to calculate the Hurst exponent H [42, 43], which is a quantification of the long-term memory in time series data. If the rescaled range analysis yields $ H > 0.5 $, it indicates the presence of positive long-range correlations in the time series data. Conversely, if $ H < 0.5 $, negative long-range correlations are presented. If $ H = 0.5 $, it suggests an absence of long-range correlations within the data.

Here we give a concise introduction of the rescaled range analysis method [44−46]. A time series of length N is divided into l non-overlapping subseries, each of length n. For each subseries $ X_m $, where $ m=1,2,...,l $, the rescaled range analysis follows the steps bellow: (a) find the mean value $ E_m $ and the standard deviation $ S_m $ for the subseries $ X_{m} $; (b) normalize the data in subseries $ X_{m} $ by subtracting the mean value to obtain a mean-adjusted time series, $ Y_{i,m}=X_{i,m}-E_m $ for $ i=1,...,n $; (c) construct a cumulative deviation time series given by $ Z_{i,m}=\sum_{j=1}^i Y_{j,m} $ for $ i=1,...,n $; (d) determine the series range $ R_m = \max \{Z_{1,m},...,Z_{n,m}\} - \min \{Z_{1,m},...,Z_{n,m}\} $; and (e) rescale the range using the standard deviation, $ R_m/S_m $. Finally we calculate the mean value of the rescaled range for all subseries of length n,

By varying the length n of subseries, we construct a series of the rescaled range, which is asymptotically following the relation

where H is the Hurst exponent, which can be derived through linear regression method in the logarithmic scale,

In this paper, we employ the R/S analysis on the time series of waiting time and energy of magnetar bursts. The public package $\textsf{NOLDS}$ [47] is used in the analysis. Figure 2 illustrates the rescaled range series as a function of n for four SGR samples, with waiting time represented in green dots and energy in magenta squares. The points correspond to the rescaled range series derived from the sample data and the straight lines represent the results of the linear regression. As we expect, the rescaled range series is a simple power-law function of n. The data points fit well with a straight line in the log-log plot, where the slope of the line gives the Hurst exponent H. The best-fitting lines obtained using the usual least-$ \chi^2 $ method are shown in Figure 2, where the shaded regions represent the $ 1\sigma $ uncertainty. The Hurst exponents H for waiting time and energy of the four SGR samples are summarized in Table 2. We also calculate the intrinsic scatter $ \sigma_\mathrm{int} $, which is defined by the square root of the reduced chi-square. For the samples of SGR 1806−20, SGR 1900+14, SGR J1935+2154 and SGR J1550−5418, the Hurst exponents are $ H=0.73 \pm 0.02 $, $ 0.64 \pm 0.03 $, $ 0.67 \pm 0.03 $ and $ 0.73 \pm 0.03 $ for waiting time, and $ H=0.56 \pm 0.01 $, $ 0.58 \pm 0.02 $, $ 0.64 \pm 0.02 $ and $ 0.57 \pm 0.02 $ for energy, respectively. The Hurst exponents are larger than $ 0.5 $ at $ \gtrsim 3\sigma $ confidence level for all the samples. Therefore, we conclude that long-term memory exists in the time series of both waiting time and energy for the magnetar bursts.

Figure 2.  (color online) The rescaled range series of waiting time and energy for SGRs. The dashed lines and solid lines are the linear regression curves of waiting time and energy, respectively. The shaded regions represent the $ 1\sigma $ uncertainty of the linear regression.

	Waiting time	Energy
SGR 1806−20	0.73 $ \pm $ 0.02	0.56 $ \pm $ 0.01
SGR 1900+14	0.64 $ \pm $ 0.03	0.58 $ \pm $ 0.02
SGR J1935+2154	0.67 $ \pm $ 0.03	0.64 $ \pm $ 0.02
SGR J1550−5418	0.73 $ \pm $ 0.03	0.57 $ \pm $ 0.02

Table 2.  The Hurst exponent H of waiting time and energy for SGRs. The uncertainty is the intrinsic scatter of the linear regression.

In this section, we explore the dynamical stability of magnetar bursts in the context of randomness and chaos. The randomness and chaos of a time series are quantified by the Pincus index (PI) [48] and the largest Lyapunov exponent (LLE) [49], respectively. We compute the PI and LLE values for the waiting time and energy in four SGRs samples and illustrate the results in randomness-chaos phase space.

The PI value is a quantification of the randomness in a dynamical system, and its calculation is based on the concept of approximate entropy. The approximate entropy of a time series $ {u(i)} $ with length N is defined by [48, 50]

Here $ x_{m}(i)=[u(i),..., u(i + m-1)] $ and $ x_{m}(j)=[u(j),..., u(j+m-1)] $ are the subseries of $ {u(i)} $, and the embedding dimension m is the length of the subseries. $ {\rm dist}[x,y] $ is the Chebyshev distance between x and y. $ \theta(x) $ is the step function, i.e. $ \theta=1 $ for $ x\geq 0 $, and $ \theta=0 $ for $ x<0 $. We use the public code $\textsf{EntropyHub}$ [51] to compute ApEn in this paper. We follow the traditional convention to take the embedding dimension $ m=2 $ and vary the distance threshold r in the range of [0.01,0.09] multiplied by the standard deviation of the data series. The maximum approximate entropy (MAE) [50, 52] is defined as the maximum value of ApEn across varying the value of r,

Through bootstrap sampling we can compute PI for a series, which is defined as the ratio of MAE of the original series to that of the randomly shuffled series,

The PI quantify the randomness of a time series, with a value of zero indicating a completely ordered system and a value of unity representing complete randomness. In this paper, we perform random shuffling on the original series 1000 times, so obtain 1000 values of $ {\rm MAE_{shuffled}} $. The bootstrap sampling also gives the uncertainty of PI, which is propagated from the uncertainty of $ {\rm MAE_{shuffled}} $,

where $ \sigma_{{\rm MAE}_ {\rm{shuffled}}} $ is defined as the standard deviation of the $ {\rm MAE} $ values of the randomly shuffled series.

We calculate the MAE and PI values for the time series of waiting time and energy for the four SGR samples. In Figure 3, we show the MAE value for the original series in magenta solid line, and the distribution of the MAE values for the 1000 shuffled series in green histogram. The 16th, 50th and 84th percentiles of the distribution are plotted in blue, red and black dashed lines, respectively. In Table 3, we list the PI values and their uncertainties. We can quantify the randomness of the samples through the distance between $ {\rm MAE}_{ \rm{original}} $ and $ {\rm MAE}_{\rm{shuffled} }$, or directly through the values of PI. As can be seen, the waiting time series of SGR 1900+14 and SGR J1550−5418 significantly deviate from a random organization, with PI values of $ 0.84 \pm 0.02 $ and $ 0.76 \pm 0.03 $, respectively. For SGR 1806−20 and J1935+2154, the waiting time series are also slightly different from a totally random organization. Hence, we conclude that the waiting time series are not totally random for all the four SGRs samples. As for the energy series, the PI value of a totally random organization ($ {\rm PI}=1 $) is inside the $ 1 \sigma $ confidence interval of the PI values of SGR 1806−20 and SGR 1900+14. In contrast, the PI values of SGR J1935+2154 and SGR J1550−5418 deviates from complete randomness ($ {\rm PI}=1 $) at more than $ 2 \sigma $ confidence level. Therefore, the energy series is totally random for SGR 1806−20 and SGR 1900+14 samples, but it is less random for SGR J1935+2154 and SGR J1550−5418 samples.

Figure 3.  (color online) The MAE distribution for SGRs. The magenta solid line is MAE of the original series. The green histogram is the distribution of MAE of the 1000 shuffled series. The blue, red and black dashed lines are the 16th, 50th, and 84th percentile of the 1000 simulations, respectively.

	Waiting time	Energy
SGR 1806−20	0.92 $ \pm $ 0.02	0.98 $ \pm $ 0.02
SGR 1900+14	0.84 $ \pm $ 0.02	0.98 $ \pm $ 0.03
SGR J1935+2154	0.92 $ \pm $ 0.03	0.91 $ \pm $ 0.04
SGR J1550−5418	0.76 $ \pm $ 0.03	0.95 $ \pm $ 0.02

Table 3.  The PI values of waiting time and energy for SGRs. The uncertainty is the standard deviation of the randomly shuffled series.

The LLE is a quantification of chaos in a non-linear dynamical system, and it quantifies the local stability features of attractors and other invariant sets in phase space. In the $ m- $dimensional phase space, an initially infinitesimal m-sphere of radius $ r_0 $ will deform into a m-ellipsoid due to the locally deforming nature of the phase flow. The length of the i-th ellipsoidal principal axis evolves exponentially with time $ r_i=r_0\exp(\lambda_it) $. The Lyapunov exponents is then defined by [49]

LLE is the maximum value of Lyapunov exponents $ \lambda_i $, which quantifies the rate of separation of two adjacent trajectories in phase space. A positive LLE value means that two infinitely close trajectories in phase space will diverge exponentially as time elapses, thus implies the existence of chaos. On the contrary, a negative value of LLE suggests a stable system. We use the algorithm of Rosenstein et al. [53] implemented in the public package $\textsf{NOLDS}$ [47] to calculate LLE, where the embedding dimension is chosen to be the default value $ m=10 $. Since the LLE is the maximum value of the whole spectrum of Lyapunov exponents, it is hard to define the uncertainty. One possible way is to calculate the LLE values by varying the embedding dimension m, and use the standard deviation of LLE values to represent the uncertainty. However, we find that the LLE is almost independent of the embedding dimension (see bellow). So we ignore the uncertainty of LLE here. Table 4 summarizes the LLE values of waiting time and energy series in the four SGRs. Although all the LLE values are positive, they are very close to zero. Therefore, we conclude that there is no strong evidence for the existence of chaos for all the four SGR samples.

In Figure 4, we plot the LLE and PI values for each SGR in the randomness-chaos phase plane. As a comparison, we also illustrate the LLE and PI values of four samples from three extremely active repeating FRBs observed by the FAST telescope, i.e. FRB 20121102A, FRB 20201124A and FRB 20220912A. The bursts from FRB 20201124A are divided into two samples, which are observed in two separate periods with about three months observational gap. The PI values of the four FRB samples were calculated in our previous work [31], and the LLE values are calculated by reanalyzing the same data samples using the method as applied in this paper. In the time domain (the left panel of Figure 4), except for SGR J1550−5418, the rest three SGR samples and all the four FRB samples concentrate on a small region (with the average values $ {\rm PI}\sim 0.9 $ and $ {\rm LLE}\sim 0.05 $). SGR J1550−5418 seems to be an outlier in the randomness-chaos phase plane, it has a smaller PI ($ \sim 0.76 $) and a larger LLE ($ \sim 0.14 $) than the other SGRs and FRBs. Hence the waiting time of SGR J1550−5418 are less random but more chaotic than other samples. In the energy domain (the right panel of Figure 4), all the data samples seem to concentrate on a small region in the randomness-chaos phase plane (with the average values $ {\rm PI}\sim 0.9 $ and $ {\rm LLE}\sim 0.06 $). Especially, the two samples from FRB 20201124A do not significantly diverge in the randomness-chaos plane, implying that there is no strong temporal variability in the burst activity. This result is consistent with our previous conclusion [31]. In summary, we conclude that there is no significant difference between SGRs and FRBs in the randomness-chaos phase plane, for both waiting time and energy.

Figure 4.  (color online) The randomness-chaos plane of waiting time (left-panel) and energy (right panel) for SGRs and FRBs. The dots and squares represent for SGRs and FRBs, respectively.

Both the PI and LLE values depend on the embedding dimension m. We use different embedding dimension m in the calculations of PI and LLE to assess the robustness of our results. The PI and LLE values as a function of m are depicted in Figure 5 for waiting time and energy of the four SGR samples. As can be seen, both PI and LLE values do not significantly affected by the choice of embedding dimension, implying the robustness of our results.

Figure 5.  (color online) The PI and LLE values vs. the embedding dimension m.

In this paper, we investigated the statistical properties of waiting time and energy of magnetar bursts using four samples. Using the rescaled range analysis to calculate the Hurst exponent, we measure the long-term memory in magnetar bursts for the first time, and found that long-term memory exists in the time series of both waiting time and energy of magnetar bursts. We investigated the dynamical stability of magnetar bursts through the measurement of randomness and chaos, which are quantified by the Pincus index and the largest Lyapunov exponent, respectively. In the time domain, all the four SGR samples show evidence for the deviation from random organization. In the energy domain, however, the result is somewhat subtle: SGR 1806−20 and SGR 1900+14 is totally consistent with random organization, but SGR J1935+2154 and SGR J1550−5418 seem to be less random. We also found that both waiting time and energy exhibit weak chaos in all the four SGR samples. Finally, we compared SGRs with repeating FRBs and found that there is no significant difference between these two astronomical phenomena in the randomness-chaos phase space, for both waiting time and energy.

Previous discussions on the statistical similarity between SGRs and repeating FRBs focus on the power-law distribution of the observed quantities such as energy, waiting time and burst duration, or the scale-invariant Tsallis q-Gaussian distribution of their fluctuations. Here we extend the scope of the investigation to the long-term memory and randomness-chaos characteristics, and we also found that SGRs share similar properties with repeating FRBs. The Hurst exponent H for repeating FRBs are found to be around 0.6 [30], consistent with the H values for SGRs obtained in this paper. Therefore the long-range correlations in SGRs and FRBs are positive but not strong. The PI values for FRBs and SGRs are also very similar, $ {\rm PI}\sim 0.8-1.0 $. These statistical similarities, along with the observation of the FRB-associated magnetar SGR 1935+2154 [17, 18], imply that there may be common emission mechanism between SGRs and FRBs.

The repeating FRBs were also compared with other physical phenomena in the randomness-chaos phase space. For instance, Zhang et al. [34] found FRBs deviate significantly from pulsars, earthquakes, and solar flares. We found that, for both waiting time and energy, SGRs do not deviate significantly from FRBs within the randomness-chaos phase space. The randomness and chaos have also been investigated for two SGRs by Yamasaki et al. [35]. In their work, the authors compared SGRs and FRBs in the randomness-chaos plane, and found that SGRs exhibit significantly lower randomness and a slightly higher degree of chaos compared to FRBs in time domain, but they exhibit a broad consistency in the energy domain. Our results, however, show that SGRs and FRBs are consistent with each other in both time domain and energy domain, which conflict with the results of Yamasaki et al. [35]. One reason may be the different data samples we used. Another reason may be that Yamasaki et al. [35] used the energy fluctuation in the analysis, while we used the energy itself. This issue needs to be further investigated with a larger data sample in the future.