Unveiling hidden variability components in accreting X-ray binaries using both the Fourier power and cross spectra

Abstract

We present a novel method for measuring the lags of (weak) variability components in neutron-star and black-hole low-mass X-ray binaries (LMXBs). For this we assume that the power and cross spectra of these sources consists of a number of components that are coherent in different energy bands, but are incoherent with one another. The technique is based on fitting simultaneously the power spectrum (PS) and the Real and Imaginary parts of the cross spectrum (CS) with a combination of Lorentzian functions. We show that, because the PS of LMXBs is insensitive to signals with a large Imaginary part and a small Real part in the CS, this approach allows us to uncover new variability components that are only detected in the CS. We also demonstrate that, contrary to earlier claims, the frequency of the type-C quasi-periodic oscillation (QPO) in the black-hole binary GRS 1915+105 does not depend on energy. Rather, the apparent energy dependence of the QPO frequency can be explained by the presence of a separate QPO component with a slightly higher frequency than that of the QPO, whose rms amplitude increases faster with energy than the rms amplitude of the QPO. From all the above we conclude that, as in the case of the PS, the CS of black-hole and neutron-star binaries can be fitted by a combination of Lorentzian components. Our findings provide evidence that the frequency-dependent part of the transfer function of these systems can be described by a combination of responses, each of them acting over relatively well-defined time scales. This conclusion challenges models that assume that the main contribution to the lags comes from a global, broadband, transfer function of the accreting system.

1 Introduction

The X-ray light curves of accreting black-hole and neutron-star X-ray binaries show a complex pattern of variability with time scales ranging from milliseconds to years (e.g. Levine et al., 2006; Belloni et al., 2011; Ingram & Motta, 2019, and references therein). In the last four decades our understanding of the properties of these sources has advanced significantly due to the realisation that their power spectra can be decomposed into several variability components with well-defined properties (e.g., van der Klis, 1994; Nowak, 2000) that correlate with each other (Belloni et al., 2002), and other source properties (see below).

An effective approach to explore the geometry of the accretion flow in these systems is through the energy- and frequency-dependent phase lags of the variability (van der Klis et al., 1987; Nowak et al., 1999a). The phase lags measure the phase angle in the complex Fourier plane of the cross vector of correlated signals measured in two energy bands as a function of Fourier frequency (Vaughan & Nowak, 1997; Nowak et al., 1999a).

In the last decade, several models have been proposed to constrain the physical and geometrical properties of the accretion flow from measurements of the rms amplitude and lag spectra of these sources (e.g. Ingram & Done, 2011; Ingram et al., 2016; Mastroserio et al., 2018; Kylafis et al., 2020; Karpouzas et al., 2020; Bellavita et al., 2022). While the rms amplitude of the different variability components can be obtained from fits to the power spectra of these sources with a combination of Lorentzian functions (Nowak, 2000; Belloni et al., 2002), the lags are obtained from measurements that do not separate the contribution of those individual components (e.g., van der Klis et al., 1987; Reig et al., 2000, but see Nowak et al. 1999b).

Specifically, until now the method of measuring the lags of a broadband noise (BBN) component or a quasi-periodic oscillation (QPO) in the power spectrum of low-mass X-ray binaries (LMXBs) was the following:

1. Compute the power spectrum of the light curve of the source (e.g., van der Klis, 1989b) in the broadest energy band available.

2. Find the interesting components in the power spectrum (e.g., Psaltis et al., 1999). In the case of a BBN component, identify the minimum and maximum frequency over which one wants to measure the lags (e.g. Reig et al., 2000; Altamirano & Méndez, 2015; Kara et al., 2019). In the case of a QPO, find the centroid frequency, $\nu_{0}$, and the full width at half maximum (FWHM), $\Delta$, of the QPO to measure the lags in the range from $(\nu_{0}-x\Delta)$ to $(\nu_{0}+x\Delta)$, where usually $x=1/2$ (e.g., Wijnands et al., 1999)

3. Compute the cross spectrum using light curves in two energy bands (e.g., van der Klis et al., 1987; Vaughan et al., 1998; Nowak et al., 1999a; Uttley et al., 2014).

4. Compute the average of the Real and Imaginary parts of the cross spectrum in the selected frequency range (e.g., van der Klis et al., 1987; Lewin et al., 1988; Reig et al., 2000; Belloni et al., 2021).

5. Compute the ratio of the Imaginary to the Real parts and take the inverse tangent function of that to get the phase lag, $\Delta\phi$.

6. If one is interested in the time lag, $\Delta\tau$, in that frequency range, take some representative frequency in the range, e.g., $\langle\nu\rangle=(\nu_{\rm min}+\nu_{\rm max})/2$ in the case of the BBN, or $\nu_{0}$ in the case of the QPO, and divide $\Delta\phi$ by $2\pi$ times that frequency to get $\Delta\tau$ (e.g., Nowak et al., 1999a; de Avellar et al., 2013; Barret, 2013). Here we call this the traditional method.

When using this method, the (underlying) assumption is that the component of interest dominates in the power and the cross spectra over the frequency range of interest. This is usually true for kilohertz QPOs in neutron-star systems (kHz QPOs; Méndez & Belloni, 2021) and high-frequency QPOs in black-hole systems (Morgan et al., 1997; Belloni et al., 2012; Méndez et al., 2013), since these QPOs appear in a frequency range in which the only other component present in the power spectrum is the Poisson noise and the QPO dominates. In many occasions this is also the case for low-frequency QPOs (see Belloni & Stella, 2014; Ingram & Motta, 2019, and references therein) that appear in a region of the power spectrum in which other components are present (e.g., a BBN component, sub-harmonics and second or higher harmonics of the QPO, or other timing components not harmonically related to the QPO of interest), when the QPO itself is significantly stronger than those other components (e.g., Zhang et al., 2017; Zhang et al., 2020). However, this is not always the case (e.g., Nowak et al., 1999a; van Straaten et al., 2002), and if another component contributes significantly to the power and cross spectra in the frequency range of interest, the above procedure fails to give the lags of the component of interest. The situation is even worse when this other component dominates the variability in that frequency range (e.g., Ma et al., 2021; Alabarta et al., 2022).

In the past this problem has been tackled in different ways; for instance Ma et al. (2021, 2023b) find the “net” lag spectrum of a QPO by subtracting from the lags of the QPO the average lags over a frequency range just outside the QPO. This method is mathematically incorrect: The two measured lags give the angles of the cross vectors in Fourier space, whereas the angle of the net cross vector is not the difference of the angles of the two individual cross vectors, but depends also on their magnitudes. Alternatively, one could subtract the Real and Imaginary parts of the cross spectrum outside the QPO instead of the lags, which would be mathematically correct (as far as we are aware this has not been done in the literature). However, unless the Real and Imaginary parts of the unwanted components are constant with frequency, their contribution over the frequency range of interest is not necessarily the same as what one measures over the adjacent frequency range. Furthermore, this procedure would over-correct the lags whenever the wings of the component of interest contribute significantly to the power in the adjacent frequency range.

Another approach that has been used in the literature is to ignore the frequency range in which the unwanted components contribute significantly to the variability (e.g., Wang et al., 2021, 2022). This, however, does not work when parts of the unwanted components contribute significantly to the frequency range of interest (see, for instance, Fig. 4 in Wang et al., 2021, and the discussion in §4.2). In all these cases, because the power and cross spectra of these sources contain many overlapping components (e.g., Nowak, 2000; Zhang et al., 2020), it is often impossible to find a “clean” frequency interval in which one can measure the variability of the component of interest without contamination from other unwanted components.

If one has a model of the cross spectrum both as a function of energy and Fourier frequency, instead of extracting the lags one can fit the model directly to the Real and Imaginary parts of the cross spectrum as a function of energy in a number of frequency bands (e.g., Mastroserio et al., 2018) or vice versa (e.g., Rapisarda et al., 2014). For this, one uses the transfer function of the system (e.g., between the corona and the accretion disc, Reynolds et al., 1999, we expand further on the transfer function in §2.3) as a function of photon energy and Fourier frequency, which is a mathematically valid approach as long as the model of that transfer function includes all the components that contribute to the variability. An advantage of this is that one can build the possible interactions of the individual components in the model of the transfer function (e.g., Ingram & van der Klis, 2013; Rapisarda et al., 2014). Some models, however, only consider the global (broadband in Fourier frequency) transfer function of the system, and explicitly exclude the frequency range in which the variability is dominated by narrow components, like QPOs (e.g., Mastroserio et al., 2018; Ingram et al., 2019).

It is traditional to fit power spectra of LMXBs with a linear combination of Lorentzian functions (usually called a multi-Lorentzian model; see, e.g., Nowak, 2000; Belloni et al., 2002). The centroid frequency, $\mathrm{FWHM}$, rms fractional amplitude, and rms and lag spectra of each of these Lorentzians correlate with each other (e.g., van Straaten et al., 2002, 2003; Casella et al., 2004; Altamirano et al., 2005, 2008; Motta et al., 2011), with properties of other Lorentzians in the power spectrum (e.g., Méndez et al., 2001), with the total intensity and different hardness ratios of the source (e.g., Méndez et al., 1999; Jonker et al., 2002; van Straaten et al., 2005), with the spectral parameters of the different components used to fit the total energy spectrum of the source (e.g., Di Matteo & Psaltis, 1999; Homan et al., 2001; Remillard et al., 2002; Pottschmidt et al., 2003; Vignarca et al., 2003; Zdziarski et al., 2005; Shaposhnikov & Titarchuk, 2009; Motta et al., 2009; Motta et al., 2010, 2011; Reig et al., 2013; Stiele et al., 2013; Shidatsu et al., 2014; Grinberg et al., 2014; Altamirano & Méndez, 2015; Kalamkar et al., 2015; Reig et al., 2018; De Marco et al., 2015; De Marco et al., 2017; De Marco et al., 2021), with the fluxes of those components (e.g., Markwardt et al., 1999; Sobczak et al., 2000) and the total flux (e.g., Sobczak et al., 2000; Remillard et al., 2002) or luminosity (Ford et al., 2000) and with radio properties of the source (e.g., Muno et al., 2001; Fender et al., 2004; Méndez et al., 2022; García et al., 2022). These correlations (the list of references in this paragraph is certainly incomplete, but it is impossible to do justice to all the papers that discuss this topic) offer compelling evidence that those Lorentzian components are not just a useful empirical description of the power spectrum, but they represent physical phenomena in the system with (rather) well-defined characteristic time scales.

In the case of GX 339-4, Nowak et al. (1999b) observed that whenever one Lorentzian component dominates the power spectrum (PS) coherence function (see §2) approaches unity and the phase lags exhibit a relatively flat ‘shelf’; whenever two Lorentzian components intersect in the PS, there is a decline in the coherence function and a transition from one characteristic phase-lag shelf to another (c.f., for instance, Figs. 1 and 10 in Nowak et al., 1999a). This result indicates that, at least in this source, there is a relation between the individual components that fit the PS and the lags in the cross spectrum (CS).

Here we propose a new method to measure the lags of these components in the cross spectrum of LMXBs. In §2 we describe the mathematical foundation of this approach. In §3 we give examples of the application of the new proposed method to a number of cases, and demonstrate the advantages of this method in comparison with the traditional one. In fact, in doing this we unveil new variability components that have gone undetected before, and we show that ignoring these components led to wrong conclusions regarding some characteristics of the QPOs. In §4 we summarise our findings and discuss some unexpected consequences of our results that have significant impact upon models of the variability that have been presented in the literature. Finally, in §5 we sum up our conclusions.

2 Mathematical formalism

Suppose we have two noiseless time series, $x(t)$ and $y(t)$, that represent the X-ray intensity of a variable source measured in two different energy bands, with corresponding complex Fourier transforms $X(\nu)$ and $Y(\nu)$. We define the power spectrum of each of these series, $G_{xx}(\nu)=\langle X(\nu)X^{*}(\nu)\rangle=\langle|X(\nu)|^{2}\rangle$, $G_{yy}(\nu)=\langle Y(\nu)Y^{*}(\nu)\rangle=\langle|Y(\nu)|^{2}\rangle$, and the cross spectrum between the two series, $G_{xy}(\nu)=\langle X(\nu)Y^{*}(\nu)\rangle=\langle|X(\nu)||Y(\nu)|e^{i\Delta\phi_{xy}(\nu)}\rangle$, where $\Delta\phi_{xy}(\nu)$ is the phase lag between the two series at frequency $\nu$, and the angle brackets indicate averaging over an ensemble of measurements of $X(\nu)$ and $Y(\nu)$ (see Bendat & Piersol, 2010, for details). We can then define the coherence function between the two series (Bendat & Piersol, 2010; Vaughan & Nowak, 1997; Nowak et al., 1999a):

$$\gamma^{2}_{xy}(\nu)=\frac{|G_{xy}(\nu)|^{2}}{G_{xx}(\nu)G_{yy}(\nu)}.$$ (1)

From the above definition it is clear that if the two time series are related by a linear transformation, $G_{yy}(\nu)=|H(\nu)|^{2}G_{xx}(\nu)$ and $G_{xy}(\nu)=H(\nu)G_{xx}(\nu)$, the coherence function is unity at all frequencies. The complex function $H(\nu)=|H(\nu)|e^{-i\Delta\phi(\nu)}$ is the frequency response function, sometimes also called the transfer function, of the system (c.f., Bendat & Piersol, 2010, see §2.3 for more details of the possible interpretations of the transfer function).

Let us consider that the time series $x(t)$ and $y(t)$ can be decomposed in a finite number of individual components that are significant over a limited frequency range. It is traditional to fit the power spectrum of LMXBs with a linear combination of Lorentzian functions (Nowak, 2000; Belloni et al., 2002),

$$\;L(\nu;\nu_{0},\mathrm{\Delta})=\frac{\Delta}{\pi+2\tan^{-1}{\left(\frac{2\nu_{0}}{\Delta}\right)}}\frac{1}{\left(\nu-\nu_{0}\right)^{2}+\left(\frac{\Delta}{2}\right)^{2}},$$ (2)

where $\nu_{0}$ and $\Delta$ are, respectively, the centroid frequency and the FWHM of each Lorentzian. (With the above definition the integral of the Lorentzian function from zero to infinity is one.) Some of these Lorentzians have a central frequency equal to zero and are therefore called zero-centred Lorentzians, while others have a central frequency that is different from zero and, depending on their quality factor, $Q=\nu_{0}/\Delta$, they are called either peaked noise if $Q<2$ or QPOs if $Q>2$ (e.g. van der Klis, 1994). Since this is an arbitrary definition (Belloni et al., 2002), we call these components QPOs regardless of their $Q$ values.

We will now work out the expression of the Real and Imaginary parts of the cross spectrum assuming, as in Nowak et al. (1999b), that:
(i) there are a number of input processes that are coherent with the corresponding output processes;
(ii) the individual input processes (as well as the individual output processes) are incoherent with one another. We note that we only require that two input components are incoherent with each other if they overlap in Fourier frequency. This assumption is not necessary if the components are very far apart from each other in frequency, as is for instance the case of harmonics/sub-harmonics of a QPO.

Assumption (i) is justified by the fact that, as far as this was measured, strong QPOs, which dominate the variability over a certain frequency range, have a coherence that is consistent with unity. The cleanest case is that of the kHz QPOs, which appear in a part of the power spectrum in which no other component contributes to the variability. This is important because the coherence drops if two or more components with different amplitudes and phases of their cross vectors contribute to the variability over the same frequency range (Vaughan & Nowak, 1997). For instance, de Avellar et al. (2013) showed that the coherence of one of the strongest kHz QPO in the neutron-star system 4U 1608–52 is consistent with unity across the QPO profile (see their Fig. 1). On the other hand, Troyer et al. (2018) showed that for the kHz QPOs in 14 neutron-star LMXBs the total rms amplitude in the power spectrum and the cross amplitudes in the cross spectra are consistent with being the same. For this to be the case the coherence at the QPO frequency must be unity.

