<figure><img src="image_1.png" /></figure>

# A Systematic Spectral-timing Analysis of Kilohertz Quasi-periodic Oscillations in the Rossi X-Ray Timing Explorer Archive

Jon S. Troyer<fnref n="1" />, Edward M. Cackett<fnref n="1" />, Philippe Peille<fnref n="2" />, and Didier Barret<fnref n="3" />,<fnref n="4" />

Received 2018 February 12; revised 2018 May 4; accepted 2018 May 10; published 2018 June 25

## Abstract

Kilohertz quasi-periodic oscillations or kHz QPOs occur on the orbital timescale of the inner accretion flow, and may carry signatures of the physics of strong gravity ($c^2 \approx GM/R$) and possibly clues to constraining the neutron star equation of state. Both the timing behavior of kHz QPOs and the time-averaged spectra of these systems have been studied extensively, yet no model completely describes all the properties of kHz QPOs. Here, we present a systematic study of spectral-timing products of kHz QPOs from low-mass X-ray binary systems through the use of archival *Rossi X-ray Timing Explorer/Proportional Counter Array* data. For the lower kHz QPOs in 14 objects and the upper kHz QPOs in six objects, we were able to obtain correlated time lags as a function of QPO frequency and energy, as well as energy-dependent covariance spectra and intrinsic coherence. For the lower kHz QPOs, we find a monotonic decrease in lags with increasing energy, rising covariance to ~12 keV, and near unity coherence at all energies. For the upper kHz QPOs, we find near zero lags, rising covariance to ~12 keV, and less well-constrained coherence at all energies. These results suggest that while kHz QPOs are likely produced by similar mechanisms across the population of LMXBs, the lower kHz QPOs are likely produced by a different mechanism than the upper kHz QPOs.

**Key words:** accretion, accretion disks – stars: neutron – X-rays: binaries

## 1. Introduction

Kilohertz quasi-periodic oscillations or kHz QPOs (see e.g., van der Klis 1998) are the highest-frequency oscillations observed in the X-ray power spectrum of neutron star low-mass X-ray binary systems (LMXBs), and were first detected in 1995 after the launch of NASA's *Rossi X-ray Timing Explorer* (*RXTE*) (Bradt et al. 1993). The presence of two simultaneous kHz QPOs is common in many neutron star LMXBs, and has led to the classification of kHz QPOs as either lower or upper, depending on their frequency (see e.g., van der Klis et al. 1997). The short timescale of kHz QPOs is consistent with the dynamical timescales of the inner accretion flow around neutron stars (≈100 µs). Thus, the physical origins of kHz QPOs are of great interest, because they offer a possible window to the study of strong gravity ($GM/R \approx c^2$) and constraints on the physics of neutron stars (see e.g., Psaltis 2008).

Timing analysis of neutron star LMXBs has mostly been phenomenological, characterizing QPO properties through fitting Lorentzians, and observing how those properties change with source state. Many phenomenological correlations have been for kHz QPOs. These include various dependences of the properties of all types of QPOs (see e.g., Méndez et al. 1999; van Straaten et al. 2000; Di Salvo et al. 2001, 2003; Altamirano et al. 2008) on the position of the source on a color–color diagram (see e.g., Hasinger & van der Klis 1989), which describes the state of the source. Any model that we hope to apply to kHz QPOs must explain in a global sense, the many correlations between kHz QPO frequency and other lower-frequency oscillations, as well as the correlations between source state and kHz QPO appearance and properties. While many models of kHz QPOs have been developed over the years

since their discovery (see van der Klis 2000, 2006 for a review), none satisfactorily explains all the properties of kHz QPOs.

Studies involving analyzing the spectral and timing properties of kHz QPOs together can provide a fuller picture of their properties. However, until relatively recently, timing and spectral analysis approaches have largely been separate, and spectral analysis of neutron star LMXBs has usually focused on fitting physical models to the time-averaged energy spectrum, including the disk, boundary layer (BL), corona etc. Performing spectral-timing analyses allows one to connect which components of the energy spectrum are related to the QPO, and how those components are connected to each other. For instance, this can be done through studying the energy dependence of the variability with a root mean square (rms) spectrum (or similarly a covariance spectrum). Moreover, time lags between the variability at different energies provide further constraints on the emission processes involved. Both the covariance spectrum and time lags can be isolated for just the Fourier frequencies covered by the kHz QPOs (see Uttley et al. 2014 for a review of spectral-timing techniques).

Previous studies have shown that the fractional rms amplitude of kHz QPOs increases steadily from ~3 keV to 12 keV in nearly all sources (see e.g., Berger et al. 1996; Zhang et al. 1996; Wijnands et al. 1997; Méndez et al. 2001), suggesting a source with thermal Comptonization rather than soft, blackbody emission. Additionally, Gilfanov et al. (2003) and Revnivtsev & Gilfanov (2006) showed that the rms energy spectrum of neutron star LMXBs are consistent with a thermal Comptonization spectrum, which they suggest is due to the neutron star BL. Conservative arguments place the emission power of the BL as at least half of the total accretion powered emission (e.g., Sunyaev & Shakura 1986; see Gilfanov &

<footnote n="1">¹Department of Physics & Astronomy, Wayne State University, 666 W. Hancock Street, Detroit, MI 48201, USA; jon.troyer@wayne.edu</footnote>

<footnote n="2">²Centre National d'Études Spatiales, Centre spatial de Toulouse, 18 Av, Edouard Belin, F-31401 Toulouse cedex 9, France</footnote>

<footnote n="3">³Université de Toulouse; UPS-OMP; IRAP; Toulouse, France</footnote>

<footnote n="4">⁴CNRS; Institut de Recherche en Astrophysique et Planétologie; 9 Av, colonel Roche, BP 44346, F-31028 Toulouse cedex 4, France</footnote>
Sunyaev 2014 for a review of neutron star BL physics). Specifically, Gilfanov et al. (2003) showed that the rms energy spectrum of kHz QPOs in 4U 1608−52 is consistent with Comptonized emission. Additionally, Peille et al. (2015) (for 4U 1608−52 and 4U 1728−34) and Troyer & Cackett (2017) (for Aql X-1) both showed that the energy-dependent covariance spectra (see Section 4.2) are also consistent with Comptonized emission. Thus, kHz QPOs are most likely associated with thermal Comptonization.

Several physical scenarios can lead to energy-dependent time lags. For instance, where thermal accretion disk photons are Compton up-scattered by a hot corona or BL to higher energies, delays between variations of higher-energy photons and variations of lower-energy photons are expected to arise. This is because higher-energy photons must have undergone more Compton scatterings, and hence spend a longer time within the corona compared with lower-energy photons, which undergo fewer scatterings. A detailed analysis of Comptonization models suggest that this delay should scale with $\sim \log \frac{E_{high}}{E_{low}}$ (Nowak & Vaughan 1996), which implies that higher-energy photon oscillations occur after lower-energy oscillations. See Lee & Miller (1998), Lee et al. (2001), Kumar & Misra (2014), and Kumar & Misra (2016) for detailed discussions of Comptonization lags in QPOs. The study of kHz QPOs has shown that in some instances, high-energy oscillations lead the low-energy oscillations, and in some instances that the reverse is true (see e.g., Vaughan et al. 1998; Kaaret et al. 1999; Barret 2013; de Avellar et al. 2013; Peille et al. 2015; Troyer & Cackett 2017).

