© The Authors 2026

1. Introduction

Black hole X-ray binaries (BHXBs) consist of a stellar-mass black hole and a companion star. Transient BHXBs spend most of their lifetimes in quiescence, occasionally undergoing X-ray outbursts that typically last for weeks to months (see e.g. Remillard & McClintock 2006; Bahramian & Degenaar 2023, for reviews). The soft X-ray emission is generally attributed to the accretion disc (Shakura & Sunyaev 1973), which is characterised by a multi-temperature blackbody spectrum (Mitsuda et al. 1984; Makishima et al. 1986), peaking at around 0.1–2.5 keV (e.g. Zhang et al. 1997; Belloni et al. 1997). A fraction of these soft photons is Compton up-scattered in a plasma of energetic electron (Sunyaev & Titarchuk 1980), the so-called corona, producing a hard X-ray power-law-like spectrum in the hard X-ray band (e.g. Méndez & van der Klis 1997; Done et al. 2007; Gilfanov 2010).

During an outburst, BHXBs typically trace an anticlockwise ‘q’-shaped track in the hardness–intensity diagram (HID; e.g. Homan et al. 2001; Fender et al. 2004; Belloni et al. 2005), transitioning through four spectral-timing states (see Belloni & Motta 2016, for a review). The outburst begins in the low–hard state (LHS), dominated by Comptonised emission, with a weak thermal disc (e.g. Sharma et al. 2018; Shidatsu et al. 2019). As accretion increases, the source brightens and passes through the hard- and soft–intermediate states (HIMS and SIMS, respectively), where the thermal and Comptonised components contribute comparably (e.g. Méndez & van der Klis 1997; Shidatsu et al. 2019). As the outburst continues, the source reaches the high–soft state (HSS), and the X-ray spectrum is dominated by the thermal component with a weak power-law component (e.g. Méndez & van der Klis 1997; Done et al. 2007). During the decay of an outburst, the source returns to quiescence through the SIMS, HIMS, and LHS. Some outbursts are failed-transition outbursts (Alabarta et al. 2021) remain only in the LHS and HIMS.

During outbursts, BHXB X-ray light curves exhibit variability over milliseconds to years (e.g. Motta 2016; Ingram & Motta 2019, and references therein). The Fourier power density spectra (PDS) can be modelled with Lorentzians, representing both broadband noise (BBN) and narrow quasi-periodic oscillations (QPOs; e.g. van der Klis 1989; Belloni et al. 2002). Low-frequency QPOs (LFQPOs) are the most common, with frequencies ranging from a few megahertz to ∼30 Hz (e.g. Belloni et al. 2002; Remillard & McClintock 2006). Low-frequency QPOs are further classified into three types: A, B, and C (Wijnands et al. 1999; Remillard et al. 2002; Casella et al. 2004, 2005), based on the spectral state, quality factor¹ (Q), the PDS shape, and the total fractional root-mean-square (rms) amplitude of the BBN and the QPOs. With centroid frequencies of ∼6 − 8 Hz and low rms (a few percent), type-A QPOs are rarely detected, occasionally appearing in the HSS and displaying relatively broad peaks (Q  ≲  3) (Wijnands et al. 1999; Casella et al. 2004; Belloni & Stella 2014). Type-B QPOs occur only in the SIMS characterised by frequencies ≲10 Hz, low rms amplitudes (≲5%), narrow peaks (Q ≳ 6), and weak (≲10%) red noise (Casella et al. 2004; Belloni & Stella 2014; Belloni & Motta 2016). Type-C QPOs typically appear during the LHS and HIMS, spanning a wide range of centroid frequencies (∼0.01 − 30 Hz) and exhibiting high rms amplitudes (up to ∼20%) and narrow peaks (Q ≳ 8). They are accompanied by strong BBN of up to 40% (Casella et al. 2004; Motta et al. 2011; Belloni & Stella 2014).

The complex Fourier cross spectrum (CS) and the coherence function provide further insights into X-ray timing. The CS argument gives the frequency-dependent phase lags between correlated light curves observed simultaneously in two energy bands (van der Klis et al. 1987; Nowak et al. 1999a; Uttley et al. 2014), while the coherence function quantifies, as a function of frequency, the degree of linear correlation between the light curves measured in two energy bands. If multiple variability components contribute to a certain frequency range in both energy bands, the coherence function may drop below unity, even if each component individually produces perfectly coherent variability (Vaughan & Nowak 1997). This behaviour applies to both narrow QPOs and broad noise components and leads to diverse coherence patterns observed across different sources and observations (e.g. Nowak et al. 1999a; Muno et al. 2001; Cui et al. 2000; Rapisarda et al. 2017a,b; Alabarta et al. 2025).

Quasi-periodic oscillations are typically identified in the PDS, but sometimes the PDS appears smooth, described by two to three broad Lorentzians (e.g. König et al. 2024). Méndez et al. (2024) introduced a novel method for measuring variability component lags in X-ray binaries. Their approach involves simultaneously fitting the PDS as well as the real and imaginary parts of the CS using a linear combination of Lorentzian functions. The model assumes that each individual component is coherent across energy bands, whereas different components are mutually incoherent. This technique enables the detection of weak signals that may be suppressed in the PDS but remain significant in the CS. This method has been successfully applied to several BHXBs, including Cygnus X–1 (König et al. 2024; Fogantini et al. 2025), MAXI J1820+070 (Méndez et al. 2024; Bellavita et al. 2025), and Swift J1727.8–1613 (Brigitte et al. 2025). These studies unveiled a new type of QPO, called ‘imaginary’ QPOs because they are only detected in the imaginary part of the CS. This results in the QPOs having a large phase lag. Therefore, the occurrence of an imaginary QPO is often accompanied by a sharp increase in phase lags and a drop in coherence functions at QPO centroid frequencies. König et al. (2024) reported similar narrow features in the phase-lag and coherence-function frequency spectrum in one of the NICER observations of AT2019wey at a frequency where no apparent narrow component was seen in the PDS. Motivated by these findings, we aim to investigate the presence of imaginary QPOs in AT2019wey and study their properties in detail.

AT2019wey was initially discovered as an optical transient by the Asteroid Terrestrial-impact Last Alert System (ATLAS; Tonry et al. 2018) optical survey in December 2019 (Tonry et al. 2019). X-ray activity was first detected in March 2020 by the Astronomical Roentgen Telescope - X-ray Concentrator (ART-XC; Pavlinsky et al. 2021) and the extended Roentgen Survey with an Imaging Telescope Array (eROSITA; Predehl et al. 2021) on board the Spektrum-Roentgen-Gamma (SRG) satellite (Mereminskiy et al. 2020). Since then, multi-wavelength observations have been conducted. Yao et al. (2021b) estimated that the distance to AT2019wey is between ∼1 and ∼10 kpc and proposed that the compact object was a black hole candidate from its locations in the L_radio − L_X (Bahramian et al. 2018) and L_opt − L_X diagrams (Russell et al. 2006). AT2019wey is a low-inclination (≲30°) system with a low-mass companion star (≲1 M_⊙) and a short orbital period (≲16 h; Yao et al. 2021b). Observations over ∼5 years show that the source is mostly in the LHS and HIMS (Yao et al. 2021a) and may occasionally enter the HSS (Negoro et al. 2023; Yang et al. 2024; Rout et al. 2025b).

The structure of this paper is as follows. In Sect. 2 we describe the NICER observations and data reduction. We provide details of the methods and techniques used in the timing analysis in Appendix A. Our results, including features in the phase lags, coherence function, and timing properties of the imaginary QPO, are presented in Sect. 3. Finally, in Sect. 4 we discuss our findings.

2. Observation and data reduction

Since August 2020, the Neutron Star Interior Composition Explorer (NICER; Gendreau et al. 2016) carried out a long-term monitoring campaign using its onboard X-ray Timing Instrument (XTI). In this work, we used all publicly available NICER archival observations of AT2019wey obtained up to date. The archival dataset comprises more than 600 observations performed between August 4, 2020, and May 7, 2025. We processed the data using HEASOFT (v6.35.2) and NICERDAS (2025_06_11_V014). We used the nicerl2 task with CALDB version xti20240206 to produce the cleaned event files. We only used orbit night data, which is largely unaffected by the light leak, as reported by the NICER team². To minimise potential instrumental effects that could affect the reliability of the timing analysis, we kept only the events with an undershoot rate below 200 counts/s and an overshoot rate lower than 5 counts/s, and a cut-off rigidity (COR_SAX) above 1.5 GeV/c. After filtering, over 400 observations remained for further analysis.

