The Chaotic Long-term X-ray Variability of 4U 1705--44

The low-mass X-ray binary 4U1705-44 exhibits dramatic long-term X-ray time variability with a timescale of several hundred days. The All-Sky Monitor (ASM) aboard the Rossi X-ray Timing Explorer (RXTE) and the Japanese Monitor of All-sky X-ray Image (MAXI) aboard the International Space Station together have continuously observed the source from December 1995 through May 2014. The combined ASM-MAXI data provide a continuous time series over fifty times the length of the timescale of interest. Topological analysis can help us identify 'fingerprints' in the phase-space of a system unique to its equations of motion. The Birman-Williams theorem postulates that if such fingerprints are the same between two systems, then their equations of motion must be closely related. The phase-space embedding of the source light curve shows a strong resemblance to the double-welled nonlinear Duffing oscillator. We explore a range of parameters for which the Duffing oscillator closely mirrors the time evolution of 4U1705-44. We extract low period, unstable periodic orbits from the 4U1705-44 and Duffing time series and compare their topological information. The Duffing and 4U1705-44 topological properties are identical, providing strong evidence that they share the same underlying template. This suggests that we can look to the Duffing equation to help guide the development of a physical model to describe the long-term X-ray variability of this and other similarly behaved X-ray binary systems.


