Supergiant Fast X-ray Transients uncovered by the EXTraS project: flares reveal the development of magnetospheric instability in accreting neutron stars

TL;DR

This study analyzes X-ray flares from Supergiant Fast X-ray Transients using XMM-Newton data, revealing that their flaring activity aligns with magnetospheric instability models, and provides detailed flare characteristics.

Contribution

It offers a comprehensive analysis of SFXT flares using Bayesian methods, linking observed properties to magnetospheric instability theories and accretion processes.

Findings

144 flares identified across a luminosity range of 10^32 to 10^36 erg/s.

Flares' temporal properties match predictions from Rayleigh-Taylor instability models.

Rapidly rotating neutron stars may experience sporadic accretion from temporary disks.

Abstract

The low luminosity, X-ray flaring activity, of the sub-class of high mass X-ray binaries called Supergiant Fast X-ray Transients, has been investigated using XMM-Newton public observations, taking advantage of the products made publicly available by the EXTraS project. One of the goals of EXTraS was to extract from the XMM-Newton public archive information on the aperiodic variability of all sources observed in the soft X-ray range with EPIC (0.2-12 keV). Adopting a Bayesian block decomposition of the X-ray light curves of a sample of SFXTs, we picked out 144 X-ray flares, covering a large range of soft X-ray luminosities (1e32-1e36 erg/s). We measured temporal quantities, like the rise time to and the decay time from the peak of the flares, their duration and the time interval between adjacent flares. We also estimated the peak luminosity, average accretion rate and energy release in the flares. The observed soft X-ray properties of low-luminosity flaring activity from SFXTs is in qualitative agreement with what is expected by the application of the Rayleigh-Taylor instability model in accreting plasma near the neutron star magnetosphere. In the case of rapidly rotating neutron stars, sporadic accretion from temporary discs cannot be excluded.

Tables

Table 1. Table 1: Logbook of the X ​ M ​ M 𝑋 𝑀 𝑀 XMM - N ​ e ​ w ​ t ​ o ​ n 𝑁 𝑒 𝑤 𝑡 𝑜 𝑛 Newton  observations used here.

Target	Obs. ID	Start Date	Duration
			(ks)
IGR J08408-4503	0506490101	2007-05-29	45.7
IGR J11215-5952	0405181901	2007-02-09	22.2
IGR J16328-4726	0654190201 (a)	2011-02-20	21.9
IGR J16328-4726	0728560201 (b)	2014-08-24	36.2
IGR J16328-4726	0728560301 (c)	2014-08-26	23.0
IGR J16418-4532	0206380301 (a)	2004-08-19	23.2
IGR J16418-4532	0405180501 (b)	2011-02-23	39.6
XTE J1739-302	0554720101	2008-10-01	43.1
IGR J17544-2619	0148090501	2003-09-11	11.2
IGR J18410-0535	0604820301	2010-03-15	45.9
IGR J18450-0435	0306170401 (a)	2006-04-03	19.2
IGR J18450-0435	0728370801 (b)	2014-10-13	22.9
IGR J18483-0311	0694070101	2013-04-18	57.6

Table 2. Table 2: Source properties (see Sidoli & Paizis 2018 and referencees therein).

Name	Dist	$P_{\mathrm{orb}}$	orbital	$P_{\mathrm{spin}}$
	(kpc)	(d)	eccentricity	(s)
IGR J08408-4503	2.7	9.54	0.63	$-$
IGR J11215-5952	7.0$\pm 1.0$	164.6	$>$0.8	187
IGR J16328-4726	7.2$\pm 0.3$	10.07	$-$	$-$
IGR J16418-4532	13	3.75	0.0	1212
XTE J1739-302	2.7	51.47	$-$	$-$
IGR J17544-2619	3.0$\pm 0.2$	4.93	0.4	71.49(a)
IGR J18410-0535	3$\pm 2$	6.45	$-$	$-$
IGR J18450-0435	6.4	5.7	$-$	$-$
IGR J18483-0311	3.5$\pm 0.5$	18.52	$\sim$0.4	21.05

Table 3. Table 3: Properties of the SFXT flares. The asterisk (*) denotes “unresolved flares”.

Flare	Peak Luminosity	Energy	Waiting time $\Delta T$	Duration $\Delta$tf	Average Luminosity	$\delta$trise	$\delta$tdecay
ID	(1034 erg s-1)	(1037 erg)	(s)	(s)	(1034 erg s-1)	(s)	(s)
IGR J08408-4503
1	1.2$\pm 0.9$	$-$	$-$	$-$	$-$	$-$	240
2*	0.6$\pm 0.5$	0.29$\pm 0.22$	7900	470	0.6$\pm 0.5$	70	80
3	0.9$\pm 0.7$	1.2$\pm 0.4$	2700	3500	0.33$\pm 0.13$	700	2500
4	0.27$\pm 0.20$	$-$	14200	$-$	$-$	400	$-$
IGR J11215-5952
5	83$\pm 24$	31$\pm 7$	$-$	570	54$\pm 12$	520	15
6*	79$\pm 23$	27$\pm 8$	200	340	79$\pm 23$	17	190
7*	83$\pm 24$	12$\pm 3$	700	140	83$\pm 24$	120	32
8	112$\pm 32$	27$\pm 5$	600	304	87$\pm 18$	180	19
9*	114$\pm 33$	14$\pm 4$	180	125	114$\pm 33$	18	14
10	142$\pm 41$	38$\pm 8$	160	330	113$\pm 23$	40	227
11	132$\pm 38$	15$\pm 3$	360	130	115$\pm 24$	12	60
12*	123$\pm 35$	16$\pm 5$	210	132	123$\pm 35$	13	17
13	111$\pm 32$	24$\pm 5$	180	341	70$\pm 14$	26	212
14	25$\pm 7$	27$\pm 6$	2400	1480	18.4$\pm 3.7$	880	70
15	22$\pm 6$	24$\pm 5$	900	1900	12.7$\pm 2.6$	100	1400
16	34$\pm 10$	24$\pm 3$	3100	1400	17.2$\pm 2.5$	900	323
17	15$\pm 4$	23$\pm 5$	1500	2390	9.5$\pm 2.0$	130	1860
18*	10$\pm 3$	13$\pm 4$	3900	1220	10$\pm 3$	130	230
19	31$\pm 9$	27$\pm 4$	4900	2200	12.4$\pm 2.1$	1900	50
20*	35$\pm 10$	10$\pm 3$	300	290	35$\pm 10$	35	50
21	35$\pm 10$	$-$	400	$-$	$-$	50	$-$
IGR J16328-4726 (a)
22	9.1$\pm 0.8$	$-$	$-$	$-$	$-$	$-$	46
23	9.7$\pm 0.8$	8.8$\pm 0.5$	1100	2100	4.2$\pm 0.2$	39	1600
24	12.7$\pm 1.1$	7.0$\pm 0.3$	3400	920	7.6$\pm 0.4$	331	430
25	6.1$\pm 0.5$	4.3$\pm 0.3$	1200	1100	3.8$\pm 0.2$	70	700
26*	3.3$\pm 0.3$	1.3$\pm 0.1$	1800	410	3.3$\pm 0.3$	210	210
27	9.5$\pm 0.8$	10.1$\pm 0.6$	3500	1330	7.6$\pm 0.4$	190	650
28	15.6$\pm 1.3$	7.5$\pm 0.5$	1200	1290	5.8$\pm 0.4$	15	1170
IGR J16328-4726 (b)
29	39$\pm 3$	39$\pm 2$	$-$	1274	30.2$\pm 1.8$	28	867
30*	25$\pm 2$	32$\pm 3$	4300	1282	25.2$\pm 2.1$	43	77
31	74$\pm 6$	64$\pm 3$	3400	1800	34.9$\pm 1.8$	1450	230
32	102$\pm 9$	67$\pm 4$	1300	900	71.9$\pm 3.8$	500	253
33	40$\pm 3$	42$\pm 3$	1500	1180	35.1$\pm 2.6$	100	286
34	96$\pm 8$	210$\pm 10$	4100	3610	58.7$\pm 2.9$	990	1410
35	30$\pm 2$	28$\pm 2$	5800	1390	20.0$\pm 1.2$	32	870
36	19$\pm 2$	9.2$\pm 0.6$	2400	900	10.7$\pm 0.6$	80	700
37	147$\pm 12$	$-$	4790	$-$	$-$	750	$-$
IGR J16328-4726 (c)
38	24$\pm 2$	$-$	$-$	$-$	$-$	$-$	480
39	34$\pm 3$	20$\pm 1$	5500	860	23.5$\pm 1.5$	460	44
40	70$\pm 6$	29$\pm 2$	1300	540	54.1$\pm 2.9$	136	160
41*	32$\pm 3$	12$\pm 1$	1600	380	31.9$\pm 2.7$	90	70
42	101$\pm 8$	93$\pm 5$	1200	1703	54.6$\pm 3.0$	126	1009
43	121$\pm 10$	144$\pm 6$	4200	1870	76.8$\pm 3.3$	790	744
44*	33$\pm 3$	22$\pm 2$	2300	680	32.5$\pm 2.7$	130	80
45	43$\pm 4$	15.2$\pm 0.9$	700	519	29.2$\pm 1.8$	35	391
46	21$\pm 2$	$-$	1200	$-$	$-$	60	$-$

Table 4. Table 1: Properties of the SFXT flares. (continued)