We defined several energy bands of interest for extracting light curves and Fourier products. These include broadbands (0.3 − 2 keV, 2 − 12 keV and 0.3 − 12 keV) as well as narrower sub-bands (0.3 − 0.8 keV, 0.8 − 1.2 keV, 1.2 − 1.6 keV, 1.6 − 1.9 keV, 1.9 − 2.4 keV, 2.4 − 4 keV, 4 − 5 keV and 5 − 12 keV). We extracted light curves in the defined energy bands for each observation using the nicerl3-lc task and estimated background light curves with the Scorpeon model. Prior to running nicerl3-lc, we used nicerl3-spect to generate the response files required for background estimation. We then used the resulting background files to produce background-subtracted light curves and apply background corrections in the computation of the rms normalisation of the power density and cross spectra.

For the X-ray timing analysis we used GHATS³ to compute the fast Fourier transform (FFT) for each energy band of interest of each observation ID, based on the cleaned event files. We computed the FFTs with a segment length of 16.384 s and a time resolution of 0.5 ms. GHATS generates a Leahy-normalised (Leahy et al. 1983) PDS for each 16.384 s segment of each observation ID.

In this work, one of our objectives is to search for narrow features in the phase lags and coherence function. Only a few individual observation IDs of AT2019wey exhibit narrow features with low signal-to-noise ratio (S/N). Previous studies in other sources have shown that the shape of the PDS and CS and their counterparts – characteristic structures in the phase-lag and coherence-function frequency spectra – depend on both spectral hardness and source intensity (Bellavita et al. 2025; Fogantini et al. 2025; Rout et al. 2025a). Based on this, we assumed that observations within a narrow hardness ratio (HR) interval share similar timing properties. To enhance S/N and enable a more detailed investigation of these features, we selected 16.384 s segments in a range of HR and averaged the PDS and real and imaginary parts of CS of these segments. Afterwards, we used the averaged PDS and real and imaginary parts of CS in two energy bands to compute the phase-lag frequency spectrum and coherence-function frequency spectrum for each selection. We corrected each PDS for Poisson noise level by subtracting the average power in the 200–800 Hz range where no other timing features are present. Following Belloni & Hasinger (1990), we normalised the PDS and CS to units of fractional rms squared per hertz considering the background contribution. We checked that the background light curve is consistent with being constant and therefore adopted the average background count rate of the entire selection in the corresponding energy band for the correction. When computing the CS, phase-lag, and coherence-function frequency spectra, we used the 0.3 − 2.0 keV band as the reference band for comparisons between the 0.3 − 2 keV and 2 − 12 keV ranges. For the energy-dependent analysis, we took the total 0.3 − 12 keV band as the reference, while the narrower sub-bands were treated as subject bands. To correct for partial correlations arising from photons that are both in the subject and reference bands, we subtracted the average real part of the CS measured over a frequency range free of source contributions (Ingram 2019; Belloni et al. 2024; Méndez et al. 2024; Bellavita et al. 2025; Fogantini et al. 2025). We rotated⁴ all the CS by 45° to have approximately equal real and imaginary parts; as shown by Méndez et al. (2024) this is equivalent to a change of coordinates and therefore does not affect the final results. We applied a logarithmic rebinning to all Fourier products such that the size of each frequency bin increases by 10^1/100 with respect to the previous one to increase the S/N. For all Fourier products, the segment length and time resolution determine the accessible frequency range, yielding a minimum Fourier frequency of approximately 1/16.384 Hz and a Nyquist frequency of 1000 Hz. Considering the frequency range of the variability components of interest, we restricted our analysis to Fourier products below 20 Hz.

3. Results

3.1. Grouping HID

Figure 1 presents the HID of AT2019wey, where the HR is defined as the ratio of the count rate in the 2 − 12 keV and 0.3 − 2 keV energy ranges and the intensity is defined as the count rate in the 0.3 − 12 keV energy range. Each data point corresponds to a single 16.384 s segment. As mentioned in Sect. 2, only a few observations display possible characteristic narrow timing variability features in the phase-lag and coherence-function frequency spectra extracted from individual observation IDs with low S/N. To ascertain the existence of these features, we combined 16.384 s segments with similar HR to enhance the S/N. In the inset of Fig. 1, we show in black the points along the main track of the HID for the segments that we retain, while in grey, we show those that deviate from the main track. These latter points are excluded from the rest of the analysis. We first divided the HID into narrow intervals according to HR width. These intervals allow us to examine the gradual evolution of the PDS and CS shapes, as well as the characteristic features in the phase-lag and coherence-function frequency spectra. Adjacent intervals with similar spectral timing properties were then merged into ten zones, with widths in HR of 0.015 for zones 1–5, 0.02 for zone 6, 0.025 for zone 7, 0.065 for zone 8, 0.125 for zone 9, and 0.2 for zone 10. These zones are colour-coded in Fig. 1. Three groups of grey data points represented with squares (MJD 59075), stars (MJD 59083), and triangles (MJD 59112) in the main plot represent the three groups of data from Yao et al. (2021a) during the initial outburst, before MJD 59122. They find that the source begins to enter the HIMS at the position marked by the grey stars. Therefore, we divided the data around the sampling point at HR = 0.21, as indicated by the vertical dashed black line, to separate the states. Data with HR < 0.21 are classified as belonging to the HIMS, while those with HR > 0.21 correspond to the LHS. To verify the robustness of the hardness-resolved analysis, we divided the data in each interval shown in Fig. 1 into high- and low-count-rate subsets and find no significant differences in their timing properties. We further checked the stability of the PDS and cross-spectral products on shorter timescales by examining them over a number of consecutive NICER orbits within individual observation IDs. We find no significant changes over these shorter timescales.

Fig. 1. HID constructed from the NICER observations of AT2019wey between August 9, 2020, and May 7, 2025. Each data point represents a 16.384 s segment. Main panel: Ten colour-coded regions corresponding to the data groups defined by hardness intervals and used in this work. The grey squares, grey stars, and grey triangles mark the three epochs analysed by Yao et al. (2021a), identified as the LHS (grey squares) and HIMS (grey stars and triangles), respectively. The vertical dashed line at HR = 0.21 roughly separates the LHS from the HIMS. Inset panel: Data used in our analysis (black points; same as that in the main panel), along with data excluded from our analysis (grey points).

3.2. The imaginary QPO

We simultaneously fitted the two PDS (0.3 − 2.0 keV and 2.0 − 12 keV) together with the real and imaginary parts of the CS for each zone, using the method described in Sect. A.2. We detect a QPO component in zones 2–9. The best-fit models for the remaining zones are presented in Appendix C.1. We use zone 5 as a representative example to illustrate our results. As shown in the top panels of Fig. 2, the two PDS and the real and imaginary parts of the CS are well described by six Lorentzian components, yielding a reduced χ² of 0.969 (410.99/424 d.o.f). Both PDSs in these energy bands are smooth, showing no narrow peaks that typically indicate a QPO. However, a narrow Lorentzian component at ∼2.88 Hz is clearly detected in the CS, where it exhibits a much stronger imaginary than real part in the rotated CS, with a quality factor of Q ∼ 3. In the non-rotated CS, the centroid frequency of this component aligns with the frequency of the peak, as displayed in the middle right panel of Fig. 2. Across zones 2–9, we find that this component is not always significantly detected in the PDS of either the 0.3 − 2 keV or 2 − 12 keV energy bands, but it is consistently significant and strong in the imaginary part of the CS. On this basis and given that it shares several properties (i.e. the sudden changes in phase lags and coherence function; see the next section) with the imaginary QPOs reported in other sources (e.g. Méndez et al. 2024; Bellavita et al. 2025; Fogantini et al. 2025; Brigitte et al. 2025), we hereafter call this component the imaginary QPO. The bottom panels of Fig. 2 shows the comparison between the derived model (solid black curve) and the observed phase lags and coherence function, indicating agreement between them. At ∼2.88 Hz, the phase lags exhibit a sharp increase (‘cliff’), accompanied by a simultaneous drop (‘dip’) in the coherence function. These two behaviours are coincident with the frequency of the imaginary QPO. We note that the solid black line in the bottom panels is the prediction of the phase lags and coherence function derived from the fits to the two PDS and the real and imaginary parts of the CS, and that we did not fit the phase lags and coherence function. We mark the frequency of the imaginary QPO with a vertical dotted black line in each panel. As illustrated by the dashed black curve in the bottom panel of Fig. 2, we also show the model without the imaginary QPO and without refitting. This comparison demonstrates that the imaginary QPO is responsible for producing the maximum in the phase lags and the minimum in the coherence function. When we refit the spectrum without the imaginary QPO, one of the other components shift and broaden to occupy its frequency range, indicating that the imaginary QPO is a necessary component. However, the resulting fit is significantly worse. Comparing both fits, an F-test yields an F-statistic of 16.85 and a probability of 1.79 × 10⁻¹⁷, indicating that the QPO is highly significant.

