© D. Bogensberger et al. 2020

1. Introduction

Low-mass transient black hole binaries usually spend most of their time in a quiescent state, with a low accretion rate, and hence a low luminosity. These periods of quiescence are interspersed by outbursts in which the accretion rate and source flux increase by many orders of magnitude (see e.g. Chen et al. 1997, Remillard & McClintock 2006). Outbursts typically last for a few months, with the source passing through several well-defined states (see e.g. Belloni et al. 2011). These are distinguished and classified via their intensity, spectral, and timing properties (Homan & Belloni 2005). These different states can be visualised on a hardness intensity diagram (HID), which shows the X-ray intensity against the flux ratio between a hard and a soft X-ray energy band. Black hole transients (BHTs) evolve through a hysteresis “q”-shaped path traversed in an anti-clockwise fashion (see e.g. Fender et al. 2004).

At the start and end of the outburst, the system is located in the lower right region of the HID. It has a low luminosity and a hard X-ray spectrum, which is dominated by a power-law component produced via Compton up-scattering of soft thermal X-ray photons from the accretion disc by high energy electrons in the corona (Gilfanov 2010). In the standard truncated disc model (Esin et al. 1997, and Poutanen et al. 1997), the geometrically thin, optically thick accretion disc is terminated at a large radius, so that the disc black body component only makes a small contribution to the X-ray spectrum. The system also features a steady, compact jet (see e.g. Corbel et al. 2001). This is known as the low hard state (LHS). As the accretion rate and luminosity increase, the truncation radius of the accretion disc starts to decrease. This causes the emission of the disc to become stronger, and the spectrum overall to be characterised by a combination of a multi-temperature black body, as well as a power-law component. This part of the outburst is classified into two types: the hard intermediate and soft intermediate states (HIMS and SIMS, respectively), which can, for example, be distinguished by their variability. The transition from the HIMS to the SIMS is marked by the steady compact jet switching off, being quenched by more than 3.5 orders of magnitude (Russell et al. 2019) after a short-lived, bright and rapidly flaring transient jet is launched (see e.g. Fender et al. 2004). Sometimes a BHT temporarily reaches an anomalous state (AS), which branches off of the HID at the SIMS, towards even higher luminosities and greater hardness.

The truncation radius continues to decrease as the luminosity starts to drop slightly. When the accretion disc extends all the way to the innermost stable orbit, the spectrum is dominated by black body emission from the accretion disc, whereas the Comptonised power-law component becomes relatively insignificant. This part of the outburst is characterised as the high soft state (HSS), during which no jet is launched (see e.g. Fender et al. 1999). Eventually, after the source has faded enough, the hardness starts to increase again, and the BHT returns to the LHS, via the SIMS and HIMS, before finally returning to quiescence. This hardening in the intermediate states occurs at a lower luminosity than the initial softening and features the re-establishment of the compact jet over a few weeks (see e.g. Russell et al. 2014). Multiple transitions between the HIMS and the SIMS usually occur in both of these horizontal regions of the HID.

The above describes the standard picture of the evolution of a BHT in an outburst, but it is by no means universal. For example, some outbursts fail to reach some of the softer states, before returning to the LHS, and then to quiescence (see e.g. Belloni & Motta 2016).

Temporal analysis of the X-ray light curves of many BHTs has revealed the existence of strong stochastic noise, along with quasi-periodic oscillations (QPOs). These appear as Lorentzian peaks in the power density spectrum (PDS), which is the distribution of squared Fourier frequencies of the light curve. Greater variability and QPO amplitude at higher X-ray energies suggests that the QPOs originate very close to the black hole; in the corona or the inner accretion flow (see e.g. Belloni et al. 1997). Accordingly, QPOs provide an additional useful window for studying this region and are essential in distinguishing different parts of an outburst. They come in two classes: low (LFQPOs), and high-frequency QPOs (HFQPOs) (Motta 2016). The LFQPOs have frequencies in the range 0.01–30 Hz, but are typically found below 10 Hz, and have been detected and studied much more than HFQPOs, which commonly have frequencies of several hundreds of Hz. The LFQPOs have themselves been grouped into three different types: A, B, and C (Wijnands et al. 1999, and for a detailed analysis, see e.g. Casella et al. 2005, or Motta 2016). The three types are distinguished via their PDS properties. PDSs are usually fitted with a sum of several Lorentzian functions (Belloni et al. 2002), which describe both the QPO and the stochastic noise continuum:

$$\begin{array}{c}\hfill L(\nu )=\frac{r^2{\nu }_0}{2\pi Q}{({(\nu -{\nu }_0)}^2+{\left(\frac{{\nu }_0}{2Q}\right)}^2)}^{-1}.\end{array}$$

Here, ν₀ is the centroid frequency, at which the Lorentzian reaches its maximum; Q = ν₀/FWHM is the quality factor, which is used to characterise the width of the Lorentzian; r describes the variability of a single Lorentzian. The temporal variability in the light curve is parameterised by the fractional root mean square variability (rms). If the PDS is normalised according to Belloni & Hasinger (1990), the rms is equal to the square root of the PDS integral: $r=\sqrt{{\int }_{-\infty }^{\infty }L(\nu )d\nu }$. The total rms of a PDS is determined by taking the square root of its properly normalised integral between two frequencies: ν₁, and ν₂, and is typically expressed as a percentage: $\mathrm{rms}=\sqrt{{\int }_{{\nu }_1}^{{\nu }_2}\left(\sum _i^n{L}_i(\nu )\right)d\nu }$, where n is the number of Lorentzian components used in the fit.

Type C QPOs are observed most frequently and are usually detected in the LHS and HIMS, although they can appear in all spectral-timing states (Motta 2016). They have frequencies spanning the entire LFQPO frequency range: 0.01 Hz ≲ ν₀ ≤ 30 Hz, are very narrow: Q ≳ 10, and have a large rms variability: 5%≲rms ≲ 20%. The strong broad-band continuum associated with type C QPOs contributes significantly towards their large rms. The continuum appears as a flat top noise extending up to a break frequency, above which it drops steeply. The break frequency has a similar value to, and correlates with the QPO frequency (Wijnands & van der Klis 1999). The type C frequency also roughly correlates with the X-ray flux, and anti-correlates with the broad-band rms (Motta et al. 2011).

Type B QPOs are the defining characteristic of the SIMS, distinguishing it from the HIMS (Belloni & Motta 2016). They are typically located at ν₀ ∼ 6 Hz in the early bright phase of the outburst but have also been found at lower frequencies, especially in the low flux horizontal branch traversed in the return to the LHS (Motta et al. 2011). They are also quite narrow: Q ≳ 6, and have low variability amplitude: rms ≲ 4%. The rms is predominantly due to the QPO peak, with a comparatively small contribution from the power-law continuum noise (see e.g. Motta 2016).

Type A QPOs are observed even less frequently than type Bs and are the hardest to detect. This is due to their broad and very shallow peak, with Q ≲ 3, at frequencies of 6 Hz ≲ ν₀ ≲ 8 Hz. Type A QPOs have the lowest variability amplitude, of rms ≲ 3%, a result of the weak QPO and a small contribution from the continuum noise. Type A QPOs do not feature harmonics, unlike types B and C. They are commonly found in the HSS (Homan & Belloni 2005).

Despite the considerable analysis of LFQPOs, there is no consensus as to their physical origin. It is thought that QPOs originate close to the black hole, but LFQPO frequencies are substantially smaller than the Keplerian orbital frequencies at these radii, which are several hundred Hz. Motta et al. (2011, 2012) argue that QPO types B and C have different physical origins, whereas types A and C could be due to the same process. They point to the different correlations between QPO frequency and power-law flux of the three QPO types, along with the observation of simultaneous type B and C QPOs in the PDS of GRO J1655-40.

The relativistic precession model (Stella & Vietri 1998) identifies the centroid frequency of LFQPOs with the frequency of Lense-Thirring precession at a single radius in the accretion disc. Within this model, pairs of HFQPOs are identified as corresponding to periastron precession and Keplerian frequency at the same radius (Stella et al. 1999).

Ingram et al. (2009, 2016), and Ingram & van der Klis (2013) describe the precessing inner flow model. In this model, the entire inner flow undergoes Lense-Thirring precession as a rigid body. A variation in inclination of the inner flow leads to a modulation of the Doppler boost, which subsequently generates the QPO.

The accretion ejection instability model (Tagger & Pellat 1999, Varnière et al. 2002, Varnière & Tagger 2002, and Varniere & Vincent 2016) hypothesises an instability with which energy and angular momentum are transported from a magnetised accretion disc to the corona, via a spiral density wave and a Rossby vortex. In this model, QPO types A, B, and C are distinguished by being produced in a relativistic, semi-relativistic, and non-relativistic regime of the instability, respectively.

Another model for the QPOs is the transition layer model (Titarchuk & Fiorito 2004). In this model, a transition layer separates the outer parts of the accretion disc, which orbit at a Keplerian frequency, from the inner regions, which have a sub-Keplerian orbital frequency. The orbital frequency at this transition layer would equal the QPO frequency.

A fourth hypothesis was developed by Cabanac et al. (2010). It considers a magneto-acoustic wave in the corona, which causes a variation in the efficiency of Comptonisation, thereby generating a QPO.

In summary, QPO properties provide an additional mode for ascertaining the physical processes occurring in the vicinity of the black hole and are an effective way of classifying the state of the source, providing information which the flux and energy spectra alone cannot supply. However, the physical origin of the QPOs is still debated, with several hypotheses positing vastly different mechanisms. While observations suggest a clear distinction between QPO types A, B, and C, it remains unknown whether these are the result of different physical processes or merely different regimes of the same fundamental mechanism.

Flip-flops, first noted by Miyamoto et al. (1991) (hereafter M91), are rare phenomena which appear as top-hat like flux variations in the light curve of some black hole transients in outburst. Flip-flops are distinguished from absorption dips by their longer duration, their top-hat like shape, and the positive correlation between flux and hardness during transitions (in contrast, dips have softer spectra at higher fluxes). Flip-flops are also seen to accompany changes in the PDS and are considered to be associated with some state transitions. Most state transitions, however, do not involve a flip-flop. Due to the marked and abrupt change in flux which defines the flip-flop, we refer to the higher and lower flux levels as the ‘bright state’ and ‘dim state’ respectively, and analyse them separately.

Out of a population of ∼60 galactic BHTs (Corral-Santana et al. 2016), only seven systems exhibited properties fitting the description of a flip-flop: GX 339-4 (M91, Nespoli et al. 2003, hereafter N03), GS 1124-683 (Takizawa et al. 1997, hereafter T97), XTE J1550-564 (Homan et al. 2001, hereafter H01, Sriram et al. 2016), XTE J1859+226 (Casella et al. 2004, hereafter C04, Sriram et al. 2013), H1743-322 (Homan et al. 2005, hereafter H05), XTE J1817-330 (Sriram et al. 2012), and MAXI J1659-152 (Kalamkar et al. 2011, hereafter K11, Kuulkers et al. 2013). There are other systems in which similar properties have been observed, such as in GRS 1915+105 (Soleri et al. 2008), or MAXI J1535-571 (Huang et al. 2018), but for which an identification as a flip-flop is not certain.

Flip-flops exhibit a large variety of different properties, not only between different systems, or different outbursts of the same system but also within a single outburst of a BHT. Some flip-flops are seen at, or very close to the peak of the outburst (M91, T97, H01, N03, and C04), while others are observed somewhat later, when the flux was dropping back down again (C04, H05, K11, and S12).

The observed flux change during a flip-flop also differs from system to system. Bright states can have an X-ray flux of between 3% and 33% greater than the flux of neighbouring dim states. Transition times have also been observed to span a large range, from fractions of a second, up to more than 1 ks. A large spread of values is also seen in the duration of each state. The time between adjacent flip-flop transitions can be anywhere between a few tens of seconds and several ks. There seems to be a relation between those two parameters, in that short-lived flip-flop states have fast transitions (M91, T97, and S12), and long-lived states have slow transitions (N03, H01, C04, and H05).

All flip-flops observed so far have involved a type B QPO in either the bright or the dim state. Therefore, flip-flops seem to require transitions between the SIMS and either the HIMS or the AS.

C04 noticed a hierarchy of QPO states in the flip-flops of XTE J1859+226, with type As occurring at the highest luminosities, type Cs at the lowest luminosities, and type Bs in between. The limiting fluxes separating the different QPO types were found to decrease exponentially with time. All observed flip-flops between QPO types B and C agree with this hierarchy (T97, H01, C04, and K11), but flip-flops between QPO types A and B disagree with it (M91, N03, and S12) more often than they agree with it (C04, and H05). Despite the fascinating phenomenology of flip-flops, a clear physical interpretation is still lacking.

At 05:37:23 UTC on February 16, 2018 (MJD 58165), the burst alert telescope on the Neil Gehrels Swift observatory (Swift/BAT, Gehrels et al. 2004) was triggered by an X-ray burst from a point source in the galactic plane, (Barthelmy et al. 2018). This X-ray transient was subsequently called Swift J1658.2-4242. This first detection was quickly followed up with an observation by the Swift X-Ray telescope (Swift/XRT), by D’Avanzo et al. (2018), and Lien et al. (2018), as well as one observation with the Mobile Astronomical System of TElescope Robots, at the Observatorio Astronomico Felix Aguilar (MASTER-OAFA, Lipunov et al. 2010), by Lipunov et al. (2018). Grebenev et al. (2018) performed subsequent examinations of data gathered by Swift/BAT (Lien et al. 2018), and by the INTErnational Gamma-Ray Astrophysics Laboratory (INTEGRAL, Winkler et al. 2003), which revealed that the source had started flaring on February 13, 2018 (MJD 58162). Lien et al. (2018) determined a more accurate location of it, at RA = 16h 58m 12.64s, Dec  =   −42°  41′ 54.4″, and found its spectrum to be well described by a highly absorbed power-law. Radio observations with the Australia Telescope Compact Array (ATCA) strongly supported the notion that Swift J1658.2-4242 is a black hole X-ray binary, rather than a neutron star X-ray binary (Russell et al. 2018).

Xu et al. (2018) described the first observation of the source by the Nuclear Spectroscopic Telescope ARray (NuSTAR, Harrison et al. 2013). They fitted the spectrum using a highly absorbed power-law with a relativistic reflection component, and a Gaussian component to fit an absorption line at 7.03 keV. They found the properties of the source to be consistent with known properties of black hole binaries and therefore described it as a black hole candidate. They also noted a type C QPO, whose centroid frequency increased continuously along with the flux, from 0.14 Hz to 0.21 Hz throughout this observation.

Jin et al. (2019) analysed the spectra detected by the Chandra X-ray observatory (Chandra, Weisskopf et al. 2000), and by the X-ray multi-mirror mission (XMM-Newton, Jansen et al. 2001). They detected a strong X-ray dust scattering halo (DSH), which significantly affected the observed source spectrum, and could be fitted with three dust layers at different distances from the BHT. Multiple methods were used to determine the distance to Swift J1658.2-4242, which was estimated to be ∼10 kpc.

Xu et al. (2019) examined the first flip-flop transition at the brightest part of the outburst, using data from both NuSTAR and XMM-Newton. They described a 45% flux decrease in ∼40 s. The lower flux level contained a type C QPO, but no QPO was detected at the higher flux level. Spectral fitting revealed only small differences in the high and low flux spectra. The strong relativistic reflection effect detected in the LHS by Xu et al. (2018) was no longer present in these observations.

Detailed timing analysis of the HIMS of Swift J1658.2-4242, using observations from the hard X-ray modulation telescope (Insight-HXMT, Zhang et al. 2018), the Neutron star Interior Composition ExploreR (NICER, Gendreau et al. 2012), and Astrosat (Singh et al. 2014) was performed by Xiao et al. (2019). They traced the increase in the QPO frequency alongside a decrease in the rms, determined the independence of QPO frequency and photon energy, and found the phase lag to be close to 0. Their results are consistent with a high inclination for Swift J1658.2-4242.

