Owens Valley Radio Observatory monitoring of LS I +61 ◦ 303 completes three cycles of the super-orbital modulation

Context. The high-mass X-ray binary LS I + 61 ◦ 303 is composed of a Be-type star and a compact object in an eccentric orbit. The emission from this source is variable and periodic across the electromagnetic spectrum, from radio to very high-energy gamma rays. The orbital period has been determined as P 1 ≈ 26 . 5 d, and the source also features a super-orbital period with a value of P long ≈ 4 . 6 years. Long-term monitoring of the binary by the Owens Valley Radio Observatory (OVRO) at 15 GHz has now completed 13.8 years, which corresponds to three full cycles of the super-orbital period. This is exactly one long-term cycle more than in the previous publication about OVRO observations of this source. Aims. Our aim is to investigate the presence and the stability of periodic signals in the radio data and to test if they are in agreement with previous results. This will contribute to the understanding of the physical processes behind the non-thermal emission from this source. Methods. We performed a timing analysis of the OVRO radio light curve and made use of the generalized Lomb-Scargle periodogram. We also combined the OVRO data with the full archive of previous radio observations and computed the discrete autocorrelation function. Results. The most powerful features in the periodogram of the OVRO data are two peaks at P 1 = 26 . 49 ± 0 . 05d and P 2 = 26 . 93 ± 0 . 05d, which are well separated from each other and clearly stand out above the very low noise level. The previously detected long-term period is still present in these new radio data, and our measurement is P long = 1698 ± 196d. Dividing the OVRO data into three segments of equal length showed that the two periods, P 1 and P 2, are present in the periodogram of each of the consecutive long-term cycles. Our analysis of the full radio archive resulted in the detection of the same three periods, and the autocorrelation function showed a regular pattern, proving the continuity of the decades-spanning stability of the super-orbital modulation. In addition, we report a possible systematic modulation of the radio flux density with a timescale of approximately 40 years that has so far remained unnoticed. Conclusions. The physical model of a relativistic jet whose mass loading is modulated with the orbital period P 1 and is precessing with the slightly larger period P 2 , giving rise to a beating with period P long , had previously been able to reproduce the radio and gigaelectron volt emission from this source. The ongoing presence and the stability of the periodic signals imply that this model is still the most plausible explanation for the physical processes at work in this source.