INTRODUCTION
X-ray binaries exhibit periodicities on multiple timescales, which provide information about the physical mechanisms at play.High magnetic fields and disc-magnetosphere interactions create pulsations on timescales from milliseconds to seconds (van der Klis 2006).Orbital modulations are seen from minutes to tens of days (Levine et al. 2011).A number of X-ray binaries show evidence of long-term periodicities, or super-orbital periods, on timescales much longer than their orbital periods (Clarkson et al. 2003a, Clarkson et al. 2003b).Some X-ray binaries have variability on timescales of over a hundred days that are not strictly periodic (Boyd et al. 2000).
For low-mass X-ray binaries (LMXBs), it is generally accepted that the accretion is primarily due to a substantial accretion disc, which is unstable to irradiation-driven warping (Foulkes et al. 2010).The precession of the accretion disc is potentially the mechanism underlying the observed long periods.Warped or twisted E-mail: rebecca.a.phillipson@drexel.edu† E-mail: patricia.t.boyd@nasa.govaccretion discs in general are also invoked to explain phenomena observed across many types of systems including high-mass Xray binaries, cataclysmic variables, proto-planetary discs and active galactic nuclei (Ogilvie & Dubus 2001).One relevant example is SS433, which has a measured 160-day precession of the relativistic jets identified to be associated with precession of the accretion disc about the central neutron star (Whitmire & Matese 1980;Collins & Scher 2002).Other well-known examples in X-ray binary systems include Her X-1, with a varying long-term period of about 35 days associated with its observed anomalous low states (Still & Boyd 2004), and SMC X-1, with confirmed 60-day variations in the X-ray (Wojdowski et al. 1998).
In this study, we consider the low-mass X-ray binary, 4U1705-44, which exhibits the high-amplitude transitions and non-periodic long-term variability of interest ( Hasinger & van der Klis 1989;Muno et al. 2002).4U1705-44 is a neutron star with mass 1.1 − 1.6 M (Olive et al. 2003) at a distance of 7.4 kpc (D'Aì et al. 2010).It is considered to be part of the class of 'atoll' sources (Wang et al. 2012), however it has sampled the upper parts of a Z in the X-ray colour-colour diagram during its soft-to-hard state transitions (Barret & Olive 2002).An infrared counterpart has been observed suggesting a dwarf star companion with a 1-10hr orbital period (Homan et al. 2009), corresponding to a mass less than 0.5 solar masses.
The RXTE (Rossi X-ray Timing Explorer) All Sky Monitor obtained approximately 14 years of daily monitoring in the 2-12 keV energy range of 4U1705-44 with its scan of 80 per cent of the sky every 90 minutes (Bradt et al. 1993).MAXI (Monitor of All-Sky Xray Image) on-board the ISS, in operation since August 2009, scans the sky every 96 minutes in the 2-20 keV band and continues to provide daily monitoring on 4U1705-44 (Matsuoka et al. 2009).With MAXI and RXTE combined, we have a nearly 20 year continuous light curve revealing over 50 cycles of the high amplitude transitions that occur on timescales on the order of hundreds of days.It is only recently that such continuous, high signal-to-noise, uninterrupted time series have been available at the timescales of interest for investigating long term variability in X-ray binaries on the order of years.The long-term modulations, which have also been observed in some other LMXBs, are not strictly periodic (Boyd & Smale 2004).In this paper, we present evidence that the long-term variability is reminiscent of a nonlinear double-welled oscillator evolving in the chaotic regime.We find that the Duffing oscillator is a good candidate to describe the system.Since the parameter ranges that best represent the 4U1705-44 light curve are indeed in the chaotic regime of the Duffing oscillator, we use analyses appropriate for non-linear and chaotic time series.
There are two broad approaches to understanding chaotic behaviour in dynamical systems.The metric approach generally involves computing the Lyapunov exponents and various dimensions, such as the Correlation and Minkowski dimensions, giving a tight range for the fractal dimension (Ott 2002;Anishchenko et al. 2007).However, these methods require very large datasets and degrade rapidly with noise and are therefore rarely applicable to astronomical data (Mindlin et al. 1990).The topological approach involves the identification of the two mechanisms that are responsible for the creation of a strange attractor.These two mechanisms are the stretching and squeezing mechanisms correlating to the sensitive dependence on initial conditions and recurrence phenomenon, respectively (Mindlin & Gilmore 1992).The stretching and squeezing mechanisms can be characterised by computing how the unstable periodic orbits are uniquely organized topologically (Letellier et al. 2006).In fact, it is possible to determine how the unstable periodic orbits embedded in the attractor are organized in terms of a set of integers (Solari & Gilmore 1988a, Tufillaro et al. 1990).Extracting these integers from a time series is robust against noise and can be performed for smaller datasets (Gilmore & Lefranc 2002).
Thus, we will follow the topological approach.Two of the topological invariants are called the linking numbers, the unique set of integers aforementioned, and the relative rotation rates (RRR), a set of fractions which when summed over all initial conditions equal the linking numbers (Gilmore 1998).It is straightforward and sufficient to compute the RRR (Gilmore 1998), which we will do for both 4U1705-44 and the Duffing oscillator (Solari & Gilmore 1988b).By extension of the Birman-Williams theorem (Birman & Williams 1983a;Birman & Williams 1983b), two chaotic attractors are equivalent if they are described by branched manifolds that can be smoothly deformed, one into the other.In other words, the RRR and linking numbers will be invariant under transformations and control-parameter changes.Conversely, the 4U1705-44 and Duffing oscillator time series share the same underlying template if their RRR are identical.Identifying the underlying template of 4U1705-44 via the RRR allows for the prediction of the RRR of all other possible orbits in the time series and provides a qualitative model for the flow that uniquely generates the chaotic times series (Mindlin et al. 1990;Letellier et al. 2006).If both the 4U1705-44 and Duffing oscillator time series share the same underlying template, then we can look to the Duffing oscillator equation and its template to make qualitative predictions of the behaviour of 4U1705-44 and infer possible physical parameters of the system (Gilmore & Lefranc 2002).
To follow the topological analysis procedure we first must determine the lowest order period in the data in question.We can then use the close returns method as used in Boyd et al. (1994) to locate and extract regions in the time series that can be used as surrogates for the unstable periodic orbits.The topological invariants of RRR of all pairs of unstable periodic orbits extracted can then be computed from the time series and compiled into a matrix or table, as exemplified by Solari & Gilmore (1988b) for the Duffing oscillator.This procedure was laid out generically for the computation of RRR by Mindlin & Gilmore (1992) and Gilmore (1998).
This paper is organized as follows.In Section 2 we will review the data obtained from 4U1705-44 using the RXTE ASM and MAXI and our data reduction techniques, as well as the generation of numerical times series from the Duffing oscillator equations of motion.In Section 3 we will determine the lowest order period of an unstable periodic orbit for both the 4U1705-44 and Duffing data using the close returns method and via dynamical power spectra and time-delay embedding analyses.In Section 4 we will use the lowest order period and close returns method to locate and extract regions in the time series that can be used as surrogates for the unstable periodic orbits.In Section 5 we will compute the topological invariants of RRR of all pairs of periodic orbits extracted from the time series and compile these RRR into matrices for the 4U1705-44 and Duffing data.Finally, in Section 6, we will compare the resulting intertwining matrices containing the RRR followed by a discussion and concluding remarks.

4U 1705-44
The Rossi X-ray Timing Explorer (RXTE) satellite was launched December 30, 1995 and was decommissioned January 5, 2012 (Bradt et al. 1993).The RXTE satellite embodied three instruments to make observations in the X-ray bandwidth (2-200 keV) with a large effective collecting area ( 0.8 m 2 ).One of the instruments was the All-Sky-Monitor (ASM), which scanned 80 per cent of the sky every 90 minutes in the 2-12 keV energy band, with a timing resolution of 90 minutes or more and sensitivity of 30 mCrab (Levine et al. 1996).The low-mass X-ray binary system 4U1705-44 was one of the 75 sources continuously observed by RXTE's ASM, and was immediately noted by Levine et al. (1996) to have highamplitude variability ranging from <25 mCrab up to 300 mCrab.
The Monitor of All-sky X-ray Image (MAXI) telescope is aboard the International Space Station (ISS) orbiting Earth as part of the Japanese Experiment Module with the goal of obtaining a map of the X-ray sky with greater sensitivity than previous allsky monitors.MAXI has been in operation since August 2009 and scans nearly the entire sky every 96 minutes in the 2-20 keV energy band and with a sensitivity of 20 mCrab (Matsuoka et al. 2009).MAXI obtained observations of 4U1705-44 since its operation start.We therefore have continuous, daily observations of 4U1705-44 and other bright X-ray sources from December 1995 through the present, sampling more than two decades.For this in-vestigation, we consider the light curve of 4U1705-44 from 6 January 1996 through 16 May 2014.Furthermore, because MAXI and RXTE ASM were in operation simultaneously for 3 years, the light curves of 4U1705-44 can be combined in a seamless manner.
The construction of the 4U1705-44 time series consisted of combining RXTE ASM and MAXI data.For both datasets, we rebinned the data to every 4 days, since our primary timescale of interest was over hundreds of days, thereby improving our signalto-noise.The online tool1 based on Kirsch et al. (2005) converts from a flux expressed in Crab flux units in a specific energy band to physical units of erg cm −2 s −1 for ease of reporting the flux from various instruments in terms of the Crab nebula.Given the respective average count rate from the Crab nebula observed by each instrument over the same observation time as 4U1705-44 (75.3 counts/second for RXTE ASM and 3.25 counts/second for MAXI) and the corresponding physical units from both instruments in each bandwidth (1.5 -12 keV and 2 -20 keV) corresponding to 1 Crab, we can therefore report the light curve in terms of the physical units taking into consideration the differing energy bands of RXTE ASM and MAXI.The same method of normalising to Crab units to form an extended light curve in physical units was used for the multiwavelength, multi-instrument analysis of the X-ray binary, LMC X-3, by Torpin et al. (2017).
Our methods will also require an evenly sampled dataset.Following the same procedure as Smale & Boyd (2012), any gaps or values that contained errors greater than 5σ were replaced by a random value containing the same mean and standard deviation as the entire dataset.This treatment was applied to 4.8 per cent of the data.
The final year of observations from the RXTE ASM had consistently higher noise, which was seen in many light curves from various sources as also noted by Smale & Boyd (2012).In order to avoid potential contamination of the time series, we determined a cut-off point for the ASM data.As both the Crab Nebula and Cassiopeia A are bright sources with a consistently high signal-tonoise, we chose the cut-off point to be where the errors in the daily measurements of these two sources consistently exceeded the 3σ range in the same daily monitoring of 4U1705-44, which occurred in October of 2010, in agreement with the cut-off date determined by Smale & Boyd (2012).
MAXI started its observations of 4U1705-44 in 2009, before the October 2010 ASM cut-off date, and we could therefore use the overlap to appropriately scale the MAXI data to the ASM data.That is, starting with the two light curves displayed in Fig. 1, we used the light curve from RXTE ASM as a reference and scaled the entire light curve from MAXI such that the amplitude range (maximum and minimum flux difference) in the region in which both light curves overlapped were equal.Finally, we normalised the entire combined dataset by the maximum amplitude of flux present in the light curve before applying all methods described below.
The resulting time series, showing both the raw data with errors and a low-pass filtered version (in order to smooth the highfrequency oscillations) is plotted in Fig. 2. It is clear by visible inspection of Fig. 1 and Fig. 2 that, while the light curve displays large amplitude oscillations where the flux increases by more than a factor of ten from minimum to maximum, these oscillations are far from strictly periodic.In addition to the large amplitude oscillations, the flux from 4U1705-44 displays smaller amplitude oscillations around the low flux level and the high flux level.In general, the longer the system stays near the average high(low) flux level, the more smaller amplitude oscillations are seen.This is reminiscent of the one-dimensional time series from the classic Duffing equation, a canonical nonlinear oscillator shown in Fig. 3.This similarity motivates us to more quantitatively compare the 4U 1705-44 light curve to the Duffing oscillator.

Chaos and Topology
2.2.1 Is the Time Series of 4U 1705-44 Chaotic?
The low-mass X-ray binary, 4U 1705-44, has previously been found to have high-amplitude transitions not easily fitted by a simply periodic model but with a power spectrum that displays clusters of multiple periodicities (Durant et al. 2010), characteristics often exhibited by both chaotic and stochastic nonlinear systems.Signatures of deterministic chaos have also been detected as the driver of observed variability in five black hole X-ray binaries (Suková et al. 2016) on short timescales relating to quasi-periodic flares and oscillations.Unlike stochastic systems, low-dimensional chaotic systems display characteristic behaviour when cast into a phase space (for example, by plotting the flux-like variable and its first time derivative; Takens (1981); Ott (2002)).The phase space embedding of the Duffing oscillator time series (Fig. 3) is shown in Fig. 4. The two centres in phase space about which the trajectory evolves are the hallmark of this 'double-well potential' (Holmes & Rand 1980).We discuss the Duffing oscillator in more detail in Sec.2.3.
The phase space trajectory (the flux versus its first time derivative) of 4U1705-44 is shown in Fig. 5. Upon initial visual inspection of the phase-space embedding (Fig. 5), we see the first indication of a chaotic structure, i.e. an axis about which the flow evolves, or a 'hole' in the middle of the projection (Letellier et al. 2006), which indicates the presence of (unstable) periodic orbits.The 4U1705-44 phase space embedding also shows evidence of two centres of rotation suggestive of a double-well potential.
Distinguishing between chaotic or stochastic driven time series is an ongoing endeavour in time series analysis in general and particularly challenging in the astrophysical context due to the often short and unevenly sampled data.Chaotic time series arising from highly nonlinear systems produce similar trends to time series arising from stochastic processes, namely, the long-term unpredictable behaviour and broad-band power spectra (Osborne & Provenzale 1989).Second-order metrics, such as the power spectrum or autocorrelation function, do not contain enough information to distinguish between the underlying processes; a nonlinear stochastic process and a deterministic chaotic process can have the same statistics (Zunino et al. 2012).Furthermore, non-stationarity in a time series can also result in falsely identifying nonlinear structures (Schreiber & Schmitz 2000).A new trend of utilising higherorder metrics, such as time reversibility and third-order autocorrelation functions, has therefore developed in order to explore the presence of nonlinearity and determinism (Collis et al. 1998;Gao et al. 2007;Zunino et al. 2012).
For data driven analysis, the application of nonlinear time series methods should be justified by establishing nonlinearity in the time series.Using the method of surrogate data (Theiler et al. 1992), we can determine to a significance level of 95 per cent that nonlinearity is present in the 4U1705-44 time series by testing for time reversal asymmetry in the 4U1705-44 light curve (Diks et al. 1995;Schreiber & Schmitz 2000).Furthermore, the local vs.global linear prediction test by Casdagli (1992) and recurrence analysis (Casdagli 1997;Marwan et al. (2007)) also indicate that the 4U1705-44 time series is likely deterministic and chaotic.See  the Appendix for a full description of the surrogate data method, nonlinearity tests, and recurrence analysis.We will therefore use nonlinear analysis and chaos topology methods to characterise the variability of the 4U1705-44 light curve.

Topology Methods
At the center of chaos theory are two fundamental properties of an attractor.First is the 'stretching' mechanism, or the sensitive dependence on initial conditions.This can be mathematically characterised by the exponential divergence between nearby trajectories in the phase space of an attractor, where the exponent of the diverging term is called the Lyapunov exponent and must be positive.The 'squeezing' mechanism relates to the characteristic that the trajectories in the phase space are bounded, corresponding to negative Lyapunov exponents.The constant stretching and squeezing of trajectories constitutes the fractal behaviour of a chaotic attractor (Gao et al. 2007, Ott 2002).
There are various ways to visualise an attractor beyond determining the Lyapunov exponents or fractal dimension.Poincaré developed topology (the compilation of his original papers from 1892 through 1904, translated into English, are available in Poincaré (2010)) in order to study differential equations geometrically.This resulted in Poincaré sections, bifurcations, and, with the entrance of Lorenz, more complex return maps, all of which are related to the phase space of an attractor (Lorenz 1963;Ott 2002;Anishchenko et al. 2007;Poincaré 2010).In the past thirty years, topological methods from Lie group theory have also been used to study lowdimensional chaotic attractors (Gilmore 1998;Gilmore & Lefranc 2002;Letellier et al. 2006).Poincaré identified chaotic attractors as recurrent but not periodic (the Poincaré recurrence theorem).That is, a chaotic trajectory will return to the neighbourhood of an initial position given enough time.Trajectories that are then perturbed just enough to close on themselves are periodic orbits, but are not stable.The mechanism(s) that creates the attractor simultaneously creates and organizes all the unstable periodic orbits in it.The advantage of the topological methods drawn from group theory is in their ability to then apply an integer invariant associated with each pair of closed, unstable orbits, called a linking number.Topological analysis using templates and bounding tori (Birman & Williams 1983a;Birman & Williams 1983b;Mindlin et al. 1990;Mindlin et al. 1991;Mindlin & Gilmore 1992) extends the geometrical methods of Poincaré to characterise and classify chaotic attractors embedded in three-dimensional spaces (Letellier et al. 2006).The first branched manifold for describing a chaotic attractor was drawn in 1963 by Lorenz using 'isopleths' (Lorenz 1963).
The connection between the embedded phase space of a chaotic attractor and a template (branched manifold, or 'knotholder') depends on the Birman-Williams theorem, which states that topological invariants of the periodic orbits are the same in the attractor as in its so-called 'caricature', the two-dimensional branched manifold (Birman & Williams 1983a;Birman & Williams 1983b;Gilmore 1998).The transformation from a chaotic attractor embedded in a three-dimensional phase space to a branched manifold or template and back preserves the topological invariants (e.g.linking numbers, relative rotation rates).Furthermore, these invariants are unique to specific branched manifolds.For example, the canonical Lorenz attractor corresponds to a variation of the Smale horseshoe template (Rössler 1977) whilst the Rössler system, considered the simplest chaotic attractor from the topological point of view, corresponds to a simple stretched and folded ribbon (Rössler 1976).
The crucial point relevant to our analysis is the fact that these two chaotic attractors are topologically inequivalent because their underlying templates are different and thus their unstable periodic orbits are organized differently, as described by the difference in their linking numbers and relative rotation rates.In fact, the three most widely cited examples of low-dimensional dynamical systems exhibiting chaotic behaviour, the Lorenz, Rössler, and van der Pol-Shaw attractors, are all associated with different branched mani-folds, and are therefore intrinsically inequivalent (Gilmore 1998;Gilmore 2007).
The Birman-Williams theorem becomes especially powerful in the comparison of topological invariants for periodic orbits extracted from data with the invariants of corresponding orbits in a branched manifold.The equivalence of two underlying topological templates does not necessarily dictate that they look identical when projected onto two dimensions, such as in a time series or embedded phase space.Rather, by identifying the underlying template through the determination of its topological invariants, one can make a mathematically rigorous statement about whether two systems are related, or not (Solari & Gilmore 1988a;Mindlin et al. 1990;Tufillaro et al. 1990).
In this paper, we seek to extract and compare the topological invariants of two different time series, one from the light curve of 4U1705-44 and the other from the numerical solution of the Duffing oscillator.By the Birman-Williams theorem, the equivalence of the invariants (or lack thereof) will directly determine the relationship between the two underlying templates that organize the unstable periodic orbits on the attractor.If we find that the templates are identical (or, more accurately, not inequivalent) then we can make the case that the equations of motion producing the two different time series are related.

The Duffing Oscillator
The phase space trajectory of 4U1705-44 displays a double-welled behaviour where one well is more favoured and tightly wound than the other.Such behaviour is shared by some damped and driven harmonic oscillators.We therefore chose to explore various nonlinear oscillators that could display double-welled behaviour similar to 4U1705-44 and compare their characteristics.This led us to study the Duffing oscillator, which within certain parameter regimes remarkably resembles the light curve and phase space trajectory of 4U1705-44.With the goal of performing a parallel topological analysis on the 4U1705-44 data and generated solutions of the Duffing equation in mind, we systematically solved the Duffing equation with randomly varying parameter ranges in order to produce many solutions for comparison.
A generic form involving five parameters of the Duffing equation was used: (1) Solutions were generated using a 4th-order Runge-Kutta method within randomly sampled parameter ranges for which chaotic behaviour arises.Those solutions with qualitatively similar frequency of low-order, low-high amplitude transitions and doublewelled phase space features whereby one well was favoured over the other predominated an increasingly narrow range in parameter space.The size of the damping is controlled by δ, the non-linearity of the restoring force by β, the stiffness of the oscillator by α, the amplitude of the forcing by γ and the driving frequency by ω.Of particular interest, the ω driving parameter frequented the range between 3.91 and 5.15, corresponding to driving periods of 122 to 160 days.The Duffing solution plotted in phase space in Fig. 4 corresponds to the parameter values (β = 4.99, α = 8.18, δ = 0.31, γ = −6.81,and ω = 4.01), with a driving period of 140 days.These representative parameter values allow for a double-welled oscillator with moderate chaotic behaviour and is the solution we used as the reference Duffing solution for our topological analysis in this paper.(For a comparison, the parameter values often used to study the chaos of the Duffing oscillator are (β = 1, α = 1,  Note that our goal is not to generate a numerical time series from the Duffing oscillator that exactly reproduces the X-ray binary light curve; since sensitive dependence on initial conditions is a hallmark of chaos, this would not be possible.Rather, our goal is to investigate whether the two time series share the same topological invariants, as would be expected if they both are low-dimensional chaotic systems evolving according to the same family of equations of motion.To do this comparison requires extracting unstable pe-riodic orbits from both the real and the numerical time series, in order to identify the underlying templates, and compare them to each other topologically.For this, we chose a segment of the simulated time series that qualitatively reproduces the same behaviour in phase space, sufficient to perform this analysis.

FINDING THE LOWEST ORDER PERIOD OF 4U1705-44
A chaotic attractor is characterised by its periodic orbits, both stable and unstable.A topological analysis of the data can be performed once these periodic orbits are determined.Since the lowest-order driving period for 4U1705-44 is unknown and we would like to compare its value to the range of driving periods in the Duffing solutions, we use the method of close returns as utilised by Mindlin & Gilmore (1992) and Boyd et al. (1994), ideal for small, chaotic datasets to identify nearly closed orbits near unstable periodic orbits of low period.In Fig. 6 a blue point is plotted at (i, p) when, for that time i, there is a point p days later with flux that is 'close' to the flux value at time i.We define 'close' as within 10 per cent of the maximum amplitude of the time series, represented by the intensity of the colour map, .The maximum amplitude is the difference between the maximum and minimum flux present within the time series.Regions of the time series that are close to repeating themselves appear as horizontal lines in blue.These correlate to sections in the light curve in which the time series comes close enough to an unstable periodic orbit such that it remains close for at least one period and is therefore within a region close to a periodic orbit.
One example of a horizontal stretch is highlighted in Fig. 6, corresponding to the regions in the light curve outlined in blue (portions of the light curve at position i) and red (the light curve at position p) in Fig. 14.We locate and extract all such regions from the light curve and find the lowest order orbits have periods between 110 and 200 days.For comparison, we computed the close returns map for the Duffing solution time series in Fig. 7.The Duffing close returns map more distinctly shows regions of the time series that repeat itself, appearing as horizontal lines throughout, and without the noise that is present in the 4U1705-44 time series.We can constrain the lowest order period obtained from the close returns analysis with the dynamic power spectrum.Fig. 8 shows the total power spectrum of the light curve of 4U1705-44 plotted against the period in the left pane.In order to determine the significance of the spectral peaks, we follow the Shimshoni (1971) version of the Fisher test of significance in harmonic analysis (Fisher 1929).The spectrum is normalised to the total power minus the first power term in the Fourier series.A vertical, dotted line is plotted in the total power spectrum in the left pane at a threshold value of g = 0.00646 (in the notation of Shimshoni (1971) and Fisher (1929)) above which all peaks are considered significant at a 99 per cent confidence level or above.
Immediately, we observe all of the significant power in the spectrum is at periods longer than 100 days, corresponding to low frequencies, indicative of longer-term variability.The right pane is the dynamic power spectrum showing the power spectrum of a window size of 2048 days, centred on the middle day in the window as it evolves through the time series.That is, for a particular time in the time series, the power spectrum is plotted with the window centred at that particular time and the power at each frequency represented by the colour map.The ends are padded with random noise, out to half a window size, with the same standard deviation as the entire dataset in order to extract the periodicities in the ends of the light curve.
The lowest order period varies between 120 and 150 days, particularly in the first half of the time series, in agreement with the periods found in the close returns analysis.If we consider the total power spectrum in the left pane, we observe a cluster of peaks isolated right around 125 days.We also potentially see the appearance of modes of higher order periods, which could relate to period-2 and period-3 orbits in the data (twice and three times the lowestorder period).For a first-order period of 120 days, we would expect a period-3 to be at 360 days and we do in fact observe the strongest peak in the total power spectrum to occur near this value.The very long period near the end of the light curve is likely related to 4U1705-44 remaining in a high-state for hundreds of days, a deviation from its previous behaviour, in which longer wavelength Fourier modes on the order of the size of the window could be favoured.For comparison, we have also computed the dynamic power spectrum for the reference Duffing solution in Fig. 9.Most notably, we observe similar modal behaviour at similar frequencies to the 4U1705-44 power spectrum.
Finally, we can consider the periodic information we obtain from a time-delay embedding of the 4U1705-44 time series.A strange attractor is embedded in a 3-dimensional space and topological and geometrical analyses in chaos theory are performed on the embedded n-dimensional phase space of the attractor.One can reconstruct a strange attractor from scalar data, x(t), like a light curve, by constructing vectors y(t) with n components determined by the scalar x(t).One such embedding is determined as the differential phase space projection, as in Fig. 5, where y(t) is merely the derivative of x(t).Another is a time-delay embedding.This involves creating an n vector by the map where τ k are called the time delays.The time delays are typically evenly spaced multiples of the time delay: τ k = (k − 1)τ .According to Takens (1981), we can reconstruct a mapping of the strange attractor from the time series via Eq. 2 by merely introducing a time-related shift.By adjusting the time delay component, we can determine the value at which it approaches the differential phase space embedding, which is related to the lowest order period of the time series.A simple example of a time-delay embedding and its relationship to its derivative is that of sine and cosine.That is, cosine is the derivative of sine, which when shifted by a quarter of its period, perfectly overlaps sine, as exemplified in Fig. 10.The shift that overlaps sine with its derivative is related to its lowestorder period.In a time-delay phase-space, this becomes apparent when sine is plotted against its time-delayed self.As the delay coordinate is increased, the phase-space projection progresses from a diagonal line (for a time-delay of 1 unit) and expands, or unfolds, towards the image of its differential phase-space projection, which would be a perfect circle for sine.Time-delay embeddings are traditionally examined via visual inspection (for other examples see Takens 1981;Casdagli 1992;Sauer 1994).
From Fig. 11, we look for the optimal unfolding with increasing time delay of 4U1705-44 towards its image in the differential phase space projection of Fig. 5, which occurs roughly between 30 and 35 days.As the time delay is a quarter of the signal's period, the lowest-order period from this embedding is between 120 and 140 days.If we were to continue to increase the time-delay coordinate, we would see the phase space projection continue to expand away from the image of the differential phase space embedding and then come back towards it at the next underlying period (e.g towards the next large peak in the power spectrum).This is similar to encountering the successive peaks and zero-crossings of the autocorrelation function of the light curve.For our purposes, we are merely looking for a range of values relating to the lowest-order driving period consistent with the power spectrum and close returns analysis.

EXTRACTING UNSTABLE PERIODIC ORBITS
We have determined a range for the lowest-order period in the 4U1705-44 data.We have also produced numerical solutions of the Duffing oscillator.Its lowest order period is simply the driving period, ω, from the Duffing equation, which we know to be   6 for the reference Duffing solution.Note that the Duffing solution was created to be ten times the length of the 4U1705-44 time series.This close returns map is zoomed in to a shorter portion of the time series in order to highlight the regions of period-1 orbits.The dotted, red line corresponds to p = 140 days, from which we can extract a 1D close returns plot for a single period, as in Fig. 12. 140 days for the solution plotted in Fig. 4. Starting with the numerically generated Duffing time series, we can extract low-order, unstable periodic orbits from the time series by using a 1-D version of the close returns method.For a particular period we can compare each position in the time series to each successive period later, calculate the relative distance and, if this distance falls within a small value, , and evolves within that neighbourhood for at least one pe-riod, then we have located an unstable periodic orbit (Mindlin & Gilmore 1992).This is equivalent to extracting the row located at p = 140 days (highlighted by the red, dashed line) in Fig. 7 for the Duffing time series and plotting its position, i, against .Such a determination is exemplified in Fig. 12 for a period-1 orbit and is distinct from the full close returns mapping as it searches for periodicities of a specific period one at a time.Once we have located the  region in the time series where a periodic orbit exists (e.g. a period-1 orbit is located about 55800 MJD identified by Fig. 12(b)), and do so for as many periods as we can discover in the time series, then we can save each of these periodic orbits in order to compute the relative rotation rates.
Using the method of close returns on the numerically generated Duffing solutions, we were able to extract unstable periodic orbits up to period-6 by varying the parameter n in Fig. 12.That is, we were able to locate and extract 5 period-1, 4 period-2, 3 period-3, 4 period-4, 2 period-5, and 2 period-6 orbits.
Given that we determined the lowest order period in the 4U1705-44 data to be in the rough range of 120 to 170 days, we chose the Duffing period-1 length of 140 days for optimal comparison.Using the same method, we thereby reduced the close returns plot in Fig. 6 to a one-dimensional version as done in the Duffing Figure 10.Plot of sin(x) and cos(x), which are phase-shifted by π/2, a quarter of the lowest-order period of sin(x) of 2π.If we were to apply a time-delay embedding on sin(x), we would recover its derivative, cos(x) using a time-delay of π/2.We seek to find the corresponding time-delay embedded in the 4U1705-44 data that would recover its phase space trajectory.The time-delay is related to the lowest-order period in the 4U1705-44 time series.
case.As a result, we located regions of unstable periodic orbits and successfully extracted three orbits for each period-1, -2 and -3.

RELATIVE ROTATION RATES
A strange attractor can be characterised by the invariants that it possesses.The traditional classification of a strange attractor is by the determination of its dynamical and metric invariants, e.g. the Lyapunov exponents and various dimensions (Correlation or Minkowski).The other type of invariant is topological.Where the first two types are invariant only under coordinate transformations, the topological type is also invariant under control-parameter changes.For experimental conditions in which the control param-Figure 11.Time-delay embedding of 4U1705-44, for time-delays of 10, 30, and 50 days.The point here is to recover an embedding that closely resembles the differential phase space embedding of the time series, Fig. 5, which will be related to an intrinsic lowest order period in the data.The time-delay closest to the derivative phase-space embedding is traditionally determined via visual inspection and, in this case, occurs for time delays of 30-35 days, corresponding to driving periods of 120-140 days.This can be compared to the relationship between sine and cosine, which are derivatives of each other and resemble each other with an offset of some phase factor, or time-delay, related to the driving period.
eters experience perturbations the determination of the topological invariants is appropriate (Solari & Gilmore 1988a) for up to threedimensional attractors.Furthermore, the topological approach can be applied to noisy systems and shorter datasets.
Topological invariants draw dependence on the periodic orbits existing in a strange attractor.The mechanisms that drive the behaviour of a strange attractor uniquely organize all the unstable periodic orbits embedded in that strange attractor.The 'stretching' mechanism is correlated to the sensitivity to initial conditions (nearby points essentially repel each other at an exponential rate), whilst the 'squeezing' mechanism is responsible for keeping the trajectory bounded in phase space.Identifying the organisation of the unstable periodic orbits can be used to identify these underlying mechanisms and thereby provide us with the 'fingerprints' on which the attractor is built (Mindlin et al. 1990).The underlying structure described by these mechanisms is completely responsible for the organisation of all periodic orbits in a flow.This becomes a powerful mathematical tool because if we can make a statement about the equivalence of the underlying mechanisms driving two systems, then we can qualitatively compare these two systems.Most notably, in the case of the Duffing oscillator and 4U1705-44, we can look to the Duffing equation to garner insights into the long-term behaviour of 4U1705-44.
One of the topological invariants introduced by Solari & Gilmore (1988b), intended to describe periodically driven two- x(i + n)] as a fraction of the total amplitude of the time series, plotted as a function of i, for a fixed period n.In this example, we are searching for period-1 orbits, which are of length 140 days; (a) is of a large sample of the time series, (b) is a zoomed in portion highlighted by the red box in (a) in which we find a few candidate period-1 orbits.We are hunting for regions in which the fractional difference between a point and one 140 days later is under 0.1 (or 10 per cent of the maximum flux amplitude) and persisting for at least one period (or 140 days).This shows us the location of potential unstable periodic orbits, identified as a region where the time series is close to repeating itself for at least one period.We can then locate the same regions within the Duffing time series (e.g. about days 66000 MJD), and extract and label these sub-sections of the time series as surrogate period-1 orbits.
dimensional dynamical systems, such as nonlinear oscillators, is called the Relative Rotation Rates (RRR).Generally speaking, RRR describe how one periodic orbit rotates around another; or, more specifically, the average value, per period, of this rotation rate.
The RRR can be computed only after the orbits have been embedded in a three-dimensional space (Gilmore 1998).We chose a three-dimensional differential phase space embedding (Mindlin & Gilmore 1992): Let A be an orbit of period pA which has intersections (a1, a2, . . ., ap A ) with a Poincaré section t = const, and similarly for orbit B. The relative rotation rate Rij(A, B) of A around B for the initial conditions (i, j) is defined as: where ∆r = [xB(t) − xA(t), yB(t) − yA(t)] is the difference vector between points on the two orbits, n is the unit vector orthogonal to the plane spanned by ∆r and d∆r, and the integral extends over pA × pB periods.The initial conditions are the points ai, bj on the Poincaré section and the resulting Rij(A, B) is a single number represented as a fraction for the orbits of A and B. All of the relative rotation rates for a system can be assembled into a table, or 'intertwining matrix' (Solari & Gilmore 1988a).An example of the 3D embedding of an extracted period-4 orbit plotted against another period-4 orbit in the numerical Duffing time series is plotted in Fig. 13.To compute the RRR, we located the intersection of each extracted periodic orbit with the same Poincaré section.Next, we connected each pair of points in the two orbits by a directed line segment.This line segment will evolve in time under the flow and will have undergone an integer number of full rotations (2π radians) in the plane perpendicular to the flow in pA×pB periods.The relative rotation rates can be computed in four equivalent ways as described by Solari & Gilmore (1988a).The first is the mathematical description as defined in Eq. 4 as an integral for continuous datasets.The line integral essentially describes the length of the rotating difference vector per period through phase space.The line integral can be described using graphical means, much like the trajectory of a time series can be plotted and the area under the curve can be computed by eye.One of these other methods, a graphical approach and therefore very straightforward for our means, is as follows: Whenever ∆r 2 = 0 and ∆r 1 0, define σ(t): In summary, we sum the number of times the difference vector crosses the positive x-axis from above or below per total period.Fig. 15 shows each time the difference ∆r crosses the half line ∆r 2 = 0, ∆r 1 > 0, whereby the crossing direction is counted positive (σ = +1) if d∆r 2 /dt > 0 and negative if d∆r 2 /dt < 0. For the period-4 compared against a different period-4, we find the RRR is 3/4 for the Duffing solution.
Following the same procedure with all extracted unstable periodic orbits, we collect the RRR for both 4U1705-44 and the numerically generated Duffing time series into their respective intertwining matrices, where we have given each extracted orbit a different identifying name according to its behaviour in phase space.An example of extracted surrogate, unstable period-2 orbits (akin to Fig. 13) from the 4U1705-44 light curve in a 2D phase-space embedding is plotted in Fig. 14 and the graphical computation of the RRR for the same period-2 orbits is shown in Fig. 16.
In computing the RRR, we noted that the time at which we began collecting data on 4U1705-44 is arbitrary to its mathematical representation.We therefore assumed that the initial point of the first orbit appearing in the time series intersected a Poincaré section and then chose the initial point of the second orbit such that it intersected the same Poincaré section.We then varied the initial point of the first orbit within the range that it was still close to an unstable periodic orbit, matched the second orbit to the new Poincaré section, and computed the RRR again.We found no changes in the computation of the RRR for each pair of orbits when varying the initial point.This was important for self-verification as the numerical Duffing time series is much smoother than the 4U1705-44 time series.Up through period-3 (the longest period extracted from the Table 1.Orbits are labeled as (p.n), i.e. the n th orbit of period p; x and y signify the orbit's presence in each lower-or upper-well, γ signifies symmetric orbits, and α, β asymmetric orbits.Only distinct orbits are included in the matrix.
4U1705-44 data) the matrices are identical.Given that these matrices are suggestively identical, it is therefore possible that the flows are equivalent.This method is essentially a means to falsify a set of differential equations and we have shown that the Duffing oscillator, and its family of differential equations, cannot yet be ruled out as related to the equations of motion generating the 4U1705-44 light curve.As such, if the RRR do continue to agree, then the underlying template of the Duffing oscillator can serve as a geometric model for the dynamics which generate the chaotic time series of the 4U1705-44 system.13 and normalised) of a pair of extracted unstable periodic orbits of period 2 from the 4U1705-44 embedded time series, highlighted within the original time series of 4U1705-44 (Top Panel).These two unstable periodic orbits of period 2 were identified using the Close Returns method, highlighted in Fig. 6.
Table 2. RRR of numerically generated Duffing time series with the same identifying naming convention of orbits as 4U1705-44 RRRs.Orbits of the same type as in 4U1705-44 are included here.

DISCUSSION
In the previous sections we have presented evidence that the time evolution of X-ray binary 4U1705-44 shares key features with the nonlinear Duffing equation in a chaotic regime.The mechanism in 4U1705-44 and other X-ray binaries showing super-orbital variability is not well understood.Our results suggest that a doublewelled, non-linear oscillator is a strong candidate to describe these systems.Specifically, we have found that the low-order driving pe- riod lies near 125 days, as is seen in the power spectra, time-delay embedding and close returns analysis.Similarly, the driving frequency (ω) of the reference Duffing solution has a similar driving period of around 140 days.The relative rotation rates analysis suggests that 4U1705-44 and the Duffing equation share the same underlying template.This means we can look to the Duffing equation to provide insight into the behaviour of the 4U1705-44 system.For example, we can consider the five parameters of the Duffing equation and how each might relate to the dominant physical processes in the neutron star binary system.
The Duffing oscillator can be described by a flexible metal beam attached to a rigid, oscillating frame whereby the beam moves between two magnets (Fig. 17).Returning to the Duffing equation, Eq. 1, we see the flexibility (or stiffness) of the metal beam is described by α, while β describes the non-linearity of the metal beam material, much like a position-varying spring constant.The rigid frame is driven harmonically, described by frequency ω, where the driving has an amplitude of γ.The rigid frame drives both the metal beam and the two magnets on either side, which means that the driving frequency correlates the motion of the metal beam to that of the magnets.Furthermore, strong enough magnets would force the metal beam to oscillate at the driving frequency.
With a picture of a desk-sized oscillator in mind, we can hypothesise the relationship of each component of the toy model to physical mechanisms in the neutron star binary system.Perhaps the most important is the source of the driving frequency.We know that the observed period in the time series of 4U1705-44 is over a hundred days, most likely occurring at approximately 120 days, and is therefore substantially longer than the probable orbital period of 1-10 hours (Homan et al. 2009).We therefore should consider precessional periods that might be present in the system, which would explain longer, super-orbital periods.2013) found compelling evidence that there exists variable neutron star free precession in Hercules X-1 relating to its observed 35-day cycle.Free precession of a body that is not perfectly spherical results in an angle between a point on the body's surface and the total angular momentum vector and thus the location of the magnetic poles.The ellipticity of the neutron star body dictates the length of the precession period (Landau & Lifshitz 2013).For Her X-1, the neutron star non-sphericity should be ≈ 10 −6 resulting in a 35-day free precession period.In order to obtain a period of 120 days, as in the 4U1705-44 system, the nonsphericity would be approximately a third of Her X-1, if we assume the same parameters.That is, the 4U1705-44 neutron star would be more spherical than that of Her X-1.
For Her X-1, the time series is strongly periodic about the 35day period.This suggests that the strength of the driving force is much higher in Her X-1 than in 4U1705-44.That is, if γ determines the strength of the magnets in the Duffing toy model and thus the strength of the driving force, then the strength of the neutron star magnetic field, where the magnetic poles are misaligned, could slave the accretion disc to the precession period in Her X-1.In the case of 4U1705-44, where we observe more chaotic behaviour, the neutron star magnetic field might therefore be weaker, allowing for other parameters in the system to compete with the free precession period.A cyclotron line is present in the spectra of Her X-1, corresponding to a magnetic field of approximately 10 12 G (Shapiro & Teukolsky 1983).In contrast, the non-detection of a cyclotron line by Fiocchi et al. (2007), Piraino et al. (2007), and Lin et al. (2010) and the estimated 10 8 G magnetic field (Lin et al. 2010) of the 4U1705-44 neutron star confirms our theory that 4U1705-44 has a weaker magnetic field than Her X-1.
A weaker magnetic field in 4U1705-44 allows for other physical mechanisms to compete with the free precession period of the neutron star driving the variability.For example, a disc-magnetosphere interaction, also considered by Postnov et al. (2013), occurs when a magnetised neutron star is surrounded by a thin diamagnetic accretion disc.The magnetic field of the neutron star will exert a torque on the accretion disc, driving it out of the equatorial plane (Lai 1999).As a result, the spin axis of a freely precessing neutron star in a binary system can remain tilted with respect to the orbital angular momentum vector thereby maintaining a tilt in the accretion disc out of the equatorial plane which can then precess.
The companion to 4U1705-44 is a low-mass, dwarf star (Homan et al. 2009), a class of stars which are known to exhibit powerful flares on timescales from minutes to hours (Hilton 2011).If the companion star is variable with a high-frequency of flare activity then it could have a significant impact on the observed variability.
We therefore have a cocktail of strong effects potentially producing the observed chaotic behaviour: free precession of the neutron star providing the clock of a 125 day driving period with the neutron star magnetic field modulating its strength; the discmagnetosphere interaction dictating the non-linearity of the accretion disc; the variation in the companion star damping the supply of material in the accretion disc as well as the well-known magnetorotational instability responsible for the viscosity mechanism in the accretion disc dictating the 'stiffness' of the oscillator.
There are other possible mechanisms to explain the long-term variability.For example, Lense-Thirring precession, which is a general relativistic frame-dragging effect due to the spin of a central massive object, has been examined for various parameters of lowmass neutron star binaries by Morsink & Stella (1999) and is found to relate to the observed QPOs in such systems.However, the computed precession frequencies (20-35 Hz) are still many orders of magnitude too short to be related to the driving period of over a hundred days.
Another promising source for the long-term driving period could be due to an irradiation-driven warped accretion disc.According to Pringle (1996), an accretion disc that is initially flat is unstable to radiation-driven warping due to a central illuminating source.The resulting warp can then precess on the order of hun-dreds of days.This effect has been studied using smooth-particlehydrodynamics (SPH) simulations of low-mass X-ray binaries by Ogilvie & Dubus (2001), Foulkes et al. (2006), andFoulkes et al. (2010).The warp that develops due to a central source is predicted to occur in the outer portions of the accretion disc, as opposed to the central flow warp that develops due to the Lense-Thirring effect, and therefore is a more likely source of a longer period driving force.
Any physical model chosen should also be in agreement with the observed spectral state transitions that occur in 4U1705-44 and the corresponding physical explanations.For example, Piraino et al. (2016) used a double Comptonisation model to describe the hard and soft state transitions of 4U1705-44 observed with BeppoSAX.Their model suggests an alternating dominance of soft Comptonisation from a hot plasma surrounding the neutron star in the soft state and a hard Comptonisation arising from the disc in the hard state.
We can verify the template identification of 4U1705-44 by computing the intertwining matrices and linking numbers for other numerically generated Duffing time series as well as from other external computations of the Duffing oscillator (Tufillaro et al. 1990).In a future work, we will extend this analysis to other X-ray binaries to investigate the presence of non-linearity and chaos in their light curves and whether the underlying template of the Duffing equation, or its family of differential equations, continues to be topologically equivalent to those generating the observed variability of other accretion powered systems.Already, the intertwining matrix for another chaotic system, the Belousov-Zhabotinsky chemical reaction (Mindlin & Gilmore 1992), is not identical to 4U1705-44 and differs by a non-integer (thus, we can predict that the flows and the equations that describe them are inequivalent).We also intend, in a future publication, to compare our analysis with the results of other methods aimed at exploring the source of the long periods in X-ray binaries such as simulations of accretion disc systems subjected to a central irradiating source or magnetically driven instabilities.

