A Precessing Jet Scenario for the Multi-Wavelength Long-Term Modulation of LS I +61°303

Abstract

The high-mass X-ray binary LS I +61°303 is detected across the electromagnetic spectrum from radio until the very high energy $\gamma$-ray regime. The emission is not only highly variable on many time scales, but is also periodic at all observed wavelengths. Periodic modulation was observed on different time-scales, ranging from hours, over months to several years. The subject of this article is a super-orbital, long-term modulation of ∼4.6 years. We review the observation of this periodic modulation at multiple wavelengths and investigate systematic relationships between them. IN particular, radio observations reveal that the long-term modulation is a very stable feature of the source. Observations at other wavelengths result in a phase-shift of the modulationpattern that is a systematic function of energy. The stability of this period favors a scenario in which the long-term modulation is the result of a precessing jet giving rise to periodic changes in the Doppler factor, beating with the orbital modulation of the accretion rate. We explain the phase-shifts across energy bands in a scenario with shorter wavelengths originating closer to the base of the presessing jet. A significant deviation of the TeV emission from this trend possibly requires a different explanation related to magnetic reconnection events.

1. Introduction

The high-mass X-ray binary LS I +61°303 is composed of a B0 Ve star [1] and a black hole candidate [2]. It was discovered in 1978 during a search for variable radio sources [3]. In the course of its observational history, spanning more than four decades now, this source has proven variable and periodic on different time scales not only in radio but all across the electromagnetic spectrum into the very high energy $\gamma$-ray regime. Variability time scales which have been observed in this source range from hours [4,5,6,7], over approximately one month ([8,9] and references therein), up to several years ([8,9,10] and references therein). The variability of hour time scales was observed in the radio and X-ray emission of the source and is superimposed on large outbursts which occur once per orbital cycle. This short-term (intra-day) variability is sometimes characterized by periodic oscillations, but these are only of quasi-periodic nature [4,5,6]. The occurrence of the larger outbursts themselves, in contrast, is periodic and precisely predictable [11], clearly related to the orbital period $P_1=26.4960\pm 0.0028$ d of the binary system [8]. However, the emission from LS I +61°303 is periodic on an even longer time scale: both amplitude and orbital phase occurrence of the radio outbursts are modulated with a period of $P_{\mathrm{long}}=1667\pm 8$ d [8], and there is observational evidence that the same long-term modulation is active at other wavelengths too, in particular at X-rays [12], GeV [13], and TeV [14]. This long-term modulation (LTM hereafter) and its behavior at multiple wavelengths of the electromagnetic spectrum is the subject of this article.

Since the discovery of the source LS I +61°303 the nature of the compact object in this system and the physical processes behind its non-termal emission had been debated. The compact object has been suggested to be either a neutron star or a black hole. And also for the physical processes two scenarios were discussed: pulsar binary and microquasar. In the first scenario the compact object has to be a neutron star and the physical processes are powered by the spin-down of a millisecond pulsar, the non-thermal emission being produced in the collision between the winds of that putative pulsar and the Be star ([15] and references therein). This hypothesis is part of a generalization of the processes observed during the periastron passages of the pulsar binary PSR B1259-63 [16], occurring once every ∼3.4 years ([17] and references therein), to a larger class of binary systems and is based exclusively on the $\gamma$-ray loudness of these sources [18]. As far as LS I +61°303 is concernced, recent observations of the system’s X-ray properties clearly show that this source fits into the class of accreting black hole binaries ([2] and references therein). Furthermore, its radio-loudness and the structure of its periodically occurring radio outbursts [19] as well as the transition from optically thick to optically thin of the radio spectral index during these outbursts [20] support the presence of a jet in LS I +61°303 and complete the microquasar scenario. The afore mentioned quasi-periodic intra-day variabliy at radio and X-rays has also been explained in this scenario [6,9]. Microquasars are stellar binary systems in which a compact object accretes matter from a companion star, an accretion disk is formed, and a jet is launched, very similar to the processes occurring in the nuclei of radio-loud active galaxies, but time-scales of the intrinsic physical processes are much shorter because they scale with the mass of the accretor [21]. In 2020, Massi et al. [22] reported an observational campaign on LS I +61°303 at multiple wavelengths and show how the evolution of the emission traces accretion and ejection during one orbital cycle of the binary system.

The presence of a jet in LS I +61°303 is not only indirectly proven by the properties of its radio emission, but there is also direct evidence provided by radio interferometric observations. Massi et al. [23] analyzed a sequence of VLBA observations, clearly showing an elongated structure which changes its morphology from one-sided to two-sided and which also changes its position angle as time evolves. This not only shows that there is a visible radio jet in LS I +61°303 but is also an indication that this jet precesses. These authors determine a precession period of 27–28 d, which is not equal to but larger than the orbital period of the binary system. The value of this precession period has since been detected by timing analysis of radio [24], X-ray [25], and GeV [9,26] data. The most precise value of this period comes again from VLBI observations: In 2018, Wu et al. [27] revisited LS I +61°303 with VLBI astrometry and determined a value of $P_2=26.926\pm 0.005$ d.