Fig. 2. Simultaneous fitting of the PDS and CS of AT2019wey with a multi-Lorentzian model assuming constant phase lags, using the combined NICER data in zone 5 as an example. Top panels: Best fits (solid black lines) and residuals for the PDS (left panel) in the 0.3 − 2 keV (blue) and 2 − 12 keV (red) bands, and to the real (blue) and imaginary (red) parts of the CS rotated by 45° (right panel). Middle panels: Real (left) and imaginary (right) parts of the non-rotated CS, showing a linear scale on the y-axis. Each model consists of six Lorentzian components; the one corresponding to the imaginary QPO is highlighted with thicker lines. Bottom panels: Phase-lag (left) and coherence-function (right) frequency spectra, along with their residuals with respect to the model derived (solid black curves) from the simultaneous fitting of PDS and non-rotated CS. The dashed black curves show the derived model without the imaginary QPO and without refitting. We rebinned the plots by a factor of 5 to improve the visibility of the dip in the coherence function.

Table 1 lists the parameters (with 1σ errors) of the imaginary QPO in each zone. The frequency of the imaginary QPO decreases from 4.95 Hz to 1.11 Hz as the source spectrum hardens, with the HR increasing from HR ≈ 0.1 (zone 2) to HR ≈ 0.4 (zone 9). The quality factor of the imaginary QPO is between ∼1.5 − 3.0. The phase lags are hard for all zones, increasing from ∼0.3 rad to ∼0.7 rad as the source hardens. The $\sqrt{{\mathrm{A}}_{\mathrm{i}}{\mathrm{B}}_{\mathrm{i}}}$ of the imaginary QPO in each zone is consistent within errors with C_i, in accordance with the assumptions described in Sect. A.1, suggesting that the corresponding Lorentzian components are coherent in the two energy bands.

Motivated by the QPO classification based on the relations between total fractional rms and QPO frequency (Casella et al. 2005; Motta et al. 2011, 2012), Fig. 3 displays the correlation between the frequency of the imaginary QPO and the broadband (0.1 − 20 Hz) fractional rms in the soft (0.3 − 2 keV), hard (2 − 12 keV) and full (0.3 − 12 keV) energy bands. The broadband fractional rms amplitude in the soft band increases from 10% to 31% as the QPO frequency decreases from ∼5 Hz to ∼1 Hz and the source moves to the hard state. In the hard band, the broadband fractional rms amplitude increases from 23% to 31% as the QPO frequency decreases to ∼2 Hz. Below 2 Hz, corresponding to zones 8 and 9, the rms amplitude begins to decrease as the frequency decreases. In the full band, the evolution of the broadband rms amplitude with QPO frequency is the same as that in the soft band.

Fig. 3. Frequency of the imaginary QPO vs the 0.1 − 20 Hz broadband fractional rms in the 0.3 − 2 keV (circles), 2 − 12 keV (triangles), and 0.3 − 12 keV (squares) energy bands derived from our model. Each point represents the observation in a zone. The colour bar indicates the evolution of the HR across all zones. The two grey arrows separate the observations in the LHS (left) and HIMS (right).

3.3. The increase in phase lags and drop in the coherence function

Figure 4 presents the phase-lag and coherence-function frequency spectra between 0.3 − 2 keV and 2 − 12 keV for each zone, from the softest zone to the hardest zone. The results for different zones are shown in using the same colours as in Fig. 1. The vertical dashed black line indicates the position of the imaginary QPO, obtained from the simultaneous multi-Lorentzian fit to the PDS and CS (see Sect. 3.2 for the imaginary QPO). For zone 1, there are no apparent characteristic structures in the phase lags and coherence function. In zone 2, the phase-lag frequency spectrum exhibits a sudden increase in the phase lag at ∼2 Hz from 0 rad to a maximum value of ∼0.3 rad, followed by a gradual return to 0 rad at ∼10 Hz. We refer to this abrupt increase in the phase lag as the ‘cliff’. However, no distinct feature is visually identified in the coherence function at this frequency. From zones 3 to 8, in addition to the presence of the phase-lag cliff, the coherence exhibits a sudden drop at around the same frequency of the cliff, followed by a recovery to unity, forming a distinct ‘dip’ structure. In zone 9, a small phase lag cliff appears at ∼1.1 Hz, above which the phase lags form a plateau, while the dip in the coherence function disappears again. For zone 10, no obvious phase-lag cliff is observed. Below ∼10 Hz, the phase lags are hard and gradually increases from low frequencies and flattens above ∼1 Hz, whereas the coherence displays no evident structure. In all the zones in which the dip is present, the frequency at which the coherence is minimum coincides with that of the imaginary QPO. We note that, in most cases, the peak of the cliff in the phase lags also corresponds to the imaginary QPO frequency. However, in zones 2, 3, and 8, the peak appears at a slightly higher frequency, similar to what was observed in Cygnus X–1 (Fogantini et al. 2025). In our case, this is due to an additional QPO at an adjacent higher frequency, whose larger phase lag shifts the maximum of the phase-lag frequency spectrum to higher frequency. In zone 2, the phase lag increases from 0.33 rad at the imaginary QPO frequency to 0.49 rad at the frequency of the adjacent component. In zone 3, the lags rise from 0.38 rad to 0.81 rad, and in zone 8, from 0.47 rad to 0.55 rad. In the zones where either the phase-lag cliff or the coherence-function dip is present, we find that both features gradually shift towards lower frequencies as the source hardens.

Fig. 4. Phase lags and coherence function for each region of the HID of AT2019wey for the 0.3–2 keV and 2–12 keV bands. The vertical dashed black line marks the frequency of the imaginary QPO (see text) where the phase-lag cliff and coherence-function dip occur. The colour coding used in the frequency spectra for the different zones matches that of Fig. 1.

3.4. Energy dependence of rms amplitude, normalised covariance, and phase lags of the imaginary QPO

To study the energy-dependent properties of the imaginary QPO, we performed separate simultaneous fits, each involving the PDS in the 0.3 − 12 keV band, the PDS in one of the narrow bands defined in Sect. A.1, and the CS between them, following the procedure described in Sect. A.2. For the zones exhibiting the ‘cliff’ or ‘dip’ features, we constructed the fractional rms energy spectra (blue points) of the imaginary QPO measured from the PDS, as presented in Fig. 5. The error bars indicate 1σ uncertainties, while the upper limits of the measurements at higher energies are given at the 95% confidence level.

Fig. 5. rms (red), normalised covariance (blue), and phase-lag (orange) spectra of the imaginary QPOs in AT2019wey for different zones in the HID. All quantities are derived from the fits using our multi-Lorentzian model. The error bars represent 1σ uncertainties. The arrows in the rms spectra represent the 95% upper limit.

Across zones 2–9, the rms spectra generally show two types of behaviours. Zone 2 displays a rising-flat profile, with the rms increasing at low energies (≲2.2 keV) and remaining roughly constant at higher energies. In all other zones, the rms spectra exhibit a possible rising-falling profile: the rms increases towards ∼1.5 − 2 keV and then decreases at higher energies, although the detailed shape above ≳3 keV is often unconstrained. This profile is particularly pronounced in zone 6. The decreasing trend at higher energies can be identified at 2.2 keV and 3.2 keV. At higher energies the large uncertainties and upper limits make it difficult to assess the shape of the rms spectrum. Regarding the evolution of the maximum fractional rms amplitude, during the HIMS (zones 2–7), the maximum fractional rms amplitude decreases from 7.2% to 4.3% as the source hardens. In the LHS (zones 8–9), the maximum fractional rms amplitude initially rises back to 7.6% in zone 8, but then decreases to 5% in zone 9 as the hardness increases. Correspondingly, we also find that the energy at which the maximum fractional rms amplitude occurs shifts to lower energies with increasing hardness during HIMS: it is approximately 2.2 keV in zones 2–4 and about 1.4 keV in zones 5–7. Similarly, during the LHS, the energy of the maximum fractional rms decreases from 2.2 keV (zone 8) to 1.4 keV (zone 9).