Introduction
An X-ray binary is a stellar system composed of a normal star and a compact object, which can be either a neutron star or a black hole.Some X-ray binary systems are particularly bright in the γ ray regime (Dubus 2013), and where and how this high-energy emission is produced is still a debate.The high-mass X-ray binary LS I +61 • 303 is composed of a Be-type star and a compact object (Casares et al. 2005).The system is a source of emission across the electromagnetic spectrum, from radio to the very high-energy γ-rays (Jaron 2021, and references therein), and is also included in the search for neutrino emission (Abbasi et al. 2022, and references therein).
The timing characteristics of its electromagnetic emission make LS I +61 • 303 special among the class of γ-ray emitting stellar binary systems.Radio emission was first detected from this source during a search for variable radio sources along the galactic plane (Gregory & Taylor 1978).Since then, it has been the target of repeated radio monitoring, which has revealed that it is not only variable but also periodic on different timescales.The most accurate measurement of the orbital period was obtained by Gregory (2002), who analyzed radio data from a long-term monitoring program of the Green Bank Interferometer (GBI) at 2 and 8 GHz and reported a value of P 1 = 26.4960± 0.0028 d.In the same article, Gregory (2002) also firmly established that the radio emission from this source is subject to a super-orbital long-term modulation with a period of P long = 1667 ± 8 d.This super-orbital modulation is present across the spectrum until the teraelectron volt regime, with a systematic phase relationship between wavelengths (Jaron 2021).In particular, the large database in the radio has proven that this long-term modulation has remained stable over decades of observations (Massi & Torricelli-Ciamponi 2016).
Concerning the eccentricity of the orbit, analyses of optical data from LS I +61 • 303 have produced mixed results.Casares et al. (2005) obtained a value of e = 0.72 ± 0.15 using absorption lines that were least contaminated by emission lines from the Be circumstellar disk.Aragona et al. (2009) analyzed the He I λ6678 absorption line, which resulted in e = 0.537 ± 0.034.However, this absorption line was slightly affected by blue and red wings in emission.More recently, Kravtsov et al. (2020) reported a value of only e ≈ 0.1, which they determined by model fitting the optical linear polarization of the source.However, their reported value is at odds with other observations.This outcome can be understood considering the approach used by Kravtsov et al. (2020), as they used a simple model of variations induced by the orbital motion of the compact star with respect to the decretion disk of the Be star, and they did not include the other disk in their model, i.e., the disk around the compact object.For these reasons, we assumed a high eccentricity of the orbit for the interpretation of our results and considerations about the physical processes at work in this source.
Evidence for periodic accretion and ejection along the orbit of LS I +61 • 303 has been given in Massi et al. (2020), and the associated X-ray luminosity is related to its photon index, as established for accreting black holes at all masses (Massi et al. 2017).Evidence for the presence of a jet in LS I +61 • 303 comes from direct very long baseline interferometry (VLBI) observations (Massi et al. 2012;Wu et al. 2018).Indirect evidence has been obtained by analyzing the radio spectrum.A six-year archive of radio observations was analyzed by Massi & Kaufman Bernadó (2009), revealing the characteristics of radio jet emission, namely, a radio peak with a flat spectrum followed by an optically thin radio outburst, which is typically observed in microquasars (i.e., in accreting black holes or neutron stars with a low magnetic field) (Mirabel & Rodríguez 1999).This is also in agreement with the expected two-peaked Bondi & Hoyle (1944) accretion rate profile along the eccentric orbit of this source, as outlined by several authors (Taylor et al. 1992;Marti & Paredes 1995;Bosch-Ramon et al. 2006;Romero et al. 2007).A physical model of a self-absorbed jet, which through precession gives rise to periodic changes in the Doppler boosting of the intrinsic emission reproduces 37 years of radio observations and its spectral index (Massi & Torricelli-Ciamponi 2014).The inclusion of external inverse Compton scattering and synchrotron self-Compton into that model reproduces several years of simultaneous radio and Fermi-LAT γ-ray observations and their timing characteristics (Jaron et al. 2016).This model also explains why the radio emission peaks only once per orbit (i.e., at apastron), the reason being that at periastron the relativistic electrons of the jet suffer catastrophic inverse Compton losses in the strong UV photon field in the proximity of the Be star.Finally, the presence of the aforementioned super-orbital modulation at all observed wavelengths and a systematic phase offset between wavelengths finds a straightforward explanation in a precessing jet in which the higher energy emission is produced upstream from the lower energy emission (Jaron 2021).Nevertheless, there is also a pulsar scenario being discussed for LS I +61 • 303 as a generalization of the phenomena observed in the pulsar binary PSR B1259−63 to the larger class of γ ray binaries (Dubus 2013).
The possibility that the physical processes behind the nonthermal emission from LS I +61 • 303 are powered by the spindown of a millisecond pulsar have been discussed since soon after the discovery of the source (Maraschi & Treves 1981).However, despite dedicated searches, pulses were never detected from this source (Cañellas et al. 2012).The only signal that pointed to the presence of a pulsar in the system was a short X-ray burst that was interpreted as a "magnetar-like event" by Torres et al. (2012).This signal was observed from a region on the sky that contains several sources other than LS I +61 • 303, putting into question the association with that particular source.More recent observations with the Five-hundred-meter Aperture Spherical radio Telescope (FAST) have resulted in the detection of 42 radio pulses from the direction of LS I +61 • 303, with a pulse period of P = 269.15508± 0.00016 ms (Weng et al. 2022).The field of view of these observations is even larger than that of Torres et al. (2012), so the possibility of the pulses originating from another source cannot be ruled out.Furthermore, the pulses were detected during an orbital phase of Φ = 0.59 with an exposure time of three hours (see Table 1 in Weng et al. 2022).Two subsequent observations listed in their Table 1 made with the same telescope during similar orbital phases (Φ = 0.58, 0.62) and with similar exposure times (3 and 2 h, respectively) did not result in the detection of any pulses.This means that, so far, there is not any proof of a dependency of the occurrence of pulses on the orbital period of the system.Finally, Weng et al. (2022) did not detect any Doppler shift of the pulses, which would indicate that they do not originate from a pulsar in a binary system at all.In conclusion, there is not enough evidence to make a firm association between the observed pulses and the astrophysical object LS I +61 • 303.In their supplementary Fig. 1, Weng et al. (2022) reported a spin-down of Ṗ = (4.2± 1.2) × 10 −10 that should be regarded as tentative because it does not result from direct timing, as also pointed out by Suvorov & Glampedakis (2022).Assuming this value would imply a spin-down power of Ė = (8.5 ± 2.4) × 10 38 erg s −1 , which is indeed a relatively large value, but it is not reliable because of the questionable measurement of Ṗ itself.More observations would be needed to measure the spin-down with a higher accuracy and to probably disentangle it from the orbital motion if the detected pulsar is indeed part of a binary system.
Even if LS I +61 • 303 contained a pulsar, this would not rule out the possibility of a jet.Pulsar wind nebulae can have flat radio spectra and jets (see, e.g., Slane 2017 for a review).Extended X-ray emission close to the γ-ray-loud binary system LS 5039 has been interpreted as a pulsar wind nebula by Durant et al. (2011).For LS I +61 • 303, however, such an observation has not been reported anywhere as of yet.Only the extended structure in the aforementioned VLBI images (Massi et al. 2012) has previously been interpreted as a "cometary tail" by Dhawan et al. (2006).The reported values of P and Ṗ by Weng et al. (2022) implying a surface magnetic field on the order of 10 14 G would certainly inhibit the formation of a radio jet (Massi & Kaufman Bernadó 2008 determined an upper limit of 10 7 −10 8 G for the formation of a radio jet from a neutron star).Such a large magnetic field is therefore in contradiction to the observation of a radio jet that has a self-absorbed spectrum (Massi & Kaufman Bernadó 2009;Zimmermann et al. 2015), which is the spectrum that is typical of microquasars.Observations of radio emission from the millisecond pulsar binary PSR J1023+0038, which has a magnetic field of ∼10 8 G (Deller et al. 2012), has recently been provided by Baglio et al. (2023) and interpreted as originating from a radio jet.
The precise timing characteristics of LS I +61 • 303 make this source a unique laboratory to study the physical processes behind the non-thermal emission.In particular, the presence of the super-orbital modulation sets this source apart from other objects in the class of γ-ray emitting X-ray binaries.In this article, we present new radio observations of LS I +61 15 GHz (Richards et al. 2011).With a baseline of 13.8 years, these observations have now completed three cycles of the 4.6-year super-orbital modulation.This represents exactly one super-orbital cycle more of data than was available at the time of the previous publication about OVRO monitoring of this source presented in Jaron et al. (2018).The first aim of our analysis in this work is to investigate the presence of periodic signals in the updated OVRO data.The second aim is to test the stability of periodic signals, and the long-term modulation in particular, by combining the new OVRO observations with the radio archive from Massi & Torricelli-Ciamponi (2016).
The paper is structured as follows.In Sect.2, we present the data sets.In Sect.3, we describe the methods that we used for data analysis.We present our results in Sect. 4 and discuss them in Sect. 5. We give our conclusions in Sect.6.