Jithesh et al. (2019) made a timing and spectral analysis of the three Astrosat observations of Swift J1658.2-4242, examining the properties of the detected QPO and its harmonics, their energy dependence, and time lag. The first observation featured a type C QPO which increased in frequency from 1.56 Hz to 1.74 Hz. The second observation featured additional flux variations reminiscent of the one observed by Xu et al. (2019). A 6.6 Hz type C QPO was only observed at low fluxes. By the final observation, the QPO frequency had dropped to 4.0 Hz. The QPO fractional rms was found to increase with photon energy. Positive time lags suggested that the QPO originated in the corona, and propagated outwards.

Beri & Altamirano (2018) studied Swift/XRT observations of Swift J1658.2-4242 from August 10, and September 28, 2018. They reported that the count rates and spectral shape suggested that Swift J1658.2-4242 had returned to quiescence by the time of these observations.

2. Observations

Following the initial detection, Swift J1658.2-4242 was observed by NuSTAR, XMM-Newton, Astrosat, Chandra, Swift, Insight-HXMT, INTEGRAL, and NICER in the X-ray band. ATCA also provided coverage in the radio band. We utilised measurements from all these instruments to obtain the best possible understanding of the phenomena occurring during this outburst.

XMM-Newton, NuSTAR, Astrosat, and Chandra provided detailed high quality information of the properties of the source at a few distinct parts of the outburst, enabling comprehensive analysis of the light curve, PDS, and energy spectrum. Swift/XRT and BAT, the Insight-HXMT low, medium, and high energy X-ray telescopes (LE, ME, and HE), as well as NICER and INTEGRAL, provided short snapshots at many different times, allowing an investigation of day-to-day variations, and an overall impression of the entire outburst, thereby complementing the sparse, yet long-duration observations by other instruments.

Out of the eight NuSTAR observations, six were performed simultaneously with XMM-Newton, and one simultaneously with Chandra. There was also an overlap between the second combined NuSTAR and XMM-Newton observation analysed here, with the second observation by Astrosat. These long-duration observations by NuSTAR, XMM-Newton, Astrosat, and Chandra are summarised in Table 1.

Swift/XRT performed 78 observations of Swift J1658.2-4242, spanning more than seven months. Swift/BAT produces daily average count rates. We used 42 observations by NICER to examine the timing evolution of the BHT. INTEGRAL obtained 151 observations of the BHT during its outburst, providing good hard X-ray coverage on shorter timescales than Swift/BAT.

Swift J1658.2-4242 was also observed with Insight-HXMT in 21 individual pointing observations, from MJD 58169 to 58196, with a total net exposure time of ∼700 ks, at low, medium, and high X-ray energies. ATCA examined the radio emission from Swift J1658.2-4242 at ten different times within the outburst. These observations are summarised in Table A.1.

3. Data analysis

3.1. X-ray data

We produced Swift/XRT light curves and spectra using the University of Leicester’s online data analysis software¹, using standard input parameters to create one spectrum for each observation. Dead pixel, pileup, and vignetting corrections were applied. The BAT light curve was obtained from the Swift/BAT transient monitor website² (Krimm et al. 2013), which bins observations made throughout each day into one data point.

NuSTAR data were reduced using the standard procedure with nupipeline and nuproducts, using the NuSTAR data analysis software (NuSTARDAS) version 1.8.0, version 4.7.8 of the CALibration DataBase (CALDB) from the high energy astrophysics science archive research center (HEASARC), as well as the general functionality of the high energy astrophysics software (HEASoft) version 6.25. In all observations, the lower-left corner of focal plane module B (FPMB) of NuSTAR was contaminated by stray light from a nearby source. To prevent this from affecting the extracted source signal, we restricted the size of the source extraction region, as there was some overlap between the source point spread function and the affected region of the CCD. We chose a circular region with a radius of 155" centred at the location of the source, for FPMB. This selection excluded virtually all the stray light seen in the CCD and included almost all the source photons. Focal plane module A (FPMA) did not suffer the same problem, so we chose a much larger source extraction region of a circle with a radius of 270" instead, also centred on the source. This different selection enabled us to test whether our choice of source extraction region for FPMB affected its spectrum. There are slight differences between the spectra of the two instruments at low energies, but they are not significant to this analysis. Slight normalisation differences between the two instruments are as expected. Background extraction regions were selected manually for each observation and focal plane module, as the area of a circle of radius 30", placed far away from the source, and any other possible contaminants in the field of view. We restricted the NuSTAR energy range to 3–79 keV. NuSTAR spectra were also rebinned with GRPPHA to contain at least 20 counts per energy bin.

Observations by XMM-Newton were taken by the pn, and MOS2 detectors of the European photon imaging camera (EPIC, Strüder et al. 2001). The EPIC-pn camera was set to Timing and Burst mode, in which the pnCCD achieves high temporal resolution, but loses spatial resolution along the y direction of the CCD. Thus, instead of a circular region, we selected individual pixel columns. As Swift J1658.2-4242 was very bright throughout all XMM-Newton observations, the entire CCD was dominated by source photons. We, therefore, used all observed photons in the creation of the light curves and spectra and did not apply a background correction. The contribution from the background was estimated by comparing the spectra in the strip of 4 × 200 pixels furthest away from the source to the spectra of the remaining 60 × 200 pixels of the CCD. The four-pixel wide column had a comparatively larger soft spectral component than the 60-pixel column. We interpret this as being primarily due to background contamination. As the four-pixel column was still dominated by the source, it would be inadvisable to use this as a background estimate. Instead, we decided to ensure minimal background contamination of the source spectra, by ignoring all energies for which the four-pixel column had a flux equal to at least 5% of the flux in the remaining pixels. This limit was exceeded at all energies below 2.7 keV, but not at any energies higher than that. Therefore, we restricted the XMM-Newton spectra to the energy range 2.7–10 keV.

We used version 16.1.0 of XMM-Newton’s Science Analysis System (SAS), implementing epchain to create event files, then using the standard procedures to produce spectra and light curves. We used epiclccorr to correct the light curves for telemetry dropouts. A few soft proton flares were detected when plotting the XMM-Newton light curves at energies exceeding 12 keV. All time intervals for which the count rate at these energies exceeded 4 counts s⁻¹ were excluded from the light curves. Photons were binned into 6 ms intervals, limited to that value by the frame time for EPIC-pn in Timing mode. We compared the count rates of the observations made in Timing mode with those made in Burst mode shortly before or afterwards, to test whether the XMM-Newton observations were suffering from pile up. Such comparisons indicated a negligible level of pileup.

Lien et al. (2018) determined that Swift J1658.2-4242 is heavily obscured, with a hydrogen column density of N_H ∼ 10²³ cm⁻². Jin et al. (2019) also detected a strong DSH around the black hole candidate. This provides a challenge for spectral fitting, as the DSH causes a significant, energy-dependent modification of the point spread function of the source. Spectra extracted from insufficiently large regions are therefore distorted, in particular at low energies. Jin et al. (2017, 2019) developed the spectral fitting model dscor to account for scattering of X-rays by a DSH along the line of sight. However, these models only work for circular source extraction regions and could therefore not be used for the XMM-Newton observations in Timing mode. To fit these spectra accurately as well, we developed a new spectral fitting model to describe the effect of the DSH on rectangular extraction regions which are not centred on the source.

We used the Chandra Grating-Data Archive and Catalogue webpage³ to calibrate the data from the Chandra high energy transmission grating (HETG) observation of Swift J1658.2-4242. Then we also extracted energy spectra for the positive and negative grating directions.

The two Astrosat observations with the large area X-ray proportional counter (LAXPC) were first reduced using the laxpc_make_event code inside the laxpcsoft package. From the resulting level-2 event files, we generated light curves using the task laxpc_make_lightcurve. The Astrosat observations had a significant background count rate, which was determined by the chain of analysis tasks laxpc_make_spectra, laxpc_make_backspectra, and finally laxpc_make_backlightcurve.

Swift J1658.2-4242 was observed with Insight-HXMT in pointed observation mode starting on MJD 58169. There are 21 available individual pointing observations with a total net exposure time of ∼700 ks. All the data were reduced following standard procedures using the Insight-HXMT data analysis software package HXMTDAS v2.01. The good time intervals (GTIs) were filtered with the screening criteria: (1) Earth elevation angle >10 degrees; (2) the value of the cutoff rigidity >8; and (3) the pointing offset angle <0.04 degrees. The events taken during satellite slews and passages through the South Atlantic Anomaly were also filtered out. For the data detected by the ME, the task megti was applied for making some corrections to the GTI file.

The Insight-HXMT count rate observed by the ME detector was affected by the contribution of nearby sources in the field of view (Xiao et al. 2019). Within the 5.5–10 keV energy range of interest, we determined the contamination from the brightest background source in the field of view, GX 340+0, to be 2.6 counts s⁻¹. By comparing the light curves of Insight-HXMT ME with those of NuSTAR, and XMM-Newton for the same energy range, we however found the best possible agreement for an additional background count rate of 3.9 counts s⁻¹. As GX 340+0 is not the only uncorrected background source in the field of view, we used this empirical best fit value for the estimation of the additional background contribution.

To improve our understanding of the hard X-ray evolution during the outburst, we also created light curves from observations by the IBIS/ISGRI instrument on INTEGRAL. We extracted the count rates of individual images, each corresponding to one pointing observation, produced using the ISDC Off-line Scientific Analysis (OSA) software version 11, for five energy bands in the range 20–200 keV.

3.2. X-ray timing and spectral analysis

X-ray observations were sorted according to whether flip-flops were detected in them. When we did not detect any flip-flops, we generated one PDS and one energy spectrum for the entire observation. But when flip-flops were detected, observations were split into individual regions, which were analysed separately. We selected intervals that started or ended at least 100 s before, or after any flip-flop transition, to be used in determining the timing, and spectral properties of each of the flip-flop states.

Medium energy X-ray light curves were extracted for the 5.5–10 keV band. The lower limit of 5.5 keV was chosen because the DSH caused different telescopes to disagree on the flux densities measured at lower energies. The upper limit was chosen because of the energy range of XMM-Newton, and the desire to maintain consistency across different instruments.

We generated PDSs using the general high-energy aperiodic timing software (GHATS) v.1.1.1 developed by T. M. Belloni, and the Stingray timing software (Huppenkothen et al. 2016). We binned individual photon detections into 10.24 ms bins and generated individual PDSs for intervals of 2¹⁰ bins. These were then combined to form averaged PDSs.

We created one PDS for every orbit within each observation by NuSTAR and Astrosat, as long as the orbit did not feature a flip-flop transition. If it did, the relevant orbit was subdivided further to separate the bright from the dim state. XMM-Newton data were binned into 2 ks segments, or shorter ones if a flip-flop transition required it. We excluded all regions with less than 500 s duration, for which precise spectral fits were hampered by low statistics.

We subtracted the Poisson shot noise from all Leahy-normalised PDSs (Leahy et al. 1983), and subsequently applied the square fractional rms normalisation (Belloni & Hasinger 1990). The PDSs were rebinned logarithmically and were then exported to XSPEC. There, we fit them using at least three different Lorentzians, at least one of which was zero-centred, following Belloni et al. (2002). We fit the main QPO peak, the broad-band continuum, and all visible harmonics. We extracted the main QPO properties (ν₀, Q, r, rms) from the fits, whenever a QPO was detected. The rms was computed by integrating fitted PDSs between 0.5 Hz, and 50 Hz, and calculating the associated source and background count rates in the interval. The values of the NuSTAR fractional rms in the QPOs were corrected for the effects of the NuSTAR dead time (van der Klis 1988). We also applied the rms correction for suppression of high-frequency power, as described by van der Klis (1988).

Spectral fits were performed using XSPEC version 12.9.1u. To model the intervening absorption we used solar abundances from Wilms et al. (2000) and the photoionisation cross-sections from Verner et al. (1996). The spectrum of the very first NuSTAR observation, with ObsID 90401307002 was also analysed but is not described here, as it was already thoroughly analysed by Xu et al. (2018), and because it differs greatly from the other NuSTAR and XMM-Newton spectra we investigated. It features a dominant power-law component, and a negligible multicolour disc black body component (Xu et al. 2018).

3.3. Radio data

We obtained radio monitoring of Swift J1658.2-4242 with ATCA, under project code C3057. Data were obtained at eight epochs from MJD 58166–58235. The majority of the radio observations were taken at central frequencies of 5.5 GHz, 9 GHz, 17 GHz, and 19 GHz, where each frequency pair (5.5/9 and 17/19 GHz) was recorded simultaneously. However, our three observations after 2018-03-05 were only taken at 5.5 GHz and 9 GHz. These frequency bands were recorded with a bandwidth of 2 GHz. Additional 8.4 GHz ATCA observations were taken on 2018-02-26 and 2018-03-15 as part of observations with the Australian Long Baseline Array (project code V456). These data were recorded with a bandwidth of 64 MHz. For all radio observations, the bandpass and flux calibration was performed using PKS 1934-638. Secondary phase calibration was done using the nearby (4.65 degrees away) calibrator J1714-397 for all observations except those taken on 2018-02-26 and 2018-03-15, which used J1713-4257 (2.9 degrees away).

The radio data were flagged and calibrated following standard procedures within the common astronomy software application (CASA, version 4.7.2; McMullin et al. 2007). Imaging was carried out using CLEAN within CASA. Due to a significant amount of diffuse emission in the field we used a Briggs roust parameter of 0 to help mitigate effects from diffuse emission within the field. Typically, flux densities were determined by fitting for a point source in the image plane. However, due to either short on-source time, or a very compact telescope configuration, flux densities of observations on 2018-03-05, 2018-03-08, and 2018-03-17 (MJD 58182, 58185, and 58194, respectively) were determined by fitting for a point source in the uv-plane using UVMULTIFIT (Martí-Vidal et al. 2014). We also used UVMULTIFIT to search for intra-observational variability within the radio observations (to identify the possibility of radio flaring associated with the X-ray flip-flops).

4. Light curves

In Fig. 1 we show the light curve of the first 60 days of the 2018 outburst of Swift J1658.2-4242 in three different energy bands. We separated this light curve into six temporal regions based on their different features. In the first region, we detected Swift J1658.2-4242 in the LHS, with the typical flux increase from quiescence. We depict the total outburst in the HIDs of Figs. 2 and 3. The set of data on the lower right of Fig. 3 is from this LHS, and shows a negative correlation between intensity and hardness, with a Pearson’s correlation coefficient of r_xy = −0.71, and a two-tailed p-value of p = 7.5 × 10⁻¹¹.

Fig. 1. Light curve of the 2018 outburst of Swift J1658.2-4242 in three different energy bands. The top figure shows the flux in an intermediate energy X-ray band of 5.5–10 keV, as observed by Swift/XRT, Insight-HXMT/ME, XMM-Newton, NuSTAR, and Astrosat. In the second graph we plot the Swift/BAT, Insight-HXMT/HE, and INTEGRAL light curves in the hard X-ray band of 20–200 keV. The radio light curves at 5.5, 8.4, 9, 17, and 19 GHz detected by ATCA are plotted in the lowest graph. We split the data temporally into six regions. The initial flux increase and spectral softening in the LHS and HIMS, the flip-flops, the intermediate HIMS sandwiched between the two flip-flop intervals, the late flip-flops, and finally another interval for which we lack data to characterise accurately. The rest of the outburst was also observed, and analysed, but we primarily focused on the regions shown here. We scaled the count rates of Astrosat, Insight-HXMT, and Swift/BAT to agree with observed fluxes measured by other telescopes. For both Astrosat, and Insight-HXMT we used the count rate in the same energy band, but the Swift/BAT light curve was obtained for the 15–150 keV range, rather than the 20–200 keV interval. We optimised the conversion factors by determining the overlap between the scaled counts and the actual fluxes measured by other instruments. Due to changes in the spectrum, a single conversion factor from count rate to flux is not accurate everywhere. This is the reason for the divergence of the Insight-HXMT and Swift/XRT light curves in the region between the two flip-flop intervals. NuSTAR, XMM-Newton, Astrosat, Insight-HXMT/ME, and Insight-HXMT/HE light curves were binned into 1 ks segments. Swift/XRT, Swift/BAT, INTEGRAL, and ATCA light curves were binned to one data point per observation. We restricted the energy range of Swift/XRT, NuSTAR, and XMM-Newton to 5.5–10 keV, to minimise the effects of the DSH, which has a significant impact on the spectra below 5.5 keV, and to ensure consistency between light curves of different instruments. We added 2.6% error to all Insight-HXMT/ME, and 1.3% error to all Insight-HXMT/HE measurements. Insight-HXMT/ME was affected by a significant background. The contribution of the background to the count rate was estimated to be 3.9 counts s−1 in the 5.5–10 keV range. The X-ray flare observed by Swift/XRT on MJD 58201 only appeared after correcting for the removed bad pixels in the image.

Fig. 2. HID of all Swift/XRT observations. Colours and shapes distinguish the different states of the outburst. Hardness is defined here as the ratio of fluxes in the two energy bands 6–10 keV, and 3–6 keV. Each data point represents one entire observation. Black squares are observations for which a state classification was not possible. The errors in the flux are shown but are almost always smaller than the symbols used for the data points.

Fig. 3. NuSTAR HID. Colours and shapes distinguish QPO types, and flip-flop detections, respectively. Hardness is defined as the ratio of count rates in the two energy bands 8–79 keV, and 3–8 keV. Data were binned into 200 s intervals. Error bars have been included but are almost always smaller than the symbols used for the data points.

Next, we observed the HIMS. The time of transition from the LHS to the HIMS is not known, but spectral analyses of NuSTAR and Swift/XRT observations indicate that it probably occurred around MJD ∼58168. Within this state, we detected an increase in the 5.5–10 keV X-ray flux, in parallel with a decrease of the 20–200 keV flux, showing the significant softening of the spectrum. This softening can also be seen in the upper right region of the HIDs of Figs. 2 and 3.

On MJD 58174, we observed the start of a major radio flare with ATCA, which reached flux densities at least 8.7 times brighter than those measured at any other time during the outburst. As the radio flux was still rising when the observation finished, the flare presumably reached even greater fluxes. Merely 200 s after the end of this ATCA observation, we measured the highest NuSTAR and XMM-Newton X-ray flux of the entire outburst, 5.4 times brighter than in the previous NuSTAR observation. This can be seen in Fig. 4a. We also observed the brightest hard X-ray flux of the outburst at the same time, in an observation by Insight-HXMT HE in the 20–200 keV band (see Fig. 1). When fitting the NuSTAR energy spectrum of this brightest state, we measured a total flux in the 3–79 keV energy range, of F₃₋₇₉ = 9.75 × 10⁻⁹ erg cm⁻² s⁻¹. This corresponds to an unabsorbed and unscattered flux of F_3−79, 0 ≈ 1.64 × 10⁻⁸ erg cm⁻² s⁻¹. Using the distance estimated by Jin et al. (2019) of ∼10 kpc, we calculated the unabsorbed and unscattered source luminosity to be L₃₋₇₉ ∼ 2 × 10³⁸ erg s⁻¹.

Fig. 4. Light curves and hardness ratios of all NuSTAR, XMM-Newton, and Astrosat observations containing flip-flops. The NuSTAR and XMM-Newton light curves are shown for the energy range 5.5–10 keV in all cases. Astrosat count rates are rescaled to the fluxes observed in the other instruments. The hardness ratios are computed as the ratio of counts in the two energy bands 8–79 keV, and 3–8 keV for the NuSTAR observations only. (a) (b) and (c) mark the three almost continuous long observations of the source within the early flip-flop interval. The light curves are binned into 100 s intervals, except for the two insets, which have time bins of 40 s, and 20 s, respectively. The hardness ratios are determined for 200 s bins. Vertical lines are plotted every 2761 s, showing that almost all flip-flop transitions occur at integer multiples of a fundamental period. This is discussed in more detail in Sect. 5.

Within the same simultaneous observation by NuSTAR and XMM-Newton in which we measured the greatest intensity of the outburst, there was a sudden flux drop of 43% within 40 s. In subsequent observations by NuSTAR, XMM-Newton, Astrosat, and Insight-HXMT/LE and ME, we saw similar flux increases or decreases, with the bright states having a 49–77% higher flux than the dim states (see Figs. 4 and 5). Transitions between the two flux levels occurred on timescales of 26–800 s (see Fig. 4 and Table 2). Unlike Xu et al. 2019, and Jithesh et al. 2019, we identify these phenomena as flip-flops, as they feature most of the characteristics of flip-flops seen before in other BHTs.

Fig. 5. Light curve from Fig. 1, zoomed into the flip-flop interval. We omitted the 20–200 keV light curve as it does not provide as much useful information as the other two light curves do, in this interval.

In Table 2, we see that the greatest ratio of bright to dim state flux between neighbouring states that we observed, is 1.77, and occurred during the first observed flip-flop transition. Subsequent flux ratios indicated a decreasing amplitude of fractional flux variation between the two states over time. The second to last of the flip-flops that we observed directly, had a flux ratio of merely 1.49. This is reminiscent of a damped oscillation.

In Figs. 4a–4c, we plotted the light curves and hardness ratios of individual NuSTAR and simultaneous XMM-Newton observations containing flip-flops. Figure 4a also features the ATCA observation of a radio flare, and Fig. 4c is complemented by an Astrosat observation. These figures show that changes in the total flux occurred coincident to changes in the hardness. Higher flux corresponded to greater hardness, but the fractional increase in hardness was smaller than the fractional increase in flux. We can, therefore, deduce that a flip-flop transition from a dim to a bright state corresponds to a flux density increase for both hard and soft X-ray energies. But the hard X-ray flux must have a larger fractional increase than the soft X-ray flux, to generate the observed change in hardness. However, this trend does not extend to very high energies, as the 20–200 keV light curve of Fig. 1 does not feature as significant fractional flux variations within the flip-flop intervals, as the 5.5–10 keV light curve does.

The hardness ratios in Fig. 4 seem to be correlated with flux variations in the bright states, and anti-correlated with flux variations in the dim states. This is verified in the best fit lines of individual regions within the NuSTAR HID of Fig. 3. Remarkably, despite being separated by several days, all bright states, and all dim states lie in the same region of the HID. We combined all the bright state measurements to find a correlation coefficient of r_xy = 0.66, and associated two-tailed p-value of p = 2.5 × 10⁻⁷, indicating a strong positive correlation between intensity and hardness. In contrast, the combined dim state measurements have correlation statistics of r_xy = −0.17, and p = 0.045.

We observed six bright states, seven dim states, and nine transitions between the two, in observations by NuSTAR, XMM-Newton, and Astrosat. In Fig. 4c, we noticed one additional instance of a significant flux difference on either side of a gap within an observation by NuSTAR, caused by its low Earth orbit, at about MJD 58181.84. The differences in the flux, timing, and spectral properties on either side of the gap, are however consistent with the flip-flop behaviour, and we consider the unobserved region to also feature a flip-flop.

In Fig. 5, we zoom into the flip-flop region, to better illustrate what is happening within it. Thanks to frequent observations by Swift and Insight-HXMT, we can narrow down the start of the flip-flop interval to MJD 58174. Interestingly, this implies that the start of the flip-flops coincides with the time of highest intermediate and hard X-ray fluxes, and the major radio flare.

On MJD 58184, following an extended bright state, which appeared to last for more than two days, we observed a final transition to a significantly lower flux. None of our observations over the subsequent 20 days showed any flip-flop activity. Therefore, it seems like this final drop down transition on MJD 58184–58185 marked the end of this interval of flip-flops, which had lasted for about 11 days.

The Insight-HXMT, and Swift/XRT light curves during the flip-flop interval revealed the existence of seemingly long (∼1 day) dim state periods. These were interspersed with shorter segments involving several flip-flop transitions, to and from bright states, which we observed to last for 2.2–19.0 ks (see Table 2). It is however likely that we missed several bright states between individual Insight-HXMT, and Swift observations. In the light curve of Fig. 5, there appear to be four shorter segments within which bright states and flip-flop transitions were observed. We obtained detailed measurements with NuSTAR, XMM-Newton, and Astrosat of three of these segments, which can be examined in greater detail in Figs. 4a–4c. The duration of these short segments in which one or more bright states were observed, seems to increase from one instance to the next (lasting approximately 0.15, 0.25, 0.5, and 3.3 days, respectively), whereas the time between them, in which no bright state was observed, appears to decrease (lasting approximately 2.5, 2.0, and 1.6 days, respectively).

In the 20–200 keV range, we detected similar up-and-down behaviour in the flux, which corresponds to, but is less pronounced than the light curve evolution observed in the 5.5–10 keV interval. Individual states are harder to identify in the hard X-ray light curve, so we focused on the 5.5–10 keV one. However, as the hardness ratios of Fig. 4 indicate, another way to single out flip-flop transitions is to monitor the ratio of the disc black body to power-law emission.

The flux change during transitions decreased gradually over time, as can be seen by the dotted lines in Figs. 1 and 5. But there is only very little variation between the flux changes of adjacent transitions. The gradual flux difference decrease observed in the flip-flops is reminiscent of the amplitude decay of a damped oscillator.

We determined the time it takes to transition between the two states, by fitting this function:

