A new algorithm for detecting X-ray shots in Cyg X-1

The Galactic X-ray binary Cygnus X-1 (Cyg X-1) consists of a black hole with a mass of $21.2\pm 2.2\,M_{\sun}$ and a high-mass companion star. The binary system has an orbital period of $5.6~{}\mathrm{d}$ and an inclination of $27\fdg 5^{+0\fdg 8}_{-0\fdg 6}$, and is located at a distance of $2.2^{+0.2}_{-0.2}\,\mathrm{kpc}$ (Miller-Jones et al., 2021). The X-ray emission is persistent and transitions between two major spectral states, the low/hard and high/soft state (Tananbaum et al., 1972).

Cyg X-1 shows strong X-ray variability over a variety of timescales from milliseconds to years, among which the non-periodic variability known as “shots” resembles second-duration mini-flares (Terrell, 1972). The shots have been observed in the full X-ray band, from $\sim$$0.1~{}\mathrm{keV}$ up to $\sim$$200~{}\mathrm{keV}$ (Yamada et al., 2013; Bhargava et al., 2022), and also in different spectral states (Feng et al., 1999). Similar shots have also been observed in other black hole X-ray binaries, e.g., GX 339–4, GS 2023+338, and GS 1124–68 (Negoro, 1999). Since the X-ray is emitted in the immediate proximity to the black hole, from the inner accretion disk, corona, or jet base, studying the shots can help us further understand the physical processes of the hot plasma around the black hole.

Shots in Cyg X-1 have been intensively studied in the past. Terrell (1972) found that the Cyg X-1 shots were consistent with random pulses that resembled the observed power spectra. While the rough shape of power spectrum can be reproduced by the shot model with appropriate shot parameters (e.g., Miyamoto & Kitamoto, 1989; Belloni & Hasinger, 1990; Lochner et al., 1991), it is difficult to explain other details and the amplitude variations of the power spectrum, e.g., the rms-flux relation (Uttley & McHardy, 2001). The average shot profile was revealed to be slightly asymmetric with a slow rise and a fast decay (Miyamoto & Kitamoto, 1989; Feng et al., 1999; Negoro et al., 2001). The profile can be modeled by the sum of two exponential components, a fast one with a timescale of $\sim$$0.1~{}\mathrm{s}$ and a slow one with a timescale of $\sim$$1~{}\mathrm{s}$ (Negoro et al., 1994; Feng et al., 1999; Negoro et al., 2001). The shot profile was found to be narrower and more asymmetric in the higher energy band, and vary in different spectral states (Feng et al., 1999). Both the shot amplitude and waiting time (time between shots) follow an exponential distribution (Negoro et al., 1995); whether or not there is a deficit of shot occurrence at short time separations is controversial (Negoro et al., 1995; Focke et al., 2005). Negoro & Mineshige (2002) further suggested that the two distributions could be log-normal, implying a physical link between shots and gamma-ray bursts due to similar behaviors (Li & Fenimore, 1996).

The hardness or spectral parameters were found to vary with time within a shot (Negoro et al., 1994; Feng et al., 1999; Yamada et al., 2013; Bhargava et al., 2022). In particular in the hard state, the hardness ratio gradually decreases before the peak but abruptly returns to the average around the peak (Negoro et al., 1994; Feng et al., 1999; Yamada et al., 2013). Below the frequency of a few Hz, the Fourier power and cross spectra, at least the shapes, can be interpreted as due to shot variability (Negoro et al., 2001). A near-unity coherence between different RXTE energy bands was obtained in this frequency range (Nowak et al., 1999; Negoro et al., 2001). Also, a hard time lag was found (Miyamoto et al., 1988; Feng et al., 1999).

Based on these observations, several possible physical mechanisms for the shot have been proposed. The magnetic flare model suggests that the energy release of magnetic reconnections on the accretion disk is responsible for the shots (Galeev et al., 1979; Pudritz & Fahlman, 1982; Poutanen & Fabian, 1999), and is able to explain the prompt spectral hardening (Yamada et al., 2013). The disturbance propagation is another physical model to explain the shot behavior (Manmoto et al., 1996). With numerical simulations, Manmoto et al. (1996) demonstrated that the propagation and reflection of the disturbance wave can naturally reproduce the observed shot profile. When the disturbance propagates toward the innermost region where the radius is comparable to the characteristic scale of the disturbance, the disturbance will evolve from the initial thermal mode into a mixture of thermal and acoustic modes because the plasma is inhomogeneous. While the inward modes are swallowed by the black hole, the acoustic mode can propagate outward, producing the decaying tail of the shot (Kato et al., 1996; Manmoto et al., 1996). Although the disturbance propagation model can successfully reproduce the shot profile, it has difficulties in explaining the prompt hardening around the peak. On the other hand, the magnetic flare model can reproduce the prompt hardening but is difficult to explain the gradual hardness decrease before the peak. Therefore, some authors proposed that both mechanisms might have been involved (Negoro et al., 2001; Yamada et al., 2013).

