A systematic study of the phase difference between QPO harmonics in black hole X-ray binaries

We perform a systematic study of the evolution of the waveform of black hole X-ray binary low-frequency QPOs, by measuring the phase difference between their fundamental and harmonic features. This phase difference has been studied previously for small number of QPO frequencies in individual sources. Here, we present a sample study spanning fourteen sources and a wide range of QPO frequencies. With an automated pipeline, we systematically fit power spectra and calculate phase differences from archival Rossi X-ray Timing Explorer (RXTE) observations. We measure well-defined phase differences over a large range of QPO frequencies for most sources, demonstrating that a QPO for a given source and frequency has a persistent underlying waveform. This confirms the validity of recently developed spectral-timing methods performing phase resolved spectroscopy of the QPO. Furthermore, we evaluate the phase difference as a function of QPO frequency. For Type-B QPOs, we find that the phase difference stays constant with frequency for most sources. We propose a simple jet precession model to explain these constant Type-B QPO phase differences. The phase difference of the Type-C QPO is not constant but systematically evolves with QPO frequency, with the resulting relation being similar for a number of high inclination sources, but more variable for low-inclination sources. We discuss how the evolving phase difference can naturally arise in the framework of precession models for the Type-C QPO, by considering the contributions of a direct and reflected component to the QPO waveform.