$$\begin{array}{c}\hfill y(t)=\{\begin{array}{cc}{y}_0+m_1(t-t_0)\hfill & \text{if}t\le t_0\hfill \\ y_0+\frac{y_1-y_0}{\mathrm{\Delta }t}(t-t_0)\hfill & \text{if}{t}_0<t\le t_0+\mathrm{\Delta }t\hfill \\ y_1+m_2(t-t_0-\mathrm{\Delta }t)\hfill & \text{if}t>t_0+\mathrm{\Delta }t\hfill \end{array},\end{array}$$(1)

to the region of the light curves from 100 s before the start, to 100 s after the end of each transition. In the above equation, y₀ is the flux just before the transition, and y₁ is the flux just after it; t₀ is the time when the transition starts, and the transition duration is Δt; m₁ and m₂ are the gradients of the light curve just before and after the transition. We found that the transition durations measured from observations by XMM-Newton in Burst mode were substantially longer than other measurements, which is a consequence of its larger uncertainties. Transitions from bright to dim states were found to take longer than transitions from dim to bright. Bright to dim transitions lasted for about 42–160 s, whereas dim to bright transitions only lasted for about 26–35 s.

In the inset of Fig. 4c, we indicate the presence of a short (∼90 s) drop in the light curve, with a fractional flux decrease of ∼15%, just 300 s before the transition from the bright to the dim state. We also see two long (∼600 s), almost triangular flux drops of ∼18% just before the first flip-flop transition, as can be seen in the inset of Fig. 4a. These features might be interpreted as failed transitions back down to the dim state.

After the end of the flip-flops, the variability of the light curves on timescales of days was small (see Fig. 1). We did not find any evidence of flip-flops in the two NuSTAR and simultaneous XMM-Newton observations made in this interval, and even frequent monitoring by Swift and Insight-HXMT did not reveal any significant flux changes. The flux decreased exponentially.

Both NuSTAR and XMM-Newton observations in this interval lie in almost the same region of the HID (see Fig. 3), which is however very clearly distinguished from every other observed part of the outburst. We see a negative correlation between hardness and intensity when combining the two observations, with r_xy = −0.75, and p = 2.6 × 10⁻²².

On MJD 58201, Swift/XRT detected an X-ray flare. Five days later, in the next observation by NuSTAR and XMM-Newton, we detected flip-flops again (see Fig. 6). These, however, differ quite significantly from the ones seen previously. They have much smaller flux differences, with flux ratios of between 1.22 and 1.36 between the bright and dim states in the 3–79 keV energy range.

Fig. 6. Light curves and hardness ratio of the NuSTAR and XMM-Newton observations of the late flip-flops. We used the same energy ranges and time bins as in Fig. 4. The Burst mode is not included here, as it was strongly affected by soft proton flares. Vertical lines are drawn every 2610 s, showcasing that the late flip-flops also feature a period underlying the times of transition between the two flip-flop states.

Interestingly, the X-ray flare is placed in the same region of the HID as the early flip-flops were observed in. The late flip-flops observed afterwards, however, fall into a different region, as can be seen in Fig. 2.

Unlike the early flip-flops, here we observed a positive correlation between hardness and intensity in both the bright and dim states. The same correlation applies to both, and has a correlation coefficient and p-value of: r_xy = 0.95, and p = 4.1 × 10⁻²⁰.

Due to fewer observations of Swift J1658.2-4242 during the late flip-flop interval, we cannot place stringent bounds on its start and end times. It is possible that the X-ray flare marked the start of this interval. In the days that followed the flare, we detected a greater variability between subsequent observations by Swift/XRT and BAT than in the 15 days that followed the end of the first flip-flop interval. This variability seemed to end by about MJD 58215, which is our estimate for the end of this late flip-flop interval.

We observed four bright states, four dim states, and six transitions between them. One additional transition was inferred but was not observed. Contrary to the early flip-flops, the late flip-flop transitions from bright to dim states (50–110 s) were faster than the reverse transitions (90–700 s).

The final NuSTAR observation, on MJD 58235–58236 contained a drop in the flux of about 27% during a gap in the data, followed by a gradual flux rise up to the previous level (see Fig. 7). The hardness ratio correlates with the flux for this observation (with r_xy = 0.89, and p = 4.1 × 10⁻²⁰), but the shape of the light curve is vastly different to what was seen during the early or late flip-flops. We, therefore, decided against labelling these flux changes as flip-flops.

Fig. 7. Light curve and hardness ratio of the final NuSTAR observation, starting at MJD 58235.96. We used the same energy ranges and time bins as in Fig. 4.

In the first NuSTAR observation, two absorption dips were seen (Xu et al. 2018). These are distinguished from flip-flops via their short duration, their shape in the light curve, and the increase of the hardness ratio. We did not observe dips in any of the other observation of Swift J1658.2-4242.

The major radio flare seen at the start of the early flip-flop interval died down quickly. During the next observation, ∼1.5 days after the detection of the radio flare, the flux had almost completely returned to its previous level again. We did not detect any additional radio flares at any other point during the outburst.

In Fig. 2 we show all Swift/XRT observations of the outburst, which continued until September 28, 2018 (MJD 58389). The rest of the outburst features the typical decrease in flux, followed by an increase in hardness, as can be seen in Fig. 8. As far as we can tell, Swift J1658.2-4242 never reached the HSS but remained in the intermediate states until it reached the LHS.

Fig. 8. Light curve of the entire outburst, for X-rays of 5.5–10 keV, 20–200 keV, and at radio wavelengths.

The duration of the flip-flop states seems to decrease over time. By this, we define the amount of time for which a particular bright, or dim state was maintained, from the end of the transition leading into it to the start of the transition leading out of it. In the first observation by NuSTAR, and XMM-Newton within the early flip-flop interval, only one bright and one dim state was seen. In the subsequent two observations, two bright, and two dim states were detected. And within the late flip-flop interval, a single NuSTAR observation spanned four bright, and four dim flip-flop states (see Fig. 4). In this outburst alone, we recover some of the variety of flip-flops seen in the literature. We observe a range of flip-flop amplitudes, durations, and also detect different behaviour in the time domain.

5. Discovery of an underlying clock in flip-flop transitions

Flip-flop states have vastly different durations. The time between one flip-flop transition and the next, was observed to range from anywhere between 2.56 ks to at least 65.2 ks. There does not seem to be a fixed time after which the light curve repeats itself in a periodic manner. However, in the XMM-Newton observation 0811213401, we noticed that two of the observed flip-flop states during MJD 58176–58177 had times between transitions which had a ratio very close to 2:1 (see Fig. 4b, and Table 2). The time between the first and second, and between the second and third transition in this observation, is 5590 s and 2785 s, respectively. Based on this finding, we investigated the possibility of an underlying timescale manifesting itself throughout all flip-flops. Starting from one of the transitions of this observation, we extrapolated integer multiples of the shorter of the two durations throughout the early flip-flop interval via the equation:

$$\begin{array}{c}\hfill f(n)=nP+t,\end{array}$$(2)

where t is the centre time of one particular flip-flop transition; P is a period, which we initially estimated to be P = 2785 s, and n ∈ ℤ. We found that numerous transitions lie very close to one of the f(n) times determined by this formula. Transitions do not happen at every integer n, but when they occur, they lie close to one of the f(n). The light curve does not appear to be periodic, but changes in the light curve occur at seemingly random multiples of a fundamental period, P. Due to the lack of overall periodicity, the standard methods to search for a period over which a cycle repeats itself, are not applicable here.

We used the centre times of each transition, which were determined by fitting the light curves around every transition with Eq. (1). Of the nine transitions we observed, eight are fitted well by Eq. (2) and a period of about P = 2785 s. Often, one of the f(n) times passes through the width of a transition or is just barely shy of it (see Fig. 4). However, the first transition observed by XMM-Newton in Burst mode on MJD 58181 (Fig. 4c) is not close to any of the times calculated by Eq. (2).

Excluding this anomalous transition time, and using ∼2.78 ks as our initial estimate of the period, we determined a best fit with P = 2761.0 ± 0.4 s, and t = 58174.2955 ± 0.0006 MJD, having a goodness of fit of χ² = 16.5, for 6 degrees of freedom. The results of this fit are shown in Fig. 9. The vertical lines drawn in Fig. 4 showcase the best fit results overplotted by the light curve. We also searched for a period larger than 2761.0 s, that could fit our observations better than this one, but did not find any. Shorter periods were found, but these turned out to be fractions of the 2761.0 s period.

Fig. 9. Difference between the best fit period and the observed time of the centre of each transition. We ignored the anomalous transition time in the fit, but show it on the graph nevertheless.

We next considered the significance of this detection, by investigating the likelihood of obtaining at least as good of an agreement from a sample that was not based on a fundamental period. We, therefore, simulated sets of nine random transition times within the window of our observations of the early flip-flop interval, to which we randomly assigned our measured errors in the transition times. The shorter the period is, the greater the likelihood is of finding a good agreement between randomised transition times, as the maximum delay between a transition time and an extrapolated period decreases. We, therefore, considered the probability of being able to fit for a period in the way we did, which is at least as long as the one we determined. We fitted all simulated transition times with Eq. (2), using a given input period, and having set n = 0 for the chronologically first transition. For each set of nine simulated transitions and each initial period estimate, we excluded the transition deviating the most from the best fit, leaving only eight transitions to be fitted. In doing so, we ensured that simulated and real data were treated in an identical way.

We fitted the remaining eight transitions again, to determine the closest integer n for each of them, to make sure that the ignored transition did not affect the result. Using the newly found best fitting integers n, we performed a third fit to find the optimal P and t, and to minimise the χ², which we finally compared with the best fit for the actual data. We provided a sample of initial period estimates, starting from 2761.0 s, up to the maximum period that would allow for nine transitions to occur within our set of observations. The increment in the initial period estimates was chosen in such a way as to correspond to a decrement of one in the number of multiples of the period that could be placed between the start and the end of the early flip-flop interval. In doing so, we sampled the entire range of possible periods.

Every set of simulated transition times was fit with all of these initial period estimates, and we investigated how many simulations had transition times that could be fitted equally good, or better than the actual data, using the same procedure, with at least one period in the specified range. When running 10⁵ iterations of this simulation, we found that only 3.38% of randomly generated sets of transition times fitted with a χ² of at most 16.5, to a fundamental period at least as large as 2761.0 s, when the worst fitting time was excluded from the fit.

We noticed that the single anomalous transition occurs very close to the half-way point between two of the f(n). This raises the question of whether a period of about half of our first estimate is instead the actual recurrent timescale of the flip-flop transitions. Indeed, a period of P = 1380.5 ± 0.2 s agrees well with all the early flip-flop transition times we observed, fitting with a χ² of 18.2 for 7 degrees of freedom. However, we found that in 10⁴ simulations, the probability of obtaining at least as good of an agreement between nine randomly generated transition times, for a period at least as long as this one, is 18.2%. Despite the existence of one anomaly, there is a greater significance in the existence of an underlying period of 2761.0 s in the early flip-flops. It is possible that P = 1380.5 ± 0.2 s is the actual underlying period defining the times of transition, but our data are insufficient to distinguish periodicity on these time scales from noise. We, therefore, rely on the agreement obtained for a period of 2761.0 s.

Next, we turned our attention to the late flip-flops, and the six transitions observed in them. A simple extrapolation of the best fit we had found for the early flip-flops, did not fit the transition times of the late flip-flops. Additionally, two transition times were separated by about 2560 s, which suggests that the late flip-flops cannot be fitted well with a period of 2761.0 s, and might require a different period. We investigated the possibility of a different linear relation linking the times of the late flip-flop transitions, via Eq. (2), independent of our results for the early flip-flops. The best fit to these data was obtained with P = 2610  ±  20 s and t = 58205.165  ±  0.003 MJD, and has a goodness of fit of χ² = 42.5, for 4 degrees of freedom. These best fit parameters were used to plot the vertical lines in Fig. 6. This best fit period underlying the late flip-flop transitions is about 5.5% smaller than the one found for the early flip-flops.

Using a similar simulation as before, after updating the parameters to match the observations of the late flip-flops, we found that for 5 × 10⁵ iterations, the probability that six random transition times within our late flip-flop observations, fit with a period at least as large as 2610 s, and a goodness of fit of at most 42.5, is about 7.37%. This higher probability is a consequence of the large χ² that we obtained for the best fit to the observed times of transition in this interval. The timing of transitions in the late flip-flops is less consistent than in the early flip-flops.

The detection of similar underlying periods in both intervals is intriguing. We investigated whether all transitions, in both the early and the late flip-flop intervals, could be fitted with one universal equation, of the form:

$$\begin{array}{c}\hfill f(n)=nP+A{n}^2+t,\end{array}$$(3)

where P describes the initial period at n = 0; t describes the reference time, corresponding to the time of the transition at n = 0; A describes the gradual change in the period, and is equal to $A=\frac{1}{2}\mathrm{\Pi }\dot{\mathrm{\Pi }}$, where Π(n) = P + 2nA is the period at n. We, however, do not measure the periods at any of the n, but rather the time between any two n, such as: $f(n)-f(n-1)={\int }_{n-1}^n\Pi (n)dn=P+(2n-1)A$.

By setting appropriate initial estimates for P, of the best fitting period found for the early flip-flops, and for A, such that the 5.5% difference in the periods of the early and late flip-flops could be generated within the 32-day range of early to late flip-flop observations, we discovered a fit that achieves a good agreement for all transition times in both intervals, except for the one anomalous transition in the early flip-flops. Excluding that one, the best fit parameters were found to be P = 2770.0 ± 0.6 s, t = 58174.2963 ± 0.0008 MJD, and A = (−8.70 ± 0.06) × 10⁻² s, resulting in an initial rate of period decrease equal to $\dot{\mathrm{\Pi }}(0)=(-6.28\pm 0.04)\times {10}^{-5}$. The best fit has 11 degrees of freedom and a goodness of fit of χ² = 52.6. We indicate the difference between the measured transition times, and the closest value of Eq. (3) to each transition time, using these best fitting parameters, in Fig. 10. The ability to fit both early and late flip-flop transitions could indicate a strong connection between them.

Fig. 10. Best fit of Eq. (3) to the times of transition in both the early and late flip-flops. The transition times cannot be fitted with Eq. (2), assuming a constant, and not too short period. But by adding an extra term to allow for a change in the period, we can fit the times of transition reasonably well.

We once again investigated the significance of this detection, under the assumption that the early and late flip-flops should be considered as one phenomenon, rather than two. To do so, we modified the above simulations, to instead generate 15 random transition times within our observations of early and late flip-flops. We fitted these times with Eq. (3) instead, and once again removed the worst fitting instance. We also had to specify a range of input values for A, which we limited to have an absolute value of less than 0.0870. After 10⁴ iterations, we found that the probability to obtain at least as good of an agreement in randomised data, as we had found for our observations, is 0.23%, or 2.8σ.

This analysis is based on the idea that the early and late flip-flops are the product of the same process, which evolves over time, but is not seen for about two weeks. A complementary possibility is that the early and late flip-flops are completely separate realisations of the same process, which are possibly a consequence of the major flares observed in X-ray and radio.

We, therefore, also consider the probability of obtaining a linear fit from Eq. (2) to both the early and late flip-flop intervals individually, in both instances fitting with a period of at least as large as the one we found in our observations, for equally many transitions as we observed, and having a goodness of fit smaller than the one we detected, and also having periods in both intervals that are separated by less than 5.5%. Using the results of our simulations for a constant period in the early and late flip-flops described above, we calculated the probability that any randomly chosen combination of an early and late flip-flop period, which managed to meet our fitting criteria applied to simulated data, had a separation of at most 5.5%. This probability was found to be about 26%. Therefore, the probability of randomly generated, and independent early and late flip-flop transition times being able to be fitted with a constant period in each interval, which is greater than or equal to the one we found in the data, which fit at least as well as our measurements did, and have periods separated by at most 5.5%, is P = 0.0338 × 0.0737 × 0.26 = 6.48 × 10⁻²%, or 3.2σ. This result is based on 10⁵ simulations for the early flip-flops, and 5 × 10⁵ simulations for the late flip-flops. In these simulations, we excluded the worst fitting transition, to ensure that the data and simulations were handled in the same way. We used the exact same procedure on both the real data and all simulated data. These results, therefore, suggest that this underlying period is significant.

This entire analysis is based on the 15 flip-flop transitions that we observed. There were two more instances when a flip-flop transition was inferred to have occurred in a short gap in the observations, on MJD 58181.8, and MJD 58205.5. We note that in both instances, there was one f(n) time that occurred within the gap and might have coincided with the unobserved transition. Neither of the two unobserved flip-flops contradicts the detection of a period in either the early or the late flip-flop interval. We did not use the gaps as additional measurements of times of flip-flop transition, as the associated error in their measurements is so large, that they did not aid the investigation of possible periods.

6. Power density spectra and QPOs

To learn more about the flip-flop transitions, we generated PDSs of individual flip-flop states, and entire observations without flip-flop activity, from the light curves of Astrosat, NuSTAR, and XMM-Newton. In Fig. 11 we plotted the PDSs of all the flip-flop states observed by Astrosat, as well as the PDS of the Astrosat observation in the HIMS, for comparison.

Fig. 11. Poisson noise subtracted, squared fractional rms normalised Astrosat PDSs of different intervals of the outburst. The first Astrosat observation in the LHS is marked in black, all other PDSs originate from different regions of the second Astrosat observation, which contains flip-flops (see Fig. 12). Both the LHS and Flip-flop dim state 1 PDSs have inflated widths due to being generated from a large time interval containing a variable type C QPO frequency. Generating the PDS for shorter intervals yields a narrower QPO. The same is not true for the type A QPOs seen in the bright states.

During the initial rise in flux in the outburst, we found a narrow QPO, perched on a strong broad-band continuum, featuring a second harmonic, and no subharmonics. It has properties indicative of a type C QPO (Xu et al. 2018). We noticed that the frequency of this QPO increased alongside the flux, rising from 0.13 Hz at the start of the first observation by NuSTAR, to 1.9 Hz by the end of the first Astrosat observation. This occurred alongside a flux rise from 3.9 × 10⁻¹⁰ erg cm⁻² s⁻¹ to 1.1 × 10⁻⁹ erg cm⁻² s⁻¹ in the energy range 3–79 keV.

We analysed PDSs of the flip-flop states individually and noticed a clear duality of properties, matching the duality of the flux levels. All bright states have similar PDSs, and so do all dim states. But the PDSs of the bright states differ remarkably from those of the dim states (see Fig. 11).

The dim flip-flop states detected by Astrosat, NuSTAR, and XMM-Newton all contained a reasonably narrow (4.7 ≤ Q ≤ 17) QPO at frequencies of 6–7.5 Hz (see the purple, red, and dark blue PDSs in Fig. 11). They feature a sub-harmonic, but no higher order harmonic. The rms in the QPO is 2.5–8.2%. Together with a strong, flat, broad-band continuum with a break frequency at about 3 Hz, the total rms of the dim states between 0.5 Hz and 50 Hz, was found to be 6.4–10.3%. The large effective area of Astrosat enabled a more detailed examination of the change in QPO properties on shorter timescales. We noticed that the centroid frequency of the QPO in the dim states of the second Astrosat observation scales with the count rate (Figs. 12 and 13), with r_xy = 0.78 and p = 7.3 × 10⁻¹⁹. We additionally found that the quality factor increases slightly with increasing count rate (r_xy = 0.38, p = 3.2 × 10⁻⁴). And the rms in the QPO was found to decrease with increasing count rate (r_xy = −0.51, p = 3.8 × 10⁻⁷), thereby indicating that the rms also correlates negatively with the centroid frequency (r_xy = −0.64, p = 2.0 × 10⁻¹¹).

Fig. 12. Spectrogram of the second Astrosat observation, starting on MJD 58180.786, with an overplotted light curve. Colours denote the Leahy normalised, Poisson noise subtracted power. The light curve is plotted in green, the blue curve represents the same light curve, but with a count rate reduced by 650 counts/s in the bright states. Gaps in the observations due to the low Earth orbit of Astrosat have been left out for display clarity. The time is measured in seconds of observing time since the start of the observation. The light curve is split into 200 s segments, for which the average count rate and the PDS were computed. The PDS frequencies are rebinned on a logarithmic scale. The first three state transitions shown here all lie within gaps in the data. We do however see the entire fourth transition, near the end of the observation.

Fig. 13. Correlation between the best QPO fitting parameters and the count rate, in each of the 200 s intervals of Fig. 12. We fitted each PDS between 4 and 10 Hz with a Lorentzian added to a straight line with a constant gradient. The properties of the Lorentzian were extracted to determine the centroid frequency, quality factor, and fractional rms in the QPO.

By comparing these results with the distinguishing characteristics of the three different QPO types (Motta 2016), we found that almost all of these properties suggest that this is a type C QPO. However, this QPO has a large range of quality factors (4.7 ≤ Q ≤ 17) and is usually found at Q <  10, which is too wide to fit the standard definition of a type C QPO (Motta 2016). Type C QPOs with smaller quality factors than the standard definition have however been seen before, by Motta et al. (2012), suggesting that the lower limit of the quality factor of type C QPOs is ill-defined. The shape of the continuum, the large rms, and the correlations between count rate, centroid frequency, and rms, are all inconsistent with type B and A QPOs.

Whereas the type C QPO of the dim states can easily be seen in the spectrograms of Figs. 12 and 14, no QPO can be identified during the bright states by eye. So the QPO seems to disappear entirely in dim to bright transitions. This is, however, a consequence of the low amplitude and broad profile of the QPO in the bright state, as well as the short time sampling used in these figures. Generating PDSs over longer time intervals instead revealed the existence of a very low amplitude, broad peak with a centroid frequency of about 5.0–7.2 Hz, a quality factor of 2.0–3.8, and an rms of 1.4–2.3% (see Fig. 11). Figure 12 also shows that the broad-band continuum at low frequencies is strongly suppressed in bright states. Figure 11 however indicates that the bright state continuum grows at low frequencies, such that bright and dim states have approximately equal power at 0.1 Hz. The continuum contributes slightly to the total integrated rms in bright states, which is about 2.0–2.7%. These broad, low amplitude QPOs are very difficult to detect, and could only be seen in the Astrosat PDSs. This QPO type is undoubtedly a type A QPO, because all of these fitting parameters, as well as the shape of the PDS, agree with the standard properties of a type A QPO (see e.g. N03). Even though we did not detect the type A QPO in the PDSs of NuSTAR and XMM-Newton, we found that the shape of their PDSs, and the low rms in the bright states, are consistent with the Astrosat result, of the bright states featuring type A QPOs.

Fig. 14. Spectrogram of the last flip-flop transition shown in Fig. 12. Colours denote the Leahy normalised, Poisson noise subtracted power. The light curve is plotted in green. The time is measured in seconds of observing time since the start of the orbit. We computed the PDS, and the average count rate, for each of the 20 s intervals into which we divided the data. We plotted a smaller frequency range for this figure, due to significant noise at low frequencies, a consequence of the short interval size.

Following Lewin et al. (1988), we determined the detection significance of the type A QPOs. We split the data into individual orbits, and only kept data firmly belonging to the bright flip-flop state. For each of these intervals, we determined the QPO detection significance and found it to be between 3.2 σ, for a very brief interval, and up to 11.6 σ. We, therefore, conclude that all bright states observed by Astrosat featured a type A QPO. Through the similarity in their properties, we infer that a type A QPO was also present in all other bright states of the early flip-flops.

Thanks to the large effective area and short time resolution of Astrosat, we were able to compute accurate PDSs over short time intervals as well, enabling the creation of the spectrograms shown in Figs. 12 and 14. In Fig. 14, we zoom into the first-ever detection of a direct transition from a type A to a type C QPO. Other groups (C04, Del Santo et al. 2008) have detected a type C QPO in one observation, and a type A in the next. But they could not rule out the possibility that the interval between the observations contained a transitional type B QPO. We, however, see multiple direct transitions between types A and C, and can, therefore, determine whether the change between the two QPO types involved a type B QPO. If a type B QPO had occurred during the transition, it would have featured as a very bright spike in the spectrograms, with next to no broad-band continuum, and would presumably have been found at a slightly different frequency than the type C QPO. Visually, we do not see any features matching this description in Fig. 14. We also fitted the individual PDS across the transition in the search for a type B QPO but did not find any. We even used intervals as short as 10 s without detecting any type B QPO. The search for QPOs in short time segments becomes increasingly difficult because many cycles are needed to crystallise out a QPO from a strong continuum, and noise. So even though we cannot rule it out entirely, we can rule out the existence of a type B QPO on any timescale longer than about 10 s during this transition.

We note that the change in the PDS is very sudden. In Fig. 14, we see that the type C QPO took less than 20 s to appear. Using smaller time intervals indicates that it took at most ∼10 s to get started. There was no indication of the impending change of the QPO in the PDS before the transition. This is significantly shorter than the flux transition, which instead took ≈140 s. Another interesting feature is that the QPO suddenly appears about 100 s after the start of the transition in the light curve, and ∼40 s before the flux has finished dropping down to the dim state level.

In Fig. 15, we plot the relation between the Astrosat total integrated rms and the QPO frequency. In the inset, we plot the dead time corrected rms in the QPO, in observations by NuSTAR, and XMM-Newton. Comparing this with the rms vs. ν plots of other transient BHTs (see e.g. Motta et al. 2011), we can distinguish QPO types, which fall into different regions of this graph. We see the typical negative gradient relation between the rms and centroid frequency of the type C QPOs (r_xy = −0.997, p = 1.0 × 10⁻¹⁶ for the Astrosat observations, and r_xy = −0.91, p = 1.0 × 10⁻²⁰ for the combined NuSTAR, and XMM-Newton observations), and the independence of the rms and centroid frequency of type A QPOs (r_xy = 0.23, p = 0.61). We indicate the region where we would expect the type B QPOs to appear, based on the observations of other BHTs. Compared to similar plots by for example (Motta et al. 2011), we see the clear lack of type B QPOs in this outburst. This is yet another indication that we only observed QPO types A and C, but no type B QPO.

Fig. 15. Total integrated rms and QPO centroid frequency of each orbit from the two Astrosat observations. The inset shows the centroid frequency, and the corrected rms in the QPO, detected in individual orbits of NuSTAR, and in 2 ks regions of the XMM-Newton light curves. Total rms values are computed in the 0.5–50 Hz range. We use colours and shapes to distinguish the type of QPO, and whether the measurement was made during the flip-flop period, or not.

We also tested whether the change in the timing properties could be explained by the addition of a separate flux component which dilutes the variation seen in the PDS. We added a component to the dim state Astrosat light curve, whose value in every time bin is given by a Poisson distribution with a mean value set equal to the mean counts per bin of this Astrosat dim state. This simulates a 100% flux increase by a separate feature with purely Poissonian variability. As expected, this hardly affected the overall shape of the PDS, and QPO properties. We retained the same broad-band continuum and comparatively narrow QPO. Hence, the increase in flux is not directly responsible for the change in the QPO. The source of the type C QPO, and the broadband continuum must therefore either be disrupted or blocked from being viewed, to obtain the observed changes to the PDS.

In the time between the early and late flip-flops, we once again detected type C QPOs. These were found at frequencies of 3.5–6.1 Hz but were stronger (with an rms in the QPO of ∼10%), and narrower (Q ∼ 10) than those observed in the dim flip-flop states.

We did not find any QPOs in the late flip-flops but do detect a change to the PDS in these transitions. NuSTAR suffers from a significant and variable deadtime of ∼2.5 ms, whose effect can, however, be mostly corrected for by cross-correlating the two light curves to generate a cospectrum (Bachetti et al. 2015). The total integrated rms obtained by this method differs slightly from the one that would be obtained from an ideal non-dead time affected PDS. But substantial changes to the rms are nonetheless significant. Despite the non-detection of a QPO, we still determined the total integrated rms between 0.01 and 20 Hz, of the NuSTAR cospectra in the bright and dim states of the late flip-flops. It was found to be about 8–12% in the dim states, and 6% in the bright states. We, therefore, still observe a significant change in the PDS in late flip-flop transitions, which cannot be obtained merely by the dilution of variations due to the higher count rate of the bright states. The value of the rms in the two states, as well as their position in the HID, seems to suggest that the late flip-flops correspond to transitions between the HIMS and the SIMS. However, the bright state, with the lower rms was found at a greater hardness than the dim state.

The final observation by NuSTAR in the SIMS (see Fig. 7) also featured some changes in flux. We also did not detect any QPO in this observation. The orbit with the lowest flux level had an rms of 7%, and the orbit with the highest flux had an rms of 10%. So there is some indication of changes to the PDS in this observation as well, though the difference is smaller and less significant than the rms differences found in any of the early or late flip-flops. No HFQPOs were detected at any point in the outburst.

Even though the fractional flux difference between bright and dim states decreased over time, adjacent flip-flop transitions showed next to no variation in flux difference. This is illustrated by the blue line in Fig. 12, which does not show any discontinuities at the transitions.

7. Energy spectra

7.1. Dust scattering halo and energy shifts

In Fig. 16 we plot the NuSTAR and XMM-Newton spectra of a simultaneous set of intervals from the long dim state observed on MJD 58174. We only include times when both instruments were observing Swift J1658.2-4242. As Fig. 16a shows, these spectra differ quite a lot, especially at low energies. It is impossible to fit them both using the same parameters unless only data at energies exceeding 5.5 keV is considered. This is why we restricted the medium energy X-ray light curve of Figs. 1, 4, 5, and 8 to the range of 5.5–10 keV. The difference at low energies is caused by the DSH, combined with distinct source extraction methods and regions.

Fig. 16. Absolutely simultaneous NuSTAR (black for FPMA, red for FPMB) and XMM-Newton (green) spectra from the long dim state observed on MJD 58174, showcasing the necessity of including both a DSH correction, and an energy calibration correction when fitting these X-ray spectra. Only when both corrections are applied do the spectra agree reasonably well. (a) Spectra without a DSH, and without an energy scale correction. (b) Spectra with a DSH, but without an energy scale correction. (c) Spectra with both a DSH, and an energy scale correction.

Using the DSH model we developed for the XMM-Newton pn CCD in Timing mode, and the standard DSH model by Jin et al. (2017, 2019) for NuSTAR, we managed to remove almost the entire difference between the NuSTAR and XMM-Newton spectra at low energies, as is shown in Fig. 16b. The remaining difference between the two spectra at low energies could be due to background contamination or indicate the limit of the usefulness of our DSH model at large hydrogen column densities.

Xu et al. (2018, 2019) describe an absorption line at an energy of 7.03 keV, very close to the iron K-edge at 7.112 keV. This line can be seen in Fig. 16b, and might be the signature of an accretion disc wind (Ponti et al. 2012, 2016, and Díaz Trigo & Boirin 2016). The simultaneous XMM-Newton spectra, however, feature an emission, rather than an absorption line at this energy. There is a contradiction here which needs to be resolved. As the absorption and emission lines are very close to the iron K-edge, we investigated the possibility that a wrong energy calibration at ∼7 keV could produce a fake emission or absorption line. By adding a non-zero redshift component to describe a slight shift in energy, having different values for NuSTAR and XMM-Newton, we achieved a good agreement between the two spectra at these energies.

Therefore, either the NuSTAR, or the XMM-Newton spectra, or both, require an energy recalibration. We examined the Chandra HETG spectrum, which was obtained simultaneous to the last NuSTAR observation, after the end of the late flip-flop region. As the Chandra spectrum is also strongly affected by the DSH at low energies, we restricted the fit to 5.5–10 keV, to ensure that the DSH correction would not bias our results. The NuSTAR spectrum contains a strong absorption line at ∼7 keV, as in all other observations. We fitted the possible absorption or emission line by including the additive gauss model in the spectral fit. For NuSTAR, this possible absorption line was best fitted with with a centroid energy of $7.{224}_{-0.056}^{+0.099}\mathrm{keV}$, a variance of (1.47 ± 0.85) × 10⁻¹ keV, and a normalisation of $(-2.{23}_{-0.53}^{+0.75})\times {10}^{-4}\mathrm{photons}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}$, with the negative sign indicating that this is an absorption line. This line has an associated equivalent width of $-28.5_{-6.8}^{+9.6}\mathrm{eV}$. However, the Chandra spectrum was best fitted without the inclusion of an emission or absorption line at this energy. When the additive Gaussian component was included, and the centroid energy allowed to vary, it did not fit to anywhere near to the value found by NuSTAR. When forcing the centroid energy and variance of the Gaussian included in the Chandra fit to be equal to their best fitted values in the NuSTAR spectrum, we found the best fit Gaussian normalisation to be: $(1.9_{-2.9}^{+3.2})\times {10}^{-4}\mathrm{photons}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}$, corresponding to an equivalent width of $(2.4_{-3.7}^{+4.0})\times {10}^{-2}\mathrm{keV}$. This does not support the existence of either an absorption or an emission line. Due to the higher resolution and better calibration of Chandra spectra, we therefore decided to use an energy shift in the NuSTAR and XMM-Newton spectra to ensure their consistency with the non-detection of an emission or absorption line at ∼7 keV by Chandra.