Flare	Peak Luminosity	Energy	Waiting time $\Delta T$	Duration $\Delta$tf	Average Luminosity	$\delta$trise	$\delta$tdecay
ID	(1034 erg s-1)	(1037 erg)	(s)	(s)	(1034 erg s-1)	(s)	(s)
IGR J16418-4532 (a)
47	38$\pm 6$	$-$	$-$	$-$	$-$	$-$	298
48	79$\pm 12$	18$\pm 1$	1970	1020	17.6$\pm 1.4$	780	168
49	31$\pm 5$	22$\pm 2$	3800	1290	17.0$\pm 1.9$	940	60
50*	51$\pm 8$	10$\pm 2$	600	195	51$\pm 8$	19	32
51	73$\pm 11$	24$\pm 3$	400	390	62$\pm 8$	140	60
52*	73$\pm 11$	13$\pm 2$	300	170	73$\pm 11$	60	80
53*	72$\pm 11$	15$\pm 2$	500	210	72$\pm 11$	90	23
54*	68$\pm 11$	13$\pm 2$	600	180	68$\pm 11$	27	59
55*	78$\pm 12$	15$\pm 2$	300	190	78$\pm 12$	36	150
56	107$\pm 16$	21$\pm 2$	300	240	86$\pm 9$	31	140
57	132$\pm 20$	33$\pm 3$	460	360	92$\pm 8$	50	179
58	15$\pm 2$	$-$	1800	$-$	$-$	160	$-$
IGR J16418-4532 (b)
59*	75$\pm 12$	10$\pm 2$	$-$	130	75$\pm 12$	40	70
60*	115$\pm 18$	9$\pm 1$	160	80	115$\pm 18$	15	40
61	110$\pm 17$	40$\pm 4$	230	530	75$\pm 8$	50	370
62*	58$\pm 9$	1.5$\pm 0.2$	800	27	58$\pm 9$	15	12
63	90$\pm 14$	11$\pm 1$	1210	180	61$\pm 8$	80	120
64	405$\pm 62$	234$\pm 15$	1220	1105	212$\pm 14$	873	172
65	180$\pm 28$	30$\pm 3$	360	250	119$\pm 12$	50	164
66	58$\pm 9$	18$\pm 2$	1900	450	41$\pm 4$	180	60
67	68$\pm 10$	18$\pm 2$	600	390	46$\pm 5$	220	9
68	37$\pm 6$	16$\pm 2$	677	770	21$\pm 2$	31	610
69	128$\pm 20$	26$\pm 2$	1510	363	70$\pm 6$	282	16
70	110$\pm 17$	24$\pm 2$	210	324	74$\pm 7$	27	208
71*	11$\pm 2$	5.0$\pm 0.8$	3500	500	11$\pm 2$	400	300
72*	19$\pm 3$	2.6$\pm 0.4$	2900	130	19$\pm 3$	60	90
73*	21$\pm 3$	5.2$\pm 0.8$	500	240	21$\pm 3$	70	57
74*	14$\pm 2$	5.1$\pm 0.8$	800	400	14$\pm 2$	190	300
75	26$\pm 4$	15$\pm 1$	1600	1600	9.0$\pm 0.8$	40	1500
76	15$\pm 2$	25$\pm 2$	7200	3200	7.9$\pm 0.6$	2300	390
77	64$\pm 10$	17$\pm 2$	1300	420	41$\pm 5$	190	21
78*	50$\pm 8$	7$\pm 1$	200	130	50$\pm 8$	40	11
79	37$\pm 6$	5.2$\pm 0.6$	820	190	27$\pm 3$	80	40
80*	26$\pm 4$	9$\pm 1$	300	350	26$\pm 4$	120	42
81	89$\pm 14$	35$\pm 5$	800	460	76$\pm 11$	90	32
82	78$\pm 12$	10$\pm 1$	300	230	46$\pm 6$	50	130
83	43$\pm 7$	13$\pm 2$	1200	340	38$\pm 5$	15	110
84	152$\pm 23$	37$\pm 5$	700	311	118$\pm 16$	96	21
85	122$\pm 19$	45$\pm 4$	300	570	78$\pm 7$	60	399
86	103$\pm 16$	44$\pm 4$	2100	970	46$\pm 4$	660	150
87	90$\pm 14$	17$\pm 2$	900	240	72$\pm 8$	105	50
88	102$\pm 16$	18$\pm 3$	300	201	92$\pm 13$	31	32
89*	43$\pm 7$	8$\pm 1$	800	181	43$\pm 7$	21	37
90*	36$\pm 6$	5.2$\pm 0.8$	500	140	36$\pm 6$	60	70
91	78$\pm 12$	12$\pm 1$	410	180	65$\pm 8$	60	40
92	102$\pm 16$	$-$	300	$-$	$-$	18	$-$

Table 5. Table 1: Properties of the SFXT flares. (continued)

Flare	Peak Luminosity	Energy	Waiting time $\Delta T$	Duration $\Delta$tf	Average Luminosity	$\delta$trise	$\delta$tdecay
ID	(1034 erg s-1)	(1037 erg)	(s)	(s)	(1034 erg s-1)	(s)	(s)
XTE J1739-302
93	1.7$\pm 1.3$	0.58$\pm 0.31$	$-$	490	1.2$\pm 0.6$	280	12
94	0.6$\pm 0.4$	0.84$\pm 0.43$	1100	3100	0.28$\pm 0.14$	80	2800
95	1.7$\pm 1.2$	1.6$\pm 0.7$	11500	2370	0.66$\pm 0.29$	1960	100
96	3.9$\pm 2.9$	1.5$\pm 0.6$	500	1200	1.3$\pm 0.5$	110	1100
97*	0.9$\pm 0.6$	0.10$\pm 0.08$	1860	120	0.9$\pm 0.6$	50	60
98*	0.7$\pm 0.5$	0.47$\pm 0.35$	1900	710	0.7$\pm 0.5$	110	120
99*	0.7$\pm 0.6$	0.47$\pm 0.35$	1400	640	0.7$\pm 0.5$	80	130
100*	2.0$\pm 1.6$	0.053$\pm 0.039$	500	25	2.1$\pm 1.6$	12	8
101	0.3$\pm 0.2$	0.86$\pm 0.49$	3000	3900	0.22$\pm 0.12$	700	1700
102	0.6$\pm 0.4$	0.47$\pm 0.25$	11000	1650	0.28$\pm 0.15$	29	1290
103	4.1$\pm 3.1$	0.44$\pm 0.23$	5300	586	0.7$\pm 0.4$	535	14
104*	3.9$\pm 2.9$	0.78$\pm 0.57$	900	210	3.9$\pm 2.9$	18	60
105	4.2$\pm 3.1$	$-$	300	$-$	$-$	42	$-$
IGR J17544-2619
106	3.6$\pm 0.5$	$-$	$-$	$-$	$-$	$-$	110
107	6.4$\pm 0.8$	1.38$\pm 0.12$	700	680	2.0$\pm 0.2$	16	570
108	0.8$\pm 0.1$	$-$	5300	$-$	$-$	400	$-$
IGR J18410-0535
109	23 ${}_{-20}{}^{}{}_{}{}^{+40}$	60 ${}_{-50}{}^{}{}_{}{}^{+100}$	$-$	3810	16${}_{-9}{}^{}{}_{}{}^{+38}$	550	2169
110*	7 ${}_{-6}{}^{}{}_{}{}^{+12}$	1.2 ${}_{-1}{}^{}{}_{}{}^{+2}$	3300	170	7${}_{-6}{}^{}{}_{}{}^{+11}$	70	33
111*	7 ${}_{-6}{}^{}{}_{}{}^{+12}$	0.7 ${}_{-0.6}{}^{}{}_{}{}^{+1.2}$	290	90	7${}_{-6}{}^{}{}_{}{}^{+11}$	27	70
112	8 ${}_{-7}{}^{}{}_{}{}^{+14}$	1.0${}_{-0.9}{}^{}{}_{}{}^{+1.8}$	450	130	8${}_{-6}{}^{}{}_{}{}^{+15}$	36	70
113	12 ${}_{-10}{}^{}{}_{}{}^{+20}$	6${}_{-5}{}^{}{}_{}{}^{+10}$	150	729	8${}_{-6}{}^{}{}_{}{}^{+11}$	11	688
114*	7 ${}_{-6}{}^{}{}_{}{}^{+12}$	1.1 ${}_{-0.9}{}^{}{}_{}{}^{+1.9}$	910	170	7${}_{-6}{}^{}{}_{}{}^{+11}$	40	70
115	8 ${}_{-7}{}^{}{}_{}{}^{+14}$	7 ${}_{-6}{}^{}{}_{}{}^{+13}$	380	3390	2.2${}_{-1.9}{}^{}{}_{}{}^{+3.8}$	33	3280
116	0.7 ${}_{-0.6}{}^{}{}_{}{}^{+1.3}$	$-$	4500	$-$	$-$	260	$-$
IGR J18450-0435 (a)
117	51$\pm 16$	$-$	$-$	$-$	$-$	$-$	100
118	77$\pm 24$	25$\pm 5$	400	477	52$\pm 10$	72	237
119*	14$\pm 4$	1.3$\pm 0.4$	510	90	14$\pm 4$	50	150
120	43$\pm 14$	6.2$\pm 1.4$	830	210	31$\pm 7$	126	70
121*	67$\pm 21$	43$\pm 13$	500	637	67$\pm 21$	8	7
122	68$\pm 21$	7.8$\pm 1.8$	700	130	58$\pm 13$	65	50
123	84$\pm 26$	29$\pm 6$	390	370	77$\pm 17$	210	130
124*	90$\pm 28$	27$\pm 8$	500	290	90$\pm 28$	60	43
125	82$\pm 26$	44$\pm 8$	500	1100	39$\pm 7$	120	900
126*	8.4$\pm 2.6$	5.6$\pm 1.7$	2500	660	8.4$\pm 2.6$	140	90
127*	7.5$\pm 2.4$	5.7$\pm 1.8$	800	760	7.5$\pm 2.4$	150	250
128	13$\pm 4$	8.1$\pm 1.8$	1400	760	10.7$\pm 2.4$	510	80
129	15$\pm 5$	4.6$\pm 1.3$	400	720	6.5$\pm 1.8$	40	680
130	15$\pm 5$	19$\pm 4$	2490	2310	8.4$\pm 1.7$	1390	830
131*	13$\pm 4$	1.8$\pm 0.6$	2310	140	13$\pm 4$	50	60
132	32$\pm 10$	8.2$\pm 1.40$	2050	600	13.7$\pm 2.3$	500	43
IGR J18450-0435 (b)
133	5.7$\pm 0.6$	$-$	$-$	$-$	$-$	$-$	250
134	42$\pm 5$	30$\pm 4$	4100	2800	10.7$\pm 1.6$	2100	500
135	39$\pm 4$	15$\pm 3$	1800	520	28$\pm 6$	310	43
136	161$\pm 18$	49$\pm 9$	600	540	92$\pm 17$	43	380
137	23$\pm 3$	29$\pm 7$	1300	4600	6.3$\pm 1.4$	160	4100
138	165$\pm 19$	$-$	7200	$-$	$-$	400	$-$
IGR J18483-0311
139	4$\pm 1$	$-$	$-$	$-$	$-$	$-$	4100
140	3.1$\pm 0.9$	4.0$\pm 0.7$	9100	2450	1.65$\pm 0.30$	1740	90
141	2.9$\pm 0.8$	4.6$\pm 1.1$	3200	3620	1.28$\pm 0.31$	110	3340
142*	1.27$\pm 0.36$	1.8$\pm 0.5$	5200	1400	1.27$\pm 0.36$	260	500
143	4$\pm 1$	7.7$\pm 1.4$	4300	5900	1.31$\pm 0.24$	820	4800
144	1.2$\pm 0.3$	8.2$\pm 1.4$	8000	12100	0.68$\pm 0.11$	600	9900

