A bright X-ray flaring event from the intermediate supergiant fast X-ray transient IGR J00370+6122 observed with NuSTAR

Introduction

High-mass X-ray binaries (HMXBs) consist of an OB-type donor star and a compact object, either a neutron star (NS) or, more rarely, a black hole, that accretes matter from its companion and produces X-ray emission (for a review, see e.g. Kretschmer et al. 2019). Within this class, supergiant X-ray binaries (SGXBs) host a compact object that orbits an OB supergiant donor and accretes from its powerful, radiatively driven stellar wind. Classical persistent SGXBs (e.g. Vela X–1 and 4U 1700–377) are characterised by continuous X-ray emission with moderate variability, typically about one order of magnitude, which is generally attributed to intrinsic clumpiness and inhomogeneities in the wind from the donor star (e.g. Fürst et al. [16]).

The discovery rate of SGXBs increased substantially following the 2002 launch of the International Gamma-Ray Astrophysics Laboratory (INTEGRAL; Winkler et al. [56]). Its hard X-ray monitoring capabilities revealed a previously unrecognised subclass of highly variable sources, subsequently classified as supergiant fast X-ray transients (SFXTs; Negueruela et al. [31]; Sguera et al. [43]; Smith et al. [47]). Unlike classical SGXBs, SFXTs spend most of their time in a faint, quiescent state ($L_{\mathrm{X}} \sim10^{32}-10^{33}$ erg s$^{-1}$), interrupted by rare, bright outbursts lasting a few tens of s, during which the luminosity typically reaches $L_{\mathrm{X}} \sim10^{36}-10^{37}$ erg s$^{-1}$ (e.g. Romano [36]; Sidoli [46]); a few exceptional cases reach up to $\sim10^{38}$ erg s$^{-1}$ (e.g. Romano et al. [37]). These outbursts are composed of multiple short flares with durations of $\sim1000$ s. On short timescales ($<10^{4}$ s), SFXTs display dynamic ranges in luminosity of $\sim10^{3}-10^{5}$, much larger than those of classical SGXBs, which makes them one of the most extreme classes of variable HMXBs known. To date, several dozen SFXTs and SFXT candidates have been identified, although orbital and spin periods have been measured for only a subset of these sources (e.g. Romano [36]; Sidoli [46]).

The physical mechanism responsible for the extreme variability of SFXTs remains debated, and more than one mechanism likely contributes across the class. Proposed explanations include (i) clumpiness of the supergiant wind, whereby accretion of individual dense clumps produces sudden X-ray flares (e.g. in ’t Zand 2005); (ii) centrifugal and magnetic gating mechanisms, in which the NS spin and its strong magnetic field can temporarily inhibit accretion (e.g. Grebenev & Sunyaev [18]; Bozzo et al. [10]); and (iii) quasi-spherical settling accretion, in which a hot, quasi-static shell forms above the magnetosphere and regulates the plasma entry rate through Rayleigh-Taylor instabilities. In this scenario, moderate flares arise when local fluctuations in the wind density or velocity shift the cooling regime near the magnetosphere, while brighter flares occur when the NS captures a magnetised wind structure. This triggers magnetic reconnection, which induces the unstable collapse of the shell onto the NS surface (Shakura et al. [45]).

Between the two extremes of classical SGXBs and highly variable SFXTs lies a small group of systems commonly referred to as ‘intermediate SFXTs’. These sources show qualitatively similar flaring activity to that of SFXTs, but with a considerably smaller dynamic range, typically of order $\sim10^{2}$ ([54]).

IGR J00370+6122 is one such intermediate SFXT, which was discovered in December 2003 with INTEGRAL in the $20$–$60\,\mathrm{keV}$ ([12]) energy range. The All-Sky Monitor (ASM) on board the Rossi X-ray Timing Explorer (RXTE) ($1.5$–$12\,\mathrm{keV}$) revealed a periodic modulation of the X-ray flux with a period of $15.665 \pm0.006$ d, which was interpreted as the orbital period of the binary system ([12]). Uchida et al. [52] later refined the orbital period to $15.6649 \pm0.0014$ d (90% c.l.) using 15 yr of Swift/Burst Alert Telescope (BAT) monitoring data. Radial velocity measurements have established the orbital parameters of the system. In particular, Grunhut et al. [19] and González-Galán et al. [17] derived eccentricities of $e = 0.48^{+0.02}_{-0.03}$ and $e = 0.56 \pm0.07$, respectively. Recent measurements based on Gaia data place the system at a distance of approximately $3.4 \pm0.2$ kpc ([8]).