OVRO monitoring
The OVRO regularly observes LS I +61 • 303 at 15 GHz as part of a monitoring program (Richards et al. 2011).The radio light curve that we used for our analysis covers the time span modified Julian day (MJD) 54908−59961 (2009March 18-2023 January 17).This represents 13.8 years of data, which is three cycles of the long-term modulation of LS I +61 • 303.The exact super-orbital phase range of the data is Θ = 6.92−9.96(as defined in Sect.3.3).A description of the data calibration can be found in Sect.2.1 of Jaron et al. (2018), and for further details, we refer to Richards et al. (2011).
Figure 1 shows the OVRO flux density plotted against time, given in MJD in the lower x-axis.The upper x-axis shows the long-term phase Θ.For the analysis of this present work, we selected data above 1σ, meaning that the value of their flux density is greater than their uncertainty.Flux variability and a long-term modulation pattern can be seen in the light curve.A quantitative analysis of the flux variability is the subject of our investigation.

Archival data
Since its first radio detection on 1977 August 11 (Gregory & Taylor 1978), LS I +61 • 303 has been the subject of repeated radio observations.Massi & Torricelli-Ciamponi ( 2016) compiled all available 4 to 15 GHz radio data from this source for a detailed timing analysis.Here, we extend this database with the new OVRO observations described above.The resulting long-term light curve is shown in Fig. 2.These observations were carried out at different radio frequencies.For a better visual display in the figure, we rounded the frequencies to half gigahertz steps and color coded them as indicated by the legend in the plot.These archival radio data span the time range from MJD 43367.480until MJD 59961.318,which is a total of 45.5 years.The sampling of the data is however impacted by irregularity from over the years and different observing programs, as is immediately visible in the plot.

Lomb-Scargle periodogram
The data we used for our analysis are not regularly sampled, which is a typical situation in radio astronomy.The average sampling rate of the OVRO data set, however, is 6.7 days, which is well below the Nyquist limit for the range of the shortest periods that we are interested in (P ∼ 26−27 d).More details about the sampling of the OVRO data can be found in Appendix A.
For our timing analysis, we made use of the implementation of the Lomb-Scargle periodogram (Lomb 1976;Scargle 1982) in the astropy Python package (Astropy Collaboration 2013, 2018, 2022).Details about the estimation of the false alarm probability of features in the periodograms can be found in Appendix B. In order to determine the center period and uncertainty of peaks, we carried out functional fits to the periodogram, as described in Appendix C. A potentially more robust method of computing the Lomb-Scargle periodogram has been introduced by Zechmeister & Kürster (2009) and implemented in Python by Czesla et al. (2019).We also applied their method to our data for comparison, and we found that the differences between the results obtained with their method and those obtained with the astropy package are marginal.

Discrete cross-correlation function
A powerful tool to quantify the self-similarity of data is the auto-correlation function, which we realized as the discrete A228, page 3 of 12 The irregular sampling of the data also had to be accounted for in this context.We computed the DCF with the help of a Python implementation of the Edelson & Krolik (1988) method by Robertson et al. (2015).