Another process that could lead to time lags is reverberation (see Uttley et al. 2014 for a review of X-ray reverberation in accreting black holes). Neutron star LMXBs often show prominent broad Fe Kα emission lines (see, e.g., Cackett et al. 2010). These lines are produced by irradiation of the accretion disk by high-energy photons from the corona and/or BL and are reflected and reprocessed by the accretion disk. Because photons that arrive directly to the observer from the corona/boundary layer travel a shorter distance than those reflected off the accretion disk, a time lag should be present between the direct and reflected emission.

The first time lags in kHz QPOs were observed by Vaughan et al. (1998) in 4U 1608−52, and Kaaret et al. (1999) in 4U 1636−53, and used 3–4 energy bins. More recently, spectral-timing analysis has been employed using significantly more energy bins, and detailed studies of four neutron star LMXBs have been performed: 4U 1608-52 (Barret 2013; de Avellar et al. 2013), 4U 1636−53 (de Avellar et al. 2013, 2016), 4U 1728−34 (Peille et al. 2015), and Aql X-1 (Troyer & Cackett 2017). These studies have shown that for the lower kHz QPOs, the higher-energy photons systematically arrive first. The magnitude of these broadband lags have all been on the order of the light-travel time from the neutron star/BL to the inner accretion disk.

In cases where the upper kHz QPO is detected, the broadband lags are generally consistent with zero, or they exhibit a slight increase with energy (see e.g., de Avellar et al. 2013; Peille et al. 2015). The interpretation of time lags is unclear, and the complete picture of the source of time lags is likely complex, involving a combination of several physical processes as previously stated. Understanding the ubiquity of time lags among different neutron star sources, as well as additionally complementary relationships between time lags

and other spectral-timing products, may lead to further insights regarding their origin.

To date, no comprehensive study exists that seeks to compute spectral-timing products of a large collection of neutron star LMXB systems. In this paper, we present the first systematic study of neutron star LMXB systems for data in the RXTE Proportional Counter Array (PCA) archive that looks at (a) time lags as a function of QPO frequency and energy; (b) intrinsic coherence; and (c) fractional rms and covariance in relative rms units.

In Section 2, we describe the data reduction and kHz QPO search methodology. In Section 3, we discuss how we classify a detected QPO as either lower or upper. In Section 4, we detail the computational process and results of various spectral-timing products. In Section 5, we discuss the implications and interpretation of our results, as well as some of the complementary nature of the various spectral-timing results. Finally, in Section 6, we summarize our results and state any conclusions that might be drawn from this work.

## 2. Data Reduction

We searched the entire RXTE/PCA archive for event mode observations of neutron star LMXB systems. In order to provide meaningful and consistent spectral-timing products, we considered only observations where the full 64 channel RXTE/PCA energy bands were available. We also required observations with as good or better than 125 µs time resolution. The list of neutron star LMXB sources we considered, as well as the nature of the observational limitations present in the RXTE/PCA archive, are shown in Table 1. In a small number of cases, objects present in the archive did not have event mode data available.

Once observations meeting the above criteria were identified, we created their associated response matrices using the FTOOL pcarsp to determine the channel-to-energy conversion table. Barycentric correction was then performed using the FTOOL faxbary. We applied good time intervals (GTI) to avoid the Earth's limb, South Atlantic Anomaly (SAA), and to remove any thermo-nuclear bursts that occurred.

We computed the discrete Fourier transform (DFT) over 4 s intervals, with a Nyquist frequency of 2048 Hz, and averaged 256 of these 4 s DFTs to produce a 1024 s average spectra. This average time was adopted to facilitate easier detection of the upper kHz QPOs. We binned the average power spectrum in 1 Hz bins. The power spectrum was initially computed over the 3–30 keV energy range. We searched the power spectra for excess power between 300 and 1300 Hz, and fit a single Lorentzian profile using the maximum likelihood estimator (MLE) method as described in Barret & Vaughan (2012). We attempted a fit with a double Lorentzian profile in cases where a second excess power was found. In this way, each 1024 s power spectrum was characterized by one, two, or no candidate kHz QPOs. Only individual 1024 s spectra where the QPO was significantly detected were considered for further analysis. We consider a QPO as significant if the ratio of the fitted Lorentzian normalization ($I_{\text{lor}}$) to the error in the normalization($\Delta I_{\text{lor}}$) is: $I_{\text{lor}}/\Delta I_{\text{lor}} \ge 3.0$ (e.g., Boutelier et al. 2010). A significant QPO spectrum is a 1024 s power spectrum in which a kHz QPO is significantly detected.

The initial detections of kHz QPOs in the sources we analyze are well documented in van der Klis (2000) with the following exceptions: Kaaret et al. (2007) for XTE J1739−285,
<table><caption>**Table 1**
Neutron Star Low-mass X-Ray Binaries Considered for Analysis</caption><thead><tr><th>Object ID</th><th>Source Classification</th><th>Event Mode Obs. Present?</th><th>Full Energy Range Available?</th><th>kHz QPOs Detected?ª</th></tr></thead><tbody><tr><td>4U 1608-52</td><td>Atoll</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>4U 1636-53</td><td>Atoll</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>4U 0614+09</td><td>Atoll</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>4U 1728-34</td><td>Atoll</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>4U 1702-43</td><td>Atoll</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>4U 1820-30</td><td>Atoll</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>Aql X-1</td><td>Atoll</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>4U 1735-44</td><td>Atoll</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>XTE J1739-285</td><td>Burster</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>EXO 1745-248</td><td>Burster</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>4U 1705-44</td><td>Atoll</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>SAXJ1748.9-2021</td><td>Atoll</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>IGR J17191-2821</td><td>Atoll</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>4U 1915-05</td><td>Atoll</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>SAXJ1808.4-3658</td><td>AMP</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>SAXJ1750.8-2900</td><td>Atoll</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>RXJ1709.5-2639</td><td>Other</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>XTE J2123-058</td><td>Atoll</td><td>Yes</td><td>Yes</td><td>No</td></tr><tr><td>KS 1731-260</td><td>Burster</td><td>Yes</td><td>Yes</td><td>No</td></tr><tr><td>Cyg X-2</td><td>Z</td><td>Yes</td><td>Yes</td><td>No</td></tr><tr><td>EXO 0748-676</td><td>Transient</td><td>Yes</td><td>Yes</td><td>No</td></tr><tr><td>MXB 1659-298</td><td>Transient</td><td>Yes</td><td>Yes</td><td>No</td></tr><tr><td>XTE J1723-376</td><td>Burster</td><td>Yes</td><td>Yes</td><td>No</td></tr><tr><td>2S 0918-549</td><td>Burster</td><td>Yes</td><td>Yes</td><td>No</td></tr><tr><td>1A 1246-588</td><td>Burster</td><td>Yes</td><td>Yes</td><td>No</td></tr><tr><td>2A 1822-371</td><td>AMP</td><td>Yes</td><td>Yes</td><td>No</td></tr><tr><td>1A 1744-361</td><td>Atoll</td><td>Yes</td><td>Yes</td><td>No</td></tr><tr><td>MXB 1730-33</td><td>Burster</td><td>Yes</td><td>Yes</td><td>No</td></tr><tr><td>GRO J1744-28</td><td>AMP</td><td>Yes</td><td>Yes</td><td>No</td></tr><tr><td>Sco X-1</td><td>Z</td><td>Yes</td><td>No</td><td>No</td></tr><tr><td>GX 5-1</td><td>Z</td><td>Yes</td><td>No</td><td>No</td></tr><tr><td>GX 17+2</td><td>Z</td><td>Yes</td><td>No</td><td>No</td></tr><tr><td>GX 340+0</td><td>Z</td><td>Yes</td><td>No</td><td>No</td></tr><tr><td>GX 349+2</td><td>Z</td><td>Yes</td><td>No</td><td>No</td></tr><tr><td>XTE J1751-305</td><td>AMP</td><td>Yes</td><td>No</td><td>No</td></tr><tr><td>XTE J1807-294</td><td>AMP</td><td>Yes</td><td>No</td><td>No</td></tr><tr><td>XTE J1814-338</td><td>AMP</td><td>Yes</td><td>No</td><td>No</td></tr><tr><td>2S 1803-245</td><td>Atoll</td><td>Yes</td><td>No</td><td>No</td></tr><tr><td>IGR J18245-2452</td><td>AMP</td><td>Yes</td><td>No</td><td>No</td></tr><tr><td>Ser X-1</td><td>Atoll</td><td>Yes</td><td>No</td><td>No</td></tr><tr><td>XTE J1701-462</td><td>Z</td><td>Yes</td><td>No</td><td>No</td></tr><tr><td>XTE J0929-314</td><td>AMP</td><td>No</td><td>No</td><td>No</td></tr><tr><td>4U 1746-37</td><td>Atoll</td><td>No</td><td>No</td><td>No</td></tr></tbody></table>