The optical counterpart, BD +60° 73, has been the subject of several spectroscopic investigations. González-Galán et al. [17] classified it as a BN0.7Ib based on high-resolution spectra and model atmosphere analysis. Hainich et al. [20] confirmed the main parameters of the optical donor through a detailed UV and optical spectral analysis using the Potsdam Wolf-Rayet atmosphere models. González-Galán et al. [17] found that BD +60° 73 showed nitrogen enrichment, indicating an advanced evolutionary stage. They also classified it as an unusual low-luminosity supergiant.

The nature of the compact object in IGR J00370+6122 has been investigated through timing analyses. in’t Zand et al. (2007) detected a possible $346 \pm6$ s periodicity during a bright X-ray flare with the RXTE/Proportional Counter Array (PCA), which they interpreted as the spin period of an accretion-powered pulsar. However, because the flare was short, the signal repeated only seven times, making this detection uncertain. Subsequently, Uchida et al. [52] analysed XMM-Newton data, finding evidence for a $674 \pm15$ s periodicity; they propose that the $\sim346$ s signal is the second harmonic of this fundamental period. Nevertheless, this revised periodicity also remains unconfirmed because the signal is transient and appears predominantly during bright flaring states, with only eleven cycles observed.

In the X-ray band, IGR J00370+6122 shows a luminosity ranging from $\sim10^{32}\,\mathrm{erg\,s^{-1}}$ to $\sim10^{36}\,\mathrm{erg\,s^{-1}}$ ([52]). The X-ray spectral properties have been studied across a broad energy range. in’t Zand et al. (2007) modelled the $3$–$60\,\mathrm{keV}$ flare spectrum with an absorbed power law ($\Gamma\approx2$, $N_{\mathrm{H}} \approx9 \times10^{22}\,\mathrm{cm^{-2}}$). Uchida et al. [52] reported a harder-when-brighter trend in the $1$–$10\,\mathrm{keV}$ band. The combined RXTE/PCA and INTEGRAL Soft Gamma-Ray Imager (ISGRI) $10$–$80\,\mathrm{keV}$ spectrum is well described by a single power law ($\Gamma\approx2$) with no clear high-energy cutoff. The emission is strongly modulated by the orbital phase: the peak occurs after periastron passage, with a delay of about $0.05$ in orbital phase, and the luminosity drops suddenly from $\sim4 \times10^{33}\,\mathrm{erg\,s^{-1}}$ to $\sim10^{32}\,\mathrm{erg\,s^{-1}}$ at an orbital phase of $0.3$ ([52]).

González-Galán et al. [17] suggested that IGR J00370+6122 is an intermediate SFXT because its X-ray flares resemble those of SFXTs, but its flux varies over a considerably smaller dynamic range ($\sim10^{2}$) than that of canonical SFXTs ($\gtrsim10^{3}$; e.g. [36, 46]). Intermediate SFXTs are generally interpreted as bridging the gap between highly variable, flare-dominated SFXTs and persistently bright, wind-fed classical SGXBs. Fewer than ten intermediate SFXTs are known; two other examples are IGR J18483−0311 and IGR J17354−3255 ([35, 44]).

In this paper, we present an analysis of a long NuSTAR observation of IGR J00370+6122 that was carried out in 2026 around the expected X-ray peak of the orbital modulation. The observation lasted about 1.8 d, with a net exposure time of about 93 ks. We report all uncertainties at the 1$\sigma$ confidence level unless otherwise stated.

Data analysis

NuSTAR

The Nuclear Spectroscopic Telescope Array (NuSTAR) satellite is equipped with two identical co-aligned telescopes, each with a focal plane module (FPMA and FPMB). Both modules operate in the $3$–$79\,\mathrm{keV}$ energy band ([21]). The NuSTAR satellite observed IGR J00370+6122 on 1 March 2026 for a net exposure time of approximately 93 ks (ObsID 31101014002). We reduced the data using NuSTAR DAs v2.1.5, which is part of HEASOFT v6.35.1, and the calibration files distributed with CALDB v20260309.