In the previous studies, shots were detected by techniques based on their relatively high fluxes (e.g., Negoro et al., 1994; Focke et al., 2005; Bhargava et al., 2022). These techniques perform well for strong shots but are less efficient in detecting weak ones. Moreover, they may have difficulties in distinguishing proximate shots and lead to biases in some statistical results (see Focke et al., 2005). In this paper, we propose a new technique that is unbiased for shots with small amplitudes and temporal proximity, enabling us to obtain a more complete and unbiased shot sample. We describe the observations and introduce our new shot detection algorithm in Section 2, including comparisons with previous techniques, and present the results in Section 3. The shot mechanisms are discussed in Section 4.

Neutron Star Interior Composition Explorer (NICER, Gendreau et al., 2012) observations of Cyg X-1 with an effective exposure ${>200~{}\mathrm{s}}$ as of 2023 January 1 are adopted. The level 2 data are downloaded from HEASARC and reprocessed with the nicerl2 tool in HEASoft 6.33.2 with the NICER CALDB version 20240206. We extracted the ${0.25\textendash 12~{}\mathrm{keV}}$ light curves with a time resolution of $0.05~{}\mathrm{s}$ using the tool nicerl3-lc, normalized to 52 FPMs. Cyg X-1 may show dips in its light curve, possibly due to absorption of clumpy winds or clumps in the accretion flow (Feng & Cui, 2002). We identified such dips with a duration of $\sim$$1~{}\mathrm{s}$ in 13 observations¹¹1ObsIDs 1100320101, 1100320102, 1100320103, 1100320104, 1100320106, 1100320107, 1100320119, 1100320120, 1100320121, 1100320122, 2636010102, 5100320129, and 5100320134. and discarded them, as dips with a duration similar to that of the shot may cause problems in shot detection. We also found that when the count rate is high, the shots occur so frequently that the typical waiting time is similar to the shot duration, which may affect the shot extraction. We therefore excluded another 17 observations²²2ObsIDs 0100320102, 0100320110, 1100320108, 1100320109, 1100320118, 2636010201, 5100320120, 5100320124, 5100320125, 5100320126, 5100320127, 5100320128, 5100320130, 5100320131, 5100320132, 5100320139, and 5100320140. and only used those with an average count rate below $12000~{}\mathrm{counts~{}s^{-1}}$. In summary, 30 observations are discarded and 70 observations with a total exposure of 402 ks are used in this work.

We fit the ${3\textendash 12~{}\mathrm{keV}}$ NICER spectra with a single power-law model to determine the spectral state, and found that the photon index is always less than 2.2 in all observations, corresponding to the hard or intermediate state, but not the soft state (Grinberg et al., 2013). Following Grinberg et al. (2013), we also checked the MAXI (Matsuoka et al., 2009) and Swift/BAT (Krimm et al., 2013) data to cross-check the spectral state of Cyg X-1 during the NICER observations, and obtained consistent results.

The new shot detection technique is based on the baseline detection and template fitting. Shots are identified from a baseline-subtracted light curve by fitting with an average shot profile iteratively. The light curve is in the unit of counts.

(a) Baseline subtraction. The baseline $b_{i}$ is subtracted from the original light curve $c_{i}$ to get the new one ${n_{i}=c_{i}-b_{i}}$, in which the long-term variations are removed. We adopted the Asymmetric Least Squares algorithm (asls; Eilers, 2003) for baseline detection. This technique calculates the baseline by minimizing the object function

where $N$ is the number of bins in the light curves, $\Delta^{d}$ is the $d$th order finite-difference operator, and $\lambda$ is a penalty scale factor. The first term in Equation 1 represents the contribution from residuals, and a small $p$ ($<0.5$) tends to generate positive residuals. The second term is associated with the roughness, which requires the baselines to vary smoothly. The baseline is calculated using the Python library pybaselines (Erb, 2022), with ${\lambda=100},{p=0.05}$ and ${d=2}$.

(b) Profile construction. An average shot profile is constructed using strong and isolated shots identified from the baseline-subtracted light curve. A time bin is marked as the potential peak position of a shot if it has the maximum count over the neighboring 21 bins ($1.05~{}\mathrm{s}$). Initially, shots with a peak count greater than a certain threshold $n_{\mathrm{th}}$ are selected and summed to generate the shot profile by aligning the peak position. The shot profile is normalized to have a unity peak amplitude. To exclude superimposed shots, we fit each individual shot with the average profile,

where $n_{i}$ and $p_{i}$ are the counts of the data and normalized shot profile at bin $i$, respectively, $A$ is the baseline-subtracted peak amplitude of the shot, and $\sigma^{2}(c_{i})$ is the variance of the data. Shots with ${\chi^{2}/20<C}$ are used to construct a new profile, where $C$ is an empirical cut for the goodness of fit. This procedure is repeated until the input and output shot samples become identical.