Note. Source classification from Liu et al. (2007). While all objects considered are present in the RXTE/PCA archive, additional columns indicate the presence of observation files in the event mode necessary for analysis. Additionally, where event mode files are present, we indicate when these have the full RXTE/PCA energy range available. Finally, we also show whether we were able to detect kHz QPOs from these observations in the manner described in Section 2. In some cases, while we were able to detect kHz QPOs, there was not sufficient detected QPO observation time to obtain spectral-timing results. Those sources shown in bold are the ones where we present analysis.

Mukherjee & Bhattacharyya (2011) for EXO 1745-248, Altamirano et al. (2010b) for SAXJ1748.9-2021, and Altamirano et al. (2010a) for IGR J17191-2821.

Within a single event mode file, each 1024 s spectra with a significantly detected QPO was shifted-and-added (see Méndez et al. 1998) to the mean QPO frequency. In cases where we detect both an upper and lower kHz QPO in a single 1024 s spectrum, the lowest-frequency QPO was used for shifting and adding. In all cases, the QPO frequency of each individual 1024 s spectrum within an event file was within 25 Hz of the mean QPO frequency. This narrow frequency range demonstrates that we are not mixing lower and upper kHz QPOs with

this shift-and-add procedure. The resulting spectrum was again fit with a single or double Lorentzian profile using the MLE. We used these profiles to classify the QPOs of each event file. This provided cleaner profiles with higher signal-to-noise-ratios (S/Ns) over individual spectra, while still maintaining temporal continuity of the observation.

kHz QPOs are not detected in all sources. We emphasize that we only searched event mode observations with the previously mentioned constraints for kHz QPOs, and attempted no detailed reconstruction of the QPO. Additionally, we only performed single detections on 1024 s observation segments, while other studies might average over longer periods of time

<footnote>ª kHz QPOs detected means detection using the analysis method in this paper and does not imply that there are not detectable kHz QPOs in these sources.</footnote>
to detect kHz QPOs. Finally, in some cases, kHz QPOs were detected in only a small number of averaged spectra and thus too few data points were available for meaningful spectral-timing analysis.

### 3. kHz QPO Classification

We wish to partition the analysis products by QPO type (lower or upper kHz QPO). Therefore, following Peille et al. (2015) in the case where a single kHz QPO is significantly detected in an event file, we classify that QPO using the quality factor (Q) parameter. Quality factor is defined as the ratio of the QPO frequency to its FWHM. In cases where two kHz QPOs are simultaneously detected, the quality factor is not needed, as the classification is obvious. For a single kHz QPO detection, we look to the relationship between quality factor and QPO frequency. The distinction between lower and upper kHz QPO can be generally inferred from this relationship. Barret et al. (2006) demonstrated that the lower and upper kHz QPOs follow distinct tracks on this diagram. Using this relationship for each source, we set quality factor limits to classify single kHz QPO detections. We compute the quality factor of the QPOs detected in each event file. These plots are shown in Figure 1. Looking at each object individually, we set frequency ranges for lower and upper kHz QPOs. We also conservatively set allowed ranges of Q to ensure unambiguous classification.

Once the kHz QPOs were classified, we were able to organize the event files for each object in order to combine observations with lower kHz QPOs and observations with upper kHz QPOs. These event files were used in the computation of various spectral-timing products described in Section 4.

### 4. Spectral-timing Products

For each object, we combined all the event files for each QPO type across the entire RXTE/PCA archive. This is necessary because in many cases the QPO is only weakly detected, and spectral-timing analysis requires a significant number of counts in order to constrain the results. To allow meaningful comparisons from source to source, we performed an identical analysis on all objects. Having identified all kHz QPOs present in the data, we next produce a range of spectral-timing products of interest. See Nowak et al. (1999) and Uttley et al. (2014) for a review of spectral-timing analysis methods. We calculate all spectral-timing products following the prescription given in Uttley et al. (2014).

The detector response varied over the multi-year range of data we averaged. Therefore, it was necessary to rebin the channel-to-energy conversion tables to accommodate both small fluctuations in the energy boundaries of observations with identical channels, and allow averaging of observations with different channel-to-energy conversion matrices. All spectra were normalized to the total count rate using fractional rms normalization (Miyamoto et al. 1991). We subtracted Fourier amplitudes where no source signal is nominally present (1350–1700 Hz) from the cross spectra in order to correct for dead time induced cross-talk (van der Klis et al. 1987).

For the frequency-dependent lags, we compute the cross spectrum between the 3–8 keV band and the 8–20 keV for each 1024 s spectrum with a significantly detected kHz QPO, averaging over the FWHM of the QPO profile. In this way,

we obtain a lag measurement as a function of kHz QPO frequency.

For the energy-dependent spectral-timing products, we shifted-and-added (Méndez et al. 1998) both the power spectra and the cross spectra of the reference band and the CI to the mean frequency of the classified QPO frequency range of each object. We did this separately for the lower and upper kHz QPOs, and only for spectra where a QPO was detected significantly as defined in Section 2. After shifting all the spectra, we fit a single Lorentzian to the shift-and-added spectra in order to characterize the FWHM of the shift-and-added QPO. Averaging across one FWHM, we then computed the intrinsic coherence (Vaughan & Nowak 1997), covariance (Wilkinson & Uttley 2009), as well as the fractional rms, and energy-dependent lags (Nowak & Vaughan 1996; Vaughan et al. 1998). See Uttley et al. (2014) for details of each calculation. The details of the intrinsic coherence results are shown in Section 4.1, the rms and covariance results in Section 4.2, a look at the frequency dependence and energy dependence of the lags are shown in Sections 4.3 and 4.4, respectively. Unless otherwise stated, all errors reported are at the $1\sigma$ level.

#### 4.1. Intrinsic Coherence

The intrinsic coherence provides a quantitative measure of how one signal is related to another by linear transformation (e.g., Vaughan & Nowak 1997). Generally speaking, if we wish to interpret the various spectral-timing products, especially time lags, as a result of some type of physical process, i.e., energy-dependent Comptonization delays or light-travel time delays due to spatial separation of one or more emitting regions, it is crucial that the underlying signals have near unity coherence. In other words, if intrinsic coherence is not near unity, we lose the underlying attachment of the two signals to a direct physical process. Thus, understanding the coherence of the kHz QPOs across the energy range of interest is an important place to start. Both Vaughan et al. (1998) and Kaaret et al. (1999) showed that kHz QPOs in 4U 1608–52 and 4U 1636–53 have high (near unity) coherence. Additionally, de Avellar et al. (2013) nicely showed the frequency dependence of the coherence across the kHz QPO profile. In Figure 1 of that paper, the unity coherence of the QPOs is illustrated.