. QPOs have subsequently been observed in many BHXRBs, especially after the launch of the X-ray timing observatory the Rossi X-ray Timing Explorer (RXTE). In BHXRBs, QPOs are observed at both high (50 to hundreds of Hz) and low (up to 10-20 Hz) frequencies. Here, we focus on low-frequency QPOs, where three different classes of low-frequency QPOs have been identified, based on their frequency, Q-value, amplitude, noise, and phase lag properties (Wijnands, Homan & van der Klis 1999;Remillard et al. 2002;Casella, Belloni & Stella 2005). The Qvalue is defined as the frequency of the QPO divided by its fullwidth half-maximum (FWHM), which is a measure for the width of the QPO peak in the power spectrum. In this paper, we will focus on the analysis of Type-B and Type-C QPOs. The main difference between these two types lies in the shape and amplitude of the broadband noise component in the power spectrum: this noise component has significantly lower amplitude (∼1 per cent rms versus > 10 per cent rms) when associated with Type-B QPOs than with Type-C QPOs. Most QPOs analysed in this work are Type-C QPOs, since these are the most commonly observed, with a fundamental frequency ranging between 0.1 and 12 Hz, while the Type-B QPOs tend to cluster more between 4 and 6 Hz. We will ignore QPOs in the mHz and kHz range, which have also been observed (Wijnands et al. 1997;Linares et al. 2010;Belloni, Sanna & Méndez 2012;Huppenkothen et al. 2017;Rapisarda, Ingram & van der Klis 2017). QPO rms spectra indicate that the emitting region of Type-C QPOs appears to be the corona: an optically thin (optical depth τ ∼ 1) cloud near the black hole where cool disc photons are Compton up-scattered by energetic electrons (Sobolewska & Zycki 2006). Different mechanisms have been suggested to explain the existence of QPOs. The fundamental distinction between these mechanisms lies in the origin of the observed variability, which can be divided into two categories: geometric or intrinsic. A geometric model assumes changes in (apparent) accretion geometry as the QPO origin, while the emission is intrinsically constant. An intrinsic model bases the QPO origin on modulations of intrinsic emission, e.g. driven by changes in accretion rate. As we will discuss below, recent evidence strongly favours geometric models, in particular, the relativistic precession model (Ingram, Done & Fragile 2009;Veledina, Poutanen & Ingram 2013). The key aspect of this model is the precession of a radially extended region of the hot inner flow of the accretion disc. If this inner flow -which is hot and geometrically thick inside the inner radius of a truncated accretion disc -has a spin that is misaligned with the black hole, it will start to precess due to (relativistic) Lense-Thirring torques. Changes in the observed solid angle of the flow and changes in Doppler boosting can then explain the quasi-periodic variations in flux.
Recently, evidence for a general geometric QPO origin has been found in the inclination dependence of QPO amplitudes (Motta et al. 2015;Heil, Uttley & Klein-Wolt 2015) and the sign of lags between soft and hard photons (van den Eijnden et al. 2017). Recent research also shows that the reflection of photons from the inner flow of the disc depends on the phase of the QPO cycle (Ingram & van der Klis 2015;Ingram et al. 2016Ingram et al. , 2017, which is a direct prediction of the relativistic precession model. While the relativistic precession model works well for Type-C QPOs, there are hints that Type-B QPOs also have a geometric origin (Motta et al. 2015;Stevens & Uttley 2016). Finally, for the BHXRB GRS 1915+105, the energydependent behaviour of its low-frequency QPO can be explained by an extension of the solid-body precession model, where the inner flow precesses differentially (van den Eijnden, . For a better understanding of the origin of QPOs, spectral-timing methods provide a powerful approach. Such methods combine both the observed variability and spectral information to reach beyond the information available in each individually. For QPO studies, such methods naturally focus on studying spectral changes throughout the QPO oscillation cycle (Ingram & van der Klis 2015;Stevens & Uttley 2016). However, this approach fundamentally assumes that one can define a phase within the QPO cycle -in other words, such methods assume the presence of a well-defined and stationary waveform throughout an observation. Interestingly, if such a waveform exists, it not only validates phase-resolved spectroscopy of the QPO but might also directly encode information about the origin of the QPO. Therefore, in this work, we set out to characterize QPO waveforms across a sample of BHXRBs over a wide range of their QPO frequencies. We use a statistical approach to determine whether a consistent phase relation between the QPO fundamental and harmonic exists and thereby (i) confirm whether a meaningful underlying waveform is present and (ii) search for signatures of Figure 1. PSD for observation 20402-01-16-00 from GRS1915+105 in the 2-13 keV energy range. The Poisson noise, incorporated in all models, is subtracted from the power which is then multiplied by frequency by convention. The reduced χ 2 /ν-values for models one to three, which allow for increasingly more harmonic components, respectively, are 5.929, 2.079, and 1.569. the origin of the QPO harmonic. Fourteen different sources are analysed, so the effects of inclination, QPO type, and frequency on the waveform can be determined. To do so, the process of extracting the QPO properties from the original light curve is automated and the phase difference is calculated as a function of QPO frequency for a large sample of RXTE observations.

Rationale: a quasi-periodic waveform?
A periodic signal x(t) consisting of two harmonics can, in the most general mathematical case, be written as: where subscripts f and h correspond to fundamental and harmonic, respectively, A is the variability amplitude, ν the oscillation frequency (ν h = 2ν f ), and φ is the phase offset. The absolute phase of the fundamental can be defined arbitrarily; in this work, we set φ f = 0. If all amplitudes, frequencies, and phase-offsets are constant, this waveform will remain the same over time.
The waveform of a QPO can be thought of in the same way as the deterministic waveform introduced above, consisting of multiple harmonic waves combining to a single waveform. The harmonic feature with the largest amplitude determines the main QPO frequency and is referred to as the fundamental in this paper. Following common practice, the second harmonic, at twice the fundamental frequency, will simply be referred to as the harmonic from here on. This harmonic can be much smaller in amplitude. An example of a QPO power spectrum in Fourier space is given in Fig. 1. From left to right, the spectrum shows three bumps: the subharmonic, fundamental, and harmonic, respectively. By expressing the QPO waveform in terms of equation 1, we will effectively ignore any sub-harmonic and higher order harmonic terms, approximating the QPO signal as just the sum of the fundamental and harmonic. While visible in Fig. 1, such features are rarely present in BHXRB observations.
Considering a QPO following equation (1), the quasi-periodic nature arises from the fact that the frequency/phase-offset and/or amplitude of each harmonic does not remain constant with time. However, for the waveform to be well-defined and stationary, each of these parameters should have a preferred value during the observation. The QPO frequency and amplitude have such a value, measurable from the power spectrum (see, e.g. Fig. 1). The phase difference, , cannot however be determined from the power spectrum. This phase difference is defined as: This phase difference is the parameter of interest in this study, as we aim to test whether it has a preferred value in each observation; if so, together with the fundamental and harmonic frequency and amplitude, a well-defined waveform exists. If not, the waveform would change continuously and erratically during observations. In this section, we will introduce each step of the method to calculate the phase difference in depth. Afterwards, Section 3 will describe the data we used to obtain our results and their automatic extraction from the RXTE data archive. A schematic summary of the pipeline, applied to each observation, can be found in Fig. 2, while each of the following subsections corresponds to one step in this figure. Schematically, the methodology is set-up in such a way that first, starting from individual light curves, the power spectrum is created. This power spectrum is then fitted with different models to determine the QPO and harmonic parameters. Using these parameters, a Fourier analysis of individual small segments of the light curve is performed. This new power spectrum is finally used to find the phase difference, which is defined as in equation (2). We stress that it is not possible to calculate the phase difference directly from the cross-spectrum; we refer the reader to the start of Section 3 in Ingram et al. (2017) for a detailed explanation.
The analysis method introduced here is based on the one developed by Ingram & van der Klis (2015), but generalized and automated in such a way that a single pipeline can be used for a large set of light curves. As stated, the next sections will provide a more detailed description of the steps of the analysis method.

Creating the power spectrum
We calculate the power spectral density (PSD) in fractional rms normalization (Belloni & Hasinger 1990) for all individual light curves using a fast Fourier transform (FFT) algorithm. We average over 8 s segments with a time step of dt = 1/128 s, resulting in ensembleaveraged power spectra in the 0.125 to 64 Hz frequency range.

Fitting the power spectrum
We look for QPO fundamental and harmonic features in the PSDs by fitting with a multi-Lorentzian model (following e.g. Belloni & Hasinger 1990). We use two zero-centred Lorentzians for the broadband noise (BBN) (following van den Eijnden et al. 2016) and a narrow Lorentzian for each QPO harmonic. We consider a maximum of three harmonics (fundamental, first overtone, and either second overtone or sub-harmonic) and use F-tests to determine how many harmonics to include in our best-fitting model. We use a threshold of p = 0.01 to determine whether the addition of a harmonic Lorentzian component is required by the data. Only observations with at least one significant overtone are useful for our analysis. An example of these three models fitted to a PSD is given in Fig. 1 for RXTE ObsID 20402-01-16-00 (source GRS1915+105), which clearly shows the improvement of the fit as the model complexity increases.
To be able to fit the power spectrum to the different models, an initial estimate for the QPO fundamental frequency is required. We obtain such an initial estimate by linearly rebinning the power spectrum by a factor three, averaging out statistical outliers. The frequency with the maximum jump in power (after subtracting the Poisson noise) in this rebinned power spectrum is our initial estimate for the QPO frequency. We use P noise = 2/ x where x is the mean count rate and the background is ignored; e.g. Van der Klis (1989); Uttley et al. (2014)). At this frequency, the power change is the most extreme and therefore represents an estimate of the QPO fundamental frequency at the resolution of the PSD. This method only yields the frequency of the sharpest peak and therefore isn't affected by other bumps in the PSD. Averaging over three data points ensures that there is a real peak building instead of one data point that is higher at random. A complete list of initial guesses and boundaries on fit parameters can be found in Appendix B.

Phase difference
In order to measure the phase of the QPO fundamental and harmonic components, we first take the Fourier transform 1 of many short segments of the light curve. Following Ingram & van der Klis (2015), we set the segment length for each observation to ensure To show a clear peak, we repeat the data from ψ/π = 0 to 1 to show two cycles. The probability that this distribution is randomly drawn from a uniform distribution is 1×10 −24 .
that the QPO fundamental lies in a single frequency bin. The harmonic feature is therefore spread across two frequency bins. The phase of the fundamental for a single segment is calculated as the argument of the Fourier transform for that segment at the frequency bin containing the QPO fundamental. For the harmonic, we use the frequency bin at exactly twice the frequency of the bin used for the fundamental. We then calculate the phase difference using equation (2). Note that throughout the remainder of this work, we will consider the phase difference as a fraction of π radians (which we denote ψ/π ), to express it more clearly as a fraction of a cycle of possible values.
It is essential to calculate the phase difference for each segment individually, since the phase of the QPO fundamental and harmonic is random, and only their phase difference is constant. Averaging over the phases of the QPO fundamental and harmonic first and calculating the difference afterwards would result in a phase difference of zero (averaging over the difference of two random quantities). Since we calculate this phase difference for each segment, we can visualize the results in a histogram. An example of such a phase difference histogram is given in Fig. 3 for observation 20402-01-16-00 from GRS 1915+105 in the 2-13 keV energy range. We have repeated the ψ/π = 0 to 1 interval twice in the figure, to show the peak clearly. This histogram shows a clear peak around /π = 0.1, which is repeated at /π = 1.1. The histogram is normalized, such that the relative frequency of measurements at the largest peak is set to 1. To check whether the histogram shows a significant peak, we use a Kuipers test. 2 This Kuipers test is similar to a KS-test, but adapted for cyclic quantities, and returns the probability that the phase difference histogram is randomly drawn from a uniform distribution. Histograms are labeled significantly peaked using a significance level of 3σ .
For phase difference histograms with a significant peak, the peak mean is calculated by summing all distances for all m as follows: and minimizing the result with respect to mean . mean has to be calculated in this way to deal with the cyclic nature of the phase difference. This can most easily be visualized by imagining a circle with circumference 2π that describes all possible phases. Given two points on the circle, one for the fundamental and the other for the harmonic, we want to know the shortest distance between the two. The above calculation takes into account that this distance can be either clockwise or counterclockwise. This visualization of phases on a circle also explains why the phase differences range between 0 and π rad instead of 2π rad.

Confidence interval on the phase difference
The final step in our analysis is to determine an error margin on the mean phase difference mean . Via bootstrapping, a statistical test that relies on repeated iterations of random sampling with replacement, we can assign an uncertainty to the phase difference. This is done by randomly drawing with replacement from the phase difference histogram. Using the mean of this new sequence and repeating this process multiple times, we build a new histogram that shows the bootstrapped mean phase differences. An example (using the same observation as shown in Figs. 1 and 3) is shown in Fig. 4, which shows the distribution of inferred phase difference obtained from bootstrapping. The 1σ uncertainty of the phase difference is determined as the standard deviation of this bootstrapped distribution. For another example of how bootstrapping can be used to estimate parameter accuracy in spectral-timing studies, we refer the reader to Stevens & Uttley (2016). The number of bootstrap realizations of the data is determined using: N boot = 50 − int(number of phase differences/1000). This implies that the maximum number of bootstrapped mean phase differences used to determine the 1σ uncertainty is 50 and will become smaller for long observations. The minimum number of bootstrap realizations is set to 15. This simple formula for the number of bootstrap realizations also optimizes the tradeoff between a high number of bootstrap realizations and short computing time. To account for the cyclic property of the phase difference, the bootstrapped measurements mean,bootstrapped /π ∼ 0 are shifted up by 1 if the actual value of mean /π ∼ 1. This way mean,bootstrapped /π = 0.1 is, for example, shifted to 1.1 such that we obtain a distribution around mean /π = 1, without enormous gaps. The same principle applies to values of mean,bootstrapped /π ∼ 1, which are shifted down by 1 if the phase difference measured from the data itself ∼0.

RXTE O B S E RVAT I O N S A M P L E
For our analysis, we use the sample of Rossi X-ray Timing Explorer (RXTE) observations compiled by Motta et al. (2015) and further  Motta et al. (2015) in their identification of QPO type, which is based on Casella et al. (2005), and considers: the total rms variability in the power spectrum ( 10 per cent for Type-C QPOs, ∼5 − 10 per cent for Type-B QPOs), the QPO frequency and its evolution during the outburst, the Q factor, and the shape of the noise component (flattop-dominated for Type-C QPOs, red-noise for Type-B QPOs).
The previous analyses of this data sample have shown that it is representative of QPO behaviour in BHXRBs and covers a range of source inclinations and possible QPO characteristics. A complete overview of the basic QPO and BBN characteristics of these observations, such as QPO frequencies, amplitudes, and phase lags, can be found in the online materials of Motta et al. (2015) and van den Eijnden et al. (2017). In Table 1, we list the number of Type-C and Type-B observations in the initial sample per source. We reiterate that, since QPOs are not present in the soft state or in dim hard states, only a small fraction of RXTE observations of compact objects shows statistically significant QPOs, depending on the observing strategy as well as the relatively rapid nature of the state transitions which show easily detectable QPOs. Therefore, the percentage of observations that contains QPOs differs per source and can range from less than one to ∼20 per cent. Table 1 also shows the number of observations where we were able to successfully define a phase difference. The rest of the observations either lack a significant harmonic component or do not show a welldefined phase difference . The lack of significant harmonic could follow from short exposures or low fluxes, causing any harmonic features to become undetectable in the noise -especially QPOs with fundamental frequency 5 Hz are prone to such effects, as their harmonics lie in lower signal-to-noise parts of the power spectrum. The lack of a measurable phase difference could also arise from short observations, resulting in few segments to calculate individual phase differences and therefore no significant peak in their distribution. Alternatively, a harmonic could simply be absent -any model for its origin and properties should therefore be able to account for a weak power or even its absence. To extract light curves for every observation, we use the automated CHROMOS pipeline 3 (Gardenier & Uttley 2018). This pipeline automatically extracts light curves from raw RXTE data in a specified energy range and with a pre-defined time resolution, independent of observation mode. We extract all light curves with a 1/128 s time resolution and a 2-13 keV energy range. As the gain of RXTE's Proportional Counter Array (PCA) changes with time and the energy resolution depends on observing mode, CHROMOS selects the absolute PCA channels most closely matching the energy range above for each individual observation. We also extract soft and hard light curves in the 2-7 and 7-13 keV ranges, matching those used in the QPO phase lag analysis of this sample by van den Eijnden et al. (2017), to search for any energy dependence in the measured phase difference.

R E S U LT S
Figs. 5 and 6 show the phase difference between the fundamental and harmonic as a function of the fundamental QPO frequency for Type-C and Type-B QPOs, respectively. To account for the cyclic property of the phase difference, we subtract 1.0 from data points with /π > 0.8 such that these points are mapped to the [ − 0.2, 0] range. For clarity in Fig. 5, only half of the sources are shown in colour per panel. The other half of the sources are plotted in grey in the background, and shown in colour in the other panel. From Fig. 5, it is clear that all sources have a well-defined phase difference i.e., a well-defined QPO waveform, that evolves with QPO frequency. The phase difference decreases with frequency for all sources except for Swift J1753.5-0127 (magenta open squares in the left-hand panel) and XTE J1650-500 (yellow crosses in the left-hand panel). Furthermore, most sources have a phase difference of /π ∼ 0.25-0.4 at 0 Hz, and this phase difference decreases to around 0 at higher frequencies. Exceptions to this behaviour are found in GX 339-4 (light blue leftwards pointing triangles in the left-hand panel) and GRS 1915+105 (yellow dots in the right-hand panel), which start at /π ∼ 0.7-0.8 at 0 Hz. Finally, there are two data points that seem to be outliers at /π ∼ 0.55 and ∼5.5 Hz in the right-hand panel (XTE J1550-564 and GRS 1915+105).
Type-B QPOs show a completely different behaviour, as shown in Fig. 6. In general, the phase difference is constant with frequency at /π ∼ 0.55-0.6. Obvious exceptions are found in H1743-322 (green upwards pointing triangles) and XTE J1817-330 (purple dots) that show a wide variety of phase differences.
The fact that sources show a systematic and well-defined phase difference for a given QPO frequency and QPO type, implies that for a given source's QPO type and frequency, there exists a simple underlying waveform of the QPO, which remains well-defined over each full observation. The presence of such an underlying waveform is a fundamental assumption of many recently developed QPO phase-resolving methods (Stevens & Uttley 2016;Ingram & van der Klis 2015): such methods aim to measure changes in the energy spectra as a function of phase within the QPO cycle, which evidently requires such a phase to exist throughout the analysed observation.
For both the Type-B and Type-C QPOs, we made a preliminary analysis of the energy dependence of the phase difference by analysing light curves in the 2-7 and 7-13 keV bands. This analysis does not reveal any obvious dependence on energy. Therefore, we do not plot these extra energy bands in Figs. 5 and 6. The lack of energy dependence might be explained by the broad nature of the two bands,  averaging out any more subtle changes with energy. Alternatively, as discussed in Section 5.4, it might be more fundamentally related to the origin of the QPO harmonic in relation to the energy range probed by RXTE PCA.

D I S C U S S I O N
In this section, we will first discuss the search for any inclination dependence of our results. Afterwards, we introduce a new approach to identify QPO types based on the phase difference between the fundamental and harmonic. We then continue by discussing geometric models for the Type-B and Type-C QPO and whether these can explain the observed evolution of the QPO phase difference, and finish by recommending future steps for the observational and modelling side.

Waveform inclination dependence
Different QPO properties -such as the amplitude of Type-B and C QPOs (Heil et al. 2015;Motta et al. 2015) and the Type-C fundamental energy-dependent phase lag (van den Eijnden et al. 2017)have been shown to depend on the system binary orbit inclination, favoring a geometric QPO origin. The question therefore arises as to whether an inclination dependence can be found in the evolution of the QPO waveform. Our sources show a variety of relations between the QPO frequency and the phase difference between harmonic and fundamental. To show these differences more quantitatively, we fit a straight line to the data in Figs. 5 and 6 for each source individually. The fitted slopes and inclinations for each source are summarized in Tables C1 and C2 for Type-C and Type-B QPOs, respectively, in Appendix C.
Although there is no single mapping of phase difference evolution to system inclination, it is worth noting two apparent patterns for the Type-C QPOs. First, there is a wide scatter in relations for low-inclination BHXRBs, with phase-difference either increasing or decreasing with frequency, or showing no obvious frequencydependence. Meanwhile, the phase difference decreases with frequency for all high-inclination BHXRBs, with data from four sources (i.e. except GRS 1915+105) clustered around a similar relation. The inclination of XTE J1859+226 is ill-defined (see Appendix C) but its energy-dependent phase-lag evolution is consistent with a high-inclination source (van den Eijnden et al. 2017) and its phase difference evolution is also consistent with that of the majority of high-inclination sources. Since the overall sample size is small and there is a wide range of patterns for the lowinclination sources, we cannot state that the differences are formally significant. However, they are certainly suggestive that there may be an inclination effect on the consistency of the phase difference evolution with frequency, with high-inclination sources showing a more consistent pattern.
Although detailed modelling is beyond the scope of this paper, we can speculate how this behaviour may be interpreted in the context of the precessing inner flow model (Ingram et al. 2009). In this model, the precession period of the inner flow modulates the X-ray light curve via a number of mechanisms: variations during each precession cycle in the solid angle subtended by the inner flow, in Doppler boosting, and in the luminosity of seed photons from the disc intercepting the inner flow. The solid angle peaks at the point in the precession cycle when the flow is seen maximally face-on and Doppler boosting peaks when it is seen maximally edge-on (Veledina et al. 2013;. The observed flux as a function of the angle between the flow rotation axis and the observer's line of sight is then a trade-off between these two effects. Crucially, if there are no more modulation effects, the observed flux is a monotonic function of this angle, such that the waveform we see from a given precessing inner flow system depends only on the angle between the black hole spin axis and the observer's line of sight.
This azimuthal symmetry is broken if the outer disc is misaligned with the black hole spin axis by some angle β and the inner flow precesses around the spin axis maintaining a constant misalignment of β, as suggested in Ingram et al. (2009). The angle between the disc and flow rotation axes therefore varies over the precession cycle from a minimum of 0 to a maximum of 2β (e.g. Veledina et al. 2013;. As the angle between disc and flow changes, the number of seed photons from the disc intercepted by the flow changes, therefore varying the overall luminosity budget of the inner flow (Zycki, Done & Ingram 2016;You, Bursa & errorZdotycki 2018). The QPO waveform now depends on the (polar) angle between the spin axis and the observer's line of sight and on the azimuthal angle, defined, for instance between the projections on the black hole equatorial plane of the observer's line of sight and the line of nodes between the black hole and disc rotation axes (e.g. see the schematic in . Perhaps for high inclinations, the solid angle and Doppler effects dominate over the variable seed photon effect. In this case, we would expect a roughly similar QPO waveform for a given geometry of the inner flow, with the waveform changing as the flow outer radius moves inwards. If the variable seed photon effect becomes more important for low-inclination objects, then different low inclination sources would be expected to display different QPO waveforms, since we are presumably viewing them from a range of azimuthal angles.
For Type-B QPOs there is also no clear relation between inclination and phase difference. From Fig. 6 it seems that the phase difference is constant with frequency at around /π ∼ 0.6 for 6 sources with Type-B QPOs, with the exceptions being H1743-322 and XTE J1817-330.

A quantitative method to classify QPOs
Examining the results of our analysis, shown in Figs 5 and 6, two distinct types of QPO behaviour can be distinguished: on one hand, a large subset shows a decreasing phase difference with QPO frequency, as seen in Fig. 5. A smaller set of observations show QPOs in a small region of the QPO frequency-phase difference parameter space, where /π > 0.5 and ν > 4 Hz. This result represents a clear, qualitative method to classify QPOs into different types, based on only the QPO's properties and directly connected to its physical origin (which sets both the QPO frequency and phase difference).
An obvious question is then how this classification relates to the existing Type-B and Type-C classes of QPOs. The observations  shown in Figs. 5 and 6 are divided based on their prior separation into Type-C and Type-B QPOs, respectively. Therefore, our quantitative classification based on QPO properties alone almost perfectly overlaps with the existing classification, which relies on a more qualitative evaluation of a combination of QPO and BBN properties (see e.g. Casella et al. 2005;Motta et al. 2015).
Considering Fig. 5 in detail however, there are two observations identified as a Type-C QPO, which fit into our second category based on their phase difference (i.e. with /π > 0.5 and ν > 4 Hz). In Fig. 7 and in Fig. 8, the power spectra of these two outliers are shown, together with power spectra of observations of the same source at similar QPO frequency. As shown in Fig. 7, the power spectrum of the outlier XTE J1550-564 observation is very similar to another Type-B QPO power spectrum (shown in red). It does not, however, look like the Type-C QPO power spectra shown as well in grey. Therefore, we conclude that this observation was most likely mis-classified as a Type-C QPO in Motta et al. (2015) and van den Eijnden et al. (2017), showing the power of our more quantitative approach. Similarly, Fig. 8 also appears different from the Type-C QPO power spectra shown: while the QPO itself has similar width and frequency, the noise level at lower and higher frequencies is significantly lower. While an unambiguous re-classification is difficult, it appears more similar to the Type-A QPO power spectra seen in for instance XTE J1550-564 (Wijnands et al. 1999) or MAXI J1535-571 (Stevens et al. 2018).

Interpretation of the Type-B phase difference
In this section we describe a simple geometric model that models the phase difference for Type-B QPOs. In our work we have not considered whether non-detections of a phase difference imply that in some cases there is no meaningful waveform (implicitly assuming that non-detections are due to signal to noise limitations alone), or the implications of the non-detections of harmonic components (which would also require consideration of the upper limits on harmonic amplitudes). However if, in future studies, cases arise where there is clearly no meaningful waveform, or very weak harmonics, these would also need to be accounted for in any physical framework. The phase difference for the Type-B QPOs seems to be approximately constant at around ψ/π 0.5-0.6 for most sources. Recent studies suggest that Type-B QPOs could be related to a black hole with precessing jets (Stevens & Uttley 2016;Sriram, Rao & Choi 2016). Therefore, we present a phenomenological model of a simple jet geometry in this section and analyse the resulting phase difference.
In this model, we set up the jet simply as a narrow rotating cylinder (a 'stick'), where the relativistic boosting and solid angle effects of an optically thick or thin emission region are taken into account (note that optically thin and thick refer to the X-ray properties and not the radio emission, and are therefore not related to the presence of either a steady compact jet or transient ejecta). This is a very simple representation of reality since the intrinsic angular dependence of emission from the jet may also play an important role. The observed luminosity for an optically thin jet can be approximated as (e.g. Ellis 1971): where β is the jet speed in units of c and γ is the jet Lorentz factor. For an optically thick jet, the whole jet is not observed all the time, so the solid angle has to be taken into account: Here, θ is the angle between the line of sight of the observer and the jet. θ can simply be calculated by taking the dot product of the jet and observer vector cos(θ ) = sin(θ incl ) · sin(θ jet ) · cos(φ) + cos(θ incl ) · cos(θ jet ).
Here, θ jet is the angle between the jet axis and the precession axis (i.e. the half opening angle of the precession cone), θ incl is the angle between the observer's line of sight and the precession axis and φ is the precession phase. In the following simulations, we set θ jet = 10 • , β = 0.15 (Miller-Jones, Fender & Nakar 2006), either θ incl = 40 • or θ incl = 70 • , and φ varies from 0 to 2π during each precession cycle. We assume the intrinsically emitted luminosity is constant. Fig. 9 shows the result of this model: the left panels show the observed waveform for an optically thick (red line) and an optically thin (black dashed line) jet. The right panels shows the power spectrum of the waveform. For both inclinations the harmonic arises due to the harmonic content required to describe the non-sinusoidal, symmetric waveform, instead of a physical process occurring at twice the precession frequency. Furthermore, from this simulation it follows that independent of inclination the phase difference for an optically thick jet is /π = 0.5, close to the frequently-observed values. However, for the optically thin (i = 40 • and i = 70 • ) scenarios in Fig. 9 one might not find a significant harmonic component in an observation due to the low harmonic contribution.
In this model, only the simplest geometric effects of a 'stick-like' optically-thick emitting jet are taken into account, so that a more realistic geometry (perhaps a composite of jet-like and disc-like), angle-dependent emissivity and general relativistic effects (likely significant in the part of the corona closer to the disc-plane) are not considered. Combining these effects with the simpler geometric effects described here could produce deviations from /π = 0.5 and perhaps explain the observed Type-B phase differences of /π = 0.5-0.6.

Modeling
The phase difference for the Type-C QPO systematically decreases with QPO frequency, for all except two sources. Such an evolving phase difference can be explained if the observed QPO consists of different contributions. For instance, in geometric precession models, the observed QPO waveform can be a combination of the direct emission from the precessing corona and reflected emission off the disc. The harmonic feature in the direct emission could be caused by general-relativistic effects distorting the waveform into a non-sinusoidal shape. The reflected component, on the other hand, could contain a feature at twice the QPO frequency because the precessing flow illuminates each disc azimuth twice per QPO cycle, with its back and front side, as implied by phaseresolved spectroscopy of H1732-322 . The two contributions to the net QPO waveform would likely not have the same phase difference.
Mathematically, this idea corresponds to introducing two QPO waveforms, both with a different value of that is constant as a function of QPO frequency. If the relative strength of these two waveforms varies with QPO frequency, the net phase difference of the observed waveform will evolve with QPO frequency as well. We can define a waveform x c (t) and x r (t) for the direct coronal and reflection QPO contributions, respectively, as: and where the subscripts f and h refer to fundamental and harmonic, while the c and r correspond to coronal and reflected, respectively. Here, φ r is the phase difference between the fundamental components of the direct and reflected emission. The waveforms in equations (7) and (8) are equivalent to equation (1), using the definition of the phase difference in equation (2). The superposition of the above two waveforms yields a full waveform with the well-defined phase difference , that depends on seven parameters of the original waveforms (namely all four amplitudes, the two phase differences, and the phase φ r in equations (7) and (8)) in a non-trivial fashion. By fine-tuning any of these parameters (or more specifically, their dependence on QPO frequency), one can in principle reproduce any relation between and QPO frequency.
However, it is difficult to match physical scenarios for the QPO origin on to a single phase difference. For instance, both the phase difference between the fundamental and harmonic variation in the reflected flux, and the exact non-sinusoidal shape of the waveform of a precessing inner flow, are unknown and non-trivial to determine: they depend on for instance on the balance between modulation mechanisms discussed in Section 5.1. Detailed examination of the different scenarios will require simulation of the relevant relativistic effects and radiative transfer, which will be an extensive effort and well beyond the scope of this work. In our work we have not considered whether non-detections of a phase difference imply that in some cases there is no meaningful waveform (implicitly assuming that non-detections are due to signal to noise limitations alone), or the implications of the non-detections of harmonic components (which would also require consideration of the upper limits on harmonic amplitudes). However if, in future, cases arise where there is clearly no meaningful waveform, or very weak harmonics, these would also need to be accounted for in any physical framework.

Observations
We speculate that multiple mechanisms underlie the harmonic feature in the QPO, changing in relative strength as the QPO frequency and the outburst evolve. Both a reflection and a GR-based origin of the harmonic signal, discussed in the previous sections, could show a strong dependence on energy: reflection spectra are strongest below 2 keV, at the Fe Kα line, and above 15 keV (the Compton hump). Alternatively, the GR effects might dominate in the power law emission that originates from the innermost regions of the accretion flow. Therefore, if the harmonic contributions of one or both of these two possible mechanisms, this could be visible by considering the phase difference (evolution) in different energy ranges.
In this work, we have considered the QPO waveform in a broad RXTE energy band between 2 and 13 keV. A preliminary test of any energy dependence, repeating the analysis for the 2-7 and 7-13 keV bands, did not reveal any clear effect of the choice of energy band. However, that is not particularly surprising considering the range of energies; given the RXTE response, both energy bands are dominated by the power law emission. In addition, the main feature of reflection, which might be related to the harmonic origin, in the RXTE band -the iron line -is located right at the division of the two bands.
Instead of analysing the light curves in different energy bands, the phase difference as a function of energy can alternatively be calculated using simply the energy dependence of both the fundamental and harmonic phase lag (see the introduction to the method in Ingram et al. 2017). Given RXTE's energy range and response, this approach is again not expected to show clear differences for RXTE observations. However, more novel observatories with different energy responses might reveal such an energy dependence; for instance, the recently employed Neutron star Interior Composition ExploreR (NICER) aboard the International Space Station is an Xray spectral-timing instrument with a soft response (peaking below 2 keV) and excellent timing capabilities. Alternatively, the Nuclear Spectroscopic Telescope ARray (NuSTAR) can detect photons up to 79 keV and has been used previously to measure phase differences between the QPO fundamental and harmonic . Observing campaigns with these observatories tracking the QPO evolution and probing different parts of the reflection spectrum (beyond the reach of RXTE PCA), can further test whether multiple mechanisms underly the Type-C QPO harmonic.

C O N C L U S I O N S
We have presented the first systematic analysis of the phase difference between the fundamental and harmonic of the QPOs from X-ray binaries. This phase difference gives important information about the physical processes that underlie the QPO emission mechanism, which is still unknown. We find that most of the studied X-ray binaries show a well-defined phase difference for a given QPO frequency and type, implying that they have welldefined waveforms. For Type-C QPOs, we find that the phase difference decreases with QPO (fundamental) frequency except for Swift J1753.5-0127 and XTE J1650-500. For Type-B QPOs, we find a completely different behaviour with the phase difference in general being constant with frequency, at around /π ∼ 0.55-0.6. Exceptions are found in H1743-322 and XTE J1817-330, which show a wide variety in phase differences.
A detailed analysis of the phase difference as a function of inclination does not yield a unique relation, although the data are suggestive that inclination plays a role in the consistency of the relation that is observed. Type-C QPOs in high-inclination sources show decreasing phase difference with frequency, with most sources clustered around a similar relation. Low-inclination sources show both increasing and decreasing phase difference versus frequency relations. For Type-B QPOs, the behaviour is similar regardless of the inclination of the source. The different behaviours for Type-B and Type-C QPOs lead to a new, quantitative distinction between Type-B and C QPOs, improving on the more qualitative classification based on the overall power spectral shape and broadband noise strength. For the sources analysed in this study, we draw the conclusion that when /π > 0.5 rad and ν > 4 Hz, the QPO should not be classified as a Type-C. Using this new scheme, we have found two incorrectly classified QPOs in our sample, that were originally classified as Type-C QPOs by Motta et al. (2015) and van den Eijnden et al. (2017).
To interpret the observed phase differences of the Type-B QPO, we propose a jet toy-model: the constant Type-B QPO phase difference observed in most sources can be explained by assuming that the QPO originates from the precession of an optically thick, X-ray emitting jet (Sriram et al. 2016;Stevens & Uttley 2016); the harmonic feature simply arises from the symmetric but nonsinusoidal nature of the resulting fundamental waveform.
For the Type-C QPO phase difference, we observe a more scattered dependence on QPO frequency for low-inclination sources. This larger scatter can be explained in the framework of geometric (Lense-Thirring) precession models. Such models, through the presence of both a direct and reflected component to the QPO waveform, might also be able to account for the evolution of the phase difference. However, detailed modelling is required to test our explanations of both the inclination-dependent scatter in, and the evolution of, the Type-C QPO phase difference. From the observational side, monitoring of the QPO with the X-ray observatories NICER and NuSTAR could further our understanding by searching for any energy dependence in the shape of the QPO waveform. Table B1. Initial guesses and lower and upper bounds for the fit parameters to models 1, 2, and 3. ν fund, 0 is the initial QPO fundamental frequency guess. 1.9 · ν fund, 0 2 · ν fund, 0 2.1 · ν fund, 0 K 5 [rms 2 /Hz] 0 10 −2 ∞ σ 5 [Hz] 0 1 10 ν center, 5 [Hz] 0.2 · ν fund, 0 0.5 · ν fund, 0 0.8 · ν fund, 0 C [rms 2 /Hz] 0 10 −3 max power short for the initial QPO fundamental frequency guess. C represents the constant power level that is added to the sum of Lorentzians. For model 1 only the BBN parameters (K 1 , σ 1 , K 2 , σ 2 ) and the parameters for Lorentzian for the fundamental peak (K 3 , σ 3 , ν center, 3 ) are taken into account. For models 2 and 3, another Lorentzian is added each time, and the bounds on these Lorentzian parameters can be found in K 4 , σ 4 , ν center, 4 and K 5 , σ 5 , ν center, 5 .

A P P E N D I X C : I N C L I NAT I O N D E P E N D E N C E O F T H E P H A S E D I F F E R E N C E
Tables C1 and C2 show the results of the straight line fits to Figs. 5 and 6 and inclination for the Type-B and C QPOs, respectively. The table shows number of data points used for the fit (data points), the slope of the fit, the fit intercept with the y-axis, the inclination of source, the inclination sample the source belongs to, and the reference for the inclination data. The errorbars on the slope and intercept represent a 1σ interval. For some sources, only two data points were available, so the straight line fit has zero error. Sources with only one data point aren't included in this table.  Motta et al. (2015). References for the specific inclinations: