An Atlas of Exotic Variability in IGR J17091-3624: A Comparison with GRS 1915+105

TL;DR

This study analyzes RXTE data of IGR J17091-3624's 2011-2013 outburst, identifying variability classes, comparing them with GRS 1915+105, and exploring the source's properties and evolution.

Contribution

It introduces a new set of variability classes for IGR J17091-3624 and compares them with GRS 1915+105, challenging the idea that Eddington-limited accretion is necessary for certain variability.

Findings

Identified 9 variability classes in IGR J17091-3624.

First detection of IGR J17091-3624 above 150 keV.

Constraints on mass, distance, and accretion rate of the source.

Abstract

We performed an analysis of all RXTE observations of the Low Mass X-ray Binary and Black Hole Candidate IGR J17091-3624 during the 2011-2013 outburst of the source. By creating lightcurves, hardness-intensity diagrams and power density spectra of each observation, we have created a set of 9 variability `classes' that phenomenologically describe the range of types of variability seen in this object. We compare our set of variability classes to those established by Belloni et al. (2000) to describe the similar behaviour of the LMXB GRS 1915+105, finding that some types of variability seen in IGR J17091-3624 are not represented in data of GRS 1915+105. We also use all available X-ray data of the 2011-2013 outburst of IGR J17091-3624 to analyse its long-term evolution, presenting the first detection of IGR J17091-3624 above 150 keV as well as noting the presence of `re-flares' during the latter stages of the outburst. Using our results we place new constraints on the mass and distance of the object, and find that it accretes at <33% of its Eddington limit. As such, we conclude that Eddington-limited accretion can no longer be considered a sufficient or necessary criterion for GRS 1915+105-like variability to occur in Low Mass X-Ray Binaries.

Figures

Tables

Table 1. Table 1: Chandra observations log covering the three observations considered in this paper. I 𝐼 I refers to Imaging mode, C 𝐶 C refers to CC33_Graded mode and T 𝑇 T refers to Timed Exposure Faint mode. HETG refers to the High Energy Transmission Grating.

OBSID	Instrument	Grating	Exposure (ks)	Mode	MJD
12505	HRC-I	NONE	1.13	$I$	55626
12405	ACIS-S	HETG	31.21	$C$	55774
12406	ACIS-S	HETG	27.29	$T$	55840

Table 2. Table 2: Lower and upper quartile count rates, fractional RMS and median colour averaged across all observations belonging to each class. Count rates and fractional RMS are taken from the full energy range of RXTE/PCA , and fractional RMS values are 2–60 keV taken from lightcurves binned to 0.5 s. Count rates are normalised for the number of PCUs active during each observation. All values are quoted as 1 ​ σ 1 𝜎 1\sigma ranges.

Class	LQ Rate	UQ Rate	Frac. RMS	Median C1
	(cts s-1)	(cts s-1)
I	84–108	106–132	0.13–0.19	0.4–0.68
II	43–57	59–71	0.15–0.23	0.4–0.68
III	64–84	80–110	0.17-0.23	0.35–0.45
IV	63–81	92–122	0.27–0.37	0.32–0.4
V	49–67	88–134	0.44–0.54	0.28–0.46
VI	64–98	111–155	0.29–0.47	0.33–0.61
VII	65–79	128–140	0.45–0.57	0.32–0.42
VIII	62–88	142–178	0.42–0.52	0.36–0.49
IX	87–111	114–144	0.16–0.24	0.42-0.6

Table 3. Table 3: A tally of the number of times we assigned each of our nine Variability Classes to an RXTE  orbit. We have also calculated the amount of observation time corresponding to each class, and thus inferred the fraction of the time that IGR J17091 spent in each class. Note: the values in the Total Time column assume that each orbit only corresponds to a single variability Class.

Class	Orbits	Total Time (s)	Fraction
I	31	69569	14.8%
II	26	50875	10.8%
III	14	26228	5.6%
IV	31	69926	14.9%
V	35	72044	15.3%
VI	29	54171	11.5%
VII	11	19241	4.1%
VIII	16	26553	5.7%
IX	50	81037	17.3%

Table 4. Table 4: Results from the IBIS/ISGRI analysis of the 2011–2013 Outburst of IGR J17091. The 20–40 keV flux is given in units of mCrab and (10 -11  ergs s -1  cm -2 ). Conversion between counts and mCrab was obtained using an observation of the Crab taken during Revolution 1597 between MJD 57305.334 and 57305.894 and the conversion factors of Bird et al. ( 2016 ) and Bazzano et al. ( 2006 ) .

Energy	Intensity	Significance	Exposure	Flux	Flux
(keV)	(cts/s)	$\sigma$	(ks)	(mCrab)	(10-10ergs s-1 cm-2)
20–40	12.39$\pm$0.05	247	115	93.5$\pm$0.38	7.08$\pm$0.03
40–100	7.06$\pm$0.05	157	163	83.5$\pm$0.60	7.87$\pm$0.06
100–150	1.05$\pm$0.03	40	173	66.9$\pm$1.91	2.14$\pm$0.06
150–300	0.23$\pm$0.03	7.6	179	46.6$\pm$5.96	2.24$\pm$0.29

Table 5. Table 5: The nine variability classes of IGR J17091-3624, showing the name of the closest corresponding variability class in GRS 1915+105. The names of GRS 1915+105 classes are taken from Belloni et al. ( 2000 ) , where more detailed descriptions can be found. Eight additional classes of GRS 1915+105 have been described; we do not find analogies to these classes in IGR J17091-3624.

IGR J17091-3624 Class	GRS 1915+105 Class
I	$\chi$
II	$\varphi$
III	$\nu$
IV	$\rho$
V	$\mu$
VI	$\lambda$
VII	None
VIII	None
IX	$\gamma$

Table 6. Table 6: The six OBSIDs explicitly classified in Altamirano et al. ( 2011a ) . We also present the GRS 1915 class with which we implicitly label each OBSID in this paper.

OBSID	Altamirano et al. Class	Court et al. Class
		(implied)
96420-01-04-03	$\alpha$	$\rho /\nu$
96420-01-05-00	$\nu$	$\rho /\nu$
96420-01-06-00	$\rho$	$\rho /\nu$
96420-01-07-01	$\rho$	$\mu$
96420-01-08-03	$\beta /\lambda$	$\lambda$
96420-01-09-06	$\mu$	$\lambda$

Table 7. Table 7: Here is listed the Observation IDs for every R ​ X ​ T ​ E 𝑅 𝑋 𝑇 𝐸 RXTE observation that was used in this analysis, along with the variability class which has been assigned to it. Orb. is the orbit ID (starting at 0) of each observation segment, Exp. is the exposure time in seconds and X is the prefix 96420-01. This table is continued overleaf in Table 8 .

MJD	OBSID	Orb.	Class	Exp.	MJD	OBSID	Orb.	Class	Exp.	MJD	OBSID	Orb.	Class	Exp.
55622	X-01-00	0	I	1840	55676	X-09-06	0	VI	3540	55741	X-18-05	0	VII	782
55622	X-01-000	0	I	3480	55677	X-09-01	0	V	1676	55743	X-19-00	0	VII	1412
55622	X-01-000	1	I	1656	55678	X-09-04	0	V	2090	55744	X-19-01	0	VIII	1938
55622	X-01-000	2	I	3384	55679	X-09-02	0	V	2306	55745	X-19-02	0	VII	2172
55622	X-01-000	3	I	3400	55680	X-10-02	0	V	952	55747	X-19-03	0	VIII	1691
55622	X-01-000	4	I	3384	55681	X-10-00	0	V	3725	55748	X-19-04	0	VI	1283
55623	X-01-01	0	I	1240	55682	X-10-03	0	V	1157	55749	X-19-05	0	VIII	1417
55623	X-01-01	1	I	752	55684	X-10-01	0	III	1504	55751	X-20-05	0	VI	1726
55623	X-01-01	2	I	992	55686	X-10-04	0	III	1127	55752	X-20-01	0	VIII	1079
55623	X-01-01	3	I	1184	55686	X-10-05	0	II	2179	55753	X-20-02	0	VIII	1433
55623	X-01-01	4	I	1056	55687	X-11-00	0	II	3537	55754	X-20-03	0	VII	1122
55623	X-01-010	0	I	2080	55688	X-11-01	0	II	1153	55756	X-20-04	0	VIII	1486
55623	X-01-010	1	I	1832	55690	X-11-02	0	II	1408	55757	X-21-00	0	VIII	3372
55623	X-01-010	2	I	1648	55691	X-11-03	0	II	886	55758	X-21-01	0	VIII	3383
55623	X-01-010	4	I	1424	55692	X-11-04	0	II	3566	55759	X-21-02	0	VI	1938
55623	X-01-010	5	I	400	55693	X-11-05	0	II	1817	55761	X-21-04	0	VII	1497
55623	X-01-02	0	I	3056	55694	X-12-00	0	II	2761	55762	X-21-05	0	VII	1548
55623	X-01-02	1	I	2792	55695	X-12-01	0	II	1374	55763	X-21-06	0	VII	2202
55623	X-01-02	2	I	2432	55695	X-12-02	0	II	2041	55764	X-22-00	0	VII	1682
55623	X-01-020	0	I	3456	55696	X-12-03	0	II	1456	55765	X-22-01	0	VII	1221
55623	X-01-020	1	I	3464	55698	X-12-04	0	II	1916	55766	X-22-02	0	V	720
55623	X-01-020	2	I	3512	55698	X-12-05	0	II	3139	55767	X-22-03	0	V	1801
55623	X-01-020	3	I	3520	55700	X-12-06	0	II	1189	55768	X-22-04	0	VIII	1983
55623	X-01-020	4	I	3512	55701	X-13-00	0	II	1214	55769	X-22-05	0	VIII	999
55623	X-01-020	5	I	464	55702	X-13-01	0	II	980	55770	X-22-06	0	VIII	667
55624	X-02-00	0	I	1758	55704	X-13-02	0	II	732	55771	X-23-00	0	VIII	2075
55626	X-02-01	0	I	1380	55705	X-13-03	0	III	1217	55772	X-23-01	0	VII	3385
55628	X-02-02	0	I	3305	55706	X-13-04	0	III	1161	55773	X-23-02	0	VII	2218
55630	X-02-03	0	I	1876	55707	X-13-05	0	IV	2763	55774	X-23-03	0	V	1811
55632	X-03-00	0	I	1712	55708	X-14-00	0	IV	1188	55775	X-23-04	0	V	3356
55634	X-03-01	0	III	3590	55709	X-14-01	0	IV	3342	55776	X-23-05	0	V	2603
55639	X-04-00	0	IV	3099	55710	X-14-02	0	IV	1094	55777	X-23-06	0	IV	912
55642	X-04-02	0	IV	2972	55712	X-14-03	0	IV	1404	55777	X-23-06	1	IV	1544
55643	X-04-01	0	III	1190	55713	X-14-04	0	V	871	55778	X-24-00	0	IV	1309
55644	X-04-03	0	III	2903	55714	X-14-05	0	V	1311	55779	X-24-01	0	IV	3599
55645	X-05-02	0	I	3578	55715	X-15-00	0	IV	1241	55779	X-24-02	0	IV	2013
55647	X-05-00	0	IV	2872	55716	X-15-01	0	IV	1262	55782	X-24-03	0	V	1761
55647	X-05-000	0	IV	3472	55717	X-15-02	0	III	1557	55782	X-24-04	0	V	1725
55647	X-05-000	1	IV	3520	55718	X-15-03	0	III	1334	55784	X-24-05	0	V	3144
55647	X-05-000	2	IV	3512	55720	X-15-04	0	IV	1486	55784	X-24-06	0	V	2591
55647	X-05-000	3	IV	3520	55721	X-15-05	0	IV	1500	55785	X-25-00	0	V	2366
55647	X-05-000	4	IV	3512	55722	X-16-00	0	IV	900	55786	X-25-01	0	V	1804
55647	X-05-000	5	IV	648	55723	X-16-01	0	III	1004	55787	X-25-02	0	V	1951
55649	X-05-03	0	IV	2409	55724	X-16-02	0	II	1923	55788	X-25-03	0	V	1619
55650	X-05-01	0	IV	1473	55725	X-16-03	0	II	1919	55789	X-25-04	0	V	2601
55651	X-05-04	0	IV	2954	55726	X-16-04	0	III	1935	55790	X-25-05	0	V	1473
55653	X-06-00	0	IV	2723	55727	X-16-05	0	II	730	55791	X-25-06	0	V	922
55654	X-06-01	0	IV	3388	55728	X-16-06	0	II	1953	55792	X-26-00	0	V	2336
55656	X-06-02	0	IV	2908	55729	X-17-00	0	II	2735	55794	X-26-01	0	V	1385
55657	X-06-03	0	V	1842	55730	X-17-01	0	II	3556	55795	X-26-02	0	VIII	1458
55661	X-07-00	0	V	1754	55731	X-17-02	0	II	3605	55796	X-26-03	0	VI	1325
55662	X-07-01	0	V	3365	55732	X-17-03	0	II	1647	55798	X-26-04	0	VI	2075
55663	X-07-02	0	V	3373	55733	X-17-04	0	II	1459	55799	X-27-00	0	VI	1396
55666	X-08-00	0	V	3338	55734	X-17-05	0	III	1736	55800	X-27-01	0	VI	2684
55669	X-08-01	0	V	3368	55735	X-17-06	0	III	3653	55801	X-27-02	0	VI	1016
55670	X-08-03	0	VI	2489	55736	X-18-00	0	III	2317	55802	X-27-03	0	VI	1179
55671	X-08-02	0	VI	2609	55737	X-18-01	0	IV	1387	55803	X-27-04	0	VI	1304
55673	X-09-03	0	VI	1011	55738	X-18-02	0	V	1291	55805	X-27-05	0	VI	1663
55674	X-09-00	0	VI	1386	55739	X-18-03	0	V	2178	55806	X-28-00	0	VI	1456
55675	X-09-05	0	IX	1148	55740	X-18-04	0	V	1478	55808	X-28-01	0	VIII	577