The coherence function of strong low-frequency QPOs in black-hole LMXBs is also consistent with being unity. For instance, in XTE J1550–564, while over a broad range of Fourier frequencies the coherence is less than one, at the frequency of the type-C QPO the coherence always climbs back to one (see, e.g., Fig. 1 of Cui et al., 2000, and Fig. 6 of Rapisarda et al. 2017a). The same is true for the type-C QPO in GRS 1915+105 (Muno et al., 2001).

This holds true not only for narrow QPOs, but also for broader variability components. For instance, the power spectrum of Cyg X–1 consists (basically) of two broad Lorentzian components (e.g. Grinberg et al., 2014), while the coherence is about unity over the whole frequency range in which either of these components dominates, and it drops slightly when the two components cross (Nowak et al., 1999a). In fact the power spectrum of Cyg X-1 consists of more than just two components (see, e.g., Fig. 3 in Nowak, 2000), and the coherence drops whenever two of those components cross (see, e.g., Fig. 5 in Nowak et al., 1999a, and Fig. 5 in Rapisarda et al. 2017b). Interestingly, in Cyg X–1 the frequencies at which the coherence drops coincide with the frequencies at which the power spectrum shows breaks, and the lag spectrum shows shelves of more or less constant lags (compare, for instance, the bottom-left panel of Fig. 1, the bottom-right panel of Fig. 5 and the bottom-right panel of Fig. 10 in Nowak et al., 1999a).

Assumption (i) implies that we can write, for each Lorentzian, $G_{xy,i}=H_{i}(\nu)G_{xx,i}(\nu)$, and therefore each Lorentzian will have its own lag spectrum, $\arg{[H_{i}(\nu)]}=\Delta\phi_{xy,i}(\nu)$. We can therefore write that:

$$G_{xx}(\nu)=\sum_{i=1}^{n}G_{xx,i}(\nu)\coloneq\sum_{i=1}^{n}A_{i}L(\nu;\nu_{0,i},\Delta_{i})$$ (3)

$$G_{yy}(\nu)=\sum_{i=1}^{n}G_{yy,i}(\nu)\coloneq\sum_{i=1}^{n}B_{i}L(\nu;\nu_{0,i},\Delta_{i}),$$

where $A_{i},B_{i}\in\mathbb{R}$ are the integrated power, from zero to infinity, of each Lorentzian component in each of the two energy bands.

Assumption (ii), which allows us to decompose the power spectra as in eq. 3, is that the $n$ Lorentzians are uncoupled and their contribution to the power spectrum are additive. While it has been used repeatedly in the literature to fit the power spectra of LMXBs (e.g. Nowak, 2000; Belloni et al., 2002; van Straaten et al., 2003; Altamirano et al., 2008, etc.), this assumption would appear to contradict the linear relation between the broadband absolute rms amplitude and the X-ray flux (e.g. Uttley et al., 2005) in LMXBs and Active Galactic Nuclei; this rms-flux relation suggests that the variability should be due to a multiplicative rather than an additive process, favouring models of propagating mass accretion rate fluctuations (Arévalo & Uttley, 2006, see Ingram & van der Klis 2013 and Rapisarda et al. 2016 for the mathematical description and an implementation of the model of propagating fluctuations to produce the power and cross spectra of LMXBs). However, Heil et al. (2011) showed that the low-frequency QPO in the black-hole candidate XTE J1550–564 does not follow a single linear rms-flux relation, but that the slope of this relation depends on the frequency of the QPO. Rapisarda et al. (2017a), on the other hand, demonstrated that propagation models fail to reproduce the data, and concluded that the discrepancies between data and model are a generic issue inherent to these models, rather than being specific to any particular implementation of that scenario. Finally, since one computes power spectra of stationary time series, there must be a linear, additive, model that describes the data (Scargle, 2020). In any case, we can test the predictions of our hypothesis that the Lorentzians are mutually incoherent, as we explain below.

If we accept that assumption (ii) above holds, it is enough to work out the case of a signal with a power spectrum consisting of only one Lorentzian, and combine the results to get the power and cross spectra of a combination of Lorentzians at the end:

$$G_{xx}(\nu)=A\;L(\nu;\nu_{0},\Delta)$$ (4)

$$G_{yy}(\nu)=B\;L(\nu;\nu_{0},\Delta).$$

From eq. 1, since for each Lorentzian $\gamma^{2}_{xy}=1$, we can write:

$$|G_{xy}(\nu)|^{2}=\mathrm{Re}[G_{xy}(\nu)]^{2}+\mathrm{Im}[G_{xy}(\nu)]^{2}=$$ (5)

$$=G_{xx}(\nu)G_{yy}(\nu)=ABL^{2}(\nu;\nu_{0},\Delta).$$

The frequency-dependent phase lag between the two signals represented by the Lorentzian functions can be written in the usual way:

$$\Delta\phi_{xy}(\nu)=\tan^{-1}{\left(\frac{\mathrm{Im}[G_{xy}(\nu)]}{\mathrm{Re}[G_{xy}(\nu)]}\right)}=g(\nu;p_{j}),$$ (6)

where $g(\nu;p_{j})$ is a function of frequency with $m$ parameters $p_{j}$. It is easy to see that we can combine eqs 5 and 6 to get

$$\mathrm{Re}[G_{xy}(\nu)]=C\;L(\nu;\nu_{0},\Delta)\cos{[g(\nu;p_{j})]}$$ (7)

$$\mathrm{Im}[G_{xy}(\nu)]=C\;L(\nu;\nu_{0},\Delta)\sin{[g(\nu;p_{j})]},$$

where $C=\sqrt{A\;B}$.

In the general case, if we assume that the power spectra of two signals $x(t)$ and $y(t)$ can be decomposed into a finite number, $n$, of Lorentzian functions that are uncorrelated with each other, because the individual Lorentzians are coherent between the two energy bands at all Fourier frequencies, the Real/Imaginary part of the cross spectrum between the two signals is also a linear combination of Lorentzian functions multiplied by a cosine/sine function of the, in principle frequency-dependent, phase lag between the two corresponding Lorentzian signals:

$$\mathrm{Re}[G_{xy}(\nu)]=\sum_{i=1}^{n}C_{i}\;L(\nu;\nu_{0,i},\Delta_{i})\cos{[\Delta\phi_{xy,i}(\nu)]}$$ (8)

$$\mathrm{Im}[G_{xy}(\nu)]=\sum_{i=1}^{n}C_{i}\;L(\nu;\nu_{0,i},\Delta_{i})\sin{[\Delta\phi_{xy,i}(\nu)]},$$

with $\Delta\phi_{xy,i}(\nu)=g_{i}(\nu;p_{j,i})$ and $C_{i}=\sqrt{A_{i}B_{i}}$.

In practice, one fits the PS in two energy bands using eq. 3 and the CS between the same two bands using eq. 8, for an arbitrary choice of $g(\nu;p_{j,i})$ (more on this below). The goal is to find, for each Lorentzian in the model, the parameters $A_{i}$, $B_{i}$, $\nu_{0,i}$, $\Delta_{i}$ and $p_{j,i}$, $j=1...m$. Instead of the PS in two bands one can use the full-band PS, $G_{zz}$, which can also be written as a a combination of the same Lorentzian functions, $G_{zz}(\nu)=\sum_{i=1}^{n}D_{i}L(\nu;\nu_{0,i},\Delta_{i})$ (because, except for the normalisations, the Lorentzians that fit the QPOs in the individual bands also fit the QPOs in the full band); in that case, for each Lorentzian, the parameters of the model are $D_{i}$, $\nu_{0,i}$, $\Delta_{i}$, $C_{i}$ and $p_{j,i}$.

From here we can write the total, frequency-dependent, phase lags, $\Delta\phi(\nu)$, as:

$$\Delta\phi(\nu)=\tan^{-1}{\left(\frac{\sum_{i=1}^{n}C_{i}\;L(\nu;\nu_{0,i},\Delta_{i})\sin{[g_{i}(\nu;p_{j,i})]}}{\sum_{i=1}^{n}C_{i}\;L(\nu;\nu_{0,i},\Delta_{i})\cos{[g_{i}(\nu;p_{j,i})]}}\right)}.$$ (9)

(Notice that the total phase lags given in eq. 9 are different from the phase lags of each Lorentzian, $\Delta\phi_{xy,i}(\nu)=g_{i}(\nu;p_{j,i})$ defined in eq. 6.) One prediction of our proposal is that the same model that fits the power and cross spectra should correctly reproduce the phase-lag vs. frequency spectra, and in particular yield the same “shelves” observed in the plots of the phase-lag vs. Fourier frequency (see, for instance, Fig. 10 of Nowak et al., 1999a, especially the bottom-right panel).

If we use two PS and the CS, we can also write the total intrinsic coherence function of the two signals as a function of frequency:

$$\gamma^{2}_{xy}(\nu)=\frac{|G_{xy}(\nu)|^{2}}{G_{xx}(\nu)G_{yy}(\nu)}=$$ (10)

$$\left[\left(\sum_{i=1}^{n}C_{i}\;L(\nu;\nu_{0,i},\Delta_{i})\cos{[g_{i}(\nu;p_{j,i})]}\right)^{2}+\right.$$

$$\left.\left(\sum_{i=1}^{n}C_{i}\;L(\nu;\nu_{0,i},\Delta_{i})\sin{[g_{i}(\nu;p_{j,i})]}\right)^{2}\right]\times$$

$$\frac{1}{\left(\sum_{i=1}^{n}A_{i}\;L(\nu;\nu_{0,i},\Delta_{i})\right)\left(\sum_{i=1}^{n}B_{i}\;L(\nu;\nu_{0,i},\Delta_{i})\right)}.$$

(This is the same as eq. 10 in Vaughan & Nowak 1997 but for $n$ components.) From this, the second prediction of our proposal is that, if indeed the Lorentzian functions are incoherent with one another, the same model that fits the power and cross spectra should correctly reproduce the coherence function vs. Fourier frequency. Specifically, if two or more Lorentzian functions contribute significantly at some Fourier frequency, and the Lorentzians have different phase lags and different ratios of the Fourier amplitude in each energy band, the observed intrinsic coherence function should drop at those frequencies. Conversely, if the observed intrinsic coherence function drops at some Fourier frequency, the best-fitting model of the power and cross spectra should contain two or more Lorentzians with different properties of the cross vectors that contribute significantly to the variability at that Fourier frequency (see Vaughan & Nowak, 1997, for details). Finally, if a strong Lorentzian dominates the variability over a frequency range, and if our assumption (i) is correct, the coherence should tend to unity in that frequency range.

2.1 Phase-lag model

Since any function $g(\nu;p_{j})$ will satisfy eqs. 5 and 7, as we mention above, to fit data we need to specify these functions for each Lorentzian $i$ in the model. The most basic case is when each of the functions depends only on one parameter, $p_{i}=k_{i}$, and from that family of functions the two simplest cases are: (i) The phase lags are constant, $\Delta\phi_{xy,i}(\nu)=g_{i}(\nu;p_{i})=2\pi k_{i}$. In this case, for all Fourier frequencies, the cross vector of each individual Lorentzian function points in the same direction, at an angle $2\pi k_{i}$ radians with respect to the Real axis. (ii) The time lags are constant, and hence the phase lags increase linearly with frequency, $\Delta\phi_{xy,i}(\nu)=g_{i}(\nu;p_{i})=2\pi k_{i}\nu$. This could be the case if the lags are produced by delays of photons propagating in the accretion flow of LMXBs. In this case the cross vector rotates, always in the same direction, at a constant rate of $2\pi k_{i}$ radians per Hz in the Fourier plane as a function of Fourier frequency.

From the above it is apparent that, unless $g(\nu)$ is constant, the cross vector will rotate in some more or less complex way, perhaps changing the direction of the rotation, as a function of frequency in the Fourier plane. If this happens, the phase lags may wrap as the angle changes abruptly from $\pi$ at some Fourier frequency to $-\pi$ at the next frequency, or vice versa. For the phase lags not to wrap it must be true that $-\pi\leq g(\nu)<\pi$ over the frequency range in which the signal contributes significantly to the variability. (We would not be able to observe a phase wrap if it happened at frequencies where the power spectrum was dominated by noise.) For instance, if a component in the power spectrum has a time lag of $\sim 100$ ms that is constant with Fourier frequency, the cross vector associated with that component will rotate in the Fourier plane at $\sim 0.6$ radians per Hz, such that if the component extends over a frequency range of $\gtrsim 10$ Hz its phase lags will wrap at least once. Since, until now, such a wrap has not been reported, and since in LMXBs the time lags go as $\nu^{\sim-0.7}$ over a relatively broad frequency range (at least $\sim 100$ Hz; e.g., Nowak et al., 1999a), such that the phase lags increase very slowly, $\Delta\phi\sim\nu^{\sim 0.3}$, from the two simplest cases mentioned above the one of constant phase lags appears more natural. In what follows we will assume that the phase lags of individual components in the power spectrum are independent of Fourier frequency (constant phase lags) or increase linearly with frequency (constant time lags).

2.2 Average phase and time lags over a frequency range

In the above discussion the functions $g_{i}(\nu;p_{j,i})$ are the frequency-dependent phase lags of the individual Lorentzian components. A quantity sometimes given in the literature (e.g. Reig et al., 2000) is the total phase lag over the frequency range $\Delta\nu=\nu_{2}-\nu_{1}$:

$$\Delta\phi_{\Delta\nu}=\tan^{-1}{\left(\frac{\int_{\nu_{1}}^{\nu_{2}}\mathrm{Im}[G_{xy}(\nu)]d\nu}{\int_{\nu_{1}}^{\nu_{2}}\mathrm{Re}[G_{xy}(\nu)]d\nu}\right)}=$$ (11)

$$= ~{}\tan^{-1}{\left(\frac{\int_{\nu_{1}}^{\nu_{2}}\sum_{i=1}^{n}C_{i}\;L(\nu;\nu_{0,i},\Delta_{i})\sin{[g_{i}(\nu;p_{j,i})]}d\nu}{\int_{\nu_{1}}^{\nu_{2}}\sum_{i=1}^{n}C_{i}\;L(\nu;\nu_{0,i},\Delta_{i})\cos{[g_{i}(\nu;p_{j,i})]}d\nu}\right)}.$$

A related, but different, quantity is the average phase lag over the same frequency range:

$$\langle\Delta\phi\rangle_{\Delta\nu}=\frac{1}{\Delta\nu}\int_{\nu_{1}}^{\nu_{2}}\Delta\phi(\nu)d\nu=$$ (12)

$$\frac{1}{\Delta\nu}\int_{\nu_{1}}^{\nu_{2}}\tan^{-1}{\left(\frac{\mathrm{Im}[G_{xy}(\nu)]}{\mathrm{Re}[G_{xy}(\nu)]}\right)}d\nu=$$

$$=\frac{1}{\Delta\nu}\int_{\nu_{1}}^{\nu_{2}}\tan^{-1}{\left(\frac{\sum_{i=1}^{n}C_{i}\;L(\nu;\nu_{0,i},\Delta_{i})\sin{[g_{i}(\nu;p_{j,i})]}}{\sum_{i=1}^{n}C_{i}\;L(\nu;\nu_{0,i},\Delta_{i})\cos{[g_{i}(\nu;p_{j,i})]}}\right)d\nu}.$$

