© The Authors 2025

1. Introduction

A large fraction of the transient sources can be associated with systems containing compact objects. Being able to study the emission of these systems simultaneously at multiple wavelengths and on multiple timescales has proven to be of fundamental importance to correctly interpreting the physical phenomena that are still poorly understood and estimating the fundamental properties of the compact objects that power them.

MAXI J1820+070 is a bright X-ray black hole transient that was discovered on 6 March 2018 in the optical band (ASASSN-18ey; Denisenko 2018; Tucker et al. 2018) by the All-Sky Automated Survey for SuperNovae (ASAS-SN; Shappee et al. 2014) and six days later in the X-ray band (Kawamuro et al. 2018) by the Monitor of All-sky X-ray Image (MAXI; Matsuoka et al. 2009)¹. The main outburst lasted several months, roughly until the end of 2018, and it was followed by a series of rebrightenings that were recorded until June 2020 (plus some weak activity reported in March and April 2021; Baglio et al. 2021a,b). During the long outburst, the source was widely observed at different wavelengths, from the radio band up to the X-rays. No detection was reported at γ-ray energies (Abe et al. 2022).

The high luminosity attained by MAXI J1820+070, due in part to the relatively small distance of the source from Earth (d = 2.96 ± 0.33 kpc; Atri et al. 2020), and the complex phenomenology displayed during the outburst made this source a favoured target of many observational campaigns at different wavelengths (Kalemci et al. 2022). In the first months of the outburst, the source went through all the typical states and transitions of a black hole accreting X-ray binary: the low/hard state (LH), high/soft state (HS), and intermediate state (IM; as reported in detail by Shidatsu et al. 2019). For a review on the different accretion states of X-ray binaries, we refer to Remillard & McClintock (2006) and Done et al. (2007), for example.

The mass of the black hole was estimated from the spectral shifts and the ellipsoidal modulation of the optical emission of the companion star during quiescence ($M_{\mathrm{BH}}=8.{48}_{-0.72}^{+0.79}{\mathrm{M}}_{\odot }$; Torres et al. 2020). Some of the jet properties, such as the inclination angle (θ_{_JET} = 63 ± 3 deg) and the velocity (v = 0.89 ± 0.09c), have been estimated from European VLBI radio measurements (Atri et al. 2020). Other estimates of the jet properties derived from modelling the highly correlated emission between different bands confirmed the existence of a high relativistic and confined jet ($\mathrm{\Gamma }=6.{81}_{-1.15}^{+1.06}$, $\varphi =0.{45}_{-0.11}^{+0.13}$ deg) that carries away a good fraction of the total power (∼0.6L_{_1-100keV}; Tetarenko et al. 2021).

The properties of these systems can also be studied using fast timing features observed in the Fourier domain. The power density spectra (PDS) often reveal broad and narrow features related to particular phenomena occurring in the accretion flow. Variability in narrow frequency ranges are usually referred to as quasi-periodic oscillations (QPOs), and they can be divided into low-frequency QPOs (LFQPOs; ∼0.005 − 40 Hz) and high-frequency QPOs (HFQPOs; ∼40 − 450 Hz). The former are further divided into different types based on their properties (type-A, type-B, and type-C QPOs). For a description of the different types of QPOs, we refer to Ingram & Motta (2019), for example, and references therein. When MAXI J1820+070 was in the LH state, a type-C LFQPO with characteristic frequencies below 1 Hz was detected in the X-ray band with Swift/XRT and NICER (Neutron star Interior Composition Explorer) as reported by Stiele & Kong (2020). Around the time of the first state transition (from the LH to the HS state), the LFQPO switched quickly from type-C to type-B together with the emission of a strong radio flare. These facts were interpreted in terms of the launch of a superluminal jet (Bright et al. 2020; Homan et al. 2020). Broad-band noise features in the PDS of MAXI J1820+070 were also reported and interpreted again in terms of a truncated accretion disc (Dziełak et al. 2021; Yang et al. 2022).

The LFQPOs synchronous to the X-ray QPOs have also been observed at optical wavelengths with various instruments (Zampieri et al. 2019; Mao et al. 2022; Thomas et al. 2022a). A time delay of ∼165 ms was found cross-correlating the X-ray and optical light curves, with the optical lagging the X-ray (Paice et al. 2019). These findings are compatible with the hypothesis that the emission in the two bands is produced in close-by regions, in the inner accretion flow, and/or at the base of a precessing jet (Paice et al. 2019; Thomas et al. 2022a). In a later work, analysing further simultaneous optical/X observations, Paice et al. (2021) found that the emission can be separated into two distinct synchrotron components originating from a compact jet and a hot flow, respectively.

Looking at the variability at longer timescales, Patterson et al. (2018) reported a ∼16.87 h optical modulation appearing some days before the transition to the HS state (corrected to ∼16.57 h in a later work; Patterson 2019). This modulation was interpreted as a super-orbital motion (superhump) emission since it differs from the orbital ones (∼16.45 h, Torres et al. 2019). Superhumps are often observed in cataclysmic variables (CVs) and are large-amplitude photometric oscillations with a period that is longer than the orbital period of the binary system (the phenomenon of superhumps in CVs is explained in detail, for example, in Kato et al. 2009, 2017). Also, in Niijima et al. (2021), where a detailed photometric study of MAXI J1820+070 is presented, these oscillations are interpreted as being a process similar to the superhumps seen in many CVs. Thomas et al. (2022b) proposed that the large optical modulation in MAXI J1820+070 could be caused by a warped precessing accretion disc (Ogilvie & Dubus 2001).

In principle, the quality of the data available for MAXI J1820+070 allows us to perform rather accurate estimates of the spin of the black hole. However, diverse methods bring us to quite different results. Using the continuum-fitting method with a thin disc model (Zhang et al. 1997) and the X-ray soft-state spectra taken with Insight-HXMT, (Hard X-ray Modulation Telescope; Zhang et al. 2020), Guan et al. (2021), and Zhao et al. (2021) found a relatively slow-spinning black hole with spin parameter a_* equal to $0.2_{-0.3}^{+0.2}$ and 0.14 ± 0.09, respectively. On the other hand, using some characteristic frequencies in the PDS (Motta et al. 2014) and the relativistic precession model (RPM; Stella & Vietri 1998), Bhargava et al. (2021) estimated a fast-spinning black hole with $a_{\ast }=0.{799}_{-0.015}^{+0.016}$. Further uncertainties on the validity of the measurements of the spin of the black hole in MAXI J1820+070 also come from optical polarimetric observations, which led a lower limit of 40° to be placed on the angle between the spin and orbital axes (Poutanen et al. 2022).

The purpose of this work is to try to provide the most comprehensive view of the variability in two different bands (optical and X-ray) and on different timescales of the 2018–2020 outburst of MAXI J1820+070. To this end, we made use of public and proprietary optical and X-ray data from various telescopes and resources with different time resolutions.

The paper is structured as follows. In Sect. 2, we briefly describe the observations and the data reduction process. In Sect. 3, we show the analysis made on the data, starting from the analysis of the low-cadence data then moving to the sub-second timescale and finally to the spectroscopic data. In Sect. 4, we discuss our results, and in Sect. 5 we report our conclusions.

2. Observations and data reduction

For the low-cadence photometric data, we used observations taken with the Schmidt telescope in Asiago, the telescopes of the ANS Collaboration² (Munari et al. 2012), and observations from the AAVSO³ (American Association of Variable Stars Observers, Kafka 2021). For the sub-second optical variability, we used our own observations collected with the fast photon counters IFI+Iqueye (Naletto et al. 2009) and Aqueye+ (Zampieri et al. 2015) in Asiago. Finally, for the X-ray light curve and fast variability, we exploited the rich dataset collected with the NICER satellite (Gendreau et al. 2012, 2016). We also used a series of low resolution spectroscopic data taken with the Asiago 1.22m telescope. The following subsections report a summary of the entire optical and X-ray dataset used in this work.

2.1. Optical photometry

In Table A.3 we have summarised the log of the photometric observations of MAXI J1820+070. The observations span almost 2.5 years (from March 2018 to August 2020) and the dataset comprises ∼2 × 10⁵ photometric measurements. We have obtained the optical photometry of MAXI 1820+070 in the UBVRI system defined by the equatorial standards of Landolt (1992) and in the g′ r′i′ Sloan bands as given in the APASS All-Sky Survey (Henden & Munari 2014; Henden et al. 2018). The observations come from multiple facilities in Italy, Slovakia, Crimea and Russia. The telescopes involved in the observations are reported in Table A.2.