Table 8. Table 8: A continuation of Table 7 . Orb. is the orbit ID (starting at 0) of each observation segment, Exp. is the exposure time in seconds and X is the prefix 96420-01.

MJD	OBSID	Orb.	Class	Exp.	MJD	OBSID	Orb.	Class	Exp.	MJD	OBSID	Orb.	Class	Exp.
55810	X-28-02	0	VI	1251	55836	X-32-02	0	IX	1591	55857	X-35-01	0	IX	1912
55811	X-28-03	0	VI	2000	55837	X-32-03	0	IX	2155	55859	X-35-02	0	IX	200
55813	X-29-00	0	VIII	1309	55838	X-32-04	0	IX	2641	55859	X-35-02	1	IX	1296
55819	X-29-04	0	VIII	1686	55838	X-32-05	0	IX	2077	55860	X-35-03	0	IX	1372
55820	X-30-00	0	VI	1488	55840	X-32-06	0	IX	3392	55861	X-35-04	0	IX	836
55821	X-30-01	0	VI	1503	55840	X-32-06	1	IX	3512	55862	X-36-00	0	IX	1145
55822	X-30-02	0	VI	1417	55840	X-32-06	2	IX	3934	55863	X-36-01	0	IX	1322
55823	X-30-03	0	VI	1290	55840	X-32-06	3	IX	3880	55865	X-36-03	0	IX	1485
55824	X-30-04	0	VI	1489	55840	X-32-06	4	IX	1896	55866	X-36-04	0	IX	1795
55825	X-30-05	0	VI	2581	55841	X-33-00	0	IX	1188	55867	X-36-05	0	IX	1732
55826	X-30-06	0	VI	2747	55842	X-33-01	0	IX	855	55868	X-36-06	0	IX	1657
55827	X-31-00	0	VI	1559	55843	X-33-02	0	IX	1156	55871	X-37-00	0	IX	815
55828	X-31-01	0	VI	2954	55845	X-33-04	0	IX	1713	55871	X-37-02	0	IX	1460
55829	X-31-02	0	IX	3005	55846	X-33-05	0	IX	934	55872	X-37-03	0	IX	1683
55830	X-31-03	0	IX	1472	55847	X-33-06	0	IX	717	55873	X-37-04	0	IX	1402
55830	X-31-03	1	IX	288	55848	X-34-00	0	IX	1159	55874	X-37-05G	0	IX	1536
55831	X-31-04	0	IX	1586	55849	X-34-01	0	IX	973	55875	X-37-06	0	IX	1536
55832	X-31-05	0	VI	3812	55851	X-34-02	0	IX	2261	55876	X-38-00	0	IX	1497
55833	X-31-06	0	IX	3675	55852	X-34-03	0	IX	1092	55877	X-38-01	0	IX	1134
55834	X-32-00	0	IX	1217	55853	X-34-04	0	IX	741	55878	X-38-02	0	IX	1289
55835	X-32-01	0	IX	1445	55856	X-35-00	0	IX	797	55879	X-38-03	0	IX	1433

Equations

$$HSS=\frac{2(AD-BC)}{(A+B)(B+D)+(A+C)(C+D)}$$

Full text

††pubyear: 2017††pagerange: An Atlas of Exotic Variability in IGR J17091-3624: A Comparison with GRS 1915+105–B

1 Introduction

X-ray binaries are systems in which a black hole or neutron star accretes matter from a stellar companion, and they provide us with opportunities to test how accretion takes place in the most extreme physical regimes. Some X-ray Binaries are believed to be accreting at very close to the Eddington limit, the limit at which the radiation pressure on accreting material is equal to the force due to gravity. As such, these objects can also provide a laboratory with which to explore accretion in radiation pressure-dominated systems (White & Zhang, 1997).

Low Mass X-ray Binaries (hereafter LMXBs) are a subclass of X-ray binary in which the compact object accretes matter transferred to it due to a Roche-lobe overflow from the companion star (e.g. Paczynski, 1979). In general, accretion in LMXBs is a variable process, with variability seen on timescales from milliseconds to decades. On the shortest timescales the X-ray lightcurves of these objects can show band-limited noise and low-frequency quasi-periodic oscillations (QPOs) at frequencies from $\sim$mHz to $\sim$200 Hz (e.g. van der Klis, 1989). Black hole binaries also show so-called ‘high frequency QPOs’ (e.g. Remillard et al., 1999a; Remillard et al., 1999b; Belloni et al., 2002; Belloni & Motta, 2016), thought to be caused by motion of matter in the innermost region of the accretion disk (e.g. Stefanov, 2014).

Three sources – GRS 1915+105, IGR J17091-3624 and the neutron star ‘Rapid Burster’ (MXB 1730-335) – also show a variety of exotic variability on timescales of seconds to minutes in addition to the kinds of variability seen in other LMXBs. This exotic variability consists of quasi-periodic flares, dips and other high-amplitude behaviours (e.g Belloni et al., 2000; Altamirano et al., 2011a; Bagnoli & in’t Zand, 2015). The second-to-minute-scale lightcurve profiles of these sources change over timescales of days. In GRS 1915+105 and IGR J17091-3624, this behaviour can be described as a set of ‘variability classes’. These classes themselves vary widely in terms of flux, structure, periodicity and spectral properties.

GRS 1915+105 (Castro-Tirado et al., 1992), hereafter GRS 1915, is a black hole LMXB which accretes at between a few tens and 100% of its Eddington Limit (e.g. Vilhu, 1999; Done et al., 2004; Fender & Belloni, 2004; Reid et al., 2014). Most LMXBs go through periods of low-intensity ‘quiescence’ and high-intensity ‘outbursts’, the latter consisting of black-body dominated ‘soft’ and power-law dominated ‘hard’ spectral states. However, GRS 1915 has been in outburst since its discovery in 1992 (Castro-Tirado et al., 1992). GRS 1915 is also notable for the incredible variety and complexity of variability classes it exhibits (e.g. Yadav et al., 2000; Belloni et al., 2000) in addition to the less exotic variability seen in other black hole binary systems. GRS 1915 additionally shows high-frequency and low frequency QPOs similar to those seen in other black hole LMXBs (Morgan et al., 1997). In total, 15 distinct variability classes have been observed (Belloni et al., 2000; Klein-Wolt et al., 2002; Hannikainen et al., 2007; Pahari & Pal, 2009). This remarkable range of behaviour is believed to be caused by instability in the inner accretion disc (e.g. Janiuk et al., 2000; Nayakshin et al., 2000), which is in turn caused by the existence of a radiation pressure dominated regime in the inner disc (e.g. Done et al., 2004). Accounting for this complexity could be key to our understanding of radiation-dominated accretion regimes.

One of the best-studied variability classes of GRS 1915 is the highly regular flaring $\rho$, or ‘heartbeat’, class, so named for the similarity of its lightcurve to an electrocardiogram. It has been shown that hard X-ray photons tend to lag soft ones in this class (e.g. Janiuk & Czerny, 2005; Massaro et al., 2010). Numerical models derived from Shakura & Sunyaev, 1973 which reproduce this lag can also reproduce other flaring classes seen in GRS 1915 (e.g Nayakshin et al., 2000; Massaro et al., 2014). These numerical models predict that GRS 1915-like variability should be seen in systems accreting with a global Eddington fraction of $\gtrsim 0.26$ (Nayakshin et al., 2000). However, other LMXBs (e.g. GX 17+2, Kuulkers et al., 2002, and V404 Cyg, Huppenkothen et al., 2016) have been observed to exceed this Eddington fraction without displaying GRS 1915-like variability.

Neilsen et al. (2011) proposed a physical scenario, based on the mathematical model proposed by Nayakshin et al. (2000), to explain the presence of the hard lag in the flaring classes of GRS 1915. This is outlined schematically in Figure 1. First, an overdensity of matter forms via the thermal-viscous Lightman-Eardley Instability (Lightman & Eardley, 1974) and propagates inwards through the accretion disc. This destabilises the disc, collapsing its inner radius and vastly increasing photon emission. If the local Eddington Limit in the inner accretion disc is then approached, extreme outflows are triggered that deplete the inner accretion disc and allow the cycle to begin again. As the matter ejected from the disc collides with the non-thermal ‘corona’ above the central object, a flash of hard Bremsstrahlung radiation is produced. This causes a hardening of the spectrum and an apparent lag between soft and hard photons. Janiuk & Czerny (2005) instead propose that the lag is caused by the corona smoothly adjusting to the changing brightness of the disc after a light travel time.

