An underlying clock in the extreme flip-flop state transitions of the black hole transient Swift J1658.2-4242

Aims: Flip-flops are top-hat-like X-ray flux variations which occur in some transient accreting black hole binary systems and feature simultaneous changes in the spectral hardness and the Power Density Spectrum (PDS). They occur at a crucial time in the evolution of these systems, when the accretion disk emission starts to dominate over coronal emission. Flip-flops have only rarely been observed and are poorly understood. Methods: We detect 15 flip-flops in the 2018 outburst of Swift J1658.2-4242, in observations by XMM-Newton, NuSTAR, Astrosat, Swift, Insight-HXMT, INTEGRAL, and ATCA. We analyse their light curves, search for periodicities, compute their PDS, and fit their X-ray spectra, to investigate the source behaviour during flip-flop transitions, and how the interval featuring flip-flops differs from the rest of the outburst. Results: The flip-flops of Swift J1658.2-4242 are of an extreme variety, exhibiting flux differences of up to 77% within ~100s, much larger than has been seen so far. We observe radical changes in the PDS simultaneous with the sharp flux variations, featuring transitions between the Quasi-Periodic Oscillation types C and A, which have never been observed before. Changes to the PDS are delayed, but more rapid than changes in the light curve. Flip-flops occur in two intervals, separated by two weeks in which these phenomena were not seen. Transitions between the two flip-flop states occurred at random integer multiples of a fundamental period, of 2.761ks in the first interval, and 2.61ks in the second. Spectral analysis reveals the high and low flux flip-flop states to be similar, but distinct from intervals lacking flip-flops. A change in the inner temperature of the accretion disk is responsible for most of the flux difference in the flip-flops. We highlight the importance of correcting for the influence of the dust scattering halo on the X-ray spectra.


Spectral-timing states of a black hole transient
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 visualized 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 disk by high energy electrons in the corona (Gilfanov 2010). In the standard truncated disk model (Esin et al. 1997, andPoutanen et al. 1997), the geometrically thin, optically thick accretion disk is terminated at a large radius, so that the disk 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 disk starts to decrease. This causes the emission of the disk to become stronger, and the spectrum overall to be characterized by a combination 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). QPOs 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). 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 HFQ-POs, which commonly have frequencies of several hundreds of Hz. The LFQPOs have themselves been grouped into three different types: A, B, and C , 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: Here, ν 0 is the centroid frequency, at which the Lorentzian reaches its maximum. Q = ν 0 /FWHM is the quality factor, which is used to characterize the width of the Lorentzian. The temporal variability in the light curve is parameterized by the fractional root mean square variability (rms). r describes the variability of a single Lorentzian. If the PDS is normalized according to Belloni & Hasinger 1990, the rms is equal to the square root of the PDS integral: r = ∞ −∞ L(ν) dν. The total rms of a PDS is determined by taking the square root of its properly normalized integral between two frequencies: ν 1 , and ν 2 , and is typically expressed as a percentage. rms = ν 2 ν 1 n i L i (ν) dν, 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 ν 0 ≤ 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 . The type C frequency also roughly correlates with the X-ray flux, and anti-correlates with the broad-band rms .
Type B QPOs are the defining characteristic of the SIMS, distinguishing it from the HIMS. (Belloni & Motta 2016). They are typically located at: ν 0 ∼ 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 . 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 ν 0 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 considerable analysis of LFQPOs, there is no consensus as to their physical origin. QPOs are thought to 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, andMotta et al. 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 disk. 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, Ingram & van der Klis 2013, and Ingram et al. 2016 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, and Varniere & Vincent 2016 hypothesizes an instability with which energy and angular momentum are transported from a magnetized accretion disk 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 nonrelativistic 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 disk, 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 by Cabanac et al. 2010 considers a magneto-acoustic wave in the corona, which causes a variation in the efficiency of comptonization, 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
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.
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 flipflops 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-flops 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.

Swift J1658.2-4242
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, . 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, andLien 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 Article number, page 3 of 30 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 high 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.

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 utilized 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, Insight-HXMT LE, ME, and HE, 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 summarized 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. 42 observations by NICER were used 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 10 different times within the outburst. These observations are summarized in Table A.1