Massi & Jaron [24] in 2013 were the first authors to point out that the interference between the orbital period $P_1$ and the precession period $P_2$ results in a beating with a period identical to the period of the LTM. In this scenario, the accretion rate onto the compact object is periodically modulated with the orbital period because of the eccentricity of the orbit [28,29,30,31]. The precession of the relativistic jet gives rise to periodic changes in the amplification of its intrinsic emission as a result of variable Doppler boosting because of periodic changes of the angle between the jet bulk velocity and the line of sight. This scenario has been quantitatively confirmed by a physical model of a self-absorbed jet which precesses and is periodically refilled with relativistic electrons producing synchrotron emission [32]. These authors show that this model reproduces the flux densities and the spectral index, both as a function of time, of almost four decades of observational radio data. By extending this model to include inverse Compton scattering of both internal synchrotron photons (synchrotron self-Compton) as well as external ultra-violet photons from the companion star (external inverse Compton), Jaron et al. [33] successfully reproduce several years of simultaneous radio and GeV data resulting from long-term monitoring by the Owens Valley Radio Observatory (OVRO) at 15 GHz and the Fermi Large Area Telescope (Fermi-LAT) in the energy range 0.1–3 GeV, respectively. Analyzing continued observational data from these two ongoing monitoring programs, Jaron et al. [9] not only showed that the presence of the two periodicities $P_1$ and $P_2$ is still evident at both wavelengths, but furthermore detected a phase-offset in the LTM patterns between the radio and GeV emission. Offsets between the phases of the LTM pattern had previously been reported for other wavelengths too, as explained in the following Section 2, but in 2018 Jaron et al. [9] were the first to propose a physical scenario to explain such an offset: for the radio and GeV emission, the LTM phase-offset fits into the scenario of a precessing jet in which the higher energy emission (GeV) is emitted in a jet region upstream from the position of the lower energy (radio) emission. The reason for this interpretation is that it is the precession profile of the emission at these two wavelengths that is shifted in phase with respect to each other by an amount which agrees with the shift in the interference (i.e., LTM) pattern. The locations of these two emission regions fit into a cooling model, as e.g., used by Lisakov et al. [34] to explain the regions of the emission at these two wavelength in the jet of the active galactic nucleus 3C 273. However, in addition, opacity effects can play a role similar to the core-shift effect at radio frequencies [35,36].

The aim of this paper is to investigate how LTM phase-offsets at other wavelengths of the emission of LS I +61°303 could fit into a similar scenario. For this purpose, we are going to review how long-term observations of this source resulted in the detection of the LTM at multiple wavelengths and investigate systematic relationships between these observations. The article is structured as follows. In Section 2 we review the observational evidence for a LTM being active in the source LS I +61°303 at multiple wavelengths, namely at radio, optical, X-rays, GeV, and TeV. Systematic phase-offsets between the LTM pattern at different wavelengths are quantitatively investigated in Section 3. We discuss our results in the context of a physical scenario of periodic accretion and ejection in combination with jet precession in Section 4. We give our conclusions in Section 5.

2. Observations of the LTM at Multiple Wavelengths

Since its first radio detection in 1978 [3], LS I +61°303 has been subject of monitoring programs at different wavelengths across the electromagnetic spectrum, from radio until the very high energy $\gamma$-ray regime. At all of these wavelengths a periodic modulation of approximately 4.6 years has been detected. This section provides a review of the observational evidence for the LTM at each wavelength.

The following two definitions are conventionally used. The orbital phase as a function of time t is defined as

$$\mathsf{\Phi }=\frac{t-T_0}{P_1}-\mathrm{int}\left(\frac{t-T_0}{P_1}\right),$$

(1)

where $P_1=26.4960\pm 0.0028$ d is the orbital period of the system and $T_0$ = 43,366.275 d is the MJD of the first radio detection of the source [8]. Periastron passage occurs at orbital phase $\mathsf{\Phi }=0.23$ [1]. The phase of the LTM is defined in an analogous way, i.e.,

$$\mathsf{\Theta }=\frac{t-T_0}{P_{\mathrm{long}}}-\mathrm{int}\left(\frac{t-T_0}{P_{\mathrm{long}}}\right),$$

(2)

where $P_{\mathrm{long}}=1667\pm 8$ d is the period of the LTM [8]. In both cases $\mathrm{int}\left(x\right)$ takes the integer part of x. If this term is omitted, then the phase has the meaning of orbital or long-term cycles elapsed since $T_0$, which is often useful for calculations.

2.1. Radio

The LTM of LS I +61°303 was first discovered in its radio emission. Already in 1987, Paredes [37] suspected in his PhD thesis that there could be a four year LTM. Then in 1989, based on the analysis of ten years of observational radio data in the centimeter regime, Gregory et al. [38] reported that periodic radio outbursts, which had been known to occur once per orbit [39], are strongly modulated in amplitude on a super-orbital, long-term time-scale. These authors confirm an approximate value of four years for the period of this LTM, and they already discussed the possibility that this long-term amplitude modulation could be the result of variable Doppler boosting caused by a similar long-term precessing jet. However, they discarded this option in favour of “variable accretion due to quasi-cyclic Be star envelope variations” (their Sect. V b). The reason they disfavor the precessing jet scenario is an asymmetry of the LTM pattern that they observe. Such an asymmetry is indeed unexpected in a scenario which assumes the precession period of the jet to be identical to the ∼4 years period of the LTM. However, as we will see later on in this artice, a precessing jet model with a precessing period close to the orbital one reproduces the observed flux densities. The presence of the LTM was then confirmed again in 1995 by Martí & Paredes [28], who added the observation that the orbital phase of the periodic outbursts is also affected by the LTM. From 1994 until 2000, LS I +61°303 was part of a monitoring program of the Green Bank Interferometer (GBI) by the National Radio Astronomical Observatory (NRAO). By observing the source several times per day, simultaneously at 2.25 and 8.3 GHz, a database was created which remains to be highly valuable for timing analysis and the investigation of changes in the radio spectral index. In 2002, Gregory [8] applied Bayesian analysis, based on hypothesis testing, to the GBI data and determined a refined value of $P_{\mathrm{long}}=1667\pm 8$ d (4.6 years) for the period of the LTM1. In 2016, Massi & Torricelli-Ciamponi [10] combined the entirety of the radio observations of LS I +61°303 that were available at that time. Timing analysis of this archive of 36.8 years of observational data not only confirmed the value of the long-term period, determined by Gregory [8] before, but most importantly proved that the value of the long-term period had remained stable since the first radio detection of the source, which is eight full cycles of the LTM. Massi & Torricelli-Ciamponi [10] point out that variations in the circumstellar disks of Be stars, as had been proposed by Gregory et al. [38] but also other authors, are never observed to remain stable over such a long period of time and refer to the review paper by Rivinius et al. [40]. That the LTM is still active in the radio emission from LS I +61°303 at 15 GHz (OVRO) has been shown by Jaron et al. [9].