We extracted source events from a circular region centred on the source. We optimised the extraction radius to maximise the signal-to-noise ratio in the hard part of the X-ray spectrum, where the source flux is lower and the background contribution is more significant. To achieve this, we used the nustar-gen-utils package (Madsen et al. 2023)². The analysis yielded an extraction radius of $50''$ centred on the source coordinates. We extracted background events from a circular region with a radius of $80''$ centred at RA = $00^\mathrm{h}36^\mathrm{m}46.57^\mathrm{s}$ and Dec = $+61^\circ24'45''.6$ (J2000), located on the same detector in a source-free area. We also tested different extraction radii for the source and background and verified that the spectral results are consistent. We corrected the event times from the satellite frame to the Solar System barycentre using the barycorr tool and the JPL EPHEM.440 ephemeris.

The NuSTAR data are not continuous. Regular gaps, typically lasting about 21 minutes, are caused by Earth occultations of the source during the satellite’s low-Earth orbit, whereas shorter, irregular gaps are introduced by standard GTI screening, which removes passages through the South Atlantic Anomaly and intervals affected by other instrumental constraints.

We searched for the spin period in the $3$–$70\,\mathrm{keV}$ events of modules A and B (combined and separated) using a $Z^2$ test (e.g. [11]) with up to three harmonics. We restricted the search to periods between 1 s and $4 \times10^{4}$ s. We detected no statistically significant pulsations, including in the previously reported periods of $\sim674$ s and $\sim346$ s. In addition, given previous reports that the tentative spin-periods were transient and lasted only several cycles ((in’t Zand et al. 2007; Uchida et al. 2021)), we employed the weighted wavelet transform (WWT; Foster [15]) to track transient modulations that may not remain stationary throughout the entire observation. For the WWT, we used barycentre-corrected FPMA and FPMB data in the $3$–$20\,\mathrm{keV}$ range (PI channels 35–460) with 8 s binning. We obtained the transforms for both detectors independently and combined them using their geometric mean; we generated the WWT maps in the 160–800 s range (Figure 1). The resulting WWT map indicates that the strongest variabilities are concentrated mainly during the bright part of the observation ($\sim70$–$85 ks after the start). Transiently enhanced WWT powers are also present around 674 s and 346 s; however, they appear to form part of the formation of several-hundred-second-long, evolving structures rather than a horizontal ridge, as would be expected from coherent pulsations. Thus, the WWT map does not reveal a stable coherent periodicity at $\sim346$ and $\sim674$ s. These broad structures can more naturally be attributed to sporadic variabilities due to flaring episodes.

To better highlight the variability properties of the source, we re-binned the background subtracted light curve using the optimal segmentation technique based on the Bayesian-block representation described by Scargle et al. [42].³ We applied this algorithm to the event lists to partition the data into an optimal number of blocks, where each block represents a statistically significant change in the count rate. We adopted a false-positive rate probability of $p_0 = 0.01$ for the Bayesian-blocks analysis to remain sensitive to the rapid rate fluctuations exhibited by the flares. Because the algorithm accounts for the exposure of each good time interval, it does not merge blocks separated by a gap, even when their count rates are consistent within the uncertainties. For this reason, the light curve shown in Figure 3 contains consecutive, but not contiguous, blocks with comparable rates.

We extracted the FPMA and FPMB spectra from the same regions used for the timing analysis. We rebinned the NuSTAR spectra to a minimum number of 30 counts per bin so that $\chi^2$ could be used as the fit statistic. We performed the spectral analysis using XSPEC v12.15.0 (Arnaud [2]). We included a free multiplicative constant in the spectral fitting to account for calibration uncertainties between the FPMA and FPMB modules. We modelled the column density ($N_{\rm H}$) with the tbabs model in XSPEC and we assumed the interstellar-medium abundances (wilm in XSPEC; [55]) and vern cross sections (Verner et al. [53]).

Swift/BAT