APPENDIX A: THE METHOD OF SURROGATE DATA
The application of non-linear time series methods should be motivated by evidence of non-linearity in the time series in question.The method of surrogate data testing introduced by Theiler et al. (1992) attempts to find the 'least interesting' explanation that cannot be ruled out based on the data, namely that some type of linear stochastic process is responsible for the observables of interest.In summary, we do this by first formulating a suitable null hypothesis for the observed structures in the time series.We then generate 'surrogate' data that preserves the temporal correlations and probability distribution of the original data using a Fourier-based Monte Carlo resampling technique.We can look for additional structure that is present in the data but not in the surrogates via statistical tests, thereby indicating non-linearity and possibly determinism, depending on the nature of the test statistic used.
The method of surrogate data is different than the statistical approach of choosing model equations fit to the data, producing simulated time series from the chosen model, and then comparing to the original dataset, called the 'typical realisations' approach (Theiler & Prichard 1996).Instead, we follow a 'constrained realisations' approach (Theiler & Prichard 1996) where we produce simulated time series, or surrogates, which are generated by the original dataset itself, thereby imposing the features of interest onto the surrogates whilst eliminating dynamical information, rather than fitting model equations to the data.
We seek to reproduce the bimodal behaviour (the low-and high-flux states) and the strong recurrence phenomena (the nonperiodic recurrence of high and low states) present in the 4U1705-44 light curve and phase-space projection.Our null hypothesis is that a variation of a linear stochastic process produces these observed features (which variation is specified by the chosen test statistic).Our original and unfiltered 4U1705-44 light curve, normalised to unity, is plotted in Fig. A1 where τ is the time delay or lag, N is the total number of data points, y is the average, σ is the standard deviation, and x is the time series (Hegger et al. 1999).That is, we have plotted the correlation of the time series to a delayed copy of itself as a function of the delay.The first peak, which is also the most positive correlation, corresponds to a time delay of 364 days.This timescale is on the order of the period-3 extracted from the 4U1705-44 light curve in the Close Returns analysis in this paper.It is worth noting that strong period-3 occurrences are indicators of chaos (Li & Yorke 1975).We also observe a very slow decay overall in the autocorrelation, with multiple significant correlated and anti-correlated peaks (above or below the dotted red lines of 95 per cent confidence levels in Fig. A1(c)), signifying the presence of both longterm memory and periodic components.