The black hole candidate LMXB IGR J17091-3624 (hereafter J17091) was discovered in outburst by INTEGRAL in 2003 (Kuulkers et al., 2003). In 2011, it again went into outburst (Krimm & Kennea, 2011). GRS 1915-like variability was discovered in its lightcurve, as well as high-frequency QPOs which behave much like the QPOs seen in GRS 1915 (Altamirano et al., 2011a, c; Altamirano & Belloni, 2012). As IGR J17091 is around a factor of 20 fainter at 2–25 keV than GRS 1915, this object has either a lower black hole mass $M$, a lower accretion rate $\dot{m}$ or lies at a larger distance $D$ than GRS 1915. Assuming by analogy with GRS 1915 that IGR J17091 is accreting at its Eddington rate, the black hole must have a mass of $M\lesssim 3$M⊙ or lie at a distance of $D\gtrsim 20$ kpc (Altamirano et al., 2011a). However, independent estimates based on empirical relationships between black hole mass and high-frequency QPOs have suggested values of $M$ between $\sim 8.7$ and $15.6$ M⊙ (Rebusco et al., 2012; Iyer et al., 2015a, b), while multi-wavelength observations of the hard-to-soft state transition have suggested values of $D$ between $\sim$11 and $\sim$17 kpc (Rodriguez et al., 2011a). This implies that IGR J17091 may have an accretion rate $\dot{m}$ that is significantly below its Eddington rate.

The suggestion that IGR J17091 accretes at significantly below the Eddington Limit raises several questions. Despite evidence of disc winds in IGR J17091 faster than 0.03c (King et al., 2012), the wind generation procedure described in Neilsen et al. (2011) cannot take place without near-Eddington-limited accretion. Additionally, it makes it increasingly unclear what differentiates IGR J17091 and GRS 1915 from the vast majority of LMXBs which do not display such complex behaviour, and what physical system properties are required for GRS 1915-like variability to be observed. This latter point was further complicated by the observation of GRS 1915-like variability in 2 out of 155 RXTE observations of the Rapid Burster (Bagnoli & in’t Zand, 2015). As thermonuclear (Type I) X-ray bursts are also seen in the Rapid Burster (Hoffman et al., 1978), it is known to contain a neutron star. As such, the presence of variability classes in this object rules out any black hole-specific effects as the cause of the complex variability. In addition to this, the persistent luminosity of the Rapid Burster during periods of GRS 1915-like variability is known to be no greater than $\sim 10\%$ of its Eddington limit (Bagnoli et al., 2015).

Accounting for GRS 1915-like variability is required for a complete understanding of the physics of accretion in LMXBs. As such, Belloni et al. (2000) performed a complete model-independent analysis of variability classes in GRS 1915. This work highlighted the breadth and diversity of variability in GRS 1915, and allowed these authors to search for features common to all variability classes. For example, Belloni et al. (2000) also found that every variability class can be expressed as a pattern of transitions between three quasi-stable phenomenological states.

Previous works have noted that some of the variability classes seen in IGR J17091-3624 appear very similar to those seen in GRS 1915 (e.g. Altamirano et al., 2011a; Zhang et al., 2014). However, although $\rho$-like classes in the two objects both show lags between hard and soft X-rays photons, these lags appear to possess different signs (Altamirano et al., 2011a). Additionally, at least two variability classes have been reported in IGR J17091 which have not yet been reported in GRS 1915 (Pahari et al., 2012). Previous works have described some of the behaviour seen in IGR J17091 in the context of the variability classes described by Belloni et al. 2000 for GRS 1915 (e.g. Altamirano et al., 2011a; Pahari et al., 2014). To further explore the comparison between GRS 1915 and IGR J17091, here we perform the first comprehensive model-independent analysis of variability classes in IGR J17091 using the complete set of RXTE (Bradt et al., 1993) data taken of the 2011-2013 outburst of the object. We also use data from all other X-ray missions that observed the source during this time to analyse the long-term evolution of the outburst.

2 Data and Data Analysis

In this paper, we report data from RXTE, INTEGRAL, Swift, Chandra, XMM-Newton and Suzaku covering the 2011-2013 outburst of IGR J17091. Unless stated otherwise, all errors are quoted at the 1$\sigma$ level.

In Figure 2 we present long-term lightcurves from RXTE, INTEGRAL and Swift with one datapoint per observation, as well as indicating when during the outburst Chandra, XMM-Newton and Suzaku observations were made.

2.1 RXTE

For our variability study, we focus on the data from the Proportional Counter Array (PCA, Jahoda et al., 1996) aboard the Rossi X-Ray Riming Experiment (RXTE, Bradt et al., 1993). We analysed all 2011 PCA observations of IGR J17091, corresponding to OBSIDs 96065-03, 96103-01 and 96420-01. The observations taken for proposals 96065-03 and 96103-01 were contaminated by the nearby X-ray source GX 349+2 (Altamirano et al., 2011a; Rodriguez et al., 2011b). As such we only use observations performed for proposal 96420-01, corresponding to a total of 243 orbits from 215 separate observations. These were offset by 25’ such that GX 349+2 was not in the $1^{\circ}$ PCA field of view. RXTE was decommissioned during a period of Sun constraint centred on MJD 55907, and hence the last observation of IGR J17091 was taken on MJD 55879.

We extracted data from the native FITS format using our own software111https://github.com/jmcourt/pantheon. To perform medium- to high-frequency ($\gtrsim 1$ Hz) timing analysis, we merged files formatted in PCA ’s ‘Good Xenon’ data mode and extracted their data at the maximum time resolution ($\sim 9.5\times 10^{-7}$ s) without accounting for the background. We divided these data into 128 s segments as this allowed us to reach frequencies below $\sim 0.015$ Hz, partly sampling the high amplitude quasi-periodic flaring behaviour seen in many classes. Using the Fast Fourier Transform (FFT), we produced the power spectrum of each segment separately. We then averaged these spectra to create a one co-added Power Density Spectrum (PDS) for each observation.

For low-frequency ($\leq 1$ Hz) timing and correlated spectral/timing analysis, we rebinned the data to 0.5 s and normalised count rates by the number of proportional counters (PCUs) active in each observation. Our choice of 1 Hz allows us to analyse high amplitude ‘flaring’ behaviour (seen at frequencies $\lesssim 0.5$ Hz) separately from the lower-amplitude behaviour seen at $\gtrsim 5$ Hz.

We split the data into three energy bands: A (PCA channels 0–14, $\sim 2$–$6$ keV), B (PCA channels 15–35, $\sim 6$–$16$ keV) and C (PCA channels 36–255, $\sim 16$–$60$ keV). We chose these energy bands to be consistent with the energy bands used by the model-independent classification of variability classes of GRS 1915 in Belloni et al. (2000). For each of the energy-filtered lightcurves produced we estimated background using pcabackest from the FTOOLS package (Blackburn, 1995) with the PCA faint source background model222http://heasarc.gsfc.nasa.gov/FTP/xte/calib_data/pca_bkgd/Faint/pca_bkgd_cmfaintl7_eMv20051128.mdl. In all observations, we found that counts in the C band were consistent with background. We then created Lightcurves $L_{A}$ and $L_{B}$ from background-subtracted photons counted in the A and B bands respectively. We used these lightcurves to define the full-band lightcurve ($L_{T}=L_{A}+L_{B}$) and the soft colour ($C_{1}=L_{B}/L_{A}$) of each observation. To complement the Fourier spectra, we also constructed Generalised Lomb-Scargle Periodograms of $L_{T}$ from each dataset, a modified version of the standard Lomb-Scargle periodogram (Lomb, 1976; Scargle, 1982) that takes into account errors in the dataset (Irwin et al., 1989). Using the Lomb-Scargle periodogram instead of the Fourier periodogram here allows us to sample the low-frequency behaviour of lightcurves with data gaps. This is important, for example, in lightcurves which show two populations of flares, as it allows each population to be studied independently by cropping the other from the lightcurve.

We also used data from Altamirano et al. 2011a to sample the long-term colour evolution of IGR J17091. We use 2 hardness ratios defined by Altamirano et al.: $H_{A1}$ and $H_{A2}$, corresponding to the ratios of the 2–3.5 keV band against the 3.5–6 keV band and the 6–9.7 keV band against the 9.7–16 keV band respectively.

When possible, if low-frequency peaks were present in the Lomb-Scargle spectrum of an observation, we used the position of the highest peak to define a value for a period. This period was then used to rebin data by phase (or ‘fold’ the data) to search for reccurent hysteretic patterns in the hardness-Intensity diagram (hereafter HID1, a plot of $L_{T}$ against $C_{1}$). We found that quasi-periodic oscillations in our observations tended to show significant frequency shifts on timescales shorter than the length of the observations. As such, we employed the variable-period folding algorithm outlined in Appendix A where appropriate. For cases in which this algorithm was not appropriate, we considered small sections of each lightcurve, with a length equivalent to small number of periods, before performing folding.

Additionally, in observations which showed a pattern of high-amplitude X-ray flaring in $L_{T}$, we used our own algorithm to find individual flares (this algorithm is described in Appendix A) and collect statistics on the amplitude, duration and profile of these events.

A list of all observations used in this study can be found in Appendix B.

2.2 Swift

In this paper, we consider data from the Burst Alert Telescope (BAT, Barthelmy 2000), and the X-ray Telescope (XRT, Burrows et al. 2003) aboard Swift Gamma-ray Burst Mission (Gehrels, 2004). IGR J17091-3624 was observed with XRT for a total of 172 pointed XRT observations between MJDs 55575 and 56600, corresponding to Target IDs 31921, 34543, 30967, 30973, 31920, 35096, 67137, 81917, 522245, 677582 and 677981. These observations were interrupted during sun constraints centred on MJDs 55907 and 56272. We created a long-term 0.3–10 keV Swift/XRT light curve, with one bin per pointed observation, using the online light-curve generator provided by the UK Swift Science Data Centre (UKSSDC; Evans et al., 2007). We have also created a long-term 15–50 keV lightcurve using the publicly available Swift/BAT daily-averaged light curve333http://swift.gsfc.nasa.gov/results/transients/weak/IGRJ17091-3624/. These are shown in Figure 2 Panels (b) and (c) respectively.

2.3 INTEGRAL

The INTErnational Gamma Ray Astrophisical Laboratory (INTEGRAL, Winkler et al., 2003) is a medium-sized ESA mission launched in 2002. Unique hard X-ray (15–1000 keV with the ISGRI detector plane) sensitivity and wide field of view make INTEGRAL ideally suited to surveying the hard X-ray sky (see Bird et al. (2016) and Krivonos et al. (2015) for recent surveys).

INTEGRAL observations are divided into short ($\sim$2 ks) pointings called Science Windows (ScWs). We analyse all available observations of IGR J17091 with INTEGRAL/IBIS (Ubertini et al., 2003) between MJD 55575–55625 where the source is less than 12 degrees from the centre of the field of view and where there is more the 1 ks of good ISGRI time per ScW. This corresponds to the spectrally hardest period of the 2011-2013 outburst. The filtering of observations results in a total of 188 Science Windows which were processed using the Offline Science Analysis (OSA) software version 10.2 following standard data reduction procedures444http://www.isdc.unige.ch/integral/analysis in four energy bands (20–40, 40–100, 100–150, 150–300 keV). These bands were selected as they are standard energy bands used in the surveys of Bird et al. (2016) and Bazzano et al. (2006) and allow comparison to these previous works. Images were created at the ScW level as well as a single mosaic of all Science Windows in each energy band.