Intrinsic periods and epoch folding
The orbital phase of the binary system LS I +61 • 303 is defined as where t is the time of observation, t 0 = 43366.275MJD is the epoch of first radio detection of the source, P 1 is the orbital period, and int(x) takes the integer part of x, so Φ is a real number between zero and one.Periastron has been estimated to occur at the orbital phase Φ = 0.23 ± 0.02 (Casares et al. 2005).Gregory (2002) estimated the value of the long-term modulation with high accuracy as P long = 1667 ± 8 d.The long-term phase is obtained by replacing P 1 with P long in Eq. (1).It is occasionally useful to keep the integer part.The phases then represent the number of cycles elapsed since MJD 43366.275.Phases of other periodic signals from LS I +61 • 303 are defined in an analogous way by inserting the corresponding period value.
In order to investigate the shape of the periodic signals, we made use of the epoch folding technique.For this purpose, we plotted the data against the phase of the corresponding period.
A variant of this is to put the data into phase-bins for averaging first.

Spectral analysis
The Lomb-Scargle periodogram of the OVRO data (plotted in Fig. 1) is shown in Fig. 3.The left panel shows the full spectrum.The strongest feature is in the position of the orbital period (∼26.5 days), and a harmonic of this is also present at half the period.And finally there is also a peak in the position of the long-term period (∼1667 days).The right panel of Fig. 3 shows a zoom into the region around the orbital period.Two very distinct peaks are present that are well detached from each other and clearly stand out above the very low noise level.The presence of these two peaks is in agreement with previous findings, as is explained in Sect.5.1.
In order to quantify the positions of the two peaks, we fit the periodogram in the range of interest with the superposition of two Gaussian functions, as described in Appendix C. The estimates of the two periods are P 1 = 26.49± 0.05 d and P 2 = 26.93 ± 0.05 d.Details of the resulting fit parameters are given in the upper part of Table C.1, and the fitted function is plotted as the solid purple curve in Fig. 3.The agreement between the data and the fit is remarkable.
We performed a fit of a single Gaussian to the peak at the long-term period (visible in the left panel of Fig. 3) and obtained the results shown in the lower part of   2008) method has been used, which is the default option of the astropy package.The dotted horizontal line is at the height of a probability of p = 0.1, which is conventionally the upper limit for a periodogram feature to be considered significant (Linnell Nemec & Nemec 1985;Baluev 2008).All three peaks have p-values well below this limit, which means that they are significant.The probabilities for P 1 and P 2 are numerically zero.For P long , we obtained p = 10 −4 , which is still very significant.Besides the Baluev (2008) method, three other options are available in astropy, all of which gave consistent results concerning the significance of the three periods.

Signal profiles
To investigate the profiles of the signals that are periodic with P 1 , P 2 , and P long , we folded the OVRO data over these periods.We put the data into ten phase-bins and averaged the flux values in each bin.The folded light curves are shown in Fig. 4, where the error bars represent the 1σ standard deviation of the distribution in each phase-bin.The phase of the x-axes range from zero to two, with values greater than one repeated from the unit interval.The top panel of the figure shows the data folded over the orbital period P 1 .In the middle panel, the data are folded over the precession period P 2 , and the bottom panel shows the data folded over the long-term period P long .All three profiles show a distinct one-peaked pattern.The orbital profile peaks at Φ ≈ 0.7, the precession profile peaks at Φ 2 ≈ 0.7, and the long-term profile peaks at Θ ≈ 1.0.This is in agreement with previous findings reported in the literature (Jaron et al. 2018) and is an indication of the stability of these patterns.In particular, we find it worth noting that the peak of the orbital profile is still at apastron, with periastron occurring at Φ = 0.23 (Casares et al. 2005).This is not only in agreement with previous observations but also with the physical modeling of Massi & Torricelli-Ciamponi (2014) and Jaron et al. (2016).Folding the data over an arbitrary period does not result in any distinct variability profile, as is shown in Appendix D.

Time-resolved analysis of the OVRO data
The OVRO monitoring of LS I +61 • 303 has now completed three full cycles of the long-term period of the source.In order to analyze the spectral characteristics of each of the three longterm cycles, we divided the full OVRO data set into three segments of equal length.Each segment is 1684.2d long.We performed generalized Lomb-Scargle timing analysis on each one of these segments.The resulting periodograms are shown in Fig. 5.In each of the periodograms, the two-peaked profile of P 1 and P 2 is clearly present, and the period values are in agreement with the ones found in the analysis of the entire OVRO  2. The upper panel shows the entire spectrum, with a prominent feature at the position of the orbital period (and a harmonic at ∼13 d) and a feature at the position of the long-term period and a subharmonic of it.The lower panel shows a zoom into the region around the orbital period, revealing a two-peaked feature.data, as indicated by the two vertical lines, which correspond to the values reported in Sect.4.1.No significant deviation from these values was detected.