Equations

$$S=\frac{2}{9}\frac{{\left(R_2-R_1\right)}^3}{{\left(\delta R_1\times \sqrt{T_1}+\delta R_2\times \sqrt{T_2}\right)}^2}$$

$$S\cong 2{\left(\frac{n}{3}\right)}^2\frac{R_2-R_1}{T_1+T_2}$$

$$E_j={\sum }_{B.b{.}_i}L(B.b{.}_i)\times dur(B.b{.}_i)$$

$$⟨L_j⟩=E_j/\Delta t_{f_j}$$

$$t_{\mathrm{rad}}\approx \frac{3k_{\mathrm{B}}T}{n_e\Lambda (T)}$$

$$t_{\mathrm{rad}}\approx 7\times 10^3[\mathrm{s}]\frac{R_9}{{\dot{M}}_{16}}f(u{)}_{\mathrm{rad}}\sim 700[\mathrm{s}]{\zeta }^{2/9}{\mu }_{30}^{2/3}{\dot{M}}_{16}^{-1},$$

$$R_{\mathrm{m}}\approx 10^9[\mathrm{cm}]{\zeta }^{4/81}{\mu }_{30}^{16/27}{\dot{M}}_{q,16}^{-6/27}$$

$$f(u{)}_{\mathrm{rad}}\approx 0.1{\zeta }^{14/81}{\mu }_{30}^{2/27}{\dot{M}}_{16}^{6/27}$$

$$Z\sim \alpha Ag{t}^2,$$

$$Z^{\prime }=Z-\int u(t)dt$$

$$u(t)=\frac{g{t}^2}{2t_{\mathrm{rad}}}\cos \chi$$

$$Z^{\prime }=\alpha Ag{t}^2-g\frac{t^3}{6t_{\mathrm{rad}}}.$$

$$\Delta T\approx 6\alpha A{t}_{\mathrm{rad}}\approx 0.18\left(\frac{\alpha }{0.03}\right)A{t}_{\mathrm{rad}}.$$

$$\Delta T\approx 130[\mathrm{s}]\left(\frac{\alpha }{0.03}\right)A{\zeta }^{2/9}{\mu }_{30}^{2/3}{\dot{M}}_{16}^{-1}.$$

$$v_b=a\sqrt{Ag\mathcal{R}}-g{t}^2/2t_{\mathrm{rad}}$$

$${\nu }_t=\frac{1}{3}{v}_t{l}_t$$

$${\lambda }_{max}=4\pi \left(\frac{{\nu }^2}{Ag}\right)$$

$$t_{k_{max}}\simeq 1.5(A{\alpha }_t^2{)}^{1/6}\sqrt{t_{\mathrm{ff}}{t}_{\mathrm{rad}}}\simeq 30[\mathrm{s}]{\zeta }^{4/27}{\mu }_{30}^{7/9}{\dot{M}}_{16}^{-2/3}.$$

$$\Delta t\sim t_{\mathrm{ff}}(R_{\mathrm{B}})=\frac{R_{\mathrm{B}}^{3/2}}{\sqrt{2G{M}_x}}\approx 400[\mathrm{s}]{\left(\frac{v_w}{1000[{\mathrm{km}\mathrm{s}}^{-1}]}\right)}^{-3}$$

$$\Delta M\approx 1.3\times 10^{18}[\mathrm{g}]\left(\frac{\alpha }{0.03}\right)A{\zeta }^{2/9}{\mu }_{30}^{2/3}.$$

$$\Delta E\approx 1.3\times 10^{38}[\mathrm{erg}]\left(\frac{\alpha }{0.03}\right)A{\zeta }^{2/9}{\mu }_{30}^{2/3},$$

$$⟨L⟩\equiv \frac{\Delta E}{\Delta t}\approx 3\times 10^{35}[{\mathrm{erg}\mathrm{s}}^{-1}]\left(\frac{\alpha }{0.03}\right)A{\zeta }^{2/9}{\mu }_{30}^{2/3}{v}_8^3.$$

$$\frac{\Delta E}{\Delta T}=10^{36}[{\mathrm{erg}\mathrm{s}}^{-1}]{\dot{M}}_{16}$$

$$\frac{{\dot{P}}_x}{P_x}=K_2\frac{{\mu }^2{P}_x}{4{\text{piup}}^{2}I{R}_{\mathrm{m}}^3}\simeq 2\times 10^{-12}(P_x/10s){\mu }_{30}^2{R}_{\mathrm{m},9}^{-3}$$

Full text

††pagerange: Supergiant Fast X-ray Transients uncovered by the EXTraS project: flares reveal the development of magnetospheric instability in accreting neutron stars–References††pubyear: 2019

1 Introduction

Supergiant Fast X–ray Transients (SFXTs) are a kind of high mass X–ray binaries (HMXBs) where a neutron star (NS) accretes a fraction of the wind of an early-type supergiant donor (see Sidoli 2017, Martínez-Núñez et al. 2017, Walter et al. 2015, for the most recent reviews). They were recognized as a new class of massive binaries thanks to rare, short and bright flares, reaching a peak luminosity LX$\sim$1036–1037 erg s*-1*, caught during $INTEGRAL$ observations (Sguera et al. 2005, 2006; Negueruela et al. 2006).

Low level X–ray flaring activity characterizes also their emission outside outbursts, down to LX$\sim$1032 erg s*-1*. Since SFXT flares usually display a complex morphology, it is somehow difficult to disentangle multiple flares (or structured flares) from the quiescent level and measure interesting quantities, such as, e.g., their duration, rise and decay times, time interval between flares. The investigation of SFXT flares requires a twofold approach: firstly, high throughput and uninterrupted observations are needed, both to detect flares even at very low X–ray fluxes and to measure the flare timescales without data gaps (within each observation) that might bias them; secondly, an efficient, systematic procedure to pick out flares and determine their observational properties, to be compared with the theory. Both conditions are met by the database of products made available to the community by the EXTraS project.

