EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 683 (March 2024) A&A, 683 (2024) A228 Full HTML
Open Access
Issue
A&A
Volume 683, March 2024
Article Number A228
Number of page(s) 12
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/202347871
Published online 22 March 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. 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 P1 = 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 Plong = 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 non-thermal emission from LS I +61°303 are powered by the spin-down 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) × 1038 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 1014 G would certainly inhibit the formation of a radio jet (Massi & Kaufman Bernadó 2008 determined an upper limit of 107 − 108 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 ∼108 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°303 carried out with the Owens Valley Radio Observatory (OVRO) at 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.

2. The data

2.1. 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 (2009 March 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.

thumbnail Fig. 1.

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.

2.2. 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.480 until 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.

thumbnail 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.

3. Methods

3.1. 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.

3.2. 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 cross-correlation function (DFC) of the data with themselves. 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).

3.3. Intrinsic periods and epoch folding

The orbital phase of the binary system LS I +61°303 is defined as

Φ = t − t 0 P 1 − int ( t − t 0 P 1 ) , $$ \begin{aligned} \Phi = \frac{t - t_0}{P_1} - \mathrm{int}\left(\frac{t - t_0}{P_1}\right), \end{aligned} $$(1)

where t is the time of observation, t0 = 43366.275 MJD is the epoch of first radio detection of the source, P1 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 Plong = 1667 ± 8 d. The long-term phase is obtained by replacing P1 with Plong 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.

4. Results

4.1. 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.

thumbnail Fig. 3.

Lomb-Scargle periodogram of the OVRO data (shown in Fig. 1). Left: full spectrum showing a peak at the long-term period, a feature at the orbital period, and a harmonic of it at ∼13 d. Right: zoom into the orbital range revealing a two-peak feature.

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 P1 = 26.49 ± 0.05 d and P2 = 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 Table C.1. Our estimate of the long-term period is Plong = 1698 ± 196 d. Again, the fitted function is plotted as the purple solid curve in the figure.

To determine the significance of peaks in the Lomb-Scargle periodogram, we made use of the functionality in the astropy package to compute the false alarm probability p of peak powers. Figure B.1 shows the false alarm probability plotted against peak power. The peak powers of the periods P1, P2, and Plong are marked by the vertical dashed lines. For the data shown in Fig. B.1, the Baluev (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 P1 and P2 are numerically zero. For Plong, 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.

4.2. Signal profiles

To investigate the profiles of the signals that are periodic with P1, P2, and Plong, 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 P1. In the middle panel, the data are folded over the precession period P2, and the bottom panel shows the data folded over the long-term period Plong. 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.

thumbnail Fig. 4.

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

4.3. 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 long-term cycles, we divided the full OVRO data set into three segments of equal length. Each segment is 1684.2 d 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 P1 and P2 is clearly present, and the period values are in agreement with the ones found in the analysis of the entire OVRO 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.

thumbnail 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 P1 and P2 in the full OVRO data set.

4.4. 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 Plong 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 P1 = 26.50 ± 0.12 d and P2 = 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 Plong = 1601 ± 152 d. The detailed results of the fit parameters are listed in Table C.2.

thumbnail Fig. 6.

Lomb-Scargle periodogram of the entire light curve shown in Fig. 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.

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 ⋅ Plong, where n is an integer number and Plong = 1601 d. With this plot, we confirmed 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 long-term modulation has by now remained stable for ten full cycles.

thumbnail 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.

5. 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, P1 and P2, 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.

5.1. 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 P1 = 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 P2 ≈ 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., P1, P2, and Plong)1. Secondly, as shown by Massi & Jaron (2013), the periodic feature P2 has a lower power in the periodogram of the GBI data set than P1, so if one only tests for one period in that range, the analysis will detect the stronger P1 and not the weaker P2.

Massi & Jaron (2013) interpreted the periodic feature at P2 = 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 Pprecession = 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 P1 and P2 in the power spectra. A physical model of a self-absorbed relativistic jet precessing with a period P2 and periodically refilled with a population of relativistic electrons with period P1 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 P2 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 P1 and P2, as discussed above. The third prominent feature is a peak at the position of the known long-term period Plong. The value of this period was determined with the highest precision by Gregory (2002) as Plong = 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.

5.2. 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 P2, that is, a few tens of days. For this reason, the possibility that P2 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 Plong 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, bright-hard, 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, 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 “low-luminosity 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 quasi-periodic 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 Plong 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 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 well-studied 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. (1999, 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 P2 = 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 P1 and P2.

5.3. The long-term period as a beating between orbit and precession

Massi & Jaron (2013) were the first authors who reported the presence of P1 and P2 in the radio emission from LS I +61°303 and who interpreted the long-term period Plong as the result of the interference between the orbital period P1 and the precession period P2 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,

sin ω 1 t + sin ω 2 t = sin ω 1 + ω 2 2 ⏟ ω avg t · cos ω 1 − ω 2 2 ⏟ ω long / 2 t , $$ \begin{aligned} \sin \omega _1t + \sin \omega _2t = \sin \underbrace{\frac{\omega _1 + \omega _2}{2}}_{\omega _{\rm avg}}t\cdot \cos \underbrace{\frac{\omega _1 - \omega _2}{2}}_{\omega _{\rm long}/2}t, \end{aligned} $$(2)

where ωi = 2π/Pi. 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 beating 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 long-term 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,

ν outburst = ν 1 + ν 2 2 = 1 2 ( 1 P 1 + 1 P 2 ) = 1 P outburst , $$ \begin{aligned} \nu _{\rm outburst} = \frac{\nu _1 + \nu _2}{2} = \frac{1}{2}\left(\frac{1}{P_1} + \frac{1}{P_2}\right) = \frac{1}{P_{\rm outburst}}, \end{aligned} $$(3)

where ν i = P i −1 $ \nu_i = P_i^{-1} $. Inserting the most precise measurements of P1 and P2 (i.e., P1 = 26.4960 ± 0.0028 d, Gregory 2002 and P2 = 26.926 ± 0.005 d, Wu et al. 2018), this results in an expected period of Poutburst = 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.02 d. 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 Plong). The effect of this is that the phase of the radio outburst with respect to Poutburst performs a phase jump of 0.5 every Plong 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 Poutburst = 26.709 d. In addition, we divide the data into phase intervals with respect to

Θ 2 = t − t 0 2 P long − int ( t − t 0 2 P long ) , $$ \begin{aligned} \Theta _2 = \frac{t - t_0}{2P_{\rm long}} - \mathrm{int}\left(\frac{t - t_0}{2P_{\rm long}}\right), \end{aligned} $$(4)

thumbnail Fig. 8.

Entire radio data set (shown in Fig. 2) folded over Poutburst = 26.709 d and separated into double long-term (2 ⋅ Plong) phase intervals. 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 P1 and P2, as explained in Sect. 5.3.

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 P1 and P2 (cf. Jaron & Massi 2013).

5.4. Systematic flux modulation on timescales longer than Plong?

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 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”.

6. Conclusions

In this article we have analyzed the updated database of OVRO radio observations of the γ-ray-loud X-ray binary LS I +61°303 at 15 GHz. In addition, we have combined these data with decades-spanning radio observations of this source in order to investigate the long-term stability of periodic signals in the emission. We conclude the following:

  1. We measure the orbital period as P1 = 26.49 ± 0.05 d from analysis of the full OVRO data set alone. Analysis of the entire radio archive gives P1 = 26.50 ± 0.12 d. The larger uncertainty is most likely the result of combining different radio frequencies.

  2. The presence of a period at P2 = 26.93 ± 0.05 d is confirmed from the OVRO data alone, and analysis of the full radio archive results in P2 = 26.91 ± 0.11 d. By combining different frequencies, we obtained a larger uncertainty for the full data set as well.

  3. By analyzing each long-term cycle of the OVRO data separately, we showed that the two periods P1 and P2 are present in each of the three consecutive long-term cycles and that their values remain stable.

  4. The period of the long-term modulation is Plong = 1698 ± 196 d in the OVRO data, and it is Plong = 1601 ± 152 d in the full radio data.

  5. Analysis of the updated long-term radio archive of LS I +61°303 using the auto-correlation function revealed that the long-term modulation has now completed ten full cycles with remarkable stability.

The presence and the stability of the two close periods, P1 and P2, and the agreement of Plong with the beat period of the two combined with the stability of the long-term modulation over ten full cycles strongly favors the interpretation that P2 is the precession period of a relativistic jet. This model has already been successful in reproducing the radio (Massi & Torricelli-Ciamponi 2014) and gigaelectron volt emission (Jaron et al. 2016) from this source. The interpretation of the long-term modulation being the result of changes in the Be star’s decretion disk is unlikely because these variations are not expected to remain so stable over such a long period of time, and this has never been observed (Rivinius et al. 2013; Massi & Torricelli-Ciamponi 2016).

From inspection of the full radio archive spanning almost 46 years, we obtained the strong impression that there is a systematic modulation of the amplitude of the long-term maxima. However, this data set is still not long enough to make a firm statement. Only the availability of more radio data will allow us to find out if there is another signal in the radio emission from this source. Knowledge about this may improve our understanding of the physical processes at work in this intriguing source. Therefore, continued radio monitoring of LS I +61°303 is strongly encouraged.

1

We note that in his publication, Gregory (2002) used the symbol P2 for the long-term period, which we denote as Plong.

Acknowledgments

We thank Gunther Witzel for carefully reading the manuscript and for providing very useful feedback. We thank Maria Massi for reading the manuscript and for valuable input during very useful discussions. This research has made use of data from the OVRO 40-m monitoring program (Richards et al. 2011), supported by private funding from the California Institute of Technology and the Max Planck Institute for Radio Astronomy, and by NASA grants NNX08AW31G, NNX11A043G, and NNX14AQ89G and NSF grants AST-0808050 and AST-1109911. F.J. acknowledges funding by the Austrian Science Fund (FWF) [P35920]. S.K. acknowledges support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 771282.

References

  1. Abbasi, R., Ackermann, M., Adams, J., et al. 2022, ApJ, 930, L24 [NASA ADS] [CrossRef] [Google Scholar]
  2. Aragona, C., McSwain, M. V., Grundstrom, E. D., et al. 2009, ApJ, 698, 514 [Google Scholar]
  3. Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  5. Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167 [NASA ADS] [CrossRef] [Google Scholar]
  6. Baglio, M. C., Coti Zelati, F., Campana, S., et al. 2023, A&A, 677, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Baluev, R. V. 2008, MNRAS, 385, 1279 [Google Scholar]
  8. Bondi, H., & Hoyle, F. 1944, MNRAS, 104, 273 [Google Scholar]
  9. Bosch-Ramon, V., Paredes, J. M., Romero, G. E., & Ribó, M. 2006, A&A, 459, L25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Cañellas, A., Joshi, B. C., Paredes, J. M., et al. 2012, A&A, 543, A122 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Casares, J., Ribas, I., Paredes, J. M., Martí, J., & Allende Prieto, C. 2005, MNRAS, 360, 1105 [NASA ADS] [CrossRef] [Google Scholar]
  12. Chernyakova, M., Neronov, A., Molkov, S., et al. 2012, ApJ, 747, L29 [NASA ADS] [CrossRef] [Google Scholar]
  13. Czesla, S., Schröter, S., Schneider, C. P., et al. 2019, Astrophysics Source Code Library [record ascl:1906.010] [Google Scholar]
  14. D’Aì, A., Cusumano, G., La Parola, V., Segreto, A., & Mineo, T. 2016, MNRAS, 456, 1955 [CrossRef] [Google Scholar]
  15. Deller, A. T., Archibald, A. M., Brisken, W. F., et al. 2012, ApJ, 756, L25 [NASA ADS] [CrossRef] [Google Scholar]
  16. Dhawan, V., Mioduszewski, A., & Rupen, M. 2006, Proceedings of the VI Microquasar Workshop: Microquasars and Beyond, 52.1 [Google Scholar]
  17. Dubus, G. 2013, A&ARv, 21, 64 [Google Scholar]
  18. Durant, M., Kargaltsev, O., Pavlov, G. G., Chang, C., & Garmire, G. P. 2011, ApJ, 735, 58 [NASA ADS] [CrossRef] [Google Scholar]
  19. Edelson, R. A., & Krolik, J. H. 1988, ApJ, 333, 646 [Google Scholar]
  20. Fender, R. P., Russell, D. M., Knigge, C., et al. 2009, MNRAS, 393, 1608 [NASA ADS] [CrossRef] [Google Scholar]
  21. Gregory, P. C. 2002, ApJ, 575, 427 [NASA ADS] [CrossRef] [Google Scholar]
  22. Gregory, P. C., & Taylor, A. R. 1978, Nature, 272, 704 [NASA ADS] [CrossRef] [Google Scholar]
  23. Hutchings, J. B., & Crampton, D. 1981, PASP, 93, 486 [NASA ADS] [CrossRef] [Google Scholar]
  24. Jaron, F. 2021, Universe, 7, 245 [NASA ADS] [CrossRef] [Google Scholar]
  25. Jaron, F., & Massi, M. 2013, A&A, 559, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Jaron, F., & Massi, M. 2014, A&A, 572, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Jaron, F., Torricelli-Ciamponi, G., & Massi, M. 2016, A&A, 595, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Jaron, F., Massi, M., Kiehlmann, S., & Hovatta, T. 2018, MNRAS, 478, 440 [NASA ADS] [CrossRef] [Google Scholar]
  29. Kravtsov, V., Berdyugin, A. V., Piirola, V., et al. 2020, A&A, 643, A170 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Linnell Nemec, A. F., & Nemec, J. M. 1985, AJ, 90, 2317 [NASA ADS] [CrossRef] [Google Scholar]
  31. Lomb, N. R. 1976, Ap&SS, 39, 447 [Google Scholar]
  32. Maraschi, L., & Treves, A. 1981, MNRAS, 194, 1P [Google Scholar]
  33. Marti, J., & Paredes, J. M. 1995, A&A, 298, 151 [NASA ADS] [Google Scholar]
  34. Massi, M., & Jaron, F. 2013, A&A, 554, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Massi, M., & Kaufman Bernadó, M. 2008, A&A, 477, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Massi, M., & Kaufman Bernadó, M. 2009, ApJ, 702, 1179 [NASA ADS] [CrossRef] [Google Scholar]
  37. Massi, M., & Torricelli-Ciamponi, G. 2014, A&A, 564, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Massi, M., & Torricelli-Ciamponi, G. 2016, A&A, 585, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Massi, M., Ros, E., & Zimmermann, L. 2012, A&A, 540, A142 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Massi, M., Migliari, S., & Chernyakova, M. 2017, MNRAS, 468, 3689 [Google Scholar]
  41. Massi, M., Chernyakova, M., Kraus, A., et al. 2020, MNRAS, 498, 3592 [NASA ADS] [CrossRef] [Google Scholar]
  42. Mirabel, I. F., & Rodríguez, L. F. 1999, ARA&A, 37, 409 [Google Scholar]
  43. Paredes-Fortuny, X., Ribó, M., Bosch-Ramon, V., et al. 2015, A&A, 575, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Rajoelimanana, A. F., Charles, P. A., Meintjes, P. J., et al. 2017, MNRAS, 464, 4133 [NASA ADS] [CrossRef] [Google Scholar]
  45. Ray, P. S., Foster, R. S., Waltman, E. B., Tavani, M., & Ghigo, F. D. 1997, ApJ, 491, 381 [NASA ADS] [CrossRef] [Google Scholar]
  46. Richards, J. L., Max-Moerbeck, W., Pavlidou, V., et al. 2011, ApJS, 194, 29 [Google Scholar]
  47. Rivinius, T., Carciofi, A. C., & Martayan, C. 2013, A&ARv, 21, 69 [Google Scholar]
  48. Robertson, D. R. S., Gallo, L. C., Zoghbi, A., & Fabian, A. C. 2015, MNRAS, 453, 3455 [NASA ADS] [Google Scholar]
  49. Romero, G. E., Okazaki, A. T., Orellana, M., & Owocki, S. P. 2007, A&A, 474, 15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Scargle, J. D. 1982, ApJ, 263, 835 [Google Scholar]
  51. Slane, P. 2017, in Handbook of Supernovae, eds. A. W. Alsabti, & P. Murdin, 2159 [CrossRef] [Google Scholar]
  52. Suvorov, A. G., & Glampedakis, K. 2022, ApJ, 940, 128 [NASA ADS] [CrossRef] [Google Scholar]
  53. Taylor, A. R., Kenny, H. T., Spencer, R. E., & Tzioumis, A. 1992, ApJ, 395, 268 [NASA ADS] [CrossRef] [Google Scholar]
  54. Torres, D. F., Rea, N., Esposito, P., et al. 2012, ApJ, 744, 106 [NASA ADS] [CrossRef] [Google Scholar]
  55. Štefl, S., Rivinius, T., Carciofi, A. C., et al. 2009, A&A, 504, 929 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Weng, S.-S., Qian, L., Wang, B.-J., et al. 2022, Nat. Astron., 6, 698 [NASA ADS] [CrossRef] [Google Scholar]
  57. Wu, Y. W., Torricelli-Ciamponi, G., Massi, M., et al. 2018, MNRAS, 474, 4245 [NASA ADS] [CrossRef] [Google Scholar]
  58. Zamanov, R. K., Martí, J., Paredes, J. M., et al. 1999, A&A, 351, 543 [NASA ADS] [Google Scholar]
  59. Zamanov, R., Stoyanov, K., Martí, J., et al. 2013, A&A, 559, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  60. Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577 [CrossRef] [EDP Sciences] [Google Scholar]
  61. Zimmermann, L., Fuhrmann, L., & Massi, M. 2015, A&A, 580, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

Appendix A: OVRO monitoring sampling rate

Figure A.1 shows a histogram of the time gaps Δt between adjacent observations within the OVRO monitoring. We only show values until 45 days for better visual display. There are five larger gaps extending until 110 days, which means that gaps of that length are an exception.

thumbnail Fig. A.1.

Histogram of the time gaps between two adjacent observations within the OVRO data. We only show values until 45 days. There are five longer gaps, the longest being 110 days.

Appendix B: Lomb-Scargle false alarm probability

There are four methods for estimating the false alarm probability of periodic features in the Lomb-Scargle implementation in astropy. We tried all four of them, and they gave consistent results. In Fig. B.1, we show the results for the Baluev (2008) method, plotting false alarm probability against periodogram peak power. The powers of all three periodic features, P1, P2, and Plong, are shown by the vertical dashed lines. The period Plong has a value of ∼10−4, which is still highly significant. The two periods P1 and P2 have values that are practically in agreement with zero. We are aware that the use of the false alarm probability is only valid in the case of white noise. Indeed, the periodograms obtained from the data analyzed in this article feature only the discussed peaks and otherwise show a very flat noise spectrum.

thumbnail Fig. B.1.

False alarm probability plotted against Lomb-Scargle peak power. The vertical dashed lines mark the peak powers of the three periods found in the periodogram. The dotted horizontal line marks a probability of p = 0.1. For the data shown here, the Baluev (2008) method has been used as implemented in astropy.

Appendix C: Fitting the periodogram

In order to quantify the positions of peaks in the periodograms and their widths, we performed least squares fitting of functions to the periodograms. In general, periodogram peaks are best fit with a Gaussian function of the form

f ( P ) = A e − 1 2 ( P − P i ) 2 σ 2 , $$ \begin{aligned} f(P) = A\,\mathrm{e}^{-\frac{1}{2}\frac{(P - P_i)^2}{\sigma ^2}}, \end{aligned} $$(C.1)

with amplitude A, center period Pi, and standard deviation σ. We took the center period Pi as the period value and its standard deviation σ as its uncertainty in the period values that we report throughout this article. We find that the two-peaked profile of the Lomb-Scargle periodogram is best fit with the function

f ( P ) = A 1 e − 1 2 ( P − P 1 ) 2 σ 1 2 + A 2 e − 1 2 ( P − P 2 ) 2 σ 2 2 , $$ \begin{aligned} f(P) = A_1\,\mathrm{e}^{-\frac{1}{2}\frac{(P - P_1)^2}{\sigma _1^2}} + A_2\,\mathrm{e}^{-\frac{1}{2}\frac{(P - P_2)^2}{\sigma _2^2}}, \end{aligned} $$(C.2)

which is the superposition of two Gaussian functions with amplitudes A1, 2; standard deviations σ1, 2; and center periods P1, 2. In this way, we also prove that the two-peaked profile is significantly fitted by a function that is intrinsically two peaked. The results are listed in Table C.1 for the OVRO data alone and in Table C.2 for the periodogram of the entire archive of radio observations.

Table C.1.

Parameters resulting from fitting functions to the periodogram of the OVRO data alone. The upper part of the table shows the results of fitting the superposition of two Gaussians to the double-peak feature. The lower part shows the results from fitting the long-term period with a single Gaussian.

Table C.2.

Same as Table C.1 but for the periodogram of all radio data.

Appendix D: Folding over a random period

In Sect. 4.2, we folded the OVRO data over the periods found by the timing analysis. The results of this are shown in Fig. 4. For comparison, we folded the data over an arbitrary period of 30 days, shown in Fig. D.1. Compared to the data folded over the period peak values, the data folded over the arbitrary period does not show any distinct profile. Figures 4 and D.1 share the same scale of the y-axis, which makes it obvious that the variability pattern of the arbitrarily folded data is very flat and does not deviate significantly from the fitted constant, which appears as the solid horizontal line in the plot.

thumbnail Fig. D.1.

OVRO data folded over the arbitrary period of 30 days. The profile does not show any regular pattern and is reasonably well fitted by a constant, as shown by the solid line.

All Tables

Table C.1.

Parameters resulting from fitting functions to the periodogram of the OVRO data alone. The upper part of the table shows the results of fitting the superposition of two Gaussians to the double-peak feature. The lower part shows the results from fitting the long-term period with a single Gaussian.

In the text
Table C.2.

Same as Table C.1 but for the periodogram of all radio data.

In the text

All Figures

thumbnail Fig. 1.

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.
In the text
thumbnail 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.
In the text
thumbnail Fig. 3.

Lomb-Scargle periodogram of the OVRO data (shown in Fig. 1). Left: full spectrum showing a peak at the long-term period, a feature at the orbital period, and a harmonic of it at ∼13 d. Right: zoom into the orbital range revealing a two-peak feature.
In the text
thumbnail Fig. 4.

Variability profiles of the periodic signals in the OVRO data (shown in Fig. 1). Top: data folded over P1, showing the orbital profile. Middle: data folded over P2, showing the precession profile. Bottom: data folded over Plong, showing the profile of the long-term modulation.
In the text
thumbnail 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 P1 and P2 in the full OVRO data set.
In the text
thumbnail Fig. 6.

Lomb-Scargle periodogram of the entire light curve shown in Fig. 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.
In the text
thumbnail 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.
In the text
thumbnail Fig. 8.

Entire radio data set (shown in Fig. 2) folded over Poutburst = 26.709 d and separated into double long-term (2 ⋅ Plong) phase intervals. 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 P1 and P2, as explained in Sect. 5.3.
In the text
thumbnail Fig. A.1.

Histogram of the time gaps between two adjacent observations within the OVRO data. We only show values until 45 days. There are five longer gaps, the longest being 110 days.
In the text
thumbnail Fig. B.1.

False alarm probability plotted against Lomb-Scargle peak power. The vertical dashed lines mark the peak powers of the three periods found in the periodogram. The dotted horizontal line marks a probability of p = 0.1. For the data shown here, the Baluev (2008) method has been used as implemented in astropy.
In the text
thumbnail Fig. D.1.

OVRO data folded over the arbitrary period of 30 days. The profile does not show any regular pattern and is reasonably well fitted by a constant, as shown by the solid line.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.