The orbital period of IGR J00370+6122 is known to be $\sim15.66$ d, although slightly different values have been reported in the literature. The most recent determination, based on approximately 15 yr of Swift/BAT observations, is 15.6649 $\pm$ 0.0014 d at a 90% confidence level (Uchida et al. [52]). We refined the orbital period using the full Swift/BAT monitoring available to date, covering approximately 21 yr. We retrieved the 15–50 keV orbital light curve of IGR J00370+6122 from the Swift/BAT Transient Monitor⁴ (Krimm et al. [27]). The light curve is averaged over the Swift orbital period of $\sim90$ minutes and covers the interval between MJD 53416 and 61068. We excluded all measurements with quality flags 1, 2, and 3. We searched for periodic signals using the Lomb–Scargle periodogram (Lomb [28]; Scargle [41]) implemented in the astropy Python package.⁵ We adopted the fast $\chi^2$ algorithm of Palmer [33] with three Fourier terms. We searched over the 2–1000 d period range, oversampling the independent frequency spacing by a factor of 20 to improve the determination of the periodogram peak. We estimated the uncertainty on the orbital period through 1000 Monte Carlo simulations. We generated synthetic light curves by randomly sampling each data point according to its measured statistical uncertainty, assuming Gaussian errors. We adopted the standard deviation of the recovered orbital periods as the uncertainty.

Results

Orbital period

We determined the orbital period of IGR J00370+6122 by applying the Lomb–Scargle periodogram to the full Swift/BAT light curve. The resulting value is $P_{\rm orb} = 15.6654 \pm0.0008$ d, where the uncertainty is the $1\sigma$ standard deviation. This is currently the most accurate orbital-period measurement for this source and improves on the previous measurement by Uchida et al. [52]. We folded the Swift/BAT light curve using this orbital period and adopted a reference epoch $T_0 = 55084.02$ MJD (phase 0), which corresponds to the periastron passage calculated by González-Galán et al. [17]. Figure 2 shows the folded light curve. The horizontal red lines in the same figure mark the orbital-phase interval covered by the NuSTAR observation.

X-ray variability

The top panel of Figure 3 shows the $3$–$79$ keV light curve, obtained with the optimal segmentation technique based on the Bayesian-block representation (see Sect. 2). This light curve displays a strong flaring activity. During the first approximately $40$ ks of the observation, IGR J00370+6122 shows mild flaring activity, with a maximum flux ratio of $F_{\max}/F_{\min}$ of approximately $40$. After this period, the flux decreases significantly, dropping by a factor of approximately ten for roughly $20$ ks. This decline is followed by a brighter phase of flaring activity, as shown in the bottom panel of Figure 3. Hereafter, we refer to this episode as the outburst.

Spectral variability

To investigate the spectral variability of IGR J00370+6122 during the NuSTAR observation we constructed hardness-ratio (HR) light curves by extracting source counts in the soft ($3$–$7$ keV) and hard ($7$–$79$ keV) energy bands. The top panel of Figure 4 shows the total ($3$–$79$ keV), soft, and hard light curves together with the evolution of the HR as a function of time, while the bottom panel displays the HR as a function of the total X-ray flux. The HR analysis reveals mild spectral variability throughout the observation, with the source generally becoming harder at higher flux levels.

To investigate this behaviour further, we performed HR-resolved spectroscopy by extracting two separate spectra based on the hardness ratio: one comprising data with HR $\leq0.7$ (hereafter, the soft spectrum) and the other with HR $>0.7$ (the hard spectrum). We chose the threshold value of HR $=0.7$ to divide the data into two spectra with comparable statistics.