EXTraS (acronym of “Exploring the X–ray Transient and variable Sky") is a project funded within the EU/FP7 framework (De Luca et al., 2017), aimed at extracting from the $XMM$-$Newton$ public archive the temporal information (periodic and aperiodic variability) of all sources observed by the EPIC cameras in the 0.2–12 keV energy range.

In the first part of the paper we report on the behaviour of some essential flare quantities that we have extracted from the EXTraS database, by means of a Bayesian blocks analysis of the light curves of a sample of SFXTs.

In the second part of the paper, we discuss the behaviour of the flare properties and discovered dependences in terms of the interchange instability of accreting matter near NS magnetosphere. We show that for slowly rotating NSs, the development of Rayleigh-Taylor instability in a quasi-spherical shell above NS magnetosphere can qualitatively describe the observed properties of the SFXT flares. For rapidly rotating NSs accreting from the stellar wind, the propeller mechanism could lead to the formation of an equatorial dead disc that may trigger the magnetospheric instability once the centrifugal barrier at its inner edge is overcome.

This paper is organized as follows: Sect. 2 introduces the EXTraS project, for what is relevant here (for more details, we refer the reader to De Luca et al. 2017); Sect. 3 reports on the $XMM$-$Newton$ observations we have considered in our study; Sect. 4 reports on the automatic procedure used to pick out flares and to measure the parameters of the flares, together with the observational results. Sect. 5 discusses the results in the context of the interchange instability near the NS magnetosphere. In Sect. 6, we summarize our findings and conclusions.

2 EXTraS

The EXTraS project dived into the public soft X–ray data archive of $XMM$-$Newton$, building on top of the 3XMM-DR4 catalog (Rosen et al., 2016), and aimed at characterising the variability of as many EPIC point sources as possible. It provided the community with tools, high-level data products, catalogues, documentation, all available through the EXTraS web-site111http://www.extras-fp7.eu and the LEDAS astronomical data archive222https://www88.lamp.le.ac.uk/extras/archive. One of these products is a set of adaptively binned light curves, one for each exposure (observation segment specific to a single EPIC camera), particularly suited for the systematic analysis carried out in this paper. In particular, these light curves make it possible to identify easily flares and provide some parameters (like the slope between two blocks, see below) useful to their characterisation.

EXTraS addressed the characterisation of variability of point-like X-ray sources under various aspects, each one dealing with specific problems and dedicated techniques (De Luca et al., 2016). In particular, aperiodic variability has been treated separately on the short term (within a single observation segment or exposure) and on the long term (combining different pointings to the same source, including slew data). We will focus in this paper on short-term variability products. Dealing with separated exposures avoids the problem of combining data collected with different EPIC cameras, with different filters and operating modes. Although EXTraS provides also a characterisation of aperiodic variability in the frequency domain, we use here only light curves, in the form of count rates versus time. All the light curves considered in this paper always remain in the Poisson regime, and are not affected by significant instrumental pile-up.

An important aspect of the EXTraS approach is the technique of background modelling and subtraction. EPIC data can be affected for a large fraction of the observing time by strong and rapid background flares due to soft protons (up to 35%, Marelli et al. 2017). Standard recommendations for data screening lead to discard 21% of the data (Rosen et al., 2016).

EXTraS instead takes care of modelling the background, by disentangling its steady and variable components, and evaluating separately their distribution on the detector. Background maps for both components are built in order to take into account the background distribution on the detector. In this way it is possible to effectively subtract the variable EPIC background from source regions and use all the exposure time (see Marelli et al. 2017 for more details). The source region has a circular shape, while the background region covers most of the detector, excluding circular regions around contaminating sources. All radii are chosen to maximise the signal to noise ratio while minimising the contamination from other sources (the region files are provided through the EXTraS archive).

One of the algorithms implemented within EXTraS for the characterisation of short term (within a single, uninterrupted exposure) variability is based on Bayesian blocks (hereafter B.b., Scargle 1998; Scargle et al. 2013). This is a segmentation technique, often applied to astronomical time series, which aims to split the data into the maximum number of adjacent blocks, in such a way that each block is statistically different from the next. The B.b. algorithm starts from a fine segmentation of a time series in cells, that are subsequently merged in a statistically optimal way. We define a time cell as containing at least 50 source photons or 50 background photons expected in the source region. This algorithm for defining the initial segmentation allows us to identify features both in the source and background light curves, which can be considered constant within each cell. As a figure of merit we consider the likelihood of an average net source rate within each block, reduced by a fixed amount. This represents a cost for each added block, implementing an Occam’s razor: the higher the cost, the more significant the difference between two neighboring blocks. EXTraS provides two sets of Bayesian blocks light curves (nominally at 3$\sigma$ and 4$\sigma$, respectively): one with a lower cost, more sensitive to small variations in rate; the other with an higher cost, more robust in its segmentation. In this work we use the latter.

If we consider two neighboring blocks, we can be confident that the rate of the source has changed between blocks, while it is consistent with a constant within each block. However, this is not sufficient for us to tell whether the rate of the source has changed sharply or smoothly and to this extend we introduce a parameter that we call slope (S). This is the minimum rate of change in the counts rate of the source between two neighboring blocks. To find S, we shrink each of the blocks until their associated rates, $R_{1}$ and $R_{2}$, are compatible within 3$\sigma$, assuming that the uncertainty in the rates, $\delta R_{1}$ and $\delta R_{2}$, decreases with time as $T^{-\frac{1}{2}}$, as expected for Poisson events. Then, we assume that the rate of the source has changed linearly for the duration of the 2 blocks, $T_{1}+T_{2}$, compatibly with the 2 rates, and obtain:

[TABLE]

For similar blocks that are $n\sigma$ apart (as expected from a source that undergoes a linear trend in flux, with no background flares), this relation reduces to:

[TABLE]

.

3 $XMM$-$Newton$ observations

We searched for the SFXTs in the EXTraS database, selecting the B.b. EPIC light curves from observations pointed on the following members of the class: IGR J08408-4503, IGR J11215-5952, IGR J16328-4726, IGR J16418-4532, XTE J1739-302, IGR J17544-2619, IGR J18410-0535, IGR J18450-0435 and IGR J18483-0311. Multiple pointings were available for three sources, as reported in Table 1, where we list the observations considered in our study.

The B.b. light curves extracted from observations performed after 2012 are not present in the public EXTraS database, but were produced for this work by the team, using the same techniques explained in Sect. 2. For the sake of completness, note that we excluded from this investigation two observations present in the EXTraS database and targeted on IGR J16479-4514 and on IGR J18483-0311, because no flares were present (according to the definition assumed below).

All pointings are Guest Observer observations, except two Target of Opportunity ones (targeted on IGR J11215-5952 and the 2011 observation targeted on IGR J16418-4532) which were triggered at the occurrence of an outburst.

We have considered only the EPIC pn light curves (all in Full Frame mode), except for IGR J11215-5952 and the second observation of IGR J18450-0435, where the EPIC MOS light curves were considered. The reason is because in these observations EPIC pn operated in Small Window mode, where the background cannot be appropriately treated using the EXTraS techniques.

A systematic analysis of the $XMM$-$Newton$ observations considered in this work, but focussed on the X-ray spectroscopy, can be found in Giménez-García et al. (2015), Bozzo et al. (2017), Pradhan et al. (2018) (and references therein).

We list the source properties in Table LABEL:tab:prop, adopting the same values published by Sidoli & Paizis (2018). For sources with no published range of variability of the distance (the ones with no uncertainty present in Table LABEL:tab:prop), an error of $\pm{1}$ kpc has been assumed, when needed.

4 Results

In Fig. 1 we show the comparison between a SFXT light curve with a uniform binning and the one obtained with a B.b. segmentation. In Appendix (Figs. 11 and 12) the B.b. light curves of the other SFXTs in our sample are reported.

The adoption of the B.b. segmentation of the light curves offers an efficient and systematic way to select a flare without the need of assuming any specific model for its profile. Indeed, we have considered a “flare" as a statistically significant peak (i.e., a B.b. containing a local maximum), with respect to the surrounding, adjacent emission. We also used the local minima (“valleys") to calculate the duration of each flare (see below). This automated procedure picked out 144 SFXT flares.

4.1 Measuring observational quantities

After the selection of the local maxima and minima in the light curves, we have estimated the following temporal quantities (which are also explained graphically in Fig. 2, for clarity):

• •

Waiting time ($\Delta T$): we define it as the time interval between the peaks of subsequent flares. Given our B.b. segmentation of the SFXT light curve, for each B.b. containing a flare peak (a local maximum), the waiting time is calculated as the difference between the midtime of the B.b. containing the peak of the flare, and the midtime of the B.b. containing the peak of the previous flare. For all the first flares of the light curves the waiting time could not be calculated. The observations are uninterrupted so that, within each EPIC observation, $\Delta T$s do not suffer any bias, as well as other parameters in Fig. 2. In Fig. 2, $\Delta T$=[t2 + 0.5(t3 – t2) – t0].

• •

Rise time ($\delta t_{rise}$) to the peak of the flare: for unresolved flares (those made of a single B.b.), the rise time is calculated as $\delta t_{rise}$= $\Delta$R/S, where $\Delta$R is the difference in count rate between the two adjacent B.b. (the one containing the flare, and the B.b. located immediately before it), and “S" is the positive slope measured before the peak (in units of counts s*-1* ks*-1*). For the resolved flares (those spanning more than one B.b.), it is calculated as the time interval between the stop time of the B.b. containing the valley before the flare, and the start time of the block containing the peak of the flare. In Fig. 2, $\delta t_{rise}$=t2-t1. These definitions return the best estimates for the flare rise time. However note that in case of unresolved flares, the true rise time to the peak might be formally (although unphysically) zero. On the other hand, the rise time defined above for resolved flares can be considered a minimum value, by definition. The same is valid for the decay time, as defined below.

• •

Decay time ($\delta t_{decay}$) from the peak of the flare: for the unresolved flares, it is calculated as $\delta t_{decay}$= $\Delta$R/S, where $\Delta$R is the difference in count rate between the two adjacent B.b. (the one containing the peak, and the B.b. just after it), and “S” is the (negative) slope measured at the peak (in units of counts s*-1* ks*-1*). Note that we have considered the absolute value of the decay times. For the resolved flares (those spanning more than one B.b.), it is calculated as the time interval between the stop time of the block containing the peak of the flare and the start time of the block containing the valley next to it. In Fig. 2, $\delta t_{decay}$=t5-t3.