2.2. Optical

First Zamanov et al. [41] in 1999 and then Zamanov & Martí [42] in 2000 reported a LTM of EW(H$\alpha$), i.e., the equivalent width of the H$\alpha$ emission line. Zamanov & Martí [42] also report a long-term phase-shift of about 0.25 with respect to the radio emission2. In 2015 Paredes-Fortuny et al. [43] reported on a phase-shift of the orbital phase occurrence of peaks of HW(H$\alpha$), a phenomenon that is known to be part of the LTM at radio and GeV wavelengths ([26] and references therein). However, the conclusions by Paredes-Fortuny et al. [43] are based on only 1.5 years of observational data, i.e., only about one third of the 4.6 years LTM.

A direct comparison of the 0.25 LTM phase-offset reported by Zamanov & Martí [42] to the other phase-shifts reported in our present investigation is complicated by the fact that they used a different value of 1584 d for the period of the LTM. Therefore, before these data could be used for our investigation they would have to be re-analyzed taking into account current knowledge about the timing characteristics of the emission from LS I +61°303, in particular the value [8] and the stability [9,10] of the LTM period. However, there is also a more fundamental reason that leads us to the decision not to include any EW(H$\alpha$) data: it would require a thorough investigation to disentangle the contributions from different emission regions within the stellar binary system to EW(H$\alpha$). Such an investigation is beyond the scope of this article, but a discussion of this issue in connection with optical emission regions will be given in Section 4.3.

2.3. X-rays

X-ray emissions from LS I +61°303 have been detected since 1981 [44], but the first authors to report the LTM in this energy range were Li et al. [12] in 2012. They analyzed four years of RXTE-PCA data of the energy range 3–30 keV. In the same publication they report a shift between the LTM at X-rays compared to the radio emission. They quantify this shift to $281.8\pm 44.6$ d, which translates to $0.17\pm 0.03$ in terms of long-term phase-offset with respect to the radio emission, meaning that the X-ray modulation lags the radio. The left panel of their Figure 1, showing the raw peak count rates put into long-term phase bins, reveals that the envelope of the LTM pattern in the X-rays has the same asymmetry as has also been observed in the radio. It has the characteristic shape of a steep rise and a slower decline. The phase-offset cited here was obtained by the authors by fitting a sine wave to the “modulated flux fraction” (MFF). The effect of this MFF is that indeed the long-term modulation of it looks more sinusoidal than the original X-ray count peaks. In reality the original count peaks look much more similar to the radio data. Indication that the mechanism behind the LTM in the X-ray emission is the beating between orbit and precession comes from the fact that D’Ai et al. [25] detected the corresponding spectral features of orbital ($P_1$) and precession ($P_2$) periods in the light curves of Swift/BAT survey data.

In 2014, Li et al. [45] analyzed INTEGRAL X-ray data from a slightly higher energy interval (18–60 keV). They detect the LTM also in this energy range and report that the LTM pattern has a peak at $\mathsf{\Theta }\approx 0.2$. Since the LTM at radio has its maximum at $\mathsf{\Theta }=0.94\pm 0.02$ [9], this corresponds to a phase-shift of ∼0.26. However, because this is only a rough estimate without reported uncertainty and because the energy range overlaps with the previously cited Li et al. [12] we do not include this data point in our quantitative analysis later on in this article.

2.4. GeV

Since the source had been detected in the radio [3], the possibility of $\gamma$-ray emission from LS I +61°303 was discussed as early as in 1979 [46]. A firm association with a GeV source, previously detected by Hermsen et al. [47] in 1977, could only be established in 2009, when Abdo et al. [48] detected the orbital period in the emission observed by the Fermi-LAT. The Fermi-LAT has also been the first mission to provide continuous monitoring in the GeV regime, and in this way, it enabled the investigation of the long-term behavior of the source in this energy range.

Other than in the radio, where the LTM is already evident by an eye-inspection of the raw observational data, detecting the same modulation in the GeV required a closer look at the data. By analyzing several years of Fermi-LAT data but selecting only certain orbital phase ranges, in 2013 Ackermann et al. [13] detected for the first time a significant long-term amplitude modulation which, however, only affects apastron orbital phases ($\mathsf{\Phi }$ = 0.5–1.0). Jaron & Massi [26] found out that the GeV emission has indeed two distinct peaks along the orbit. The first peak occurs at periastron and remains stable in orbital phase. The second peak occurs at a later orbital phase and is subject to the same phase-shift as observed for the radio outbursts. In combination, the results obtained by Ackermann et al. [13] and Jaron & Massi [26] completed the picture that the LTM has the same period and the same effects, in terms of amplitude modulation and orbital phase-shift of the outbursts, as is the case for the radio emission. In addition to that, there is a GeV outburst at periastron, which remains stable in amplitude and orbital phase, while the radio emission is lacking a distinct periastron outburst.