2.4 XMM-Newton

In this paper we only consider data from the European Photon Imaging Camera (EPIC Bignami et al., 1990) aboard XMM-Newton (Jansen et al., 2001). EPIC consists of one telescope with a pn-type CCD (EPIC-pn, Strüder et al., 2001) and two telescopes with MOS CCDs (EPIC-MOS1 and -MOS2, Turner et al., 2001).

XMM/Newton observed IGR J17091 thrice during the period from 2011–2013 (represented by the blue arrows in Figure 2). One of these (OBSID 0721200101) was made on 12 September 2013; we do not consider this observation further as IGR J17091 had returned to quiescence by this time (Altamirano et al., 2013). The remaining two observations, corresponding to OBSIDs 0677980201 and 0700381301 respectively, were taken on March 27 2011 (MJD 55647) and September 29 2012 (MJD 56199).

During observation 0677980201, EPIC-pn was operating in burst mode and EPIC-MOS was operating in timing mode. Given the low efficiency of burst mode, we only consider data from EPIC-MOS for this observation. During observation 0700381301, EPIC-pn was operating in timing mode, and thus we use data from EPIC-pn for this observation.

We used the XMM-Newton Science Analysis Software version 15.0.0 (SAS, see Ibarra et al., 2009) to extract calibrated event lists from EPIC in both observations. We used these to construct lightcurves to study the X-ray variability, following standard analysis threads555http://www.cosmos.esa.int/web/xmm-newton/sas-threads.

2.5 Chandra

In this paper we consider data from the Advanced CCD Imaging Camera (ACIS, Nousek et al., 1987) and the High Resolution Camera (HRC, Murray et al., 1987) aboard Chandra (Weisskopf, 1999). Chandra made 7 observations of IGR J17091 during the period 2011–2013. Four of these observations were taken after IGR J17091 returned to quiescence, and we do not consider these further in this paper. The Chandra observations log is reported in Table 1.

We analysed these data using CIAO version 4.8 (Fruscione et al., 2006), following the standard analysis threads. In order to apply the most recent calibration files (CALDB 4.7.0, Graessle et al., 2006), we reprocessed the data from the three observations using the chandra_repro script666See e.g. http://cxc.harvard.edu/ciao/ahelp/chandra_repro.html, and used this to produce data products following standard procedures.

The first Chandra observation(OBSID 12505) of this source was made shortly after it went into outburst in February 2011. It was a 1 ks observation performed to refine the position of the X-Ray source, using the High-Resolution Camera in Imaging mode (HRC-I). We created the 0.06–10 kev light curve accounting for the Dead Time Factor (DTF) to correct the exposure time and count rate, due to the deviation of the detector from the standard detection efficiency, using the dmextract tool in the CIAO software.

Two additional observations (OBSIDs 12405 and 12406) were performed within 214 days of this first observation, using the High Energy Transmition Grating Spectrometer (HETGS) on board Chandra. The incident X-Ray flux was dispersed onto ACIS using a narrow strip array configuration (ACIS-S). Continuous Clocking and Time Exposure modes were use in each observation respectively (see King et al., 2012 for further details). We exclude any events below 0.4 keV, since the grating efficiency is essentially zero below this energy. In the case of the OBSID 12405 observations we also excluded the Flight Grade 66 events in the event file, as they were not appropriately graded. We extracted the 0.5-10 kev HEGTS light curves, excluding the zeroth-order flux, adopting standard procedures.

2.6 Suzaku

In this paper, we only consider data from the X-ray Imaging Spectrometer (XIS, Koyama et al., 2007) aboard Suzaku (Mitsuda et al., 2007). Suzaku observed IGR J17091 twice during the period 2011–2013; a 42.1 ks observation on October 2–3, 2012 (MJD 56202–56203, ObsID: 407037010) and an 81.9 ks observation on February 19–21, 2013 (MJD 56342–56344, ObsID: 407037020). XIS consists of four X-ray CCDs (XIS 0, 1, 2 and 3), and all them except for XIS 2 were operating in the 1/4 window mode which has a minimum time resolution of 2 seconds.

We analysed the Suzaku data using HEASOFT 6.19 in the following standard procedures after reprocessing the data with aepipeline and the latest calibration database (version 20160607). We extracted XIS light curves in the 0.7–10 keV range, and subtracted background individually for XIS 0, 1 and 3 and then summed these to obtain the total background. We created Power density spectra (PDS) using powspec in the XRONOS package.

3 Results

3.1 Outburst Evolution

The 2011-2013 outburst of IGR J17091-3624 was first detected with Swift/XRT on MJD 55595 (3 Feb 2011) (Krimm & Kennea, 2011), and was observed by RXTE , Swift/XRT, Swift/BAT and INTEGRAL/ISGRI (see Figure 2, Panels a, b, c, d respectively). There were also pointed observations by XMM-Newton, Chandra and Suzaku during this time (denoted by coloured arrows in Figure 2).

RXTE data were taken within the first week of the outburst, but they were heavily contiminated by the nearby source GX 349+2 (Rodriguez et al., 2011b). Thus these data are not considered here.

The onset of the outburst can be seen in the Swift/BAT lightcurve (Figure 2 Panel c). In a 22 day period between MJDs 55584 and 55608, the 15–50 keV intensity from IGR J17091 rose from $\sim 9$ mCrab to a peak of $\sim 110$ mCrab. This onset rise in intensity can also be seen in 0.3–10 keV Swift/XRT data and 20–40 keV INTEGRAL/ISGRI data.