It is easy to show that when the cross spectrum consists of a single Lorentzian with $g(\nu;k)=2\pi k$ (constant phase lags), $\langle\Delta\phi\rangle_{\Delta\nu}=2\pi k=\Delta\phi_{\Delta\nu}$, but for any other form of $g(\nu,p_{j})$, or if the cross spectrum consists of a linear combination of Lorentzians, $\langle\Delta\phi\rangle_{\Delta\nu}\neq\Delta\phi_{\Delta\nu}$, and hence $\langle\Delta\phi\rangle_{\Delta\nu}$ should not be reported.

It is not possible to define the total time lag over the frequency range $\Delta\nu$, $\Delta\tau_{\Delta\nu}$, in a similar way as we defined $\Delta\phi_{\Delta\nu}$ in eq. 11, but a quantity commonly given in the literature is $\Delta\tilde{\tau}_{\Delta\nu}=\Delta\phi_{\Delta\nu}/(2\pi\langle\nu\rangle)$, where $\langle\nu\rangle=(\nu_{1}+\nu_{2})/2$. In the case of a narrow component, this is more or less similar to the usual definition for a pure sinusoidal function with frequency $\nu_{0}$, $\Delta\tau(\nu_{0})=\Delta\phi(\nu_{0})/(2\pi\nu_{0})$, but for a broad component, e.g., the BBN (e.g. Kara et al., 2019; Wang et al., 2022), these two expressions are not mathematically equivalent and hence $\Delta\tilde{\tau}_{\Delta\nu}$ is not a well-defined quantity and should not be used.

As with the phase lags, we can define the average time lags over a frequency range $\Delta\nu$:

$$\langle\Delta\tau\rangle_{\Delta\nu}=\frac{1}{\Delta\nu}\int_{\nu_{1}}^{\nu_{2}}\Delta\tau(\nu)d\nu=\frac{1}{\Delta\nu}\int_{\nu_{1}}^{\nu_{2}}\frac{\Delta\phi(\nu)}{2\pi\nu}d\nu,$$ (13)

which is equal to $\langle\Delta\phi\rangle_{\Delta\nu}/(2\pi\langle\nu\rangle)$ only when the time lags are constant with frequency, $\Delta\tau(\nu)=k$. For instance, if $\Delta\phi(\nu)=2\pi k$ (constant phase lags), $\langle\Delta\tau\rangle_{\Delta\nu}=\left(k/\Delta\nu\right)\ln{\left(\nu_{2}/\nu_{1}\right)}$, and the average time lags will reflect the (arbitrary) boundaries of the selected frequency range. Since in LMXBs it is generally true that the time lags extending over a broad frequency range are not constant (see, e.g., Fig. 1 in De Marco et al. 2021, or Fig. 3h in Wang et al., 2022), contrary to what is sometimes mentioned in the literature (e.g. Belloni & Bhattacharya, 2022) one should not report $\langle\Delta\tau\rangle_{\Delta\nu}$ to characterise those lags either. (Notice that rebinning the time-lag spectrum as, e.g., in Wang et al., 2022, is equivalent to this, and should therefore not be done.)

Table 1: Definitions and acronyms

Acronym/text shortcut	Explanation
PS	Power density spectrum
CS	Either the cross spectrum, or the Real and Imaginary parts of the CS, depending on the context
CV	Cross vector in the Fourier plane with phase $\Delta\phi(\nu)$ and modulus $|\mathrm{CV}(\nu)|$
Phase-lag frequency spectrum	Phase lags between light curves in two energy bands vs. Fourier frequency
Derived model	Model of the phase-lag frequency spectrum derived from the model fitted to the Real and Imaginary parts of the CS
	Model of the intrinsic coherence function derived from the model fitted to PS in two energy bands and the Real and
	Imaginary parts of the CS in those same two bands.
	These models are not fitted to the lags or the intrinsic coherence function
Phase lags of a QPO	Phase lags of the Lorentzian that fits the QPO, computed as $\tan^{-1}{\left(N_{\rm IM}/N_{\rm RE}\right)}$, where $N_{\rm RE}$ and $N_{\rm IM}$ are the integrals
	from zero to infinity of the Lorentzian functions used to fit, respectively, the Real and Imaginary parts of the CS
Traditional phase lags of a QPO	$\tan^{-1}{\left(\langle\mathrm{IM}\rangle/\langle\mathrm{RE}\rangle\right)}$, where $\langle\mathrm{RE}\rangle$ and $\langle\mathrm{IM}\rangle$ are the averages of, respectively, the Real and Imaginary parts of the
	CS over a fixed frequency range across the QPO profile
Constant phase-lags model	Model of the CS assuming constant phase lags with Fourier frequency for each individual Lorentzian component,
	$\Delta\phi_{i}(\nu)=g_{i}(\nu;p_{j,i})=2\pi k_{i}$ (see §2.1)
Constant time-lags model	Model of the CS assuming constant time lags with Fourier frequency for each individual Lorentzian component,
	$\Delta\tau_{i}(\nu)=k_{i}$ such that $\Delta\phi_{i}(\nu)=g_{i}(\nu;p_{j,i})=2\pi k_{i}\nu$ (see §2.1)

2.3 Transfer/response function

Any linear time-invariant system is characterised by an impulse response function, $h(\tau)$, which gives the output of the system at any time, $y(t)$, to an input, $x(t)$, applied a time $\tau$ before, $y(t)=\int_{-\infty}^{\infty}h(\tau)x(t-\tau)d\tau$. The Fourier transform of the impulse response function, $H(\nu)$, is the frequency response or transfer function of the system (Bendat & Piersol, 2010). If $H(\nu)$ is the same for all realisations of the two processes, the processes are said to be coherent at frequency $\nu$ (see §2).

In a manner similar to Bendat & Piersol (2010), one can consider $x(t)$ and $y(t)$ as correlated light curves of a source in two energy bands, as discussed earlier. Here, $H(\nu)$ represents the transfer function that yields the light curve in band $y$ based on the light curve in band $x$. Alternative, one can think of $x(t)$ and $y(t)$ as the signals of two physical components of the system, e.g., the accretion disc and the corona, such that the transfer function gives the signal of one in terms of the signal of the other (e.g., Reynolds et al., 1999; Mastroserio et al., 2018; Ingram et al., 2019). In the first case, for light curves in two energy bands, the transfer function depends only upon Fourier frequency (for the given energy bands), and it refers to two observed quantities, each of them (possibly) being produced by a combination of more than one physical mechanism (e.g., the accretion disc, the corona, etc.). In the second case, the transfer function refers to the physical mechanisms themselves, and it depends both on Fourier frequency and energy; if one wants to compute light curves or power and cross spectra in certain energy bands to compare with observations, one needs to integrate the transfer function over energy (see Mastroserio et al., 2018, 2019, for details).

In this paper we will mostly use the term transfer function to refer to the first case. Specifically, as stated in §2, $\arg{[H(\nu)]}$ gives the phase-lags of light curve $y(t)$ with respect to light curve $x(t)$ at frequency $\nu$. Hence, when we fit Lorentzian functions to the power and cross spectra of a source, as explained above, we will be talking of the frequency-dependent lags of the Lorentzian that represents the variability in the power and cross spectra in one energy band with respect to the Lorentzian that represents the same variability in another energy band. In some passages, when we talk about global transfer or response functions, we refer to the second case.

3 Application to data

In this section we apply the formalism presented in §2 to five observations of three black-hole X-ray binaries: an observation of GX 339–4 and two observations of GRS 1915+105 with the Proportional Counter Array (PCA; Jahoda et al., 2006) on board the Rossi X-ray Timing Explorer (RXTE; Bradt et al., 1993), and two observations of MAXI J1820+070 with the Neutron Star Interior Explorer (NICER; Gendreau et al., 2012). Each of the examples highlights certain aspects of the method or provide useful insight into aspects of the variability of accreting X-ray binaries. (We have already presented results obtained with an initial version of this method in Alabarta et al. 2022 and Peirano & Méndez 2022.)

In each of the subsections we first give very briefly the necessary information of how we process the data and produce power, cross and lag spectra and coherence function of the source, and we then give the results of fitting those data. In all cases, when we fit simultaneously the power spectrum and the Real and Imaginary parts of the cross spectrum , we link the frequency and FWHM of each Lorentzian in the model so that these parameters are the same in the PS and CS. In some cases we fix the frequency and FWHM of all the Lorentzians to the values we obtain from an initial fit to the PS alone, while in other cases we leave these parameters free to vary, but we always link them across the PS and the Real and Imaginary parts of the CS. In each case we indicate whether we do one or the other. On the other hand, we always leave free the normalisations of all the Lorentzians and allow them to vary independently in the PS and the CS. If we plot the phase-lag frequency spectrum or the intrinsic coherence function together with the best fitting model, we do not fit the model to the phase-lag frequency spectrum or the intrinsic coherence function, but we derive the model from the best-fitting model to the Real and Imaginary parts of the CS. We call this the “derived model” of the phase-lag spectrum/intrinsic coherence function.

To simplify the reading we give the definition of some acronyms and language shortcuts that we use often in the text in Table 1. For completeness, we also give the acronym in the text the first time we define it.

We use the recommended tools for each mission to extract calibrated clean event files for each observation. When extracting events in certain energy bands, we always use the channel space of the detector, and convert those values into energies taking into account the channel to energy conversion for each instrument. In all cases we give both the channel and the (approximate) energy range that we use to analyse the data.

Once we have obtained the final event files, we use GHATS to compute the Fast Fourier Transform (FFT) in each band of interest. Except in one case (see below), we always compute the FFT of the data from the full observation (we give the exact ObsIDs that we use in each subsection). For this we compute FFTs on segments of a given duration, $T_{\rm FFT}<T_{\rm exp}$, where $T_{\rm exp}$ is the total duration of the observation, with a time resolution, $\Delta t$, that allows us to reach (and exceed) the frequency of interest in the PS and CS, and we average the PS and CS of the individual segments (van der Klis, 1989b). The number of segments that we average for each PS and CS is therefore $M=\mathrm{integer}(T_{\rm exp}/T_{\rm FFT})$. From this, the minimum frequency and the frequency resolution of the PS and CS are $\nu_{\rm min}=\Delta\nu=1/T_{\rm FFT}$, while the Nyquist frequency is $\nu_{\rm Nyq}=1/(2\Delta t)$. The averaging of the PS and CS per observation increases the signal-to-noise ratio by a factor $\sim\sqrt{M}$. We further rebin the average PS and CS of an observation logarithmically in frequency such that the size of a bin increases by a factor $\approx 1.023(=10^{1/100})$ compared to the size of the previous bin to increase the signal-to-noise ratio further.

We use Xspec v.12.13.0 (Arnaud, 1996) to fit the PS and CS. We always start by first fitting the full-band PS with a model with only one Lorentzian, and we add a new Lorentzian until the reduced $\chi^{2}$ is about one and the fit shows no significant structured residuals. The final model of the PS, consisting of $n$ Lorentzians (eq. 3), is $\sum_{i=1}^{n}D_{i}L(\nu;\nu_{0,i},\Delta_{i})$. In §3.1-3.4 we fit simultaneously the full-band PS, $G_{zz}(\nu)$, and the the Real and Imaginary parts of the CS in two bands, $\mathrm{Re}\left[G_{xy}(\nu)\right]$ and $\mathrm{Im}\left[G_{xy}(\nu)\right]$, while in §3.1, §3.4 and §3.5 we also use the PS in the two bands, $G_{xx}(\nu)$ and $G_{yy}(\nu)$, to be able to compute the coherence function.

When we fit the PS and the CS simultaneously, we need to assume the frequency dependence of the phase lags (see §2.1). Here we only consider two cases: either the phase lags or the time lags are constant with frequency (constant phase-lags or constant time-lags model, respectively); although it is possible that the phase lags of different variability components depend differently upon frequency, here we always assume that the same functional dependence applies to all the Lorentzian components. We mention in the text and in the Figure and Table captions the functional dependence of the lags that we use when we fit the data.

Specifically, when we use the constant phase-lags model we fit functions of the form $\sum_{i=1}^{n}C_{i}L(\nu;\nu_{0,i},\Delta_{i})\cos{(2\pi k_{i})}$ to the Real part of the CS and $\sum_{i=1}^{n}C_{i}L(\nu;\nu_{0,i},\Delta_{i})\sin{(2\pi k_{i})}$ to the Imaginary part, where $C_{i}$ is the integral from zero to infinity of the modulus squared of the cross vector (CV) and $2\pi k_{i}$ are the phase lags for each component (eq. 8 with $\Delta\phi_{xy,i}(\nu)=2\pi k_{i}$). The $C_{i}$ and $k_{i}$ are free parameters that need to be fitted.

When we use the constant time-lags model we fit the Real part of the CS with $\sum_{i=1}^{n}C_{i}L(\nu;\nu_{0,i},\Delta_{i})\cos{(2\pi k_{i}\nu)}$ and the Imaginary part with $\sum_{i=1}^{n}C_{i}L(\nu;\nu_{0,i},\Delta_{i})\sin{(2\pi k_{i}\nu)},$ where the $k_{i}$ are now the time lags of each component (eq. 8 with $\Delta\phi_{xy,i}(\nu)=2\pi k_{i}\nu$).

Since we normalise the PS and CS to units of fractional rms squared per Hz (Belloni & Hasinger, 1990), the integrated fractional rms amplitude of a Lorentzian is the square root of its normalisation. We consider that a Lorentzian is significantly needed in the model of the PS if the normalisation of that Lorentzian divided by its 1-$\sigma$ error is larger than 3. For the CS we take the modulus squared of the CV and its error. If the errors are asymmetric, we take the negative error to calculate the significance. Unless otherwise indicated, we give the 1-$\sigma$ error for one parameter and, if a Lorentzian is absent, we give the 95% confidence upper limit to the rms fractional amplitude of that Lorentzian in the power or the cross spectrum.

In the literature, authors call “the phase lags of a QPO” to the inverse tangent function of the ratio of the average of the Imaginary part to the average of the Real part of CS; these averages are computed over a fixed frequency range across the QPO profile, usually covering a FWHM, from $\nu_{\rm QPO}-\Delta/2$ to $\nu_{\rm QPO}+\Delta/2$. From this, “the time lags of the QPO” are defined as the phase lags of the QPO divided by $2\pi$ times the centroid frequency of the QPO (but see §2.2 for the problems of this definition). We call these quantities the “traditional phase/time lags of a QPO”.

On the contrary, following the explanation in §2.1, here we call “the phase lags of a QPO” to the quantity $2\pi k_{i}$ for the Lorentzian that fits the QPO in the CS for the constant phase-lags model, and “the time lags of a QPO” to the parameter $k_{i}$ for the Lorentzian that fits the QPO in te CS for the constant time-lags model. As we explained in §2.2, this means that we cannot compute the time lags of a QPO when we use the constant phase-lags model, and vice versa. Withal, the phase and time lags of a QPO are not related via the centroid frequency of the QPO as in the case of the traditional definition of the phase and time lags of the QPO (see §2).

Finally, as is customary, we say that lags are hard or positive if the high-energy photons lag the low-energy ones, and soft or negative if the opposite happens.

Figure 1: Top left panel: Full-band PS of the RXTE observation 92035-01-03-06 of GX 339–4 fitted with a model (solid thick line) consisting of seven Lorentzian functions (thin dotted lines). Middle left panel: Phase-lag vs. Fourier frequency (middle panel) together with the derived model (thick solid line) from the CS (not shown). During the fit, we let the centroid frequency and FWHM of each component free to vary but we link them to be the same in the PS and the CS, further assuming the constant phase-lags model (see Table 1 and §2.1). Bottom left panel: Residuals of the phase lags with respect to the derived model. The vertical line indicates the centroid frequency of the narrow QPO. Right column: Real (top) and Imaginary (middle) parts of the CS with the best-fitting model assuming constant time lags. The bottom panel shows the intrinsic coherence together with the derived model. For each Lorentzian in the model, the vertical lines mark $\nu_{\rm max}$, the frequency at which the Lorentzian function peaks in the $\nu\times P_{\nu}$ representation (Belloni et al., 2002).