The average coherence as a function of energy for the lower and upper kHz QPOs are shown in Figures 2 and 3, respectively. Here, we adopt confidence limits for the high-signal, high-coherence regime as defined in Nowak & Vaughan (1996) for all values of coherence due to the difficulties in quantifying the coherence confidence limits outside of this regime (e.g., Vaughan & Nowak 1997). Thus, confidence limits for low values of coherence should be considered unreliable.

General speaking, we find that both kHz QPOs show near unity coherence across a majority of the energy range we use. Objects with larger amounts of data generally show more well-constrained coherence. Additionally, at higher energies, the QPO signal can become dominated by noise, seen most convincingly in 4U 0614+09, which can reduce the coherence.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Quality factor vs. kHz QPO frequency for the 14 objects studied. Quality factor is the ratio of the QPO frequency to the FWHM ($\nu/\Delta\nu$). Each point represents a single shift-and-added event mode file using 1024 s spectra with a significantly detected kHz QPO. The different quality factor tracks (see e.g., Barret et al. 2006) are used to classify the QPOs as lower or upper in cases where the QPOs are ambiguous. Limits in quality factor ($Q$) and frequency ($\nu$) were adopted from a visual-inspection of the data in the case of a single kHz QPO detection. In all cases, a minimum quality factor of five was employed to ensure a good detection. No maximum $Q$ was used to filter detections of lower kHz QPOs. The vertical red dashed lines indicate the frequency boundaries of the lower kHz QPOs. The vertical blue dashed lines indicate the frequency boundaries of the upper kHz QPOs. The horizontal red dashed line indicates the minimum quality factor for lower kHz QPOs. The horizontal blue dashed line indicates the maximum quality factor for upper kHz QPOs.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Average intrinsic coherence as a function of energy for the lower kHz QPOs. Intrinsic coherence measures the degree to which two signals are related by linear transformation. In this case, one signal is the CI and the other is the reference band. In general, objects with larger numbers of observations tend to have a more well-constrained coherence. This is likely simply due to more signal being available.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Average intrinsic coherence as a function of energy for the upper kHz QPOs. See the caption in Figure 2 for additional details.</figcaption></figure>

## 4.2. Covariance

The covariance measures the variability power in an energy bin that is correlated with the reference band (see Wilkinson & Uttley 2009; Uttley et al. 2011, 2014). For high-coherence signals, the covariance is essentially the rms spectrum with intrinsically higher S/N (Uttley et al. 2014). It is well known that the rms spectrum of kHz QPOs increases with increasing energy up to about ≈15 keV, where it levels off (see e.g., Berger et al. 1996; Zhang et al. 1996; Wijnands et al. 1997; Méndez et al. 2001). Peille et al. (2015) showed similar behavior in the fractional covariance for 4U 1728−34, and this trend was recently extended below 3 keV through NICER observations of an upper kHz QPO in 4U 0614+09 (Bult et al. 2018). Our results are generally consistent with this behavior of increasing fractional covariance with energy, flattening off above 15 keV.

In order to illustrate the usefulness and limitations of the covariance statistic, we have computed both the rms and covariance in fractional rms units. These data for the lower and upper kHz QPOs are shown in Figures 4 and 5, respectively. The errors on the rms and covariance calculation are valid in the limit where the intrinsic coherence is unity (Uttley et al. 2014).

In general, we note increasing rms and covariance with increasing energy for both the lower and upper kHz QPOs. In a few cases, a drop in rms and covariance is seen in the highest-energy bin. In 4U 0614+09, we note decreasing rms and covariance at energies above ~15 keV. We also note that there is sometimes a divergence between the rms and covariance. These are both likely due to a drop in coherence seen at those energies. For example, see Figure 2. At higher energies, the QPO signal can become dominated by noise, which can result in the reduction of both the rms and covariance.

## 4.3. Lags versus Frequency

To establish the presence of any lags, and explore any possible frequency dependence, we computed lags between two broad-energy bins: 3–8 keV and 8–20 keV. We note here that the broad-energy lags are not calculated using the channel-of-interest and reference band convention. These lags are produced directly between the two broad-energy bins. We considered lower kHz QPOs separately from uppers. We were able to derive lower kHz QPO lag versus frequency relationships for 14 objects and upper kHz QPO lag versus frequency relationships for six objects. The lag versus frequency relationships for the lower kHz QPOs and the upper kHz QPOs are shown in Figures 6 and 7, respectively. A positive lag indicates that the higher-energy band photon variability arrives before the lower-energy band photon variability. We binned the lag/frequency data into a maximum of eight bins in an effort to highlight any frequency dependence.

In general, the lags for the lower kHz QPOs show positive average lags, i.e., so-called “soft lags” where the lower-energy photons lag the higher-energy ones. This result is seen in all previously published work where such lags have been computed (Barret 2013; de Avellar et al. 2013, 2016; Peille et al. 2015; Troyer & Cackett 2017). The upper kHz QPOs generally have reduced data quality. In general, we find lags consistent with zero or slightly negative average lags, except in two cases. Zero or slightly negative average lags are also consistent with previously published work and illustrates that the nature of the lags is fundamentally different between the lower and upper kHz QPO (de Avellar et al. 2013, 2016; Peille et al. 2015; Troyer & Cackett 2017). Moreover, the values of the average lags are consistent with all comparable previous analyses which include: 4U 1608−52, 4U 1636−53, 4U 1728−34, and Aql X-1. Both 4U 1608−52 and 4U 1702−43 show
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Average fractional rms and covariance in fraction rms units as a function of energy for the lower kHz QPOs. Fractional rms is shown in red squares and covariance in black dots. In general, objects show increasing fractional variability with increasing energy. The divergence of fractional rms and covariance implies a drop in coherence, seen here at the lowest energy bins and the highest. This behavior is consistent with the behavior of the coherence seen in Figure 2. At energies above ~15 keV, the covariance falls off generally, while the fraction rms levels off or increases while becoming less well constrained. The drop in covariance at high energies generally tracks the drop in intrinsic coherence, which reduces the degree to which one can interpret the covariance as conveying the variability of the QPO. Poorly constrained rms values due to poor S/N are not shown.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Average fractional rms and covariance as a function of energy for the upper kHz QPOs. The fractional rms is shown in red squares and covariance in black dots. Objects generally show increasing variability with increasing energy, with the behavior of upper kHz QPOs similar to that of the lower kHz QPOs</figcaption></figure>

positive average lags for the upper kHz QPOs. In both cases, there are relatively few data and the upper kHz QPOs are very weakly detected as shown by their poorly constrained coherence.

While in some cases significant upper kHz QPOs are found in event file shift-and-added spectra shown in Figure 1, an individual 1024 s spectrum with a significant upper kHz QPO may not exist and thus not appear in the lag versus frequency analysis.

In general, the dependence of lag on kHz QPO frequency is likely complex. Barret (2013) showed the binned frequency-dependent lags for 4U 1608−52 increase monotonically to ~700 Hz and then decrease monotonically at higher frequencies. de Avellar et al. (2013) showed no dependence of lags on frequency for 4U 1608−52 using a smaller data set than Barret (2013). Additionally, de Avellar et al. (2013) showed a slight decrease in lower kHz QPO lags for 4U 1636−53 above 850 Hz, with no frequency dependence for the upper kHz QPO. de Avellar et al. (2016), who computed phase lags, find no significant frequency dependence for either the lower or upper kHz QPOs in 4U 1636−53. When fitting the lag/frequency relationship, both de Avellar et al. (2013) and de Avellar et al. (2016) note that linear or quadratic fits are not favored over a constant fit. In 4U 1728−34, Peille et al. (2015) find no frequency dependence for the lower kHz QPO lags, and the possibility of a slight decrease of the lags at higher frequency for the upper kHz QPO. In Aql X-1, Troyer & Cackett (2017) find no significant frequency dependence on the lags for the lower kHz QPO.