A1 The Surrogate Data
Using the time series analysis code, TISEAN (Hegger et al. 1999, Schreiber & Schmitz 2000), we create a set of surrogate data whereby an iterative algorithm (the routine called surrogates) ensures the surrogates obtain the same autocorrelations as the data and the same probability distribution (Schreiber & Schmitz 1996), the features of interest as in Fig. A1.To produce our confidence level, we use a rank-order test, as introduced by Theiler et al. (1992) and summarised by Schreiber & Schmitz (2000), avoiding the likely misplaced assumption that our signal has a single Gaussian distribution.We can obtain a level of significance of (1 − α) × 100%, where α is the residual probability of a false rejection, by generating M = 2K/α − 1 surrogate datasets where K is an integer resulting in a probability α that the data gives either one of the K smallest or largest values (see Hope (1968) for an extensive proof).
We desire a significance level of at least 95 per cent and choose K = 2 (generally chosen to be a small integer in order to reduce computation time), requiring 79 surrogate time series.Whichever test statistic we use will require the original time series statistic to fall in the top or bottom K = 2 values for a 95 per cent level of significance.An example surrogate time series is plotted in Fig. A2 alongside its distribution and correlogram where we note that both are conserved between the original time series and the generated surrogate.
The test design for generating Fourier based surrogates and other methods and their associated confidence levels are described by Schreiber & Schmitz (2000) and the corresponding C and Fortran based computational methods for surrogate generation and testing are available in the TISEAN package (Hegger et al. 1999).