Stability of the long-term modulation
In this section, we investigate the stability of the long-term modulation of the radio emission from LS I +61 • 303 by analyzing the long-term radio light curve shown in Fig. 2.This light curve is the combination of the new OVRO data with the entire radio archive from Massi & Torricelli-Ciamponi (2016).
We computed the Lomb-Scargle periodogram of the full data set with the same method that we applied to the OVRO data alone.The result is shown in Fig. 6, where the upper panel shows the full spectrum and the lower panel shows a zoom into the region around the orbital period.The result is very similar to the periodogram of the OVRO data alone (cf.Fig. 3) except that the subharmonic of P long is more pronounced here, gaining almost the same power as the long-term period itself.In the lower panel, the zoom into the orbital period range shows a double-peaked profile.We fit this profile with the superposition of two Gaussian functions, which is analogous to what we did with the OVRO data above.The resulting periods are P 1 = 26.50 ± 0.12 d and P 2 = 26.91 ± 0.11 d.The fitted function appears as the solid purple curve in the figure.We performed a fit of a single Gaussian to the peak at the position of the long-term modulation and obtained P long = 1601 ± 152 d.The detailed results of the fit parameters are listed in Table C.2.
Following up on the analysis by Massi & Torricelli-Ciamponi (2016), we computed the discrete autocorrelation function of the long-term radio data.The result is shown in Fig. 7, where the autocorrelation coefficient is plotted against the time lag.The plot shows an oscillatory pattern with peaks at integer multiples of the long-term period, as highlighted by the red arrows.This shows that the radio data become self-similar when they are shifted against themselves by n • P long , where n is an integer number and P long = 1601 d.With this plot, we confirmed  2).The autocorrelation coefficient is plotted against the lag.Distinct peaks occur at integer multiples of the long-term period, marked by the red arrows.
the result previously obtained by Massi & Torricelli-Ciamponi (2016) in the left panel of their Fig. 5, where they showed that the long-term modulation had remained stable for eight full cycles.Our result now shows that this trend is ongoing and that the longterm modulation has by now remained stable for ten full cycles.

Discussion
The values of the periodic features found by the timing analysis presented here are in very good agreement with previously published results.In particular, the periodogram of the third segment of the OVRO data (i.e., the new data only) shows that the two periods, P 1 and P 2 , are still present in the most recent radio data, and that their values have remained stable (bottom panel of Fig. 5).In this section, we first discuss which physical processes are the possible mechanisms behind these periodic signals.We then show how these periodicities are mathematically connected through interference.And finally, we also briefly discuss properties of the full radio data set on even longer timescales.