3.1 Case study 1: The complex lag spectrum of the type-C QPO in GX 339–4

We use the RXTE observation 92035-01-03-06 of GX 339–4 that was analysed by Zhang et al. (2017, see their Fig. 1) and Altamirano & Méndez (2015, see their Fig. 3). This observation contains a prominent type-C QPO (see Wijnands et al., 1999; Remillard et al., 2002; Casella et al., 2005, for the definition of the three types of low-frequency QPOs, types A, B and C, in black-hole X-ray binaries) at $\sim 5.5$ Hz that shows a rather complex phase-lag frequency spectrum: The hard phase lags around the QPO are double-peaked, with the minimum centred on the QPO centroid frequency (see top and middle panels of Fig. 1)

Although no systematic study of this type of lag-frequency spectra exists, this is not an isolated case, and a similar trend is apparent in several other observations of GX 339–4 (see Fig. 1 of Zhang et al., 2017, for other examples) and GRS 1915+105 (Zhang et al., 2020, see also §3.2 below).

We take the Binned and Event Mode files of this observation to compute the full-band (channels 0 to 249) PS and CS of photons in the $5.4-115$ keV band (channels 14 to 249) with respect to those in the $1.95-5.4$ keV band (channels 0 to 13). Both for the PS and the CS we use $T_{\rm FFT}=16$ s, yielding $\nu_{\rm min}=\Delta\nu=1/16$ Hz, at a time resolution $\Delta t=1/2048$ s such that $\nu_{\rm Nyquist}=1024$ Hz.

We initially fit the PS in the range $0.5-20$ Hz adding one Lorentzian at a time as described in §3. A fit with six Lorentzians gives $\chi^{2}=111.8$ for 87 degrees of freedom (dof), while that with seven Lorentzians gives $\chi^{2}=93.4$ for 84 dof, and all Lorentzians are at least 4-$\sigma$ significant. Adding an eighth Lorentzians does not improve the fit significantly, so we stop at seven.

Constant phase-lags model: We next fit simultaneously the full-band PS and the CS with the same number of Lorentzians, fixing the frequency and FWHM of each Lorentzian to the values that we obtain from the fits to the PS; for this fit we assume the constant phase-lags model. The fit gives $\chi^{2}=305$ for 294 dof, with structured residuals around 5 Hz (not shown). We next let the frequency and FWHM free but linked in the PS and CS for each Lorentzian. This fit, shown in Figure 1, gives $\chi^{2}=287.9$ for 280 dof; while it is not statistically better than the one with fixed frequency and FWHM of the Lorentzians, this fit shows no structured residuals around 5 Hz. Because of this, and because later on we find that in other cases the fit with fixed frequencies and FWHM does not work, we adopt this one as the final model.

As it is apparent from Figure 1, the QPO is actually fitted with two Lorentzians, a relatively broad one that fits the wings, and a relatively narrow one that fits the central part of the QPO profile. The broad QPO component is reminiscent of the shoulder mentioned by Belloni et al. (1997, later on called the “hump” in ; see also ), and we therefore call this the QPO shoulder. Both the QPO and the QPO shoulder are significantly required by the fit, with a significance of 12-$\sigma$ and 20-$\sigma$, respectively. Since the QPO and the QPO shoulder have slightly, but significantly ($\sim 5.5$-$\sigma$), different frequencies, $\nu_{\rm QPO}=5.50\pm 0.01$ Hz and $\nu_{\rm shoulder}=5.78\pm 0.05$ Hz, the QPO profile appears to be slightly asymmetric (top panel of Fig. 1).

Table 2: Phase and time lags ($5.4-115$ keV vs. $1.95-5.4$ keV) of the QPO and the QPO shoulder in GX 339–4 with our method

Component	$\nu_{0}$ (Hz)	Phase lags${}^{\dagger}$ (rad)	Time lags${}^{\ast}$ (ms)
QPO	$5.50\pm 0.01$	$-0.004\pm 0.024$	$0.9\pm 0.7$
Shoulder	$5.78\pm 0.05$	$0.51\pm 0.03$	$10.8\pm 0.7$
${}^{\dagger}$Using the constant phase-lags model.
${}^{\ast}$Using the constant time-lags model.

In principle, it is possible that the asymmetric QPO profile is due to slight changes of the QPO frequency during the course of the observation, with the QPO spending most of the time at around $\sim 5.5$ Hz, and a small fraction of the time at $\sim 5.8$ Hz. To check this, we compare the full-band rms amplitude and phase lags between the two bands given above, $\Delta\phi$, of the two components, the QPO and the shoulder; we find rms${}_{\rm QPO}=4.9\pm 0.2\%$, rms${}_{\rm shoulder}=5.5\pm 0.3\%$, $\Delta\phi_{\rm QPO}=-0.004\pm 0.024$ rad and $\Delta\phi_{\rm shoulder}=0.51\pm 0.03$ rad. While the rms amplitudes of the QPO and the QPO shoulder are consistent with being the same, the phase lags are significantly ($\sim 13.4$-$\sigma$) different. Since the lags of the QPO do not change so rapidly with QPO frequency (Zhang et al., 2017; Zhang et al., 2020), we conclude that the QPO and the QPO shoulder are two different components.

Constant time-lags model: We also compute the time lags of the QPO and the QPO shoulder by fitting the PS and the CS assuming the constant time-lags model, and we find that time lags of the QPO are significantly different from those of the QPO shoulder. The fit yields $\chi^{2}=284.0$ for $290$ dof. We give the frequencies, phase and time lags of the QPO and QPO shoulder in Table 2.

As we explained in §2.1, in our method the phase (time) lags of the Lorentzians are parameters of the constant phase-lags (time-lags) model. Under the assumption that the phase (time) lags are constant, the time (phase) lags change with Fourier frequency across the QPO profile and, therefore, the time lags of a QPO are not exactly equal to the phase lags divided by $2\pi$ times the QPO centroid frequency.

Traditional lags: The case of the shoulder of the QPO highlights the problem of measuring lags in the traditional way, which we explained in §1. To show this we compute the traditional phase-lag, $\Delta\phi_{\rm trad}$, of the QPO in the $5.285-5.715$ Hz frequency range (one FWHM), $\Delta\phi_{\rm QPO,trad}=0.23\pm 0.01$ rad, which is significantly larger than the lag using our method (see Tab. 2), $\Delta\phi_{\rm QPO}=-0.004\pm 0.024$ rad. Similarly, the traditional phase lag of the QPO shoulder within one FWHM ($4.795-6.765$ Hz) is $\Delta\phi_{\rm shoulder,trad}=0.17\pm 0.02$ rad, significantly smaller than the phase lag using our method, $\Delta\phi_{\rm shoulder}=0.51\pm 0.03$ rad.

Likewise, the time lags computed using the traditional method over the same frequency ranges, $\Delta\tau_{\rm QPO,trad}=4.9\pm 0.5$ ms and $\Delta\tau_{\rm shoulder,trad}=6.5\pm 0.3$ ms, are significantly different from those obtained with our method, $\Delta\tau_{\rm QPO}=0.9\pm 0.7$ ms and $\Delta\tau_{\rm shoulder}=10.8\pm 0.7$ ms.

Figure 2: Left panel: Imaginary part of the CS of the RXTE observation 30703-01-34-00 of GRS 1915+105 fitted with six Lorentzian functions. The upper panel shows the data and the best-fitting model (thick solid line) and the Lorentzian functions (thin dotted lines), while the lower panel shows the residuals with respect to the best fitting model. Structured residuals are visible at $\sim 2$ Hz and $\sim 4$ Hz, just above the frequency of, respectively, the QPO fundamental and the second harmonic (see Table 3). Middle panel: PS of the same observation, but now fitted with eight Lorentzian functions. The upper panel shows the data and the best-fitting model (thick solid line) and the Lorentzian functions (thin dotted lines), while the lower panel shows the residuals with respect to the best fitting model. (The eighth Lorentzian is the high-frequency bump, Trudolyubov, 2001; Zhang et al., 2022b, which is only visible above 10 Hz.) Right panel: Imaginary part of the CS fitted with eight Lorentzians (upper panel), phase lags vs. Fourier frequency (middle panel) and residuals of the phase lags with respect to the derived model (bottom panel). In all cases, during the fits we let the centroid frequency and FWHM of each component free to vary but we link them to be the same in the PS and the CS, further assuming the constant phase-lags model (see Table 1 and §2.1).

Finally, we note that the second harmonic of the QPO is also fitted by two Lorentzian components (see Fig. 1), a narrow one at $11.07\pm 0.06$ Hz and a broad one at $11.4\pm 0.3$ Hz, both of them significant (11.2-$\sigma$ and 3.4-$\sigma$, respectively). The frequencies of both these components are consistent with being in a 2:1 relation with, respectively, the QPO fundamental and the QPO shoulder. The phase lags of these two components are, respectively, $0.25\pm 0.11$ rad and $-0.07\pm 0.04$ rad, and the time lags are, respectively, $3\pm 2$ ms and $-2.3\pm 0.4$ ms, marginally consistent within errors; we therefore cannot assess whether the QPO second harmonic and its shoulder are different components or the same component that drifts slightly during the observation period. Guided by this, in the next section we explore the presence of a shoulder of the QPO second harmonic in an observation of GRS 1915+105.

We subsequently compute PS of this same observation for the channels $0-13$ and $14-249$, using the same parameters mentioned at the start of this subsection. We fit these two PS together with the CS in the same bands above with a model consisting of nine Lorentzians, and assuming the constant time-lags model. In the right panels of Figure 1 we plot the Real (top) and Imaginary (middle) parts of the CS and the intrinsic coherence with the derived model (bottom). The derived model reproduces the data rather well, with the largest deviations appearing above $\sim 10$ Hz, where the signal-to-noise ratio drops very quickly. The vertical lines in the bottom-right panel mark the characteristic frequency of each Lorentzian, $\nu_{\rm max}=\sqrt{\nu^{2}_{0}+\left(\frac{\mathrm{FWHM}}{2}\right)^{2}}$ (Belloni et al., 2002). We note that the changes in the behaviour of the coherence happen at $\nu_{\rm max}$ and not at $\nu_{0}$ of the corresponding Lorentzians.

The coherence is $\sim 1$ at low frequencies where a single, broad, Lorentzian with $\nu_{\rm max}=1$ Hz dominates (see top-right panel of Figure 1 for the Lorentzians mentioned in this part of the text); the coherence then drops at $\sim 1.5$ Hz when a new Lorentzian with $\nu_{\rm max}=1.41$ Hz appears. This Lorentzian is not strong enough to lead to a peak in the coherence, but its presence interferes with the previous Lorentzian and leads to a drop. After that, the coherence drops more or less steadily over a frequency range in which the Lorentzian with $\nu_{\rm max}=1$ Hz and two other Lorentzians with, respectively, $\nu_{\rm max}=2.74$ Hz and $\nu_{\rm max}=4.04$ Hz, alternate dominance. The coherence goes up again at $\sim 5$ Hz, where the strong QPO starts to dominate the variability. As the frequency increases further, in the model the coherence decreases and increases again going from the QPO to the shoulder, but the errors in the data are too large to see that. At higher Fourier frequency the model is consistent with the data, but there, the coherence has large errors.

The results of this part show that the best-fitting model of the PS and the CS correctly describes both the phase-lag spectrum and the coherence function. These two outcomes align with the model’s predictions, which are based on the assumptions outlined in Section §2.

3.2 Case study 2: Shoulders of the fundamental and second harmonic of the type-C QPO in GRS 1915+105

For the next example we use the RXTE observation 30703-01-34-00 of GRS 1915+105 with a strong type-C QPO at $\sim 1.8$ Hz, which was analysed by Zhang et al. (2020). We select this observation because the type-C QPO shows a significant second harmonic at $\sim 3.6$ Hz that displays a double peak in the phase-lag frequency spectrum. Motivated by this, and after our findings in §3.1, we decided to study the potential presence of a shoulder to the QPO fundamental and the second harmonic.

We take the Single-Bit and Event Mode data of this observation to compute the full-band (channels 0 to 249) PS and the CS of photons with energies in the $5.7-100$ keV band (channels 14 to 249) relative to those with energies in the $2-5.7$ keV band (channels 0 to 13). Both for the PS and the CS we use $T_{\rm FFT}=256$ s, which yields $\nu_{\rm min}=\Delta\nu=1/256$ Hz, at a time resolution $\Delta t=1/8192$ s such that $\nu_{\rm Nyquist}=4096$ Hz.

We initially fit the full-band PS in the range $0.1-64$ Hz adding one Lorentzian at a time as described in §3, and find that we need a model consisting of six Lorentzians. The fit gives $\chi^{2}=144.7$ for 143 dof, and all Lorentzians are at least 3-$\sigma$ significant.

Constant phase-lags model: We next fit simultaneously the PS and the CS with the same number of Lorentzians, fixing the frequency and FWHM of each Lorentzian to the values that we obtain from the fit to the PS; for this fit we assume the constant phase-lags model. The fit gives $\chi^{2}=684.0$ for 619 dof, with structured residuals in the Imaginary part of the CS around 2 Hz and 4 Hz, close to the frequencies of the QPO fundamental and second harmonic. We next let the frequency and FWHM of each Lorentzian free but linked in the PS and CS. This fit gives $\chi^{2}=675.3$ for 607 dof, but the residuals around 2 Hz and 4 Hz in the Imaginary part of the CS remain. We show the Imaginary part of the CS and the residuals of the fit with six Lorentzians in the left panel of Figure 2.

Table 3: Phase and time lags ($5.7-100$ keV vs. $2-5.7$ keV) of the QPO fundamental and second harmonic and the corresponding QPO shoulders in GRS 1915+105 with our method

Component	$\nu_{0}$ (Hz)	Phase lags${}^{\dagger}$ (rad)	Time lags${}^{\ast}$ (ms)
QPO fund	$1.822\pm 0.006$	$-0.076\pm 0.006$	$-7.1\pm 0.7$
Should. fund	$2.11\pm 0.06$	$1.9^{\,\,+0.7}_{-1.0}$	$114\pm 20$
2nd harm	$3.56\pm 0.02$	$0.20\pm 0.2$	$11.2\pm 1$
Should. 2nd	$4.00\pm 0.03$	$0.49\pm 0.08$	$21.7\pm 5$
${}^{\dagger}$Using the constant phase-lags model.
${}^{\ast}$Using the constant time-lags model.

We therefore add two new Lorentzians to the model at, respectively, $\sim 2$ Hz and $\sim 4$ Hz with, as for the other Lorentzians, the centroid frequency and FWHM free and linked in the PS and CS. This fit, shown in the middle and right panels of Figure 2, gives $\chi^{2}=627.1$ for 598 dof, significantly better than the one with six Lorentzians (and also better than a model with seven Lorentzians not shown here), and the fit no longer shows the structured residuals at around 2 Hz and 4 Hz. In fact all Lorentzians are at least 3-$\sigma$ significant in either the PS, the CS or both (see below). We therefore adopt this one as the final model.

Figure 3: Left panel: PS of the second segment of the RXTE observation 10408-01-27-00 of GRS 1915+105 in two energy bands (see legend). The upper panel shows the data and the best-fitting model (thick solid line), consisting of five Lorentzian functions (thin solid and dashed lines). We fitted the two PS simultaneously with the centroid frequency and FWHM of all Lorentzians linked so that they are the same in both bands. The two Lorentzians plotted using thin solid blue lines represent, respectively, the main QPO at 5.8 Hz and the QPO shoulder at 6.3 Hz in the $5.1-6.9$ keV band; similarly, the Lorentzians plotted using the thin dashed red lines represent, respectively, the main QPO at 5.8 Hz and the QPO shoulder at 6.3 Hz in the $13-18.1$ keV band. The bottom panel shows the residuals with respect to the best fitting model. Middle panel: Fractional rms amplitude of the main QPO at 5.8 Hz (black) and the QPO shoulder at 6.3 Hz (purple) vs. energy. Right panel: Time lags of the main QPO and the QPO shoulder vs. energy. The reference band for the lags is the full band (channels 0-249).