A2 Test Statistics
There are several test statistics that can detect non-linearity and tracers of chaos to provide motivation for pursuing the non-linear topological approach in this paper.A test often used for detecting non-linearity in a time series is the time reversal asymmetry statistic.That is, a system that is non-linear (deterministic or stochastic) will lead to time irreversibility, or poor forecasting of the timereversed data.Diks et al. (1995) showed that if reversibility can be rejected, all static transformations of linear stochastic processes can be excluded as a model for the original time series.Such a distinction is important because the power spectrum does not contain information about the direction of time and is therefore insufficient for distinguishing linear stochastic and other related processes.The TISEAN package includes the function timerev for the purpose of testing time reversibility, which Schreiber & Schmitz (2000) defines quantitatively as: where τ represents the time delay and N the length of the time series, sN .We compute the time reversibility asymmetry statistic for the data and the 79 surrogates with the results plotted in Fig. A3, in which the time asymmetry of the original data is found to be significantly different from that of the surrogates.That is, the null hypothesis of any static transformation of a linear stochastic process generating the time series of 4U1705-44 can be rejected at a 95 per cent significance level, since the resulting statistic for the data is found to be either the first or the last two values against the surrogates.It should be noted, however, that time reversal asymmetry detects non-linearity, but not the nature of the non-linearity, i.e. deterministic or stochastic.Another method to test the null hypothesis of a stationary, possibly rescaled, linear Gaussian random process generating the observed time series is via a simple non-linear prediction algorithm (Kantz & Schreiber 2004), which was originally proposed as a forecasting method by Lorenz (1963) referred to as the "Lorenz' method of analogues."This method proposes that, if we have observed a system for some time, there will be states in the past which are arbitrarily close to the current state, in which a 'state' is consid- ered to be some segment of the time series, embedded in phase space.We can therefore use a localised segment of the embedded time series to predict, or forecast, future behaviour.
Given a time series y of length N and its corresponding mdimensional embedding in phase space, one searches for all embedding vectors, sn, in the past which are closer than some distance to the current state, sN (i.e. for the last segment in the time series, we search for every other previous segment in the time series which follows a similar time evolution in phase space).The m × region   (Farmer & Sidorowich 1987, Hegger et al. 1999) of the 4U1705-44 data against its 79 surrogates.The original data maintains one of the smallest two errors and is therefore significant at the 95 per cent level to rule out the null hypothesis of a linear stochastic process.centred about sn in phase-space forms a neighbourhood, U (sN ).One then takes the average of all regions, sn, found to be within the neighbourhood U (sN ) and uses it as the prediction for the future evolution of the time series, a distance ∆n away from the end of the time series.In the TISEAN package, we follow the Kantz & Schreiber (2004) definition of the resulting predicted state: The described algorithm above is considered a zeroth-order prediction of the dynamics whereby we are using a locally constant predictor to forecast the future of the time series.Given that this method exploits deterministic structures in the signal, the expectation is that, in applying a prediction algorithm in surrogate testing, the resulting rms error in the prediction algorithm (Farmer & Sidorowich 1987) would be smallest for the original time series than in its surrogates if the null hypothesis is to be rejected (Hegger et al. 1999).For the case of the 4U1705-44 time series and its 79 surrogates, we require the simple non-linear prediction test error to be one of the two smallest errors in order to reject the null hypothesis with a 95 per cent level of significance.Per Fig. ??, this is the case.
The third non-linear statistic that is used for detecting nonlinearity present in the 4U1705-44 time series but which can also indicate determinism is the local versus global linear prediction test as described by Casdagli (1992).By defining a neighbourhood using a local linear ansatz about the data in phase space, one can then compute the average forecast error of a locally linear model as a function of the neighbourhood size.This method is algorithmically similar to the previous simple non-linear prediction test, but with using a local linear first-order model rather than a zeroth order, constant value predictor.In this case, if the agreement with the model is best at large neighbourhood sizes (that is, the error in the forecast is smallest), then the data is best described by a linear stochastic process.If instead the agreement with the model is best at the smallest neighbourhood sizes, as exemplified by the test Given that the time series has been normalised, this means that for a neighbourhood size of 5 × 10 −1 , the distance between predictor and predicted states is more than 70 per cent the maximum distance between any two given points in the time series.in Fig. A5, then a non-linear, possibly deterministic, description is preferred.