After peak intensity, the 15–50 keV flux (Swift/BAT) began to steadily decrease, until returning to a level of $\sim$20 mCrab by MJD 55633. A similar decrease in flux can be seen in the data obtained by INTEGRAL at this time (Figure 2 Panel (d). However, there was no corresponding fall in the flux at lower energies; both the long-term 2–16 keV RXTE data and Swift/XRT data (Panels a and b respectively) show relatively constant fluxes of 45 mCrab between MJDs 55608 and 55633.

The significant decrease in high-energy flux during this time corresponds to IGR J17091 transitioning from a hard state to a soft intermediate state (Pahari et al., 2014). This transition coincides with a radio flare reported by Rodriguez et al. (2011a) which was observed by the Australian Telescope Compact Array (ATCA).

Altamirano et al., 2011b first reported a 10 mHz QPO in RXTE data on MJD 55634 , evolving into ‘Heartbeat-like’ flaring by MJD 55639 (Altamirano et al., 2011c). Between MJDs 55634 and 55879, the global RXTE lightcurve shows large fluctuations in intensity on timescales of days to weeks, ranging from a minimum of $\sim 1$ mCrab on MJD 55731 to a maximum of $\sim 66$ mCrab on MJD 55756. The Swift/XRT lightcurve shows fluctuations that mirror those seen by RXTE during this period, but the amplitude of the fluctuations is significantly reduced.

Swift/XRT was unable to observe again until MJD 55952. Between this date and MJD 55989, Swift/XRT observed a gradual decrease in intensity corresponding to a return to the low/hard state (Drave et al., 2012).

Between MJD 55989 and the end of the outburst on MJD 56445, we see secondary peaks in the Swift/XRT, Swift/BAT and INTEGRAL/ISGRI lightcurves that evolve over timescales of $\lesssim 100$ days. Similar humps have been seen before in lightcurves from other objects, for example the black hole candidate XTE J1650-500 (Tomsick et al., 2003) and the neutron stars SAX J1808.4-3658 (Wijnands et al., 2001) and SAX J1750.8-2900 (Allen et al., 2015). These humps are referred to as ‘re-flares’ (also as ‘rebrightenings’, ‘echo-outbursts’, ‘mini-outbursts’ or a ‘flaring tail’, e.g. Patruno et al., 2016). We see a total of 3 apparent re-flares in the Swift/BAT data, centred approximately at MJDs 56100, 56220 and 56375.

The observation with XMM-Newton/EPIC-pn on MJD 56547 (12 September 2013) recorded a rate of 0.019 cts s*-1*. An observation with EPIC-pn in 2007, while IGR J17091 was in quiescence (Wijnands et al., 2012), detected a similar count rate of 0.020 cts s*-1*. Therefore we define MJD 56547 as the upper limit on the endpoint of the 2011-2013 outburst. As such the outburst, as defined here, lasted for $\lesssim$952 days.

After the end of the 2011-2013 outburst, IGR J17091 remained in quiescence until the start of a new outburst around MJD 57444 (26 February 2016, Miller et al., 2016).

3.2 RXTE

Using the data products described in Section 2, we assigned a model-independent variability class to each of the 243 RXTE/PCA orbits. To avoid bias, this was done without reference to the classes defined by Belloni et al. (2000) to describe the behaviour of GRS 1915.

Classes were initially assigned based on analysis of lightcurve profiles, count rate, mean fractional RMS (Vaughan et al., 2003), Fourier and Lomb-scargle power spectra and hardness-intensity diagrams. For observations with significant quasi-periodic variability at a frequency lower than $\sim 1$ Hz, we also attempted to fold lightcurves to analyse count rate and colour as a function of phase. When flares were present in the lightcurve, we used our algorithm (described in Appendix A) to sample the distribution of parameters such as peak flare count rate, flare rise time and flare fall time. All parameters were normalised per active PCU, and fractional RMS values were taken from 2–60 keV lightcurves binned to 0.5 s. We identify nine distinct classes, labelled I to IX; we describe these in the following sections.

Although the criteria for assigning each class to an observation was different, a number of criteria were given the most weight. In particular, the detection, $q$-value and peak frequency of a QPO in the range 2 Hz–10 Hz were used as criteria for all classes, as well as the presence or absence of high-amplitude quasi-periodic flaring with a frequency between 0.01–1 Hz. The folded profile of these flares, as well as the presence of associated harmonics, were also used as classification diagnostics in observations. Additionally, the presence or absence of low count-rate ’dips’ in a lightcurve was used as a criterion for Classes VI, VIII and IX. Detailed criteria for each individual class are given below in Sections 3.2.1 to 3.2.9.

For hardness-intensity diagrams, we describe looping behaviour with the terms ‘clockwise’ and ‘anticlockwise’; in all cases, these terms refer to the direction of a loop plotted in a hardness-intensity diagram with colour on the $x$-axis and intensity on the $y$-axis.

In Appendix B, we present a list of all orbits used in the study along with the variability classes we assigned to them.

In Figure 3, we show global 2–16 keV lightcurves of IGR J17091 during the 2011-2013 outburst. In each panel, all observations of a given class are highlighted in red. A characteristic lightcurve is also presented for each class. In Figure 4 panel (a), we show a plot of average hardness $H_{A2}$ against $H_{A1}$ for each observation, showing the long-term hysteresis of the object in colour-colour space. Again, observations belonging to each variability class are highlighted. In Figure 4 panels (b) and (c), we show global hardness-intensity diagrams for $H_{A1}$ and $H_{A2}$ respectively.

In Figure 4 Panel (a), we see that IGR J17091-3624 traces a two branched pattern in colour-colour space corresponding to a branch which is soft ($\sim 0.9$) in $H_{A1}$ and variable in $H_{A2}$ and a branch which is soft ($\sim 0.5$) in $H_{A2}$ and variable in $H_{A1}$. The ‘soft’ HID shown in Figure 4 Panel (b) is dominated by a branch with a wide spread in $H_{A1}$ and intensities between $\sim 40\mbox{--}60$ mCrab. A second branch exists at lower intensities, and shows an anticorrelation between intensity and $H_{A1}$. Finally, the ‘hard’ HID shown in Figure 4 Panel (c) shows an obvious anticorrelation between $H_{A2}$ and intensity, but there is also a secondary branch between $H_{A2}\approx 0.7\mbox{--}0.9$ at a constant intensity of $\sim 40$ mCrab.

For characteristic count rates and colours in each class, we quote the upper and lower quartile values (Kenney, 1939) instead of the mean. This is due to the presence of high-amplitude but short-lived flares in many of the classes we describe. Using the upper and lower quartiles as our measure of average and distribution means that our values will be less susceptible to outlier values of count rate and colour present in these flares. All count rates have been background corrected (see Section 2.1).

We have obtained mean values for these count rate quartiles, as well as values for colour $C_{1}$ and fractional RMS, by calculating these values individually for each orbit. Histograms were then constructed from these datasets for each class, such that the mean and standard deviation of these values could be measured for each class. These values are presented in Table 2.

We describe QPOs in terms of their $q$-value; a measure of coherence defined by the ratio of peak frequency and full-width half-maximum of each QPO. We collected these values by fitting our power spectra with Lorentzians.

For each class, we present three standard data products; a 500 s lightcurve, a variable-length lightcurve where the length has been selected to best display the variability associated with the class and a Fourier PDS. Unless otherwise stated in the figure caption, the 500 s lightcurve and the Fourier PDS are presented at the same scale for all classes. In Table 3 we present a tally of the number of times we assigned each Variability Class to an RXTE orbit.

3.2.1 Class I – Figure 5

In the 2 s binned lightcurve of a Class I observation, there is no structured second-to-minute scale variability. The Fourier PDS of all observations in this class shows broad band noise between $\sim 1$–$10$ Hz, as well as a weak QPO (with a $q$-value of $\sim 5$) which peaks at around 5 Hz.

3.2.2 Class II – Figure 6

Class II observations are a factor of $\sim 2$ fainter in the $L_{T}$ band than Class I observations. They also occupy a different branch in a plot of hardness $H_{A2}$ against flux (see Figure 4c). The PDS shows no significant broad band noise above $\sim 1$Hz unlike that which is seen in Class I. The $\sim$5Hz QPO seen in Class I is absent in Class II.

3.2.3 Class III – Figure 7

Unlike Classes I & II, Class III lightcurves show structured flaring, with a peak-to-peak recurrence time of $42$–$80$ s. Most flares consist of a steady $\sim 60$ s rise in count rate and then an additional and sudden rise to a peak count rate at $\gtrsim 200$ cts s*-1* PCU*-1*which lasts for $\lesssim$0.5 s before returning to continuum level (we have magnified the y-scaling in the lightcurve of Figure 7 to emphasise this behaviour). This sudden rise is not present in every flare; in some observations it is absent from every flare feature. No 5Hz QPO is present in the PDS and there is no significant variability in the range between $\sim 1\mbox{--}10$Hz.

As this class has a well-defined periodicity, we folded data in each observation to improve statistics using the best-fit period obtained from generalised Lomb-Scargle Periodogram Analysis; we show a representative Lomb-Scargle periodogram in Figure 8. We find an anticlockwise hysteretic loop in the folded HID1 of all 15 Class III orbits. In Figure 9 we show an example of one of these loops.

3.2.4 Class IV – Figure 10

The lightcurves in this class show regular variability with a peak-to-peak recurrence time of $25$–$39$ s. We performed peak analysis (see Appendix A) on observations belonging to this class, finding that each peak has a rise time with lower and upper quartile values of $19.5$ and $33.5$ s, a fall time with lower and upper quartile values of $4.6$ and $13.5$ s and a peak count rate of $159$–$241$ cts s*-1* PCU*-1* . There are no prominent significant QPOs in the Fourier PDS above $\sim 1$Hz.

We folded individual Class IV lightcurves and found anticlockwise hysteretic loops in the HID1 of 14 out of 30 Class IV observations. In the top panel of Figure 11 we show an example of one of these loops. However, we also find clockwise hysteretic loops in 6 Class IV observations, and in 10 orbits the data did not allow us to ascertain the presence of a loop. We provide an example of both of these in the lower panels of Figure 11. We note that the structure of clockwise loops are more complex than anticlockwise loops in Class IV, consisting of several lobes777In HIDs with multiple lobes, the loop direction we assign to the observation corresponds to the direction of the largest lobe. rather than a single loop (Figure 11, bottom-left).

Compared with Class III, the oscillations in Class IV occur with a significantly lower period, with a mean peak-to-peak recurrence time of $\sim 30$ s compared to $\sim 60$ s in Class III.

In Figure 4 we show that Classes III and IV can also be distinguished by average hardness, as Class III tends to have a greater value of $H_{A2}$ than Class IV.

3.2.5 Class V – Figure 12

The lightcurves in this class, like in Classes III and IV, show flaring behaviour, with flares separated by a few tens of seconds. At higher frequencies, the PDS shows a prominent QPO centred at $\sim 4$Hz with as $q$-value of $\sim 3$. There is also significant broad band noise between $\sim 0.1$–$1$Hz

In Figure 13 we show that the flaring in this class is more complex than that seen in Classes III and IV. Class V lightcurves consist of short strongly peaked symmetrical flares (hereafter Type $V_{1}$) and a longer more complex type of flare (hereafter Type $V_{2}$). The Type $V_{2}$ flare consists of a fast rise to a local maximum in count rate, followed by a $\sim 10$ s period in which this count rate gradually reduces by $\sim 50\%$ and then a much faster peak with a maximum count rate between 1 and 2 times that of the initial peak. In both types of flare, we find that the increase in count rate corresponds with an increase in soft colour. The two-population nature of flares in Class V can also clearly be seen in Figure 14, where we show a two-dimensional histogram of flare peak count rate against flare duration.

We folded all individual Class V lightcurves, in each case cropping out regions of $V_{2}$ flaring. We find clockwise hysteretic loops in the HID1 of 30 out of 33 Class V observations, suggesting a lag in the aforementioned relation between count rate and soft colour. In the upper panel Figure 15 we present an example of one of these loops. In one observation however, we found an anticlockwise loop in the HID1 (shown in Figure 15 lower-left panel). We were unable to ascertain the presence of loops in the remaining 2 orbits; for the sake of completeness, we show one of these in the lower-right panel of Figure 15.

3.2.6 Class VI – Figure 16

The lightcurves of observations of this class show large dips in count rate; this can be seen in Figure 16 at, for example, $t\approx 125$–$150$ s . These dips vary widely in duration, from $\sim 5$ to $\sim 50$ seconds, and the count rate in both $L_{A}$ and $L_{B}$ fall to a level consistent with background. The dips’ rise and fall times are fast, both lasting no longer than a second. They do not appear to occur with any regular periodicity.

Aside from the dips, Class VI observations show other structures in their lightcurves. Large fluctuations in count rate, by factors of $\lesssim 3$, occur on timescales of $\sim 1\mbox{--}5$ s; no periodicity in these oscillations could be found. This behaviour is reflected in the PDS, which shows high-amplitude broad band noise below $\sim 0.5$Hz with RMS-normalized power (Belloni & Hasinger, 1990) of up to $\sim 1.1$Hz*-1*. As can be seen in Figure 16, this feature takes the form of a broad shoulder of noise which shows a either weak peak or no clear peak at all. The $\sim 5$Hz QPO seen in the PDS of other classes is not present in Class VI observations.

We attempted to fold all individual Class VI lightcurves, ignoring the sections of data corresponding to the large count rate dips described above. In general, folding lightcurves belonging to this class is difficult; many orbits showed low-amplitude oscillations which were difficult to fold using our flare-finding algorithm (see Appendix A), while many others only showed oscillatory behaviour for a small number of periods between each pair of dips. As such, we only succesfully folded 23 of the 40 Class VI orbits. Of these, 19 showed clockwise loops in the HID1 (top panel, Figure 17), 3 showed anticlockwise loops (bottom-left panel, Figure 17). In the remaining 1 observation, the data did not allow us to ascertain the presence of loops (bottom-right panel, Figure 17).

Like in Class VI, we note that the clockwise loops in Class VI appear more complex than clockwise loops. Again, the clockwise loop shown in Figure 17 appears to have a 2-lobe structure; this is repeated in all clockwise loops found in this class.

3.2.7 Class VII – Figure 18

Class VII shows high-amplitude flaring behaviour with a peak-to-peak recurrence time of $6$–$12$ s. In Figure 19 we show a dynamical Lomb-Scargle spectrogram of a Class VII observation, showing that the fast flaring behaviour has a frequency which moves substantially over time. This in turn accounts for the large spread in the value of the flare peak-to-peak recurrence time.

In Figure 19 we show that the peak frequency of the QPO also varies in a structured way. We also suggest that the variabilitity of the frequency is itself a QPO with a period of $\sim 150$.

At higher frequencies, the PDS shows a weak QPOs centred at $\sim 8$Hz, with a $q$-values of $\sim 2$.

We used our flare-finding algorithm (see Appendix A) to perform variable-frequency folding of Class VII orbits. We find clockwise loops in 9 out of 11 Class VII orbits. In the remaining two observations, the oscillations were extremely fast. As a result, the errors in the HID1 of these too observations were too large to succesfully select peaks, and we are unable to confirm or reject the presence of loops.

3.2.8 Class VIII – Figure 20

The lightcurve of this variability class shows the dipping behaviour seen in Class VI, as can be seen in Figure 20 at $t\approx 125$–$150$ s. The dips are less frequent than in Class VI. The behaviour outside of the dips is dominated by highly structured high-amplitude oscillations consisting of flares with a peak to peak separation of $3.4\pm 1.0$ s. The PDS shows this behaviour as a very significant ($q$-value > 20) QPO; two harmonics of this QPO are also visible. The PDS also shows a strong ($q$-value = 4.7) QPO at $\sim 9$Hz.

We attempted to fold Class VIII lightcurves, ignoring the portions of data corresponding to dips, using our flare-finding algorithm. The high frequency of the dominant oscillation in Class VIII resulted in large errors in the peak times of individual flares, which translated to large errors in all HID1s; however, we were able to ascertain the presence in loops in 8 out of 16 orbits. All 8 of these loops are clockwise.

3.2.9 Class IX – Figure 21

The 1 s lightcurve of a Class IX observation is superficially similar to the lightcurve of a Class I observation, with little obvious structured variability at timescales larger than 2 s; however, large count rate dips like those seen in Classes VI and VIII (e.g. the feature at $t\approx 410$ s in the lightcurve of Figure 21) are very occasionally observed. These dips may in turn be coupled to short second-scale flares in which count rate briefly increases by a factor of 2–3.

Outside of these dips and flares, the lightcurve of a Class IX observation is indistinguishable from the lightcurve of a Class I or Class II observation. However, in Figure 4, we show that Class IX occupies a very different part of the global $H_{A2}$/$H_{A1}$ colour-colour diagram. Class IX observations show a significantly larger $H_{A2}$ than Class I and II observations, but a significantly lower $H_{A1}$.

The PDS reveals significant broad band noise peaked at $\sim$0.3 Hz, and the $\sim 5$Hz QPO seen in other classes is absent. Altamirano & Belloni (2012) discovered high frequency ($\sim 66$Hz) QPOs in observations corresponding to this variability class.

3.3 Swift

Observations with Swift took place throughout the 2011-2013 outburst of IGR J17091-3624. Between MJDs 55622 and 55880, 17 Swift/XRT were at least partly simultaneous with an RXTE observation, corresponding to at least one observation of all 9 classes. In each case, the Swift and RXTE lightcurves were similar. The remainder of the Swift/XRT observations during this time were also consistent with belonging to one of our nine classes. Given that the RXTE data have higher count rate and time resolution, we do not further discuss the Swift observations taken before MJD 55880. A more detailed comparison of RXTE and Swift data is beyond the scope of this paper.

Between MJD 55952 and 56445, Swift observations showed IGR J17091-3624 decreasing in flux. For all observations longer than 500 s, we rebinned the lightcurves to 10 s and calculated the RMS. We find the lower and upper quartiles of the fractional RMS in these measurements to be 18.3% and 21.7% respectively. INTEGRAL observations taken as part of a scan programme of the Galactic Plane (Fiocchi et al., 2012) and reported by Drave et al. (2012) suggest that IGR J17091-3624 returned to the hard state between MJDs 55952 and 55989. Therefore these observations sample IGR J17091-3624 the hard state.

3.4 INTEGRAL

The results of the INTEGRAL/IBIS analysis are presented in Table 4. We see clear detections of IGR J17091-3624 in all energy bands during the hardest period (MJD 55575–55625) of the 2011–2013 outburst. Conversion from detected counts to flux was achieved using an INTEGRAL/IBIS observation of the Crab taken between MJD 57305.334 and 57305.894. Conversion from Crab units to standard flux units was obtained by conversion factors listed in Bird et al. (2016) and Bazzano et al. (2006).

Comparing these results with those of Bazzano et al. (2006), we see that IGR J17091 is detected for the first time above 150 keV with a detection significance of 7.6 $\sigma$, corresponding to a flux of $2.24\pm 0.29\times 10^{-10}$ ergs s*-1* cm*-2* (Figure 22).

3.5 Chandra

In Figure 23, we present lightcurves from the three Chandra observations considered in this paper (see also Table 1 for details of these observations).

Observation 12505 was performed within 24 hours of RXTE observation 96420-01-02-01, which showed Class I variability. No structured variability is seen in the lightcurve of OBSID 12505 (Figure 23, upper panel), which is consistent with Class I. Note that we consider the energy range 0.06-10 keV for this observation but 0.5-10 keV for observations 12405 and 12406.

Observation 12405 was performed within 24 hours of RXTE observation 96420-01-23-03, which showed Class V variability. The two observations were not simultaneous; OBSID 12405 began $\sim 8.4$ ks after OBSID 96420-01-2303 finished. The lightcurve of Chandra OBSID 12405 (shown in Figure 23, middle panel) shows a mean count rate of 41 cts s*-1*. The lightcurve shows fast flaring behaviour (with a recurrence time on the order of 10s of seconds) in which the frequency changes widely on timescales of $\sim 1000$ s. This observation strongly resembles a Class VII lightcurve, but with its characteristic timescales increased by a factor of $\sim 4$. This leads to the possibility that the low number of Class VII RXTE observations we identify is due to a selection effect; we would not have been able to see this observation’s long-term Class VII-like behaviour if the observation had been shorter than $\sim 2$ ks.

Observation 12406 was performed within 24 hours of RXTE observation 96420-01-32-06, which showed Class IX variability. The lightcurve presented for Chandra OBSID 12406 shows a mean count rate (36 cts s*-1*), which is consistent with IGR J17091 being harder in this observation than in Observation 12505. This, combined with the lack of variability seen in its lightcurve, suggests that Observation 12505 is consistent with Class IX.

3.6 XMM-Newton

In Figure 24 we show lightcurves from two XMM-Newton observations. The lightcurve of XMM-Newton observation 0677980201, shown in the upper panel of Figure 24, shows the regular flares characteristic of Class IV variability. A simultaneous RXTE observation (OBSID 96420-01-05-000) also showed Class IV variability.

XMM-Newton observation 070038130, shown in the lower panel of Figure 24, was made after the end of RXTE observations IGR J17091-3624. As such it cannot be compared with contemporaneous RXTE data. The 5 s binned lightcurve shows no apparent variability, but a Fourier PDS of the observation (shown in Figure 25) reveals a QPO centred at around $\sim 0.15$Hz and a broad band noise component at lower frequencies. Drave et al. (2012) reported that IGR J17091 transited to the hard state in February 2012, seven months before this observation was taken. As such, we find that observation 0677980201 samples the hard state in IGR J17091 and is thus beyond the scope of our set of variability classes.

3.7 Suzaku

The two Suzaku observations of IGR J17091-3624 considered, OBSIDs 407037010 and 407037020, were performed during the 2nd and 3rd re-flares of the hard state phase of the 2011–2013 outburst. OBSID 407037010 was taken simultaneously with XMM-Newton observation 0700381301. The XIS 0 count rates are 7.8 cts s*-1* and 2.5 cts s*-1* respectively.

Neither lightcurve shows ‘heartbeats’ or any other type of GRS 1915-like variability. However, we find evidence of a low frequency QPO feature at $\sim$0.15 Hz in the OBSID 407037010; this QPO is also seen in XMM-Newton observation 0700381301 (Figure 25). The presence of a QPO below 1 Hz and flat-topped power density spectrum confirm that IGR J17091 was in the hard state at this time.

4 Discussion

Using observations from XMM-Newton, RXTE and Chandra, we describe the complex variability seen in IGR J17091 as a set of nine variability ‘classes’, labelled I to IX. These classes are distinguished from each other by values of upper and lower quartile (i.e. 25th and 75th percentile) count rates, mean RMS, the presence of QPOs in Fourier PDS, the shape of flare and dip features in the lightcurve and the presence of loops in the 6–16/2–6 keV hardness-intensity diagram HID1. See Section 3 for a full description of these classes.

The classification of some observations is clearer than others. Some orbits were too short to definitively quantify the behaviour of the source, whereas some other orbits contain a transition between two classes. An example lightcurve showing a transition from Class III to Class IV is presented in Figure 26.

Our set of classes is analogous to, but not based upon, the set of variability classes defined by Belloni et al. 2000 to describe the behaviour of the similarly complex LMXB GRS 1915. This ensures that our set of classes is not biased by an a priori assumption that the two objects are similar. However if we do assume that wide range of variability seen in these two objects are driven by the same physical processes, a direct comparison between the variability classes in the two systems can further our understanding of the physics that drive these exotic objects.

We also use all 2011-2013 IGR J17091-3624 data from RXTE, XMM-Newton, Chandra, Swift, INTEGRAL and Suzaku to analyse the long-term evolution of the 2011–2013 outburst. This in turn corresponds to all available X-ray data taken during this outburst.

4.1 Variability Classes: IGR J17091 vs. GRS 1915

As observations of IGR J17091 and GRS 1915 suffer from different values of interstellar absorption $N_{H}$, we cannot directly compare the absolute colours of these two objects. However, we can compare the evolution of colour both over time and as a function of count rate. We therefore use these parameters, along with power spectra and lightcurve morphology, when comparing GRS 1915 with IGR J17091.

For seven of our classes, we were able to assign the closest matching class described by Belloni et al. 2000 for GRS 1915 (see Table 5). We are unable to find analogues to our classes VII and VIII in observations of GRS 1915, and we suggest that these classes are unique to IGR J17091.

Below, we evaluate our mapping between GRS 1915 and IGR J17091 classes, and interpret the differences between each matched pair.

4.1.1 Classes I and II – Figures 5, 6

Classes I and II both show low count rates and little structure in their lightcurves. The two classes in GRS 1915 that also show this lightcurve behaviour are Class $\chi$888Note that, in GRS 1915+105, Class $\chi$ is further subdivided into four classes based on hard colour (Belloni et al., 2000; Pahari et al., 2013b). As we cannot obtain hard colour for IGR J17091, we treat $\chi$ as a single variability class here. and Class $\phi$. Belloni et al. 2000 differentiate between Classes $\phi$ and $\chi$ based on the hard colour (corresponding to $C_{2}$), as Class $\chi$ has a significantly higher value for this colour than Class $\phi$.

Data from RXTE  indicates that the transition from the hard state to the soft intermediate state between MJDs 55612 and 55615 (Drave et al., 2012). This was confirmed by a radio spectrum taken on MJD 55623 which was consistent with an observation of discrete ejecta (Rodriguez et al., 2011a). This observation of discrete ejecta at the transition between the hard state and the intermediate state has been reported in other LMXBS (e.g. XTE J1550-564, Rodriguez et al., 2003), and has also been associated with transitions to the $\chi$ Class in GRS 1915 (Rodriguez et al., 2008, see also review by Fender, 2006).

Using Fourier PDS, we conclude that Class I is analogous to Class $\chi$ in GRS 1915, while Class II is analogous to Class $\phi$. In Class $\chi$ observations of GRS 1915, broad band noise between $\sim 1-10$Hz and a QPO at around 5Hz are seen in the PDS. We find that both of these are present in Class I observations of IGR J17091. On the other hand, we find that Class $\phi$ observations of GRS 1915 do not show this broad band noise, and show either a weak ($q$-value $\lesssim 3$) QPO at $\sim 5$Hz or no QPO at all. We find that the weak QPO and lack of broad band noise are also seen in the PDS of Class II observations.

4.1.2 Classes III and IV – Figures 7, 10

Classes III and IV both show highly regular flaring activity in their lightcurves, but they differ in terms of timescale and pulse profile. As can be seen in lightcurves in Figure 10, flares in Class IV occur every $\sim 32$ s and are nearly identical to each other in shape. On the other hand, as can be seen in Figure 7, flares in Class III occur every $\sim 61$ s and may or may not end in a much faster sharp peak which is never seen in Class IV. In Figure 27 we show a two-dimensional histogram of flare peak count rate against flare duration, showing all flares in all observations classified as Class III or Class IV. In this figure, we can see that flares tend to group in one of two regions in count rate-duration space; a region between $\sim 90\mbox{--}110$  cts s*-1* PCU*-1and $\sim 35\mbox{--}55$ s, corresponding to flares seen in Class III, and a region between $\sim 150\mbox{--}250$  cts s-1* PCU*-1*and $\sim 20\mbox{--}55$ s, corresponding to flares seen in Class IV. From this plot, we conclude that the flares seen in Class III exist in a different population to the flares seen in Class IV.

The GRS 1915 classes that show behaviour most similar to these are $\rho$ and $\nu$; both produce similar structures in their lightcurve, but Class $\nu$ is differentiated from Class $\rho$ by the presence of a secondary count rate peak which occurs $\sim 5$ s after the primary (Belloni et al., 2000).

The secondary peak is present in most Class III observations and some Class IV observations (Figure 28), suggesting that both classes consist of a mix of $\rho$-like and $\nu$-like observations. However, the poor statistics sometimes make the presence of this secondary peak difficult to detect. As such, we do not use the presence or absence of this peak as a criterion when assigning classes. Instead we choose to separate Classes III and IV based on the larger-scale structure in their lightcurves (see Section 3.2.4). Due to the aforementioned difference in burst populations between the two classes, we suggest that classes III and IV do represent two distinct classes rather than a single class with a period that drifts over time. We suggest that Classes $\rho$ and $\nu$ in GRS 1915 could also be re-partitioned in this way.

However, HID1 loops are found to generally execute in an anticlockwise direction in Classes III and IV (previously noted by e.g. Altamirano et al., 2011a); the opposite direction to the clockwise loops in Classes $\rho$ and $\nu$ reported by e.g. Belloni et al., 2000 and repeated by us using the same method we apply to data from IGR J17091-3624 (see Section 2). This suggests that Classes III and IV could be generated by a different physical mechanism to Classes $\rho$ and $\nu$. Alternatively, Classes III and IV could be generated by the same mechanism as $\rho$ and $\nu$ if some other unknown process was able to alter the spectral evolution of flares in these classes.

4.1.3 Class V – Figure 12

The lightcurve of a Class V observation appears similar to that of a Class $\mu$ observation of GRS 1915, as both are characterised by rapid $\rho$-like flares which occur less regularly than in Class $\rho$. In addition to this, flares in Class $\mu$ fall into two clear populations, as do the flares in Class V. However, significant differences exist between Class V and Class $\mu$. Class $\mu$ observations are characterised by long ($\sim 100$ s) excursions to plateaus of high count rate, a behaviour which is not seen in any Class V observation thus far.

We note that the HID1 in Class V observations displays a loop in the clockwise direction; the opposite direction to the looping seen in Classes III and IV but the same direction seen in Class $\mu$.

Regarding the two-population nature of flares seen in this class (see Section 3.2.5), we suggest that V2 flares may simply be two V1 flares that occur close together in time, such that the second flare starts during the decay of the first flare. This would result in an apparent two-peaked flare structure, as we see in type V2 flares. This interpretation also accounts for the bimodal distribution of flare duarations shown in the 2D histogram of Figure 14, as this could be caused by the misinterpretation of two-flare V2 events as a single event. This also accounts for the Gaussian distribution of peak flare intensities seen in Figure 14), as the constituents of each V2 event would be from the same population as V1 flares.