To verify the shape of the rms spectra over the entire energy range, we also plot the normalised covariance spectrum (see Appendix A.3 for definition), shown as the red points in Fig. 5, to cross-check the shape of the rms spectrum at the higher energies. The normalised covariance is equal to the rms amplitude of the imaginary QPO when the signals in the two energy bands are perfectly coherent. The normalised covariance spectrum shows good overall agreement with the rms spectrum, particularly below ∼2 keV. Considering the 95% upper limits, this agreement extends above ∼2 keV as well. In zone 2, the normalised covariance follows the same shape as the rms spectrum, increasing from ∼1% at ∼0.55 keV to ∼8% at ∼2 keV and above. For the other zones, above ∼2 keV, the normalised covariance decreases with increasing energy. We observe that the energy of maximum of the normalised covariance spectrum moves towards lower energies with increasing HR, from zone 2 to zone 9. We conclude that the normalised covariance of the imaginary QPO exhibits a rising-flat trend with energy in zone 2, whereas in the remaining zones it shows a rising-falling behaviour.

In Fig. 5, the orange points show the phase-lag spectra. We fitted a quadratic polynomial to the data around the minimum of the phase-lag spectra of zones 2–7, which show the possible U-shape. We find that the energy at which the phase lags are minimum, E_min, decreases as the source hardens and the frequency of the imaginary QPO decreases, with E_min = 1.46 ± 0.17 keV (zone 2), 1.24 ± 0.09 (zone 3), 1.26 ± 0.05 (zone 4), 1.13 ± 0.06 (zone 5), 1.03 ± 0.04 (zone 6), and 0.98 ± 0.14 (zone 7). We used the E_min − ν_imQPO relation with a linear function, E_min = aν_imQPO + b and find a = 0.18 ± 0.05 and b = 0.60 ± 0.13. This shows that the trend is 3.6σ significant. For the remaining zones (8 and 9), the typical U-shaped feature in the phase-lag spectrum disappears as the source hardens further. In zone 8, the phase lag increases monotonically with energy below ∼4.5 keV and then decreases slightly or remains constant above this energy. In zone 9, the phase lags show a similar trend to that in zone 8, but the phase lags remain nearly constant at higher energies within 1σ errors. Based on the state classification in Sect. 3.1, we find that, within the current energy range, the U-shaped phase-lag spectra occur only during the HIMS and disappear in the LHS.

4. Discussion

We performed a simultaneous analysis of the PDS, CS, phase lags, and coherence function of AT2019wey using all available NICER observations, following the technique proposed by Méndez et al. (2024). We detect a QPO that is not always significant in PDS but is significant in the imaginary part of the CS. Following previous works, we call it an imaginary QPO. The phase-lag frequency spectrum exhibits a sudden increase and the coherence-function frequency spectrum shows a drop in the QPO frequency. We refer to these two narrow features as the ‘cliff’ and the ‘dip’, respectively. The fits to the PDS and CS correctly reproduce the cliff of the phase-lag frequency spectrum and the dip of the coherence-function frequency spectrum, lending further support to the assumption that the PDS and CS can be fitted with a linear combination of Lorentzian functions that are coherent in two energy bands but incoherent with one another. The imaginary QPO has a large hard phase lag (> 0.3 rad). We find that the centroid frequency of the imaginary QPO decreases from ∼5 Hz to ∼1 Hz as the source spectrum hardens. The fractional rms amplitude and normalised covariance of the imaginary QPO increase at energies below ∼2 keV but, interestingly, remain more or less constant or show a possible decreasing trend at higher energies. Additionally, we observe a typical U-shaped phase-lag spectrum, and a continuous change in the shape of the phase-lag spectrum as the source transitions from the HIMS to the LHS.

4.1. The implication of the phase-lag cliff

At the frequency of the imaginary QPO, the phase lags exhibit a sharp increase, forming a pronounced phase-lag cliff. This behaviour can be naturally attributed to the fact that the imaginary QPO has a larger phase lag than the other components in the same frequency range. Several models have been proposed to explain the lags of QPOs. In the jet precession model (Ma et al. 2021), if lower-energy photons are emitted from higher altitudes of a precessing jet, soft lag would be produced because the jet base comes into view first. The soft QPO lags predicted by this model are inconsistent with our results. In the time-dependent Comptonisation model, vKompth (Karpouzas et al. 2021; García et al. 2021; Bellavita et al. 2022), photons from the disc that are inverse-Compton scattered in the corona and escape later produce hard lags, while reprocessing of returning Comptonised photons (i.e. feedback) in the disc introduces soft lags. This model provides a plausible explanation for the hard phase lags observed in the imaginary QPO.

Interestingly, compared to other sources where the imaginary QPO was detected, the significance of the phase-lag cliff and coherence-function dip in AT2019wey exhibit a different dependence on the chosen reference energy bands. A clear phase-lag cliff remains visible in AT2019wey when the 2 − 4 keV band is used as the reference band (Fig. D.1), in contrast with other sources (e.g. König et al. 2024; Fogantini et al. 2025; Bellavita et al. 2025). Assuming the phase lag can still be interpreted within the vKompth model, this finding may indicate that the accretion disc of AT2019wey contributes non-negligibly to the emission above 2 keV. However, this interpretation may not be unique. As mentioned in Appendix D, the phase-lag cliff above 2 keV is primarily associated with a broad variability component at frequencies higher than that of the imaginary QPO. This broad component has a larger imaginary part than real part, which causes the hard lag. The high-frequency broad variability component in the PDS could be directly associated with the Comptonising medium (König et al. 2024). Therefore, the phase-lag cliff above 2 keV could originate solely from the Comptonising region, in which higher-energy photons undergo more scatterings and are emitted later, naturally producing the hard lag observed between photons above 2 keV.

Under the assumption that the time lag associated with the QPO reflects the light-travel time across the Comptonising region, one can describe the evolution of the coronal size. For a 10 M_⊙ black hole, the size of the corona of AT2019wey would increase from ∼3.2 × 10³ km to 1.5 × 10⁴ km (215 − 1023 R_g), as the source hardens, corresponding to an increase in the phase lag of the imaginary QPO from ∼0.3 rad to ∼0.7 rad during the HIMS. In this case, the inferred values should be regarded as upper limits for a Comptonising cloud with small optical depth (τ ≳ 0). The actual size would be much smaller for τ ≫ 1. From the fitting of the energy spectrum, Yao et al. (2021b) found that during the HIMS the inner disc radius remained almost constant at ∼100 km (1000 km) assuming a distance of ∼10 kpc (1 kpc). The coronal sizes inferred from the time lags are systematically larger than the inner disc radii derived from the spectral fitting. Furthermore, while the spectral analysis suggests that the inner disc radius remains nearly constant during the HIMS, the time lags indicate that the characteristic coronal size increases by approximately an order-of-magnitude.

4.2. The implication of the coherence-function dip

The coherence can be lower than unity when multiple, physically independent emission regions that are incoherent with one another contribute to the variability in both energy bands, even if each region individually produces perfectly coherent signals (Vaughan & Nowak 1997). In the PDS of AT2019wey (see Fig. 2), multiple Lorentzian components overlap within the narrow frequency range where the coherence dip is observed. If these Lorentzian components correspond to physically independent variability processes, the observed coherence dip is naturally accounted for within the framework proposed by Vaughan & Nowak (1997).

Hua et al. (1997) examined the physical origins of the reduced coherence function in hard X-ray bands above 2 keV. They showed that evolution of the macroscopic parameters of the Comptonising region and/or in the energy of the seed photons during an observation can lead to a loss of coherence. In their scenario, multiple variability components within a given energy band were not considered. Instead, the light curves in two different energy bands were assumed to originate from two distinct coronal configurations, effectively representing two independent time series. As a result, the coherence between the two bands naturally drops below unity. However, in contrast to our findings, such a scenario would lead to a loss of coherence over a broad range of frequencies (see their Fig. 2).