(c) Shot detection. The shots are detected by slide fitting with the shot profile in the light curve one after another. By moving the average shot profile through the whole light curve, one fits the shot amplitude at each time bin and finds the strongest shot that leads to the largest $\Delta\chi^{2}$. A shot component is then added into the model light curve and the search is repeated to find the next strongest one until $\Delta\chi^{2}$ due to the new shot component is lower than a given threshold $\Delta\chi^{2}_{\mathrm{th}}$.

When applying the technique on the NICER data, in each observation, we chose an $n_{\mathrm{th}}$ such that the stronger half of shots are selected. We set $C=5$ in step (b). Different $n_{\mathrm{th}}$ and $C$, which only determine the shot profile, have almost no influence on the scientific results. The stopping criterion is set as $\Delta\chi^{2}_{\mathrm{th}}=70$. Figure 1 shows $\Delta\chi^{2}$ and the shot amplitude $A$ vary as a function of the iteration step. As one can see, $\Delta\chi^{2}$ is roughly scaled with the amplitude of detected shots. We obtained consistent results if $\Delta\chi^{2}_{\mathrm{th}}=60$ is used.

We examined the shot profiles generated in different observations and found that they have a small variation. They have almost the same normalized shape, although the amplitudes may vary. Thus, a combined shot profile from all observations is used, shown in Figure 2. We notice that there is convex (deficit in count rate) on both wings. This is simply because we have used isolated shots to construct the profile. Shots outside the time window raise the outer wings. If we do not require isolated shots for the profile, the composite shot profile decreases monotonically away from the peak, just like those shown in previous studies. In step (c), we only use the central part of the shot profile (see Figure 2), and the convex has no impact on the results. We note that our shot profile shows only a single temporal component with an exponential timescale of 0.07–0.08 s (fitted with an exponential function). It is hard to determine whether or not there exists a second, slow component ($\sim$1 s) as reported in the literature (Negoro et al., 1994; Feng et al., 1999; Negoro et al., 2001), due to shot superpositions.

With this algorithm, we successfully detected shots in the 70 NICER observations. A segment of NICER light curve is displayed in Figure 3 for illustration.

Previous shot detection techniques are mainly based on flux comparison (e.g., Negoro et al., 1994; Feng et al., 1999; Yamada et al., 2013; Bhargava et al., 2022). For example, they required that the shot peak should have the maximum counts among a certain time interval, and also be significantly higher than the local average (Negoro et al., 1994; Yamada et al., 2013; Bhargava et al., 2022). In Figure 3, we also marked shots detected using the technique described in Negoro et al. (1994) and Bhargava et al. (2022). As one can see, the previous techniques cannot effectively detect shots with small amplitudes and those superimposed on each other, e.g., the shots near 39 s and 52 s in Figure 3.

In this paper, we introduced a new algorithm for shot detection. How the detection sensitivity affects the results and how to restore the intrinsic shot amplitude and waiting time is presented in Section 3.1. Our technique allows us to find shots with a time separation down to $\sim$1 s, close to the shot duration. This is consistent with the minimum $\Delta t$ shown in Figure 7, contradicting the conjecture of a deficit of shots below a time separation of $\sim{5\textendash 10~{}\mathrm{s}}$ (Negoro et al., 1995).

As shown in Figure 6 and Figure 7, the distributions of shot amplitude and waiting time are both exponential. The lower cutoff is consistent with the sensitivity of the technique. In particular, the lower cutoff in the $A$ distribution in Figure 6 is scaled with $\bar{c}$. For the waiting time, the exponential slope of the distribution is a strong function of the mean count rate. Therefore, if one plots the distribution in a wide range of count rate, one may see a log-normal distribution as reported in the literature (Negoro et al., 1995; Negoro & Mineshige, 2002).

The mean count rate can be regarded as the mean mass accretion rate. The correlations between the mean count rate and the shot properties may imply that the generation of shots is simply modulated by the accretion rate. At a higher accretion rate, both the average shot amplitude and shot recurrence rate are higher (Figure 4). However, the dependence of the shot amplitude on the waiting time is weak (Figure 5). These mean that the mass accretion rate has a higher impact on the shot recurrence rate than on the shot amplitude. These results can be useful in constraining the shot models.

Negoro et al. (2001) and Yamada et al. (2013) proposed that both magnetic flares and disturbance propagation are needed in the shot generation to account for both the long timescale ($\sim$1 s) flux evolution and short timescale ($\sim$0.1 s or less) hardening, with the disturbance propagation for the longer one and the magnetic reconnection for the shorter. As we demonstrate here that the slightly different flux rising timescales in the two energy bands can lead to a prompt rise of the hardness ratio, there is no need to invoke any short timescale physical processes. An energy dependence of the flux rising timescale is also consistent with the disturbance propagation from the cooler outer disk into the inner hotter disk (Feng et al., 1999). Therefore, the disturbance propagation alone can explain both the shot profile and hardness evolution.