We obtained further support for this idea when fitting for the energy of the iron K-edge. We not only found that the best fit energy of the edge differs from the value it should have, of 7.112 keV, but that in using the best fit edge energy, we no longer see absorption or emission lines at these energies.

To measure the energy shift of the iron K-edge, we fitted NuSTAR and XMM-Newton spectra between 5.5–9 keV using the combination of XSPEC models: edge*diskbb. By comparing the fitted edge energy with the value it should have, of 7.112 keV, we can determine a redshift which can be implemented into a model for the entire spectrum, to correct for the energy offset observed at the iron K-edge, and ensure consistency in the spectra of NuSTAR, XMM-Newton, and Chandra.

The best fit results of the shift in iron K-edge energy from its expected value, in the NuSTAR and XMM-Newton spectra, are shown in Fig. 17. As expected, NuSTAR spectra require a negative redshift energy correction, and XMM-Newton spectra require a positive redshift energy correction. When applying these energy corrections via a redshift component, we obtain the spectra shown in Fig. 16c. We still retain slight differences between the spectra, but they are now undoubtedly more consistent than they were without the DSH correction, and without the energy calibration correction.

Fig. 17. Difference between the energy at the expected iron K-edge, at 7.112 keV, and the best fitted edge energy, in NuSTAR and XMM-Newton spectra. The energy shift was determined by fitting a simple zedge*diskbb model to spectra confined to the energy range of 5.5–9 keV.

7.2. Spectral Analysis

Comparing the spectra of the bright, dim, and late flip-flop states, we found that they are remarkably similar (see Fig. 18). Visually, besides having a different normalisation, these spectra are more similar to each other than to either of the three non-flip-flop spectra shown in this figure. The difference between the bright and dim flip-flop spectra is larger at higher energies. The hardness ratio, therefore, increases simultaneously with the flux, as seen in Fig. 4.

Fig. 18. Comparison of the NuSTAR FPMA spectra and the varying strengths of the individual components they are fit with, at different times during the outburst. The spectra here are all fitted with Model 1.

The initial LHS spectrum differs noticeably from the other spectra we measured (see Fig. 18). We also found that the spectrum in the interval between the two sets of flip-flops, and the spectrum after the completion of the late flip-flops, are visually rather similar, besides having a different normalisation. This would indicate that the spectra themselves distinguish the flip-flop intervals from other parts of the outburst. This behaviour is reminiscent of state transitions, which typically are observed to occur at a very similar hardness.

To obtain a better understanding of the differences between the bright and dim states of the flip-flop, and what distinguishes the flip-flop intervals from other parts of the outburst, we fit the spectra of individual flip-flop states, and entire observations which did not feature any flip-flops. We used two sets of XSPEC models in the spectral fitting, to obtain two independent estimates of various parameters, which we describe below.

In Model 1, we fit the spectral energy contribution from the multicolour black body component of the accretion disc, using the model diskbb. The high energy power-law component is represented by the cutoffpl model, which includes a cutoff energy, beyond which the power-law becomes steeper. We also add a diskline component to describe the iron K_α line at 6.4 keV. The effect of absorption by interstellar dust is modelled through the multiplicative ztbabs component. We shift the energy of the spectrum, to ensure the iron K-edge is located at 7.112 keV, through freezing the redshift at the energy shift found previously for each spectrum. We correct for the effects of the DSH through the multiplicative model dscor. We also used the constant multiplicative model to apply a cross-calibration correction constant between the two focal plane modules of NuSTAR. Combining all these elements, our Model 1 in XSPEC jargon is: dscor*ztbabs*constant*(diskbb+diskline+cutoffpl).

The inclination of the accretion disc around the black hole is unknown, but the dips seen at the start of the outburst, and the shape of the HID (compare Fig. 2 with the HIDs of Muñoz-Darias et al. 2013), suggest that the accretion disc around the black hole has a high inclination relative to the line of sight. To maintain consistency in the spectral fits, we set the inclination equal to 70° for all spectra. The choice of inclination does not affect the other spectral fitting parameters significantly. Unless a close to face-on inclination is chosen, the best fit parameter values agree within their respective errors. The energy of the added diskline was set equal to 6.4 keV. We were unable to fit for the inner and outer disc radii using this model, so we set these parameters equal to 6 GM/c², and 1000 GM/c², respectively. The remaining spectral parameters were left free, and their values determined when fitting the data using this Model 1. These fitting parameters are: the hydrogen column density (N_H), the inner accretion disc temperature (T_in), the emissivity power-law index of the iron line (B₁₀), the power-law photon index (Γ), the high energy cutoff (E_cut), as well as the normalisation of the three additive spectral components, diskbb, diskline, and cutoffpl.

We also fitted the energy spectra using a more physical Comptonisation model by Poutanen & Svensson (1996), compPS. As this model describes both the disc black body and the power-law component of the spectrum, we replaced both these components of Model 1 by compPS, unlike (Xu et al. 2019). Therefore, Model 2 is: dscor*ztbabs*constant*(diskline+compPS). We used the same fixed values for the iron K_α line energy, the inner and outer accretion disc radii, and the inclination of the accretion disc, which were already used in Model 1. We used a spherical corona geometry and specified a multicolour disc black body component by confining the disc temperatures specified in the model to have negative values. The absolute magnitudes of those values then represent the temperatures at the inner edge of the geometrically thin, optically thick accretion disc. We fit the spectra with a free hydrogen column density (N_H), emissivity power-law index (B₁₀), coronal electron temperature (T_e), inner accretion disc temperature (T_in), coronal optical depth (τ_y), relativistic reflection normalisation (R_r), and normalisation of the compPS, and diskline components. Spectra fitted with both Models 1 and 2 always had a lower χ² at the same number of degrees of freedom when fitted with Model 2. This indicates that the more physical model provides a better description of the data.