4.1.4 Class VI – Figure 16

Class VI is dominated by long flaring periods which separate periods of low count rate, as can be seen in the lightcurve presented in Figure 16. Similar behaviour is seen in the lightcurves of observations of GRS 1915 belonging to Classes $\lambda$ and $\omega$ (Klein-Wolt et al., 2002). However, the long count rate ‘dips’ are far less regular in Class VI than in Classes $\lambda$ and $\omega$, and we also note long periods of medium count rate during which neither flares nor dips occur. This variability class is noted by Pahari et al. (2012) who suggest that this class is unique to IGR J17091999Pahari et al. (2012) refers to Class VI as Class C2.. However, Pahari et al. (2013a) show that, in a plot of burst decay time against burst rise time, Classes VI and $\lambda$ fall in a straight line, suggesting a similar physical origin for both.

While it is cetainly true that Class VI is not a perfect analogue of either Class $\lambda$ or Class $\omega$, Class VI only differs noticeably from Class $\lambda$ during the extended low-variability portions of its lightcurves. As such, we associate Class VI with Class $\lambda$.

4.1.5 Class VII – Figure 18

We are unable to find an analogue of Class VII in observations of GRS 1915. This class, and its apparent uniqueness, have previously been noted by Pahari et al., 2012101010Pahari et al. (2012) refers to Class VII as Class C1.. Pahari et al. found that the $C_{2}$ hard colour in this class increases during count rate dips and decreases during count rate peaks. Here we reproduced the results of Pahari et al. and found that the anti-correlation between hard-colour and intensity is not physical, but due to the definition of $C_{2}$: the count rate in band $L_{C}$ is approximately constant and consistent with background, and therefore $C_{2}=L_{C}/L_{A}\propto L_{A}^{-1}$, which will naturally anticorrelate with intensity.