The middle panel of Figure 2 shows the PS (upper panel) and the residuals (bottom panel) of the fit with the adopted model with eight Lorentzians. As is apparent from the Figure, both the QPO fundamental and the QPO second harmonic are fitted with two Lorentzians, a relatively broad one that fits the bulk of the QPO profile, and a relatively narrow one at slightly higher frequency than the other one. While both the fundamental QPO at $\nu_{\rm QPO}=1.822\pm 0.006$ Hz and the second harmonic at $\nu_{\rm 2nd}=3.56\pm 0.02$ Hz are significantly required by the fit (24.6-$\sigma$ and 6.8-$\sigma$, respectively), the shoulder of the fundamental QPO at $\nu_{\rm shoulder,1}=2.11\pm 0.06$ Hz and that of the second harmonic at $\nu_{\rm shoulder,2}=4.00\pm 0.03$ Hz are not detected significantly in the PS ($<1$-$\sigma$ and $2.6$-$\sigma$, respectively). However, both shoulders are significantly required (4-$\sigma$ and $\sim 5.5$-$\sigma$, respectively) to fit the CS. This is remarkable, since a component that is not required to fit the PS is significantly required to fit the CS. As this occurs again in §3.4, we expand on this there.

The right panel of Figure 2 shows the Imaginary part of the CS (upper panel), the phase-lag spectrum (middle panel) and the residuals of the phase-lag spectrum with respect to the derived model with eight Lorentzians (bottom panel; we do not plot the Real part of the PS because it looks very similar to the PS ). We give the frequencies and phase lags of the fundamental, second harmonic and the shoulders of the QPO fundamental and second harmonic in Table 3. The phase lags of the shoulder of the second harmonic are significantly (3.5-$\sigma$) larger than those of the second harmonic itself, but the phase lags of the shoulder of the fundamental have very large errors, therefore we cannot conclude whether the phase lags of the shoulder and the fundamental QPO are different.

Using the traditional method, the lags of the QPO, the QPO shoulder, the second harmonic and the shoulder of the second harmonic are, respectively, $-0.056\pm 0.005$ rad, $-0.001\pm 0.009$ rad, $0.123\pm 0.008$ rad and $0.16\pm 0.01$ rad. Given that the error bars in our method are relatively large compared to those in the traditional method (see Appendix A for a discussion of this), the lags of the QPO fundamental, the QPO shoulder and the second harmonic are all consistent with being the same using either the traditional or our method (see Table 3). On the contrary, the lags of the shoulder of the second harmonic using our and the traditional method are $4\sigma$ different. In Appendix A we use simulations to explore the difference in the results of our and the traditional method.

Constant time-lags model: When we fit the data assuming the constant time-lags model (this fit is marginally worse, $\chi^{2}=633.7$ for 598 dof, than the one with the constant phase-lags model), the time lags of the fundamental QPO and those of its shoulder are 6.1-$\sigma$ different, whereas the time lags of the second harmonic and its shoulder are only 2.1-$\sigma$ different. We give the time lags from this fit in Table 3.

To summarise, we find that in this observation of GRS 1915+105 either the phase or the time lags of the shoulders are significantly different from those of the corresponding QPO, and hence the shoulders are not due to small drifts of the QPO frequency during the observation, but are separate components in the PS or the CS of the source.

3.3 Case study 3: No energy dependence of the QPO frequency in GRS 1915+105

Table 4: Phase and time lags (given band vs. the full band) of the QPO and QPO shoulder in GRS 1915+105 with our method

Component	$\nu_{0}$ (Hz)	Phase lags${}^{\dagger}$ (rad)	Time lags${}^{\ast}$ (ms)	Phase lags ${}^{\dagger}$ (rad)	Time lags${}^{\ast}$ (ms)
		($5.1-6.9$ keV)		($13-18.1$ keV)
QPO	$5.78\pm 0.02$	$0.16\pm 0.03$	$4.9\pm 0.5$	$-0.61\pm 0.06$	$-20\pm 1$
QPO shoulder	$6.34\pm 0.06$	$0.24\pm 0.08$	$5.5\pm 0.6$	$-0.36\pm 0.08$	$-12.5\pm 0.8$
${}^{\dagger}$Using the constant phase-lags model.
${}^{\ast}$Using the constant time-lags model.

Now that we have established that the shoulders of the QPOs are separate components, we explore whether previous claims that the QPO frequency changes with energy (Qu et al., 2010; Li et al., 2013a; Li et al., 2013b; Yan et al., 2018) are connected to the presence of a shoulder of the QPO. For this we analyse the RXTE observation 10408-01-27-00 of GRS 1915+105 for which Qu et al. (2010) reported a significant change of the centroid frequency of the QPO with energy, from $\nu_{\rm QPO}=5.9$ Hz for energies below 5 keV to $\nu_{\rm QPO}=6.7$ Hz for energies above 20 keV (see their Figure 2). We speculate that this apparent change of the QPO frequency with energy could be due to the presence of a shoulder if the rms spectrum of the QPO and the shoulder are different. We will therefore call the QPO reported by Qu et al. (2010) the “QPO feature”, since it could be a combination of a QPO and a QPO shoulder.

We take the Binned and Event Mode data of this observation to compute the full-band (channels 0 to 249) spectrogram (also known as dynamical PS; e.g., Markwardt et al., 1999) of the observation, and find that the frequency of the QPO feature changes with time (not shown), from above $\sim 6.5$ Hz at the start to $\sim 5$ Hz at the end of the observation; the frequency of the QPO feature appears to be more or less constant at $\sim 6$ Hz in the second segment of the observation, therefore we use only that segment for the rest of the analysis. We then compute PS in the same six energy bands used by Qu et al. (2010, channels $0-13$, $14-18$, $19-25$, $26-35$, $36-49$ and $50-103$; see their Table 2 for the energy ranges) to study the rms spectrum of the QPO feature, plus the CS of the photons in the individual bands with respect to those in the full band (channels 0 to 249). Both for the PS and the CS we use $T_{\rm FFT}=128$ s, yielding $\nu_{\rm min}=\Delta\nu=1/128$ Hz, at a time resolution $\Delta t=1/512$ s that yields $\nu_{\rm Nyquist}=256$ Hz.

The left panel of Figure 3 shows the PS of the data during the second segment of the observation in two energy bands, $5.1-6.9$ keV and $13-18.1$ keV. We fit the two PS simultaneously in the range $0.5-15$ Hz with five Lorentzians, with the frequency and FWHM of each Lorentzian free but linked to be the same in both energy bands. As it is apparent from the Figure, the QPO feature is fitted with two Lorentzians, one at $\sim 5.8$ Hz and the other at $\sim 6.3$ Hz. Notably, the Lorentzian at 5.8 Hz is more or less equally strong in the two bands whereas the Lorentzian at 6.3 Hz is stronger in the $13-18.1$ keV band than in the $5.1-6.9$ keV band. (It is also true that all the other Lorentzian components in the PS are stronger in the high- than in the low-energy band.) As before, we call these two Lorentzian components the QPO and the shoulder of the QPO, respectively.

Figure 4: PS in the $0.3-12$ keV energy band (top left), Real (bottom left) and Imaginary (bottom right) part of the CS for the $2-12$ keV band with respect to the $0.3-2$ keV band of the NICER observation 1200120120 of MAXI J1820+070 fitted with a model (red solid line) consisting of four Lorentzians (dotted lines). For the fitting and the plot we rotate the cross vector by $45^{\circ}$ (see text for details), we assume the constant phase-lags model (see Table 1 and §2.1) and we fix the frequency and FWHM of the Lorentzians to the values we obtain from the best-fitting model of the PS. The top right panel shows the phase-lag spectrum with the derived model. In each case the top sub-panels show the data and the model and the middle sub-panels show the residuals. The bottom sub-panels show the residuals when we fit the data letting the frequency and FWHM of the Lorentzians free, but linked to be the same in the PS and the CS.

In the $5.1-6.9$ keV band the rms amplitude of the QPO and the shoulder are rms${}_{\rm QPO}=3.6\pm 0.3\%$ (5-$\sigma$ significant) and rms${}_{\rm shoulder}=2.4\pm 0.3\%$ (3.7-$\sigma$ significant). In the $13-18.1$ keV band the rms amplitude of the QPO and the shoulder are rms${}_{\rm QPO}=4.2\pm 0.4\%$ (4.4-$\sigma$ significant) and rms${}_{\rm shoulder}=5.1\pm 0.7\%$ (4.2-$\sigma$ significant). We give the frequencies, phase lags (obtained using the constant phase-lags models) and time lags (using the constant time-lags model) of the QPO and the shoulder in Table 4. While the phase lags ($13-18.1$ keV with respect to $5.1-6.9$ keV) of the QPO and the QPO shoulder are consistent with being the same within 3-$\sigma$, their time lags are 4.9-$\sigma$ different. (See the discussion about the difference between phase and time lags at the end of §3.)

Figure 5: Same as Figure 4 but now fitted with seven Lorentzians instead of four, assuming the constant phase-lags model. The bottom sub-panels show the residuals with respect to the best-fitting model.

The middle panel of Figure 3 shows the rms spectrum of the QPO and the QPO shoulder. From this Figure it is apparent that the rms amplitude of the QPO increases with energy up to about 10 keV and at that energy it levels off or decreases slightly, whereas the rms amplitude of the QPO shoulder increases steadily with energy, also above $\sim 10$ keV. Because of this, at energies below $\sim 8$ keV the QPO is stronger than the QPO shoulder, whereas above that energy the opposite is true. The right panel of Figure 3 shows the time-lag spectrum of the QPO and the QPO shoulder obtained using the constant time-lags model and taking the full band as reference. (The phase-lag spectrum, from the fit with the constant phase-lags model, looks almost exactly the same as this one, with slightly larger error bars.) The time lags of both the QPO and the QPO shoulder become more negative as the energy increases (compare with Fig. 2b in Qu et al., 2010, in which they plot the phase lags with respect to the lowest energy band); both time-lag spectra are consistent with being the same up to $\sim 10$ keV, at which point the time lags of the QPO continue decreasing while the time lags of the QPO shoulder level off.

The fits of the PS in two energy bands (Fig. 3, left panel), the rms spectrum (Fig. 3, middle panel) and the time-lag spectrum (Fig. 3, right panel) of the QPO and the QPO shoulder show that the data are consistent with the presence of two separate components, the QPO and the QPO shoulder, with different rms and time-lag spectra. In our model, the centroid frequency of both the QPO and the QPO shoulder are linked to be the same in the different energy bands, it is only that the change of the relative strength of the two mimics a frequency dependence of the QPO feature with energy.

In fact, our model fits the data better and with less free parameter than the model of a QPO with frequency that changes with energy. Our model of the QPO feature in $N$ energy bands with two Lorentzians that have each the same frequency and FWHM in all the bands has $4+2N$ parameters: Two centroid frequencies, two FWHM and $2N$ normalisations. The model with two separate QPOs in $N$ energy bands has $3N$ parameters: $N$ centroid frequencies, $N$ FWHM and $N$ normalisations. When we fit the PS of the six bands simultaneously with a model in which the QPO feature is a single Lorentzian with the frequency and FWHM of the Lorentzian free to change with energy, like in Qu et al. (2010), we get $\chi^{2}=601.9$ for 528 dof. The same fit with a model that has a Lorentzian for the QPO and another one for the QPO shoulder with the frequency and FWHM of each Lorentzian linked to be the same in all bands gives a better fit, $\chi^{2}=566.1$ for 530 dof.

In conclusion, in this observation of GRS 1915+105 the data are better fit with a model in which the QPO feature consists of two components with frequencies that do not change with energy rather than with a single Lorentzian with centroid frequency and FWHM that do change with energy. The rms spectra of the two Lorentzian components combine in such a way that the frequency of the QPO feature appears to depend upon energy.

3.4 Case study 4: The broadband variability in MAXI J1820+070

As a final example, in this and the next section we study the PS and CS of the transient black-hole X-ray binary MAXI J1820+070. We select this source because it was very bright and highly variable during outburst (Wang et al., 2022) and because, as we explain below, it provides a stringent test to the hypothesis that we put forward here.

From 2018 March 6 to 2018 November 21 MAXI J1820+070 was observed almost daily with NICER. During the rising part of the outburst, in the low-hard and hard-intermediate states, the source showed strong broadband X-ray variability (e.g. De Marco et al., 2021). Notably, except for a few cases (e.g., Homan et al., 2020; Ma et al., 2023a) the power spectra of most of these observations were relatively featureless (e.g., Kawamura et al., 2022) and showed no strong and narrow QPOs. Furthermore, the magnitude of the time lags was rather small and the time-lag spectra were smooth (e.g., Kara et al., 2019). All these characteristics can potentially challenge our proposal that, as for the PS, the CS consists of a combination of Lorentzian components.

Here we use ObsID 1200120120, which was previously analysed by Kara et al. (2019) and Wang et al. (2022). We process the data with the tool nicerl2 to produce clean event files; following the recommendations on the NICER website, we discard the data of detectors 14 and 34 that show episodes of increased detector noise. We compute a PS in the $0.3-12$ keV band and a CS of photons in the $2-12$ keV band with respect to those in the $0.3-2$ keV band. (In this case we do not give the channels because for NICER the channel number can be calculated directly as the energy in keV multiplied by 100.) Both for the PS and the CS we use $T_{\rm FFT}=13.1072$ s, yielding $\nu_{\rm min}=\Delta\nu=0.07629$ Hz, at a time resolution $\Delta t=1/2500$ s such that $\nu_{\rm Nyquist}=1250$ Hz.

We initially fit the PS in the range $0.07629-60$ Hz following the procedure described in §3, and find that we need a model consisting of four Lorentzians. The fit gives $\chi^{2}=153.9$ for 146 dof, and all Lorentzians are at least 3-$\sigma$ significant.

Before we discuss the joint fits of the PS and CS we note that, because the lags of MAXI J1820+070 are close to zero over a broad frequency range (e.g., Kara et al., 2019), at all Fourier frequencies the Real part of the CS is much larger than the Imaginary part. Since all fitting routines are more stable when the free parameters are of the same order (see, for instance, the Levenberg–Marquardt algorithm in the Xspec code, Arnaud, 1996), we rotated all the cross vectors by $45^{\circ}$ such that those with a zero phase lag would end up having more or less equal Real and Imaginary parts. Since the rotation does not change the modulus of the cross vector, and the rotation angle is known precisely, the rotation has no effect on the parameters of the fit, while it makes the fit more stable. Because of this, in the Figures in this section we plot the components of the rotated cross vector, $\mathrm{Re}\,\cos{(\pi/4)}-\mathrm{Im}\,\sin{(\pi/4)}\propto(\mathrm{Re}-\mathrm{Im})$ and $\mathrm{Re}\,\sin{(\pi/4)}+\mathrm{Im}\,\cos{(\pi/4)}\propto(\mathrm{Re}+\mathrm{Im})$, which also allows us to plot the components of the (rotated) cross vector vs. Fourier frequency using logarithmic axes, because both rotated quantities are always positive. Regardless of this, we always report the best-fitting phase lags minus $\pi/4$ so that the lags are given again with respect to the chosen energy band for the non-rotated cross vector.