We fitted these two models to the NuSTAR FPMA and FPMB spectra on their own, as well as to the NuSTAR FPMA, FPMB, and XMM-Newton spectra in absolutely simultaneous intervals of each observation, or flip-flop state. After correcting for the effects of the DSH and applying an energy correction, there was a good agreement in the best fit parameter values of both sets of fits (see Fig. 16). Due to the shorter XMM-Newton exposure times, the latter selection contained a smaller number of flip-flop states and shorter intervals of simultaneous spectra. To remain consistent, and describe all the observed flip-flops seen by either of the two instruments, we only plot (Figs. 19 and 20) and list (Tables B.1, and B.2) the best fit parameter values found when fitting the NuSTAR FPMA and FPMB spectra together. In such fits, we only allow the cross-calibration constant to differ between the spectra from the two focal plane modules, along with the difference in the source extractions regions used, an important element of the DSH model. Errors are quoted at a 90% confidence interval for each parameter.

Fig. 19. Comparison of the best fit NuSTAR spectral parameters of individual observations and flip-flop states fitted with Model 1, dscor*tbabs*constant*(diskline+diskbb+cutoffpl). The spectra are distinguished by their corresponding fluxes, timing properties, and their flip-flop state classifications. The first NuSTAR observation has been excluded, as its spectrum had a negligible black body component, and therefore differed too greatly from the spectra shown here. The blue dashed line in the graph of diskline normalisation as a function of flux, depicts a line of constant equivalent width, assuming no change in spectral shape at an energy of 6.4 keV. This is a reasonably good assumption, as is demonstrated by Fig. 18. The two blue dashed lines in the diskbb normalisation graph indicate the average values of this parameter for flip-flop and non-flip-flop spectra. One dim flip-flop observation featuring a type C QPO has been omitted from the graph of the high energy cutoff for display clarity, as this parameter was not well constrained in the spectral fit of this particular observation, having a best fit of ${110}_{-50}^{+220}\mathrm{keV}$. We also excluded one bright state, and three late flip-flop data points from the graph of Betor10, as the fits were insensitive to this parameter, and the errors in the measurements could not be determined.

Fig. 20. Comparison of the best fit NuSTAR spectral parameters of individual observations and flip-flop states fitted with Model 2, dscor*tbabs*constant*(diskline+compPS). We excluded three late flip-flop data points from the graph of Betor10, because the fits were insensitive to this parameter, and the errors in the measurements could not be determined.

In Fig. 21, we analyse the differences between the bright and dim state spectra. All the spectra in this figure are fitted with Model 2. Initially, the parameters of both the bright and dim state spectra were set equal to the best fit values of the dim state spectrum, in (a). The dim state spectrum is fitted very well, but the bright state spectrum has huge residuals. Next, we changed the inner disc temperature of the spectral fit describing the bright state, to the best fit value that had been found when this spectrum was fitted on its own, in (b). This one change on its own already accounts for most of the difference between the bright and dim state spectra. So, the dominant cause of the change in flux between the dim and bright states is an increase in disc temperature, and the effect this has on the power-law component.

Fig. 21. Comparison of the bright (red) and dim (black) flip-flop spectra fitted with Model 2, with the spectral parameters initially set to the best fit values of the dim state (a). By setting the inner disc temperature of the bright flip-flop spectral fit to its best fit value, we obtain fit (b) and close most of the gap between the two spectra. In plot (c) we additionally modify the compPS normalisation of the bright state spectrum to its best fitting value. Most of the remaining difference between the two spectra is corrected for by the change in NH (d). Finally, in plot (e) we show the residuals of the best fits to both bright and dim spectra, with completely independent parameters.

In plot (c), we also change the compPS normalisation of the bright state spectral fit to its best fit value. This has now removed most of the remaining residuals of the fit, except at low energy, where the greater N_H in the bright state has a noticeable impact. When the N_H itself is changed from the dim to the bright state best fit value, in (d), the bright state spectrum is fitted very well, with only minor differences remaining. By changing all the other remaining parameters, we obtain the final graph, (e).

In Figs. 19 and 20, we distinguished the results of the spectral fits, by whether flip-flops were detected in the observations, and what QPO type was observed, if any. Interestingly, the observations without a detectable QPO are the least bright, and the bright flip-flop states with a type A QPO are the brightest. All regions in which a type C QPO was observed lie between these two extremes. All identical or comparable spectral parameters involved in the two different fitting models showed similar, though not identical values and flux dependencies when applied to the same set of data.

In the results of the spectral fits, we found N_H to correlate with flux. When using all the observations shown in Figs. 19 and 20, we found r_xy, 1 = 0.93, p₁ = 1.0 × 10⁻⁸, and r_xy, 2 = 0.83, p₂ = 1.2 × 10⁻⁵ for the two spectral fitting models, respectively. There is however a partial degeneracy between N_H and T_in, as both parameters push the diskbb spectral component to higher energies. We, therefore, analysed whether the observed change in N_H could also be reproduced by a larger change in T_in. To do so, we fit all the spectra together, each spectrum with independent fitting parameters as described above, except for the hydrogen column density, which was tied together for all spectral fits. Then for comparison, we fit all spectra again, but this time with each spectrum having its own independent hydrogen column density. The second set of fits had 18 fewer degrees of freedom but also had a χ² which was smaller by 693. We, therefore, conclude that the observed correlation between hydrogen column density and flux is a real effect. This implies that a significant fraction of the N_H must be local to the system, and able to vary within ∼2 ks. However, this goes against the assumption of Jin et al. (2019), that all the absorption occurred in the DSH.

We measured a lower average hydrogen column density for the dim flip-flop states than for the bright states. The effect of the different N_H measured in the bright and dim states can be seen in Figs. 21c and d. The low p-value for the null hypothesis which assumes a constant N_H in both states, of p = 6.4 × 10⁻⁴, and p = 1.9 × 10⁻⁴ for the two models, indicates that the hydrogen column density changes during the flip-flops. Both Xu et al. (2019) and Jithesh et al. (2019) also detected a greater N_H in the bright flip-flop states, though their values differ from the ones we found, as they did not include the DSH model in their spectral fits. As another investigation of whether the N_H changes during the flip-flops, we repeated the above test, of comparing the best fit χ² of spectral fits of flip-flop spectra, having an identical N_H, with the χ² of fits with completely independent values of N_H. We found that the decrease of 9 degrees of freedom was accompanied by a drop of the total χ², of 103. This also supports the notion that the N_H changes during flip-flop transitions.

In both models, we observed an increase in the inner disc temperature with X-ray flux, most notably between the bright and dim flip-flop states. As Fig. 21 indicates, this change in temperature is the cause of most of the difference between the bright and dim state spectra. This increase is enshrined in the detected correlation between temperature and flux of the early flip-flop bright and dim states, with r_xy, 1 = 0.92, p₁ = 1.6 × 10⁻⁴, and r_xy, 2 = 0.91, p₂ = 2.4 × 10⁻⁴ for the two models. Interestingly, the two states by themselves do not show any significant correlation at all, with p₁ = 0.51, p₂ = 0.81 for the bright states, and p₁ = 0.98, p₂ = 0.59 for the dim states. At comparable fluxes, the disc has a lower inner temperature in both of the flip-flop intervals, compared to other parts of the outburst. This is a consequence of the power-law component being stronger during the flip-flops than in other regions of the outburst whose spectra we obtained, and can be seen in Fig. 18.

On the two graphs of the dependence of the iron K_αdiskline normalisation on X-ray flux, we plotted a straight line depicting a constant equivalent width, under the assumption of a uniform spectral shape. Figure 18 indicates that even though the spectra do have some differences in their spectral shape, this is still a reasonably good assumption. We notice that the equivalent width of the bright and dim flip-flop states is very similar (46 ± 12 eV). This suggests that the disc intercepts a similar solid angle of light emitted by the primary X-ray source in both states. The interval between the early and late flip-flops features a greater equivalent width (83 ± 20 eV), and both the late flip-flops (20 ± 5 eV) and the final NuSTAR observation in the SIMS (13 ± 7 eV) have lower equivalent widths than the early flip-flops. the largest equivalent width (166 ± 8 eV) was measured in the first NuSTAR observation when the BHT was in the LHS.

The emissivity power-law index, B₁₀, does not change significantly throughout the range of observations we analysed. This was found for both sets of spectral fitting models that we used.

The normalisation of the diskbb component, is given by the equation: $N_{\mathit{dbb}}={\left(\frac{R_{\mathrm{in}}}{D_{10}}\right)}^{2}cosi$. It is proportional to the square of the accretion disc inner truncation radius (R_in), and the cosine of the disc inclination (i). The inclination of the accretion disc is not known, and the distance to the binary (D₁₀, measured in units of 10 kpc) has a significant uncertainty, so absolute measurements of the truncation radius from the best fitted values of N_dbb are very uncertain. We can, however, determine fractional changes in the inner radius, assuming a constant inclination. In Fig. 19 we observe that the diskbb normalisation is almost constant in all flip-flop and late flip-flop states, and significantly larger (by a factor of 2–4) than in observations not involving any flip-flops. All the flip-flop observations are consistent with variance around a common mean of 67.3 ± 7.1, having r_xy, 1 = −0.35, and p₁ = 0.18. In contrast, the observations without flip-flops differ too much to be considered as having the same normalisation. This, therefore, suggests that all flip-flops correspond to a similar geometry, which however differs significantly from the geometry responsible for the spectra of other observations, which also differ a lot amongst themselves.

The photon index of the power-law spectral component, Γ, does not change a lot within the spectra we considered. The early flip-flops in particular have an almost constant Γ ≈ 2.05 ± 0.12, with p = 0.39.

The high energy cutoff increases slightly with increasing flux during the early flip-flops, having p₁ = 0.02, but there is not sufficient evidence to support a change in the cutoff within the range of bright or dim state observations individually. In the last two NuSTAR observations, in which we did not detect any QPO, we found that their spectra have a higher average cutoff (83 ± 40 keV) than any of the spectra of intervals with a measurable type C or A QPO (52 ± 21 keV).

The normalisation of the cutoffpl component increases almost linearly with flux for all observations, with no clear distinction between the flip-flop and non-flip-flop states, or the nature of the QPO. Using all best fit cutoffpl values, we obtain a correlation coefficient of r_xy = 0.94, and a p-value of p = 3.6 × 10⁻⁹.

We noticed an almost constant optical depth for the bright and dim flip-flop states of τ_y ∼ 0.512 ± 0.042. But it nonetheless features a very gradual rise from dim to bright state, with r_xy = 0.67, and p = 0.032. A greater optical depth of τ_y ≈ 0.74 ± 0.16 was found in the interval between the two sets of flip-flops. The remaining observations of the late flip-flops and the SIMS, for which we did not manage to identify any QPO, had a lower average optical depth of τ_y ≈ 0.242 ± 0.057.

We observe a similar trend of the compps electron temperature, and of the cutoffpl high energy cutoff as a function of flux. This was expected, as these two parameters are related. The flip-flops states have $T_e\approx 80.2\pm 4.5k_{\mathrm{B}}^{-1}\mathrm{keV}$, and only feature a very small negative dependence of electron temperature on flux, with r_xy = −0.55, and p = 0.097. A lower electron temperature is found between the two sets of flip-flops, with an average of $T_e\approx 68\pm 11k_{\mathrm{B}}^{-1}\mathrm{keV}$. Observations of the late flip-flops, and the SIMS were found to have a higher average electron temperature of $T_e\approx 140\pm 14k_{\mathrm{B}}^{-1}\mathrm{keV}$.

Finally, the relativistic reflection component of the compps model does not indicate a dependence on flux for the early flip-flop states, which have an average of R_r ≈ 0.30 ± 0.15. This also supports the idea that the area of the disc, as seen from the central X-ray source, remains constant during flip-flop transitions. Between the early and late flip-flop intervals, we detected an increase in this component, up to R_r ≈ 0.98 ± 0.28. This is followed by a drop down to close to 0 for the late flip-flops, and the SIMS.

To study the temporal evolution of the spectral features and to determine whether these change immediately at the transition in the light curve and spectrogram, we fitted the spectra of 2 ks segments on either side of every observed flip-flop. The resulting fits indicated an immediate change of the inner disc temperature and the power-law normalisation. There was a slight suggestion of delayed hydrogen column density changes in these spectra, but changes were small, and occurred within the error margin, so a possible early or delayed change in N_H could not be corroborated.

The final NuSTAR observation still had a strong power-law component, in addition to the multicolour disc black body component. Neither of them could be excluded from the fit or was significantly dimmer than the other. This is indicative of the intermediate states. It does not agree with the HSS definition, in which the disc black body component dominates. The remaining observations by Swift/XRT needed to be fitted with both a black body and a power-law component. Therefore, it seems like Swift J1658.2-4242 never reached the HSS.

8. Discussion

We have investigated the differences between the bright and dim flip-flop states, as well as the distinctions between the flip-flop intervals, and other parts of the outburst, to enhance our knowledge of the changes during, and the causes of the extreme flip-flops we observed. Unfortunately, there is a distinct lack of models describing these phenomena.

The flip-flops we observe seem to be of an extreme and rare variety, which has never been seen before. Previously observed flip-flops had flux ratios of between 1.03 and 1.33 between the bright and dim states, and corresponded to transitions between QPO types A and B, or between types B and C. The flip-flops of Swift J1658.2-4242, however, had flux ratios of up to 1.77 and corresponded to transitions between QPO types A and C.

Flip-flops have been observed a few times before, but they are not very common, having only been seen in seven other systems so far, out of ∼60 known transient black hole binaries (Corral-Santana et al. 2016). However, it is possible that a lot of flip-flops were missed, due to short exposure times, sparse monitoring, long flip-flop state durations, short intervals of flip-flop activity, or small flip-flop amplitudes.

The duration, and frequency of XMM-Newton, and NuSTAR observations in the critical period of transitions between intermediate states in the outburst of a BHT is unusual, and what enabled the discovery of 17 flip-flop transitions, 15 of which were directly observed, and two of which were inferred. This is more than has ever been seen before within a single outburst of a BHT, of a comparable flip-flop duration. The collection of these observations are therefore ideal for examining these phenomena within one outburst, without worrying about the differences between different outbursts, and different sources. As flip-flop states of the kind we observed for Swift J1658.2-4242 last for at least 2.7 ks, but can also last for tens of ks, or much longer, short observations risk missing flip-flops entirely. Part of the apparent scarcity of flip-flops in the literature may be merely the product of shorter and fewer observations than would be required to detect flip-flops. Daily averaged count rates seem to indicate the presence of flip-flops in several other outbursts but could not be verified without direct detection through long continuous observations. Flip-flops are possibly much more common than they appear to be.

In the standard definition of spectral-timing states used to describe the outburst of a BHT (Belloni et al. 2011), the dim flip-flop states we observed are classified as HIMS, as they have both strong black body, and strong power-law spectral components, feature a type C QPO and have an rms in the range typical of the HIMS. The bright flip-flop states are classified into the AS (Belloni 2010, Motta et al. 2012), due to their location in the HID. The AS has similar timing properties as the SIMS but is found at a greater hardness, and a significantly higher flux, than the dim state HIMS. The late flip-flops are classified into the HIMS (dim state), and SIMS (bright state), due to their location in the HID, and their rms, and PDS properties.

Interestingly, we always found a negative correlation between count rate and hardness in all observations with a type C QPO (see Fig. 3). In contrast, the bright flip-flop states, which feature a type A QPO, as well as later observations in which no QPO was detected, have a positive correlation between count rate and hardness. This is predominantly a result of the relation between the path traced by a BHT in outburst in the HID, and the QPO types observed at particular points of the outburst.

Compared to previously observed flip-flops, those in Swift J1658.2-4242 fit into the sub-sample with long-duration states (of order ∼0.1–10 ks), and long transition times (of order ∼10–100 s). We did not observe any short duration flip-flops, which have durations of order ∼10 s, and transition within fractions of a second. This supports the idea that although the duration and transition time of a flip-flop can change quite a lot within one outburst of a system, they remain within about one order of magnitude of each other. No system has been found yet, in which both long and short duration flip-flips were found. This could, however, be the result of an insufficiently large sample of flip-flops.