A3 Recurrence Analysis
The recurrence analysis tool introduced by Eckmann et al. (1987) and further developed by Marwan et al. (2007) has been used in time series analysis of various experimental data sets for the purpose of distinguishing recurrent phenomena not easily extrapolated by other correlation methods.Recurrence analysis was recently used to characterise deterministic versus stochastic structure in the light curves of black hole X-ray binaries (Grzedzielski et al. 2015, Suková et al. 2016).The production of a recurrence plot (RP) is similar to the Close Returns method used in Sec. 3. Every position in the time series is compared to every successive position in the time series and, if the flux-like value at those two positions are within some distance of each other in phase space, then a black dot is plotted on the RP.
Two of the parameters that must be chosen when constructing an RP are the time delay coordinate and the distance, , used as a threshold for comparing two points in the time series phase space.From the correlogram of 4U1705-44 in Fig. A1, we choose a delay coordinate of 91 (or 364 days).Secondly, we follow the 'rule of thumb' outlined by Marwan et al. (2007) whereby the threshold should not exceed 10 per cent of the mean or the maximum phase space diameter.The resulting RP for the 4U1705-44 light curve is plotted in Fig. A7.
To briefly explore some of the prominent characteristics of RPs, we provide four examples of RPs from time series generated by distinctly different mechanisms.Fig. A6 displays, in order, the RPs of uniform gaussian white noise, a periodic signal constructed from the addition of two sinusoidal waves, the logistic map corrupted by a linearly increasing term (xi+1 = 4xi(1 − xi) + 0.01i), and the chaotic Duffing oscillator (with the same parameters as the simulated time series used throughout this paper).Summarising the extensive analysis from Marwan et al. (2007), Homogeneous RPs, like those from gaussian white noise, are typical from stationary random processes, where no characteristic structure is distinguished perhaps because the timescales of interest are outside the bounds of the RP.Periodic or quasi-periodic systems display shorter diagonal lines parallel to the main diagonal throughout the entire RP.As noted by Marwan et al. (2007), even for systems whose oscillations are not easily recognisable (like in the case for unstable periodic orbits), the presence of short diagonal lines parallel or perpendicular to the main diagonal can be useful.Nonstationary systems give rise to the fading out at the corners of the RPs, as in the case of the modified logistic map (which includes a linearly increasing term).Typical behaviour of laminar states, or intermittency, are the appearance of vertical (horizontal) lines in RPs, reflecting a period in the time series in which the state evolves slowly or does not change.Chaotic systems tend to show a combination of all of these features, except for possibly non-stationarity, as exemplified by the RP of the Duffing oscillator.
The RP of 4U1705-44 shows multiple features of interest, namely, regions of diagonal lines parallel to the main diagonal indicative of fundamental periods, short horizontal and vertical lines indicative of unstable periodic orbits and intermittency, and persistence of intensity throughout the entire RP, reflecting stationarity within the length of the time series.The qualitative similarity to the Duffing oscillator is also remarkable.At a minimum, we can state from visual inspection that there likely exists simultaneous components of noise, periodicity, and intermittent, likely chaotic components in the RP of 4U1705-44.
A more thorough and extensive analysis of the presence of non-linearity and possibly deterministic chaos can be continued on the time series 4U1705-44 using the surrogate data and recurrence analysis methods.For our purposes, we seek only a motivation to pursue non-linear analysis on the time series of 4U1705-44 using a topological comparison to the template of the Duffing oscillator, which can potentially add further justification of the presence of non-linearity and chaos if a positive correlation is found.The surrogate data testing and recurrence plot visual analysis provides ample motivation to pursue a topological analysis of the 4U1705-44 light curve as compared to the Duffing oscillator.