X-ray data
We produced Swift/XRT light curves and spectra using the University of Leicester's online data analysis software 1 , 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 2 (Krimm et al. 2013), which bins observations made throughout each day into one datapoint.
NuSTAR data were reduced using the standard procedure with nupipeline and nuproducts, using NuSTARDAS version 1.8.0, CALDB version 4.7.8 and the general functionality of 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 normalization 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. We therefore 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 −1 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 23 cm −2 . Jin et al. 2019 also detected a strong DSH around the black hole candidate. This provides a challenge for spectral fitting, as the DSH results effectively in 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, andJin et al. 2019 developed the spectral fitting model dscor to account for scattering of Xrays 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. In order 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 3 to calibrate the data from the Chandra High Energy Transmission Grating (HETG) observation of Swift J1658.2-4242, and extract an energy spectrum 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 Medium Energy X-ray telescope (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 −1 . By comparing the light curves of Insight-HXMT ME with those of NuSTAR, and XMM-Newton for the same en-ergy range, we however find the best possible agreement for an additional background count rate of 3.9 counts s −1 . 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 was chosen because different telescopes disagreed on the flux densities measured at lower energies, as a result of the DSH. 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 GHATS software 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 powerspectra for intervals of 2 10 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 are hampered by low statistics.
We subtracted the Poisson shot noise from all Leahynormalized PDSs (Leahy et al. 1983), and subsequently applied the square fractional rms normalization (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 (ν 0 , 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 photoionization 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 disk black body model (Xu et al. 2018).

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 UVMUL-TIFIT to search for intra-observational variability within the radio observations (to identify the possibility of radio flaring associated with the X-ray flip-flops).

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 detect Swift J1658.2-4242 in the LHS, with the typical flux increase from quiescence. We depict the total outburst in the HIDs of Fig. 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 −11 .
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 Fig. 2 and 3.  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 characterize 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 exactly 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 optimized 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 datapoint per observation. We restricted the energy range of Swift/XRT, NuSTAR, and XMM-Newton to 5.5-10 keV, to minimize 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.
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 3−79 = 9.75 × 10 −9 erg cm −2 s −1 . This corresponds to an unabsorbed and unscattered flux of F 3−79,0 ≈ 1.64 × 10 −8 erg cm −2 s −1 . Using the distance estimated by Jin et al. 2019 of ∼ 10 kpc, we find the unabsorbed and unscattered source luminosity to be L 3−79 ∼ 2 × 10 38 erg s −1 .
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 see similar flux increases or decreases, with the bright states having a 49-77% higher flux  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. than the dim states (see Fig. 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, andJithesh 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.
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 Fig. 4a, 4b, and 4c, we plotted the light curves and hardness ratios of individual NuSTAR and simultaneous XMM-Newton observations containing flip-flops. 4a also features the ATCA observation of a radio flare, and 4c is complemented by an As-   trosat 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 state, and anti-correlated with flux variations in the dim state. 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 −7 , 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 flip-flops started at the same time as the highest intermediate, and hard X-ray flux was reached, and the major radio flare was detected.
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 each lasted for 2.2-19.0 ks, or longer (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 Figure 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 Fig. 4a, 4b, and 4c. The duration of these short segments in which one, or more bright states were observed, increases 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 disk black body to power law emission.
The flux change during transitions decreased gradually over time, as can be seen by the dotted lines in Fig. 1, and 5. But there is only very little variation between the flux changes of adjacent transitions. This is illustrated by the blue line in Fig. 12, which does not show any discontinuities at the 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: to the region of the light curves from 100 s before the start, to 100 s after the end of each transition. y 0 is the flux just before the transition, and y 1 is the flux just after it. t 0 is the time when the transition starts, and the transition duration is ∆t. m 1 and m 2 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-170 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 (∼ 600s), 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, 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 in an exponential manner.
Both NuSTAR and XMM-Newton observations in this interval lie in almost exactly 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 −22 . 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. 7). 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.
Interestingly, the X-ray flare is placed in exactly 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 observe 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 −20 .
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 is 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 flux of about 27% during a gap in the data, followed by a gradual flux rise up to the previous level (see Fig. 8). The hardness ratio correlates with the flux for this observation (with r xy = 0.89, and p = 4.1 × 10 −20 ), 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 a flip-flop. In the first NuSTAR observation, two absorption dips are seen (Xu et al. 2018), which 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. 6. As far as we can tell, Swift J1658.2-4242 never reached the HSS, but remained in the intermediate states for the rest of the Swift/XRT observations. 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.

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 Article number, page 11 of 30 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: 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 ∈ Z. 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 Equation 1. Of the nine transitions we observed, eight are fitted well by Equation 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 Equation 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 χ 2 = 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 different period that could fit our observations better than this one, but did not find any other period larger than 2761.0 s, which fitted better than it. Shorter periods were found, but these turned out to be merely fractions of the 2761.0 s period. Here we plotted the 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 randomized 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 Equation 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 exactly the same 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 minimize the χ 2 , 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 flipflop 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 at least one instance where a better, or equally good fit was obtained as we had found for the actual data, using the same procedure.
When running 10 5 iterations of this simulation, we found that only 3.38% of randomly generated sets of transition times fitted with a χ 2 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 χ 2 of 18.2 for 7 degrees of freedom. However, we found that in 10 4 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 flipflops. 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 are 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 Equation 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 χ 2 = 42.5, for 4 degrees of freedom. These best fit parameters were used to plot the vertical lines in Fig. 7. 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 5 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 χ 2 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: 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 = 1 2 ΠΠ, 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: 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 flipflops 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 −2 s, resulting in an initial rate of period decrease equal toΠ(0) = (−6.28 ± 0.04) × 10 −5 . The best fit has 11 degrees of freedom, and a goodness of fit of χ 2 = 52.6. We indicate the difference between the measured transition times, and the closest value of Equation 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.
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 flipflops. We fitted these times with Equation 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 4 iterations, we found that the probability to obtain at least as good of an agreement in randomized 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 flipflops 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 Equation 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, 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 fit 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 −2 %, or 3.2σ. This result is based on 10 5 simulations for the early flip-flops, and 5 × 10 5 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 exactly 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 gaps contradict 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.

Power density spectra and QPOs
In order to learn more about the flip-flop transitions, we generated PDSs of individual flip-flop states, and of entire observations without flip-flop activity, from the light curves of Astrosat, NuS-TAR, 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. 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. Producing 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 −10 erg cm −2 s −1 to 1.1 × 10 −9 erg cm −2 s −1 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 ( Fig. 12 and 13), with r xy = 0.78 and p = 7.3 × 10 −19 . We additionally found that the quality factor increases slightly with increasing count rate (r xy = 0.38, p = 3.2×10 −4 ). And the rms in the QPO was found to decrease with increasing count rate (r xy = −0.51, p = 3.8 × 10 −7 ), thereby indicating that the rms also correlates negatively with the centroid frequency (r xy = −0.64, p = 2.0 × 10 −11 ).
By comparing these results with the distinguishing characteristics of the three different QPO types (Motta 2016), we find 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 a smaller quality factor 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 Fig. 12 and 14, no QPO can be identified during the bright states. 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). Fig. 12 also shows, that the broad-band continuum at low frequencies is strongly suppressed in bright states. Fig. 11 however indicates that the continuum grows at low frequencies. It 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, is consistent with the Astrosat result, of the bright states featuring type A QPOs.
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 flipflop 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 small 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 Fig. 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 power spectra 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 crystallize 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. 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.
We note that the change in the PDS is very sudden. In Fig. 14, we see that the type C QPO takes less than 20 s to appear. Using smaller time intervals indicates that it takes at most ∼ 10 s to get started. There was no indication before then, that a type C QPO was being established there. This is significantly shorter than the flux transition, which instead takes ≈ 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 −16 for the Astrosat observations, and r xy = −0.91, p = 1.0×10 −20 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 e.g. 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.
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 models 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 of 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 Article number, page 15 of 30 Fig. 13: Correlation between the best QPO fitting parameters and the count rate, in each of the 200 s intervals of Figure 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, as described in the section on Quasi-periodic oscillations. Different ranges of values for each of these parameters were found when binning over longer time periods. ∼ 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 PSD 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 a 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. 8) 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  Fig. 12. Colours denote the Leahy normalized, 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. 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.

Dust Scattering Halo and Energy Shifts
In Fig. 16 we plot the NuSTAR and XMM-Newton spectra of an absolutely 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 exactly 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 Fig.  1, 4, 5, and 6 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. 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 andJin et al. 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 andXu et al. 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 disk wind (Ponti et al. 2012, Ponti et al. 2016, and Díaz Trigo & Boirin 2016. The simultaneous XMM-Newton spectra however feature an emission, rather than an absorption line at this energy. Clearly 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 for 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.099 −0.056 keV, a variance of (1.47 ± 0.85) × 10 −1 keV, and a normalization of −2.23 +0.75 −0.53 ×10 −4 photons cm −2 s −1 , with the negative sign indicating that this is an absorption line. This line has an associated equivalent width of −28.5 +9.6 −6.8 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 normalization to be: 1.9 +3.2 −2.9 × 10 −4 photons cm −2 s −1 , corresponding to an equivalent width of 2.4 +4.0 −3.7 × 10 −2 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 find 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.
(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. 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 Figure 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.

Spectral Analysis
Comparing the spectra of the bright, dim, and late flip-flop states, we find that they are remarkably similar (see Fig. 18). Visually, besides having a different normalization, 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. This results in the increase in hardness ratio simultaneous to flux increases, as seen in Fig. 4.
The initial LHS spectrum differs noticeably from the other spectra we measured (see Fig. 18). We also find 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 normalization. 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. Fig. 17: Graph of the 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.
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, using two sets of XSPEC models.
In Model 1, we fit the spectral energy contribution from the multicolour black body component of the accretion disk, 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 individual 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 disk 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 disk 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 disk radii using this model, so we set these parameters equal to 6 GM/c 2 , and 1000 GM/c 2 , 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 disk temperature (T in ), the emissivity power law index of the iron line (B 10 ), the power law photon index (Γ), the high energy cutoff (E cut ), as well as the normalization of the three additive spectral components, diskbb, diskline, and cutoffpl.
We also fitted the energy spectra using a more physical comptonization model by Poutanen & Svensson 1996, compPS. As this model describes both the disk 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 disk radii, and the inclination of the accretion disk, which were already used in Model 1. We used a spherical corona geometry and specified a multicolour disk black body component by confining the disk 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 disk. We fit the spectra with a free hydrogen column density (N H ), emissivity power law index (B 10 ), coronal electron temperature (T e ), inner accretion disk temperature (T in ), coronal optical depth (τ y ), relativistic reflection normalization (R r ), and normalization of the compPS, and diskline components. Spectra fitted with both Models 1 and 2 always had a lower χ 2 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 (Fig. 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.
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 disk 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, with all parameters described above being set free to vary, 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 disk temperature, and the effect this has on the power law component.
In plot (c), we also change the compPS normalization 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 Fig. 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 Fig. 19, and 20, we find r xy,1 = 0.93, p 1 = 1.0 × 10 −8 , and r xy,2 = 0.83, p 2 = 1.2 × 10 −5 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 χ 2 which was smaller by 693. We therefore conclude that the observed correlation between hydrogen absorption 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 Fig. 21 (c), and (d). The low p-value for the null hypothesis which assumes a constant N H in both states, of p = 6.4 × 10 −4 , and p = 1.9×10 −4 for the two models, indicates that the hydrogen column density changes during the flip-flops. Both Xu et al. 2019 andJithesh 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 test of whether the N H changes during the flip-flops, we repeated the above test, but this time only for the flip-flop spectra. We found that the decrease of 9 degrees of freedom was accompanied by a drop of the total χ 2 , by 103. This also supports the notion that the N H changes during flip-flop transitions.
In both models we find an increase in the inner disk 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 = 1.6 × 10 −4 , and r xy,2 = 0.91, p 2 = 2.4 × 10 −4 for the two models. Interestingly, the two states by themselves do not show any significant correlation at all, with p 1 = 0.51, p 2 = 0.81 for the bright states, and p 1 = 0.98, p 2 = 0.59 for the dim states. At comparable fluxes, the disk 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 normalization on X-ray flux, we plotted a straight 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 normalization 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 normalization graph indicate the average value 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 +220 −50 keV. We also excluded one bright state, and three late flip-flop data points from the graph of Betor 10 , as the fits were insensitive to this parameter, and the errors in the measurements could not be determined.
line depicting a constant equivalent width, under the assumption of a uniform spectral shape. Fig. 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 disk 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 flipflops (20 ± 5 eV), and the final NuSTAR observation in the SIMS We excluded three late flip-flop data points from the graph of Betor 10 , because the fits were insensitive to this parameter, and the errors in the measurements could not be determined.
(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 10 , does not change significantly throughout the range of observations we analysed, in either of the two models.
The normalization of the diskbb component, is given by the equation: N dbb = R in D 10 2 cos i. It is proportional to the square of the accretion disk inner truncation radius (R in ), and the cosine of the disk inclination (i). The inclination of the accretion disk is not known, and the distance to the binary (D 10 , 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 normalization 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- 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 disk 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 normalization 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 N H (d). Finally, in plot (e) we show the residuals of the best fits to both bright and dim spectra, with completely independent parameters. 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 1 = 0.18. In contrast, the observations without flip-flops differ too much to be considered as having the same normalization. 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 1 = 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 normalization of the cutoffpl component increases almost linearly for all observations, with no clear distinction between the flip-flop and non flip-flop states, or the nature of the QPO. When using the results in all spectra, we obtain a correlation coefficient of r xy = 0.94, and a p-value of p = 3.6 × 10 −9 .
We notice a near 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 see a similar trend in the electron temperature as a function of flux, of the compps model, as in the high energy cutoff of the cutoffpl model. This was expected, as these two parameters are related. The flip-flops states have T e ≈ 80.2 ± 4.5 k B −1 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 ≈ 68 ± 11 k B −1 keV. Observations of the late flip-flops, and the SIMS have a higher average electron temperature of T e ≈ 140 ± 14 k B −1 keV. Finally, the relativistic reflection component of the compps model does not indicate a dependence on flux for the early flipflop states, which have an average of R r ≈ 0.30 ± 0.15. This also supports the idea that the area of the disk, as seen from the central X-ray source, remains constant during flip-flop transitions. Between the early and late flip-flop intervals, we detect 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 disk temperature, and the power law normalization. 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 delayed change in N H could not be corroborated.
The final NuSTAR observation still had both a strong multicolour disk black body component, in addition to a strong power law 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 disk black body component dominates. The remaining observations by Swift/XRT needed to be fit with both a black body, and a power law component. Therefore, it seems like Swift J1658.2-4242 never reached the HSS.

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 2 of which were inferred. This is more than have 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. It is possible that part of the apparent scarcity of flip-flops in the literature is 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 , the dim flipflop 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. This AS has similar timing properties as the SIMS, is however 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 subsample with long duration states (of order ∼ 0.1-10 ks), and long transition times (of order ∼ 10-100 s). We did not observe any of the short duration flipflops, 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 1 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 involve 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 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 fits 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 find 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 find a clear gap between the two. A noticeable difference is however that we detected different decay rates of these two exponentials, whereas C04 found parallel lines in this log-linear light curve. 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 datapoints during transitions between the two states to show the limiting fluxes of each states 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 longlived 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. It therefore 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 separate 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, that types A and C originate from the same physical process, but type Bs have an entirely different cause.
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 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 comptonized 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 flipflop 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 flipflop 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 flipflops, 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 normalization implies that the flipflops have very similar values of R 2 in cos i, 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 exactly 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. QPOs were not detected in any of the late flip-flop states, but a noticeable difference in the rms was nonetheless observed. 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 alternatively 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 mecha-nism at play during it. These different properties also suggest that this interval really did not feature any unobserved flip-flops. However, of the two NuSTAR and XMM observations in this interval, the latter observations shows 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 Fig. 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 this 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, is seen in the fitted truncation radius of the accretion disk. Fig. 19 points to a consistent value of the diskbb normalization 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 normalization, but these are all a factor of 2-4 smaller than what is found within the flip-flop intervals. The diskbb normalization is defined as: N dbb = R in D 10 2 cos i, where R in is the inner radius of the accretion disk, D 10 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 normalizations is due to a difference in the truncation radii. This result would therefore indicate that the disk is truncated at a very specific radius whenever flip-flops occur. Along the same reasoning, other times of the outburst, which do not feature flip-flops, have no consistency in their truncation radii.
The fitted diskbb normalization 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 disk 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 normalizations to correspond to accretion disk 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 normalization of the disk 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 disk. 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 disk 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 disk temperature, at a constant inner radius. An additional increase of the power law normalization, 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 on 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 initiates an instability that pushes 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 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 disk, 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 disk.
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 such an 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 ionized, 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 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 comptonization of disk and corona emission whenever the outflow points into the line of sight. The bulk comptonization 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 PSD 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 of time. 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 polarization 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.

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 simultaneous 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 the two 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 disk temperature. The remaining disparity was caused by an additional increase of the power law normalization and hydrogen column density. We found that the inner accretion disk 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.