The LTM remains active in both the radio and the GeV regime, as shown by Jaron et al. [9]. Furthermore, these authors detect a phase-shift between the LTM at radio and GeV, which they quantify to $0.26\pm 0.03$3. As an explanation for this phase-offset they suggest a scenario in with the GeV emission is produced in a precessing jet upstream from the radio emission. We will come back to this scenario and extend it to the multi-wavelength picture further down in this article.

2.5. TeV

The source LS I +61°303 is also detected in the very high energy (VHE) $\gamma$-ray regime. Here we first summarize the history of observations of the source with Cherenkov telescopes. As we will see, the LTM has been detected in the TeV emission, but a quantitative determination of its phase-offset and uncertainties of the fit parameters has so far been missing from the literature. In order to obtain these values for our investigation here, we repeat a previously published analysis in the second part.

2.5.1. Observational History

A VHE counterpart to LS I +61°303 was first detected by Albert et al. [49]. The source has since been included in the monitoring by MAGIC and VERITAS see references in [14]. The LTM of LS I +61°303 was detected in the VHE $\gamma$-ray regime by Ahnen et al. [14] in 2016 by combining observations from both of these VHE observatories. Selecting only peak flux densitities from the orbital phase interval $\mathsf{\Phi }$ = 0.50–0.75, they showed that these data are best fit with a sine function. They report a period of this sine of $1610\pm 58$ d, which agrees well within uncertainties with the value of $1667\pm 8$ d determined from radio observations [8]. However, Ahnen et al. [14] do not mention the values of the other fit parameters.

2.5.2. Determination of the LTM Phase-Offset

To determine the value of the long-term phase-offset in the TeV regime, we repeat here the analysis of Ahnen et al. [14]. We take the data that they present in their Figure 2, plotted in Figure 1 here. We fit these data with a sine function of the form

$$f\left(\mathsf{\Theta }\right)=asin2\pi (\mathsf{\Theta }-{\mathsf{\Theta }}_0)+b,$$

(3)

which is a function of the long-term phase $\mathsf{\Theta }$, with amplitude a, offset b, and long-term phase-origin ${\mathsf{\Theta }}_0$. The long-term phase $\mathsf{\Theta }$ is defined in Equation (2). We keep the value of the long-term period $P_{\mathrm{long}}=1667$ d fixed to the value from the radio. This approach is justified by the fact that the value determind by Ahnen et al. [14] for the TeV agrees with the value for the radio. In addition to that, this enables the comparison between a phase-offset between TeV and radio, i.e., the subject of this article.

We confirm that the data are well fitted by the sine function of the form of Equation (3). The reduced ${\chi }^2$ from our analysis is $1.6$, and the values of the fit parameters are listed in

Table 1. Values of the parameters resulting from fitting the TeV data of Figure 1 with a sine function of the form given in Equation (3).

Parameter	Value
a	$(8.18\pm 1.11){10}^{-12}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}$
b	$(10.36\pm 1.11){10}^{-12}{\mathrm{cm}}^{-2}{\mathrm{s}}^{-1}$
${\mathsf{\Theta }}_0$	$0.20\pm 0.03$

. We obtain a value of ${\mathsf{\Theta }}_0=0.20\pm 0.03$ for the phase-origin of the LTM in LS I +61°303 at TeV energies. Relative to the phase-origin of the LTM at radio ${\mathsf{\Theta }}_0=0.69$ determined by Jaron et al. [9] (called ${\varphi }_0$ in the last row of their Table 2), this translates to a phase-offset of $0.51\pm 0.03$ between the LTM at TeV and radio. At this point there is the question of ambiguity of the phase-offset. The same offset could also have been interpreted as $-0.49$, i.e., the TeV being earlier than the radio by that amount. However, as evident from Figure 2, a positive TeV-radio offset leads to the best alignment with the offsets at the other wavelengths. As we will see in Section 4.1, this value of $0.51$ for the TeV phase-offset also fits into a physical scenario in which these LTM phase-offsets correspond to energy-dependent emission regions in a precessing jet. However, as we will see in Section 4.2, the TeV photons may be produced by a different mechanism, and the ambiguity resolution of this phase-offset value would then not be so obvious any more.

This concludes our review of observations of the LTM in LS I +61°303 at multiple wavelengths. In the next section, we are going to put the relative phase-offsets in context and investigate systematics between them.

3. Phase-Offsets across Wavelengths

All the known values of the LTM phase-offsets at different wavelengths of the emission from LS I +61°303 are listed in

Table 2. Phase-offsets of the long-term modulation, as observed at different wavelengths of the emission from LS I +61°303. All phase-offset values are given relative to the long-term phase of the radio emission. The higher energy emission lags the radio in terms of LTM.

Name	Energy Range [eV]	Phase-Offset	Reference
Radio	(5.0–7.4) $\times {10}^{-5}$	$0\pm 0.02$	Jaron et al. [9]
X-rays (soft)	(3–30) $\times {10}^3$	$0.17\pm 0.03$	Li et al. [12]
X-rays (hard)	(18–60) $\times {10}^3$	∼0.2	Li et al. [45] (not used here)
GeV	(0.1–3.0) $\times {10}^9$	$0.26\pm 0.02$	Jaron et al. [9]
TeV	(0.3–10) $\times {10}^{12}$	$0.51\pm 0.03$	This work, based on [14]