In all previous works, various energy boundaries, power spectral averaging times, and/or data sets are used. Barret (2013) uses a 128 s averaged power spectra to calculate the

lags, and a smaller data set than this work. Both de Avellar et al. (2013, 2016) use different energy bin boundaries and also a smaller data set than this work as well as a 16 s power spectra average time. While Peille et al. (2015) and Troyer & Cackett (2017) use the same energy boundaries and similar data sets, Troyer & Cackett (2017) uses a 256 s power spectra average time.

Our use of a 1024 s spectra averaging time, while necessary to maximize the kHz QPO detections in many sources, could mask the underlying lag/frequency relationship. Thus, we choose to apply a simple linear fit of the unbinned data in an effort to characterize the lag-frequency relationship, while noting the underlying relationship could be more complex. While not perfect, this does allow for meaningful comparisons of any trends in the lag/frequency relationship between sources.

In order to determine if a linear trend is favored over a constant, we computed the $\chi^2$ statistic for both a linear and constant model of the data and then used an F-test to determine the significance of the statistically favored model. The average lag; best-fit slope and intercept; and $\chi^2$ values with the statistically favored model are shown in Table 2 for the lower kHz QPOs and Table 3 for the upper kHz QPOs. For the lower kHz QPOs in four objects: EXO 1745−248, 4U 1705−44, SAXJ1748.9-2021, and 4U 1915−05, and the upper kHz QPOs in two objects: 4U 1608−52 and 4U 1702−43; we assessed the number of data points (<5) too small to provide a meaningful $\chi^2$ comparison of the data.

For the lower kHz QPOs, in three objects, 4U 1608−52, 4U 1820−30, and Aql X-1, the linear model is favored with >3$\sigma$ significance over the constant model. In all other models, there is no statistically significant difference between a linear fit
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Broad-energy lags as a function of kHz QPO frequency for the lower kHz QPO. The low-energy band is 3–8 keV, and the high-energy band is 8–20 keV. Each small dot represents a single 1024 s spectrum with a significantly detected lower kHz QPO. The data are binned in a maximum of eight equally spaced frequency bins shown in large blue squares. The average lag with $1\sigma$ error is shown in red and the values of the average lags are shown in Table 2. A positive lag indicates that variability in the higher-energy band (8–20 keV) arrives before variability in the lower-energy band (3–8 keV). In general, lower kHz QPOs show positive average lags. In all cases, there is either a slight decrease in lags as a function of QPO frequency or no significant frequency dependence. In objects with more than five data points, the linear fit parameters are also shown in Table 2.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Broad-energy lags as a function of kHz QPO frequency for the upper kHz QPO. See the caption of Figure 6 for additional details. The average lag with $1\sigma$ error is shown in red and the values of the average lags are shown in Table 3. The upper kHz QPOs show lags which are consistent with 0, except in the cases where the number of 1024 s segments with lags is small, i.e., 4U 1608–52 and 4U 1702–43. In all cases, there is no significant dependence of the lag on frequency. The linear fit parameters are also shown in Table 3.</figcaption></figure>

<table><caption>**Table 2**
Lower kHz QPO Lag/Frequency Data</caption><thead><tr><th>Object ID</th><th>Average Lag (µs)</th><th>Slope (µs/Hz)</th><th>Intercept (µs)</th><th>χ<sup>2</sup> (dof) Constant</th><th>χ<sup>2</sup> (dof) Linear</th><th>Favored Model</th></tr></thead><tbody><tr><td>4U 1608–52</td><td>23.7 ± 1.3</td><td>-0.09 ± 0.02</td><td>96.3 ± 1.3</td><td>115.0(92)</td><td>97.9(91)</td><td>Linear 3.8σ</td></tr><tr><td>4U 1636–53</td><td>16.5 ± 0.8</td><td>-0.02 ± 0.02</td><td>36.0 ± 0.8</td><td>767.3(742)</td><td>765.3(741)</td><td>Linear 1.4σ</td></tr><tr><td>4U 0614+09</td><td>46.8 ± 4.0</td><td>0.07 ± 0.1</td><td>-0.8 ± 4.1</td><td>117.4(112)</td><td>116.9(111)</td><td>Linear 0.69σ</td></tr><tr><td>4U 1728–34</td><td>9.5 ± 1.6</td><td>0.01 ± 0.03</td><td>1.9 ± 1.6</td><td>157.4(175)</td><td>156.8(174)</td><td>Linear 0.81σ</td></tr><tr><td>4U 1702–43</td><td>35.6 ± 2.8</td><td>0.03 ± 0.04</td><td>14.7 ± 2.9</td><td>112.3(114)</td><td>111.9(113)</td><td>Constant 0.63σ</td></tr><tr><td>4U 1820–30</td><td>5.4 ± 3.1</td><td>-0.26 ± 0.07</td><td>193 ± 3</td><td>175.2(153)</td><td>162.0(152)</td><td>Linear 3.4σ</td></tr><tr><td>Aql X-1</td><td>29.2 ± 3.0</td><td>-0.15 ± 0.04</td><td>154 ± 3</td><td>96.2(88)</td><td>86.1(87)</td><td>Linear 3.1σ</td></tr><tr><td>4U 1735–44</td><td>6.9 ± 3.1</td><td>0.12 ± 0.08</td><td>-82.5 ± 3.2</td><td>108.9(109)</td><td>106.7(108)</td><td>Linear 1.5σ</td></tr><tr><td>XTE J1739–285</td><td>9.6 ± 9.3</td><td>-0.3 ± 0.2</td><td>267 ± 10</td><td>10.0(21)</td><td>8.3(20)</td><td>Linear 1.9σ</td></tr><tr><td>EXO 1745–248</td><td>18.9 ± 10.1</td><td>0.06 ± 0.43</td><td>-22.0 ± 11.9</td><td>0.2(3)</td><td>0.2(2)</td><td>...</td></tr><tr><td>4U 1705–44</td><td>-84.3 ± 35.0</td><td>-2.4±2.5</td><td>1875 ± 44</td><td>0.9(1)</td><td>0(0)</td><td>...</td></tr><tr><td>SAXJ1748.9–2021</td><td>26.6 ± 23.0</td><td>8.7 ± 4.8</td><td>-6842 ± 31</td><td>5.2(3)</td><td>1.7(2)</td><td>...</td></tr><tr><td>IGR J17191–2821</td><td>-1.4 ± 9.5</td><td>3.4 ± 2.5</td><td>-73 ± 10</td><td>11.4(18)</td><td>11.0(17)</td><td>Linear 0.77σ</td></tr><tr><td>4U 1915–05</td><td>27.2 ± 32.8</td><td>3.4 ± 2.5</td><td>-2670 ± 33</td><td>3.4(2)</td><td>1.5(1)</td><td>...</td></tr></tbody></table>

**Note.** This includes the weighted mean of the lags; the best-fit slope and intercept of the linear model we fit to the unbinned data; the $\chi^2(\text{dof})$ of both the linear model and the constant lag model; and the result of an F-test indicating which of these two models is preferred and the statistical significance of the conclusion. In the cases of EXO 1745–248, 4U 1705–44, SAXJ1748.9-2021, and 4U 1915–05, we assess the number of data points (<5) too small to provide meaningful $\chi^2(\text{dof})$ analysis.