For the continuum, we tested several models commonly used to describe the spectra of accreting compact objects in HMXBs, including a cutoff power law and a power law with a high-energy exponential cutoff. We find that the best fit, based on the $\chi^2$ statistic and the residuals, is obtained with a double-Comptonisation model (comptt + comptt in XSPEC; Titarchuk [49]). This model has been used successfully to describe the broadband spectra of several X-ray pulsars observed at low luminosities ($L_{\mathrm{X}} \approx10^{33}$–$10^{35}$ erg s$^{-1}$), including X Per (Doroshenko et al. [13]), GX 304−1 (Tsygankov et al. [51]), A 0535+262 (Tsygankov et al. [50]), 2RXP J130159.6−635806 (Salganik et al. [38]), MAXI J0655−013 (Malacaria et al. [29]), RX J0440.9+4431 (Salganik et al. [39]), and 1A 1118−61 (Salganik et al. [40]). Interestingly, in the latter two sources, the double hump is seen up to X-ray luminosities of $\sim10^{36}$ erg s$^{-1}$. In all the low-luminosity accreting pulsars mentioned above, the double-comptt model is interpreted as arising from two distinct emission components: a soft component originating from thermal Comptonisation in the accretion-heated polar cap region, and a hard component associated with either cyclotron emission reprocessed by magnetic Compton scattering in the overheated upper layers of the NS atmosphere or thermal Comptonisation in a hotter plasma (Tsygankov et al. [51], and references therein).

Figure 5 shows the best-fitting spectra for the soft and hard HR states, and Table 1 shows the corresponding best-fit parameters.

Parameter	HR $>0.7$	HR $\leq0.7$
$N_{\mathrm{H}}$ ($10^{22}$ cm$^{-2}$)	$4.9^{+0.9}_{-1.0}$	$4.8^{+1.2}_{-0.9}$
$kT_0$ (keV)	$0.29^{+0.08}_{-0.08}$	$0.31^{+0.09}_{-0.10}$
$kT_{\mathrm{e},1}$ (keV)	$3.0^{+0.3}_{-0.3}$	$1.92^{+0.20}_{-0.10}$
$\tau_1$	$11.8^{+1.4}_{-0.8}$	$15.4^{+1.1}_{-2.0}$
$\mathrm{norm}_1$	$9.8^{+3.5}_{-1.8} \times10^{-3}$	$5.8^{+1.7}_{-1.4} \times10^{-3}$
$kT_{\mathrm{e},2}$ (keV)	$44^{+7}_{-10}$	$30.5^{+2.3}_{-2.4}$
$\tau_2$	$5.5^{+0.8}_{-0.7}$	$5.0^{+0.7}_{-0.6}$
$\mathrm{norm}_2$	$1.3^{+0.5}_{-0.2} \times10^{-4}$	$5.2^{+1.6}_{-1.2} \times10^{-5}$
$\mathrm{const}_{AB}$	$1.045^{+0.010}_{-0.013}$	$1.037^{+0.012}_{-0.014}$
$F_{\mathrm{X}}$ (erg cm$^{-2}$ s$^{-1}$)	$(1.022^{+0.024}_{-0.010}) \times10^{-10}$	$(2.69^{+0.07}_{-0.06}) \times10^{-11}$
Total $\chi^2$/d.o.f.	$630.19/588$	$456.58/510$

Notes. The quantity $F_{\mathrm{X}}$ is the $3$–$70$ keV absorbed flux. Both comptt components share the same $kT_0$ value.

eters. We modelled the absorption with a single tbabs component, so the column densities reported in Table 1 include both Galactic and local contributions. Because the Galactic column density in the direction of the source is $N_{\mathrm{HGal}} \approx4.8 \times10^{21}$ cm$^{-2}$ (HI4PI Collaboration et al. [22]), approximately an order of magnitude lower than the measured values, the observed absorption is mainly intrinsic to the source. Intrinsic column densities of this order are commonly observed in HMXBs, including SFXTs, where they are usually ascribed to the dense stellar wind of the supergiant donor in which the compact object is embedded (e.g. Kretschmar et al. [26]).

An alternative continuum model consisting of an absorbed blackbody plus a power law also provides a reasonable description of the data, although with worse $\chi^2$ (hard spectrum: $\chi^2$/d.o.f. $=647.44/591$; soft spectrum: $\chi^2$/d.o.f. $=465.43/514$) and less satisfactory residuals. For this model, we obtained a blackbody temperature of $kT_{\mathrm{bb}} \approx1.4$–$1.8$ keV, a blackbody radius of $R_{\mathrm{bb}} \sim130$–$160$ m, an absorption column of $N_{\mathrm{H}} \approx2$, and a column density varying from $N_{\mathrm{H}} \sim2.3 \times10^{22}$ cm$^{-2}$ for the hard spectrum to essentially zero for the soft spectrum. Given the better fit quality of the double-comptt model, we adopted it as our preferred description of the continuum.