. The energy range reported there is the one of the emission used to determine the phase-offset at that wavelength. The LTM pattern of the higher energy emission lags the pattern of the radio emission, the offset of which we define as zero. Figure 2 presents a plot of the values given in Table 2. Plotted is long-term phase-offset against energy of the observed emission, the energy axis appearing in a logarithmic scale. It is evident from this plot that the phase-offset monotonically increases as a function of energy, meaning that the lag of the LTM pattern, with respect to radio, increases with energy. At this point it has to be mentioned that there is an ambiguity affecting all of these offsets. For example, the value for the X-rays could also be interpreted as $-0.83$, i.e., leading the radio emission by that amount, and analogous for the other wavelengths. However, resolving these ambiguities to the values reported here results in the smallest LTM phase-offsets in between adjacent wavelengths, and it also agrees with the conclusions given in previous publications [9,12], with Jaron et al. [9] even providing a physical explanation.

To quantify the relationship between LTM phase-offset and the logarithm of the energy, we fit the data presented in Figure 2 with different functions. The first function is a straight line of the form

$$f_1\left(E\right)=a_1·{log}_{10}E+b_1,$$

(4)

where the energy E is given in eV. This fit results in the parameters $a_1=0.02\pm 0.01$, $b_1=0.09\pm 0.05$, and a reduced ${\chi }^2/\mathrm{d}.\mathrm{o}.\mathrm{f}.=21.6/(4-2)=10.8$. This linear fit is plotted as the red line in Figure 2. This shows that the slope a is clearly positive, quantitatively confirming the increasing offset with increasing energy. However, the large ${\chi }^2$ means that the relationship is certainly not that simple.

To account for a possible curvature in the trend, the second function is a parabola of the form

$$f_2\left(E\right)=a_2·{\left({log}_{10}E\right)}^2+b_2·{log}_{10}E+c_2,$$

(5)

and the parameters are estimated to $a_2=0.002\pm 0.001$, $b_2=0.013\pm 0.010$, $c_2=0.031\pm 0.058$, and a reduced ${\chi }^2/\mathrm{d}.\mathrm{o}.\mathrm{f}.=7.8/(4-3)=7.8$. The blue line in Figure 2 shows that this fit represents the data actually quite well despite the rather high ${\chi }^2$. The reason for this is the small number of available data points.

Finally, we note that the data points for the radio, X-ray, and GeV emission seem to lie on a straight line, while the TeV point is apparently offset from this otherwise linear trend. In order to quantify this impression, we restrict the linear fit of Equation (4) to the first three data points. This results in fit parameters $a_1=0.0195\pm 0.0004$, $b_1=0.0828\pm 0.0027$, and a reduced ${\chi }^2/\mathrm{d}.\mathrm{o}.\mathrm{f}.=0.04/(3-2)=0.04$. The corresponding orange line in Figure 2 shows that the three first data points happen to perfectly lie on this line, while the TeV point is clearly offset. The value of the fitted line at the position of the TeV is $f_1\left(5.15\mathrm{TeV}\right)=0.33$, and the TeV measurement data point is with $0.51\pm 0.03$ offset from that value at the $6\sigma$ level. If the phase-ambiguity is resolved to a value of $-0.49$ instead, then this deviation has a significance of even $27\sigma$.

Concerning the function fitting described above it is important to point out that the small number of data points greatly restricts the possibilities here. In the time of writing, we only have four data points to work with. In order to draw more solid conclusions from this, Table 2 and the corresponding plot in Figure 2 have to be populated with more data points. The quantitative findings of systematic LTM phase-offsets as a function of energy will be discussed in a physical scenario in the next section.

4. Discussion

The LTM in the emission from LS I +61°303 has been detected at multiple wavelengths. As evident from the observational facts presented in Section 2, the phase of the modulation pattern is not the same at all wavelengths, but is significantly offset from each other. Putting these phase-offsets into context, as we did in Section 3, reveals that this phase-offset is monotonically increasing as the energy of the observed emission increases, as evident from Figure 2 and the functional fits to these data. A physical scenario that explains the observed long-term phase-offset between the emissions at radio and GeV wavelengths has already been descibed by Jaron et al. [9]. In the following we first show that this scenario can be straightforwardly extended to also include the phase-offsets in the X-ray and TeV energy regimes. However, the divergence of the TeV data point in Figure 2 from the otherwise very simple straight line through the data points at the other three wavelengths deserves an extra comment, given in the second part of this section. Finally, in the third part we will further comment on the exclusion of optical data and describe how their inclusion could help to further discriminate between emission models for this source.

4.1. Plasma Cooling and Opaticity in a Precessing Jet