and a constant fit. While the data do generally show lags that decrease slightly with frequency (though at low statistical significance), because of the paucity of high-quality data for most of the sources and the fact that splitting into multiple frequency bins spreads the number of photons available for analysis across multiple spectral-timing products, we chose to average over all frequencies for all energy-dependent analyses (coherence, covariance, and lag-energy spectra). This allows for increased S/N when averaging all data together instead of

splitting into multiple bins, and enables a systematic method of analysis across all sources.

### 4.4. Energy-dependent Lags

Energy-dependent lags offer a way to begin describing the role of various time-dependent physical process involved in the creation of kHz QPOs. These data illustrate the time lag in each energy bin relative to the reference band. They should be read as relative lags, with the most negative lag being the energy
<table><caption>**Table 3**
Upper kHz QPO Lag/Frequency Data</caption><thead><tr><th>Object ID</th><th>Average Lag (μs)</th><th>Slope (μs/Hz)</th><th>Intercept (μs)</th><th>χ<sup>2</sup> (dof) Constant</th><th>χ<sup>2</sup> (dof) Linear</th><th>Favored Model</th></tr></thead><tbody><tr><td>4U 1608-52</td><td>38.2 ± 9.3</td><td>-0.13 ± 0.65</td><td>144 ± 9</td><td>3.3(4)</td><td>3.2(3)</td><td>...</td></tr><tr><td>4U 1636-53</td><td>2.6 ± 4.8</td><td>-0.01 ± 0.08</td><td>12.4 ± 4.8</td><td>81.9(103)</td><td>81.9(102)</td><td>Equal</td></tr><tr><td>4U 0614+09</td><td>6.8 ± 8.3</td><td>-0.10 ± 0.10</td><td>90.0 ± 8.5</td><td>18.1(32)</td><td>16.9(31)</td><td>Linear 1.4σ</td></tr><tr><td>4U 1728-34</td><td>-2.9 ± 2.8</td><td>-0.04 ± 0.04</td><td>30.0 ± 2.7</td><td>113.3(127)</td><td>112.1(126)</td><td>Linear 1.1σ</td></tr><tr><td>4U 1702-43</td><td>53.2 ± 19.6</td><td>0.10 ± 0.41</td><td>-25 ± 21</td><td>1.9(4)</td><td>1.8(3)</td><td>...</td></tr><tr><td>4U 1820-30</td><td>-11.4 ± 24.0</td><td>0.13 ± 0.16</td><td>-136 ±</td><td>3.7(8)</td><td>3.1(7)</td><td>Linear 1.1σ</td></tr></tbody></table>

Note. For additional details, see the caption of Table 2. For the upper kHz QPO, in the cases of 4U 1608–52 and 4U 1702–43, we assess the number of data points (<5) too small to provide meaningful χ²(dof) analysis.

from which the photons arrive first. We show lags as a function of energy for the lower and upper kHz QPOs in Figures 8 and 9, respectively. To further characterize the data, we performed linear fits and show the best-fit slope and intercept, as well as then number of significant QPO spectra we averaged to obtain the lag/energy spectrum in Tables 4 and 5 for the lower and upper kHz QPOs, respectively. The best-fit slope and intercept are used as a simple way to characterize the data. Additionally, the shape of the lag/energy spectrum depends on the strengths of the components of the energy spectrum, and thus would not necessarily fit a linear function.

For the lower kHz QPOs, the trend is generally monotonically decreasing lags with increasing energy. The lag of each energy bin channel of interest (CI) is measured with respect to the CI-subtracted reference band energy range (3–20 keV). This trend shows that the variability at lower energies lags the variability at higher energies, with the highest-energy photons arriving first. The upper kHz QPOs generally show lags that are flat, and in some cases show an increase in the lags at higher energies, albeit with poorer statistics than lower kHz QPO. This is consistent with the results of de Avellar et al. (2013) and Peille et al. (2015), which show different behavior of energy-dependent lags of the lower and upper kHz QPO for 4U 1608–52, 4U 1636–53, and 4U 1728–34. This implies that the physical mechanisms responsible for the creation of lower and upper kHz QPOs are likely different.

## 5. Discussion

We have performed a large-scale systematic study of the spectral-timing products of kHz QPOs for neutron star LMXBs using data in the entire RXTE/PCA archive. We examined the coherence, covariance, and lags in the lower kHz QPOs for 14 objects and in the upper kHz QPOs of six objects. This is a significant increase compared with all previous studies where lags in the lower kHz QPOs have been studied in four objects (Barret 2013; de Avellar et al. 2013, 2016; Peille et al. 2015; Troyer & Cackett 2017), and in the upper kHz QPOs for 3 objects. We also find similar results between the various sources where we were able to compute spectral-timing products. We summarize the main observational results:

1. both the lower kHz QPOs and upper kHz QPOs are highly coherent over most of the energy range studied (3–20 keV);

2. the fractional covariance of all kHz QPOs generally increases with energy and levels off around 15 keV;

3. three objects show a small decrease in lags with frequency between the 3–8 keV and 8–20 keV band for

the lower kHz QPO, but, otherwise no statistically significant trend is found; and

4. the lower kHz QPO lags generally show a decrease with energy, with the highest-energy photons arriving first, while the upper kHz QPO lags are generally consistent with zero lag.

We discuss the main findings in additional detail in the following sections.

### 5.1. Coherence

While generally consistent with unity, our results show the intrinsic coherence is not well constrained in objects with poor statistics across all energy bands. Generally, the extreme low- and high-energy bins of all objects lack statistics and thus also show poorly constrained intrinsic coherence.

The lower kHz QPOs have more and generally higher-quality data, and produce more well-constrained coherence. The confidence limits of the intrinsic coherence are difficult to quantify (see e.g., Nowak & Vaughan 1996; Vaughan & Nowak 1997; Uttley et al. 2014). The drop in intrinsic coherence in the highest-energy bin of a few sources and above ~15 keV in 4U 0614+09 is likely due to the loss of QPO signal at these low S/N energies.

### 5.2. Covariance

The covariance statistic is dependent on the coherence, thus we expect the covariance to drop when the coherence drops. In the regime of unity coherence, the covariance spectrum becomes equivalent to the rms spectrum with higher S/N, and we generally find good agreement between the covariance and rms spectra.

The fractional covariance of all kHz QPOs generally follow a similar trend, rising with increased energy and then leveling off around 15 keV, and in a few cases dropping in the highest-energy bin. An increase in fractional rms/covariance with energy is indicating that the fraction of kHz QPO signal that is variable is increasing with energy. This is consistent with the behavior of a Comptonized emission source. Gilfanov et al. (2003) showed that the rms spectrum of kHz QPOs was consistent with Comptonized emission. This idea is further supported by both Peille et al. (2015) and Troyer & Cackett (2017), who show the response-folded covariance spectra of kHz QPOs are consistent with Comptonized emission and do not require a disk component. Moreover, the first application of phase-resolved spectroscopy to the lower kHz QPOs in 4U 1608–52 also shows that it is the Comptonized component that is modulating (Stevens 2018).
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Energy-dependent lags for the lower kHz QPO. A positive value of lag in an energy bin indicates that the variability in that bin arrives after the variability in the reference band. In general, lower kHz QPOs show lags with a decreasing trend as energy increases. In order to further characterize the lags, we performed a linear fit in all cases. The best-fit line is shown in red. The best-fit intercept and slope with confidence limits are shown in Table 4.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Energy-dependent lags for the upper kHz QPO. In general, upper kHz QPOs show lags which are consistent with zero lag. In order to characterize the lags, we performed a linear fit in all cases. The best-fit line is shown in red. The best-fit intercept and slope with confidence limits are shown in Table 5.</figcaption></figure>