Figure 1 .
Figure 1.Flux of 4U1705-44, 6 January 1996 to 16 May 2014, in MJD.The blue curve represents the data obtained from RXTE ASM, while the orange curve is from MAXI with data taken in the 2-12 keV and 2-20 keV flux bands, respectively.The region in which the blue and orange curves overlap (inset) provides the basis for scaling the MAXI time series to that of RXTE.

Figure 2 .
Figure 2. The mean-subtracted flux of 4U1705-44, 6 January 1996 to 16 May 2014 (in MJD), using the combined RXTE ASM and MAXI data in physical flux units via normalisation to Crab units of each instrument.The grey curve is the raw light curve and associated error.The green curve was produced by applying a low-pass filter.

Figure 3 .
Figure 3. Duffing Equation solution time series segment spanning the same time window as the 4U1705-44 extended light curve.

Figure 4 .
Figure 4. Numerical Flux of Duffing Equation solution against its first derivative, normalised to unity.We note the similar double-welled behaviour in phase space similar to 4U1705-44, and the slight asymmetry in tightness in the two wells.

Figure 5 .
Figure 5.The mean-subtracted flux of 4U1705-44 against its first derivative, over the time period of December 1995 to May 2014 using combined RXTE ASM (1.5-12 keV) and MAXI (2-20 keV) data.The green curve represents a cubic-spline interpolation of the data (blue crosses).