• •

Flare Duration ($\Delta$tf): it is defined as the time interval comprised between the two local minima surrounding the flare, subtracting the time interval covered by the B.b. containing these same minima. Therefore, it can be calculated as the sum of the time intervals covered by the B.b. in-between two local minima. For flare peaks which do not lie between two local minima (this might occur at the beginning and/or at the end of an observation), the total duration cannot be calculated. In Fig. 2, $\Delta$tf=t5-t1. The only exception to this rule is the first peak in the light curve of the SFXT IGRJ18410–0535: since its profile is very well defined (fast rise and exponential decay, hereafter FRED), we could measure its duration, although formally a valley is not present before it.

Besides the above timescales, we have calculated the following quantities, for each flare $j$:

• •

Flare peak Luminosity (Lj): the EPIC count rate of the B.b. containing the peak of the flare has been converted to unabsorbed flux (1–10 keV) using WebPIMMS333https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3pimms/w3pimms.pl assuming a power law spectrum with a photon index $\Gamma$=1 and a column density NH=1.5$\times$1022 cm*-2*. 444This spectral shape represents an average between somehow harder X-ray emission observed in SFXT bright flares (e.g. Sidoli et al. 2007), and softer emission from outside outbursts (Sidoli et al., 2008). The energy range 1–10 keV has been assumed to better compare with the literature.

This resulted in a conversion factor of 10*-11* erg cm*-2* count*-1* (EPIC pn). Then, the flare peak X-ray luminosities have been calculated assuming the source distances reported by Sidoli & Paizis (2018) and listed in Table LABEL:tab:prop, for clarity. The luminosity is also reported in the y-axis, on the right side of the graphs in Figs. 11 and 12, to enable a proper comparison between different sources.

• •

Energy (Ej) released in a flare $j$: it is calculated summing the products (L (B.b.i)$\times$ dur(B.b.i), where L(B.b.i) is the luminosity reached by a single block $i$, while dur(B.b.i) is the time duration of the block $i$) over the blocks covered by a single flare:

[TABLE]

• •

**Average Luminosity ($\langle{L}_{j}\rangle$) during a flare $j$: **

For each flare $j$, we calculate the average X-ray luminosity as:

[TABLE]

where $\Delta~{}t_{f_{j}}$ is the total duration of the flare $j$ (as defined before).

The last quantity relevant for this investigation is the pre-flare X–ray luminosity (or accretion rate; L${}_{X_{q}}$). With this term we mean the X–ray luminosity level displayed by the valley just before each flare.

We report in Appendix the Table 1 with the values of the flare quantities defined above. We collected 144 SFXT flares, which can be “unresolved" or “resolved", depending on whether they extend on just one or more B.b., respectively. About one third of the flares are unresolved. We have clearly marked them with an asterisk in Table 1.

5 Discussion

The automated procedure we have adopted to select the SFXT flares from $XMM$-$Newton$ light curves has lead to the measurement of observational quantities which can be now compared with the theory. In Sect. 5.1 we will discuss the phenomenology of the flares in the framework of the development of the Rayleigh-Taylor instability (RTI) in accreting plasma trying to enter the magnetosphere of slowly rotating neutron stars.

In order to perform this comparison, we will highlight the behaviour of flares from single sources, within the global behaviour of SFXT flares taken as a whole. In this way it is somehow possible to identify trends which are not found in flares occurred in a single source, but are due to the superposition of flares from different sources, lying in different regions of the parameter space, and vice versa: trends in single sources might be in principle mixed up when all SFXT flares are considered together.

In Fig. 3 we show the overall distributions of flare temporal properties (flare duration, waiting times, rise and decay times), while in Fig. 4 (left panel), the decay time is plotted against the rise time, to investigate the flare shape. In order to quantify the degree of asymmetry in the flare profiles, we calculated the parameter $\xi$ = ($\delta t_{decay}$ - $\delta t_{rise}$) / ($\delta t_{decay}$ + $\delta t_{rise}$). This implies that flares where the rise to the peak is much faster than the decay have \xi$$\sim1 (as in FRED-like profiles), while flares with a much slower rise than the decay show \xi$$\sim-1.

The result is reported in the right panel of Fig. 4, where it is evident that all range of values is covered by the SFXT flares analysed here.

We show in Fig. 5 the histograms of the peak luminosity and of the energy emitted during flares. From the latter histogram, a bimodality in the energy released in flares might be present, above and below $\sim$2$\times$1037 erg, but the relatively low statistics do not permit to draw a firm conclusion. Note that most of the flares with energies below 1037 erg are contributed by the source XTE J1739–302.

In Figs. 6,  7,  8 and 9 we show flare global behaviours which will be compared with the theory in the next subsections. From Fig. 6 an apparent anticorrelation is present between the waiting time between two adjacent flares and the X-ray luminosity in-between flares. The plot of energy released in flares versus the flare duration might indicate a positive trend, while the plot of the rise time of all flares, versus their pre-flare luminosities does not apparently show any correlation. A correlation is shown over several orders of magnitude by the ratio between the flare energy and the waiting time versus the pre-flare luminosity. Similar plots are reported in Figs. 13,  14,  15 and 16, where the flares from single sources are highlighted. They are discussed in Sect. 5.2.

We note that the same $XMM$-$Newton$ observations reported here have been analyzed by Giménez-García et al. (2015), Bozzo et al. (2017), Pradhan et al. (2018) (and references therein). Although their temporal-selected spectra of flares were extracted from intervals much longer than the B.b.s adopted here (so that they cannot be directly compared), some spectra of SFXT flares showed a column density larger than 1.5$\times$1022 cm*-2* (the value we assume for all flares), implying, in principle, a larger conversion factor from count rate to flux. However, even a very high absorption of $\sim$2$\times$1023 cm*-2* would imply only a $\sim$3 times larger unabsorbed flux (1-10 keV), with no impact on the conclusions of our work, where observational facts and theory are compared over four orders of magnitude in X-ray luminosity (and emitted energy).

5.1 Rayleigh-Taylor instability

It has long been recognized that plasma entry in NS magnetosphere in accreting X-ray binaries occurs via interchange instability – Rayleigh-Taylor (RTI) in the case of slowly rotating NSs (Arons & Lea, 1976; Elsner & Lamb, 1977) or Kelvin-Helmholtz (KH) in rapidly rotating NSs (Burnard et al., 1983). In the case of disc accretion, the plasma penetration into magnetosphere via RTI was compellingly demonstrated by multi-dimensional numerical MHD simulations (Kulkarni & Romanova, 2008). However, global MHD simulations of large NS magnetospheres ($\sim 10^{9}$ cm) have not been performed yet, and information about physical processes near NS magnetospheres should be inferred from observations.

During quasi-spherical wind accretion onto slowly rotating NSs, there is a characteristic X-ray luminosity $L^{*}\simeq 4\times 10^{36}$  erg s*-1* that separates two physically distinct accretion regimes: the free-fall Bondi-Hoyle supersonic accretion occurring at higher X-ray luminosity, when the effective Compton cooling time of infalling plasma is shorter than the dynamical free-fall time (Elsner & Lamb, 1984), and subsonic settling accretion at lower luminosities, during which a hot convective shell forms above the NS magnetosphere (Shakura et al., 2012, 2018). In the latter case, a steady plasma entry rate is controlled by plasma cooling (Compton or radiative) and is reduced compared to the maximum possible value determined by the Bondi-Hoyle gravitational capture rate $\dot{M}_{B}$ from the stellar wind by a factor $f(u)^{-1}\approx(t_{\mathrm{cool}}/t_{\mathrm{ff}})^{1/3}>2$.

The necessary conditions for settling accretion are met at low-luminosity stage in SFXTs. Low X-ray luminosities make it difficult to detect X-ray pulsations and therefore to answer the question of where the observed X-rays are actually produced. They can be either generated near the NS surface (if the inefficient plasma entry rate into the magnetosphere is provided by diffusion, cusp instabilities, etc., as discussed e.g. by Elsner & Lamb 1984), or be a thermal emission from magnetospheric accretion, like in the model developed for $\gamma$ Cas stars (Postnov et al., 2017). It is quite possible that at low-luminosity states of SFXTs, no RT-mediated plasma penetration into the NS magnetosphere occurs at all. This may be the case if the plasma cooling time is longer than the time a plasma parcel spends near the magnetosphere because of convection: $t_{\mathrm{cool}}>t_{\mathrm{conv}}\sim t_{\mathrm{ff}}(R_{\mathrm{B}})\sim 300-1000$ s. Once this inequality is violated, RTI can start to develop.

At this stage, the magnetospheric instability can occur for different reasons. For example, it was conjectured (Shakura et al., 2014) that bright flares in SFXTs are due to sudden break of the magnetospheric boundary caused by the magnetic field reconnection with the field carried along with stellar wind blobs. This can give rise to short strong outbursts occurring in the dynamical (free-fall) time scale during which accretion rate onto NS reaches the maximum possible Bondi value from the surrounding stellar wind.

Another reason for the instability can be due to stellar wind inhomogeneities which can disturb the settling accretion regime and even lead to free-fall Bondi accretion episodes.