Similar coherence-function dips have been reported during the decay phases of MAXI J1820+070 (Bellavita et al. 2025), MAXI J1348–630 (Alabarta et al. 2025), Swift J1727.8–1613 (Brigitte et al. 2025), and Cygnus X–1 (Fogantini et al. 2025), all of which appear in the lower branch of the q-shaped HID and hence at relatively low accretion rates (König et al. 2024). These results suggest that the sudden coherence dips preferentially occur at low accretion rates, possibly associated with a coronal geometry similar to that inferred for MAXI J1348–630. For AT2019wey, fitting the brightest NICER spectrum with vphabs*(diskbb + nthComp + relxillCp) yields a 1.5 − 10 keV luminosity of (1.2 − 2.5)×10⁻² L_Edd for a distance of 10 kpc, or (1.2 − 2.5)×10⁻⁴ L_Edd for 1 kpc, assuming a black hole mass of 5 − 10 M_⊙ (Yao et al. 2021b). These values would place AT2019wey on the lower branch of the HID (König et al. 2024).

4.3. Energy dependence of the normalised covariance and phase lag of the imaginary QPO

As described in Sect. 3.4, the fractional rms spectra do not provide a good constraint at higher energies, but the covariance spectra give a statistical significant trend with energy. As shown by Wilkinson & Uttley (2009), the normalised covariance is equivalent to the fractional rms amplitude (See Appendix A.3). Therefore, we discuss the energy dependence of the variability using the results obtained from the covariance spectra. The covariance spectra of the imaginary QPO in AT2019wey exhibit two types of profiles: a rising-flat profile observed only in zone 2, and a rising-falling profile seen in the other zones. The rising trend of the normalised covariance of the imaginary QPO in AT2019wey could indicate that the Comptonised component dominates the variability (Zhang et al. 2020b). A similar rising trend in the fractional rms spectra of the imaginary QPO has been observed in Cygnus X–1 (Fogantini et al. 2025), MAXI J1820+070 (Bellavita et al. 2025), and Swift J1727.8–1613 (Brigitte et al. 2025).

Despite these similarities, AT2019wey exhibits notable differences compared to these sources. In Cygnus X–1, MAXI J1820+070, and Swift J1727.8–1613, the rms spectra of the imaginary QPOs generally show a rising-flat profile, whereas in AT2019wey the energy dependence of the normalised covariance predominantly follow a rising-falling trend. In addition, the energy at which the normalised covariance peaks is systematically lower in AT2019wey than in MAXI J1820+070 (Bellavita et al. 2025) and Swift J1727.8–1613 (Brigitte et al. 2025). Specifically, the rms maximum occurs at ∼1.4 − 2.2 keV, comparable to that observed in Cygnus X–1 (Fogantini et al. 2025), while in MAXI J1820+070 and Swift J1727.8–1613 the fractional rms amplitude continues to increase up to ∼10 keV.

The declining normalised covariance may be associated with the coronal electron temperature. Bellavita et al. (2022, see their Fig. 3) demonstrates that in their model a lower coronal electron temperature leads to signs of suppression of variability at higher energies and causes the energy of the maximum of fractional rms amplitude to shift towards lower values. This theoretical prediction appears to be consistent with our observational results, suggesting that AT2019wey hosts a relatively cold corona. Supporting this interpretation, Sahu et al. (2026) report the presence of two Comptonising regions in AT2019wey: a hot (kT_e ∼ 26 − 370 keV) and a cold corona (kT_e ∼ 1.2 − 2.7 keV). We therefore suggest that the imaginary QPO in AT2019wey may originate from variability within the cold corona, which might explain the normalised covariance peaking at low energies and declining towards higher energies.

Gierliński & Zdziarski (2005) investigated the pattern of the rms spectra of the BBN in the BHXBs, XTE J1650–500 and XTE J1650–500, over the 2 − 50 keV. They found that in the LHS and intermediate state the rms amplitude either decreases or remains approximately constant with increasing energy, while in the HSS it increases with energy. They suggested that a flat rms spectrum arises from variations in the normalisation of the Comptonised component in the energy spectrum, with no change in the spectral shape. They interpreted the decreasing rms amplitude as resulting from the variation in the seed-photon input, together with some required variation in the power released in the Comptonised component. We note, however, that their analysis concerned the rms spectrum of the BBN component, integrated in the broad frequency band, whereas in this work we focus on the rms spectrum of an individual variability component, the imaginary QPO. It therefore remains unclear whether the idea proposed by Gierliński & Zdziarski (2005) for the BBN can also be applied to QPOs. It remains to be seen how the decomposition of the BBN into multiple Lorentzian components, each characterised by its own rms spectrum, is reflected in the rms spectrum of BBN.

We note that the decrease in the variability amplitude at higher energies could be due to the partial loss of the imaginary QPO signal in the low S/N regime. As explained by Troyer et al. (2018) for the kilohertz QPO in the neutron star 4U 0614+09, at high energies the QPO signal can become dominated by noise, leading to a reduction in the variability amplitude.

The phase-lag spectra of AT2019wey exhibit a U-shaped pattern during the HIMS, which disappears in the LHS. A similar behaviour is also observed for the imaginary QPO in MAXI J1820+070 (see Fig. B.1 of Bellavita et al. 2025), where the U-shape in the phase-lag spectrum also starts to disappear as the source hardens, even if in that case the source did not reach the LHS. In both cases, the energy at the minimum of the U-shaped phase-lag spectra shifts towards lower energies as the source spectrum hardens. A similar behaviour has been reported in the imaginary QPO in Cygnus X–1 (Fogantini et al. 2025), where the energy of the minimum of the phase-lag spectrum was found to scale with the accretion disc temperature. They suggest that this phenomenon can be explained by the time-dependent Comptonisation model, vKompth (Karpouzas et al. 2021; García et al. 2021; Bellavita et al. 2022). In this model, the feedback produces the characteristic U-shaped phase-lag spectrum, with the energy of the phase-lag minimum being positively correlated with the temperature of the soft-photon source (i.e. accretion disc). Qualitatively, in AT2019wey, a similar trend could be seen in Fig. 5 if, as expected in the transition from the HIMS to the LHS, the disc temperature decreases as the source hardens from zone 2 to zone 9. The disappearance of the U shape in the LHS might indicate that the minimum of the phase-lag spectrum has moved to even lower energies below the NICER band. Further investigation of the application of the vKompth model to the data is beyond the scope of this paper and will be addressed in future work.

4.4. Comparison of imaginary QPOs with type-B/C QPOs

The rising-flat rms spectrum of the imaginary QPO in AT2019wey shown in zone 2 (see Fig. 5) is reminiscent of that of either the type-C or type-B QPOs in other BHXBs, such as MAXI J1348–630 (e.g. Alabarta et al. 2022; Liu et al. 2022), GRS 1915+105 (e.g. Zhang et al. 2020a; Karpouzas et al. 2021), MAXI J1535–571 (e.g. Zhang et al. 2022), MAXI J1820+070 (e.g. Ma et al. 2023), and GX 339–4 (e.g. Peirano et al. 2023). On the other hand, the rising-falling rms spectrum of the imaginary QPO in AT2019wey has also been observed for the type-C QPO during the first reflare of MAXI J1348–630 (Alabarta et al. 2022). Moreover, the U-shape phase-lag spectrum of the imaginary QPO in AT2019wey during the HIMS (zones 2–7 of Fig. 5) has often been seen in both type-B and type-C QPOs (see references above). These comparisons suggest that the energy-dependent rms amplitude and phase-lag characteristics of the imaginary QPO is similar to that of both type-B and type-C QPOs.

While these similarities point to a close phenomenological connection between the imaginary QPO and type-B/C QPOs, additional criteria are required to identify the specific QPO class. Given that type-B QPOs are generally observed during the SIMS, whereas the imaginary QPO in AT2019wey appears in the HIMS and LHS, it is more likely that the imaginary QPO is associated with type-C rather than type-B QPOs. Moreover, during the decay of the outburst in MAXI J1348–630, Alabarta et al. (2025) reported a narrow Lorentzian component significant both in the PDS and CS, which they identified as a type-C QPO. In that source, the centroid frequency of the type-C QPO coincides with a cliff-like structure in the phase-lag frequency spectrum and a dip in the coherence function, providing strong evidence that the imaginary QPO in that source, and by extension in AT2019wey, is a type-C QPO.