We used the spectral results presented above to estimate the minimum and maximum luminosity reached by IGR J00370+6122 during the NuSTAR observation. We converted the lowest and highest count rates of the light curve shown in Figure 3 into fluxes, using the unabsorbed fluxes of the two best-fitting models reported in Table 1. Because the HR correlates with flux (bottom panel of Figure 4), we adopt the HR $\leq0.7$ model for the minimum count rate and the HR $>0.7$ model for the maximum count rate. For an assumed distance of $3.4$ kpc (Bailer-Jones et al. [8]), we obtain $3$–$70$ keV luminosities of $\sim10^{34}$ erg s$^{-1}$ and $\sim1.8 \times10^{36}$ erg s$^{-1}$, corresponding to a dynamic range of approximately $200$. When extrapolated to the $0.1$–$100$ keV band, the luminosities derived from the same models are $\sim1.6 \times10^{34}$ erg s$^{-1}$ and $\sim3 \times10^{36}$ erg s$^{-1}$.

Discussion

Throughout the NuSTAR observation, IGR J00370+6122 remained below $L_{\mathrm{x}} \simeq4 \times10^{36}$ erg s$^{-1}$, the limit below which the shell above the magnetosphere cools inefficiently and quasi-spherical settling accretion is expected to operate ([45, 34]). The luminosity reached during the brightest flares approaches this limit, but only for a total of about $\sim40$ s over the entire exposure, and it remains below the limit even then. Throughout the observation, IGR J00370+6122 therefore remained in the luminosity regime where a quasi-static shell can form above the magnetosphere and regulate the accretion onto the NS. Using the orbital and stellar-wind parameters of IGR J00370+6122 from González-Galán et al. [17] and Hainich et al. [20], we therefore checked whether the observed flaring properties are consistent with episodes of enhanced plasma entry, triggered by instabilities within the quasi-static shell that surrounds the NS magnetosphere. We also note that the low-state luminosity of IGR J00370+6122 (a few $10^{34}$ erg s$^{-1}$) is of the order expected for the non-flaring state of an SFXT in the model by Shakura et al. [45], and that the dynamic range of about two orders of magnitude spanned by IGR J00370+6122 during the NuSTAR observation is consistent with the flare amplitudes that can be produced when the magnetosphere tospheric plasma entry rate switches from the settling value to free fall.

In the quasi-spherical settling accretion model, a flare develops on the free-fall time-scale from the outer boundary of the shell, $t_{\mathrm{ff}} = R_{\mathrm{acc}}^{3/2}/\sqrt{G M_{\mathrm{NS}}}$, with

where $v_{\mathrm{w}}$ is the radial wind velocity at the orbital separation of the NS, and $v_{\mathrm{r}}$ and $v_{\varphi}$ are the radial and azimuthal components of the orbital velocity [32]. For $v_{\mathrm{w}}$ we adopted the double-$\beta$ velocity law of Hainich et al. [20]. The orbital terms cannot be neglected, since the observation covers periastron passage, where $v_{\varphi}$ is comparable to the local wind velocity (see also Hainich et al. [20]). Two interpretations are then possible.

1. The outburst as a single shell collapse. The interpretation of the entire outburst as resulting from shell accretion requires $t_{\mathrm{ff}} \approx5 - 7 \times10^{4}\ \mathrm{s}$ and hence a low relative velocity. We used a terminal wind velocity of $v_{\infty} = 900\ \mathrm{km\ s^{-1}}$, corresponding to the $1\sigma$ lower bound from Hainich et al. [20], together with the orbital solution of González-Galán et al. [17] to obtain $R_{\mathrm{acc}} \approx3.1 \times10^{11}\ \mathrm{cm}$ and $t_{\mathrm{ff}} \approx1.3 \times10^{4}\ \mathrm{s}$, within a factor of approximately five of the observed duration. However, the stellar masses and wind parameters of IGR J00370+6122 are poorly constrained, and the adopted wind law neglects the feedback of the X-ray emission on the wind. Photoionisation suppresses line driving and would lower $v_{\mathrm{w}}$ further below the value assumed here (e.g. [14, 9]). Because $R_{\mathrm{acc}} \propto v_{\mathrm{rel}}^{-2}$ and $t_{\mathrm{ff}} \propto v_{\mathrm{rel}}^{-3}$, a modest reduction in $v_{\mathrm{w}}$ would suffice to reduce the discrepancy between the observed duration of the outburst and $t_{\mathrm{ff}}$.