Constant phase-lags model: The top-left and the two bottom panels of Figure 4 show the PS and the CS of MAXI J1820+070 in the range $0.07629-60$ Hz with the same four Lorentzians assuming the constant phase-lags model. During the fit of the CS, we initially fix the frequency and FWHM of the four Lorentzian components to the values obtained from the fit to the PS. The top-right panel shows the phase-lag spectrum with the derived model. In each case the top sub-panels show the data and the model and the middle sub-panels show the residuals.

It is apparent from this Figure that, while the PS is well fitted with four Lorentzians, the Real and Imaginary parts of the CS are not, with structured residuals in the full frequency range. The joint fit of the PS and the CS yields $\chi^{2}=593.7$ for 438 dof, while the derived model of the lags gives $\chi^{2}=412.1$ for 150 dof.

We next fit the model letting the frequency and the FWHM of each of the four Lorentzian components free but linked across the three spectra. The bottom sub-panels of Figure 4 show the residuals in this case. While the fits of the CS and of the derived model of the phase lags improve slightly, significant residuals remain, and this time there are also structured residuals in the fit of the PS because in the model the frequency and FWHM of the Lorentzians change. The joint fit in this case gives $\chi^{2}=529.6$ for 430 dof, while the derived model of the lags gives $\chi^{2}=417.4$ for 150 dof.

Table 5: Parameters of the best-fitting model to the PS and CS of NICER observation 1200120120 of MAXI J1820+070 with seven Lorentzians assuming the constant phase-lags model

Component	$\nu_{0}$ (Hz)	FWHM (Hz)	rms${}_{\rm PS}$ (%)	significance${}^{*}$	phase lags (rad)	rms${}_{\rm CS}$ (%)	significance${}^{*}$
Lorentzian 1	$0.054\pm 0.002$	$0.245\pm 0.002$	$25.4\pm 0.1$	$97.3$	$0.10\pm 0.01$	$26.6\pm 0.1$	$124.7$
Lorentzian 2	$0.85\pm 0.06$	$2.0\pm 0.1$	$<9.5^{\dagger}$	$<1$	$1.5\pm 0.2$	$12.5\pm 1.4$	$5.3$
Lorentzian 3	$0.99\pm 0.04$	$2.1\pm 0.09$	$11.2\pm 1.2$	$2.4$	$-0.9\pm 0.2$	$13.5\pm 1.3$	$4.0$
Lorentzian 4	$1.83\pm 0.08$	$5.68\pm 0.08$	$10.7\pm 0.3$	$13.8$	$-0.03\pm 0.02$	$13.2\pm 0.3$	$16.0$
Lorentzian 5	$7.8\pm 0.4$	$8.3\pm 1.3$	$<1.2^{\dagger}$	$<1$	$-0.91\pm 0.07$	$2.5\pm 0.3$	$5.2$
Lorentzian 6	$18.5\pm 0.8$	$15.7\pm 3.2$	$<11.2^{\dagger}$	$<1$	$-0.73\pm 0.07$	$1.9\pm 0.3$	$3.9$
Lorentzian 7	$24.1\pm 1.2$	$31.7\pm 1.6$	$1.68\pm 0.09$	$8.2$	$0.7\pm 0.2$	$1.8\pm 0.1$	$6.4$
$\chi^{2}/\mathrm{dof}$					$424.9/415$
${}^{*}$ In units of $\sigma$.
${}^{\dagger}$ 95% upper limit.
The rms amplitude of the power and cross spectra of each Lorentzian are integrated from zero to infinity.

The previous result shows that, even if the fit to the PS with a model consisting of four Lorentzians is statistically acceptable, the fits of the PS and the CS with that model are not. We therefore add extra Lorentzians to the PS and the CS until the simultaneous fit to the three spectra, and the derived model of the phase lags are statistically acceptable, even if some of the Lorentzian components are not significant in either of the three spectra, as long as they are significant in at least one of them. In doing the fits we let the centroid frequency and FWHM of the Lorentzian components free but link them so that they are the same for each component in the PS and the CS.

Figure 6: Intrinsic coherence function of the same observation of MAXI J1820+070 shown in Figures 4 and 5 with the derived model obtained from the fit to the PS and CS assuming the constant phase-lags model.

Figure 5 shows the data fitted with a model consisting of seven Lorentzians assuming the constant phase-lags model. As usual, the upper panels show the data and the model, and the lower panels show the residuals. The top-right panel shows the same for the phase-lags frequency spectrum with the derived model. The joint fit of the PS and CS gives $\chi^{2}=424.9$ for 415 dof, while the derived model of the lags gives $\chi^{2}=148.5$ for 150 dof. We give the best-fitting parameters in Table 5. Lorentzians number 2, 5 and 6 in this model are less than 3-$\sigma$ significant in the PS, but the three of them are significant in the CS.

We subsequently compute the PS of this observation of MAXI J1802+070 in the $0.3-2$ keV and $2-12$ keV bands using the same parameters mentioned at the start of this subsection, and fit them, together with the CS of those two bands, with the model consisting of seven Lorentzians over the same frequency range as before. In Figure 6 we plot the observed intrinsic coherence with the derived model for the constant phase-lags model (see Appendix B for the fit using the constant time-lags model.) The fit gives $\chi^{2}=586.9$ for 558 dof. The derived model reproduces the data rather well, with the largest deviations appearing at $0.2$ Hz, where the PS shows an unresolved QPO (see, e.g., Fig. 5) that, when studied at a higher frequency resolution (not shown), breaks into two separate QPO peaks.

Once more, similar to what we discuss in §3.1, the most appropriate model for the PS and the CS in the case of this MAXI J1820+070 observation accurately characterises both the phase-lag spectrum and the coherence function. This alignment with the two predictions we detail in §2 remains consistent.

Figure 7: PS in two energy bands (top left; see legend for the energies of each band), Real (bottom left) and Imaginary part (bottom right) of the CS for those two same bands of the NICER observation 1200120268 of MAXI J1820+070 fitted with a model (thick solid line) consisting of seven Lorentzians (dotted lines). During the fitting we let the frequency and FWHM of each component free to vary but keep them linked to be the same in the PS and the CS. In the model of the CS we assume the constant phase-lags model. The top-right panel shows the phase lags vs. Fourier frequency together with the derived model. The bottom sub-panels show the residuals with respect to the best-fitting model.

3.5 Case study 5: A QPO in the Imaginary part of the cross spectrum of MAXI J1820+070

We searched the NICER archive for an observation in which the coherence function showed a significant drop at some Fourier frequency, and came across ObsID 1200120268 (O. König, priv. comm.), which was previously analysed by De Marco et al. (2021). In this observation MAXI J1820+070 was in the decline of the outburst, in the lower branch of the ‘q’ traced by the source in the hardness-intensity diagram (see, e.g., Fig. 1 of De Marco et al., 2021, or the right panel of Fig. 1 of Ma et al. 2023a). We follow the procedures in §3.4 to process the NICER data, and compute two PS in the $0.3-2$ keV and $2-12$ keV bands and a CS of photons between those same two bands using the same $\Delta t$ and $T_{\rm FFT}$ of §3.4 to compute the FFT

Constant phase-lags model: We fit the two PS and the CS in the range $0.07629-60$ Hz with seven Lorentzians using the constant phase–lags models, which gives $\chi^{2}=505.1$ for 558 dof. All Lorentzians are at least 3-$\sigma$ significant either in one of the PS or the CS. Figure 7 shows the data and the best fit. While we again rotate the cross vector by $45^{\circ}$ during the fits, this time we plot the Real and Imaginary parts of the CS without rotation, to highlight the presence of a strong QPO at $\sim 2.1$ Hz in the Imaginary part of the CS. The phase lag of this QPO is $\Delta\phi=1.06\pm 0.12$ rad, such that the Imaginary part of the cross vector is about twice larger than the Real part. The 2.1-Hz QPO is not significant in either of the two power spectra (see below). Close to the 2.1-Hz QPO, the Imaginary part of the CS shows also QPOs at $\sim 2.6$ Hz and $\sim 4.3$ Hz, but these QPOs have smaller phase lags ($\sim 0.58$ rad and $\sim 0.36$ rad, respectively) and both are significantly detected in the PS of the two bands.

In Figure 8 we plot the observed intrinsic coherence together with the model derived from the simultaneous fit of the two PS and the CS. The observed intrinsic coherence shows a significant drop at $\sim 2.1$ Hz which, as shown in the plot, is perfectly described by the derived model. It is noteworthy that this drop occurs precisely at the frequency of a QPO that is present only in the Imaginary part of the CS (bottom-right panel of Figure 7). This “imaginary” QPO, at $2.09\pm 0.02$ Hz, is marginally present in the Real part of the CS (bottom left), where it is overshadowed by several other components, and is not significantly detected in any of the two PS (top left; the QPO is $\lesssim 2.9\sigma$ significant in the $0.3-2$ keV band and $\lesssim 2.5\sigma$ in the $2-12$ keV band). This underscores the fact that this QPO is only detected significantly in the Imaginary part of the CS, and that the Lorentzian component used to fit this imaginary QPO is incoherent with the other Lorentzians in that frequency range. This aligns with the assumptions made in Section §2, further bolstering our assertion that the variability is comprised of multiple incoherent components.

Figure 8: Intrinsic coherence function of the same observation of MAXI J1820+070 shown in Figure 7 with the derived model obtained from the fit to the PS and CS assuming the constant phase-lags model.

4 Discussion

We show, for the first time, that if the power spectrum (PS) of accreting neutron-stars and black-holes can be fitted with a combination of Lorentzian functions that are coherent in different energy bands but incoherent with each other, the same is true for the Real and Imaginary parts of the cross spectrum (CS). Based on this, we propose a novel method to measure the lags of variability components that is especially useful when those components are weak and overlap with other (stronger) components in the PS and CS.

Surprisingly, using this method we discover new variability components that are detected significantly only in the CS and not in the PS of these sources. This happens because the PS is insensitive to signals with a large Imaginary part and a small Real part in the CS when such a signal overlaps in frequency with other variability components that have a large Real part in the CS. Because in the last 40 years we have used exclusively the PS to identify variability components in these sources, we have so far missed signals with large positive or negative lags, such that the cross vector has a significant component along the Imaginary axis.

We also show that, contrary to what has been previously claimed, the frequency of a type-C quasi-periodic oscillation (QPO) in the black-hole binary GRS 1915+105 does not depend upon energy. The apparent energy dependence of the QPO frequency can be explained by the presence of a second, significant, component in the PS and CS of the source at a frequency very close to that of the QPO, but with a different rms and lag spectra. This alternative interpretation requires fewer model parameters and is statistically favoured. We propose that the same applies to other QPOs in which a similar energy dependence of the QPO frequency was observed.

Finally, we demonstrate that, in accordance with the predictions derived from the assumptions underpinning our method, the model that fits both the PS and CS reproduces both the phase-lag spectrum and the coherence function.

4.1 Weak variability components in the presence of other, stronger, components

In §2 we show that if the PS of X-ray binaries can be fitted with a linear combination of Lorentzian functions that are coherent in different energy bands and incoherent with each other, the same is true for the CS of linearly correlated light curves of the source in two energy bands. The centroid frequency and FWHM of each Lorentzian are the same in the PS and CS, with the Lorentzians in the Real (Imaginary) part of the CS being multiplied by the cosine (sine) of a function of Fourier frequency, $\Delta\phi_{xy}(\nu)=g(\nu;p_{j})$ with parameters $p_{j}$. These functions are the frequency-dependent phase lags of each Lorentzian. If $g(\nu)$ is constant, the phase lags of that component are independent of Fourier frequency, while if $g(\nu)\propto\nu$ the time lags are independent of Fourier frequency.

While we show that, under the assumptions discussed in §2, the same Lorentzians should appear both in the PS and the CS, there is no guarantee that all the Lorentzians that are significantly detected in one will be significantly detected in the other. To test this, and to explore the potential of our method, in §3 we fit the PS and CS of five observations of three black-hole binaries. We unveil weak QPO signals at frequencies very close to those of previously detected QPOs in the black-hole binaries GX 339–4 (Zhang et al., 2017; Altamirano & Méndez, 2015) and GRS 1915+105 (Zhang et al., 2020). The weak signals, which we call QPO shoulders, appear at a slightly higher frequency than that of the QPO fundamental or second harmonic. QPO shoulders have been previously reported (Belloni et al., 1997; Belloni et al., 2002; van Doesburgh & van der Klis, 2020), but until now it remained unclear if those were really separate QPOs or whether those shoulders were due to the QPO frequency drifting slightly during the time over which one calculated the PS (see, e.g., Belloni et al., 1997). Our method allows us, for the first time, to measure a significant difference of the lags and rms spectra of the QPO and the QPO shoulder, which strongly suggests that the QPO and the shoulder are truly different components. In GRS 1915+105 the shoulder of the QPO fundamental and the second harmonic are also consistent with having a harmonic relation. It is therefore possible that all harmonics of a QPO, if studied using this method, will eventually show shoulders that, as the QPOs, are also harmonically related.

An advantage of the method that we propose here is that it allows to measure the lags of components that appear in a frequency region of the PS or CS where other equally strong or stronger components are present. One cannot use the traditional method to measure the lags of such a component in those cases (see §3.1). A way used in the literature to overcome this limitation is to measure the lags in the traditional way within one FWHM of the component of interest and subtract the lags measured outside that range (Ma et al., 2021, 2023b). This method is mathematically incorrect because it ignores the modulus of the cross vector (CV) of the two components. Alternatively, one could try and subtract the Real and Imaginary parts of the CS in a frequency range range outside the QPO from the same quantities within a FWHM of the QPO, however this is rarely possible (see for instance Fig. 1 above or Fig. 1 in Nowak, 2000).

4.2 Hidden variability components

Initially, we assumed that we had to fit the PS with a number of Lorentzians first, and subsequently fit the CS with the same Lorentzians with centroid frequencies and FHWM fixed at the values obtained from the fit to the PS. While we find a good fit to the PS with a number of Lorentzians, when we fit the CS with those same Lorentzians with only the normalisation of the Lorentzians free the fit is bad, with significant structured residuals over the full frequency range. This is most clearly seen in the case of the observation of the black-hole binary MAXI J1820+070 that we present in §3.4. As shown in Figure 4 for the constant phase-lags model and Figure 9 for the constant time-lags model, while the PS is well fitted with four Lorentzians, the CS is not; because of this, the phase-lag frequency spectrum cannot be fitted with the derived model either (cf. the residuals in the middle panels of those two Figures). If we let the frequency and FWHM of the Lorentzians free but link each of them to be the same in the PS and CS, the fit improves but remains statistically unacceptable, with significant structured residuals in the CS and the lags. Not only that, but now there are significant residuals also in the PS. The reason for this is easy to understand: In the joint fits the frequency and FWHM of the Lorentzians change to reduce the residuals in the CS, but this degrades the fit of the PS.

In the case of MAXI J1820+070, in order to get a good fit we need to add three extra Lorentzians to the model. While some of these Lorentzians are not significant in the PS, they are (in some cases very) significant in the CS. This can be seen in Figures 5 and 10 and Tables 5 and 6. For instance, Lorentzians 2, 3 and 6 in the constant phase-lags model (Figure 5 and Table 5) and Lorentzian 5 in the constant time-lags model (Figure 10 and Table 6) are less than 1-$\sigma$ significant in the PS, but they are $3.9-6.5\sigma$ significant in the CS.