At typical X-ray luminosities $L_{q}\simeq 10^{33}-10^{34}$ erg s*-1*, SFXTs occasionally demonstrate less pronounced flares with phenomenology as described in previous Sections. Below we present a possible scenario of development of such flares based on the consideration of non-linear growth of RTI. During these flares, an RTI-mixed layer is advected to the NS magnetosphere, but X-ray power generated is still below $\sim 10^{36}$  erg s*-1* to provide effective plasma cooling for steady RTI. In a sense, the observed short flares during low-luminosity SFXT state are due to ‘failed’ RTI.

These will enable us to explain, without making additional assumptions, the main observed properties of SFXT flares inferred from the statistical analysis presented in this paper. In all numerical estimates below, we assume the NS mass $M_{x}=1.5M_{\odot}$ and normalize the NS magnetospheric radius as $R_{\mathrm{m}}=10^{9}[\mathrm{cm}]R_{9}$, the mass accretion rate onto NS as $\dot{M}_{x}=10^{16}[\mathrm{g\,s}^{-1}]\dot{M}_{16}$ and the NS magnetic moment as $\mu=10^{30}[\mathrm{G\,cm}^{3}]\mu_{30}$.

5.2 Non-linear RTI growth

At the settling accretion stage, in a subsonic convective shell around the NS magnetosphere, external wind perturbations gravitationally captured from stellar wind of the optical components at the Bondi radius $R_{\mathrm{B}}=2GM_{x}/v_{w}^{2}$ travel down to the NS magnetosphere $R_{\mathrm{m}}$ with the convective motions. Therefore, the response time of the magnetosphere to the external perturbations is not shorter than about free-fall time from the Bondi radius, $t_{\mathrm{ff}}(R_{\mathrm{B}})=\sqrt{R_{\mathrm{B}}^{3}/2GM_{x}}$, ($R_{\mathrm{B}}=2GM_{x}/v^{2}\simeq 4\times 10^{10}[\mathrm{cm}]v_{8}^{-2}$  cm is the Bondi radius for the relative wind velocity $v=10^{8}[\mathrm{cm\,s}^{-1}]v_{8}$), typically a few hundred seconds for the stellar wind velocity from OB-supergiant $v_{w}\sim 1000$  km s*-1*.

In an ideal case with constant boundary conditions, the development of RTI occurs via production of a collection of bubbles (or, rather, flutes) with different size (Arons & Lea, 1976), and the mean plasma entry rate $u$ into magnetosphere is determined by the slowest linear stage of RTI in a changing effective gravity acceleration determined by plasma cooling (Shakura et al., 2012, 2018): $u\approx f(u)v_{\mathrm{ff}}(R_{\mathrm{m}})$, where $v_{\mathrm{ff}}(R_{\mathrm{m}})=\sqrt{2GM_{x}/R_{\mathrm{m}}}$ is the free-fall velocity at the magnetospheric radius.

In the quiescent state, the mass accretion rate is determined by the average plasma velocity $u$ at the magnetospheric boundary, $\dot{M}_{z}=4\pi R_{\mathrm{m}}^{2}\rho u$. At low X-ray luminosity $\lesssim 10^{35}$ erg s*-1*, radiative plasma cooling dominates the Compton cooling. The characteristic plasma cooling time is

[TABLE]

where $n_{e}$ is the electron number density, $\Lambda(T)\approx 2.5\times 10^{-27}\sqrt{T}$ [erg cm3 s*-1*] is thermal cooling function dominated by bremsstrahlung at the characteristic temperatures of the problem (1-10 keV). Taking into account that the plasma temperature at the basement of the shell near the magnetosphere is about the adiabatic value, $T\approx 2/5(GM_{x})/{\cal R}R_{\mathrm{m}}\simeq 10^{10}R_{9}$ K, and by expressing $n_{e}$ from the mass continuity equation, we obtain the radiative cooling time

[TABLE]

In the last equation, the magnetospheric radius and factor $f(u)_{\mathrm{rad}}$ are derived as (Shakura et al., 2013, 2018):

[TABLE]

[TABLE]

(the dimensionless parameter $\zeta\lesssim 1$ characterizes the size of the RTI region in units of the magnetospheric radius $R_{m}$).

Increase in the density $\rho_{m}$ near the magnetosphere would shorten the plasma cooling time and lead to an increase in the X-ray photon production from the NS surface, which in turn would enhance the Compton plasma cooling and increase the plasma entry rate $f(u)$. This would result in an X-ray flare (or a collection of flares) on top of the quiescent X-ray luminosity level.

5.2.1 Flare waiting time

We start with the estimate of the flare waiting time $\Delta T$. Consider the spreading of the RTI layer during the instability development. If there is no plasma penetration into the NS magnetosphere, the thickness of RTI mixing layer at the late non-linear stage grows as

[TABLE]

where $g=GM_{x}/R_{\mathrm{m}}^{2}$ is the gravity acceleration, the dimensionless factor $\alpha\sim 0.03$, $A\lesssim 1$ is the effective Atwood number (see e.g. Carlyle & Hillier (2017) for a recent discussion of numerical calculations of the non-linear growth of magnetic RTI).

However, in our problem the layer struggles against the mean plasma flow with velocity $u(t)$ determined by the plasma cooling that further slows down the RTI, and therefore the net distance the RTI layer extends above the magnetosphere is

[TABLE]

As long as the time is shorter than the cooling time, $t<t_{\mathrm{rad}}$, during the linear stage of RTI development in the unstable region we can write (Shakura et al., 2012):

[TABLE]

(here $\chi$ is the latitude from the magnetosphere equator where RT modes are the most unstable; below we set $\cos\chi=1$). Therefore, Eq. (8) takes the form:

[TABLE]

With time, the second (negative) term in Eq. (10) overtakes the first (positive) one. The negative value of the RTI layer height above magnetosphere $Z^{\prime}<0$ would inhibit instability growth because there will be no room for plasma flutes to interchange with magnetic field above the magnetospheric boundary (more precisely, above the layer in which the mean plasma entry rate is sustained for a given plasma cooling rate). Therefore, the growth of the RTI mixing layer size at the nonlinear stage should be restricted by the time for the net travel distance of rising blobs above the magnetosphere to become zero. Thus, the time it takes for the RTI layer to grow is:

[TABLE]

We can identify this time with intervals between consecutive flares (the ‘waiting time’). That this time turned out to be of the order of the plasma cooling time is intuitively clear: the next portion of RT-unstable plasma is accumulated during the characteristic time needed for plasma to cool down to enable RTI.

Substituting $t_{\mathrm{rad}}$ from Eq. (3) and the expression for the magnetospheric radius $R_{\mathrm{m}}$ for radiation cooling, Eq. (5), into Eq. (11), we find

[TABLE]

This estimate shows that flare waiting time can be as long as a few thousand seconds. In this model, the comparison of flare waiting times in a particular source enables us to evaluate the dimensionless combination of parameters $A\zeta^{2/9}<1$, which is impossible to obtain from theory.

Fig. 6 displays the flare waiting time against the pre-flare luminosity. Flares from single sources are overlaid in Fig. 13. The straight line indicates the dependence $\Delta T=130[\mathrm{s}]\dot{M}_{16}^{-1}$ (where $\dot{M}_{16}$ is, in this context, the accretion rate measured before each flare). Most flares from single sources follow this relation, with some scatter. Flares from XTE J1739-302 are notable in following this anticorrelation but with a significantly lower normalization. In this source, we can derive \left(\frac{\alpha}{0.03}\right)A\zeta^{2/9}\mu_{30}^{2/3}$$\sim0.03 (the total range covered is 0.01-0.1). A similar situation may be valid for IGR J08408-4503 (although with only two flares from this target it is impossible to draw any conclusion). Also flares from IGR J18483-0311 appear sistematically shifted to a lower value of \left(\frac{\alpha}{0.03}\right)A\zeta^{2/9}\mu_{30}^{2/3}$$\sim0.3. On the other hand, no trend is apparent from IGR J18410-0535, probably because all these flares come from the decaying part of a FRED-like flare, a unique behaviour among SFXTs studied here.

5.2.2 Flare rising time $\delta t$, duration $\Delta t$ and energy $\Delta E$

Flare rising time. Consider a plasma blob rising due to RTI. When the instability starts, the rising blob struggles against the flow velocity towards the magnetosphere, so the net blob velocity is

[TABLE]

where $g=GM_{x}/R_{\mathrm{m}}^{2}$ is the gravity acceleration, the dimensionless factor $a\sim 0.1$, $A\lesssim 1$ is the effective Atwood number, ${\cal R}=2\piup/k$ is the blob curvature radius that we will associate with the instability wavelength $\lambda=2\pi/k$, $k$ is the wavenumber.

In the convective shell above the magnetosphere, plasma is likely to be turbulent (Shakura et al., 2012). In this case, the effective viscosity in the plasma is, according to the Prandtl rule,

[TABLE]

where $v_{t}$ and $l_{t}$ is the characteristic turbulent velocity and scale, respectively. Below we shall scale these quantities with the free-fall velocity and magnetospheric radius, respectively: $v_{t}=\alpha_{v}v_{\mathrm{ff}}$, $l_{t}=\alpha_{l}R_{\mathrm{m}}$, so that the turbulent viscosity can be written in the form $\nu_{t}=(\alpha_{t}/3)v_{\mathrm{ff}}R_{\mathrm{m}}$, where $\alpha_{t}=\alpha_{v}\alpha_{l}\lesssim 1$ is the effective turbulent viscosity parameter (Shakura, 1973).