Data reduction has involved all the usual steps for bias, dark and flat, with calibration images collected during the same observing nights. We adopted aperture photometry because the sparse field around MAXI 1820+070 did not require PSF-fitting procedures. The transformation from the local to the standard system was carried out via nightly colour equations calibrated on a photometric sequence recorded on the same frames and extracted from the APASS survey, ported to the Landolt system via the transformations calibrated by Munari et al. (2014). The final errors are the quadratic sum of the Poissonian error on the variable and the error in the transformation to the standard system via the instantaneous colour equations. Typical errors are of the order of few hundredths of magnitudes.

We included in the analysis a part of the 2018 data from the AAVSO (Kafka 2021). These data were used to complement those taken with our instrumentation, especially for filling periods when the sampling was scanty. These data are already reduced and calibrated and can be freely downloaded from the AAVSO website.

Using the entire cleaned and calibrated dataset, we searched the periodicities in the light curve of MAXI J1820+070 as described in Sect. 3.1.

2.2. Optical spectroscopy

Low resolution spectroscopy of MAXI J1820+070 was obtained with the 1.22m telescope + B&C spectrograph operated in Asiago by the Department of Physics & Astronomy of the University of Padova. The CCD camera is a ANDOR iDus DU440A with a back-illuminated E2V 42 − 10 sensor, 2048 × 512 array of 13.5 μm pixels. It is highly efficient in the blue down to the atmospheric cut-off around 3200 Å, and it is normally not used longward of 8000 Å for the fringing affecting the sensor. The adopted 300 ln/mm grating, blazed at 5000 Å, allowed us to cover the wavelength range from ∼3200 to ∼8000 Å at a spectral dispersion of 2.31 Å/pix. The slit width was set to 2 arcsec, providing a resolution with a full-width half-maximum (FWHM) of 2.2 pix. The slit was always aligned with the parallactic angle for optimal absolute flux calibration, which was achieved by observations of the spectro-photometric standard HR 6900. This standard is conveniently located just a few degrees away from MAXI J1820+070, and has a similarly hot energy distribution and blue colours, to the benefit of a higher quality of the flux calibration. The BVR magnitudes were computed by band-profile integration on all recorded spectra, and were checked for consistency against nearly simultaneous CCD photometry collected with the aim of building the light and colour-curves of MAXI J1820+070 described in this paper. MAXI J1820+070 was observed at 10 epochs distributed through the 1st to the 5th maxima of the optical light curve, each time for about 2700 s and the recorded spectra are presented in Fig. 10 were observing dates are provided.

2.3. High timing resolution optical observations

We performed 8 observational runs of MAXI J1820+070 with the IFI+Iqueye (IQ) and Aqueye+ (AQ+) fast photon counters (Barbieri et al. 2009; Naletto et al. 2009; Zampieri et al. 2015) mounted at the 1.2 m Galileo telescope and the 1.8 m Copernicus telescope (Asiago, Italy) from April to October 2018 (summary of observations in Table A.1). The source was observed in white light (without filter). The data were reduced with the QUEST software (v. 1.1.5, see Zampieri et al. 2015). Since the instruments can only observe one source at a time (due to having a field of view of a few arcseconds), between on-target acquisitions, we also regularly observed the sky and a nearby star. The selected star GSC 00444-02282 was observed for calibration purposes (to be able to properly compare data taken on different nights) and as a reference to check for possible systematics in the PDS of MAXI J1820+070. The data taken on source during a night are divided into segments of ∼30/60 minutes each, while a typical observation of the reference star is of the order of 15 minutes.

The final output of the low-level analysis consists of a sequence of time series that are subsequently time-binned to produce the light curves that can be used for the analysis that follows. The selected value of the bin size is 1 ms.

2.4. X-ray observations

During 2018, NICER, actively monitored MAXI J1820+070 via the X-ray timing instrument (XTI; Gendreau et al. 2012), with almost one observation per day between 12 March and 21 November (210 observations in 254 days). We reduced the data between 12 March and 15 October 2018 (ObsID 1200120101-1200120278). Observations usually have exposure times greater than 1 ks (one observation reaches ∼22 ks) for a total exposure of ≳400 ks. We applied the standard data processing procedure⁴, using the script nicerl2 - part of HEASoft (v.6.29) software⁵ – with version 20210707 of the calibration files. The selected energy range for the analysis is 1 − 12 keV and the chosen time bin is 1 ms. We finally barycentred the data with the script barycorr.

X-ray data were used for comparison with optical results, both for the photometry data (rebinning the X-ray light curve with longer time bins) and for the high time resolution data.

3. Data analysis and results

3.1. Analysis of the optical photometry and the X-ray light curve

3.1.1. Optical photometry

Figure 1 shows the optical light curve of MAXI J1820+070 in all the available bands. The data start around the time of the first maximum and extend for more than 2 years. 5 subsequent rebrightenings can be seen after the main burst (the vertical dashed lines show approximately the time of the maximum for each burst). The intervals between the maxima of all bursting events are: 94 days between the first and second; 110 days between the second and third; 153 days between the third and fourth; and finally 186 days between the last two bursts. The bursting activity is observed in all bands with similar trends. Interestingly, the intervals between two subsequent bursts increase progressively.

Fig. 1. Light curves of MAXI J1820+070 in different optical bands since the epoch of the first maximum after the discovery of the source and until August 2020. The data taken in the Landolt UBVRI bands (Landolt 1992) are shown in violet, blue, light green, light red, and light orange, respectively. The darker green, red, and orange represents the g′, r′, and i′ Sloan bands of the APASS system (Henden & Munari 2014; Henden et al. 2018). The dashed black vertical lines indicate approximately the dates when the local maximum of the optical luminosity is reached during all the subsequent bursts. The red vertical lines indicate the epochs when we observed optical LFQPOs with Iqueye. The yellow vertical lines indicate the dates when Torres et al. (2019) measured the periodicity of the binary system (1.4591 1/d) spectroscopically. The grey shaded area indicates the interval in which a photometric periodicity is detected in this dataset, while the orange shaded area indicates the interval in which a periodicity can be seen also in the X-ray data.

In order to search for possible periodicities in the data, we first remove the flaring activity from the long-term light curve (which is on timescales much longer than the searched periodicities). The whole procedure was carried out using a single filter, to avoid possible systematics deriving from averaging together different bands. We relied on V-band data, as they have the best temporal coverage. We fitted a spline to the average magnitude computed on N consecutive days. After subtracting the spline (detrending), we used an implementation of the generalised Lomb-Scargle algorithm (LS; Lomb 1976; Scargle 1982; VanderPlas 2018) from the Astropy library (Astropy Collaboration 2018) to finally search for periodicities in the data (see Fig. 5). It is worth nothing that a periodicity appears in the periodogram even before applying the detrending procedure, with several spurious peaks caused by the long trend flaring activity.

We initially searched for periodicities on the entire time interval using N = 1, 2, 3, or 4 days to compute the average magnitude. Results consistently showed a clear peak at a frequency around 1.45 d⁻¹. However, folding the whole dataset with this frequency gives a very noisy profile. We therefore tried to restrict the analysis to different temporal sub-intervals to find out when the most prominent periodic signal is generated. We checked that the detrending procedure did not introduce spurious artefacts in the light curve and searched for the interval with the best signal. In the following we describe the adopted iterative procedure:

1. We first clean the data computing the average over two days.

2. We then searched for a periodicity, computing the LS diagram and looking for the most prominent peak in time windows of 100 days, from the first available date (MJD 58189) to the end (MJD 59066). We found that a peak in the LS diagram appears when we consider a time window between MJD 58285 and MJD 58385.

3. We computed again the LS diagram changing the starting date in steps of 1 day in an interval of 20 days around MJD 58285. We did the same also for the ending date. In this way, the size of all time windows is approximately the same, and we can then compare the maximum powers (next point).