Figure 3 presents a sketch visualizing how the phase-offsets between the multiple wavelengths fit into a precessing jet scenario. This figure is basically an extension of Figure 7 in Jaron et al. [9], which showed an analogous sketch, but restricted to the radio and GeV emissions. Figure 3 shown here consists of four sub-figures which represent four epochs, ordered chronologically from left to right and labelled $t_0$ to $t_3$. The interval $t_0-t_3$ is short compared to the orbital period $P_1$, with a separation of a few days between adjacent time steps. As a particular example, the results in Jaron et al. [9] would translate to a difference of about five days between $t_1$ and $t_3$. The line of sight is drawn as a dashed arrow which points into the direction of the observer. Different ejections of plasma into a conically expanding jet are marked with different colors as they appear. The ejections propagate into the direction indicated by an arrow and emit radiation at an energy regime as labelled in the same color as the plasma itself. At time $t_0$ an ejection of plasma appears and is marked in blue. It propagates into a direction that is aligned with the line of sight4. The emission from this plasma is in the TeV regime. Because it travels at a significant fraction of the speed of light ($\beta \approx 0.5$, Jaron et al. [33]), the intrinsic emission appears amplified to the observer as a result of maximal Doppler boosting. At $t_1$ the blue ejection has propagated on a ballistic trajectory, following its initial direction and expanding in a conical manner. The energy range of the emission of the blue plasma is now in the GeV regime, as a result of plasma cooling. At the same time, i.e., still at $t_1$, a new ejection of plasma appears, which is now marked in orange. Because of the precession of the jet, this ejection now propagates into a slightly rotated direction compared to the previous blue ejection and is now inclined with respect to the line of sight. The emission of the orange plasma is in the TeV regime, but the Doppler boosting is now decreased compared to the TeV emission of the blue plasma at time $t_0$. In the next step, at time $t_2$, the blue plasma has continued its propagation and expansion, and now emits in the X-ray regime. The orange plasma has done the same in its own direction, now emitting GeV. In addition, there is also a new plasma, this time marked in green, and which is ejected into a yet increased angle as a result of continuous jet precession. Finally, at $t_3$, the multi-wavelengths picture is complete with the original blue plasma having cooled down to emit in the radio. The orange and green ones have continued the process to emit in X-rays and GeV, respectively. In addition, a newly appearing magenta plasma emits in the TeV. All of these emissions at different wavelengths originate from plasmas moving into different directions with respect to the line of sight and hence appear differently Doppler boosted to an observer, who remains at a fixed position.

The process described above shows that the emission at each energy regime (i.e., TeV, GeV, X-rays, and radio) has its peak of amplification at a different moment in time in the interval $t_0$–$t_3$. Following the blue plasma, the maximum Doppler boosting for the TeV emission occurs at $t_0$, for GeV at $t_1$, for X-rays at $t_2$, and for radio at $t_3$. Looking at the temporal evolution of the direction of one wavelength alone, e.g., TeV, one can trace how the amplification at this energy changes with time. We refer to this variability in the Doppler boosting of the intrinsic emission as the precession profile at a certain wavelength. The crucial point to emphasize now is that the peak of maximal Doppler boosting occurs earlier in the precession profile with increasing energy. This is the complete opposite as for the LTM pattern, where the peak of the higher energy occurs later. In the following we will see how this can be understood by investigating mathematically the interference of two periodic processes.

Mathematically, the precession profile can be expressed, to first order, as a cosine function of the form

$$f_{\mathrm{prec}}\left(t\right)=cos2\pi \left(\frac{t-T_0}{P_2}-{\varphi }_{\mathrm{mp}}\right),$$

(6)

where $P_2$ is the precession period of the jet and $T_0$ is the conventional time-origin5 for LS I +61°303 ([9] and references therein). The phase ${\varphi }_{\mathrm{mp}}={\varphi }_{\mathrm{mp}}\left(E\right)$ of maximum (m) Doppler boosting due to precession (p) is a function of the energy E of the observed emission. This has been empirically proven for the radio and GeV emission by Jaron et al. [9], and our scenario here is based on the assumption that this is also true for the X-ray and TeV emission. We follow the convention that the phase ${\varphi }_{\mathrm{mp}}$ is a number between 0 and 1, as for the orbital phase.

Several authors showed that for the eccentric orbit of LS I +61°303 two accretion peaks are predicted [28,29,30,31]. As pointed out by Taylor et al. [29], the Bondi [50] accretion rate $\dot{M}\propto \rho /v_{\mathrm{rel}}^3$ is proportional to the density $\rho$ of the stellar wind at the position of the accretor, and inversly proportional to the cube of the relative velocity $v_{\mathrm{rel}}$ between the accretor and the wind. As a consequence, one accretion peak is at periastron, where the density $\rho$ is maximal. If the eccentricity e of the orbit is sufficienty high, there is a second accretion peak at a later orbital phase, when the smaller velocity $v_{\mathrm{rel}}$ compensates for the lower density $\rho$. As Taylor et al. [29] showed (see their Figure 7), this second accretion peak develops at an eccentricity $e=0.4$ and is well pronounced for $e=0.6$. The eccentricity of the orbit of LS I +61°303 has been determined to be $e=0.72\pm 0.15$ by Casares et al. [1], so two peaks along the orbit are clearly expected. It is, however, only the emission from the second, apastron, accretion peak which is subject to the LTM. Jaron et al. [33] explain this observational fact with a physical model in which the ejection at periastron only reaches velocities of $\beta \approx 0.01$, which is not fast enough to give rise to appreciable Doppler boosting. For this reason, we restrict the accretion profile here to apastron orbital phases. Mathematically this leads to an expression analogous to the one given in Equation (6),

$$f_{\mathrm{orb}}\left(t\right)=cos2\pi \left(\frac{t-T_0}{P_1}-{\varphi }_{\mathrm{mo}}\right),$$

(7)

where $P_1$ is the well-known orbital period of the system and $T_0$ is the same time-origin as for the precession profile above. The phase ${\varphi }_{\mathrm{mo}}$ of the maximum (m) of the orbital (o) modulation of the accretion rate does not depend on the energy, but is a constant ${\varphi }_{\mathrm{mo}}=\mathrm{const}$, extrapolating the results from Jaron et al. [9] to the X-rays and TeV. (Following Jaron et al. [33] this is ${\varphi }_{\mathrm{mo}}=0.58$, see their Table 1, called ${\mathsf{\Phi }}_0$ for injection II there.) We point out that the actual value of ${\varphi }_{\mathrm{mo}}$ is not important for the following explanation, so we are going to neglect this term from here on for the sake of clarity of the formulas.