<table><caption>**Table 4**
Lower kHz QPO Lag/Energy Data</caption><thead><tr><th>Object ID</th><th>Number of Significant QPOs</th><th>Best-fit Slope (μs/keV)</th><th>Best-fit Intercept (μs)</th></tr></thead><tbody><tr><td>4U 1608-52</td><td>89</td><td>-4.9 ± 0.6</td><td>37.3 ± 2.6</td></tr><tr><td>4U 1636-53</td><td>742</td><td>-3.2 ± 0.3</td><td>24.3 ± 1.3</td></tr><tr><td>4U 0614+09</td><td>114</td><td>-8.1 ± 1.0</td><td>54.4 ± 4.9</td></tr><tr><td>4U 1728-34</td><td>174</td><td>-2.2 ± 0.4</td><td>20.3 ± 1.4</td></tr><tr><td>4U 1702-43</td><td>112</td><td>-6.4 ± 0.7</td><td>52.9 ± 3.2</td></tr><tr><td>4U 1820-30</td><td>147</td><td>-1.1 ± 0.4</td><td>8.9 ± 1.7</td></tr><tr><td>Aql X-1</td><td>92</td><td>-5.6 ± 0.7</td><td>42.5 ± 3.3</td></tr><tr><td>4U 1735-44</td><td>110</td><td>-0.9 ± 0.5</td><td>8.1 ± 1.8</td></tr><tr><td>XTE J1739-285</td><td>22</td><td>-0.5 ± 1.5</td><td>4.0 ± 5.9</td></tr><tr><td>EXO 1745-248</td><td>4</td><td>-1.5 ± 1.1</td><td>11.6 ± 4.1</td></tr><tr><td>4U 1705-44</td><td>2</td><td>6.5 ± 3.9</td><td>-54.5 ± 17.9</td></tr><tr><td>SAXJ1748.9-2021</td><td>4</td><td>-0.4 ± 3.6</td><td>-2.4 ± 12.4</td></tr><tr><td>IGR J17191-2821</td><td>18</td><td>-0.08 ± 1.41</td><td>2.7 ± 5.2</td></tr><tr><td>4U 1915-05</td><td>3</td><td>-4.0 ± 3.7</td><td>35.4 ± 16.7</td></tr></tbody></table>

The covariance is similar between the lower and upper kHz QPOs. This result is consistent with Gilfanov et al. (2003), where they show that by subtracting a putative modeled accretion disk from the average spectrum, the result looks like the response-folded rms spectrum for all the QPOs (lower, upper, and the 45 Hz hectohertz QPO) they considered. The fractional rms/covariance spectra are also very similar for different objects, suggesting the production mechanisms of kHz QPOs across different sources may be the same.

Our results also show the expected deviation between the covariance and rms spectra when the coherence departs from unity. In some cases, we note drops in both rms and covariance near 20 keV. The drop in rms and covariance in the highest-energy bin of a few sources and above ~15 keV in 4U 0614+09 is also likely due to the loss of QPO signal at these low S/N energies. The behavior of the covariance is similar in 4U 1728-34 to that described in Peille et al. (2015) using a similar aggregate data set. Mukherjee & Bhattacharyya (2012) note a drop in fractional rms above ~10 keV in 4U 1728-34 using a much smaller data set.

### 5.3. Lags

We analyzed the lags between the 3–8 keV and 8–20 keV light curves, determining both the average lag and looking for any frequency dependence. Three objects (4U 1608–52, 4U 1820–30, and Aql X-1) all show a decreasing trend of lags with increasing frequency for the lower kHz QPOs at >3σ significance. This result is generally consistent with and is well discussed in previous work (see e.g., de Avellar et al. 2013). All other objects (for both lower and upper kHz QPOs) have no statistical significant frequency dependence, though we note that the underlying frequency dependence of the lags may be more complex than linear.

The broad-energy lags for the lower kHz QPOs show correlated variations in the high-energy band occur on average before those in the low-energy band. The average lags are on the order of a few 10s of microseconds which is consistent with the putative light-travel timescale of the inner accretion flow of a neutron star. The average lags of the upper kHz QPOs are generally consistent with zero, and generally have fewer and lower quality data than the lower kHz QPOs.
<table><caption>**Table 5**
Upper kHz QPO Lag/Energy Data</caption><thead><tr><th>Object ID</th><th>Number of Significant QPOs</th><th>Best-fit Slope (μs/keV)</th><th>Best-fit Intercept (μs)</th></tr></thead><tbody><tr><td>4U 1608-52</td><td>8</td><td>-2.3 ± 1.2</td><td>15.0 ± 5.4</td></tr><tr><td>4U 1636-53</td><td>104</td><td>0.63 ± 0.57</td><td>-5.0 ± 2.0</td></tr><tr><td>4U 0614+09</td><td>114</td><td>-0.5 ± 1.1</td><td>2.3 ± 5.5</td></tr><tr><td>4U 1728-34</td><td>130</td><td>0.8 ± 0.3</td><td>-7.2 ± 1.3</td></tr><tr><td>4U 1702-43</td><td>112</td><td>-5.7 ± 2.3</td><td>44.8 ± 8.7</td></tr><tr><td>4U 1820-30</td><td>12</td><td>-1.1 ± 1.7</td><td>10.1 ± 7.0</td></tr></tbody></table>

We also study the lags as a function of energy, averaged over all frequencies. When comparing the energy-dependent lags between the lower and upper kHz QPOs, we notice different trends. The lower kHz lags generally decrease monotonically as energy increases. Those objects that do not show a strong trend typically have far fewer lower kHz QPOs detected. The exceptions are 4U 1735-44 and 4U 1820-30, which show little energy dependence, despite relatively large numbers of detections as seen in Table 4. Additionally, these objects have the smallest average broad-energy lags of all objects with well-constrained, near unity coherence.

The upper kHz QPO lags are not as well constrained, and are generally consistent with zero lag, except in the cases of 4U 1608-52 and 4U 1702-43 where the upper kHz QPOs are weakly detected. The objects with the best upper kHz QPO photon statistics, 4U-1636-53 and 4U 1728-34, show slightly increasing lags as a function of energy.

The energy dependence of the lags has the potential of constraining the geometry and physics of the emission region of kHz QPOs. For example, in the case of reverberation, the parts of the energy spectrum that are dominated by the reflected emission should show the largest lags behind the direct emission component. The energy-dependent lags in 4U 1608-52 were initially suggested to be due to reverberation by Barret (2013). However, models of the expected reverberation lags in 4U 1608-52 found the energy-dependent lags diverged from the modeled lags above ≈8 keV (Cackett 2016), eliminating reverberation as the sole source of emission for the lower kHz QPO. Moreover, a response-folded covariance spectra, produced for lower kHz QPOs, do not show an accretion disk component that would be expected from reverberation (Peille et al. 2015; Troyer & Cackett 2017). Both these studies found the seed photon temperature of the covariance spectra, which can be thought of as the spectrum of the kHz QPO, was systematically higher than the time-averaged spectra for 4U 1728-34 and Aql X-1, respectively. This is also consistent with the phase-resolved spectroscopy of 4U 1608-52, which shows phase lags between variations in the power-law index, high-energy cut-off temperature, and seed blackbody temperature (Stevens 2018).

Lee & Miller (1998), Lee et al. (2001), and more recently Kumar & Misra (2014, 2016) have applied Comptonization models to the lag/energy and rms spectra of kHz QPOs. Most of these studies use the rms spectrum of 4U 1608-52 from Berger et al. (1996) and the lag/energy spectrum from Vaughan et al. (1998). Kumar & Misra (2014, 2016) show that the lag/energy spectra of the lower kHz QPO can be recovered by a Comptonization model in which a fraction of the Comptonized emission impinges back on the seed photon source. This includes time lags where high-energy fluctuations lead lower-energy fluctuations.