2. The outburst as a sequence of depletion and replenishment cycles. Interestingly, the value $t_{\mathrm{ff}} \approx1.3 \times10^{4}\ \mathrm{s}$ derived above is comparable to the intervals between the outburst peak and the five subsequent rebrightenings, which range between $\sim10$ and $\sim13\ \mathrm{ks}$ (see the red arrows in the middle panel of Figure 3). In this scenario, the first collapse of the shell produces the initial peak, and the subsequent enhancements correspond to further collapses following the replenishment of the shell by newly captured wind material. Such a recurrence is expected in the model by Shakura et al. [45], Postnov et al. [34], in which flares repeat for as long as the enhanced plasma entry rate is sustained. We caution that we may have missed some sporadic flares during the $\sim1.4\ \mathrm{ks}$ gaps separating the $\sim4.4\ \mathrm{ks}$ visibility windows; the pattern discussed here therefore remains only suggestive.

An important feature of the outburst observed by NuSTAR and reported in this work is that it is composed of numerous flares with durations of $\sim10 - 10^{3}\ \mathrm{s}$ (Sect. 3), shorter than $t_{\mathrm{ff}}(R_{\mathrm{acc}})$ under either of the interpretations above. We speculate that these flares could reflect a non-uniform breakup of the shell, with matter leaking onto the NS surface from a few preferential regions of the magnetospheric boundary rather than simultaneously over the $4\pi$ solid angle. These flares may also explain the $\sim674\ \mathrm{s}$ and $\sim346\ \mathrm{s}$ periodicities reported by in’t Zand et al. (2007) and Uchida et al. [52], which were found only during bright flaring episodes and over a limited number of cycles. A few consecutive flares with durations and separations of this order can mimic periodicity over a few cycles.

The NuSTAR spectrum of IGR J00370+6122 is well described by an absorbed double-comptt model, a characteristic feature of X-ray pulsars accreting at low mass accretion rates (typically $L_{\mathrm{x}} \approx10^{33} - 10^{35}\ \mathrm{erg\ s^{-1}}$). Previous spectral studies of this source were limited to single power-law or cutoff power-law descriptions (in’t Zand et al. 2007; Uchida et al. 2021). The sensitivity and broad energy coverage of the present observation reveal the double-hump continuum in IGR J00370+6122 for the first time. This adds a new member to the small group of sources in which this spectral state has been resolved and extends this phenomenology to the class of intermediate SFXTs, where, to our knowledge, it has not previously been reported. At these low luminosities, X-ray pulsars are in the subcritical accretion regime, in which the emission region is a hotspot (rather than an accretion column) at the NS surface (Mushtukov et al. [30]). This double-hump spectrum is physically interpreted as arising from a strong temperature gradient in the NS atmosphere. The soft component, which peaks around 5–10 keV originates from thermal Comptonisation in the optically thick, deep atmospheric layers of the polar cap, and is dominated by extraordinary polarised photons that escape due to their reduced scattering cross section. The hard component originates from the Comptonisation of the total photon field in the hot upper atmosphere. This field consists of thermal photons that escape from the bottom and cyclotron photons produced as infalling ions decelerate and collisionally excite electrons, which then undergo multiple scatterings in the hot atmospheric plasma ([30, 51]). In our analysis, we find that the normalisation of the hard-comptt component tends to increase with luminosity more strongly than that of the soft comptt component, although the uncertainties on the normalisations do not allow a firm conclusion. A similar trend has been reported in other accreting pulsars at comparable luminosities, for example, in 1A 1118−61 (Salganik et al. [40]). Salganik et al. [40] qualitatively interpret this behaviour as the result of resonant reprocessing of atmospheric radiation within the accretion flow above the polar cap: photons emitted near the NS surface with energies close to the local cyclotron energy can scatter resonantly off the infalling Doppler-shifted electrons, exciting them to higher Landau levels whose subsequent decay regenerates photons near the cyclotron energy. This mechanism is also known as photon spawning (Araya & Harding [1]). Part of this reprocessed radiation is scattered back toward the star and re-emitted, seeding further resonant interactions. As the mass accretion rate increases, the resonant optical depth of the infalling column also increases, making this mechanism more efficient and directing more flux into the high-energy hump, while the low-energy component grows more slowly.