Further support for a type-C association comes from the frequency-dependent timing behaviour. An anti-correlation between the frequency of the imaginary QPO and the broadband fractional rms amplitude is observed in AT2019wey (Fig. 3), with the broadband rms reaching relatively high values (30%). A similar anti-correlation has been observed for the imaginary QPO in MAXI J1820+070 (see Fig. 6 of Bellavita et al. 2025). This behaviour closely resembles that of the type-C QPOs observed in GX 339-4 (Motta et al. 2011), GRO J1655–40 (Motta et al. 2012), and MAXI J1631–479 (Rout et al. 2021).

The frequency of the imaginary QPOs in AT2019wey decreases with spectral hardening. This trend is reminiscent of the type-C QPOs observed during the hard-to-soft transition in BHXBs, whose frequency increases as the energy spectrum softens (e.g. Sobczak et al. 2000; Vignarca et al. 2003; Shaposhnikov & Titarchuk 2009; McClintock et al. 2009; Zhang et al. 2022; Rawat et al. 2023). Yao et al. (2021a) report the detection of the type-C QPO in AT2019wey during MJD 59112.24–59112.98 and 59083.85–59083.94 (marked as grey triangles and stars in Fig. 1) at 6.58 ± 0.21 Hz and 2.06 ± 0.03 Hz, respectively, whereas the corresponding imaginary QPOs in the same hardness-ratio range are detected at 4.95 ± 0.21 Hz and 2.21 ± 0.04 Hz. The discrepancy between in the second interval of Yao et al. (2021a) and our HR selection may be explained if the QPO frequency depends both upon HR and intensity, since in that time interval the source was brighter than in our selection. In the first time interval of Yao et al. (2021a), their QPO frequency are consistent with ours within errors, and the average intensity of the source in their selection and ours are similar. We also examined the data in the time intervals of MJD 59112.24–59112.98 and 59083.85–59083.94 and simultaneously fitted the PDS and CS. We detect an imaginary QPO and find that it has the same frequency as that of the type-C QPO detected by Yao et al. (2021a). Therefore, we suggest that the imaginary QPO in AT2019wey is the type-C QPO.

4.5. Possible origin of the imaginary QPO

The imaginary QPO in AT2019wey can be interpreted as a distinct signal arising from independent components within the accretion system, as suggested by the linear combination model of multiple Lorentzian functions (Méndez et al. 2024). Rout et al. (2025a) discovered a similar coherence dip in Cygnus X–1 at ∼0.05 Hz, corresponding to the frequency of a hidden QPO in the hard energy band, in Cygnus X–1. In their case the dip in the coherence function is present when they use two hard energy bands above 3 keV. They propose that the dip arises from variability produced by Comptonisation in the jet and other Comptonising components. In the case of AT2019wey, evidence for jet activity has been reported. Yao et al. (2021b) observe radio brightening as the source transitioned from the LHS to the HIMS, and Yadlapalli et al. (2021) detect a resolved radio source during the HIMS, interpreted as a steady compact jet. These observations raise the possibility that the imaginary QPO in this source may also be linked to jet activity. However, in AT2019wey we do not detect a significant coherence dip and a QPO component when we use only data above 2 keV (see Appendix D). In addition, the hidden QPO and its associated coherence dip in Cygnus X–1 reported by Rout et al. (2025a) occur at much lower frequencies than those of the imaginary QPO in AT2019wey, which may indicate a fundamentally different physical origin from the imaginary QPOs observed at frequencies of > 1 Hz. Therefore, we cannot conclusively establish a connection between the imaginary QPO in AT2019wey and the jet.

The imaginary QPO could be related to a beat between two broad Lorentzian components (König et al. 2024). However, in AT2019wey an additional additive narrow Lorentzian component is required to account for the residuals in the imaginary part of the CS, suggesting that the imaginary QPO is not the result of non-linear coupling between the other broad Lorentzian components present in the PDS.

The imaginary QPO could also be related to the interplay between the corona and the accretion disc. In MAXI J1820+070 (Bellavita et al. 2025) and Cygnus X–1 (König et al. 2024; Fogantini et al. 2025), the imaginary QPO is not detected when reference energy bands above 2 keV are used, and both the associated phase-lag cliff and coherence dip disappear. Similarly, in AT2019wey, as mentioned in Appendix D, the imaginary QPO is not significantly detected above 2 keV. These findings collectively indicate that the soft energy band plays a crucial role in revealing the imaginary QPO, strongly implying that the imaginary QPO may be produced by disc–corona radiative coupling processes (Karpouzas et al. 2021; García et al. 2021; Bellavita et al. 2022; Mastichiadis et al. 2022). Furthermore, the evolutionary behaviour of the phase-lag spectrum from the HIMS to the LHS, which is consistent with the predictions of the time-dependent Comptonisation model vKompth (Bellavita et al. 2022), as discussed in Sect. 4.3, provides additional evidence that disc–corona radiative coupling processes are a plausible mechanism for the production of imaginary QPOs.

4.6. Investigation of the inclination dependence of the phase lag of the imaginary QPO

Based on a sample of 15 black hole low-mass X-ray binaries, van den Eijnden et al. (2017) find that both in low- and high-inclination systems, the type-C QPO exhibit slightly hard lags at low QPO frequencies (≲1 Hz). Conversely, at higher frequencies, high-inclination sources transition to soft lags, whereas low-inclination sources show increasingly harder lags. van den Eijnden et al. (2017) suggest that this behaviour is consistent with the idea that the type-C QPOs reflect the Lense-Thirring precession of the inner accretion flow (Ingram et al. 2009)

We investigated whether the phase lags of the imaginary QPO exhibit an inclination dependence similar to that of the type-C QPOs. AT2019wey appears to have a low inclination, as indicated by its single-peaked hydrogen lines in the optical spectrum (Yao et al. 2021b), consistent with the low inclination (i ≲ 30°) inferred from X-ray reflection modelling (Yao et al. 2021a). For Cygnus X–1, the optical orbital modulation indicates a low orbital inclination of i ∼ 27° (Orosz et al. 2011), whereas XRISM/Resolve spectral analysis suggests a higher inner disc inclination of i ∼ 63° (Draghis et al. 2025). This difference may be explained if the accretion disc is warped (Tomsick et al. 2014) and the X-ray spectrum reflects a highly inclined disc that is not aligned with the orbit of the secondary in this source. MAXI J1820+070 has a high orbital inclination of i = 66° −81° (Torres et al. 2020) and jet inclination of i = 64° ±5° (Wood et al. 2021). Relativistic reflection spectroscopy also indicates a high disc inclination of $i=64{°}_{-4}^{+3}$ (Draghis et al. 2023a). The jet inclination of Swift J1727.8–1613 is constrained to be i < 74° (Wood et al. 2024). A moderate inclination in this source is indicated by other methods, such as estimates from X-ray reflection (∼40°; Draghis et al. 2023b; ∼48°; Peng et al. 2024) and X-ray polarisation observations (∼30 − 60°; Veledina et al. 2023). Our sample also includes MAXI J1348–630. Although the QPO detected in this source has been classified as type-C, it exhibits properties similar to those of the imaginary QPO (Alabarta et al. 2022, 2025). MAXI J1348–630 was first identified as a low-inclination (∼30 − 40°) system based on reflection studies (Chakraborty et al. 2021). Carotenuto et al. (2022) measured an inclination of ∼29° for the jet of the system with respect to the line of sight. To ensure as consistent a comparison as possible, we used the common measurements of the inclination from X-ray reflection fitting to define the high (≳60°) and low (< 60°) inclination. Figure 6 shows the phase lag of the imaginary QPO between 2 − 5 keV and 5 − 12 keV energy bands⁵ as a function of QPO frequency for all these sources. This figure shows that the imaginary QPO phase lag above 2 Hz systematically increases with increasing QPO frequency. However, due to the uncertainty in the measurements of the inclination and the limited sample size, we cannot draw a strong conclusion from the current results on whether the phase lag of the imaginary QPO depends on inclination.

Fig. 6. Phase lag of the imaginary QPO as a function of QPO frequency. The inclination (high or low) of each source is inferred from X-ray reflection fitting. The phase lags are computed between the 2 − 5 keV (reference band) and 5 − 12 keV (subject band) energy bands. For AT2019wey, MAXI J1348–630, MAXI J1820+070, and Cygnus X–1, we use NICER data, while for Swift J1727.8–1613 we use XMM-Newton data.

5. Conclusions