To understand why we do not see these components in the PS, whereas they are significant in the CS, we first need to note that the majority of the variability components in LMXBs have rather small lags (see, for instance, Fig. 2 in Zhang et al., 2020, the same is true in other sources). This means that in the CS the Real part is much larger than Imaginary part (factors $\gtrsim 10$ are usual; see for instance Figs. 1 and 7). Consider two variability components, $X_{1}$ and $X_{2}$, that overlap in frequency, each of them perfectly coherent in two energy bands (c.f., assumption (i) in §2). Furthermore, consider that $X_{1}$ is weak with phase lags close to $\pm\pi/2$ (Real part in the CS close to zero and small but significantly different from zero Imaginary part) and $X_{2}$ is strong with phase lags close to zero (in the CS large Real part, Imaginary part close to zero). It is easy to show that if the two components have similar rms spectra, the contribution of component $X_{1}$ in the PS is negligible compared to that of component $X_{2}$. On the other hand, the statistical errors in the power spectrum will be of the order of the Real part of CS of $X_{2}$ divided by the square root of the product of the number of segments and the number of frequency bins (Bendat & Piersol, 2010; Ingram, 2019). As long as this error is larger than the Imaginary part of component $X_{1}$, the latter will not be detected in the power spectrum, even if it is very significant in the Imaginary part of the CS. We note that the advantage of detecting such a weak component in the cross spectrum, and not in the power spectrum, will not be available when, instead of the “contaminating” component(s), the Poisson noise dominates the variability.

An example of this is the $2.1$-Hz QPO in MAXI J1820+070 that has a phase lag of $\sim 1.1$ rad, appears close to a strong broad QPO peaking at $\sim 0.6$ Hz (Fig. 7) with a phase lag of $\sim-0.01$ rad, and is only detected in the Imaginary part of the CS. (A similar effect is likely the cause of the drop of the coherence function at $\sim 0.2-0.3$ Hz, in the observation of GRS 1915+105 shown in Fig. 3 of Ji et al. 2003, where no significant QPO is apparent in the PS.)

All the above demonstrates that by searching for QPOs exclusively in the PS we tend to find only components with relatively small lags, and we are insensitive to signals with lags that approach $\pm\pi/2$. In fact, we would never find lags of $\pm\pi/2$ rad in the CS of any of these sources if we measure the lags over a broad frequency band that contains several components with small lags, because the Real part of the CS in that frequency range will never be zero. By using only the PS to identify and characterise the variability of these sources, and by computing the average lags over broad frequency ranges, all these years we may have missed these components altogether.

4.3 No energy dependence of the frequency of QPOs

It has been argued that the frequency of some low-frequency QPOs in black-hole X-ray binaries, especially the type-C QPO in GRS 1915+105, depends on energy (Qu et al., 2010; Li et al., 2013a; Li et al., 2013b; Yan et al., 2018); in some cases the QPO frequency appears to increase, and in others to decrease, with energy (see, e.g., Fig.2 of Qu et al., 2010, and Fig. 3 of Li et al. 2013a). This has been explained by van den Eijnden et al. (2016) in terms of differential Lense-Thirring precession (Ingram et al., 2016; Nathan et al., 2022), with parts of the inner accretion flow precessing at different rates. An issue with this interpretation is that sometimes the QPO frequency appears to increase whereas other times it appears to decrease with energy (Qu et al., 2010). Since the characteristic frequencies in the vicinity of a compact object as well as the temperature of the material in the accretion flow generally decrease with increasing distance from the central object, this requires that the outer and cooler parts of the accretion flow precess faster than the inner and hotter parts.

We find that in one of the observations of GRS 1915+105 in which Qu et al. (2010) measure an increase of the QPO frequency with energy, the QPO feature can actually be fitted better with a model that consists of two components, a QPO at $\sim 5.8$ Hz and a QPO shoulder at $6.3$ Hz. Because the rms amplitude of the shoulder increases faster with energy than that of the QPO, the QPO dominates at low energies whereas the QPO shoulder dominates at high energies, mimicking a shift of the frequency of the QPO feature with energy. The fit with our model with two Lorentzian components with energy-independent centroid frequency and FWHM is statistically better and requires less parameters than the one with a single Lorentzian with energy-dependent centroid frequency and FWHM. On the other hand, the rms and time-lag spectra of the QPO and the shoulder are significantly different, showing that the shoulder is not due to a drift of the QPO frequency during the observation, but that it is a separate component, and is the reason for the apparent energy-dependence of the QPO frequency. We propose that this holds for other cases in which the QPO frequency appears to change with energy (e.g., Li et al., 2013a; Li et al., 2013b; Yan et al., 2018), significantly challenging the ideas that explain the QPO as a single component with its frequency in different bands originating in different parts of the accretion flow.

The shoulders could be due to amplitude modulation of the QPO signal (see Fig. 3 in van den Eijnden et al., 2016). Amplitude modulation, however, should produce equally strong sidebands at frequencies below and above that of the QPO, whereas there is no apparent shoulder at frequencies below those of the QPO in the data that we present here. (The 95% confidence upper limit for the full-band rms amplitude of a possible shoulder at $\sim 0.5$ Hz below the frequency of the QPO is $\sim 1$%.) The shoulders appearing only at higher frequencies than those of the corresponding QPOs resemble the positive sidebands of the kHz QPOs observed in some neutron-star systems (Jonker et al., 2000, 2005) although, different from the case of the kHz QPOs, when the frequency of type-C QPOs appears to decrease with energy (e.g., Qu et al., 2010) the shoulder would be at lower frequencies than those of the QPOs.

In this observation of GRS 1915+105, and in the observation of GX 339–4 that we present in §3.1, the shoulder appears at a slightly higher frequency than that of the corresponding QPO. In both sources we find a shoulder of the QPO fundamental, while in the observation of GRS 1915+105 in §3.2 we also find a shoulder of the second harmonic of the QPO, with a centroid frequency that is consistent with being twice that of the shoulder of the fundamental QPO. This suggests that the shoulders of the QPO fundamental and second harmonic are physically connected.

4.4 Broadband noise or QPOs?

We show that in the black-hole binary MAXI J1820+070, while the PS can be fitted with four broad Lorentzians, the joint fits of the PS and the CS require at least seven, narrower, Lorentzian functions. This result evinces that, as it was already shown for the PS of GX339–4 and Cyg X–1 (Nowak, 2000) and several other black-hole and neutron-star X-ray binaries (Belloni et al., 2002), in MAXI J1820+070 the broadband noise (BBN) in the CS is consistent with a combination of several Lorentzians.

That the CS, like the PS, comprises multiple Lorentzian components, naturally explains the findings of Nowak et al. (1999a); Nowak et al. (1999b). They observe that in Cyg X–1 and GX 339–4, in the frequency range in which one Lorentzian dominates the PS, the phase lags of light curves in two energy bands display a more or less flat shelf. Conversely, when two Lorentzians intersect in the PS, a transition occurs from one characteristic phase-lags shelf to another (compare Figs. 1 and 10 in Nowak et al., 1999a). Nowak et al. (1999b) show that, at each Fourier frequency, the total measured phase lags are the average of the phase lags of the individual Lorentzians, weighted by the product of the amplitudes of the Fourier transforms in those two energy bands. Because over a given frequency range the Fourier amplitudes of one Lorentzian dominate, the source shows a more or less constant shelf if the phase lags of that component are constant with Fourier frequency.

All the above suggests that in MAXI J1820+070 the lags of the BBN component in the CS arise from a transfer function that can be described as the combination of several individual, narrower, responses. This idea is reinforced by the fact that the same linear combination of Lorentzian functions that fits the PS and CS predicts correctly the coherence function, in particular that the coherence is one when a single Lorentzian component dominates in both the PS and CS, and drops when two or more Lorentzians with different cross amplitudes and phase lags overlap in frequency (see Fig. 8 for a remarkable case of this effect). This, in turn, reflects in the fact that each Lorentzian has its own time or phase lags, unrelated to the lags of all other Lorentzians, possibly indicating separate resonances in the variability properties of the accretion flow (e.g., Nowak et al., 1997; Méndez et al., 2013; Méndez et al., 2022; Zhang et al., 2022b). Our findings challenge the idea that there is a smooth global transfer function of the accreting system, like the one used to explain the broadband lags in LMXBs as reverberation of corona photons that reflect off the accretion disc.

Several papers in the literature measure the time lags in MAXI J1820+070 in various frequency bands in the $\sim 0.1-80$ Hz range (e.g., Kara et al., 2019; De Marco et al., 2015; De Marco et al., 2017; Wang et al., 2022), which they interpret as the time lags of the BBN produced by the (global) transfer function of the accretion disc. This is generally justified by the fact that the PS of this source appears to be smooth, without very strong and narrow QPO peaks, and in that frequency range it can be fitted by two or three broad Lorentzian functions (e.g., Kawamura et al., 2022, 2023). Our results above, however, bring this procedure into question. We show that to fit simultaneously the PS and CS one needs seven Lorentzians. Figures 5 and 10 show that there is no unique component that dominates the PS and CS in the full $0.1-80$ Hz range, but different Lorentzians dominate the variability over different parts of that frequency range. Moreover, similar to what we show in §3.1 and Appendix A, in this case the measurements of the lags using the traditional method will clearly misrepresent the actual lags of any of these components. Figures 5 and 10 manifestly show that there is no frequency interval over which one can measure the lags of the putative BBN, without being affected by other components.

In short, our findings challenge the validity of conclusions about the geometry of these systems deduced from converting the time lags of the BBN into light travel distances of photons in the accretion flow (e.g., Kara et al., 2019; Wang et al., 2022) without a model of the frequency-dependent part of the transfer function of the system that properly accounts for the response of the individual components that comprise the BBN (e.g., Mastroserio et al., 2018; Ingram et al., 2019).

4.5 Constant phase lags or constant time lags?

As we explain in §2.1, to be able to fit the CS we need to assume the form of the functions $\Delta\phi_{xy,i}(\nu)=g_{i}(\nu;p_{j,i})$ that represent the phase-lags vs frequency of each Lorentzian component, $i$, and fit the model to get the parameters $p_{j,i}$. In this paper we assume that either the phase lags or the time lags of each Lorentzian component are constant with Fourier frequency, $g_{i}(\nu;k_{i})=2\pi k_{i}$, or $g_{i}(\nu;k_{i})=2\pi k_{i}\nu$, respectively, with $k_{i}$ constants to be determined from the fits. Also, while in principle the phase lags of different Lorentzian components could depend differently upon Fourier frequency, here we assume that the same function applies to all the components in the model.

In §3 we fit both types of models to the data and find that they are statistically equivalent. We note, however, that the parameters of some of the Lorentzians are different in one model and the other. This can be seen from a comparison of Figures 5 and 10 and Tables 5 and 6. This is not surprising, given that the form of $g(\nu;p_{j})$ affects the shape of the profile of the Lorentzian in the CS. In fact, in the constant phase-lags model the components in the CS are also Lorentzian functions, but for the constant time-lags model the components in the Real and Imaginary parts of the CS are Lorentzian functions multiplied by, respectively, $\cos{(2\pi k_{i}\nu)}$ and $\sin{(2\pi k_{i}\nu)}$. That in this case these are not Lorentzian functions is apparent in Figure 10, where some of the individual components, and the total model, show oscillations as a function of Fourier frequency with a period of $1/k_{i}$ Hz.

We initially hoped that from these data, especially those of MAXI J1820+070, we would be able to conclude whether the constant phase-lags or the constant time-lags model was most likely. Since in the constant time-lags model the CV of each component rotates in the Fourier plane as a function of Fourier frequency, each time the CV rotates by $2\pi$ radians the phase lags will wrap back to the interval $\left[-\pi,\pi\right)$ and the Real and Imaginary parts of the CS will start a new cycle. We expected that, if this was the case, we could observe in the data the oscillations predicted by the model, or that we could discard that model if the residuals showed significant oscillations with Fourier frequency. Unfortunately the current data do not allow us to choose between the two alternatives.

We stress that, while we can fit the data equally well using either the constant phase- or time-lags model for each component, these two are the simplest of an infinite choice of lag models. It remains to be seen whether a more sophisticated analysis of the data, e.g., combining several observations, would allow us to discard either of the two models, or whether a more complex model is needed to describe the data.

The models that have been proposed so far to explain the lags (and in some cases also the rms amplitude) of the X-ray variability can be roughly divided into two classes: (i) Models of the broadband variability consider that the lags come from the travel time of photons or accretion-rate fluctuations, and hence the time lags are constant with Fourier frequency. (ii) Models of the QPOs consider a sinusoidal signal (sometimes with harmonics), and hence the distinction cannot be made because, for a given frequency of the sinusoidal function, the phase and time lags are simply related to each other. The first class of models includes reverberation, propagating mass accretion rate fluctuations and Comptonisation in a relativistic jet or outflow. The second class includes Lense-Thirring precession, time-dependent Comptonisation with feedback and a precessing jet.

In the RELTRANS model (Mastroserio et al., 2018, 2019; Mastroserio et al., 2021; Ingram et al., 2019), hard photons emitted from the corona (assumed to be a point source along the spin axis of the black hole; see Lucchini et al., 2023, for the case of a corona consisting of two point sources) reach the observer first. Hard photons that illuminate and reflect off the accretion disc are reprocessed and re-emitted at lower energies than the corona photons, and reach the observer at later times. In this model the soft lags of the BBN component in the PS of these sources reflect the difference of the travel times between the hard and soft photons. (The observed time lags are not simply the difference of those travel times –see above for the shortcomings of that–, but the model computes the lags from the energy- and frequency-dependent transfer function of the disc illuminated by a point source, the so-called lamppost, at a certain height above the disc; see, e.g., Reynolds et al. 1999.)

In the several variants of the model of propagating mass accretion rate fluctuations, PROPFLUC (Rapisarda et al., 2014, 2017a, 2017b, see also Ingram & Done 2011; Mahmoud & Done 2018; Mahmoud et al. 2019; Kawamura et al. 2023; Mummery 2023), the hard lags of the BBN reflect the speed at which those fluctuations propagate in the accretion disc on the viscous time scale. Similarly, in the JED-SAD model (Ferreira et al., 2006; Marcel et al., 2019, 2020) the hard time lags would reflect the viscous time scale in the disc, given that the energy carried by mass accretion rate fluctuations will dissipate at the transition radius between the standard and the jet-dominated parts of the accretion disc (Ferreira et al., 2022).

In the model of Comptonisation in a relativistic jet/outflow (Reig et al., 2003; Giannios et al., 2004; Reig & Kylafis, 2015, 2021; Reig et al., 2018; Kylafis & Reig, 2018; Kylafis et al., 2020) the hard lags of the BBN reflect the difference in the travel time of the photons emitted from the disc that reach the observer directly and the disc photons that are inverse-Compton scattered in the jet/outflow before reaching the observer.

The Lense-Thirring precession model of the QPOs (Ingram et al., 2009, 2016; Nathan et al., 2022) assumes that, as it precesses, a hot torus located inside the inner radius of the accretion disc illuminates the surface of the disc asymmetrically and produces the variability at the QPO frequency. This precession leads to a soft delay between the modulated flux of the hard component from the precessing torus and the emission from the disc (especially the iron line).

In the time-dependent Comptonisation model, VKOMPTH (Karpouzas et al., 2020; García et al., 2021; Bellavita et al., 2022, see also Karpouzas et al. 2021; García et al. 2022; Zhang et al. 2022a, 2023a; Rawat et al. 2023; Ma et al. 2023a; Rout et al. 2023; Zhang et al. 2023b), the lags of the QPO are either hard or soft depending on the “net” delay between the photons from the disc and those that are inverse-Compton scattered in the corona, or the direct corona photons and those photons that return and are re-processed and re-emitted by the disc.

Ma et al. (2021, see also ) assume that the lags of the QPO are produced in a small-scale precessing jet. In this case, the energy of the emitted photons decreases with height in the jet, and the soft lags at the QPO frequency represent the difference of the time of arrival of photons produced at different heights of the jet.