Although a correlation between QPO frequency and count rate has been noted in the $\sim 5$Hz QPO seen in GRS 1915 (e.g. Markwardt et al., 1999; Vignarca et al., 2003), this QPO is also seen in Class VII observations at the same time as the $\sim 0.1$Hz QPO. As such, the flux-frequency relationship in the very low frequency ($\sim 0.1$Hz) QPO in Class VII is apparently unique amongst the classes of both IGR J17091 and GRS 1915.

4.1.6 Class VIII – Figure 20

We are unable to find an analogue of Class VIII in observations of GRS 1915. When it is flaring, the lightcurve waveform is similar to that seen in Class $\rho$, with rapid regular spikes in count rate. The lightcurve also shows irregular dips in count rate similar to those seen in Class VI and in Class $\lambda$ in GRS 1915.

However, the amplitude of the flares in Class VIII is much larger, and the frequency much higher, than in Classes VI or $\lambda$. The amplitude of the flares in Class VIII can approach $\sim 350$cts s*-1* PCU*-1*, while the flare separation time of 4–5 s makes Class VIII the fastest flaring activity seen in any class of IGR J17091 or GRS 1915. As such, we consider this variability class distinct from both Class VI and Class $\lambda$.

4.1.7 Class IX - Figure 21

Class IX is defined by long periods of high amplitude but unstructured variability (with a broad peaked noise component in the Fourier spectrum peaked at $\sim$0.3 Hz) punctuated with infrequent irregular short-duration ‘spikes’ in which the count rate increases by a factor of $\sim 2$–$3$. A similarity between this Class and Class $\gamma$ in GRS 1915 has been previously noted by Altamirano & Belloni (2012). However, the irregular spikes seen in some Class IX lightcurves are not reproduced in Class $\gamma$ lightcurves of GRS 1915.

4.2 General Comparison with GRS 1915+105

Overall, variability in IGR J17091 tends to be faster than structurally similar variability in GRS 1915, as can be noted in Classes III and IV compared to Classes $\rho$ and $\nu$ (see also Altamirano et al., 2011a). Additionally, IGR J17091 also displays highly structured variability unlike anything yet seen in GRS 1915, with classes VII and VIII in particular showing very fine detail in their lightcurves.

In total we find 2 variability classes which are seen in IGR J17091 but not in GRS 1915, compared with 8 that are seen in GRS 1915 but not in IGR J17091. As relatively little data exists on GRS 1915-like variability in IGR J17091, the presence of classes in GRS 1915 that are not seen in IGR J17091 could simply be an observational effect. It is unknown how long each variability class lasts for and, as such, additional variability classes could have occurred entirely while IGR J17091 was not being observed (however, see Huppenkothen et al., 2017 for a study on GRS1915 based on more than 16 years of data). However, GRS 1915 has displayed variability classes consistently since its discovery in 1992, implying that the two classes seen only in IGR J17091 are either completely absent in GRS 1915 or that they occur with a much lower probability. In either case, this implies physical differences between methods of generating GRS 1915-like variability in the two objects.

As noted in sections 4.1.1 to 4.1.7, variability classes seen in both IGR J17091 and GRS 1915 show differences in the different objects. In particular, we note the presence of irregular flares in Class IX which are not seen in the analogous Class $\gamma$. If these classes are indeed generated by the same processes in both objects, the differences between them must represent physical differences between the objects themselves.

It has previously been noted that, while the hardness ratios in IGR J17091 and GRS 1915 during $\rho$-like classes are different, the fractional hardening between the dip and peak of each flare is consistent with being the same in both objects (Capitanio et al., 2012). This suggests that the same physical process is behind the ‘heartbeats’ seen in both objects.

We note the presence of hysteretic HID1 loops in some classes of both objects. Although these loops are always clockwise in GRS 1915, they can be executed in either direction in IGR J17091. Classes in IGR J17091 that show loops all have a preferred loop direction: anticlockwise in Classes III and IV and clockwise in classes V, VI, VII and VIII. In cases where the loop direction was opposite to that expected for a given class, loop detections were generally only marginally significant. In particular, we note that Classes IV and V tend to show loops in opposite directions, despite the similarities between their lightcurves and the $\rho$, $\nu$ and $\mu$ classes in GRS 1915. The fact that IGR J17091 can show HID1 loops in both directions suggests that an increase in soft emission can either precede or lag a correlated increase in hard emission from IGR J17091. Whether soft emission precedes or lags hard emission is in turn is dependent on the variability class.

There are also non-trivial similarities between variability in the two objects. We note the presence of a $\sim 5$Hz QPO in many of the classes seen in IGR J17091, and this same 5Hz QPO is seen in lightcurves of GRS 1915. Similarly Altamirano & Belloni (2012) reported the discovery of a 66Hz QPO in IGR J17091; a very similar frequency to the 67Hz QPO observed in GRS 1915 (Morgan et al., 1997). It is not clear why these QPOs would exist at roughly the same frequencies in both objects when other variability in IGR J17091 tend to be faster.

4.3 Comparison with the Rapid Burster

In 2015, Bagnoli & in’t Zand (2015) discovered the existence of two GRS 1915-like variability classes in the neutron star binary MXB 1730-335, also known as the ‘Rapid Burster’. Specifically, Bagnoli & in’t Zand (2015) note the presence of variability similar to Classes $\rho$ and $\theta$ in GRS 1915.

Class $\theta$-like variability, seen in RXTE observation 92026-01-20-02 of the Rapid Burster, is not closely matched by any of the classes we identify for IGR J17091. However, the lightcurves of a Class $\theta$ observation feature large dips in count rate similar to those seen in Classes VI and VIII in IGR J17091.

Conversely, Class $\rho$-like variability is seen in all three objects. Bagnoli & in’t Zand (2015) note that the variability of the $\rho$-like flaring is slower in the Rapid Burster than in either GRS 1915 or IGR J17091. It has previously been suggested that the maximum rate of flaring in LMXBs should be inversely proportional to the mass of the central object (e.g. Belloni et al., 1997; Frank et al., 2002). In this case, the fact that variability is faster in IGR J17091 than in GRS 1915 could simply be due to a lower black hole mass in the former object (Altamirano et al., 2011a). However if variability in the Rapid Burster is assumed to be physically analogous to variability in these two black hole objects, then we note that a correlation between central object mass and variability timescale no longer holds.

4.4 Comparison with Altamirano et al. 2011

Altamirano et al. (2011a) identify 5 GRS 1915 variability classes in a subset of observations from the 2011-2013 outburst of IGR J17091: six of these observations are presented in Table 6 along with the best-fit GRS 1915 class that we assign it in this paper (see also Table 5).

We acknowledge differences between the classifications assigned by this paper and by Altamirano et al. (2011a). We ascribe these differences to the different approaches we have used to construct our classes. In particular while we have constructed an independent set of variability classes for IGR J17091 which we have then compared to the Belloni et al. classes for GRS 1915, Altamirano et al. applied the Belloni et al. classes for GRS 1915 directly to IGR J17091.

In general, the variability classes we find to be present in IGR J17091 are broadly the same as those noted by Altamirano et al. (2011a). We do not associate any class with Class $\alpha$ in GRS 1915, but we find examples of all of the other variability classes posited by Altamirano et al. to exist in IGR J17091.

Altamirano et al., 2011a noted the presence of an anticlockwise loop in the HID of ‘heartbeat’-like observations of IGR J17091, opposed to the clockwise loop seen in HID of $\rho$-class observations of GRS 1915. This is consistent with our finding that hysteretic loops in classes III and IV also tend to execute in an anticlockwise direction. However, we additionally find that hysteretic loops in classes V, VI, VII and VIII tend to execute in a clockwise direction. This is also different from GRS 1915, in which the loop is executed in the same direction in all classes. We also additionally report that clockwise loops tend to be more complex than anticlockwise loops seen in IGR J17091, with many showing a multi-lobed structure not seen in GRS 1915. This apparent inconsistency between the objects strengthens the suggestion in Altamirano et al., 2011a that the heartbeat-like classes in GRS 1915 and IGR J17091 may be generated by physically different mechanisms.