Physical processes behind the periodic features
The most precise measurement of the orbital period of the binary system LS I +61 • 303 is still the value P 1 = 26.4960± 0.0028 d determined by Gregory (2002), who applied Bayesian hypothesis testing to the radio data from the GBI monitoring.There are two possible reasons why his analysis did not also result in the detection of the nearby period P 2 ≈ 26.9 d, as the timing analysis by Massi & Jaron (2013) of the same data set did.First of all, Gregory (2002) did not test for the hypothesis of three periodic features (i.e., P 1 , P 2 , and P long )1 .Secondly, as shown by Massi & Jaron (2013), the periodic feature P 2 has a lower power in the periodogram of the GBI data set than P 1 , so if one only tests for one period in that range, the analysis will detect the stronger P 1 and not the weaker P 2 .
Massi & Jaron (2013) interpreted the periodic feature at P 2 = 26.92± 0.07 d as the precession period of a relativistic jet.The reason for this conclusion was that Massi et al. (2012) had previously analyzed a sequence of VLBA observations of LS I +61 • 303 and by applying the method of phase-referenced astrometry, they showed that the core component traces an ellipse and that it takes this ellipse longer than one orbital cycle to return to its initial position.This result was later confirmed by Wu et al. (2018), who revisited the source with VLBI astrometry and showed that the core still traces an ellipse that overlaps with the previous observations by Massi et al. (2012) after correcting for the proper motion of the source.By analyzing the combined data set, Wu et al. (2018) determined the precession period of the core component as P precession = 26.926± 0.005 d.This value is in remarkable agreement with the values from the timing analysis of different data sets at multiple wavelengths (Massi & Jaron 2013;Jaron & Massi 2014;Massi & Torricelli-Ciamponi 2016;D'Aì et al. 2016;Jaron et al. 2018), all of which detected the double-peaked profile of P 1 and P 2 in the power spectra.A physical model of a self-absorbed relativistic jet precessing with a period P 2 and periodically refilled with a population of relativistic electrons with period P 1 is able to reproduce decades of radio observations and their spectral characteristics, as shown by Massi & Torricelli-Ciamponi (2014).The physical process responsible for the modulation of the emission with P 2 is periodic changes in the Doppler boosting of the intrinsic emission by a jet that continuously and periodically changes the angle with respect to the line of sight of an observer on Earth.With the work presented in this article, we confirm that the precession signal is still active in the most recent radio emission from this source (as shown, in particular, in the bottom panel of Fig. 5).
In the periodograms presented in this article (Fig. 3 for the OVRO monitoring and Fig. 6 for the entire radio data set), the most powerful features are found at the positions of P 1 and P 2 , as discussed above.The third prominent feature is a peak at the position of the known long-term period P long .The value of this period was determined with the highest precision by Gregory (2002) as P long = 1667 ± 8 d through Bayesian analysis of the GBI data.The timing analysis of our work here gives results that are, within their uncertainties, in agreement with this measurement.In addition, we repeated the autocorrelation analysis first carried out by Massi & Torricelli-Ciamponi (2016) on the updated data set of radio observations since 1977, the result of which is presented in Fig. 7.This result confirms the previously published result by Massi & Torricelli-Ciamponi (2016), and it shows that the long-term modulation is still active in LS I +61 • 303 and is still stable.

Periodic variability in the Be star disk
The non-degenerate component in the binary system LS I + 61 • 303 is a Be-type star (Hutchings & Crampton 1981).These stars have a high angular momentum, leading to an equatorial outflow in the form of a decretion disk, which is the source of emission lines seen in their optical spectra (see Rivinius et al. 2013 for a review).The definition of the Be phenomenon already includes the possibility that the presence of emission lines can be variable, which implies that the disk itself can be subject to variability over time.Physical processes connected to precisely periodic Be star variability have been observed to occur on timescales of 0.5 to 2 days.To the best of our knowledge, there are not any (quasi-)periodic processes in Be star disks that occur on the timescales on the order of the period P 2 , that is, a few tens of days.For this reason, the possibility that P 2 itself may correspond to any physical process in the Be star disk is ruled out.There are, however, two processes that can occur on timescales on the order of a few years, which would be compatible with the long-term modulation P long of LS I +61 • 303 and is why these processes deserve to be mentioned here.The first process is variability of the size of the Be star disk.The second is a one-armed density wave.Although both types of variability are often observed in Be X-ray binaries, neither of these has ever been observed to occur strictly periodically, especially not over a time span of several decades, which is how the long-term modulation in LS I +61 • 303 occurs.
In X-ray binaries, the decretion disk around a star displaying the Be phenomenon is not the only possible source of line emission.Fender et al. (2009, and references therein) presented observational evidence that the accretion disk around the compact object can be a source of strong emission lines.In their Fig. 1, Fender et al. (2009) present optical spectra of the black hole X-ray binary GX 339-4, which they state is almost certainly uncontaminated by any companion star.The spectra shown there correspond to different X-ray states of the source (faint, brighthard, and bright-soft).While emission lines can be seen in all three presented spectra, they are the strongest in the faint state.In all three states, the strongest feature is the Hα line, but there are also several He lines present.Since LS I +61 • 303 remains in a low-hard X-ray state, on the boundary to the quiescent state (Massi et al. 2017(Massi et al. , 2020)), it is very possible that in this system the accretion disk is also a significant source of line emission.That the optical emission lines indeed originate from a rotating flow is shown in Fig. 7 in Fender et al. (2009), where it is highlighted that most of the observed emission lines have a double-peaked profile, which is discussed in more detail in their Sect.3. In their conclusion, the authors state that "lowluminosity accreting sources should be clearly identifiable in, for example, Hα surveys by their large [equivalent widths]".These observational facts should be kept in mind when interpreting the presence and the timing characteristics of any emission lines from LS I +61 • 303.
If the size of a Be star disk changes, then this manifests itself in a variability of the photometric light curve of the optical emission.In this context, a well-studied system is the Be X-ray binary A0538-66.Rajoelimanana et al. (2017) reported on how quasiperiodic variability in the optical emission from this source, especially pronounced in the V-band data shown in the upper panel of their Fig. 1, is related to variation in the size of the Be star disk.Concerning our target of interest, LS I +61 • 303, such an optical light curve has never been reported anywhere.Furthermore, Rajoelimanana et al. (2017) show periodigrams (their Fig. 3) that reveal that there is only one peak at the orbital period.There is no two-peaked profile, which is unlike the periodograms we report (Figs. 3,5,and 6) and those that have been reported before, as explained in the previous subsection.The absence of a second peak in A0538-66 thus demonstrates that a quasi-periodic modulation of the Be star disk in the form of a gradual buildup and decay (as suggested for LS I +61 • 303 by Chernyakova et al. 2012) does not result in a beating with the orbital period.Another probe of the Be disk size is the equivalent width (EW) of the Hα emission line.However, remembering that in X-ray binaries the accretion disk can also be a source of Hα and other line emission, detections that EW(Hα) is modulated with P long in LS I +61 • 303 (Zamanov et al. 2013;Paredes-Fortuny et al. 2015) are not proof of that variability being related to periodic changes in the Be star disk.The interpretation of EW(Hα) is further complicated by the fact that both the accretion and the jet contribute their variable optical continuum to its normalization (as been pointed out in Jaron 2021).In any case, a strictly periodic variability of a Be star disk size has A228, page 7 of 12 never been observed, so it is an unlikely physical process behind the periodic long-term modulation of LS I +61 • 303.
Cyclic variation in a Be star disk in the form of a one-armed density wave manifests itself as modulation of the V/R (violet to red) ratio of optical emission lines.As Massi & Torricelli-Ciamponi (2016) already pointed out, a wellstudied example of this type of variability is the Be binary system ζ Tau, as examined in detail by Štefl et al. (2009).The database available for this system spans a century and shows that the modulation of the V/R ratio is quasi-periodic at times but can also be completely absent for a few decades.Concerning LS I +61 • 303, a periodic modulation of the V/R ratio of any emission line has never been reported despite the searches performed by Zamanov et al. (1999Zamanov et al. ( , 2013)).The periodic behavior of the long-term modulation of LS I +61 • 303 makes it very unlikely that there is any connection to this type of Be star disk variability.
The period P 2 = 26.926± 0.005 d has been explicitly measured for the precession period of the VLBI core component by Wu et al. (2018).This is why in the following section we explain in detail how the strictly periodic long-term modulation of LS I +61 • 303 fits into the scenario of a beating between P 1 and P 2 .

The long-term period as a beating between orbit and precession
Massi & Jaron (2013) were the first authors who reported the presence of P 1 and P 2 in the radio emission from LS I +61 • 303 and who interpreted the long-term period P long as the result of the interference between the orbital period P 1 and the precession period P 2 in the form of a beating.In the following, we provide a short review of how the periods found in the power spectra, here and in previous publications, fit into the mathematical concept of a beating and how this corresponds to the characteristics of the radio emission from LS I +61 • 303.
In simple terms, the phenomenon of a beating can be understood by considering the sum of two sine functions oscillating at circular frequencies ω 1 and ω 2 , where ω i = 2π/P i .For a simple demonstration of the concept of beating, we omit any phase offsets (that these are in fact relevant in the multiwavelength context is discussed in Jaron 2021).The right-hand side of this equation shows that the sum of the two sine functions can be rewritten as the product of a sine function oscillating at a circular frequency (ω 1 + ω 2 )/2, which is the average of ω 1 and ω 2 , and a cosine term that oscillates at (ω 1 − ω 2 )/2.If the difference between ω 1 and ω 2 is small, then the interference pattern has the form of a beating (i.e., a periodic pattern that is slowly modulated).The reason why we define ω long := ω 1 − ω 2 as twice the frequency of the cosine term in Eq. ( 2) is that what is observed as the long-term modulation in a beating is the envelope of the interference pattern, which oscillates at twice the frequency of the actual cosine term.The reason for this is that the cosine term has a sign flip every half period.When inspecting the interference pattern in detail, one can see that the cosine term indeed has a sign flip during every minimum of the long-term modulation pattern, as we show later in this section.
In the following we show how all of this applies to the radio emission from LS I +61 • 303.The most obvious effect of a beat-ing is the long-term amplitude modulation.This effect is most evident in the GBI data because these data are very well sampled, with several observations per day over a time span of 6.7 years (i.e., 1.5 cycles of the long-term modulation).Analysis of this data set led to the most precise measurement of the long-term period by Gregory (2002).In the work we present here, we have analyzed the OVRO radio light curve, now spanning three cycles of the long-term modulation, and we detect a very clear signal at the long-term period in the periodogram of Fig. 3.That the longterm modulation is a decades-spanning stable periodic signal is demonstrated by the autocorrelation study presented in Fig. 7.
However, amplitude modulation is not the only effect that is expected from a beating.The other effect is that the individual radio outbursts do not occur with either the orbital or the precession frequency but with a frequency that is the average of the two, where ν i = P −1 i .Inserting the most precise measurements of P 1 and P 2 (i.e., P 1 = 26.4960± 0.0028 d, Gregory 2002 and P 2 = 26.926±0.005d, Wu et al. 2018), this results in an expected period of P outburst = 26.709± 0.003 d.That the radio outbursts of LS I +61 • 303 indeed occur with that period has first been discussed by Ray et al. (1997), who reported a period of 26.69±0.02d.Furthermore, the slowly oscillating cosine term on the right-hand side of Eq. ( 2) has zero crossings that occur with a frequency of ω 1 − ω 2 (i.e., with a period of P long ).The effect of this is that the phase of the radio outburst with respect to P outburst performs a phase jump of 0.5 every P long during the minima of the long-term modulation.This has already been pointed out by Jaron & Massi (2013), who demonstrated how this behavior can be used to predict the occurrence of the radio outbursts.In Fig. 8, we show the entire radio archive (i.e., all the data points of the long-term light curve shown in Fig. 2) plotted against the phase of the period P outburst = 26.709d.In addition, we divide the data into phase intervals with respect to which is the phase corresponding to twice the long-term period (i.e., the true period of the cosine term in Eq. ( 2)).The left panel of Fig. 8 shows the data from the interval Θ 2 = 0.75−1.25,and the right panel shows data from the complimentary interval Θ 2 = 0.25−0.75.There is a clear separation between the peak occurrence of these two intervals.The data from the former interval have a peak at a position of Φ avg ≈ 0.2, and the data from the latter interval peak at Φ avg ≈ 0.7.This shows that the radio emission from LS I +61 • 303 over the past 45.5 years has been subject to the 0.5 phase jumps in alternating long-term cycles.This is exactly the behavior that is expected from a beating between the P 1 and P 2 (cf.Jaron & Massi 2013).

Systematic flux modulation on timescales longer than P long ?
A by eye inspection of the full radio light curve shown in Fig. 2 revealed that the peak flux densities reached during each cycle of the long-term modulation are not the same.This is not only the result from the frequency dependency of the flux density (see, e.g., Massi & Kaufman Bernadó 2009), which becomes evident when comparing observations carried out at the same or similar radio frequencies.There seems to be a systematic trend with A228, page 8 of 12 For Θ 2 = 0.75−1.25, the data cluster around Φ avg ≈ 0.2, and for Θ 2 = 0.25−0.75, the data cluster around Φ avg ≈ 0.7.This is expected from the beating between P 1 and P 2 , as explained in Sect.5.3.time.In particular, considering the OVRO light curve at 15 GHz observed from MJD 54908 until the end of the time series (see also Fig. 1) suggests a decline in the amplitude of long-term maxima until about Θ ≈ 9. Considering the whole radio data set, this impression seems to fit into the big picture.The highest flux densities are reached during the long-term cycle around Θ = 4 (or Θ = 5, if one wants to put importance on the few data points at ∼350 mJy).After that, there is a systematic trend of decreasing amplitude until the cycle at around Θ = 9.The next long-term cycle has not been fully observed yet, but the flux densities toward Θ = 10 are again higher, hinting at a reversal of the trend (i.e., toward higher amplitudes).Continued radio monitoring of the source will tell whether this impression is confirmed.A caveat to this interpretation of the radio light curve is of course the missing data between Θ ≈ 5−7.The data points observed at Θ ≈ 6 reaching somewhat lower values do not necessarily contradict the hypothesis of a steadily decreasing trend because these data only cover a fraction of one long-term cycle and most probably do not include the maximum of that cycle.Unfortunately, the data base before Θ ≈ 3.5 (MJD 49200) is not as good as after that point in time, with the data being considerably sparser.Looking at the data by eye, one might have the impression that there is a systematically increasing trend that starts at a minimum at around Θ ≈ 1 (MJD 45000) and lasts until a maximum at Θ = 5, after which the amplitudes follow a decreasing trend.Considering the separation of the flux minima at ∼45 000 MJD and ∼59 500 MJD would imply that this putative systematic modulation has a characteristic timescale, or probably even period, of ∼40 years.
We conclude these considerations by stating that we see evidence for a variability in the amplitude of long-term radio maxima in the radio emission from LS I +61 • 303, with an indication for a systematic modulation with time.There seem to be hints that this pattern is repeating, but with the current database at hand, we cannot make a firm statement as to whether this is another periodicity of the source, which might be called a "super-long-term modulation".understanding of the physical processes at work in this intriguing source.Therefore, continued radio monitoring of LS I +61 • 303 is strongly encouraged.

Fig. 2 .FitFig. 3 .
Fig. 2. Concatenation of all available radio data for LS I +61 • 303 between 4 and 15 GHz.The flux density is plotted against time, and the observed radio frequency is indicated by color.

Fig. 4 .
Fig. 4. Variability profiles of the periodic signals in the OVRO data (shown in Fig. 1).Top: data folded over P 1 , showing the orbital profile.Middle: data folded over P 2 , showing the precession profile.Bottom: data folded over P long , showing the profile of the long-term modulation.

Fig. 5 .
Fig. 5. Lomb-Scargle periodograms of the OVRO data (shown in Fig. 1) divided into three individual segments.The two peaks, which are found in the full data set, are also present in each of the segments.The vertical lines indicate the positions of P 1 and P 2 in the full OVRO data set.

Fig. 7 .
Fig. 7. Autocorrelation function of the entire radio light curve (shown in Fig.2).The autocorrelation coefficient is plotted against the lag.Distinct peaks occur at integer multiples of the long-term period, marked by the red arrows.
Radio light curve from OVRO monitoring of LS I +61 • 303 at 15 GHz.The flux density is plotted against the time of observation given in MJD in the lower axis and in cycles of long-term modulation in the upper axis.Only data points above 1σ are plotted here and were included in the analysis for this work.
Lomb-Scargle periodogram of the entire light curve shown in Fig.