The emission from LS I +61°303, as received by an observer on the Earth, is modulated by two mechanism. The first mechanism is the modulation of the accretion rate, and the second one is the amplification of the intrinsic emission of the jet resulting from periodic changes in the Dopper boosting due to precession of the relativistic jet. These two mechanisms interfere with each other, giving rise to an interference pattern that has the form of a beating because of the proximity of the two intrinsic periods $P_1$ and $P_2$. This beating is identical to the LTM, as first discovered by Massi & Jaron [24]. This hypothesis has been confirmed by physical modelling of the radio [32] and GeV [33] data, resulting from long-term monitoring of the source. Here we investigate the effect of an energy-dependent phase-offset of the precession profile in a way that is analogous to the presentation in the appendix of Jaron et al. [9]. The interference between orbit and precession can be studied by considering the sum of the two,

$$\begin{array}{c}\hfill f_{\mathrm{beat}}\left(t\right)=f_{\mathrm{orb}}\left(t\right)+f_{\mathrm{prec}}\left(t\right)=cos2\pi \left(\frac{t-T_0}{P_1}\right)+cos2\pi \left(\frac{t-T_0}{P_2}-{\varphi }_{\mathrm{mp}}\right)\\ \hfill \propto cos2\pi \left(\frac{t-T_0}{P_{\mathrm{avg}}}-\frac{{\varphi }_{\mathrm{mp}}}{2}\right)cos2\pi \left(\frac{t-T_0}{2P_{\mathrm{beat}}}+\frac{{\varphi }_{\mathrm{mp}}}{2}\right),\end{array}$$

(8)

where the first cosine term oscillates at a period which is the inverse of the average of the frequencies $1/P_1$ and $1/P_2$, and the second cosine term is slowly oscillating at the twice the beat period $P_{\mathrm{beat}}=1/(P_1^{-1}-P_2^{-1})$. The envelope of this interference pattern has a period of $P_{\mathrm{beat}}$ and is identical to the LTM of LS I +61°303, i.e., $P_{\mathrm{beat}}=P_{\mathrm{long}}$, as explained in Massi & Jaron [24]. The important point here is that the precession phase-offset ${\varphi }_{\mathrm{mp}}$ propagates into the interference pattern with the opposite sign. This means that emission which has its maximum earlier in the precession profile (corresponding to negative values of ${\varphi }_{\mathrm{mp}}$) has its maximum in the LTM pattern shifted to a later time (or phase). This is exactly what we observe for the LTM phase-offsets presented in Figure 2. Remembering that the LTM pattern corresponds to the envelope of the interference pattern, the phase shift of the LTM occurs exactly by the same amount as for the precession profile but in the opposite direction, as outlined by Jaron et al. [9].

In this way, we can conclude that the sequence of the emission of the energy regimes TeV, GeV, X-ray, and radio occurs in that particular order in the precessing jet scenario. However, it is important to point out that this reflects a temporal order first of all. Because the jet is filled with plasma at its base and material propagates in the jet direction from that point on, we can also conclude that this temporal order is identical to the spatial sequence of emission regions along the jet. However, to draw conclusions about the exact locations of these regions is not possible in this scenario and with the data available. The reason for this is that there is a degeneracy of the problem in terms of distance from the base and velocity of the plasma along the jet. It is well possible that the processes probed by the emission observed for this study occur in a part of the jet where acceleration of the plasma still plays an important role. This has also been pointed out by Jaron et al. [33] to explain the slow velocity of a jet that is launched at periastron of the system, and which is not affected by the LTM because the bulk velocity does not reach relativistic speeds.

The deviation of the TeV data point in Figure 2, however, is remarkable. In the scenario presented until here there are three possible reasons for this relatively large offset:

• The distance between the location of the TeV and GeV emission is considerably larger than between the other wavelengths.
• It takes the plasma longer to travel from the TeV region to the GeV region than in between the other regions because the velocity is smaller closer to the base of the jet.
• It takes the plasma longer to lose its energy from TeV down to GeV than from that point onwards.

In the end, a mixture of all of these three options would also be possible. However, there is also an alternative explanation that we will sketch in the next subsection.

4.2. Are the TeV Photons Produced by a Different Mechanism?

The fact that the LTM phase-offset values for the radio, X-ray, and GeV emission are perfectly fit by a straight line as a function of $log\left(E\right)$, while the TeV data point is offset from that simple trend by $6\sigma$ (or even $27\sigma$ with alternative ambiguity resolution for that data point, as pointed out in Section 3) justifies the consideration of a different or at least modified explanation for the TeV emission. Furthermore, the fact that the absolute phase-offset of ∼0.5 with respect to the radio means that the LTM pattern of the TeV is in perfect anti-phase with the radio, and deserves special attention. That the VHE emission is affected by the LTM, i.e., the beating between orbit and precession, implies that these photons are produced in a mechanism that is subject to both accretion and periodic changes in the relativistic beaming. One possible explanation could be that during magnetic reconnection events, which can occur both in the accretion disk [51] or the jet itself [52,53], plasmoids are injected into the steady jet and collide with the slower material there. This results in shocks in which the magnetic field lines are compressed into one direction and force the plasmoids to move into a different direction than the steady jet. Consequently, the Doppler factor would be systematically different for the TeV emission than for the other wavelengths. In this scenario, a localization of the TeV emission relative to the other emissions, as depicted in Figure 3 is not possible any more. It is worth of note that the short-term variability in LS I +61°303, observed in the radio [9] and X-ray [6] emission, mentioned in the beginning of Section 1 of this article can be explained in the same scenario of magnetic reconnection.