4.5 New Constraints on Accretion Rate, Mass & Distance

The constraints that Altamirano et al., 2011a placed on the mass and distance of IGR J17091 assumed that the object emitted at its Eddington luminosity at the peak of the 2011–2013 outburst. They report a peak 2–50 keV flux of $4\times 10^{-9}$ ergs s*-1* cm*-2* during flares in ‘heartbeat’-like lightcurves during this time. The correction factor $C_{Bol,Peak}$ to convert 2–50 keV flux to bolometric flux is not well constrained, but Altamirano et al., 2011a suggest an order-of-magnitude estimate of $\lesssim 3$, corresponding to a peak bolometric flux of $\lesssim 1.2\times 10^{-8}$ ergs s*-1* cm*-2*.

Maccarone, 2003 performed a study of the soft to hard transitions in 10 LMXBs. They found that all but one perform this transition at a luminosity consistent with between 1% and 4% of the Eddington limit. We use Swift observation 00031921058 taken on MJD 55965 to create a spectrum of IGR J17091 during the approximate time of its transition from a soft to a hard state (Drave et al., 2012). We fit this spectrum above 2 keV with a power-law, and extrapolate to find a 2–50 keV flux of $8.56\times 10^{-10}$ ergs s*-1* cm*-2*. Assuming that the transition bolometric correction factor $C_{Bol,Tran}$ is also $\lesssim 3$, this corresponds to a bolometric flux of $\lesssim 2.5\times 10^{-9}$ ergs s*-1* cm*-2*.

By comparing this with the results of Maccarone, 2003 and Altamirano et al., 2011a, we find that IGR J17091-3624 was likely emitting at no more than $\sim 5$–20% of its Eddington Limit at its peak. This number becomes $\sim 6\mbox{--}25$% if we instead use $C_{Bol,Tran}=2.4$, or $\sim 8\mbox{--}33$% if $C_{Bol,Tran}=1.8$. With this new range of values, we are able to re-derive the compact object mass as the function of the distance (Figure 29). We find that for a black hole mass of $\sim 10$M⊙, as suggested by Iyer et al., 2015a, IGR J17091 is within the galaxy at a distance of 6–17 kpc. This is consistent with the estimated distance of $\sim 11\mbox{--}17$ kpc estimated by Rodriguez et al., 2011a for a compact object mass of 10M⊙.

4.6 Implications for Models of ‘Heartbeat’ Variability

We have found that hysteretic HID loops can execute in both directions in IGR J17091 (e.g. Section 4.4), as well as found a revised estimate that IGR J17091 accretes at $\lesssim 20$% Eddington (Section 4.5). Both of these findings have implications for physical models of GRS 1915-like variability in this source.

Firstly, we find that Eddington-limited accretion is neither necessary nor sufficient for GRS 1915-like variability. The discovery of GRS 1915-like variability in the sub-Eddington Rapid Burster (Bagnoli & in’t Zand, 2015; Bagnoli et al., 2015) provided the first evidence that Eddington-limited accretion may not be a driving factor in this type of variability. We strengthen this case by finding that IGR J17091-3624 is also likely sub-Eddington. As such, we further rule out any scenario in which Eddington-limited accretion is required for GRS 1915-like variability in black hole LMXBs specifically.

Secondly, by using the direction of hysteretic HID loops, we find that hard photon lag in ‘heartbeat’-like classes of IGR J17091 can be either positive or negative. This could mean that we must rule out the causal connection between soft and hard emission being common to all classes.

In either case, we find that scenarios that require high global accretion rates or predict a consistent hard photon lag (e.g. Neilsen et al., 2011; Janiuk & Czerny, 2005), are not able to explain GRS 1915-like variability in IGR J17091 unless they also feature geometric obscuration in a subset of variability classes. We note that simulations by Nayakshin et al., 2000 require an Eddington fraction of $\gtrsim 0.26$ before GRS 1915-like variability, a value which falls in the range $\sim 0.05\mbox{--}0.33$ that we find for the peak Eddington fraction of IGR J17091.

In addition to being near its Eddington limit GRS 1915 also has the largest orbit of any known LMXB (e.g. McClintock & Remillard, 2006). Sa̧dowski, 2016 have also shown that thin, radiation dominated regions of disks in LMXBs require a large-scale threaded magnetic field to be stable, and the field strength required to stabilise such a disk in GRS 1915 is higher than for any other LMXB they studied. We suggest that one of these parameters is more likely to be the criterion for GRS 1915-like variability. If better constraints can be placed on the disk size and minimum stabilising field strength in IGR J17091, it will become clear whether either of these parameters can be the unifying factor behind LMXBs that display GRS 1915-like variability.

5 Conclusions

We have constructed the first model-independent set of variability classes for the entire portion of the 2011–2013 outburst of IGR J17091 that was observed with RXTE. We find that the data are well-described by a set of 9 classes; 7 of these appear to have direct counterparts in GRS 1915, while two are, so far, unique to IGR J17091. We find that variability in IGR J17091 is generally faster than in the corresponding classes of GRS 1915, and that patterns of quasi-periodic flares and dips form the basis of most variability in both objects. Despite this, we find evidence that ‘heartbeat’-like variability in both objects may be generated by different physical processes. In particular, while hard photons always lag soft in GRS 1915, we find evidence that hard photons can lag or precede soft photons in IGR J17091 depending on the variability class.

We also report on the long-term evolution of the 2011–2013 outburst of IGR J17091, in particular noting the presence of 3 re-flares during the later part of the outburst. Using an empirical relation between hard-soft transition luminosity and Eddington luminosity (Maccarone, 2003), we estimate that IGR J17091 was likely accreting at no greater than $\sim 33$% of its Eddington limit at peak luminosity.

We use these result to conclude that any model of GRS 1915-like variability which requires a near-Eddington global accretion rate is insufficient to explain the variability we see in IGR J17091. As such we suggest that an extreme value of some different parameter, such as disk size or minimum stabilising large-scale magnetic field, may be the unifying factor behind all objects which display GRS 1915-like variability. This would explain why sub-Eddington sources such as IGR J17091 and the Rapid Burster do display GRS 1915-like variability, while other Eddington-limited sources such as GX 17+2 and V404 Cyg do not.

Acknowledgements

J.C. and C.B. thank the Science & Technology Facilities Council and the Royal Astronomical Society for their financial support. D.A. thanks the Royal Society, and the International Space Science Institute (ISSI) for its support during the ‘The extreme physics of Eddington and super Eddington accretion onto Black Hole’ meetings (http://www.issibern.ch/teams/eddingtonaccretion/). T.B. also thanks the ISSI for their support. R.W. is supported by a NWO Top grant, Module 1. M. Pereyra gratefully acknowledges the Committee on Space Research (COSPAR) Capacity Building and Schlumberger Foundation Faculty for the Future Fellowship Programmes for their financial support. M.Pahari acknowledges the support of the UGC/UKIERI thematic partnership grant UGC 2014-15/02.

The authors acknowledge the use of AstroPy(Astropy Collaboration et al., 2013), APLpy (Robitaille & Bressert, 2012), NumPy (Jones et al., 2001) and MatPlotLib (Hunter, 2007) libraries for Python. We also acknowledge the use of public data from the RXTE (Bradt et al., 1993) archive, as well as FTOOLS (Blackburn, 1995) for data manipulation.

This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester.

Appendix A Flare-Finding Algorithm

The algorithm used to find flares is performed as such (see also Figure 30):

Choose some threshold values $T_{L}$ and $T_{H}$. Set the value of all datapoints below $T_{L}$ to zero. 2. 2.

Retrieve the x-co-ordinate of the highest value remaining in the dataset. Call this value $x_{m}$ and store it in a list. 3. 3.

Set the value of point at $x_{m}$ to zero. 4. 4.

Scan forwards from $x_{m}$. If the selected point has a nonzero value, set it to zero and move to the next point. If the selected point has a zero value, move to step (v). 5. 5.

Scan backwards from $x_{m}$. If the selected point has a nonzero value, set it to zero and move to the previous point. If the selected point has a zero value, move to step (vi). 6. 6.

Retrieve the y-co-ordinate of the highest value remaining in the dataset. Call this $y_{m}$. 7. 7.

If $y_{m}>T_{H}$, repeat steps (ii)–(vi). If $y_{m}<T_{H}$, proceed to step (viii). 8. 8.

Restore the original dataset. 9. 9.

Retrieve the list of $x_{m}$ values found in step (ii). Sort them in order of size. 10. 10.

For each pair of adjacent $x_{m}$ values, find the x-coordinate of the datapoint between them with the lowest y-value. Call these values $x_{c}$. 11. 11.

This list of $x_{c}$ can now be used to demarcate the border between peaks.

The values $T_{L}$ and $T_{H}$ can also be procedurally generated for a given piece of data:

Select a small section of the dataset or a similar dataset (containing $\sim 20$ peaks by eye) and note the location $x_{e}$ of all peaks found by eye. 2. 2.

Let $P_{L}$ and $P_{H}$ be two arbitrary values in the range $[0,100]$. 3. 3.

Let $T_{L}$ ($T_{H}$) be the $P_{L}$th ($P_{H}$th) percentile of the y-values of the subsection of dataset. 4. 4.

Run the flare-finding algorithm up to step (ix). Save the list of $x_{m}$. 5. 5.

Split the dataset into bins on the x-axis such as the bin width $b\ll p$, where $p$ is the rough x-axis separation between peaks. 6. 6.

For each bin, note if you found any value in $x_{m}$ falls in the bin and note if any value of $x_{e}$ falls in the bin. 7. 7.

Using each bin as a trial, compute the Heidke Skill Score (Heidke, 1926) of the algorithm with the method of finding peaks by eye:

[TABLE]

Where $A$ is the number of bins that contain both $x_{e}$ and $x_{m}$, $B$ ($C$) is the number of bins that contain only $x_{m}$ ($x_{e}$) and $D$ is the number of bins which contain neither (Kok, 2000). 8. 8.

Repeat steps (iii)–(vii) for all values of $P_{H}>P_{L}$ for $P_{L}$ and $P_{H}$ in $[1,100]$. Use a sensible value for the resolution of $P_{L}$ and $P_{H}$. Save the HSS for each pair of values 9. 9.

Locate the maximum value of HSS, and note the $P_{L}$ and $P_{H}$ values used to generate it. Use these values to generate your final $T_{L}$ and $T_{H}$ values.

We show an example of Heidke skill score grid for this algorithm, applied to a Class IV observation, in Figure 31.

Appendix B Model-Independent Classification of each Observation of IGR J17091-3624

Observation IDs, and orbit IDs, for every observation and observation segment that was used in our analysis are presented in Table 7. Note that not all of every observation was used; in many cases, large spikes caused by PCA PCUs switching off or on rendered $\sim 100$ s unusable. As these often occurred very close to the beginning or end of an observation segment, small sections of data before or after these spikes was also sometimes discarded. Every observation segment is presented along with the variability class assigned to it by this study.