The 2022 Outburst of IGR J17091–3624: Connecting the Exotic GRS 1915+105 to Standard Black Hole X-Ray Binaries

After its last outburst in 2016, IGR J17091–3624 entered a new outburst in 2022 March (Miller et al. 2022). When this outburst began, we triggered our NICER and NuSTAR GO Program (PI: J. Wang). Here, we analyze all 136 NICER observations taken at a near-daily cadence from 2022 March 27 to August 21, as well as six NuSTAR observations taken over this same epoch (see the observation catalog in Table 1). We also requested (by Director's Discretionary Time) one Chandra/HETG observation during this campaign, and this observation took place on June 16, which was simultaneous with the fifth NuSTAR observation. The time evolution of the count rate, fitted disk temperature (see Section 3.1), and the fractional rms are shown in Figure 1.

Zoom In Zoom Out Reset image size

2.2.1. NICER

2.2.2. NuSTAR

2.2.3. Chandra

To decipher the spectral states, we begin by identifying the dominant spectral component in each observation. In an automated way, we fit the flux energy spectra of all the 305 NICER segments of a length of 500 s (i.e., continuous 500 s intervals), making a total exposure time of 152.5 ks. The baseline model used includes the multicolor disk emission (diskbb) and a Comptonization component (nthcomp). We use cflux to calculate the flux contribution from each component. The XSPEC syntax of the model, therefore, is TBabs(cflux∗diskbb + cflux∗nthComp).

The time evolution of the fitted disk temperature is shown in Figure 1. At the beginning of the outburst, the disk temperature is low (kT ∼ 0.5 keV) and rises to 1.5–2 keV as the luminosity increases. We attempt to place these observations in the conventional hardness–intensity diagram (HID) in order to cleanly identify hard and soft states, but because of these very high disk temperatures and the high galactic absorption column (N_H > 10²² cm⁻²), in some observations more thermal-dominated spectra actually led to larger hardness ratios (e.g., see either hardness ratio in the color–color diagram in Figure 2(b)). In other words, the conventional phenomenological HID fails to capture the corona-dominated states versus the thermal-dominated states. To overcome this, we plot the fitted disk temperature as a proxy for the spectral hardness. ¹⁴ The resulting "mimicked" HID is shown in Figure 2(a). With this approach, we can map the evolution of IGR J17091–3624 in its 2022 outburst to more typical BHXBs. Following the classical pattern, IGR J17091–3624 started the outburst in the Hard State (i.e., at low disk temperature, on the bottom right side of the mimicked HID), rose in flux and transitioned to the higher temperatures, corresponding to the HIMS, SIMS, and Soft State, and eventually went back toward the Hard State at a lower flux than the hard-to-soft state transition. This is akin to the hysteresis pattern seen in the HIDs of typical BHXBs.

Zoom In Zoom Out Reset image size

Figure 3 shows representative NICER lightcurves discovered during our observing campaign. The shapes of the lightcurves vary dramatically, and these shapes can be broadly characterized into different states (see also Figure 1 for when each state was observed).

Zoom In Zoom Out Reset image size

At the beginning of the outburst (from March 14 to 17), IGR J17091–3624 showed stochastic variability (Figure 3(a); navy curve) that coincided with the hardest spectra, when the disk temperature was at its lowest, at ∼0.5 keV (see also Miller et al. 2022). On March 18, the source flux began to increase, akin to the spectral HIMS (March 18 to 19) and SIMS (March 22 to 27; Wang et al. 2022c), but the variability remained stochastic (Figure 3(a); cyan and green curves). By March 26, exotic variability started to develop in the lightcurves (Figures 3(b) and (e); termed "At the end of SIMS" and the "Transition to Class V"; Wang et al. 2022a), and by March 28, when the disk temperature was near its highest values, between 1.5 and 2 keV (akin to a soft spectral state), the lightcurves showed very clear, structured variability (Figures 3(f)–(h)). The exotic variability seen in panel (f) is reminiscent of the "Class V" variability identified in Court et al. (2017), while the near-sinusoidal variability seen in panels (g) and (h) does not resemble any previously identified classes. Therefore, we term this new variability "Class X " (more details below). From March 28 to July 18, IGR J17091–3624 made excursions between the variability of Class V, Class X, and a more "traditional" Soft State (Figure 3(c)) that shows very little variability altogether.

Here we describe the exotic variability classes in more detail. The structured, exotic variability began gradually at the end of the SIMS, as sharp flares began arising on top of the stochastic variability (panel (b)). Then, on March 27, the lightcurves started to become more variable, with different characteristics than in the SIMS. The NICER lightcurves from March 26 to March 28 are shown in Figure A1. Comparing the lightcurves at March 27 06:10 and 18:49 UTC, the average NICER count rate decreased from ∼800 to ∼600 counts s⁻¹. Later that day, the average count rate decreased even further, to ∼500counts s⁻¹. Then, within 4.5 hr, the source went from demonstrating largely stochastic variability to showing a distinct and highly structured exotic variability pattern, having firmly transitioned to Class V variability (panel (e)).

Class V lightcurves (panel (f)) are characterized as having repeated, sharp, high-amplitude flares, although the period of those flares drifts even within a 500 s segment. Each major flare's apex can be singly peaked or multipronged, the latter of which we refer to as "mini-flares."

The lightcurves in the new Class X (panels (g)–(h)) show nearly sinusoidal variations. They are distinguished from Class V by the larger amplitudes (the mean rms is 24% compared to 14%), uniformity (see the power spectral density, or PSD, in Section 4.1), and symmetry of the flares. Sometimes there can also be additional mini-flares at the peaks of major flares. In some observations, the Class X appeared to nearly vanish, but then reappeared a few hundred seconds later (panel (h)).

Finally, at the end of our campaign (August 21), we observed the source transition back to lower disk temperatures (∼1 keV), akin to a traditional intermediate state, and also the stochastic variability reemerged (panel (d)). We term this state the Intermediate State Return (IMS Return). We note that while our observations stop on August 21, after this date IGR J17091–3624 was found to exhibit exotic variability once again. The analysis of this later behavior will be published in future work.

From our analysis of the broadband spectral shape (Figure 2), we conclude that this remarkable source transitions between the spectral and timing characteristics of typical BHXBs (e.g., the Hard State, HIMS, SIMS, Soft State, and IMS Return). Then, the shapes of the lightcurves (Figure 3) revealed that there was a phase of Transition to Class V and that sometimes, when the disk component dominates over the corona (as in traditional Soft State), instead of showing very little variability (as in most BHXBs), IGR J17091–3624 can demonstrate exotic (structured and repeated) variability. Therefore, we also identified the exotic variability classes Class V and Class X. In this way, IGR J17091–3624 can be seen as a bridge between the more typical BHXBs and GRS 1915+105, with its famously complex and exotic variability. In the next section, we delve further into the spectral and timing properties of each of these identified states.

As a note to the reader: for the remainder of this paper, we use green/blue/purple colors for observations in the more typical/stochastically varying states (e.g., the Hard State, HIMS, SIMS, Soft State, and IMS Return), while the exotic variability states (Transition to Class V, Class V, and Class X) are shown in red/orange/yellow.

To quantitatively investigate the characteristics of the lightcurves, we compute PSDs of all 305 NICER segments of a length of 500 s. We use the "rms-squared" normalization (Belloni & Hasinger 1990), and the Nyquist frequency is 500 Hz. Representative PSDs from each of our identified states corresponding to the lightcurves in Figure 3 are shown in Figure 4. The PSDs in Figure 4 were computed by averaging ≳10 segments of a length of 500 s and are binned geometrically in frequency, i.e., from frequency ν to (1 + f)ν, where f is called the f factor (see Section 2.2 in Uttley et al. 2014 for more details). We choose an f factor of 0.1 to measure the characteristic frequencies more precisely.

Zoom In Zoom Out Reset image size

As IGR J17091–3624 can evolve very quickly, on a timescale of a day, we also compute the dynamical PSD (DPSD), to show the evolution of the PSDs over time (Figure 5). The DPSD can be regarded as a matrix of the PSD of each 500 s segment. The DPSD has been color-coded by the identified states, where for a given state/color, the darker shade corresponds to higher variability power. The DPSD clearly reveals the exotic variability, with a peak in the power at around 0.1 Hz, but one can also see other characteristic frequencies popping out.

Zoom In Zoom Out Reset image size

To investigate these characteristic frequencies further, we then fit all 305 raw (Poisson noise included) single-segment PSDs with multiple Lorentzian components and a constant for the Poisson noise. The Lorentzians of varying widths describe both the broadband noise and the narrower components (including the "normal" QPOs and the ones caused by exotic variability). We show the characteristic centroid frequencies of the narrower Lorentzian components versus the fractional rms in Figure 6. Below, we describe some of the PSDs in Figure 4 as well as the characteristic timescales for each state and compare to typical BHXBs.

Zoom In Zoom Out Reset image size

In the Hard State (Figure 4(a); navy curve), we detect a QPO at ∼0.3 Hz with a Q factor ∼6 and a fractional rms ∼13%. The QPO is accompanied by a flat-top noise with both low- and high-frequency breaks. The high-frequency break is at a similar frequency to the QPO frequency. These characteristics are consistent with a Type C QPO in a typical hard state.

In the HIMS (Figure 4(a); cyan), both the frequencies of the QPO and the low-frequency break of the flat-top noise increase compared to the Hard State. The QPO is still narrow, with a Q factor ∼6, and its fractional rms is ∼6%. The QPO frequency is in the range of 2.5–5 Hz, and it also anticorrelates with the fractional rms (Figure 6), which is a characteristic of Type C QPOs in normal BHXBs (e.g., see Motta et al. 2011).

In the SIMS (Figure 4(b), green), a narrow and prominent QPO is always present at 2–3 Hz, with a weak power-law noise. No clear correlation between the QPO frequency and the rms can be seen (Figure 6). The Q factor is ∼5 and the rms is ∼5%. These features are consistent with a Type B QPO in the traditional SIMS.

As found in Section 3.2, IGR J17091–3624 transitioned gradually from the SIMS to the exotic Class V, and the PSD in the Transition to Class V has an increase of power generally across frequencies from 0.002 to ∼10 Hz (Figure 4(d)). At the end of the SIMS, because of the consistent Lorentzian centroid frequencies in the segment-based PSDs, we are able to fit the averaged PSD (Figure 4(b)) using the multi-Lorentzian model. We measure the Type B QPO frequency at ν₁ = 2.70 ±0.06 Hz and an additional peak at ${\nu }_{2}={0.016}_{-0.001}^{+0.002}$ Hz. ¹⁵ The additional peak is caused by modulations occurring with a period of $1/{\nu }_{2}={62}_{-8}^{+4}$ s (the corresponding lightcurve is shown in Figure 3(b)).

In Class V (Figure 4(e)), the lack of regularity in the flare period produces a PSD that can be fitted with a broad Lorentzian centered in the range of 0.02–0.2 Hz plus a zero-centered Lorentzian for the broadband noise. The centroid frequency anticorrelates with the rms (Figure 6), meaning that when the rms is higher (the variability amplitude is larger), the characteristic exotic variability timescale is longer.

The distinguishable feature of the new Class X (Figure 4(f)) compared to Class V is the uniformity of the flare timescale, which leads to a narrow peak in the averaged PSD at 0.0154 ± 0.0005 Hz, which does not evolve with the rms. There is also a QPO between 2 and 3 Hz, and in the averaged PSD, its centroid frequency is measured to be ${2.78}_{-0.10}^{+0.07}$ Hz. We note that the ∼0.016 Hz and ∼2.7 Hz features at the end of SIMS match, within 90% uncertainties, the frequencies of the features also seen in Class X, although there is a large difference between their rms and PSD shape. This is interesting, because it might indicate some persistent and intrinsic physical timescale in the system. We will discuss this more in Section 5.1.

In the Soft State (Figure 4(c)), the fractional rms is very low (∼6%), and the corresponding PSD is absent of any component besides a weak power-law noise.

In the data set, we discovered a highly coherent QPO with Q factors (defined as the QPO frequency divided by the full width at half maximum) ≳50. The QPO evolved over time, with its frequency ranging between 5 and 8 Hz, appearing first on April 19 and disappearing on June 26 (see Figure 5). When the QPO was present, the PSD consists of the Poisson noise (the PSD is flat over frequency and is consistent with ≃2/〈x〉, where 〈x〉 is the averaged count rate for the rms-squared normalization), red noise (PSD ∝f⁻²), and an additional noise component that has either a Lorentzian centroid frequency of zero (i.e., flat-top noise) or is in the range of 0.3–0.6 Hz (>3σ away from zero; e.g., see Figure 7). The noise component with nonzero centroid frequency appears as the low-rms extension of the Lorentzian component representing heartbeat-like exotic variability in Class V (see Figure 6). Therefore, for the data with a low total fractional rms ≲6% and a disk-dominated spectrum, we classify them as being in Class V if the centroid frequency for the noise component is nonzero and in the Soft State if the noise component centers at zero. A detailed analysis of the properties and evolution of the highly coherent QPO is presented in a separate paper (Wang et al. 2024).

Zoom In Zoom Out Reset image size

In the IMS Return (Figure 4(a), light blue), the PSD is similar to the initial HIMS and SIMS, with flat-top noise and a QPO in the range of 2–6 Hz (Figure 6).

After systematically analyzing the 305 single-segment NICER flux energy spectra (Section 3.1), we combine the NICER spectra in each of our identified states to perform a more detailed spectral analysis, focusing especially on the iron K band. We also include NuSTAR data to cover a broad energy band and to increase the constraining power of the data. Among the eight states, NuSTAR observations are available in all states besides the initial Hard State, the HIMS, and the Transition to Class V (this last state occurred for only one day; see Table 1). The NuSTAR spectra in observations 1 and 2 are combined as they are both in the SIMS.

First, we model the NICER and NuSTAR flux energy spectra in all eight states with the baseline model TBabs∗crabcorr∗(diskbb+nthComp). The model crabcorr serves a NICER and NuSTAR cross-calibration purpose, multiplying each model by a power law with corrections to both the slope by ΔΓ and the normalization (Steiner et al. 2010). The data-to-model ratios are shown in Figure 8. Below, in this section, in order to test the significance level of the iron emission/absorption line, we add a gaussian line with the normalization free to be positive or negative to the baseline model. The energy, width, and normalization of the Gaussian line and the baseline model are all free to vary.

Zoom In Zoom Out Reset image size

In the states that are akin to the states in typical BHXBs (i.e., the Hard State, HIMS, SIMS, Soft State, and IMS Return), we see a broad iron emission line in the spectrum, a canonical signature of relativistic reflection. Another signature of reflection, the Compton hump, is clearly detected in the IMS Return (when the hard Comptonized component was strongest). In the Soft State, the broad iron emission line is detected at a significance level of 6σ, measured with NICER and NuSTAR spectra combined.

In the "exotic" states (i.e., Class V, Class X, and the Transition to Class V), we detect absorption lines at energies close to the rest energies of Fe xxv (6.7 keV) and Fe xxvi (6.97 keV). The energies, widths, equivalent widths, and significance levels of the significantly detected absorption lines (>3σ) are shown in Table A1. The 90% upper limit on the blueshift is 0.08 keV, corresponding to an outflow velocity <0.01c. We note that during the Transition to Class V, the absorption lines at 6.7 and 6.97 keV are detected at significance levels of 6σ and 3σ with only a 2.2 ks NICER exposure. The Chandra/HETG observation took place when the source was in Class X. The Chandra/HETG unfolded spectrum and the data-to-model ratio using the baseline model are shown in Figure 9. We measure the Fe xxv and Fe xxvi absorption lines at ${6.66}_{-0.04}^{+0.05}$ and 7.00 ± 0.04 keV with both widths <0.08 keV. The equivalent widths are ${8.7}_{-1.2}^{+8.9}$ and ${12}_{-11}^{+7}$eV for Fe xxv and Fe xxvi, respectively. The absorption lines are therefore consistent within 90% uncertainty with NICER.

Zoom In Zoom Out Reset image size

Table A1. Best-fit Parameters when Fitting the Time-averaged Flux Energy Spectra with the Final Model crabcorr∗tbabs∗(diskbb+nthcomp+relxilllpCp+gauss+gauss+gauss)∗gabs

Parameter	Epoch 1	Epoch 2	Epoch 3	Epoch 4	Epoch 5	Epoch 6	Epoch 7	Epoch 8
	Hard State	HIMS	SIMS	Transition to Class V	Class V	Soft State	Class X	IMS Return
					(exotic)		(new exotic)
NuSTAR Obs.			Nu1+2		Nu3	Nu4	Nu5	Nu6
NH	(1.537 ± 0.002) × 1022 cm−2
a*	0.998 (f)
i (degrees)	24 ± 4
AFe	1 (f)
Tin (keV)	0.20 ± 0.02	${0.61}_{-0.02}^{+0.07}$	${1.463}_{-0.004}^{+0.005}$	${1.656}_{-0.020}^{+0.010}$	${1.644}_{-0.003}^{+0.005}$	${1.694}_{-0.016}^{+0.013}$	${1.562}_{-0.004}^{+0.003}$	${1.028}_{-0.008}^{+0.015}$
Ndiskbb	$({1.2}_{-0.6}^{+1.2})\times {10}^{4}$	${250}_{-30}^{+40}$	43.0 ± 0.4	${19.4}_{-0.8}^{+2.2}$	${23.8}_{-1.6}^{+0.3}$	27.6 ± 0.2	34.6 ± 0.2	60 ± 2
Γ	1.60 ± 0.02	2.033 ± 0.013	${2.782}_{-0.015}^{+0.012}$	${3.4}_{-0.6}^{(p)}$	${2.88}_{-0.02}^{+0.33}$	${3.00}_{-0.12}^{+0.08}$	${2.30}_{(p)}^{+0.02}$	${2.35}_{-0.04}^{+0.03}$
kTe (keV)	100 (f)	100 (f)	>182 (p)	100 (f)	${5.0}_{(p)}^{+0.3}$	${5.3}_{-0.1}^{+4.8}$	${10.5}_{-1.3}^{+1.1}$	${59}_{-11}^{+23}$
Nnthcomp	<0.02	${0.14}_{-0.03}^{+0.09}$	0.23 ± 0.01	${0.14}_{-0.04}^{+0.03}$	${0.250}_{-0.011}^{+0.012}$	0.19 ± 0.01	0.120 ± 0.003	0.23 ± 0.02
h (Rg)	${5.0}_{(p)}^{+1.2}$	${9}_{-4(p)}^{+12}$	${170}_{-100}^{+230}$	⋯	⋯	${40}_{-20}^{+10}$	⋯	${5}_{(p)}^{+2}$
Rin (RISCO)	${3.1}_{-0.9}^{+2.6}$	${5.6}_{-0.8}^{+1.5}$	${1.0}_{(p)}^{+3.6}$	⋯	⋯	1.3 ± 0.1	⋯	6.2 ± 1.6
$\mathrm{log}\xi$ (erg·cm⋯−1)	${2.91}_{-0.17}^{+1.13}$	${3.32}_{-0.12}^{+0.21}$	${0.7}_{-0.2}^{+1.2}$	⋯	⋯	${0.84}_{-0.01}^{+0.51}$	⋯	${3.08}_{-0.19}^{+0.21}$
$\mathrm{log}{N}_{{\rm{e}}}$ (cm−3)	18 (f)	18 (f)	${16.8}_{-0.3}^{+0.9}$	⋯	⋯	19.8 ± 0.1	⋯	${15.4}_{-0.4(p)}^{+1.3}$
N relxill (10−3)	${15}_{-4}^{+3}$	9 ± 3	12 ± 2	⋯	⋯	0.08 ± 0.02	⋯	${8.6}_{-1.4}^{+1.0}$
Egauss,1 (keV)	⋯	⋯	⋯	${6.66}_{-0.03}^{+0.02}$	${6.95}_{-0.02}^{+0.04}$	⋯	6.70 ± 0.03	⋯
σgauss,1 (keV)	⋯	⋯	⋯	0.02 (u)	0.02 (u)	⋯	${0.02}_{(p)}^{+0.03}$	⋯
Ngauss,1 (10−4)	⋯	⋯	⋯	${6.4}_{-1.6}^{+1.9}$	1.0 ± 0.3	⋯	1.7 ± 0.5	⋯
EW1 (eV)	⋯	⋯	⋯	${16}_{-2}^{+15}$	${2.0}_{-0.1}^{+2.4}$	⋯	6 ± 3	⋯
Significance level1	⋯	⋯	⋯	6σ	5σ	⋯	5σ	⋯
Egauss,2 (keV)	⋯	⋯	⋯	7.01 ± 0.04	⋯	⋯	6.97 ± 0.04	⋯
σgauss,2 (keV)	⋯	⋯	⋯	0.02 (u)	⋯	⋯	0.02 (u)	⋯
Ngauss,2 (10−4)	⋯	⋯	⋯	4.0 ± 1.7	⋯	⋯	1.2 ± 0.5	⋯
EW2 (eV)	⋯	⋯	⋯	${16}_{-9}^{+10}$	⋯	⋯	5 ± 3	⋯
Significance level2	⋯	⋯	⋯	3σ	⋯	⋯	3σ	⋯
Egauss,3 (keV)	⋯	⋯	⋯	7.78 ± 0.07	⋯	⋯	⋯	⋯
σgauss,3 (keV)	⋯	⋯	⋯	0.05 (u)	⋯	⋯	⋯	⋯
Ngauss,3 (10−4)	⋯	⋯	⋯	${4.2}_{-2.0}^{+1.9}$	⋯	⋯	⋯	⋯
EW3 (eV)	⋯	⋯	⋯	${32}_{-23}^{+6}$	⋯	⋯	⋯	⋯
Significance level3	⋯	⋯	⋯	3σ	⋯	⋯	⋯	⋯
Egabs (keV)	1.74 (f)
σ (keV)	>0.03	${0.05}_{-0.02}^{+0.03}$	0.05 ± 0.03	<0.07	>0.08	0.04 ± 0.03	<0.06	>0.05
Strength	$-{0.007}_{-0.005}^{+0.004}$	−0.005 ± 0.002	−0.005 ± 0.002	−0.004 ± 0.002	${0.002}_{-0.004}^{+0.002}$	$-{0.005}_{-0.002}^{+0.001}$	−0.002 ± 0.001	−0.009 ± 0.002
ΔΓNICER	0 (f)
NNICER	1 (f)
ΔΓFPMA	⋯	⋯	${0.098}_{-0.006}^{+0.004}$	⋯	0.051 ± 0.005	${0.129}_{-0.003}^{+0.006}$	0.139 ± 0.007	0.097 ± 0.004
NFPMA	⋯	⋯	${1.176}_{-0.004}^{+0.010}$	⋯	1.315 ± 0.012	${1.146}_{-0.006}^{+0.012}$	1.314 ± 0.014	1.619 ± 0.011
ΔΓFPMB	⋯	⋯	${0.098}_{-0.006}^{+0.004}$	⋯	0.048 ± 0.006	0.111 ± 0.006	0.138 ± 0.007	0.096 ± 0.004
NFPMB	⋯	⋯	${1.160}_{-0.009}^{+0.005}$	⋯	1.289 ± 0.012	1.100 ± 0.012	1.290 ± 0.014	${1.562}_{-0.015}^{+0.011}$
Disk fraction	${0.042}_{-0.022}^{+0.014}$%	5.0 ± 0.5%	${52.8}_{-1.4}^{+0.8}$%	${61}_{-22}^{+5}$%	${49.8}_{-0.7}^{+0.3}$%	${50.9}_{-0.4}^{+2.9}$%	${65.8}_{-0.6}^{+0.7}$%	${22.2}_{-0.6}^{+0.9}$%
Flux (10−9 erg cm−2 s−1)	${1.51}_{-0.03}^{+0.05}$	3.25 ± 0.02	${5.044}_{-0.007}^{+0.005}$	${3.45}_{-0.03}^{+0.04}$	${4.962}_{-0.004}^{+0.003}$	${5.140}_{-0.006}^{+0.005}$	${4.401}_{-0.006}^{+0.005}$	3.038 ± 0.008
χ2/d.o.f.	4516/3857 = 1.17

Note. Errors are at the 90% confidence level and statistical only. The index (f) means the parameter is fixed, (p) means the uncertainty is pegged at the bound allowed for the parameter, and (u) indicates the parameter is unconstrained. The line widths of the absorption lines modeled by gaussian are limited in the range of 0.02–0.05 keV. The line width is set to have an upper limit of 0.1 keV in the phenomenological gabs model. The flux and disk fraction are both measured using cflux in XSPEC in the 2–20 keV band, and the total flux is unabsorbed flux. The equivalent width is calculated using the XSPEC command eqwidth. Even though the model provides a good fit for the spectra, the exact values of the parameters constrained from reflection models that are constrained mostly with the details in the iron line need to be taken with caveats (e.g., the coronal height, black hole spin, and inclination, etc.), because the assumptions in the reflection model used are not guaranteed to hold.

Download table as:  ASCIITypeset images: 1 2

These results indicate that there is a broad iron emission line from relativistic reflection when there is no exotic variability, and there are absorption lines from highly ionized absorbing material when there is exotic variability. Therefore, in our final model, we add to the baseline model (1) a relativistic reflection model relxilllpCp ¹⁶ (García et al. 2014; Dauser et al. 2022) for the Hard State, HIMS, SIMS, Soft State, and IMS Return; and (2) absorption lines detected above the 3σ level modeled by gaussian for the Transition to Class V, Class V, and Class X. The best-fitted parameter values are shown in Table A1, and the data-to-model ratios and the unfolded spectra are shown in Figures B1 and 10, respectively. We note that IGR J17091–3624 is a peculiar BHXB with a sometimes very hot disk (T_in can reach ≳1.5 keV), so the assumptions in the reflection model that we use, such as the low-energy break of the Comptonization component, and the geometrically thin disk assumed are not guaranteed to hold (see Appendix B for more details). Therefore, the exact values of the parameters constrained by the reflection model need to be taken with caveats.

Zoom In Zoom Out Reset image size

While the parameters constrained from reflection modeling (e.g., coronal height, spin, and inclination, etc.) warrant caution, they do provide a good description of the reflection features, and therefore the continuum modeling is robust, especially when NuSTAR data are included. We show the parameters describing the properties of the disk blackbody and the corona in each state at different X-ray fluxes in Figure 11. As found in the single-segment NICER spectral fits (Figure 1), the disk temperature is ≲1 keV in the Hard State, HIMS, and IMS Return, and it varies in the range of 1.5–2 keV in the other five states, when the disk fraction is >50%. These five states also have a soft coronal spectrum with a high photon index ≳2, which further solidifies the identification of the SIMS and Soft State, if mapped to the traditional spectral states of BHXBs. We note that the electron temperature of the corona is very low in Class V, Class X, and the Soft State, around 5–10 keV. This has been found also in GRS 1915+105 (e.g., Neilsen et al. 2011).

Zoom In Zoom Out Reset image size

The global parameters that should not change in time are tied between epochs and, as such, they are constrained with more confidence than the parameters from reflection spectroscopy from individual states. The column density of galactic absorption is constrained to be N_H = (1.537 ± 0.002) × 10²² cm⁻², consistent with previous measurements (e.g., Xu et al. 2017; Wang et al. 2018). The fitted inclination angle via reflection spectroscopy is ${\left(24\pm 4\right)}^{\circ }$. Previous measurements from reflection spectroscopy using NuSTAR data in the hard state resulted in $i={({37}_{-4}^{+3})}^{\circ }$ and ${\left(45.3\pm 0.7\right)}^{\circ }$ (Xu et al. 2017; Wang et al. 2018), also suggesting a low inclination for the inner disk producing the reflection.

In Section 3, we describe the broadband continuum shape (i.e., whether it is corona- or thermal-dominated) and the shapes of the lightcurves (i.e., whether they show exotic variability) in order to classify the complex phenomenology of IGR J17091–3624 into eight states. Some of those states are akin to the states in traditional BHXBs (namely, the Hard State, HIMS, SIMS, Soft State, and IMS Return), and then when the source is thermal disk-dominated, it can sometimes take excursions into states with very complex and exotic variability (namely, the Transition to Class V, Class V, and Class X). Here, we summarize the key properties of each of these states, and particularly describe the emission-/absorption-line structure and the PSD structure of each state. The states listed below are given in the order in which they appeared for the first time in this outburst (see Table 1). The quoted measured quantities using flux energy spectra can be found in Table A1.

• 1.  
  Hard State. As in typical BHXBs, the variability in this state is stochastic, with a high averaged fractional rms of 27%, and the PSD consists of a Type C QPO at ∼0.3 Hz on top of flat-top noise. The flux energy spectrum contains a cool disk with disk temperature T_in = 0.20 ±0.02 keV, a corona with a hard spectrum (photon index Γ = 1.60 ± 0.02), and a broad iron emission line due to relativistic reflection.
• 2.  
  H IMS. The variability is still stochastic, while the fractional rms has decreased to 12%. Both the frequencies of the Type C QPO and the low-frequency break of the flat-top noise increase compared to the Hard State. The QPO frequency is between 2.5 and 5 Hz. The disk temperature increases to ${0.61}_{-0.02}^{+0.07}$ keV and the coronal spectrum is softer than it was, with Γ = 2.033 ± 0.013. The broad iron emission line is present from reflection.
• 3.  
  S IMS. The fractional rms further declines to 8% and a Type B QPO at 2–3 Hz is present. At the end of the SIMS (March 26), the lightcurve starts to show flares and modulations that are signs of emerging exotic variability on a timescale of ${62}_{-8}^{+4}$ s, corresponding to a narrow peak at ${0.016}_{-0.001}^{+0.002}$ Hz in the PSD. The disk becomes hotter than it was, with a disk temperature of ${1.463}_{-0.004}^{+0.005}$ keV, and the coronal spectrum further softens to ${\rm{\Gamma }}={2.782}_{-0.015}^{+0.012}$. The broad iron emission line is detected in both NICER and NuSTAR spectra.
• 4.  
  Transition to Class V. The lightcurve shows that Class V exotic variability developed over 4 hr in this transitional state, while the flux decreased to half that of the SIMS peak (see Figure 2(a)). Due to this exotic variability, the fractional rms increases to 16%, and the power increases generally across frequencies from 0.002 to ∼10 Hz. The disk is still very hot, with ${T}_{\mathrm{in}}={1.656}_{-0.020}^{+0.010}$ keV, and the coronal spectrum is soft, with Γ > 2.8. Strong absorption lines close to the rest energies of Fe xxv and Fe xxvi are detected at the 6σ and 3σ confidence levels.
• 5.  
  Class V. Exotic flaring variability is evident in the lightcurves. The PSD contains a broad component with centroid frequency in the range of 0.02–0.5 Hz, resulting from the irregularity of the exotic variability. The irregular variability pattern and the broad component in the PSD are similar to Class V in Court et al. (2017), corresponding to Class μ in GRS 1915+105 (Belloni et al. 2000). We will discuss this further in Section 5.1. The averaged fractional rms is 12%. The disk temperature is ${1.644}_{-0.003}^{+0.005}$ keV, with a high disk fraction of ${49.8}_{-0.7}^{+0.3}$%, and the photon index is ${2.88}_{-0.02}^{+0.33}$. An iron xxvi absorption line is detected at 5σ confidence.
• 6.  
  Soft State. No exotic variability can be seen in the lightcurves and the averaged fractional rms is only 5%, the lowest among the states. Starting from April 19, a highly coherent QPO appeared (Wang et al. 2024). The disk temperature is ${1.694}_{-0.016}^{+0.013}$ keV and the disk fraction is ${50.9}_{-0.4}^{+2.9}$%. The coronal spectrum is soft, with ${\rm{\Gamma }}={3.00}_{-0.12}^{+0.08}$, and an iron emission line is detected at 6σ significance. The low rms, hot-disk domination, and soft coronal spectrum are the features that are in agreement with a traditional soft state (see also the discussion in Section 5.1).
• 7.  
  Class X. Exotic, large-amplitude, near-sinusoidal variability is prominent in the lightcurves, and the fractional rms increases to 24%. The variability amplitude can change within 500 s (see Figures 3(g)–(h)), but the uniformity of the flare timescale leads to a narrow peak in the averaged PSD at 0.0154 ± 0.0005 Hz. This class has never been seen before, in either this source or GRS 1915+105, and we will discuss it more in Section 5.1. With a disk temperature of ${1.562}_{-0.004}^{+0.003}$ keV, the disk fraction is ${65.8}_{-0.6}^{+0.7}$%, the highest among all the states. Both Fe xxv and xxvi absorption lines are detected in NICER and Chandra/HETG spectra.
• 8.  
  I MS Return. The variability is stochastic and the PSD shape is similar to the initial HIMS. This state connects the initial HIMS and SIMS in the mimicked HID and color–color diagram (Figure 2). Compared to previous states, the disk temperature has dropped to ${1.028}_{-0.008}^{+0.015}$ keV, along with a lower disk fraction of ${22.2}_{-0.6}^{+0.9}$%. The coronal spectrum is still soft, with ${\rm{\Gamma }}={2.35}_{-0.04}^{+0.03}$. Both reflection signatures, the broad iron line and the Compton hump, are prominent in the NICER and NuSTAR spectra.

Previously, nine variability classes were defined for IGR J17091–3624 in its 2011–2013 outburst (Altamirano et al. 2011; Court et al. 2017). In this section, we will compare the variability classes we identify for this 2022 outburst with those in previous works.

In the 2022 outburst, we observed Class V identified based on the repeated, sharp, high-amplitude flares in the lightcurves and the PSD shape. We note that a QPO at ∼4 Hz with a Q factor of ∼3 was observed in Class V in the 2011 outburst (Figure 11 in Court et al. 2017). While Class V variability was observed previously, here we see for the first time that the timescale of variability can evolve dramatically, and that it scales with the fractional rms of the source (Figure 6).

We define a new Class X, because of several properties that are different from the previously identified nine classes. Although the uniformity of the exotic variability timescale is reminiscent of the heartbeat class (Class IV in IGR J17091–3624 and Class ρ in GRS 1915+105), the flare patterns are more symmetric than the typical heartbeat (slow rise and quick decay). In the PSD, there is sometimes a QPO at 2–3 Hz, while in Class IV, no QPO was previously detected (Court et al. 2017) or the QPO was at 6–10 Hz (Altamirano et al. 2011). We also find that in Class X, the exotic variability amplitude can vary within even 500 s, while maintaining the variability timescale (see Figures 3(g) and (h)). When the variability amplitude is relatively small, the variability pattern is not only similar to Class III (Class V in GRS 1915+105), but also similar to the end of SIMS, when the source first started to show flares and modulations. The modulation produces a peak in the PSD at ∼0.016 Hz, in addition to a Type B QPO at ∼2.7 Hz. The two characteristic frequencies are respectively consistent within 90% uncertainty at the end of SIMS and in Class X (see Section 4.1 and Figure 6). This suggests that Class X can be regarded as a high-rms extension of the structured variability that starts at the end of SIMS.

We identify the Soft State based upon the lack of exotic variability, the lowest fractional rms of 6%, the hot-disk domination, and a very soft coronal spectrum. There are two classes identified in Court et al. (2017) that also show no exotic variability in the lightcurves—Class I and Class II. The PSD in Class I is similar to the intermediate state here, with a broadband noise at 1–10 Hz and a QPO at ∼5 Hz. We notice the PSD in our Soft State in the NICER hard band 4.8–9.6 keV is similar to the one in Class II using RXTE data (2–60 keV), as both lack any power above ∼1 Hz. However, in our Soft State, the PSD sometimes shows a highly coherent QPO at 5–8 Hz (Wang et al. 2024), which was not detected in Class II.

In this work, we find that in the 2022 outburst, IGR J17091–3624 went through a traditional hard state and intermediate state, and then entered an exotic soft state, where it sometimes exhibited heartbeat-like variability in the lightcurves. This transition from traditional BHXB states to states showing exotic variability was also suggested in its previous two outbursts in 2011 and 2016 (Pahari et al. 2014; Wang et al. 2018).

In a seminal work, for GRS 1915+105, Belloni et al. (2000) identified 12 variability classes and three basic states and suggested that the variability classes are produced by transitions between the three basic states with certain patterns. Although they found some similarities between the three basic states and the canonical BHXB accretion states, including the soft state and intermediate state, in the variability class GRS 1915+105 spends most of its time in (Class X), the spectrum is the hardest and the lightcurve shows stochastic variability, although it is still not as hard as the canonical hard state of BHXBs. Therefore, the accretion state landscape in GRS 1915+105 is significantly more challenging to understand, due to the more complex phenomenology.

In summary, IGR J17091–3624 shares some variability classes with GRS 1915+105 that are not seen in other BHXBs, while having an outburst recurrence rate (see the summary in Section 1) and evolution pattern in outbursts similar to normal BHXBs. It can therefore be regarded as a bridge between the most peculiar BHXB GRS 1915+105 and normal BHXBs.

5.3.1. A Brief Review of Radiation Pressure Instability

Although it is generally accepted that the heartbeat-like variability is due to some disk instability, the nature of that instability is unclear. In the standard Shakura–Sunyaev disk, there are two types of instability that might lead to limit-cycle oscillations in the lightcurves: hydrogen ionization instability and radiation pressure instability (e.g., Done et al. 2007; Janiuk & Czerny 2011). The ionization instability takes place in the outer disk and may explain the outburst and quiescence behavior of BHXBs, i.e., lightcurve variability on a long timescale (hundreds of days). On the other hand, radiation pressure instability sets in at a relatively higher mass accretion rate, and thus the inner disk becomes radiation-pressure-dominated. This could then lead to thermal-viscous limit cycles, and is therefore a widely accepted physical mechanism for driving the heartbeat-like variability (Belloni et al. 1997; Janiuk et al. 2000; Nayakshin et al. 2000; Neilsen et al. 2011).

One way to understand the radiation pressure instability is through the "S-curve" when plotting the mass accretion rate $\dot{m}$ versus the disk surface density Σ at a certain disk radius (e.g., Done et al. 2007). From bottom to top, the S-curve consists of three branches: a stable branch when heating is proportional to the gas pressure balanced by radiative cooling, an unstable branch when radiation pressure dominates (the "Lightman–Eardley instability"; Lightman & Eardley 1974), and another stable branch when advective cooling is effective to balance the heating (the "slim-disk" solution; Abramowicz et al. 1988). That is to say, when the inner disk is radiation-pressure-dominated and $\dot{m}$ is on the middle unstable branch, a limit-cycle instability is expected. Over each limit cycle, the mass accretion rate oscillates at each disk radius (switching between the two stable branches), resulting in a "density wave." Observational pieces of evidence supporting this hypothesis include the correlated limit-cycle timescale, disk inner radius (Belloni et al. 1997), and extensive phase-resolved spectroscopic analysis of the heartbeats (e.g., Neilsen et al. 2011; Zoghbi et al. 2016).

We observe a persistent heartbeat timescale of ∼60 s at the end of SIMS and in Class X (see Section 4.1 and Figure 6). A back-of-envelope estimate of the limit-cycle duration is therefore the viscous timescale for the accretion disk (Belloni et al. 1997), which is given by

$$\begin{eqnarray*}&&{t}_{\mathrm{visc}}=60\,{\rm{s}}\left(\displaystyle \frac{M}{10{M}_{\odot }}\right){\left(\displaystyle \frac{R/{R}_{{\rm{g}}}}{46}\right)}^{3/2}{\left(\displaystyle \frac{\alpha }{0.03}\right)}^{-1}{\left(\displaystyle \frac{H/R}{0.1}\right)}^{-2},\end{eqnarray*}$$

where M is the mass of the black hole, R is the radius in the disk, R_g is the gravitational radius GM/c², α is the viscosity parameter, and H is the vertical scale height of the disk (e.g., Frank et al. 2002). A viscous timescale of 60 s can be explained, e.g., by the parameter combination including a black hole mass of 10M_⊙, α = 0.03, H/R = 0.1, and R = 46R_g; or, if the black hole mass is 3M_⊙, R = 102R_g.

5.3.2. Did IGR J17091–3624 Reach the Eddington Limit?

The Eddington ratio threshold for radiation pressure instability to occur is very uncertain (it could range from 6% L_Edd to near Eddington), fundamentally because the S-curve depends on a variety of physical assumptions and conditions (e.g., Honma et al. 1991; Janiuk et al. 2000, 2002). It is therefore natural to find an observational threshold of the Eddington ratio where heartbeats are seen to constrain the underlying accretion disk physics.

So far, heartbeat-like (repetitive, high-amplitude, and structured) variability has been observed in both accreting black hole and neutron star systems and can be a multiwavelength phenomenon. (1) Two BHXBs, including GRS 1915+105 and IGR J17091–3624. In GRS 1915+105, besides X-ray, heartbeat-like variability was also seen in radio and infrared (e.g., Fender et al. 1997; Pooley & Fender 1997). As the mass and distance of GRS 1915+105 are known, it reaches 80%–90% of the Eddington luminosity (Neilsen et al. 2011). (2) Three neutron star low-mass X-ray binaries, including GRO J1744–28 (known as the "Bursting Pulsar"; Kouveliotou et al. 1996; Degenaar et al. 2014), MXB 1730–335 (the "Rapid Burster"; Bagnoli & in't Zand2015), and Swift J1858.6–0814 (Vincentelli et al. 2023). The heartbeat-like (GRS 1915+105-like) variability in neutron stars is also sometimes referred to as Type II X-ray bursts. The luminosity reached near Eddington, ≲20%, and 40% L_Edd for the three sources, respectively. We note that the criterion for radiation pressure instability to occur in neutron star systems is when the magnetospheric radius is larger than the neutron star radius (Mönkkönen et al. 2019). This means that the Eddington ratio is not the only key parameter, but also the magnetic field strength. It has been shown that all three neutron star X-ray binaries satisfy this condition for heartbeat-like variability to take place (Vincentelli et al. 2023). (3) An ultraluminous X-ray source in NGC 3621 with an unclear compact object nature (Motta et al. 2020).

The challenging part of finding an observational Eddington ratio threshold then arises—we only have two BHXBs exhibiting heartbeats and for IGR J17091–3624, we do not know the black hole mass (from the mass function) or distance (e.g., from parallax) to estimate the Eddington ratio. The position of IGR J17091–3624 on the radio versus X-ray luminosity diagram suggests a distance between ∼11 and ∼17 kpc for it to lie on the track followed by other BHXBs (Rodriguez et al. 2011). However, it is also possible for IGR J17091–3624 to not follow the typical relation, as was found for GRS 1915+105. A black hole mass in the range of 8.7–15.6 M_⊙ is obtained under the two-component advective flow framework (Iyer et al. 2015). On the other hand, there are several BHXBs reaching close to the Eddington limit that do not exhibit heartbeat-like variability (e.g., XTE J1550–564 and 4U 1543–47; Rodriguez et al. 2003; Negoro et al. 2021). Therefore, before the heartbeats in IGR J17091–3624 were discovered, it was proposed that GRS 1915+105 was a unique BHXB that showed heartbeats, because it was the only one that had spent any considerable time above the Eddington limit (Done et al. 2004). Could this also be the case for IGR J17091–3624?

In Figure 12, we show the disk flux F_d versus the disk temperature T_in in the segment-based fits of the broadband continuum (see Section 3.1) in three different states. We also perform fits with the model F_d ∝ Tⁿ (i.e., the disk luminosity L ∝ Tⁿ ), where the power-law index n could indicate the nature of the accretion flow. The standard Shakura–Sunyaev disk predicts L ∝ T⁴, and the advection-dominated accretion flows including the "slim disk" follow L ∝ T² (Watarai et al. 2000). We find that in both variability classes showing heartbeat-like variability, the disk luminosity–temperature (L-T) relation is consistent with that of a thin disk (n = 3.7 ± 0.2 and 4.0 ± 0.3 in Class V and Class X). The indication that heartbeat-like variability is associated with a thin disk is consistent with previous conclusions from a theoretical point of view (Nayakshin et al. 2000). We also notice previous phase-resolved spectroscopy of GRS 1915+105 found that the disk forms a loop in the L-T diagram over heartbeat cycles (Neilsen et al. 2011). On the other hand, in the SIMS, the disk has a flatter L-T relation, with n = 1.4 ± 0.3, more in line with a slim disk that is expected at a relatively high Eddington ratio (≳30%; Abramowicz et al. 2010). It is also consistent with n = 4/3 when the mass accretion rate is constant and the disk inner radius is variable, as expected when the local Eddington limit is reached (Lin et al. 2009). In either interpretation of the flatter L-T relation, a relatively high mass accretion rate (≳30% of the Eddington rate) is suggested.

Zoom In Zoom Out Reset image size

In our observing campaign, we caught the precursor of the heartbeats: the end of the SIMS and Transition to Class V, during which the flux dropped gradually from the SIMS peak to half of that. Considering the suggested slim-disk nature in the SIMS and the idea of the "S-curve," we have a plausible explanation for this precursor behavior. In the SIMS, the inner disk lies on the stable upper branch in the S-curve, which corresponds to the slim-disk solution. Then the mass accretion rate decreased, the inner disk thus entered the middle unstable branch, and heartbeat-like variability (but the variability amplitude is relatively small at the end of the SIMS and is irregular in the Transition to Class V) started to show. Then, in Class V and Class X, the inner disk started to exhibit limit-cycle variabilities, switching between the upper and lower stable branches. The heartbeat-like variability is more structured and regular compared to the precursor phase. The flux in these two classes is also in the middle of the SIMS peak and the Transition to Class V.

The flux peak in the SIMS is near the flux peak in the entire outburst. The measured unabsorbed flux (in units of 10⁻⁹ erg cm⁻² s⁻¹) at the SIMS peak is ∼8.5 (0.1–200 keV; the flux shown in Figure 2(a) is in 1–10 keV). Assuming that IGR J17091–3624 at the SIMS peak emits at different Eddington ratios, the black hole mass and distance are shown in Figure 13. Based on the coordinates of IGR J17091–3624 (galactic longitude l = 34952 and latitude b = 221), assuming the distance to the galactic center is 8 kpc, and the radius of the "edge" of the stellar disk in our galaxy is 10–15 kpc (Bland-Hawthorn & Gerhard 2016), the upper limit of the distance to IGR J17091–3624 if it is within the stellar disk is ∼23 kpc. Then, in the case of emitting at 100% Eddington, the black hole mass needs to be <4.3 M_⊙. This would make IGR J17091–3624 one of the least massive black holes known. This limit could be brought up to <8.6 M_⊙ when adopting a bolometric correction factor of 2.

Zoom In Zoom Out Reset image size

Another empirical benchmark of the Eddington ratio is that the luminosity during the soft-to-hard state transition has a mean value of 2–3% L_Edd (Maccarone 2003; Dunn et al. 2010; Tetarenko et al. 2016). As the SIMS peak has ∼two times the flux of the soft-to-hard transition, assuming a slim disk in the SIMS is expected to be ≳30% L_Edd, the soft-to-hard transition in this outburst was at ≳15% L_Edd, which is ∼1.8σ above the mean value using the standard deviation in Tetarenko et al. (2016).

In summary, our data suggest that IGR J17091–3624 in its 2022 outburst reached a relatively high Eddington ratio because of the L-T relation in the SIMS. Depending on the accretion disk model, this can be, e.g., ≳30% L_Edd (Abramowicz et al. 2010), but depends on several parameters, such as the magnetization (e.g., Marcel et al. 2022). We cannot rule out the possibility that IGR J17091–3624 reached a super-Eddington luminosity based on the available black hole mass and distance constraints. Revisiting the conclusion in Done et al. (2004), it is therefore still possible that a super-Eddington luminosity is a necessary condition for heartbeat-like variability to be present. If super-Eddington accretion is a necessary condition, we see new evidence that it alone is not a sufficient condition. In the mimicked HID (Figure 2(a)), the Soft State reaches a similar X-ray flux as the two exotic variability classes (Class V and Class X). This might suggest that it is not just the mass accretion rate that determines whether exotic variability occurs, but rather another factor is required. This other factor could perhaps be the timescale for the inner accretion disk to respond to mass accretion rate changes in the outer disk. We also cannot rule out the possibility that in the Soft State, heartbeat variability is present at other wavelengths rather than X-ray. For example, in Swift J1858.6–0814, the neutron star system where heartbeat-like variability was discovered recently, the heartbeat is most manifest in the infrared band and is much less discernible in X-ray. On the other hand, if IGR J17091–3624 reached only sub-Eddington, it is also a mystery what factor sets IGR J17091–3624 apart from other BHXBs to exhibit heartbeat-like variabilities.

5.3.3. Beyond the Radiation Pressure Instability Model

5.3.4. Interplay between Heartbeats and Iron Emission/Absorption Lines

5.3.5. Future Directions