We analysed all available archival NICER observations of AT2019wey during the LHS and HIMS by the jointly fitting the PDS and CS. We also compared our results with those obtained for other sources and with properties of type-B and type-C QPO. Our main conclusions are as follows:

1. We detect a QPO signal in AT2019wey that is not always significant in the power density spectrum but is significant in the imaginary part of the CS. The QPO is always accompanied by a sharp increase (‘cliff’) in the phase-lag frequency spectrum and a sudden drop (‘dip’) in the coherence-function frequency spectrum. We refer to this QPO as the imaginary QPO. Our results add a new member to the emerging family of imaginary QPOs. Compared to other sources with imaginary QPOs, we extend the occurrence of this phenomenon to the LHS.

2. We find a decrease in the normalised covariance of the imaginary QPO above ∼2 keV, which may be related to a relatively low electron temperature of the Comptonised component in AT2019wey.

3. We find that the phase-lag spectra of the imaginary QPO in AT2019wey during the HIMS show a U-shaped profile, while the U-shape disappears in the LHS.

4. The imaginary QPO in AT2019wey shares several properties with type-C QPOs. We suggest that the imaginary QPO in AT2019wey could be the type-C QPO.

5. Whether the phase lags of the imaginary QPOs show an inclination dependence – as seen in type-C QPOs – remains inconclusive based on the current sample size.

While dozens of BHXBs are known, imaginary QPOs have been detected only in a few. These QPOs typically appear during the soft-to-hard transition or along the low-accretion-rate branch of the HID. However, their absence in similar transitions of sources such as GX 339–4 raises the question of whether this discrepancy may stem from factors such as system inclination and disc-corona geometry. Future work should therefore focus on: (1) expanding the sample size of sources showing these QPOs to test for inclination dependence and (2) systematically comparing disc-corona geometries across different sources to understand why these QPOs occur in some cases but not others.

Acknowledgments

We are grateful to the anonymous reviewer for the comments that helped us improve the manuscript. FG is a CONICET researcher. FG acknowledges support from PIP 0113 and PIBAA 1275 (CONICET). FG was also supported by grant PID2022-136828NB-C42 funded by the Spanish MCIN/AEI/10.13039/501100011033 and “ERDF A way of making Europe”. PY acknowledges support from the China Scholarship Council (CSC), No. 202304910059. PY also thanks Kevin Alabarta for the discussion on the ‘rising-falling’ profile of the rms/covariance spectrum.

References

Appendix A: Methodology

A.1. Multi-Lorentzian model

We adopt the method proposed by Méndez et al. (2024) to search for signals in the CS and characteristic structures in the phase-lag frequency spectrum and coherence function. This method assumes that the PDS and CS of low-mass X-ray binaries consist of an additive combination of Lorentzian functions, each of which is coherent across different energy bands but incoherent with one another (see also Nowak et al. 1999b). The technique involves the simultaneous fitting of the PDS and the real and imaginary parts of the CS. It has been shown mathematically that this method enhances the detectability of subtle QPOs, which exhibit a higher S/N in the CS than in the power spectrum. Moreover, the coherence function can reveal signals that are difficult to detect in either the power or cross spectra (e.g. Méndez et al. 2024; Bellavita et al. 2025; Rout et al. 2025a).

Under the above assumptions, if x(t) and y(t) are correlated light curves in two energy bands, the PDS in the two energy bands can be modelled as

$$\begin{array}{c}\hfill \begin{array}{cc}& G_{\mathit{xx}}(\nu )={\displaystyle \sum _{i=1}^n}{G}_{xx,i}(\nu )={\displaystyle \sum _{i=1}^n}{A}_{i}L(\nu ;{\nu }_{0,i},{\mathrm{\Delta }}_i)\hfill \\ \hfill & G_{\mathit{yy}}(\nu )={\displaystyle \sum _{i=1}^n}{G}_{yy,i}(\nu )={\displaystyle \sum _{i=1}^n}{B}_{i}L(\nu ;{\nu }_{0,i},{\mathrm{\Delta }}_i),\hfill \end{array}\end{array}$$(A.1)

where the L(ν; ν_0, i, Δ_i) is the Lorentzian functions with centroid frequency ν_0, i and FWHM Δ_i, and A_i and B_i are the normalisations in the respective energy bands. We take x as the reference band and y as the subject band. If each Lorentzian is fully coherent in two energy bands, and any two Lorentzians are incoherent, the Real and Imaginary part of the CS between the two light curves can be modelled as a linear combination of Lorentzians:

$$\begin{array}{c}\hfill \begin{array}{cc}& \mathrm{Re}[G_{\mathit{xy}}(\nu )]={\displaystyle \sum _{i=1}^n}{C}_{i}L(\nu ;{\nu }_{0,i},{\mathrm{\Delta }}_i)cos[\mathrm{\Delta }{\varphi }_{xy,i}(\nu )+\pi /4]\hfill \\ \hfill & \mathrm{Im}[G_{\mathit{xy}}(\nu )]={\displaystyle \sum _{i=1}^n}{C}_{i}L(\nu ;{\nu }_{0,i},{\mathrm{\Delta }}_i)sin[\mathrm{\Delta }{\varphi }_{xy,i}(\nu )+\pi /4],\hfill \\ \hfill \end{array}\end{array}$$(A.2)

where $C_i=\sqrt{A_i{B}_i}$ and the presumed phase lag between the two signals, Δϕ_xy, i(ν) = g_i(ν; p_j, i), is frequency-dependent. The additional constant phase shift of π/4 accounts for the rotation of the cross vector. As explained in Méndez et al. (2024), g_i(ν; p_j, i) can be any arbitrary function of frequency with parameters p_j, i for each Lorentzian; here we take the simplest case of constant phase lag, g_i(ν) = 2πk_i ∈ [ − π, π), for each Lorentzian, where the k_i are obtained from the fitting of the data. The total phase lag and the coherence function can be computed from Eq. A.2 (see more details in Méndez et al. 2024)

A.2. Simultaneous fitting of PDS and CS

Using XSPEC 12.15.0f (Arnaud 1996) we first simultaneously fitted the PDS in the 0.3 − 2 keV and 2 − 12 keV energy bands as well as the rotated Real and Imaginary parts of the CS for the selected 16.384 s segments. To obtain the optimal combination of Lorentzian components, we add them one by one until no significant residuals remain. A Lorentzian is considered to be required in the model if its normalisation in either PDS or the CS is at least the 3σ significant. In zone 8, shown in Appendix C.1, the Lorentzian for the imaginary QPO is detected at a 90% confidence level. During the fitting process, we allow the centroid frequency and FWHM of each component free to vary, but we link them to be the same in the PDS and the CS. In contrast to $C_i=\sqrt{A_i{B}_i}$ in the model, in our fitting process we let the normalisation of the CS and PDS free to vary independently, and check afterwards whether they are consistent with this relation (Table 1). Finally, from the best-fit model, we derive the models for the phase lags and the coherence function, and compare them with the observed data. As the PDS is normalised in fractional rms units, integrating the best-fit model over a fixed frequency range yields the total fractional rms. Similarly, the square root of the Lorentzian normalisation represents the fractional rms amplitude of the corresponding component.

To investigate the energy-dependent properties of the variability components, from which we derive the fractional rms and phase-lag spectrum, we simultaneously fitted the PDS in the full band (reference band), the PDS in each narrow band (subject bands), and the corresponding cross spectra, using the best-fit model obtained above. In this procedure, we assumed that the centroid frequencies and FWHMs of all Lorentzian components are energy-independent and are therefore linked to be the same.

A.3. Normalised covariance

As described by Wilkinson & Uttley (2009) and Uttley et al. (2014), when two signals are highly coherent over a certain frequency range, the covariance effectively acts as a matched filter, using variability in a high signal-to-noise band to detect much weaker correlated variability in the energy channel of interest. The linear multi-Lorentzian model assumes that each Lorentzian function in different energy bands is coherent, so

$$\begin{array}{c}\hfill {\gamma }_{xy,i}^2(\nu )=\frac{|G_{xy,i}{(\nu )|}^2}{G_{xx,i}(\nu )G_{yy,i}(\nu )}=\frac{{C_i}^2}{A_i{B}_i}=1,\end{array}$$(A.3)

where the i is for the i-th Lorentzian component. We can define a normalised the covariance of each component following Wilkinson & Uttley (2009):