## 6. Conclusion

We have analyzed all RXTE/PCA data for the objects listed in Table 1. We present a large-scale, systematic comparison of spectral-timing characteristics of kHz QPOs across many sources. We have shown that the average spectral-timing properties for all sources in the RXTE/PCA archive, which remains the best data set currently available for this kind of study, are very similar. This supports the idea that the production mechanism of kHz QPOs is the same for all sources. Additionally, we show that both the broadband lags and the lag/energy spectra for the lower and upper kHz QPOs are markedly different, which supports the idea that the emission mechanisms of the two types of kHz QPOs are different.

We computed the intrinsic coherence, rms/covariance, frequency-dependent and energy-dependent lags for the lower kHz QPO in 14 objects, with six of those yielding results for the upper kHz QPO. We have demonstrated only broad trends in our results. Both high timing and higher-energy resolution are required to model spectral-timing results in a robust fashion.

This work was only possible because of the existence of RXTE/PCA, which was uniquely capable of providing data necessary for spectral-timing work of kHz QPOs. While the energy resolution coupled with variations in the detector's response limits the energy resolution of analyses, we now understand the scope of the requirements needed to fully exploit this type of analysis. The current NICER mission can extend such analyses to lower energies than previously probed (Bult et al. 2018). Future proposed missions such as STROBE-X (Wilson-Hodge et al. 2017) and eXTP (Zhang et al. 2016), equipped with such capabilities will allow the full exploration of these analytic techniques.

We thank an anonymous referee for helpful comments that improved the clarity of the paper. J.S.T. and E.M.C. gratefully acknowledge support from the National Science Foundation through CAREER award number AST-1351222.

## ORCID iDs

Jon S. Troyer https://orcid.org/0000-0002-9191-219X
Edward M. Cackett https://orcid.org/0000-0002-8294-9281
Philippe Peille https://orcid.org/0000-0002-0695-1662
Didier Barret https://orcid.org/0000-0002-0393-9190

## References

* Altamirano, D., van der Klis, M., Méndez, M., et al. 2008, ApJ, 685, 436

* Altamirano, D., Linares, M., Watts, A., et al. 2010a, MNRAS, 401, 223

* Altamirano, D., Patruno, A., Watts, A., et al. 2010b, ATel, 2426, 1

* Barret, D. 2013, ApJ, 770, 9

* Barret, D., Olive, J.-F., & Miller, M. C. 2006, MNRAS, 370, 1140
Barret, D., & Vaughan, S. 2012, ApJ, 746, 131

Berger, M., van der Klis, M., van Paradijs, J., et al. 1996, ApJL, 469, L13

Boutelier, M., Barret, D., Lin, Y., & Török, G. 2010, MNRAS, 401, 1290

Bradt, H. V., Rothschild, R. E., & Swank, J. H. 1993, A&AS, 97, 355

Bult, P., Altamirano, D., Arzoumanian, Z., et al. 2018, ApJL, submitted

Cackett, E. M. 2016, ApJ, 826, 103

Cackett, E. M., Miller, J. M., Ballantyne, D. R., et al. 2010, ApJ, 720, 205

de Avellar, M. G. B., Méndez, M., Altamirano, D., Sanna, A., & Zhang, G. 2016, MNRAS, 461, 79

de Avellar, M. G. B., Méndez, M., Sanna, A., & Horvath, J. E. 2013, MNRAS, 433, 3453

Di Salvo, T., Méndez, M., & van der Klis, M. 2003, A&A, 406, 177

Di Salvo, T., Méndez, M., van der Klis, M., Ford, E., & Robba, N. R. 2001, ApJ, 546, 1107

Gilfanov, M., Revnivtsev, M., & Molkov, S. 2003, A&A, 410, 217

Gilfanov, M. R., & Sunyaev, R. A. 2014, PhyU, 57, 377

Hasinger, G., & van der Klis, M. 1989, A&A, 225, 79

Kaaret, P., Piraino, S., Ford, E. C., & Santangelo, A. 1999, ApJL, 514, L31

Kaaret, P., Prieskorn, Z., in't Zand, J. J. M., et al. 2007, ApJL, 657, L97

Kumar, N., & Misra, R. 2016, MNRAS, 461, 2580

Kumar, N., & Misra, R. 2014, MNRAS, 445, 2818

Lee, H. C., & Miller, G. S. 1998, MNRAS, 299, 479

Lee, H. C., Misra, R., & Taam, R. E. 2001, ApJL, 549, L229

Liu, Q. Z., van Paradijs, J., & van den Heuvel, E. P. J. 2007, A&A, 469, 807

Méndez, M., van der Klis, M., & Ford, E. C. 2001, ApJ, 561, 1016

Méndez, M., van der Klis, M., Ford, E. C., Wijnands, R., & van Paradijs, J. 1999, ApJL, 511, L49

Méndez, M., van der Klis, M., van Paradijs, J., et al. 1998, ApJL, 494, L65

Miyamoto, S., Kimura, K., Kitamoto, S., Dotani, T., & Ebisawa, K. 1991, ApJ, 383, 784

Mukherjee, A., & Bhattacharyya, S. 2011, ApJL, 730, L32

Mukherjee, A., & Bhattacharyya, S. 2012, ApJ, 756, 55

Nowak, M. A., & Vaughan, B. A. 1996, MNRAS, 280, 227

Nowak, M. A., Vaughan, B. A., Wilms, J., Dove, J. B., & Begelman, M. C. 1999, ApJ, 510, 874

Peille, P., Barret, D., & Uttley, P. 2015, ApJ, 811, 109

Psalis, D. 2008, LRR, 11, 9

Revnivtsev, M. G., & Gilfanov, M. R. 2006, A&A, 453, 253

Stevens, A. 2018, PhD thesis, Univ. Amsterdam

Sunyaev, R. A., & Shakura, N. I. 1986, SvAL, 12, 117

Troyer, J. S., & Cackett, E. M. 2017, ApJ, 834, 131

Uttley, P., Cackett, E. M., Fabian, A. C., Kara, E., & Wilkins, D. R. 2014, A&ARv, 22, 72

Uttley, P., Wilkinson, T., Cassatella, P., et al. 2011, MNRAS, 414, L60

van der Klis, M. 1998, AdSpR, 22, 925

van der Klis, M. 2000, ARA&A, 38, 717

van der Klis, M. 2006, AdSpR, 38, 2675

van der Klis, M., Hasinger, G., Stella, L., et al. 1987, ApJL, 319, L13

van der Klis, M., Wijnands, R. A. D., Horne, K., & Chen, W. 1997, ApJL, 481, L97

van Straaten, S., Ford, E. C., van der Klis, M., Méndez, M., & Kaaret, P. 2000, ApJ, 540, 1049

Vaughan, B. A., & Nowak, M. A. 1997, ApJL, 474, L43

Vaughan, B. A., van der Klis, M., Méndez, M., et al. 1998, ApJL, 509, L145

Wijnands, R., Homan, J., van der Klis, M., et al. 1997, ApJL, 490, L157

Wilkinson, T., & Uttley, P. 2009, MNRAS, 397, 666

Wilson-Hodge, C. A., Ray, P. S., Gendreau, K., et al. 2017, ResPh, 7, 3704

Zhang, S. N., Feroci, M., Santangelo, A., et al. 2016, Proc. SPIE, 9905, 99051Q

Zhang, W., Lapidus, I., White, N. E., & Titarchuk, L. 1996, ApJL, 469, L17