The non-appearance of a type B QPO is at odds with previous findings by C04 according to which a flip-flop transition always involves a type B QPO. They plot a hierarchy of QPO types as a function of flux, with type A occurring at the greatest fluxes, type C at the lowest fluxes, and type B in between the two. The flux limits separating the three QPO types decayed exponentially with time. Although all flip-flops between types B and C described in the literature follow this hierarchy, several A to B flip-flops contradict it. It was thought that type B QPOs are essential for flip-flops and have to appear in either the bright or the dim state of every flip-flop. We, however, detected the first instance of flip-flops not involving a type B in either of the two states and conclude that flip-flops do not have to involve a type B QPO.

Indeed, all hard to soft state transitions were thought to require a type B QPO. Although the bright flip-flop states fit into the properties of the AS rather than the SIMS, the observed change in rms is equivalent to the standard distinction between the HIMS and SIMS. We, however, found that these kinds of state transitions do not require a type B QPO.

In Fig. 22, we plot a graph distinguishing QPO types in flux and time, analogous to Fig. 14 of C04. Plotting the temporal evolution of the maximum type C flux, and the minimum type A flux as decaying exponentials, we also found a clear gap between the two. However, we detected different decay rates of these two exponentials, whereas C04 found parallel lines in this log-linear light curve.

Fig. 22. Swift J1658.2-4242 light curve, with flux plotted on a logarithmic scale, analogous to Fig. 14 of C04. We also plotted exponentials describing the minimum type A flux, and the maximum type C flux, as a function of time. We only included NuSTAR, and XMM-Newton data, as their spectra enable an accurate flux determination, and as the QPO type can be unambiguously identified from their light curves. We removed data points during transitions between the two states to show the limiting fluxes of each state when there is no transition.

C04 observed type B QPOs in the gap between types A and C, but we did not. It should, however, be noted that they detected long-lived bright, or dim flip-flop states in this flux interval, rather than short transition periods, as we did. Perhaps the reason we did not detect any type Bs in this interval is that they take longer to switch on than the transition lasts for. The length of the transitions we found, is however longer than the observed time it took other type B QPOs in similar systems to appear or disappear. Therefore, it seems as though longer transitions would also not have involved type B QPOs.

It seems possible that there is a QPO phase diagram of parameter values, involving more than just the luminosity, which separates the different QPO states of a transient black hole, with transitions between different states happening at specific thresholds. So far, we have known of thresholds separating types A and B, and types B and C. Now we, however, know of a threshold separating types A and C directly, without a type B in between, in analogy to the solid to gaseous transition in a phase diagram.

The direct transition between QPO types A and C also supports the argument of Motta et al. (2012). They suggest that types A and C originate from the same physical process, but that type Bs have an entirely different origin.

We note that it is possible for the type A QPOs to be present in the PDSs at all times, but remain undetected whenever a type C appears, due to its very shallow amplitude and low power compared to the type C QPO continuum, and background. Either the type C replaces the type A, or the type C appears and disappears alongside changes in flux, but the type A always remains in the PDS. We could not rule out either possibility but note that a validation of either scenario could assist the interpretation of the origin of both flip-flops and QPOs.

We classify the early flip-flops of Swift J1658.2-4242 as transitions between the HIMS and the AS. Transitions between these states have not been described before, and most flip-flops observed so far would be associated with transitions between the HIMS and the SIMS. The extreme flip-flop properties observed in this system are possibly the result of transitions from the HIMS to the AS, rather than to the SIMS.

The state of the system at the time of the jet ejection, and the start of the radio flare, is unknown. The last observation with Astrosat, four days before the start of the flare was observed, found the source to be in the HIMS. The first observation after the start of the flare, by NuSTAR, was the first detected bright flip-flop state, which is classified into the AS. Even though we did not observe a peak of the radio flux, we can infer that it almost certainly reached its peak during the first flip-flop observation.

We note that the power-law involved in the flip-flops differs from the standard power-law of the LHS. The dominant power-law of the LHS is the source of its high variability, with rms values of 10–40% (Motta 2016). But during the flip-flops, we see a relative increase in the strength of the power-law and an associated sharp drop in the rms variability. This suggests that the corona generating the Comptonised power-law spectral feature has different properties during the flip-flop intervals than at the start or end of an outburst. This idea is supported by the existence of the high energy cutoff of the power-law component in the spectra.

Within the sole detailed observation of a flip-flop transition, we found flux variation to precede changes in the PDS, but take a longer time to transition. This suggests the existence of two different, yet interrelated mechanisms: one for changes in the averaged flux, and another for differences in the PDS.

The simultaneity of the radio flare with the start of the flip-flop interval implies a possible causal connection between the two effects. Jet ejections and subsequent radio flares have been linked with state transitions before (see e.g. Fender et al. 2004), and are characteristic features of X-ray binaries in outburst. However, we see no other radio flare at any later time, despite observing numerous flip-flops. We did not even see a flare on MJD 58175, when an ATCA observation coincided with a NuSTAR and XMM-Newton measurement of a bright flip-flop state. If all flip-flop transitions are related to jet ejections, we should have detected a noticeably higher radio flux in this observation.

Due to the lack of additional radio flares, our observations suggest that they are only associated with the very first flip-flop transition. This view is supported by the observations described by C04. They detect an ejection episode just before the start of the first X-ray flip-flop interval, and another one within a second flip-flop interval, 11 days later. K11 also found a rapid drop in radio emission, possibly belonging to the end of a radio flare, just before the start of the flip-flops. Therefore, the entire flip-flop behaviour might be the settling of the system back to an equilibrium state following the disruption caused by the initial jet ejection.

One might expect that a radio flare should also have occurred at the start of the late flip-flop interval, but due to a lack of radio observations in this epoch, we cannot investigate this possibility. We encourage further analysis of similar events, to determine the connection between flip-flops observed in X-rays, and radio flares, or the lack thereof.

We pose the question of whether the early and late flip-flops are two different instances of the same phenomena, or governed by entirely distinct mechanisms. The drop in the early flip-flop amplitude and state duration over time seems to continue into the late flip-flops as well. This is reminiscent of a damped oscillation, which encompasses both the early and the late flip-flops. There is a greater similarity between the spectra of the early and late flip-flops, than with the spectra obtained in the interval between the two sets of flip-flops, or after the late flip-flops have completed. For instance, the diskbb normalisation implies that the flip-flops have very similar values of $R_{\text{in}}^{2}cosi$, which are significantly larger than corresponding values in observations outside of either of the two flip-flop intervals.

However, there are some very clear distinctions between the two sets of flip-flops as well. The late flip-flops have a smaller amplitude, a greater frequency of transitions, and a shorter fundamental period underlying the timing of transitions than the early flip-flops do. All early flip-flops are found at almost the same locations on the HID, with only minor variations over several days. The late flip-flops, however, lie in an entirely different region of the HID and feature different correlations between hardness and intensity. We did not detect any QPOs in any of the late flip-flop states, but we nonetheless observed a noticeable difference in the rms of the bright and dim states. Transition times were longer on average than their earlier equivalents, but the duration of individual flip-flop states was shorter. The late flip-flop spectra also contain a noticeably smaller iron K_α equivalent width than the early flip-flops, suggesting a difference in geometry.

In both the early and the late flip-flop intervals we observed a period relating the times of transition from one state to the other. The period found in the late flip-flops is slightly smaller than the one detected in the early flip-flops. We showed that it is possible both to fit the early and late flip-flop transition times under the assumption of a decreasing period, as well as to consider them separately. The similarity in the periods found in the two intervals suggests that either the late flip-flops are a continuation of the early flip-flops, or that they are separate events, but that flip-flops within the same outburst only allow for a limited variation of the transition period. The existence of this period implies a semi-periodic mechanism, by which the system can change from one to the other state only at specific points in time, but does not have to change at every multiple of the period. The time between one transition and the next appears to be a random integer multiple of the underlying period.

The interval between the early and late flip-flops is clearly distinguished from them via its spectral, variability, and timing properties. This indicates that there is a different physical mechanism at play during it. These different properties also suggest that this interval did not feature any unobserved flip-flops. However, of the two NuSTAR and XMM observations in this interval, the latter observations show a significantly greater similarity to the flip-flop spectral and timing properties. The QPO in the first of these intra-flip-flop observations has a frequency of 3.9 Hz. But the QPO in the second of these observations has a frequency of 5.7 Hz, much closer to the QPO frequency in the dim flip-flop states observed previously, of 6.2–8.3 Hz. The fitted spectral properties also support this interpretation, with the latter intra-flip-flop observation having similar parameter values as the late flip-flops, as can be seen in Figs. 19 and 20. This possibly indicates a gradual return to the set of conditions enabling flip-flops to occur.

We notice similarities in the light curves of some, but not all, black hole transients featuring flip-flops. Swift J1658.2-4242, XTE J1859+226 (C04, and Sriram et al. 2013), H1743-322 (H05), possibly also the 2002-2003 outburst of GX 339-4 (N03), all feature an initial flip-flop region at the brightest part of the outburst, which is followed by several days of almost constant flux, before a late flip-flop interval is initiated. Flip-flops in the latter interval were only detected and described for Swift J1658.2-4242, and XTE J1859+226, but comparable regions in the light curves of other black hole transients appear very similar and might also contain flip-flops. None of these systems exhibits a third region that appears to contain flip-flops. These similarities suggest that flip-flops could be linked to a particular type of evolution of the outburst, which only occurs for some BHTs.

An intriguing distinction between the times when flip-flops occurred, and when they did not, can be seen in the fitted truncation radius of the accretion disc. Figure 19 points to a consistent value of the diskbb normalisation for all flip-flop states, both within the early, and the late flip-flop interval. Other observations within a similar part of the outburst, which do not feature flip-flops, have a range of different values for this normalisation, but these are all a factor of 2–4 smaller than what was found within the flip-flop intervals. The diskbb normalisation is defined as: $N_{\mathit{dbb}}={\left(\frac{R_{\mathrm{in}}}{D_{10}}\right)}^{2}cosi$, where R_in is the inner radius of the accretion disc, D₁₀ is the distance to the source, in units of 10 kpc, and i is the inclination. The distance to the system can be assumed to be constant, so as long as the inclination does not vary a lot, the distinction of the diskbb normalisations is due to a difference in the truncation radii. This result would, therefore, indicate that the disc is truncated at a very specific radius whenever flip-flops occur. Following the same reasoning, other times of the outburst, which do not feature flip-flops, have no consistency in their truncation radii.

The fitted diskbb normalisation falls into the range of 17–74 within the NuSTAR spectra of observations after the initial flux rise. Even though we do not know the exact inclination of the accretion disc relative to the line of sight, we assumed it to be ∼70°. Along with the estimate of the distance to the system, of ∼10 kpc (Jin et al. 2019), we can approximate the range of diskbb normalisations to correspond to accretion disc truncation radii in the range of 12–25 km. We used the conversion by Kubota et al. (1998), to obtain the true truncation radius from the spectrally fitted one. This correction increases the radius by 19–65%, depending on the choice of κ, and we chose the maximum of κ = 2.0. The resulting radii are still too small compared with the gravitational radius of a 10 M_⊙ black hole, of 14.9 km. The small spectrally fitted inner truncation radii could be a consequence of an overestimated black hole mass, or an underestimated inclination, or distance to the black hole. But all of these quantities cannot differ a lot from their estimated values. Tomsick et al. (2009), however, points out that mixed results were obtained from using the spectrally fitted normalisation of the disc black body component to estimate the inner radius. Additionally, in the relativistic rigid body precession model of the QPO, its frequency is determined by the truncation radius of the accretion disc. A larger truncation radius should yield a lower QPO frequency. However, the highest QPO frequencies were observed in the dim flip-flop states, when the spectrally fitted truncation radius was at its largest. The accretion disc inner radius as a function of time remains uncertain but is nonetheless seemingly very close to the black hole in the brightest phase of the outburst, particularly whenever there is no indication of flip-flop activity.

We found that the major spectral difference between the dim and bright flip-flop states is a result of an increase of the inner disc temperature, at a constant inner radius. An additional increase of the power-law normalisation and the hydrogen column density are required to explain most of the remaining change in flux and spectral shape between the dim and bright flip-flop states.

The flip-flop interval features a decrease in flip-flop amplitude over time, and this decrease seems to continue into the late flip-flop phase, for which both the absolute and fractional flux difference is substantially smaller than it was in the early flip-flop interval. This suggests an association of flip-flops with a damped oscillation with an irregular period. In this picture, the mechanism responsible gradually loses its amplitude and approaches an equilibrium state. Perhaps the jet ejection triggered an imbalance in the system which initiated an instability that pushed the system from one unstable state to another one, resulting in an oscillation between the two unstable states.

The physical mechanism causing flip-flop to occur remains unknown. The association of flip-flops with alterations of QPOs suggest that an understanding of the former might require a knowledge of the latter. The search for a model describing the cause of flip-flops is complicated by the ongoing discussion about the origin of QPOs. However, as our findings of transitions between QPO types A and C demonstrated, an understanding of flip-flops could aid in the search for the origin of QPOs as well. Below, we present some ideas to provide the reader with a flavour of some of the ingredients that could play a role in shaping the flip-flop phenomenology.

Flip-flops are described by two timescales; the duration of a state before the next transition and the time the BHT takes to change from one state to the other. For Swift J1658.2-4242, these timescales are of order ∼10 ks, and ∼100 s, respectively. Other BHTs feature flip-flops with a wide variety of different timescales, but the ratio between the two seems to be reasonably similar. This could suggest that the two timescales are related to each other, and the physical parameters of the disc, and the black hole. Indeed, this ratio is consistent with the expected ratio of viscous to thermal timescales within the inner region of a standard accretion disc.

A thermal instability could also account for the observed temperature variation observed during flip-flop transitions. On the other hand, we note that thermal-viscous instabilities typically result in much larger temperature variations (see e.g. Fragile et al. 2018) than we observed in the flip-flops (of at most 20%). It also remains unclear how this instability could generate the observed change of the QPO or the underlying clock in the times of transition.

A completely different framework posits the existence of an ionised, relativistic, beamed outflow misaligned with the spin of the black hole, therefore precessing with a complex pattern. Our discovery of an underlying clock governing the times of flip-flop transitions indicates the existence of a fundamental complex periodicity in the system, which would be achieved by this model. An increase in flux from dim to bright states would be the result of bulk Comptonisation of disc and corona emission whenever the outflow points into the line of sight. The bulk Comptonisation of this outflow might explain the observed spectral changes between dim and bright states, namely the increase in temperature and power-law normalisation (eventually also the increased hydrogen column density, if the outflow is multi-phase) in the bright states. To explain the variation in the PDS we would require this outflow to have a significant radial extent, opening angle, and optical depth, to ensure that a narrow type C QPO would be broadened and weakened into a type A QPO due to scattering of radiation, and the corresponding time lags across this outflow.

This model provides a straightforward explanation for the near-consistency of the flux differences of adjacent flip-flop transitions. It also immediately justifies why the type C QPO that appears after the end of a bright state is so similar to the type C QPO observed just before the start of the bright state; it was there throughout but was scattered to a type A QPO in the bright state. The pause in flip-flop activity would be explained by the outflow never pointing into the line of sight for a prolonged period. Eventually, the complex system of rotations would point the outflow back into the line of sight once in a while, and the flip-flops would be seen yet again, but this time with a lower amplitude of flip-flop flux variation. Alternatively, the pause in flip-flop activity could be the result of the outflow being quenched, and subsequently re-established.

We tried to model this hypothesis with a precessing conical outflow having a constant opening angle. We defined the system to be in a bright flip-flop state, whenever the angle between the line of sight and the outflow centre was less than the opening angle. Larger angles of separation were instead assigned to the dim state. We used this binary model to fit the occurrence of bright and dim states in our observations and compared our results with those found in randomly generated flip-flop states, which had the same number of transitions within our observing times. We were unable to fit our observations well without adding even more complexity, suggesting that this simplistic scenario can also not explain the wealth of observational constraints. Using this model, it was possible to fit randomly generated flip-flop states equally well in 45% of our simulations. Our simulations, therefore, do not support this particular interpretation with these specific assumptions. It is, however, possible that a more complex system of rotations is required, or that the assumptions we used, in particular of a constant opening angle, are too simplistic. We also note that this model cannot account for the flip-flops observed in other systems, with transitions between different sets of QPO types.