While conceiving mechanisms that produce time delays is relatively straightforward, mechanisms that produce phase delays are less obvious. A possible phase delay is mentioned in a study of the type-C QPO in GRS 1915+105 (Lin et al., 2000), where the authors discuss the phase delay between the fundamental and the second harmonic of the QPO. Constant phase lags for each variability process, without discussing the mechanism that produces them, is also considered by Nowak et al. (1999b) to explain the phase-lag shelves and the coherence function of GX 339–4.

4.6 How should we search for variability?

In this paper we show that some variability components in LMXBs are not detected significantly in the PS, and that one needs to resort to the CS to detect them. Our findings prompt the question of the most effective method for detecting variability in these sources.

The significance of a Lorentzian component is proportional to the rms amplitude of the Lorentzian squared times the source count rate (van der Klis, 1989a). The rms amplitude of most variability components increases with energy (e.g. Berger et al., 1996; Méndez et al., 2001; Gierliński & Zdziarski, 2005; Zdziarski et al., 2005; Zhang et al., 2020), whereas the observed count rate decreases at high energies, both because of the shape of the spectrum of these sources and the drop of sensitivity of X-ray detectors. At low energies the source count rate usually decreases also because the effective area of X-ray detectors drops and because of the effect of the interstellar absorption. Based on all this, and although it is not possible to give a procedure that fits all possible cases, a generic approach would be something along these lines:

1. Fit the PS in the broadest possible energy band with a model consisting of a number of Lorentzian functions to identify the strongest variability components. One could also fit the PS in several narrow bands over the full energy range of the detector, because some components may be more significant in the full band while others may be more significant in some of the narrow bands.

2. Produce CS of two bands covering, respectively, the low- and the high-energy parts of the energy spectrum. While there is no general rule to select the bands, one could consider dividing the energy spectrum such that each band has more or less the same count rate and the combined bands cover the full band of the detector.

3. Fit the full band PS and the two-bands CS with the same Lorentzian model that fits the PS, with the frequency and FWHM of each Lorentzian free to vary but linked in the PS and CS. The Lorentzian functions in the Real and Imaginary part of the CS need to be multiplied by, respectively, $\cos{\left[g_{i}\left(\nu;p_{j,i}\right)\right]}$ and $\sin{\left[g_{i}\left(\nu;p_{j,i}\right)\right]}$ assuming a model of the phase lags to get the parameters $p_{j,i}$. If necessary add new Lorentzians to the model.

4. Repeat the previous step for (some of) the narrow-band PS and the CS of (some of) the narrow bands with respect to the full band .

5 Conclusions

We propose a new method to measure the energy-dependent phase and time lags in X-ray binaries, and use it to fit the power and cross spectra (respectively, PS and CS) of a number of sources. As we show in §3, this procedure is capable of unveiling signals that are not significantly detected in the PS alone, opening up new possibilities to understand the properties of the accretion flow around neutron-star and black-hole X-ray binaries. Specifically:

• We show that, mathematically, the PS and CS of X-ray binaries can be fitted with the same combination of Lorentzian functions, with each Lorentzian in the CS multiplied by the cosine (sine) of a function of Fourier frequency, $\Delta\phi(\nu)_{xy,i}=g_{i}(\nu;p_{j,i})$ with parameters $p_{j,i}$, with the centroid frequency and FWHM of each of the Lorentzians being the same in the PS and CS.

• We successfully fit the PS and the CS of the black-hole binaries GX 339–4, GRS 1915+105 and MAXI J1820+070 assuming that, for each Lorentzian $i$ and parameters $k_{i}$, either $g_{i}(\nu;k_{i})=2\pi k_{i}$ (constant phase lags) or $g_{i}(\nu;k_{i})=2\pi k_{i}\nu$ (constant time lags).

• We find that there is a significant shoulder at a slightly higher frequency than that of (i) the fundamental of the type-C QPO in, respectively, GX 339–4 and GRS1915+105, and of (ii) the second harmonic of the type-C QPO in GRS 1915+105. The shoulder is sometimes significant only in the CS, and would go undetected if one only analysed the PS of these sources.

• Contrary to previous reports, the frequency of the type-C QPO in an observation of GRS 1915+105 is consistent with being independent of energy. The apparent change of the QPO frequency with energy that was previously reported can be explained by a QPO feature that consists of two components, the QPO and a QPO shoulder, each of them with an rms amplitude spectrum that depends differently upon energy. This model is statistically better and requires less parameters than the model in which the QPO frequency depends upon energy.

• While the PS of an observation of the black-hole binary MAXI J1820+070 can be fitted with four Lorentzians, the simultaneous fit to the PS and CS requires seven Lorentzians, all of them significant.

• We find a narrow QPO in an observation of MAXI J1820+070 that is very significant in the Imaginary part of the CS, but is not significantly detected in the PS. This “imaginary” QPO causes a sharp drop in the coherence function at the QPO frequency.

• All the above shows that, by measuring the lags in the traditional way, one is bound to miss variability components for which the Real part of the CS is small compared to the Real part of other components in the same frequency interval.

• We argue that because of this, and because so far we have only used the PS to detect variability components, we have missed signals with large positive or negative lags a significant component of the cross vector along the Imaginary axis) in these sources.

• We conclude that, as it was previously shown for the PS, in the CS of X-ray binaries the so-called broadband noise component is in fact the combination of individual Lorentzian functions with more or less well-defined time scales. The transfer function of models used to fit the broadband variability in X-ray binaries must account for this.

Data Availability

This research has made use of data obtained from the High Energy Astrophysics Science Archive Research Center (HEASARC), provided by NASA’s Goddard Space Flight Center.

Appendix A Simulations of phase lags of the QPO: The traditional vs. our method

In this subsection we carry out simulations to compare the results that we obtain with the new and the traditional methods as a function of the strength of the QPO relative to the strength of other variability components. We note that for these simulations we assume that the PS and CS consist of a linear combination of Lorentzian functions, as described in §2. Here we only discuss the case of the constant phase-lag model, but we find similar results using the constant time-lags model.

We take the best-fitting model to the PS and CS of the observation of GRS 1915+105 in §3.2 as the basis for our simulation. To make the comparison more straightforward, we simulate the data considering a simplified version of the model shown in Figure 2; in particular, we ignore the weak shoulders of the QPO fundamental at $\sim 2$ Hz and of the second harmonic at $\sim 4$ Hz. We further assume that the QPO has a phase lag of $-0.10$ rad, and the low- and high-frequency broad Lorentzians that contribute to the variability at the QPO frequency (these are the Lorentzians in the middle panel of Figure 2 peaking at, respectively, $\sim 1$ Hz and $\sim 3.5$ Hz) have lags of $-0.20$ rad and $+0.30$ rad, respectively. Finally, for the simulations we take the normalisation of the Lorentzian components in the PS and CS that we obtain from the best-fitting model. Once we have established the base model, we carry out Monte Carlo simulations of the power and cross spectra for a range of values of $A$, the normalisation of the QPO in the PS of band 1 (channels $0-13$), such that the ratio of the QPO power to the combined power of the other two components in that band, integrated over one FWHM of the QPO, is in the range $\sim 0.25-35$. We then do the same to simulate the PS of band 2, applying the same factor, in the range $\sim 0.25-35$, to the normalisation of the QPO, $B$, in band 2 (channels $14-249$) and compute the normalisation of the cross vector, $C=\sqrt{AB}$ (see §2). For the simulations we assume that the powers in the PS and the amplitude squared in the CS are normally distributed. This is a reasonable approximation given that the simulated data are based on the average of at least 25 separate realisations of those Fourier quantities, which are also rebinned logaithmically in frequency (see §3).

We find that when the power of the QPO is $\sim 25$% of the power of the broad underlying components, the QPO is $4.2\sigma$ and $4.4\sigma$ significant in, respectively, band 1 and band 2; under these circumstances our method gives a phase lag of the QPO of $-0.047\pm 0.068$ rad, already consistent with the input value of the simulation of $-0.10$ rad. On the contrary, the lags obtained with the traditional method are $0.165\pm 0.011$ rad, more than $20\sigma$ different from the input value of the simulation, but close to the lags of the underlying high-frequency broad Lorentzian. Even when the QPO is $\sim 8$ times stronger than the broad underlying components, the lags with the traditional method are $-0.068\pm 0.005$, more than $6\sigma$ different from the input value of the simulation, whereas the lags with our method are $-0.104\pm 0.006$. Only when the QPO is $\sim 25$ times stronger than the broad underlying components, the traditional lags of the QPO are $-0.091\pm 0.003$ rad, consistent (albeit just at the $3\sigma$ level) with the input value of the simulation. This last result shows emphatically that, if the PS and CS can be decomposed into a number of Lorentzian functions, and a QPO appears in a frequency range in which there are other variability components, the traditional lags of the QPO are biased towards the lags of those components, even when the QPO is significantly stronger than the other variability component.

It is worth noting that, while the values of the lags in the traditional method are significantly different from the input value of the simulation, the errors of the traditional lags are much smaller than those obtained with our method. This is easy to understand from the fact that in our method the errors of the lags of each Lorentzian reflect the uncertainties of the parameters of all the other Lorentzians in the model. Furthermore, each individual Lorentzian component is weaker, and hence less significant, than the sum of all components, which is what one measures when one averages the CS over the FWHM of a QPO. Because of this, if the PS and CS consist truly of the sum of individual Lorentzian components, the traditional method underestimates the relative errors of the Real and Imaginary parts of the CS and hence of the lags (see also Peirano & Méndez, 2022; Alabarta et al., 2022).

All the above shows that in cases in which a QPO appears on top of a broad component in the PS, while the lags obtained using the traditional method appear to be very precise, they are very inaccurate. Both the inaccuracy and the apparent high precision of the measurements are an artefact of the contamination of the underlying components. On the contrary, while the lags obtained using our method have larger errors than those of the traditional lags, they are significantly more accurate and should therefore be used instead of the traditional lags.

Appendix B The broadband variability in MAXI J1820+070 fitted with the constant time-lags model

Figure 9: Same as Figure 4 but now assuming the constant time-lags model (see Table 1 and §2.1). Notice the difference in the scale of the residuals plot of the lags in this Figure and in Figure 4.

Constant time-lags model: As we did in §3.4, in Figure 9 we show the fit to the PS and the CS with four Lorentzians, assuming the constant time-lags model (the top left and the two bottom panels). The panel at the top right shows the phase-lag frequency spectrum with the derived model. As in Figure 4, the top sub-panels show the data and the model and the middle sub-panels show the residuals in the case in which we first fit the PS and then fix the frequency and FWHM of each Lorentzian to the values that we obtain from the PS. The joint fit to the PS and CS with this model is very bad, with $\chi^{2}=724.9$ for 438 dof, while the derived model of the lags gives $\chi^{2}=2704.5$ for 150 dof. The bottom sub-panels show the residuals when we let the frequency and the FWHM of each of the four Lorentzian components free but linked across the three spectra. In this case the fit gives $\chi^{2}=587.5$ for 430 dof, while the the derived model of the lags gives $\chi^{2}=1686.6$ for 150 dof.

Figure 10: Same as Figure 5 but now assuming the constant time-lags model.

Figure 11: Intrinsic coherence function of the same observation of MAXI J1820+070. The plot is the same as that in Figure 6, except that here we assume the constant time-lags model (right panel).

As in the case of the constant phase-lags models (§3.4, Fig. 4), when the frequency and FWHM of the Lorentzian components are fixed to the values obtained from the PS, the fit to the CS and the derived model for the phase lags show large residuals (significantly larger than those in the case of the constant phase-lags model).

Figure 10 shows the data fitted with a model consisting of seven Lorentzians assuming the constant time-lags model. The panels show the same information as the panels in Figure 5. The fit gives $\chi^{2}=417$ for 415 dof, marginally better than in the case of the constant phase-lags model. The derived model of the lags gives $\chi^{2}=205$ for 150 dof, worse than for the constant phase-lags model.

In this case, some of the components have centroid frequencies and time lags such that the phase lags wrap as explained in §2.1. For instance, the oscillations in the bottom panels of Figure 10 are due to Lorentzians 1, 2 and 3 with centroid frequencies of $\sim 0.05$ Hz, $\sim 0.165$ Hz and $\sim 0.29$ Hz that have, respectively, time lags of $\sim 242$ ms, $\sim 73$ ms and $\sim 45$ ms. For instance, the phase lags of the first of those Lorentzians rotate at a rate of $\sim 1.5$ rad per Hz such that, as explained in §2.1, the phase angle of the CV wraps $\sim 14$ times over the $0.07629-60$ Hz frequency range leading to the oscillations seen in the Figure.

We report the best-fitting parameters in Table 6. Notice that the frequency and FWHM of some of the Lorentzians are different in this case compared to the fit with the constant phase-lags model. This can be seen from comparing the values in Tables 5 and 6, and looking at the plots in Figures 5 and 10. For this reason, and because the Lorentzians in the Tables are sorted by their centroid frequency, it is not possible to compare the properties of all the Lorentzians obtained from the fits with the two types of models. As in the case of the constant phase-lags model, in this case there is again a Lorentzian (number 5) that is not significant in the PS but is significant in the CS.

In Figure 11 we plot the observed intrinsic coherence with the derived model for the constant time-lags model. The largest deviations appear again at $0.2$ Hz (see §3.4 for an explanation).

Table 6: Same as Table 5, assuming the constant time-lags model

Component	$\nu_{0}$ (Hz)	FWHM (Hz)	rms${}_{\rm PS}$ (%)	significance${}^{*}$	time lags (ms)	rms${}_{\rm CS}$ (%)	significance${}^{*}$
Lorentzian 1	$0.05\pm 0.01$	$0.063\pm 0.008$	$15.8\pm 0.8$	$10.0$	$242\pm 4$	$16.8\pm 0.8$	$11.3$
Lorentzian 2	$0.165\pm 0.005$	$0.12\pm 0.02$	$14.3\pm 1.1$	$7.1$	$73\pm 3$	$14.9\pm 1.0$	$7.8$
Lorentzian 3	$0.29\pm 0.02$	$0.39\pm 0.03$	$10.7\pm 0.9$	$5.8$	$45\pm 2$	$11.4\pm 0.9$	$6.3$
Lorentzian 4	$0.74\pm 0.06$	$2.4\pm 0.2$	$14.2\pm 0.7$	$9.9$	$-1.9\pm 0.8$	$14.4\pm 0.9$	$8.7$
Lorentzian 5	$1.3\pm 0.3$	$6.8\pm 0.3$	$<2.2^{\dagger}$	$<1$	$2.5\pm 0.5$	$11.9\pm 0.9$	$6.5$
Lorentzian 6	$2.1\pm 0.3$	$5.9\pm 0.2$	$9.2\pm 0.9$	$4.7$	$-17\pm 1$	$5.7\pm 0.5$	$4.8$
Lorentzian 7	$21.0\pm 2.1$	$35.3\pm 3.9$	$1.9\pm 0.1$	$7.0$	$-2.9\pm 0.9$	$3.3\pm 0.4$	$4.4$
$\chi^{2}/\mathrm{dof}$							$417/415$
${}^{*}$ In units of $\sigma$.
${}^{\dagger}$ 95% upper limit.
The rms amplitude of the power and cross spectra of each Lorentzian are integrated from zero to infinity.

Appendix C A QPO in the Imaginary part of the cross spectrum of MAXI J1820+070: the constant time-lags model

In Figure 13 we plot the observed intrinsic coherence with the model derived from the simultaneous fit of the two PS and the CS using the constant time-lags model. As in §3.5, the significant drop in the observed intrinsic coherence at $\sim 2$ Hz, where there is a QPO that is only present in the Imaginary part of the CS, is well described by the model derived from the fit to the two PS and the CS.

Figure 12: Same as Figure 7 but now assuming the constant time-lags model.

Figure 13: Intrinsic coherence function of the same observation of MAXI J1820+070 shown in Figure 12 with the derived model obtained from the fit to the PS and CS assuming the constant time-lags model.