One of the viscosity effect on RTI is the appearance of the fastest growing mode (Plesset & Whipple, 1974):

[TABLE]

Substituting Eq. (15) into Eq. (13) using Eq. (14), from the condition $v_{b}=0$ we find the time of the most rapidly growing mode:

[TABLE]

We assume that this is the time after which the entire RTI-mixed layer falls onto the NS in the dynamical time producing a flare. Therefore, we identify this time with the flare rising time, $\delta t_{rise}=t_{k_{max}}$ (see Figs. 7 and 14).

Flare duration. During a flare triggered by an external perturbation, the mass $\Delta M$ accumulated in the mixing layer is assumed to fall onto the NS surface over the characteristic dynamical time of the entire shell $t_{\mathrm{ff}}(R_{\mathrm{B}})=R_{\mathrm{B}}^{3/2}/\sqrt{2GM_{x}}$:

[TABLE]

As this time is most sensitive to the wind velocity, the observed dispersion in the flare duration should reflect the stellar wind velocity fluctuations, $\delta\Delta t/\Delta t=-3\delta v_{w}/v_{w}$, and cannot vary more than by a factor of 2-3. Apparently, the long duration of some flares may be a result of ‘gluing’ of several shorter flares into a longer one.

See Fig. 15 for the range of flare duration in individual sources.

Flare energy. The mass accumulated in the RTI mixing layer can be estimated as $\Delta M=\dot{M}_{q}\Delta T$:

[TABLE]

In this approach, the characteristic flare energy due to accretion of the mass $\Delta M$ onto NS, $\Delta E=0.1\Delta Mc^{2}$, turns out to be:

[TABLE]

which is very close to what is observed (see Fig. 5).

The mean mass accretion rate and hence mean X-ray luminosity during the flare is

[TABLE]

Remarkably, in this model the waiting time between flares $\Delta T$ is independent on the mass accretion rate $\dot{M}_{x}$ between flares (to within the possible dependence of the Atwood number $A$ and parameter $\zeta$ on the mass accretion rate). Therefore, these qualitative considerations suggest that the power of flares and mean accretion luminosity in flares in particular sources should be of the same value during individual $XMM$-$Newton$ observations. Moreover, the ratio of the flare energy to the waiting time

[TABLE]

does not depend on unknown RTI parameters $\zeta,\alpha,A$ (which can vary in individual sources) and is proportional only to the mean mass accretion rate between the flares (Figs. 9 and  16, for individual sources). This explains the tight correlation seen in these plots.

Clearly, the actual mass accretion rate during the flares in individual sources is determined by the RTI details (e.g., the fraction of the magnetospheric surface subject to the instability, the effective Atwood number etc.), which cannot be calculated theoretically. However, we note good agreement of the expected flare duration and mean flare energy obtained from this qualitative considerations with observations (Fig. 15). In these estimates, one should also keep in mind the inevitable dispersion, from source to source, in the NS magnetic field (see the dependence on $\mu_{30}$ in above formulas).

If the external mass accretion rate from the stellar wind increases, however, the higher average mass accretion rate can be reached automatically to enable Compton cooling to control plasma entry. We remind that once $\dot{M}_{x}\gtrsim 4\times 10^{16}$  g s*-1* , the settling regime itself ceases altogether, and free-fall supersonic flow occurs until the magnetospheric boundary with the subsequent formation of a shock above the magnetosphere, as was described and studied in more detail in earlier papers (Arons & Lea, 1976; Burnard et al., 1983).

5.3 Could the SFXT flares at their low-luminosity state be

due to propeller mechanism?

The propeller mechanism (Illarionov & Sunyaev, 1975) has been also suggested for the SFXT phenomenon (Grebenev & Sunyaev, 2007; Bozzo et al., 2008) as a mechanism for gating accretion onto a rapidly rotating magnetized NS. It is feasible for disc accretion and is likely observed in luminous transient X-ray pulsars (Tsygankov et al., 2016). In the case of quasi-spherical accretion, the propeller mechanism can be involved to explain major observational features of enigmatic $\gamma$ Cas stars (Postnov et al., 2017). For low-states of SFXTs, the propeller mechanism, which requires centrifugal barrier for accretion by the condition that the Alfvén radius $R_{\mathrm{m}}\sim\dot{M}^{-2/7}$ be larger than the corotation radius, $R_{\mathrm{c}}=(GM_{\mathrm{x}}P_{x}^{2}/4\piup^{2})^{1/3}$, would need either a fast NS rotation or a large NS magnetic field: $P_{x}\leq P_{cr}\simeq 9[\mathrm{s}]\mu_{30}^{6/7}\dot{M}_{16}^{-3/7}$. Clearly, with the quiescence X-ray luminosity $L_{x}\sim 10^{34}$ erg s*-1*, the fastest SFXT from Table 2, IGR J18483-0311, with a pulse period of $P_{x}\sim 21$ s could be at the propeller stage. If so, the low X-ray luminosity can be due to the leakage of matter across the magnetospheric surface (for example, close to the rotational axis).

In the case of quasi-spherical turbulent shell above magnetosphere, the propeller regime should correspond to the so-called ‘strong coupling’ between the magnetic field and surrounding matter, when the toroidal field component is approximately equal to the poloidal one, $B_{t}\sim B_{p}$ (Shakura et al., 2012, 2018). In this regime, the NS spins down at a rate

[TABLE]

($I\approx 10^{45}$ g cm2 is the NS moment of inertia, $K_{2}\simeq 7.6$ is the numerical coefficient accounting for the structure of quasi-spherical NS magnetosphere, Arons & Lea 1976), corresponding to a spin-down time of less than $10^{5}$ yrs. This short time suggests that a fast spinning magnetized NS rapidly approaches the critical period to become accretor, $R_{\mathrm{m}}\sim R_{\mathrm{c}}$, and this fact was stressed already in the original paper by Illarionov & Sunyaev (1975).

There is a difference between the propeller stage for disc accretion and quasi-spherical accretion. In the former case, the disc is produced by accreting matter and the material is expelled by the rotating NS magnetosphere along open magnetic field lines (Lii et al., 2014). In the disc case, sporadic accretion episodes during the transition to the accretion stage were found once $R_{\mathrm{m}}\simeq R_{\mathrm{c}}$ (e.g, for pre-outburst flares in A 0535+26, Postnov et al. 2008). However, the numerical simulations by Lii et al. (2014) were carried out for small magnetospheres, and in the case of large magnetospheres the situation remains unclear.

In the quasi-spherical case onto large magnetospheres (low accretion rates or high magnetic fields), the matter acquires the (super-Keplerian) specific angular momentum of the magnetosphere $\sim\omega_{x}R_{\mathrm{m}}^{2}$. If the cooling time of this matter (e.g., if it is expelled in the form of dense blobs) is short compared to the dynamic (convection) time in the shell, an equatorial ring with radius $R_{\mathrm{p}}\simeq R_{\mathrm{m}}(R_{\mathrm{m}}/R_{\mathrm{c}})^{3}$ and some thickness $h\ll R_{\mathrm{p}}$ (Shakura et al., 2012) is likely to form. This ring spreads over in the viscous diffusion time scale, $t_{d}\sim t_{\mathrm{K}}(R_{\mathrm{p}})(R_{\mathrm{p}}/h)^{2}$ ($t_{\mathrm{K}}$ is the Keplerian time), and may end up with an accretion episode once the inner disc radius overcomes the centrifugal barrier. Therefore, it is possible to characterize the time between flares by the ring diffusion time scale, over which the disc replenishes mass by freshly propelled matter.