Conclusions

We analysed a long NuSTAR observation of IGR J00370+6122, with a net exposure of ∼ 93 ks. The observation was performed around the peak of its orbital modulation.

The observation was caught IGR J00370+6122 in a bright flaring event, providing the most detailed view obtained so far of its short-term X-ray variability. The emission is resolved into numerous short flares with durations of ∼ 10–10$^{3}$ s.

The source exhibited a dynamic range of approximately 200 over the observation, with luminosities ranging from ∼ 10$^{34}$ to ∼ 1.8 × 10$^{36}$ erg s$^{-1}$.

Throughout the observation, IGR J00370+6122 remained below $L_{\rm X} \simeq4 \times10^{36}$ erg s$^{-1}$, the luminosity at which quasi-spherical settling accretion is expected to operate. We therefore examined the flaring activity within this framework. Depending on the adopted wind velocity, the free-fall time from the outer boundary of the shell is comparable either to the duration of the entire outburst, which would then correspond to a single collapse of the shell, or to the duration of its initial bright phase, in which case the subsequent enhancements would be produced by successive cycles of depletion and replenishment of the shell. The sampling of the present observation does not allow us to distinguish between these two possibilities. In both cases, the flares are much shorter than the free-fall time. We suggest that they reflect a non-uniform break-up of the shell, with matter entering the magnetosphere from a few preferential regions rather than simultaneously over the whole magnetospheric boundary.

We searched for the ∼ 346 s and ∼ 674 s periodicities reported in the literature, using both a $Z^{2}$ search and a weighted wavelet transform, but found no significant signal at either period. This makes their interpretation as the spin period of the NS less likely and is consistent with the transient nature already noted by in’t Zand et al. (2007) and Uchida et al. [52].

The NuSTAR spectrum is well described by an absorbed double-Comptonisation model. We report the first detection of a double-hump continuum in IGR J00370+6122. This spectral state has previously been observed only in other types of X-ray pulsars accreting at low rates, and our results extend it to the class of intermediate SFXTs. The measured column density is about an order of magnitude higher than the Galactic value in the direction of the source, indicating that the observed absorption is mostly intrinsic.

Finally, using ∼ 21 yr of Swift/BAT monitoring, we refine the orbital period of the system to $P_{\rm orb} = 15.6654 \pm0.0008$ d, which is currently the most accurate measurement available for this system. We show that the properties of the flaring activity reported here are compatible with quasi-spherical settling accretion. To discriminate this scenario from magnetic or centrifugal gating mechanisms, the corotation radius and the size of the magnetosphere, which can be inferred from measurements of the NS spin period and magnetic field, are still required. Because the orbital parameters of IGR J00370+6122 are already well determined, once these two quantities are measured, the system would become one of the few SFXTs in which most of the parameters describing these models are known, and hence one of the most suitable targets for testing them.

Acknowledgements. We thank the anonymous referee for helpful comments that improved this paper. LD acknowledges funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Projektnummer 549824807. PR and EB acknowledge financial support from the Bando di Finanziamento della Ricerca Fondamentale INAF 2023 Mini Grant “Investigating SFXTs and other wind-fed HMXBs”. CM was supported by Fondazione Cariplo/Cassa Depositi e Prestiti (Grant 2023-2560, PI: Papitto). This research has made use of the NuSTAR Data Analysis Software (NuSTARDAS) jointly developed by the ASI Space Science Data Center (SSDC, Italy) and the California Institute of Technology (Caltech, USA). This research made use of Astropy, a community-developed core Python package for Astronomy ([5, 4, 3]), and of Stingray ([23] (Bachetti et al. 2024; Bachetti et al. 2024)), a community-developed spectral-timing software package for astrophysics.