4. Figure 2 shows the maximum power of the peak for different starting and ending dates. We selected the most favourable interval considering as lower bound the date when the power starts rising and as upper bound the date when it starts decreasing. Before MJD 58283 the power is somehow constant (it varies very little) and after that date it starts rising quickly, reaching a maximum at MJD 58285. After MJD 58391 the power begins to decrease slowly and monotonically. We then selected the interval MJD 58283-58391⁶.

   Fig. 2. Power of the highest peak in the LS diagram of the V-band light curve after detrending the data and considering different starting and ending dates. This plot was computed to understand which temporal interval to use for searching the best timing solution of the system by means of fitting the sinusoid from Eq. (1) to the data. The power increases from blue to red. A vertical straight line marks the time (MJD 58283) at which the power starts to rise, while a horizontal straight line marks the time (MJD 58391) at which the power starts to decrease.

5. After selecting the most favourable interval we returned to step (1) and varied the number of days to average together the data. We detrended the data using different number of days, from 1 to 20 days. With the different detrended data we computed the LS diagram in the interval MJD 58283-58391 and determined the frequency and the power of the most prominent peak (Fig. 3). The black line indicates the frequency at which we found the maximum power in the LS after detrending with N averaged days (x-axis). The background colour chart corresponds to the value of the power at the peak. We found that peak power is maximised averaging together 3 days of data.

   Fig. 3. Frequency (black line) and power (background colour chart) of the highest peak in the LS diagram computed after detrending the V-band data using the average over N days (x-axis). The maximum power is reached for N = 3 days. The spread of the frequency is compatible with the uncertainties of the measurement.

6. After fixing N = 3 days, we repeated steps (3) and (4) to verify the choice of the interval (confirming that the best interval seems to be between MJD 58283 and MJD 58391).

7. Finally, the data were rebinned in intervals of 20 minutes to compensate for the different exposure times in the original data-set and to filter out the noise of the light curve at shorter timescales.

The steps in the procedure outlined above may in principle introduce an additional uncertainty in the estimate of the periodicity. However, this does not seem to be a serious concern as in the fitting procedure reported below the error of the fitted frequency turns out to be of the same order of the dispersion visible in Fig. 3. Figure 4 shows the detrended V-band data in the selected interval, which corresponds to the grey shaded area in Fig. 1. Even after the detrending procedure is applied, the curve shows a quite pronounced intrinsic dispersion that varies inside the interval. The amplitude of the variability is large at the beginning, while in the second half of the interval is smaller. The LS diagram in Fig. 5 shows a strong peak at a frequency of 1.4517 d⁻¹, for both all the V-band data (orange curve) and the sole data between MJD 58283-58391 (blue curve). The peak for the data in the selected interval is much higher, while that for the data outside this interval (green curve) completely disappears, meaning that the periodicity outside this interval is completely absent.

Fig. 4. Detrended optical data in the time interval selected to measure the periodicity. The colour code is the same as in Fig. 1.

Fig. 5. Comparison of LS diagrams in different time intervals. The LS diagram computed considering all the available V-band data is shown in orange, whereas the diagram computed considering only the data within (outside) the interval MJD 58283-58391 is shown in blue (green).

In order to accurately determine the modulation period, we performed a fit of the detrended light curve considering the following sinusoidal function (Δ(t)):

$$\begin{array}{c}\hfill \mathrm{\Delta }(t)=Asin(2\pi f[t-t_0]+\varphi )+C,\end{array}$$(1)

where A is the amplitude, f the frequency, ϕ the phase with respect to the initial time t₀ and C is an offset (to account for possible residuals in the detrending procedure). To fit the data to Eq. (1) and estimate the uncertainties of the fitted parameters, we implemented a Markov Chain Monte Carlo procedure (MCMC) through the python package emcee (Foreman-Mackey et al. 2013). The results of the fitting procedure are reported in Table 1. The value of the frequency is in agreement with that measured in the LS diagram computed in the same interval. The error is also compatible with the dispersion of the frequency for different choices for the number of nights used in the detrending procedure (see Fig. 3).

Figure 6 shows the V-band light curve folded using the frequency obtained with the MCMC fitting procedure together with the best-fit model (black line) and the 2σ confidence interval (grey shaded area). The overall behaviour is well described by the fitted profile, with only some outliers due to data within certain time intervals where the oscillation is not perfectly visible (or is slightly out of phase), as can be seen from Fig. 7.

Fig. 6. Data from the V-band between MJD 58283 and MJD 58391 folded with the value of the frequency f reported in Table 1. The dotted black line is the sinusoidal fitted model with the corresponding 3σ confidence uncertainty shown as a grey shaded area. Some residuals can be seen in the plot and are caused by the data within certain time intervals where the oscillation is not perfectly visible or slightly out of phase, as can be seen from Fig. 7.

Fig. 7. Evolution of the amplitude and shape of the periodical signal in the light curve of MAXI J1820+070 between MJD 58283 and MJD 58391 seen in the folded optical data in different sub-intervals. The colour code is the same as in Figs. 1 and 4. The dashed black line corresponds to the best-fitting sinusoid (Eq. (1)) for the entire interval. Except for the first panel, the phase of the signal remains constant, while its shape evolves with time. All the optical bands are in good agreement.

We then detrended the data in the other bands using the same correction applied to the V-band data, properly rescaled for the different mean value of each band. The detrended data for the other bands are also shown in Fig. 4. We tried to fit them using Eq. (1) but the measurements are too scattered to allow for the fitting procedure to converge to a single solution. Therefore, we tried to redo the fit by fixing the frequency to the value found in the V-band data and found that the phases in the different bands are in agreement with those in the V-band within the error bars.

The periodic signal shows a clear evolution of the amplitude and the shape. Looking at the light curves folded with the best-fitting frequency in different time sub-intervals, it is possible to see that the periodicity is not always significant, and in addition the shape of the profile and its amplitude change with time (Fig. 7). The first two sub-intervals (MJD 58283-58288 and MJD 58288-58301) show large amplitude oscillations – even larger than the amplitude computed using the entire time interval – and the first of these two sub-intervals seems to be slightly out of phase. In the following sub-intervals the periodicity is again clearly present, but with a lower amplitude and with some points out of phase (e.g. MJD 58301-58309 and MJD 58353-58382), while in other sub-intervals the periodicity seems to disappear completely (MJD 58326-58338).

3.1.2. X-ray light curve

Similarly to what did for the photometric data, we tried to remove the overall short-term irregular oscillations and search for a hours-to-days periodic signal in the X-ray data. Before detrending the X-ray light curve, we rebinned the data in time to 60 s and subsequently renormalised them according to the number of active NICER focal plane modules, as the intense photon flux caused saturation of the internal telemetry in certain time intervals (Homan et al. 2020). We applied the same detrending procedure described above, with the difference that we used 2 (instead of 3) days for the averaging of the data. This choice was dictated by the fact that the X-ray light curve shows more variability on a daily timescale. Moreover, we had to exclude from the detrending procedure the data in the MJD range 58302-58308 because of the intense X-ray flaring activity.

After the detrending we found a small signal in a sub-interval of those in which we detected the optical periodicity (see Fig. 8). We tested the possibility that the signal is caused by the data cleaning procedure by implementing an ad hoc bootstrap method (something similar to the method described in VanderPlas 2018): (1) we simulated many randomly distributed light curves by resampling the X-ray data but keeping the daily trend; (2) we applied exactly the same cleaning procedure described above; (3) we calculated all the LS diagrams and (4) we finally extracted the probability distribution at each frequency from which we could infer a false alarm probability level (0.01%, red curve in Fig. 8). Remarkably, the resulting frequency is out of the area where we could expect to measure spurious signals because of the cleaning process (frequencies ≤1 d⁻¹) and close to the one found in the optical light curve in a time window where the optical periodicity is very well defined (MJD 58309-58326, upper right panel in Fig. 7). However, the significant variability of the X-ray light curve prevents us from effectively performing an accurate measurement. We then fixed the frequency at the value found in the optical data, leaving the phase, the amplitude, and the offset free to vary. We reported the results of the fit in Table 1. In Fig. 9 we show the folded light curve together with the best-fit function. We also overplotted the same X-ray data binned in phase (the black points; 15 bins per phase) to better visualise the periodical trend in the X-ray data. We measured a reduced χ_ν² value of 0.77 for the fitted sinusoidal model against a reduced χ_ν² value of 1.0 for a constant model, showing that a low significance periodicity is present in the data. This periodicity, if real, shows that the X-ray light curve is approximately in anti-phase with the optical light curve, with the X-ray leading the optical by about 1.1π radians or ∼0.37887 days (∼9.1 hours).