Figure 6 .
Figure6.Close Returns of 4U1705-44: A blue point is plotted at (i, p) when, for that position in time, i, there is a point p days later that is 'close' to the flux value at time i, where closeness is defined as within 10 per cent of the maximum amplitude in flux within the light curve.The intensity of the blue dot correlates to = x(p) − x(i), the difference between the two points in the time series, divided by the maximum amplitude of the time series, i.e. the fractional 'closeness' of the two points in the time series.Regions of the time series that are close to repeating themselves appear as short horizontal lines.One such example is located within the black box aligned with the line in red, which corresponds to the regions plotted in red and blue in the time series of Fig.14.

Figure 7 .
Figure7.The same close returns map as Fig.6for the reference Duffing solution.Note that the Duffing solution was created to be ten times the length of the 4U1705-44 time series.This close returns map is zoomed in to a shorter portion of the time series in order to highlight the regions of period-1 orbits.The dotted, red line corresponds to p = 140 days, from which we can extract a 1D close returns plot for a single period, as in Fig.12.

Figure 8 .
Figure 8. Left panel: The total power spectrum of the 4U1705-44 light curve, normalised according to Shimshoni(1971).Vertical dotted red line corresponds to at least 99 per cent significance in power.Right panel: Dynamic Power Spectrum of 4U1705-44, with a window size of 2048 days, normalised to unity.The edges of the time series are padded with Gaussian noise up to half the window size in order to extract periodicities at the beginning and end of the dataset.We observe the lowest order periodicity occurs between 100 and 150 days, especially early in the time series, with a strong period-3 presence developing later.The longer periods arising near the end of the time series are likely a consequence of the persistent high state in the time series from about MJD 55000 and on.The dashed red lines correspond to periods of 120 and 350 days, for reference.

Figure 9 .
Figure 9. Dynamic Power Spectrum of the reference Duffing solution time series, with a window size of 2048 days, normalised to unity.The edges of the time series are padded with Gaussian noise up to half the window size in order to extract periodicities at the beginning and end of the dataset.The dashed lines are located at periods of 160, 260, and 350 days, for reference.

Figure 12 .
Figure12.Relative distances between the points x(i) and x(i+n), d[x(i) − x(i + n)] as a fraction of the total amplitude of the time series, plotted as a function of i, for a fixed period n.In this example, we are searching for period-1 orbits, which are of length 140 days; (a) is of a large sample of the time series, (b) is a zoomed in portion highlighted by the red box in (a) in which we find a few candidate period-1 orbits.We are hunting for regions in which the fractional difference between a point and one 140 days later is under 0.1 (or 10 per cent of the maximum flux amplitude) and persisting for at least one period (or 140 days).This shows us the location of potential unstable periodic orbits, identified as a region where the time series is close to repeating itself for at least one period.We can then locate the same regions within the Duffing time series (e.g. about days 66000 MJD), and extract and label these sub-sections of the time series as surrogate period-1 orbits.

Figure 13 .
Figure13.A 3D phase space embedding and 2D projection of a pair of extracted unstable periodic orbits of period 4 from the numerical Duffing solutions over a period 4T .

Figure 14 .
Figure14.Bottom Panel: A 2D phase space embedding (akin to the projected panel in Fig.13 and normalised) of a pair of extracted unstable periodic orbits of period 2 from the 4U1705-44 embedded time series, highlighted within the original time series of 4U1705-44 (Top Panel).These two unstable periodic orbits of period 2 were identified using the Close Returns method, highlighted in Fig.6.

Figure 15 .
Figure15.Computation of RRR for two period-4 orbits of Fig.13.Each time the difference ∆r crosses the half line, ∆r 1 , crossing direction is positive (σ = +1); oppositely for negative.This is shown over a period of 4T (rather than p A × p B ).

Figure 16 .
Figure16.Computation of RRR for two period-2 orbits of Fig.14.Each time the difference ∆r crosses the half line, ∆r 1 , the crossing direction is positive (σ = +1); oppositely for negative.We observe similar behaviour to the period-4 Duffing orbits.

Figure 17 .
Figure 17.Toy model of the Duffing oscillator (Kanamaru 2008) whereby the position x(t) of the metal beam oscillates chaotically with time between the two magnets, all attached to a sinusoidally driven rigid frame.
(a).The distribution, or histogram, of the 4U1705-44 time series, Fig.A1(b), reflects the amount of time, fractionally, that the light curve is spent in a specific flux bin and clearly distinguishes the two peaks corresponding to the high-and low-flux states.The correlogram, or autocorrelation function, of 4U1705-44 is computed up to a time delay of half the length of the time series in Fig.A1(c).The autocorrelation function is defined as:

Figure
Figure A1.(a) The original time series of 4U1705-44, unfiltered and normalised to unity.(b) The distribution of 4U1705-44 normalised flux.We have fit a double gaussian model to the histogram in order to accentuate the double-welled features in the time series and for comparison to the surrogates.The fit parameters corresponding to the mean, standard deviation, and amplitude of the first and second Gaussian peaks are (0.095 ± 0.005, 0.073 ± 0.005, 0.061 ± 0.004, 0.59 ± 0.01, 0.16 ± 0.01, 0.035 ± 0.003), where the errors are the standard deviations to the fit parameters.(c) The autocorrelation function of the time series up to a lag half the length of the original light curve.The first (also the largest) peak occurs at a delay of τ = 91, corresponding to 364 days (vertical dashed line).The dotted lines correspond to 95 per cent confidence level of significance.

Figure A3 .
Figure A3.Time reversibility asymmetry statistic as defined by Eq.A2 for the data, marked by the red line and 'x', and the 79 surrogates.As the time reversibility of the 4U1705-44 data is one of the two greatest or least, we can rule out the null hypothesis of any transformation of a linear Gaussian stochastic process with a 95 per cent level of significance.

Figure A4 .
FigureA4.The rms error of the simple non-linear prediction test(Farmer & Sidorowich 1987, Hegger et al. 1999) of the 4U1705-44 data against its 79 surrogates.The original data maintains one of the smallest two errors and is therefore significant at the 95 per cent level to rule out the null hypothesis of a linear stochastic process.

Figure A5 .
Figure A5.Local versus global prediction errors as a function of neighbourhood size, against logarithmic axes, for four different delay embeddings, corresponding to 40, 120, 364, and 460 days, respectively.The neighbourhood size is determined by the square distance about the segment of the time series used for the prediction algorithm in the embedded phase space.Given that the time series has been normalised, this means that for a neighbourhood size of 5 × 10 −1 , the distance between predictor and predicted states is more than 70 per cent the maximum distance between any two given points in the time series.

Figure A7 .
Figure A7.The recurrence plots (RP) of 4U1705-44 for a time-delay embedding of phase space using a time delay of 91 (or 364 days).A black dot is plotted when the position at t is within = 10% × D of t + ∆t, when embedded in phase space, where D is the maximum phase space diameter.