In so far as accretion through such a disc is centrifugally prohibited, its structure should be described by equations of ‘dead’ discs (Syunyaev & Shakura, 1977), with the characteristic thickness $h\sim t_{\mathrm{K}}(R_{\mathrm{p}})^{6/7}\Sigma_{0}^{3/14}$, where $\Sigma_{0}$ is the surface density at its outer edge. Assume that the mass stored in this disc in a time interval $\Delta T$ be $M_{\mathrm{d}}\sim\dot{M}_{\mathrm{p}}\times\Delta T$, where $\dot{M}_{\mathrm{p}}$ is the fraction of the mass accretion rate propelled from the magnetosphere, $\dot{M}_{\mathrm{p}}=\dot{M}-\dot{M}_{\mathrm{x}}$ ($\dot{M}_{\mathrm{x}}$ – the fraction of the mass accretion rate that reaches the NS surface and produces the inter-flare X-ray luminosity). In the simplest case, $\dot{M}_{\mathrm{x}}$ is geometrically determined by the centrifugally free fraction of the magnetosphere surface, $\dot{M}_{\mathrm{x}}=\dot{M}(1-\sqrt{1-(R_{\mathrm{c}}/R_{\mathrm{m}})^{2}}$.

Next, we use the relation $M_{\mathrm{d}}\sim R_{\mathrm{p}}^{2}\Sigma_{0}$ and note that the corotation radius $R_{\mathrm{c}}$ remains pretty much constant over short time intervals. Then by identifying the time between flares $\Delta T$ with the viscous time of such a dead disc, we arrive at the relation $\Delta T\sim R_{\mathrm{p}}^{5/4}\dot{M}_{\mathrm{p}}^{-3/10}$. As the ring radius $R_{\mathrm{p}}$ scales with mass accretion rate as $R_{\mathrm{m}}^{4}\sim\dot{M}^{-8/9}$ and the propelled mass rate $\dot{M}_{\mathrm{p}}$ scales as $\dot{M}$, we arrive at $\Delta T\sim\dot{M}^{-127/90}\approx\dot{M}^{-1.4}$. The mass accumulated in the dead disc between the flares turns out to be inversely dependent on the mass accretion rate between the flares, $M_{\mathrm{d}}\sim\dot{M}^{-0.4}$.

In a quite different setup, the magnetospheric instability could be related to perturbations in the magnetosphere. These perturbations propagates with the Alfvén velocity, $v_{A}=B_{p}/\sqrt{4\piup\rho}\sim v_{ff}(R_{\mathrm{m}})$. Therefore, the characteristic time between accretion episodes due to these perturbations would be $\Delta T\sim t_{A}\sim R_{\mathrm{m}}/v_{A}\propto R_{\mathrm{m}}^{3/2}$. For any regime (disc or quasi-spherical), $R_{\mathrm{m}}$ scales with $\dot{M}$ not stronger than $R_{\mathrm{m}}\sim\dot{M}^{-2/7}$, and thus $\Delta T\sim\dot{M}^{-3/7}$ (disc) or $\sim\dot{M}^{-1/3}$ (quasi-spherical, radiation cooling). In the last plot of Fig. 13, which outlines the case of the the fastest pulsar in a SFXT known to date (IGR J18483–0311), we show the dependences $\Delta T\propto\dot{M}^{-1/3}$ and $\Delta T\propto\dot{M}^{-1.4}$.

Thus we conclude that the propeller model for SFXT flares can also provide the qualitative inverse dependence of the flare waiting time on the pre-flare luminosity. A more detailed analysis of the propeller mechanism at low accretion rates onto large NS magnetospheres definitely deserves further investigation, which is far beyond the scope of the present paper.

6 Summary

To summarize, here we propose the following model explaining the flaring behaviour of SFXTs at their low luminosity state which is based on the statistical analysis of properties of the $XMM$-$Newton$ B.b. light curves:

• •

At the quiescent states of SFXTs with low X-ray luminosity $\sim 10^{33}$  erg s*-1* , the RTI is ineffective to enable rapid plasma penetration into the NS magnetosphere. Instead, either a steady settling accretion regime controlled by the radiative plasma cooling occurs and the mass accretion rate onto the NS is reduced by factor $f(u)\ll 1$ compared to maximum available Bondi-Hoyle-Littleton value $\dot{M}_{B}\simeq\rho_{w}R_{\mathrm{B}}^{2}/v_{w}^{3}$, $\dot{M}_{x}=f(u)\dot{M}_{B}\approx\dot{M}_{B}(t_{\mathrm{ff}}/t_{\mathrm{rad}})^{1/3}$ (Shakura et al., 2012, 2018), or plasma enters the magnetosphere via ineffective processes (e.g., diffusion or magnetospheric cusp dripping, see Elsner & Lamb 1984). In the last case, a low X-ray luminosity of $10^{32}-10^{33}$  erg s*-1* can be sustained by the thermal X-ray emission of the hot magnetospheric shell (see Postnov et al. (2017) for a more detailed discussion and possible applications to the $\gamma$ Cas phenomenon).

• •

A series of flares can be triggered by an external fluctuation of the stellar wind properties (density $\rho_{w}$ and/or velocity $v_{w}$). Individual $XMM$-$Newton$ X-ray light curves of SFXTs (see Figs. 11 and  12) suggest that in most cases a series of flares is initiated by a small- or moderate-amplitude flare, with subsequent development of more powerful flares and gradual decrease in flare amplitudes. Typically, such flare series last for about $\sim 1000$ s, a few dynamical time scales of the problem.

• •

The analysis of individual flares shows (see, e.g., Table 1, sixth column) that the mean X-ray luminosity during the flare very rarely exceeds $\sim 10^{36}$ erg s*-1*. This can explain why these flares cannot switch-on the development of Compton-cooling controlled RTI and thus does not turn the source into a steady-state wind accreting state like Vela X-1. Instead, a series of flare terminates when all matter stored in the magnetopsheric shell is exhausted by the small RTI-flares. This is to be contrasted with bright SFXT X-ray flares during which the entire magnetospheric shell can accrete onto the NS because of the magnetosphere breakage due to, for example, reconnection of the magnetic field carried out by stellar wind plasma (Shakura et al., 2014).

• •

On average, the flare energy in individual sources should be proportional to the fraction of the shell subject to RTI, $\delta M/M_{sh}\sim Z/R_{\mathrm{m}}$, where $Z$ is the size of the RTI layer. In the convective/turbulent shell, there is a turbulent viscosity that singles out a specific wavelength growing most rapidly, $Z\sim 4\pi[2\alpha^{2}/(45A)]^{1/3}$, with $\alpha<1$ being the turbulent viscosity coefficient (a la Shakura-Sunyaev in discs), $A\lesssim 1$ the Atwood number. The shell mass is (Shakura et al 2014) $M_{sh}\sim\dot{M}_{x}t_{\mathrm{ff}}(R_{\mathrm{B}})\sim\dot{M}_{x}v_{w}^{-3}$. Therefore, the mean energy of flares in individual source $<\Delta E_{f}>\sim\delta M\sim<L_{x,q}>v_{w}^{-3}$, i.e. on average linearly grows with the mean pre-flare X-ray luminosity $<L_{x,q}>$ (Fig. 10).

• •

In each individual source, $E_{f}\sim\dot{M}_{x}\times\Delta T$, where $\Delta T$ is the ‘waiting time’ between individual flares, which is $\Delta T\sim t_{\mathrm{rad}}\sim 1/\dot{M}_{x}$. Thus, in each source, $E_{f}$ must be independent on variations of $\dot{M}_{x}$ between the flares.

• •

The spread of the mean X-ray luminosities $<L_{x,q}>$ between flares in individual sources is determined by the fractional change in the mass accretion rate onto the NS due to variations in the Bondi mass accretion rate captured from the stellar wind. For example, in the settling accretion theory with radiative plasma cooling $\dot{M}_{x}=f(u)\dot{M}_{B}\sim\dot{M}_{B}\dot{M}_{x}^{2/9}$, hence $\dot{M}_{x}\sim\dot{M}_{B}^{9/7}$. Therefore, the fractional change in $L_{x,q}$ in the source between flares is $\delta L_{x}/L_{x}=9/7\delta\dot{M}_{B}/\dot{M}_{B}=(9/7)\delta\rho_{w}/\rho_{w}-(27/7)\delta v_{w}/v_{w}$, where $\delta\rho_{w}/\rho_{w}$ and $\delta v_{w}/v_{w}$ is the stellar wind density and velocity fluctuations, respectively. These variations in $L_{x}$ up to one order of magnitude can be produced during the active RTI stage.

• •

The rising time of a flare can correspond to the fastest growing RTI mode in the turbulent shell, $\delta t_{r}\sim 30\mbox{s}\,\dot{M}_{16}^{-2/3}$. The inverse dependence of the flare rising time on the X-ray luminosity between flares can be seen for some individual sources (Fig. 14).

• •

For the fastest pulsar in IGR J18483–0311, the centrifugal barrier at the magnetospheric boundary may lead to the formation of an equatorial dead cold disc which could trigger a flaring activity of the source once the centrifugal barrier at its inner edge is overcome. The waiting time between the flares in this case can be characterized by the viscous time scale of the disc evolution and is also inversely proportional to the pre-flare X-ray luminosity as $\Delta T\sim\dot{M}_{x}^{-1.4}$.

We conclude that SFXT flares observed during low-luminosity states could be qualitatively compatible with the development of the Rayleigh-Taylor instability in plasma accreted from the stellar wind of the companion which tries to enter the NS magnetosphere.

However, the full development of the RTI fails because the radiative plasma cooling during the flares turns out to be insufficient for Compton cooling to enable steady-state magnetospheric plasma penetration, as in the case of persistent wind accreting X-ray pulsars like Vela X-1.

Thus, SFXT flares offer unique possibility to probe complicated processes of plasma entering into magnetopsheres of magnetic NSs through interchange instabilities under natural conditions.

Acknowledgments

This work is based on observations obtained with $XMM$-$Newton$, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. We have made use of data produced by the EXTraS project, funded by the European Union’s Seventh Framework Programme under grant agreement no 607452. We acknowledge financial support from the Italian Space Agency (ASI) through the ASI-INAF agreement 2017-14-H.0. The EXTraS project acknowledges the usage of computing facilities at INAF - Astronomical Observatory of Catania. The EXTraS project acknowledges the CINECA award under the ISCRA initiative, for the availability of high performance computing resources and support. The work of KAP is supported by RFBR grant 18-502-12025, and by the grant of Leading Scientific Schools of Moscow University ‘Physics of stars, relativistic objects and galaxies’. The authors thank the anonymous referee for the useful and constructive comments, and suggestions for alternative explanations for the SFXT flaring activity.

Appendix A Flare properties

In Table 1 we list the properties of the X-ray flares selected from the sample of SFXTs investigated in this work. Temporal quantities have been rounded to the significant digits. The ($-$) symbol means that the value of the parameters could not be determined, according to the definitions assumed in Sect. 4.1. In Figs. 11 and  12 we report the EPIC source light curves, segmented in B.b. In Figs. 13, 14, 15 and 16 we highlight the behaviour of flares from single sources.