4.3. What Is the Origin of Optical Emission in LS I +61°303?

As mentioned in Section 2.2, there is also observational evidence that the equivalent width of the H$\alpha$ line in LS I +61°303 is subject to the LTM. As explained there, we excluded these optical observations from our analysis. In the context of the physical scenarios presented here, the origin of the H$\alpha$ line is more complicated. Moreover, the possibility of the emission line to originate from the circumstellar disk of the Be star, Fender et al. [54] point out that the accretion disk of X-ray binaries is another source of H$\alpha$ line emission. Furthermore, by definition EW(H$\alpha$) is a function of both the H$\alpha$ luminosity itself as well as the continuum emission. In addition, especially in the low/hard state of X-ray binaries, the jet can contribute significantly in the optical band [55]. The continuum optical jet emission, in this scenario, is subject to variable Doppler boosting in the same way as the other wavelengths emitted from the precessing jet, and consequently the resulting LTM of the continuum emission would propagate into the modulation of EW(H$\alpha$). Further observations at optical wavelengths and investigation of the origin of the optical emission are needed to clarify these issues.

At present the largest void in Figure 2 is between radio and X-rays. A reliable determination of the LTM phase-offset at optical wavelengths would provide data just in the center of this gap. For this reason alone, optical data would be a natural first choice to be added to a future extension of the analysis presented here. However, before that, existing databases of optical observations will have to be revisited, and continuum and line emission will have to be carefully disentangled. In the scenario suggested in Section 4.1, the optical continuum emission is expected from a steady jet originating from a region between those of the radio and X-ray emission. In this case the LTM phase-offset of the continuum optical emission alone would be expected to be ∼0.1, interpolating the linear trend that extends from radio to GeV (orange line in Figure 2). Concerning the H$\alpha$ emission line, sources of this line can be both the circum-stellar disk of the Be star and the accretion disk, so a variability in this emission line probes variability in both of these regions. Furthermore, in a feasibility study presented by Massi & Zimmermann [56], the precessing jet is the result of a precessing accretion disk. Hence, the H$\alpha$ emission from the accretion disk would be subject to the same intrinsic periodicites, i.e., accretion rate modulated by the orbital period $P_1$, and variable Doppler boosting modulated by the precession period $P_2$, resulting in the beating that we identify with the LTM. In this scenario the LTM phase-offset of this H$\alpha$ emission could be approximately anti-correlated with the TeV emission, if the TeV emission originates at the base of the steady jet, as proposed in Section 4.1. In this way, optical observations have the potential to sigificantly contribute in emission model discrimination for LS I +61°303.

5. Conclusions

In this article, we investigated the long-term modulation of the $\gamma$-ray-loud high-mass X-ray binary LS I +61°303 at multiple wavelengths. The periodic modulation of ∼4.6 years is consistently detected at radio, X-rays, GeV, and TeV energies by analysis of the past four decades of observations at these wavelengths. Subject of this article was the investigation of systematic phase-offsets of the modulation pattern between these wavelengths and the explanation in a physical scenario. These are our conclusions:

• By re-analyzing archived TeV data [14] we determined that LTM pattern at TeV is offset from the pattern at radio by $0.51\pm 0.03$. We point out that, while this value fits into the monotonically incrasing trend with the other wavelengths, this also means perfect anti-correlation with the radio emission in terms of the LTM pattern.
• The LTM of LS I +61°303 is established at radio, X-ray, GeV, and TeV wavelengths by long-term monitoring of the source. There is a systematic trend of the modulation pattern being increasingly offset from the radio pattern as the energy increases. Emission at higher energy is lagging emission at lower energy in a strictly monotonically increasing manner (Figure 2).
• We extended the physical scenario first introduced by Jaron et al. [9] to X-rays and TeV. In this scenario, the emission regions are located closer to the base of the jet as the energy increases (Figure 3). Because the LTM is the result of interference between orbit and precession, an earlier peak in the precession profile propagates to a later peak in the LTM pattern with increasing energy, as observed. We also note that while the radio, X-ray, and GeV emissions are fitted by a very simple trend, the TeV emission appears to be significantly offset from that. This is a challenge for this scenario at TeV energies.
• We consider an alternative scenario in which the TeV photons are produced in shocks resulting from the injection of plasmoids during magnetic reconnection events. It is interesting to note that this would imply a connection between the large phase-offset of the LTM at TeV energies and the hour time-scale phenomenon of quasi periodic oscillations at radio and X-rays. This would explain two phenomena, occurring at different energy bands and on time-scales of different orders of magnitude (hours vs. years) in the same coherent picture. This should be subject of future investigations.

Further long-term monitoring of LS I +61°303 at more wavelengths is needed in order to populate Figure 2 with more data points. This would enable a better quantification of the relationship between phase-offset and energy. At the moment functional fits are complicated by the small number of available data points. In addition, a better sampling at X-rays and especially TeV would allow quantitatively studying the precession profile at these wavelengths and to verify whether a shift in this pattern is indeed present. This would also help to disentangle the two scenarios for the TeV emissions and answer the question whether these photons are produced in the steady or transient jet.