Fig. 8. Lomb-Scargle diagram (black curve) computed with the detrended X-ray data in the interval MJD 58309-58326. The red curve correspond to the 0.01% false alarm probability level (see text for details). The highest peak is quite broad, but is centered at the frequency inferred from the optical photometry and well above the false alarm probability level in the range of frequencies around it.

Fig. 9. Comparison of optical and X-ray folded data in the interval MJD 58309-58326 computed using the best-fitting frequency obtained with all optical data between MJD 58283 and MJD 58391 as reported in Table 1. Data are shown together with the best-fitting sinusoid (dashed black lines). The data in the upper panel correspond to the optical data (coloured dots). Lower panel: X-ray data (grey dots) as well as the same data binned in phase (black dots with error bars, 15 bins per phase). The error bars are computed from the standard deviation of the data falling inside every phase bin.

3.2. Optical spectroscopy

The spectra presented in Fig. 10 are distributed over the first 5 optical maxima exhibited by MAXI J1820+070: the first 7 spectra cover the evolution from the 1st optical maximum to the 3rd optical maximum, the 8th spectrum is taken around the 4th optical maximum, and the last two spectra around the 5th optical maximum (the optical maxima are indicated with dashed black lines in Fig. 1). For a clearer view, the bluest part of the spectra (shortward of 3500 Å, which is much noisier than the rest) has been omitted from the figure, and similarly for the wavelength interval red-ward of 6800 Å, where no significant emission lines are observed.

Fig. 10. Ten spectra of MAXI J1820+070 recorded with the Asiago 1.22m + B&C telescope. From top to bottom, the spectra are arranged in order of observing date. They have been offset for plot clarity by the indicated amount. Next to the offset is listed the B-band magnitude exhibited by MAXI J1820+070 at the time of the spectroscopic observation. The brightest emission lines are identified as well as a sample of the interstellar features seen in absorption (CaII, NaI, and the diffuse interstellar band at 6281 Å as representative of several other well visible diffuse interstellar bands).

Irrespective of which maximum they belong, the slopes of the spectra in Fig. 10 can be clearly divided in two groups depending on the brightness of MAXI J1820+070: when the magnitude of the object is B ≤ 13.58, the spectral slopes are bluer (< (B − R) >  = +0.33) than those when the object turns fainter at B ≥ 14.26 (< (B − R) >  = +0.42). This change in the slope also corresponds to the time of the transition from the HS state to the IM state (and soon after to the LH state). The continuum at λ ≤ 3700 Å shows a great range of variability, which is unrelated to the flux in Balmer emission lines. This suggests that an interpretation in terms of a Balmer continuum going back and forth between emission and absorption seems unlikely.

Even if our spectra are characterised by a low-frequency resolution, the broad and double-peak profiles of emission lines are well resolved. The velocity separation between the blue and red peaks of the emission line profiles, and their relative intensity, vary according to the observing date and the specific line, with a mean value around ∼850 km s⁻¹ as illustrated in Fig. 11.

Fig. 11. Zoomed-in view from Fig. 10 of the MAXI J1820+070 spectra for 2018 covering hydrogen Hβ, HeII 4686 Å, and the broad 4640 Å blend attributed to NIII lines. The spectra are offset by the indicated quantity for clarity of the plot. The typical velocity separation of 850 km s−1 for the double-peaked profiles is marked.