$$\begin{array}{c}\hfill Co{v}_{i,norm}:=\sqrt{\frac{|G_{xy,i}{(\nu )|}^2}{G_{xx,i}(\nu )}}=\frac{C_i}{\sqrt{A_i}}.\end{array}$$(A.4)

Therefore, when the index i corresponds to the imaginary QPO component, the normalised covariance is formally equivalent to the rms of the imaginary QPO, $\sqrt{B_i}$. On the one hand, C_i is the normalisation of the i-th Lorentzian in the model for the CS that has better S/N than the PDS. On the other hand, A_i is the normalisation of the model for the 0.3 − 12 keV reference-band PDS that always provide high S/N.

Appendix B: Remarks on the second imaginary QPO in the CS and weak hard lags below 1 Hz

B.1. Second imaginary QPO

In Fig. 2, we identify an additional relatively narrow feature with Q = 1.7 at a higher frequency of 3.75 Hz in zone 5, located immediately next to the imaginary QPO. We refer to this Lorentzian as the second imaginary QPO. The same component is also visible in some of the other zones (see Appendix C.1). The QPO adjacent to the imaginary QPO that appears as a shoulder of that QPO also shows a relatively large imaginary component, but its centroid frequency does not align with the position of the cliff and dip. Shoulder-like features also appear in Cygnus X–1 (Fogantini et al. 2025, see their Fig. 2) and MAXI J1820+070 (Bellavita et al. 2025, see their Fig. 3), although they are not reported. In Swift J1727.8–1613 (Brigitte et al. 2025, see their Fig. 2), they modelled the PDS and the rotated CS using three Lorentzian components, but the residuals indicate that an additional component is still required at frequencies slightly higher than that of the imaginary QPO.

B.2. Weak hard lags below 1 Hz

We also note that in the phase-lag frequency spectrum, the phase lags neither reache zero nor remain constant at frequencies below ∼ 1 Hz, where weak hard lags are observed. This feature is pronounced in other zones (see Appendix C.1), and moves towards lower frequency as source hardens. The phase-lag bump over the low-frequency band is also present in other sources, such as Cygnus X–1 (Fogantini et al. 2025, see their Fig. 2), MAXI 1820+070 (Bellavita et al. 2025, see their Fig. 3) and Swift J1727.8–1613 (Brigitte et al. 2025, see their Fig. 2).

Appendix C: Best fitting for zones 2–9

We simultaneously fitted the PDS (0.3 − 2 keV and 2 − 12keV) and the CS (real and imaginary parts) of AT2019wey in zones 2–9 defined in Sect.3.1, covering the LHS and HIMS. These zones correspond to the intervals where a pronounced cliff in the phase lags and a dip in the coherence function are observed. We applied the technique introduced by Méndez et al. (2024), assuming a constant phase-lag model. Figure C.1 presents the best-fitting PDS in the two energy bands, together with the real and imaginary parts of the rotated CS for each selected zone. The corresponding model predictions for the phase lags and coherence function are also shown. Table C.1 shows the best-fit parameters for those zones. We note that the frequencies of several common components evolve in the same direction as the imaginary-QPO frequency, indicating that there is correlation between the frequency of the imaginary QPO and the other timing components.

Fig. C.1. Best fit of the PDS (0.3 − 2 keV and 2 − 12 keV) and the corresponding real and imaginary parts of the rotated CS for zones 2-9, and the derived model for the phase lags and coherence function.

Fig. C.1. continued.

Appendix D: Phase lags and coherence function between 2–4 keV and 4–12 keV energy bands

König et al. (2024) first reported that in Cygnus X–1 the increase in the phase lags and the drop in the coherence usually disappear when the reference band is above ∼ 2 keV (see also Fogantini et al. 2025). The same behaviour was also observed in MAXI J1820+070 (Bellavita et al. 2025). Fig. D.1 shows the phase-lag frequency spectrum and coherence function of AT2019wey when we use a reference band above 2 keV. In zones 2–6, an increase in the phase lags is apparent, whereas the corresponding drop in the coherence function is not clearly discernible. We also apply the multi-Lorentzian model to fit simultaneously the PDS in the 2 − 4 keV and 4 − 12 keV bands as well as the corresponding CS. In zones 2–6, we can fit both the PDS and CS with a model consisting of four Lorentzian components. The phase-lag frequency spectra can be well described by the model derived from the corresponding best-fit model for both the PDS and CS. We find that a high-frequency broad component with a large imaginary part is responsible for the increase in the phase lags. In zone 3, where the cliff–dip pair is the most apparent, adding a Lorentzian component at the possible dip frequency to both the PDS and CS reproduces the dip in the coherence function, although the component remains very weak and only marginally significant.

Fig. D.1. Phase lags and coherence function between the 2 − 4 keV and 4 − 12 keV energy bands. Note: we rebin the frequency spectrum by an factor of 101/50 to improve the S/N.

All Tables

All Figures

Fig. 1. HID constructed from the NICER observations of AT2019wey between August 9, 2020, and May 7, 2025. Each data point represents a 16.384 s segment. Main panel: Ten colour-coded regions corresponding to the data groups defined by hardness intervals and used in this work. The grey squares, grey stars, and grey triangles mark the three epochs analysed by Yao et al. (2021a), identified as the LHS (grey squares) and HIMS (grey stars and triangles), respectively. The vertical dashed line at HR = 0.21 roughly separates the LHS from the HIMS. Inset panel: Data used in our analysis (black points; same as that in the main panel), along with data excluded from our analysis (grey points).

In the text

Fig. 2. Simultaneous fitting of the PDS and CS of AT2019wey with a multi-Lorentzian model assuming constant phase lags, using the combined NICER data in zone 5 as an example. Top panels: Best fits (solid black lines) and residuals for the PDS (left panel) in the 0.3 − 2 keV (blue) and 2 − 12 keV (red) bands, and to the real (blue) and imaginary (red) parts of the CS rotated by 45° (right panel). Middle panels: Real (left) and imaginary (right) parts of the non-rotated CS, showing a linear scale on the y-axis. Each model consists of six Lorentzian components; the one corresponding to the imaginary QPO is highlighted with thicker lines. Bottom panels: Phase-lag (left) and coherence-function (right) frequency spectra, along with their residuals with respect to the model derived (solid black curves) from the simultaneous fitting of PDS and non-rotated CS. The dashed black curves show the derived model without the imaginary QPO and without refitting. We rebinned the plots by a factor of 5 to improve the visibility of the dip in the coherence function.

In the text

	Fig. 3. Frequency of the imaginary QPO vs the 0.1 − 20 Hz broadband fractional rms in the 0.3 − 2 keV (circles), 2 − 12 keV (triangles), and 0.3 − 12 keV (squares) energy bands derived from our model. Each point represents the observation in a zone. The colour bar indicates the evolution of the HR across all zones. The two grey arrows separate the observations in the LHS (left) and HIMS (right).
In the text

	Fig. 4. Phase lags and coherence function for each region of the HID of AT2019wey for the 0.3–2 keV and 2–12 keV bands. The vertical dashed black line marks the frequency of the imaginary QPO (see text) where the phase-lag cliff and coherence-function dip occur. The colour coding used in the frequency spectra for the different zones matches that of Fig. 1.
In the text

	Fig. 5. rms (red), normalised covariance (blue), and phase-lag (orange) spectra of the imaginary QPOs in AT2019wey for different zones in the HID. All quantities are derived from the fits using our multi-Lorentzian model. The error bars represent 1σ uncertainties. The arrows in the rms spectra represent the 95% upper limit.
In the text

Fig. 6. Phase lag of the imaginary QPO as a function of QPO frequency. The inclination (high or low) of each source is inferred from X-ray reflection fitting. The phase lags are computed between the 2 − 5 keV (reference band) and 5 − 12 keV (subject band) energy bands. For AT2019wey, MAXI J1348–630, MAXI J1820+070, and Cygnus X–1, we use NICER data, while for Swift J1727.8–1613 we use XMM-Newton data.

In the text

	Fig. C.1. Best fit of the PDS (0.3 − 2 keV and 2 − 12 keV) and the corresponding real and imaginary parts of the rotated CS for zones 2-9, and the derived model for the phase lags and coherence function.
In the text

	Fig. C.1. continued.
In the text

	Fig. D.1. Phase lags and coherence function between the 2 − 4 keV and 4 − 12 keV energy bands. Note: we rebin the frequency spectrum by an factor of 101/50 to improve the S/N.
In the text