In conclusion, some of the requirements of this model seem to be rather extreme, even though they are still within the realms of physical possibility. A clear verifiable prediction of this model is a high degree of polarisation in the bright flip-flop states, which could be investigated by future X-ray polarimetric observations by IXPE or eXTP. This model should be regarded as tentative, and incomplete. Finally, we note that it is possible for flip-flops to be caused by a combination of several different effects, rather than one single effect, as we have considered here.

9. Conclusions

We observed the black hole candidate Swift J1658.2-4242 throughout its outburst from February to September 2018, using a suite of X-ray and radio instruments. Swift J1658.2-4242 underwent some extreme changes in its light curve and PDS, which have never been seen before.

Flip-flops with flux ratios of up to 1.77 were observed, which occurred simultaneously to changes between QPO types A and C. We report the first detection of a direct transition between these two QPO types. The presence of a type B QPO with a duration of longer than ∼10 s is ruled out in the transitions from QPO types A to C. The start of a type C QPO during a bright to dim transition was observed to be delayed relative to the start of the transition in the light curve, but the QPO changed significantly faster.

A major radio flare was detected at the same time as the greatest X-ray flux was reached, and the flip-flop interval started, suggesting that these events are related. We found a second interval of flip-flop activity, which featured smaller flux ratios between the two states, a greater frequency of transitions, no identifiable QPO, but a significant change in total rms between the bright and dim flip-flop states. It started about 16 days after the end of the first flip-flop interval. But only one radio flare was observed throughout the entire outburst.

We found flip-flop transitions at seemingly random integer multiples of 2.761 ks after the time of the first flip-flop. A slightly lower period of 2.61 ks was observed in the late flip-flops. The existence of an underlying period in both intervals, which are separated by no more than 5.5%, was detected with a significance of 3.2σ. This suggests an inherent base timescale, which defines when flip-flop transitions can occur.

Spectral fitting of flip-flop states revealed that the major increase in flux and the slight change in their spectral shape was primarily the result of an increase in the inner accretion disc temperature. The remaining disparity was caused by an additional increase of the power-law normalisation and hydrogen column density. We found that the inner accretion disc radius remained almost constant throughout all flip-flop states and was considerably larger than at times without any flip-flops.

We demonstrated the necessity of using a DSH model to describe the spectrum of Swift J1658.2-4242 and other highly obscured sources. Energy corrections were applied to NuSTAR and XMM-Newton spectra to ensure consistency between them, and the Chandra spectrum.

We consider a misaligned, precessing, and beamed outflow as a tentative explanation for the phenomenology we observed. In this model, flip-flops would be generated at semi-periodic times whenever the outflow passes into the line of sight.

Flip-flops have rarely been observed so far, but this might reflect more on the nature of past observations, rather than their fundamental rarity, as long continuous observations are required to detect them. We encourage further multiwavelength analysis of similar phenomena in other systems and a continued search for the physical nature of flip-flops.

Acknowledgments

DB acknowledges financial support from the Bundesministerium für Wirtschaft und Energie (BMWi)/Deutsches Zentrum für Luft- und Raumfahrt (DLR) grant FKZ 50 OR 1812. TMB and GP acknowledge financial contribution from the agreement ASI-INAF n.2017-14-H.0. GP acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. [865637]). CJ acknowledges the National Natural Science Foundation of China through grant 11873054, and the support by the Strategic Pioneer Program on Space Science, Chinese Academy of Sciences through grant XDA15052100. HP acknowledges the National Natural Science Foundation of China through grant 11803047 and U1838202. TDR is supported by a Netherlands Organisation for Scientific Research Veni award. TMD acknowledges support from the Spanish grants AYA2017-83216-P and RYC-2015-18148. PV acknowledges partial funding support from the DST-INSPIRE fellowship from the Government of India. FV acknowledges support from STFC under grant ST/R000638/1. This research has made use of data and/or software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC and the High Energy Astrophysics Division of the Smithsonian Astrophysical Observatory. This research has made use of the NuSTAR Data Analysis Software (NuSTARDAS) jointly developed by the ASI Science Data Center (ASDC, Italy) and the California Institute of Technology (Caltech, USA). This publication uses the data from the AstroSat mission of the Indian Space Research Organisation (ISRO), archived at the Indian Space Science Data Centre (ISSDC). This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester. This work made use of the data from the Insight-HXMT mission, a project funded by China National Space Administration (CNSA) and the Chinese Academy of Sciences (CAS). This research has made use of data obtained from the Chandra Data Archive and the Chandra Source Catalog, and software provided by the Chandra X-ray Center (CXC) in the application packages CIAO, ChIPS, and Sherpa. This research is based on observations with INTEGRAL, an ESA project with instruments and science data centre funded by ESA member states (especially the PI countries: Denmark, France, Germany, Italy, Switzerland, Spain), and with the participation of the Russian Federation and the USA. This paper includes archived data obtained through the Australia Telescope Online Archive (http://atoa.atnf.csiro.au). The Australia Telescope Compact Array is part of the Australia Telescope National Facility which is funded by the Australian Government for operation as a National Facility managed by CSIRO. The authors thank the anonymous referee for providing useful comments which improved the paper.

References

Appendix A: Table of ATCA radio observations

Appendix B: Spectral fitting results

All Tables

All Figures

Fig. 1. Light curve of the 2018 outburst of Swift J1658.2-4242 in three different energy bands. The top figure shows the flux in an intermediate energy X-ray band of 5.5–10 keV, as observed by Swift/XRT, Insight-HXMT/ME, XMM-Newton, NuSTAR, and Astrosat. In the second graph we plot the Swift/BAT, Insight-HXMT/HE, and INTEGRAL light curves in the hard X-ray band of 20–200 keV. The radio light curves at 5.5, 8.4, 9, 17, and 19 GHz detected by ATCA are plotted in the lowest graph. We split the data temporally into six regions. The initial flux increase and spectral softening in the LHS and HIMS, the flip-flops, the intermediate HIMS sandwiched between the two flip-flop intervals, the late flip-flops, and finally another interval for which we lack data to characterise accurately. The rest of the outburst was also observed, and analysed, but we primarily focused on the regions shown here. We scaled the count rates of Astrosat, Insight-HXMT, and Swift/BAT to agree with observed fluxes measured by other telescopes. For both Astrosat, and Insight-HXMT we used the count rate in the same energy band, but the Swift/BAT light curve was obtained for the 15–150 keV range, rather than the 20–200 keV interval. We optimised the conversion factors by determining the overlap between the scaled counts and the actual fluxes measured by other instruments. Due to changes in the spectrum, a single conversion factor from count rate to flux is not accurate everywhere. This is the reason for the divergence of the Insight-HXMT and Swift/XRT light curves in the region between the two flip-flop intervals. NuSTAR, XMM-Newton, Astrosat, Insight-HXMT/ME, and Insight-HXMT/HE light curves were binned into 1 ks segments. Swift/XRT, Swift/BAT, INTEGRAL, and ATCA light curves were binned to one data point per observation. We restricted the energy range of Swift/XRT, NuSTAR, and XMM-Newton to 5.5–10 keV, to minimise the effects of the DSH, which has a significant impact on the spectra below 5.5 keV, and to ensure consistency between light curves of different instruments. We added 2.6% error to all Insight-HXMT/ME, and 1.3% error to all Insight-HXMT/HE measurements. Insight-HXMT/ME was affected by a significant background. The contribution of the background to the count rate was estimated to be 3.9 counts s−1 in the 5.5–10 keV range. The X-ray flare observed by Swift/XRT on MJD 58201 only appeared after correcting for the removed bad pixels in the image.

In the text

Fig. 2. HID of all Swift/XRT observations. Colours and shapes distinguish the different states of the outburst. Hardness is defined here as the ratio of fluxes in the two energy bands 6–10 keV, and 3–6 keV. Each data point represents one entire observation. Black squares are observations for which a state classification was not possible. The errors in the flux are shown but are almost always smaller than the symbols used for the data points.

In the text

	Fig. 3. NuSTAR HID. Colours and shapes distinguish QPO types, and flip-flop detections, respectively. Hardness is defined as the ratio of count rates in the two energy bands 8–79 keV, and 3–8 keV. Data were binned into 200 s intervals. Error bars have been included but are almost always smaller than the symbols used for the data points.
In the text

Fig. 4. Light curves and hardness ratios of all NuSTAR, XMM-Newton, and Astrosat observations containing flip-flops. The NuSTAR and XMM-Newton light curves are shown for the energy range 5.5–10 keV in all cases. Astrosat count rates are rescaled to the fluxes observed in the other instruments. The hardness ratios are computed as the ratio of counts in the two energy bands 8–79 keV, and 3–8 keV for the NuSTAR observations only. (a) (b) and (c) mark the three almost continuous long observations of the source within the early flip-flop interval. The light curves are binned into 100 s intervals, except for the two insets, which have time bins of 40 s, and 20 s, respectively. The hardness ratios are determined for 200 s bins. Vertical lines are plotted every 2761 s, showing that almost all flip-flop transitions occur at integer multiples of a fundamental period. This is discussed in more detail in Sect. 5.

In the text

	Fig. 5. Light curve from Fig. 1, zoomed into the flip-flop interval. We omitted the 20–200 keV light curve as it does not provide as much useful information as the other two light curves do, in this interval.
In the text

Fig. 6. Light curves and hardness ratio of the NuSTAR and XMM-Newton observations of the late flip-flops. We used the same energy ranges and time bins as in Fig. 4. The Burst mode is not included here, as it was strongly affected by soft proton flares. Vertical lines are drawn every 2610 s, showcasing that the late flip-flops also feature a period underlying the times of transition between the two flip-flop states.

In the text

	Fig. 7. Light curve and hardness ratio of the final NuSTAR observation, starting at MJD 58235.96. We used the same energy ranges and time bins as in Fig. 4.
In the text

	Fig. 8. Light curve of the entire outburst, for X-rays of 5.5–10 keV, 20–200 keV, and at radio wavelengths.
In the text

	Fig. 9. Difference between the best fit period and the observed time of the centre of each transition. We ignored the anomalous transition time in the fit, but show it on the graph nevertheless.
In the text

	Fig. 10. Best fit of Eq. (3) to the times of transition in both the early and late flip-flops. The transition times cannot be fitted with Eq. (2), assuming a constant, and not too short period. But by adding an extra term to allow for a change in the period, we can fit the times of transition reasonably well.
In the text

Fig. 11. Poisson noise subtracted, squared fractional rms normalised Astrosat PDSs of different intervals of the outburst. The first Astrosat observation in the LHS is marked in black, all other PDSs originate from different regions of the second Astrosat observation, which contains flip-flops (see Fig. 12). Both the LHS and Flip-flop dim state 1 PDSs have inflated widths due to being generated from a large time interval containing a variable type C QPO frequency. Generating the PDS for shorter intervals yields a narrower QPO. The same is not true for the type A QPOs seen in the bright states.

In the text

Fig. 12. Spectrogram of the second Astrosat observation, starting on MJD 58180.786, with an overplotted light curve. Colours denote the Leahy normalised, Poisson noise subtracted power. The light curve is plotted in green, the blue curve represents the same light curve, but with a count rate reduced by 650 counts/s in the bright states. Gaps in the observations due to the low Earth orbit of Astrosat have been left out for display clarity. The time is measured in seconds of observing time since the start of the observation. The light curve is split into 200 s segments, for which the average count rate and the PDS were computed. The PDS frequencies are rebinned on a logarithmic scale. The first three state transitions shown here all lie within gaps in the data. We do however see the entire fourth transition, near the end of the observation.

In the text

	Fig. 13. Correlation between the best QPO fitting parameters and the count rate, in each of the 200 s intervals of Fig. 12. We fitted each PDS between 4 and 10 Hz with a Lorentzian added to a straight line with a constant gradient. The properties of the Lorentzian were extracted to determine the centroid frequency, quality factor, and fractional rms in the QPO.
In the text

Fig. 14. Spectrogram of the last flip-flop transition shown in Fig. 12. Colours denote the Leahy normalised, Poisson noise subtracted power. The light curve is plotted in green. The time is measured in seconds of observing time since the start of the orbit. We computed the PDS, and the average count rate, for each of the 20 s intervals into which we divided the data. We plotted a smaller frequency range for this figure, due to significant noise at low frequencies, a consequence of the short interval size.

In the text

Fig. 15. Total integrated rms and QPO centroid frequency of each orbit from the two Astrosat observations. The inset shows the centroid frequency, and the corrected rms in the QPO, detected in individual orbits of NuSTAR, and in 2 ks regions of the XMM-Newton light curves. Total rms values are computed in the 0.5–50 Hz range. We use colours and shapes to distinguish the type of QPO, and whether the measurement was made during the flip-flop period, or not.

In the text

Fig. 16. Absolutely simultaneous NuSTAR (black for FPMA, red for FPMB) and XMM-Newton (green) spectra from the long dim state observed on MJD 58174, showcasing the necessity of including both a DSH correction, and an energy calibration correction when fitting these X-ray spectra. Only when both corrections are applied do the spectra agree reasonably well. (a) Spectra without a DSH, and without an energy scale correction. (b) Spectra with a DSH, but without an energy scale correction. (c) Spectra with both a DSH, and an energy scale correction.

In the text

	Fig. 17. Difference between the energy at the expected iron K-edge, at 7.112 keV, and the best fitted edge energy, in NuSTAR and XMM-Newton spectra. The energy shift was determined by fitting a simple zedge*diskbb model to spectra confined to the energy range of 5.5–9 keV.
In the text

	Fig. 18. Comparison of the NuSTAR FPMA spectra and the varying strengths of the individual components they are fit with, at different times during the outburst. The spectra here are all fitted with Model 1.
In the text

Fig. 19. Comparison of the best fit NuSTAR spectral parameters of individual observations and flip-flop states fitted with Model 1, dscor*tbabs*constant*(diskline+diskbb+cutoffpl). The spectra are distinguished by their corresponding fluxes, timing properties, and their flip-flop state classifications. The first NuSTAR observation has been excluded, as its spectrum had a negligible black body component, and therefore differed too greatly from the spectra shown here. The blue dashed line in the graph of diskline normalisation as a function of flux, depicts a line of constant equivalent width, assuming no change in spectral shape at an energy of 6.4 keV. This is a reasonably good assumption, as is demonstrated by Fig. 18. The two blue dashed lines in the diskbb normalisation graph indicate the average values of this parameter for flip-flop and non-flip-flop spectra. One dim flip-flop observation featuring a type C QPO has been omitted from the graph of the high energy cutoff for display clarity, as this parameter was not well constrained in the spectral fit of this particular observation, having a best fit of ${110}_{-50}^{+220}\mathrm{keV}$. We also excluded one bright state, and three late flip-flop data points from the graph of Betor10, as the fits were insensitive to this parameter, and the errors in the measurements could not be determined.

In the text

	Fig. 20. Comparison of the best fit NuSTAR spectral parameters of individual observations and flip-flop states fitted with Model 2, dscor*tbabs*constant*(diskline+compPS). We excluded three late flip-flop data points from the graph of Betor10, because the fits were insensitive to this parameter, and the errors in the measurements could not be determined.
In the text

Fig. 21. Comparison of the bright (red) and dim (black) flip-flop spectra fitted with Model 2, with the spectral parameters initially set to the best fit values of the dim state (a). By setting the inner disc temperature of the bright flip-flop spectral fit to its best fit value, we obtain fit (b) and close most of the gap between the two spectra. In plot (c) we additionally modify the compPS normalisation of the bright state spectrum to its best fitting value. Most of the remaining difference between the two spectra is corrected for by the change in NH (d). Finally, in plot (e) we show the residuals of the best fits to both bright and dim spectra, with completely independent parameters.

In the text

Fig. 22. Swift J1658.2-4242 light curve, with flux plotted on a logarithmic scale, analogous to Fig. 14 of C04. We also plotted exponentials describing the minimum type A flux, and the maximum type C flux, as a function of time. We only included NuSTAR, and XMM-Newton data, as their spectra enable an accurate flux determination, and as the QPO type can be unambiguously identified from their light curves. We removed data points during transitions between the two states to show the limiting fluxes of each state when there is no transition.

In the text