The lower ionisation/excitation emission lines, exemplified in Fig. 12 by Hβ and HeI 5876 Å, vary in phase with the level of continuum emission, which appear not to be the case for the higher ionisation HeII 4686 Å line. The equivalent width of the HeII line exceeds 2 times that of the Hβ line for the earliest and brightest spectra in Fig. 11, while it is 1.5 times of the Hβ line or the later and fainter ones. A similar behaviour is exhibited by the broad blend centered at 4640 Å, evolving in parallel with HeII. This suggests that the blend could be due to the Bowen (1934) excitation mechanism, as observed in the high-density environments of novae and symbiotic binaries among other types of objects: HeII Lyman-α 303.78 Å photons are absorbed by OIII in its ground state, that emits at 374.43 Å upon returning to it. The 374.43 Å photons are absorbed by NIII in its ground state, and the following de-excitation produces a trio of lines around 4640 Å (multiplet #2 at 4634.16, 4640.64, and 4641.92 Å). In support of such a scenario comes the abnormally large intensity of Hδ in Fig. 10, which violates the usual Hβ:Hγ:Hδ:Hϵ progression. The abnormal intensity would be easily explained by contribution from NIII lines (multiplet #1 at 4097.31 and 4103.37 Å) emitted during the same return to ground-state following the pumping by absorption of 374.43 Å photons emitted by OIII. The 4640 Å blend is also prominent in the spectra obtained by Muñoz-Darias et al. (2019) at epochs preceding our first spectroscopic observation.

Fig. 12. Bottom panel: Evolution of the equivalent width (in Å) of representative emission lines as a function of time, calculated from the spectra of MAXI J1820+070 in Fig. 11. Top panel: Ratio between the equivalent widths of two of these lines, HeII 4686 Å and Hβ. The colour bar on the top indicates the value of the B band magnitude of the source at the time of the observations. We show for comparison also the periods where the optical (grey shaded area) and the X-ray (orange shaded area) modulation are seen in the light curves, together with the epochs of the different accretion states that the source entered (Shidatsu et al. 2019).

3.3. Analysis of the high timing resolution data

Measuring the variability of a signal on different timescales can be achieved efficiently by moving to the Fourier space (as for the LS diagrams in the previous section) and calculating the power density spectrum (PDS) to measure the variance of a signal at each Fourier frequency and pinpoint the presence of periodic or semi-periodic features in the light curves. We then searched for fast periodic features in the NICER X-ray data and the IQ/AQ+ optical data by calculating the PDS using fast Fourier transform (FFT) algorithms.

We computed the average PDS for each observing run, both in the X-ray and optical band, with a frequency resolution that varies from dataset to dataset (typically between 5 mHz and 20 mHz) and with a constant Nyquist frequency of 0.5 kHz (all light curves are time-binned at 1 ms). The frequency resolution is not constant because it depends on the width of the time window into which the observations are segmented before performing the FFTs, typically of the order of ∼50, ∼100 or even ∼200 seconds. All segments are finally averaged to compute the final PDS for each observing run. For very long exposures it was possible to use ∼200 seconds, as an adequate number of time frames were available to obtain an average PDS of good quality. In the case of shorter exposures, however, we had to use time frames of 50/100 seconds. The PDS are normalised according to Leahy et al. (1983) so that the Poisson noise level is about 2. We verified that at higher frequencies (> 100 Hz) the PDS are all flat and close to 2.

Once the PDS were calculated, we used XSPEC (Arnaud 1996) to model the data. We followed the method described in Ingram & Done (2012) to transform the PDS into the correct format needed for XSPEC, where we could easily fit the data with some functional models. In all PDS, we fitted the data using multiple components: a power law for the Poissonian noise at high frequencies, a zero-centred Lorentzian profile for band-limited noise, and a series of wide and narrow Lorentzian profiles for the different features in the PDS. In this way, we were able to evaluate the properties of the QPOs in the X-ray and optical bands.

From the PDS of the optical observations, we found QPOs in three observations: in April, June and October 2018 (Figs. 15 and 16 and Table 2). To exclude spurious signals due to systematic effects in the optical data, we tried to search for similar components in the PDS of the reference star and we did not find any hint of similar features. It is worth noting that the total fractional rms of the reference star during the April and June sessions was always substantially lower than the total fractional rms of MAXI J1820+070 (see Table A.1).

In the X-rays, the property of the QPOs have been already extensively analysed by Stiele & Kong (2020). A series of type-C LFQPOs are detected in the Swift/XRT and NICER data soon after the beginning of the first outburst (MJD 58198). These type-C QPOs have frequencies between 30 mHz and ∼1 − 2 Hz, and are seen until the source enters the IM state (from that moment on type-B QPOs appear in the PDS). Type-C QPOs reappeared again in the NICER data in October 18 after the new transition to the LH state. We did a similar analysis of the NICER data to compare the X-ray PDS to the optical one. We then focus our analyses in the periods close to the Iqueye/Aqueye+ observations (for example we did not study in detail the phase near the state transitions and during all the HS state, roughly between the beginning of July and September 18). In Fig. 13 we show the overall evolution of the central frequencies (ν₀) of the two most prominent QPOs found in the NICER data until the first state transition⁷. We also show the epochs of the different accretion states that the source entered and the epochs where it was possible to see the modulation of the light curve in the optical data (grey shaded area) and in the X-ray data (dashed orange area). An evolution of the central frequencies of the QPOs is clearly visible.

Fig. 13. Central frequency evolution of NICER QPOs (grey points) compared to the optical QPOs (green stars) and the periods where the optical (grey shaded area) and the X-ray (orange shaded area) modulation can be seen in the light curves. We also show the epochs of the different accretion states that the source entered, as reported by Shidatsu et al. (2019).

We divided the data into seven intervals, corresponding to different evolving characteristics of the QPOs central frequencies.

March to April 2018: In the optical, in MJD 58227, we observed two prominent QPOs (Zampieri et al. 2018a) on the top of three broad-band noise components (see middle-left panel of Fig. 15 and Table 2). Instead, in almost all the NICER data in the MJD range 58198-58236, we can see two QPOs (a main component and an upper harmonic; see Fig. 14). The central frequencies of these two components evolve exponentially with time with the same trend (ν₀ ∝ 10^{0.02 ⋅ t}). The main QPO in the optical data (upper green star) perfectly fits the expected frequency found from the X-ray trend. In the left panels of Fig. 15 we show a comparison of the X-ray/optical PDS considering the X-ray observations taken before and after the optical observation. We see a close match between the central frequencies of the main components, as well as the presence of a lower frequency component that is visible only in the optical data. The upper harmonic is clearly visible in the X-ray data.

Fig. 14. Time evolution of the central frequency of the two most prominent X-ray QPOs in NICER data (grey points) together with the central frequency of the optical QPOs (green stars). On the left we show the evolution during March-April 2018, while on the right is the evolution during June 2018.

Fig. 15. Comparison of optical and X-ray PDS of MAXIJ1820+070. The top and bottom panels show the NICER PDS, while the middle panels show the Iqueye PDS. The left panels correspond to the data taken around 18 April 2018, and those on the right correspond to the data taken around 8 June 2018. As no overlapping observations are available, we used for comparison the X-ray observations taken before and after the optical observations. The vertical red lines correspond to the central frequency of the QPOs observed in the optical PDS. In both periods the frequency of the most prominent QPO in the X-rays is in good agreement with the frequency of the second QPO in the optical band.

April to May 2018: The evolution of the central frequencies in this interval (MJD 58236-58249) is similar to the previous one, but with a slightly different dependence with the time (ν₀ ∝ 10^{0.025 ⋅ t}). No optical observations are present in this interval.

May 2018: During May 2018 (MJD 58249-58268), the frequencies of the QPOs stop increasing. In MJD 58259, we did observe the source with Iqueye, but we did not see any QPOs. This is most probably caused by the bad weather, since in the NICER observations before and after our observation the two QPOs are clearly visible. We estimated that a QPO with properties similar to the QPOs found in April or June had to have a fractional rms variability smaller than ∼2% (95% confidence level, Zampieri et al. 2018a) for the QPO to not be observable.

June 2018: In the optical, in MJD 58279, we observed again two prominent QPOs (Fiori et al. 2018) on the top of two broad-band noise component, plus a small third component that can be identified only when comparing the PDS of the optical and X-ray observations (see middle-right panel of Fig. 15 and Table 2). From the X-rays, in the interval MJD 58268-58288 we see that the frequencies decrease exponentially with time (Fig. 14, right panel) and again the trends for the main component and the upper-harmonic are the same (ν₀ ∝ 10^{−0.013 ⋅ t}). The main QPO in the Iqueye data and the upper-harmonic fit perfectly the expected frequencies found in the X-rays. Looking again at Fig. 15 we see a close correspondence between the central frequency of the main QPO in the X-ray and optical bands and a lower frequency component, which is again only visible in the optical data. The upper-harmonic, clearly visible in the X-ray data, is only marginally observed in the optical data.

June to July 2018: Before the state transition to the IM state (at ∼ MJD 58303.5) the evolution of the central frequencies of the QPOs changes again. At first the QPOs almost disappear, with only the upper-harmonic marginally visible in some of the observations (as reported by Stiele & Kong 2020). We decided not to include these very low-amplitude QPOs in Fig. 13. Around MJD 58297, the two harmonically related QPOs become well visible again with an exponentially increasing central frequency. Finally, when the source enters the IM state, a quickly evolving type-B QPO is observed in the NICER data (not shown in Fig. 13; more details for this phase can be found in Homan et al. 2020). No optical observations are present in this interval.

July to September 2018: After the state transition and during all the period when the source was in the HS state no QPOs are found. This is again in agreement with the optical data, where no QPOs are present. During July, when we also observed with Aqueye+ and the meteorological conditions were good, we estimated a fractional variability smaller than ∼0.5% (95% confidence level, Zampieri et al. 2018b) for a QPO not to be observable.

October to November 2018: When the source entered again the LH state, QPOs started to reappear in the NICER data and in the optical data. In the optical, in MJD 58404, we observed other three QPOs on top of broad-band noise components (see bottom panel of Fig. 16 and Table 2). In the X-rays, the QPOs can be seen only in a few observations (as again reported by Stiele & Kong 2020). We did not observe any clear evolving trend, in contrast to what was seen at earlier times.

Fig. 16. Comparison of optical and X-ray PDS of MAXIJ1820+070 for the data taken around 13 October 2018. The top panel shows the NICER PDS, and the bottom panel shows the Iqueye PDS. As no overlapping observations are available, we used for comparison the closest X-ray observation to our optical observation. The X-ray observation after our optical observation shows no clear QPO, and therefore we do not show the corresponding PDS. The vertical red lines correspond to the central frequency of the QPOs in the optical PDS. The frequency of the most prominent optical QPO (the second one in the optical PDS) is in good agreement with the frequency of the X-ray QPO.

4. Discussion

4.1. Optical super-orbital modulation

The values of the periodicities mentioned in this section are reported in Table 3. From the optical photometry, we calculated a period (p_{58283 − 58391}) that is about 0.5% longer than the orbital period (p_orb) found from optical spectroscopy in quiescence (Torres et al. 2019) and about 2% shorter than the one reported by Patterson et al. 2018 (p_P18). Later Patterson (2019) reported a slightly shorter period (p_P19) considering a different time window. These two estimates were calculated over an interval of about 15 and 30 days, which only partially overlaps with the interval used in our analysis. We repeated the analysis in similar intervals and we measured two periods that are within 1σ to those in Patterson et al. (2018) and Patterson (2019) (p_{58275 − 58290} and p_{58275 − 58310}). We also found that adding or removing dates at the beginning or at the end of this interval, greatly changes the inferred modulation period and leads to a higher actual uncertainty on this measurement. To understand these differences we have studied the behaviour of the inferred period using different starting and ending dates. This is shown in the upper panel of Fig. 17. The y-axis shows the values of the frequencies of the peak in the corresponding LS diagrams. We can see that, when we add more and more data, the frequency always evolves to a value close to the reference value of f_{58283 − 58391} = 1.4517 1/d. In addition, as the starting date is moved forward (going from MJD 58725 for the blue curve to MJD 58288 for the red curve), the different curves tend faster and faster towards the reference value. The blue curve, that is, the one including data at earlier dates (MJD < 58283), seems to deviate more from the reference value, both considering the inferred frequency using only the early data (MJD < 58309) and the whole time window (up to MJD 59391, zoomed-in panel in Fig. 17). Looking at the blue curve, it is clear that the period computed using MJD 58275-58290 and MJD 58275-58310 is longer than the one inferred from all the data between MJD 58283-58391. This may be caused by the fast evolution of this signal in the first few days after its appearance.

To check the validity of our assumption, we tried to fit the frequency and the phase of the modulation in different time intervals, following what was done by Thomas et al. (2022b). We first fitted the frequency and the phase dividing the data in two intervals, before and after MJD 58309. Considering the data in MJD 58283-58309, we found a frequency of 1.4488 ± 0.0005 1/d and a phase of 0.72π ± 0.02 (consistent within 1σ with the value reported in Patterson 2019), while considering the data in MJD 58309-58391 we found a frequency of 1.4514 ± 0.0003 and a phase of 0.72π ± 0.03 (consistent within 1σ with the frequency found in the entire time window). The phase of the signal, varying the frequency in the two time windows, remains constant.

We then fit the phase of the signal fixing the frequency at specific values: the frequency f_{58283 − 58391} and the orbital frequency f_orb reported by Torres et al. (2019). The phases are fitted considering the same time windows of Fig. 7 and are plotted in the bottom panel of Fig. 17. Using f_{58283 − 58391} (blue points) we found a roughly constant phase in all the windows – only the data in the first window are slightly out of phase – and a good agreement with the overall phase (ϕ_{58283 − 58391} = 0.67π ± 0.03). On the other hand, using f_orb, we found that the phase monotonically decreases with time, starting from ∼1.3π down to ∼ − 0.3π. This can be easily interpreted as being caused by an incorrect value of the frequency, so that in order to compensate for it, the phase must evolve over time.

Fig. 17. Top: Frequency of the highest peak in the LS diagrams computed considering different starting dates. The x-axis correspond to the stopping date used for computing the LS diagrams. The frequency always quickly evolves to the reference value of ∼1.4517 1/d−1 (dashed black line). In addition, as the starting date is moved forward, the curves tend faster and faster towards the reference value. The curve that deviates most from this trend is the one corresponding to the starting date MJD 58275 (blue curve), which includes data when the modulation was still evolving. Bottom: Phase evolution considering the dates shown in Fig. 7 (excluding the data in MJD 58326-58338) and using two different frequencies. In blue the frequency computed in this work (f58283 − 58391 = 1.4517 1/d) and in orange the frequency related to the orbital motion (forb = 1.4588 1/d, Torres et al. 2019). The dotted lines are cubic splines fitted to the data.

From all these findings we therefore tentatively conclude that the modulation is not caused by the orbital motion (as suggested by Thomas et al. 2022b for all the data after MJD 58309) and that it probably evolved from the frequencies reported in Patterson et al. (2018) and Patterson (2019) to the frequency f_{58283 − 58391}. A similar result is also reported in Niijima et al. (2021; see the top panel of their Fig. 2), who found a period within 2σ to the one reported here (p_N21; Table 3) and interpreted it as the true period of the super-orbital modulation.

4.2. X-ray/optical super-orbital modulation

A clue that the optical modulation visible during MJD 58283-58391 and the X-ray modulation visible during MJD 58309-58326 are driven by the same physical mechanism is given by the behaviour exhibited during the different state transitions. The optical modulation started to show up around the beginning of the first optical rebrightening, as can be seen in Fig. 1. At that time MAXI J1820+070 was still in the LH state. However, around the same dates, the evolution of the QPOs central frequency in the X-ray band changed (Fig. 13). The optical modulation then lasted all along the HS state and stopped only a few days after the new transition to the LH state. However, the modulation of the X-ray light curve appeared only immediately after the state transition from the IM to the HS states and lasted only for a few days. Interestingly, the optical modulation immediately after MJD 58326 (when the X-ray modulation is no longer visible) also seems to vanish for a few days (bottom left panel of Fig. 7) before appearing again, but with a much smaller amplitude. Indeed, this could be the proof that the process that generates the optical modulation could be the same that produces the modulation that seems visible in the X-rays between MJD 58309 and MJD 58326 (Fig. 9).

As already reported, if the X-ray modulation between MJD 58309 and 58326 is real, the optical and X-ray are out of phase by about ∼1.1π, with the X-ray leading the optical. Given the large uncertainty in the X-ray modulation and that in a small time window, we can also fit the optical data using the orbital modulation and compare the phases in the two bands using the orbital solution by Torres et al. (2019) (using the orbital frequency for both bands or only for the X-ray band). Even with these assumptions the X-ray always leads the optical, again by about ∼1.1π. This evidence appears to be consistent with a common geometrical origin of the optical and X-ray modulations (e.g. a warped disc).

4.3. X-ray/optical LFQPOs

The type-C LFQPOs observed in the PDS are synchronous over at least 5 orders of magnitude in terms of energy, from the optical band up to the hard X-ray band (Zampieri et al. 2019; Ma et al. 2021; Mao et al. 2022; Thomas et al. 2022a). A close link between the optical and X-ray bands is also evident from the PDS obtained in this work (see Figs. 15 and 16), even if no strictly simultaneous X-ray/optical observations were carried out. The frequencies of the LFQPOs in the two bands are consistent once the expected evolution of the central frequency is taken into account (Figs. 15, 16, 13, and 14). This result is confirmed by other authors who analysed simultaneous optical/X-ray data (Paice et al. 2019, 2021; Thomas et al. 2022a). Remarkably, the optical QPOs are present in all our observations when the source was in the LH state. In two observations out of three it was possible to detect a lower frequency component in the optical PDS (in October 2018 the lower frequency component is probably present, but seems to be very broad), while in the X-rays this component is not observed. The frequencies of the lower component in the optical data and of the most prominent peak in the X-rays appear to be harmonically related, with a ratio of ∼1:2. On the other hand, a higher frequency component is almost always present in the X-ray PDS, while in optical data it is always less pronounced.

Optical, ultraviolet and infrared QPOs have been observed in other black hole binaries in the past (Motch et al. 1983; Imamura et al. 1990; Hynes et al. 2003; Durant et al. 2009; Gandhi et al. 2010; Veledina et al. 2015; Kalamkar et al. 2016; Vincentelli et al. 2019, 2021). In several cases the central frequency of the QPOs is the same in the different bands (Hynes et al. 2003; Veledina et al. 2015; Vincentelli et al. 2021). In other cases there was no clear identification of an X-ray counterpart, either because no simultaneous X-ray observations were available (Imamura et al. 1990) or because it was just not detectable (Durant et al. 2009; Gandhi et al. 2010; Vincentelli et al. 2019). Interestingly, in GX 339-4, the central frequency of the optical/infrared QPO is half of the frequency of the QPO in the X-ray band (Motch et al. 1983; Kalamkar et al. 2016). Usually the fundamental QPO is considered to be the one with the highest peak amplitude in the X-ray band (Ingram & Motta 2019) based on the fact that in many periodic processes the fundamental frequency has higher amplitude than its overtones. However, Veledina et al. (2013) predicted that if the inclination of the system is between 45° and 75°, the second harmonic could be stronger than the fundamental. Therefore, depending on the adopted model, the first harmonic of the QPO could be the one observed in the optical/infrared/ultraviolet energy range.

Regardless of whether or not the fundamental frequency is that observed in the optical band, how can we explain the close link between X-ray and optical emission? We have already seen that there are evidences that the regions where the photons in the two bands are emitted must be very close (Paice et al. 2019, 2021). However, the geometry of the system is yet to be fully understood, and therefore different models can be adopted to interpret the data. Assuming a lamp-post geometry for the hot corona in the innermost regions of the accretion flow, measurements of the reverberation time lags between the X-ray continuum emission and the (broad) iron K emission line (produced in the outer irradiated accretion disc) imply the existence of an accretion disc with a characteristic size 6 − 20 times smaller (15 gravitational radii, R_{_g}) than what previously seen in other similar objects (Kara et al. 2019). Given that the reverberation time lag frequency increases with time, the hot corona should also contract and decrease its height following the evolution of the outburst (Kara et al. 2019; Buisson et al. 2019). A different interpretation, however, is possible; namely, one can consider a standard disc model (Shakura & Sunyaev 1973) where the steady increase of the reverberation lag frequency is caused by the decrease of the internal truncation radius r_{_in} (De Marco et al. 2021; Kawamura et al. 2022). This interpretation is also supported by broadband X-ray spectra that show that r_{_in} should be located close to the central black hole at only a few tens of gravitational radii (Marino et al. 2021).

We can thus try to model the properties of the observed LFQPOs by considering a Lense-Thirring (LT) precession of the hot accretion flow near the black hole (Stella & Vietri 1998; Motta et al. 2014). The precession could occur either in the hot inner flow (Ingram et al. 2009; Veledina et al. 2013) or at the base of the jet (Markoff et al. 2005; Stevens & Uttley 2016). In fact, the observed LFQPOs may originate from the precession of a small-scale jet, as suggested by the gradual increase of the soft phase lag with increasing energy, up to and above 200 keV (Ma et al. 2021). The evolution of the frequency of the LFQPOs (Fig. 13) implies a variable characteristic radius and, since the frequency is small (< 1 Hz), this radius can be relatively large (larger than the innermost stable circular orbit). In the LT precession model, the precession frequency (ν_{_LT}) is generally linked to the frequency of the fundamental QPO. Assuming that the fundamental frequency corresponds to the centroid frequency of the QPO in the optical (ν_{_opt}), at half of the frequency of the most prominent peak in the X-ray data (ν_{_X}), implies that the precession frequency is half of what is inferred from the X-rays. The weak/lack of detection of the first harmonic in the X-rays may be caused by, for example, partial obscuration/eclipse of the X-ray emitting region at certain precession phases. This is in agreement with the prediction of Veledina et al. (2013) given that the inclination of the system is estimated to be ∼63° (Atri et al. 2020). Thus, considering that the range of the most prominent QPO in the X-ray band is between ∼30 mHz and ∼1570 mHz (maximum frequency reached before the transition to type-B QPO; Stiele & Kong 2020), the precession frequency should then be between ν_{_LT,min} ≃ 15 mHz and ν_{_LT,max} ≃ 785 mHz.

From the expression of the LT precession frequency of a particle in the Kerr metric (e.g. Merloni et al. 1999; Ingram et al. 2009), we can estimate the dimensionless spin of the central black hole, a_*, with:

$$\begin{array}{c}\hfill {\nu }_{\mathrm{LT}}={\nu }_{\varphi }[1-\sqrt{1-\frac{4a_{\ast }}{r^{3/2}}+\frac{3a_{\ast }^2}{r^2}}],\end{array}$$(2)

where ν_ϕ is the frequency of a circular orbit in the equatorial plane at a given radius r (in units of the gravitational radius R_{_g} = GM/c²). Here, the radius is the truncation radius of the standard disc r_{_in}. In the following, we assume that it varies in the range 10 − 40 R_{_g}, accordingly to the findings of De Marco et al. (2021) and Kawamura et al. (2022). If we assume that the minimum and maximum precession frequencies are linked to the maximum and minimum radii at which the standard disc is truncated at a certain time, from a comparison of Eq. (2) with the observed frequencies, we derive an estimate for a_*⁸. The results are shown in Fig. 18, where we considered two different scenarios: the scenario in which the true fundamental frequency is the one visible in the optical data (left panel) and the scenario in which the true fundamental frequency is the most prominent peak seen in the X-rays (right panel). The two panels show the LT precession frequency as a function of the truncation radius for different values of the black hole spin (coloured lines). The vertical lines correspond to the estimated range of radii for the truncation radius, while the horizontal lines correspond to the minimum and maximum LT precession frequency that we considered. The grey area represents the parameter space that is in agreement with the data. Taking into account the uncertainties on r_{_in}, if the fundamental frequency is the one observed in the optical band, the allowed dimensionless spin of the black hole should be ≲0.15, meaning a relatively slowly spinning black hole. This value is in good agreement with the spin estimated by Zhang et al. (2020). Conversely, if the fundamental frequency is the one observed in the X-ray band, the allowed a_* should be ≲0.3, meaning a black hole that can spin twice faster then in the previous scenario. It is worth emphasising that any estimate of the black hole spin depends strongly on both the adopted model and the choice of the inner disc truncation radius. Clearly, different assumptions may lead to different results (see e.g. Bhargava et al. 2021). In particular, although the LT precession model provides a reasonable physical framework for producing the observed QPOs frequencies, different interpretations of the reverberation frequency lag and/or the geometry of the corona are possible (Kara et al. 2019; Buisson et al. 2019; De Marco et al. 2021; Kawamura et al. 2022), making the estimate of the truncation radius somewhat uncertain. However, as mentioned above, selecting r_{_in} in the range inferred from some recent X-ray spectral analyses (De Marco et al. 2021; Kawamura et al. 2022) we obtain consistent results.

Fig. 18. Lense-Thirring precession frequency computed following Eq. (1) in Ingram et al. (2009) and as a function of the truncation radius r_in for different values of the black hole spin a* (coloured lines). The vertical lines correspond to the estimated range of radii for the truncation radius, while the horizontal lines correspond to the minimum and maximum LT precession frequency that we considered. The shaded coloured areas around the lines are the uncertainties from estimated mass of the black hole as reported by Torres et al. (2020), while the grey area represents the parameter space that is in agreement with the data. On the left, we show the solutions that correspond to the scenario in which the fundamental frequency is traced by the QPOs at lower frequencies seen in the optical data. On the right, we show the solutions that correspond to the scenario in which the fundamental frequency is traced by the most prominent peak in the X-rays.

We note that a small value of the inner disc radius (∼10 R_{_g}) is also consistent with the onset of the super-orbital modulation at the end of the LH state, because the stronger X-ray irradiation from the inner disc can more easily trigger the warp observed in the X-rays and in optical.

4.4. Optical spectroscopy

From the optical spectroscopic data we obtain additional interesting information on the system during the first phase of the outburst. The blue spectra up until MJD ∼58380 are indicative of an increased X-ray irradiation of the outer disc. Indeed, until MAXI J1820+070 returns to the LH state, the X-ray flux is high (see e.g. Fig. 1 in Prabhakar et al. 2022) and irradiation is significant. Furthermore, when the X-ray flux is maximum, the ambient density is probably high enough to lead to the activation of the Bowen mechanism, as observed in the broad blend centred at 4640 Å. Finally, the variability of the profile of some broad lines (e.g. HeII 4686 Å) may be a further indicator of a (possibly warped) precessing disc. When the super-orbital modulation starts to be detected, the equivalent width of the spectral lines in Fig. 12 tends to increase, while it decreases when the modulation vanishes. Therefore, these spectral features seem to disappear when the super-orbital modulation stops and the first phase of the outburst terminates. This is a further indication that there exists a close relationship between the emission in the X-ray and optical bands.

5. Conclusions

In this work, we have presented a detailed analysis of data obtained from an intensive multi-timescale variability observing campaign of the new black hole X-ray binary MAXI J1820+070. We made use of low-cadence and fast optical photometry data as well as NICER X-ray data. The main part of our analysis was focused on the initial phases of the outburst from March to October 2018.

We detected an optical super-orbital modulation that started around the beginning of the first optical rebrightening with a frequency of f = 1.4517 ± 0.0001 1/d and an evolving profile. After the transition from the LH to the HS state, a hint of a modulation appeared to emerge for a few days also in the X-ray band. We found evidence of an X-ray super-orbital modulation with a frequency in close agreement with the optical one and out of phase by about ∼1.1π (X-ray leading the optical). The optical super-orbital modulation lasted until the transition from the HS to the LH state.

In three epochs (April, June, and October 2018), we observed synchronous optical and X-ray type-C LFQPOs. A lower frequency component harmonically related to the most prominent X-ray QPO (with a ratio 1:2) is visible in two optical observations. If the lowest modulation frequency is that observed in the optical power density spectrum, the characteristic variability frequency of MAXI J1820+070 is lower than what is inferred from the ‘fundamental’ QPO in the X-rays. Considering the LT precession model, under the assumption that the truncation radius of the standard disc varies in a range r_{_in} ≃ 10 − 40 R_{_g} and taking as precession frequency the lowest frequency traced by the optical observations, we estimated a relatively slowly spinning black hole with a_* ≲ 0.15, which is in close agreement with other measurements in the literature. In these assumptions, the observed optical and X-ray large-scale super-orbital warp could be triggered by the self-irradiation of the outer disc induced by the inner disc when extending inward (r_{_in} ∼ 10 R_{_g}) before the source enters the HS state. Further observations of black hole binaries in outburst at different wavelengths and on multiple timescales can help explain how synchronous QPOs are produced over such a large energy range (from the optical to the hard X-rays) and whether the super-orbital modulation is directly linked to the process that generates the QPOs.

Acknowledgments

This research made use of the following PYTHON packages: MATPLOTLIB (Hunter 2007), NUMPY (van der Walt et al. 2011), SCIPY (Virtanen et al. 2020), PANDAS (McKinney 2010), ASTROPY (Astropy Collaboration 2013, 2018, 2022) and EMCEE (Foreman-Mackey et al. 2013). We acknowledge financial support from INAF Research Grant “Uncovering the optical beat of the fastest magnetized neutron stars (FANS)”. (S.S.) This work was also supported by grants from the Slovak Academy VEGA 2/0030/21 and APVV-20-0148.

References

Appendix A: Log of the observations

In this appendix, we report the summaries of the optical observations. We do not report the summary of the NICER X-ray observations since it can be retrieved from the official HEASARC archive.⁹

All Tables

All Figures

Fig. 1. Light curves of MAXI J1820+070 in different optical bands since the epoch of the first maximum after the discovery of the source and until August 2020. The data taken in the Landolt UBVRI bands (Landolt 1992) are shown in violet, blue, light green, light red, and light orange, respectively. The darker green, red, and orange represents the g′, r′, and i′ Sloan bands of the APASS system (Henden & Munari 2014; Henden et al. 2018). The dashed black vertical lines indicate approximately the dates when the local maximum of the optical luminosity is reached during all the subsequent bursts. The red vertical lines indicate the epochs when we observed optical LFQPOs with Iqueye. The yellow vertical lines indicate the dates when Torres et al. (2019) measured the periodicity of the binary system (1.4591 1/d) spectroscopically. The grey shaded area indicates the interval in which a photometric periodicity is detected in this dataset, while the orange shaded area indicates the interval in which a periodicity can be seen also in the X-ray data.

In the text

Fig. 2. Power of the highest peak in the LS diagram of the V-band light curve after detrending the data and considering different starting and ending dates. This plot was computed to understand which temporal interval to use for searching the best timing solution of the system by means of fitting the sinusoid from Eq. (1) to the data. The power increases from blue to red. A vertical straight line marks the time (MJD 58283) at which the power starts to rise, while a horizontal straight line marks the time (MJD 58391) at which the power starts to decrease.

In the text

	Fig. 3. Frequency (black line) and power (background colour chart) of the highest peak in the LS diagram computed after detrending the V-band data using the average over N days (x-axis). The maximum power is reached for N = 3 days. The spread of the frequency is compatible with the uncertainties of the measurement.
In the text

	Fig. 4. Detrended optical data in the time interval selected to measure the periodicity. The colour code is the same as in Fig. 1.
In the text

	Fig. 5. Comparison of LS diagrams in different time intervals. The LS diagram computed considering all the available V-band data is shown in orange, whereas the diagram computed considering only the data within (outside) the interval MJD 58283-58391 is shown in blue (green).
In the text

Fig. 6. Data from the V-band between MJD 58283 and MJD 58391 folded with the value of the frequency f reported in Table 1. The dotted black line is the sinusoidal fitted model with the corresponding 3σ confidence uncertainty shown as a grey shaded area. Some residuals can be seen in the plot and are caused by the data within certain time intervals where the oscillation is not perfectly visible or slightly out of phase, as can be seen from Fig. 7.

In the text

Fig. 7. Evolution of the amplitude and shape of the periodical signal in the light curve of MAXI J1820+070 between MJD 58283 and MJD 58391 seen in the folded optical data in different sub-intervals. The colour code is the same as in Figs. 1 and 4. The dashed black line corresponds to the best-fitting sinusoid (Eq. (1)) for the entire interval. Except for the first panel, the phase of the signal remains constant, while its shape evolves with time. All the optical bands are in good agreement.

In the text

	Fig. 8. Lomb-Scargle diagram (black curve) computed with the detrended X-ray data in the interval MJD 58309-58326. The red curve correspond to the 0.01% false alarm probability level (see text for details). The highest peak is quite broad, but is centered at the frequency inferred from the optical photometry and well above the false alarm probability level in the range of frequencies around it.
In the text

Fig. 9. Comparison of optical and X-ray folded data in the interval MJD 58309-58326 computed using the best-fitting frequency obtained with all optical data between MJD 58283 and MJD 58391 as reported in Table 1. Data are shown together with the best-fitting sinusoid (dashed black lines). The data in the upper panel correspond to the optical data (coloured dots). Lower panel: X-ray data (grey dots) as well as the same data binned in phase (black dots with error bars, 15 bins per phase). The error bars are computed from the standard deviation of the data falling inside every phase bin.

In the text

Fig. 10. Ten spectra of MAXI J1820+070 recorded with the Asiago 1.22m + B&C telescope. From top to bottom, the spectra are arranged in order of observing date. They have been offset for plot clarity by the indicated amount. Next to the offset is listed the B-band magnitude exhibited by MAXI J1820+070 at the time of the spectroscopic observation. The brightest emission lines are identified as well as a sample of the interstellar features seen in absorption (CaII, NaI, and the diffuse interstellar band at 6281 Å as representative of several other well visible diffuse interstellar bands).

In the text

	Fig. 11. Zoomed-in view from Fig. 10 of the MAXI J1820+070 spectra for 2018 covering hydrogen Hβ, HeII 4686 Å, and the broad 4640 Å blend attributed to NIII lines. The spectra are offset by the indicated quantity for clarity of the plot. The typical velocity separation of 850 km s−1 for the double-peaked profiles is marked.
In the text

Fig. 12. Bottom panel: Evolution of the equivalent width (in Å) of representative emission lines as a function of time, calculated from the spectra of MAXI J1820+070 in Fig. 11. Top panel: Ratio between the equivalent widths of two of these lines, HeII 4686 Å and Hβ. The colour bar on the top indicates the value of the B band magnitude of the source at the time of the observations. We show for comparison also the periods where the optical (grey shaded area) and the X-ray (orange shaded area) modulation are seen in the light curves, together with the epochs of the different accretion states that the source entered (Shidatsu et al. 2019).

In the text

	Fig. 13. Central frequency evolution of NICER QPOs (grey points) compared to the optical QPOs (green stars) and the periods where the optical (grey shaded area) and the X-ray (orange shaded area) modulation can be seen in the light curves. We also show the epochs of the different accretion states that the source entered, as reported by Shidatsu et al. (2019).
In the text

	Fig. 14. Time evolution of the central frequency of the two most prominent X-ray QPOs in NICER data (grey points) together with the central frequency of the optical QPOs (green stars). On the left we show the evolution during March-April 2018, while on the right is the evolution during June 2018.
In the text

Fig. 15. Comparison of optical and X-ray PDS of MAXIJ1820+070. The top and bottom panels show the NICER PDS, while the middle panels show the Iqueye PDS. The left panels correspond to the data taken around 18 April 2018, and those on the right correspond to the data taken around 8 June 2018. As no overlapping observations are available, we used for comparison the X-ray observations taken before and after the optical observations. The vertical red lines correspond to the central frequency of the QPOs observed in the optical PDS. In both periods the frequency of the most prominent QPO in the X-rays is in good agreement with the frequency of the second QPO in the optical band.

In the text

Fig. 16. Comparison of optical and X-ray PDS of MAXIJ1820+070 for the data taken around 13 October 2018. The top panel shows the NICER PDS, and the bottom panel shows the Iqueye PDS. As no overlapping observations are available, we used for comparison the closest X-ray observation to our optical observation. The X-ray observation after our optical observation shows no clear QPO, and therefore we do not show the corresponding PDS. The vertical red lines correspond to the central frequency of the QPOs in the optical PDS. The frequency of the most prominent optical QPO (the second one in the optical PDS) is in good agreement with the frequency of the X-ray QPO.

In the text

Fig. 17. Top: Frequency of the highest peak in the LS diagrams computed considering different starting dates. The x-axis correspond to the stopping date used for computing the LS diagrams. The frequency always quickly evolves to the reference value of ∼1.4517 1/d−1 (dashed black line). In addition, as the starting date is moved forward, the curves tend faster and faster towards the reference value. The curve that deviates most from this trend is the one corresponding to the starting date MJD 58275 (blue curve), which includes data when the modulation was still evolving. Bottom: Phase evolution considering the dates shown in Fig. 7 (excluding the data in MJD 58326-58338) and using two different frequencies. In blue the frequency computed in this work (f58283 − 58391 = 1.4517 1/d) and in orange the frequency related to the orbital motion (forb = 1.4588 1/d, Torres et al. 2019). The dotted lines are cubic splines fitted to the data.

In the text

Fig. 18. Lense-Thirring precession frequency computed following Eq. (1) in Ingram et al. (2009) and as a function of the truncation radius r_in for different values of the black hole spin a* (coloured lines). The vertical lines correspond to the estimated range of radii for the truncation radius, while the horizontal lines correspond to the minimum and maximum LT precession frequency that we considered. The shaded coloured areas around the lines are the uncertainties from estimated mass of the black hole as reported by Torres et al. (2020), while the grey area represents the parameter space that is in agreement with the data. On the left, we show the solutions that correspond to the scenario in which the fundamental frequency is traced by the QPOs at lower frequencies seen in the optical data. On the right, we show the solutions that correspond to the scenario in which the fundamental frequency is traced by the most prominent peak in the X-rays.

In the text