EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 554 (June 2013) A&A, 554 (2013) A105 Full HTML
Free Access
Issue
A&A
Volume 554, June 2013
Article Number A105
Number of page(s) 8
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/201219685
Published online 12 June 2013

© ESO, 2013

1. Introduction

The TeV emitting source LS I +61°303 has radio characteristics that make it unique not only among the small number of gamma-ray emitting systems given in the X-ray binary class, a class of binary systems where a neutron star or a black hole is orbiting around a normal star, but also among the larger group of the radio emitting X-ray binary systems (Fender et al. 1997; Massi 2005; Mirabel 2012). The peak flux density of the radio outburst, which is related to the orbital period of 26.4960 ± 0.0028 d, exibits a modulation of 1667 ± 8 d (Gregory 2002). Some double peaked outbursts when observed at two frequencies show different spectral characteristics. There is a first outburst with a flat/inverted spectrum and a second optically thin outburst associated with different conditions, as indicated by its high amplitude, the spectral index, and the Hα emission line measurements (Massi & Kaufman Bernadó 2009; Grundstrom et al. 2007). The complex spectral sequence found in LS I +61°303 finds a natural explanation in the disk-jet coupling model for microquasars: first, there is a continuous outflow with a flat or inverted spectrum, then an event triggers a shock in this slow optically thick outflow (Fender et al. 2004), and the growing shock creates the optically thin outburst (Valtaoja et al. 1992; Hannikainen et al. 2006). One of the characteristics that make LS I +61°303 unique among the other radio emitting X-ray binary systems is that this spectral evolution, between inverted and optically thin spectra, may occur twice during the orbital period (Massi & Kaufman Bernadó 2009). This agrees well with the Bondi & Hoyle (1944) accretion in an eccentric orbit (as in LS I +61°303) predicting two events along the orbit as shown for LS I +61°303 by several authors (Taylor et al. 1992; Marti & Paredes 1995; Bosch-Ramon et al. 2006; Romero et al. 2007).

The binary system LS I +61°303, for which the nature of the compact object has not yet been established (i.e., a black hole or a neutron star), shares the remarkable property of slow radio quasi-periodic oscillations (Peracaula et al. 1997, 84 min) with the two black hole microquasars V404 Cyg (Han & Hjellming 1992, 20–120 min) and GRS 1915+105 (Pooley & Fender 1997; Rodriguez & Mirabel 1997, 20–40 min).

The radio morphology of the system also shows unique characteristics. The resolved extended structure changes position angle (i.e., angle of the projection of the jet onto the plane of the sky, measured from north through east) with the surprising large variation of 60° in only one day (Massi et al. 2004; Dhawan et al. 2006). Moreover, the jet is sometimes one-sided and at other times two-sided. Because of both of these variations in position angle and morphology, the hypothesis that LS I +61°303 might be a precessing microquasar was brought forth (Massi et al. 2004). The one-sidedness of jets is usually attributed to relativistic bulk motion along a relatively small angle to the line of sight, which leads to Doppler boosting of the jet and deboosting of the counterjet emission (Urry & Padovani 1995). A variation of that angle due to precession would cause variable Doppler deboosting of the counter jet, making it appear for larger angles (double-sided jet) and disappear for smaller values of the angle to the line of sight (one-sided jet, blazar like).

In 2006, VLBA observations by Dhawan et al. (2006) measured the same large rotation of 60°/day in their images as Massi et al. (2004). Some of the VLBA images, showing again a one-sided structure were, however, interpreted by Dhawan et al. (2006) as a cometary tail pointed away from the companion Be star, that is in favor of a pulsar model rather than of the precessing microquasar model. The reanalysis of this VLBA data set and the resulting higher dynamic range of the self-calibrated maps has actually revealed a double-sided structure in several images (Massi et al. 2012). Before we illustrate how the results in Massi et al. (2012) brought us to the present investigation on the precessional period, let us consider two important points in the two following paragraphs, the first concerning self-calibration and the second the pulsar model.

Self-calibration of interferometric data is a well-established technique (Cornwell & Fomalont 1999). It may fail at SNR < 4 (Martí-Vidal & Marcaide 2008) or may create spurious symmetrization for unbalanced closure phase triangles (resulting when a very displaced telescope is added to the array) (Massi & Aaron 1999). None of these two cases apply to the used set of only VLBA data, where double-sided structures at 8–16σ are present in 6 out the 12 images. Concerning the effects on an image of strong variations of flux density during observations (Stewart et al. 2011), the source LS I +61°303 shows a low radio flux density at all orbital phases, apart from the maximum of the long-term modulation. In these epochs a large outburst lasting few days occurs around apastron. This means that during the maximum of the long-term modulation one might expect a reduction of the dynamic range of the produced maps around apastron, i.e., weak features will be lost. This is not the case for the VLBA observations of Dhawan et al. (2006) performed toward the minimum of the long-term modulation.

We are therefore faced in LS I +61°303 with a changing structure from a double-sided structure to a one-sided structure. Are such variations possible also in the pulsar scenario? Simulations (Moldón et al. 2012) show that the emission from the cometary tail of a pulsar for a particular orientation and inclination of the orbit, after almost one orbital cycle and at a particular orbital phase may look like a double-sided nebula at a fixed position angle. In these conditions, the fading and expanding last part of the cometary tail may appear detached from the brightest part of the cometary tail which is closer to the orbit. This possibility is clearly ruled out for LS I +61°303, where along one orbital cycle the VLBA images show a double-sided jet at several different orbital phases and at different position angles (Massi et al. 2012).

We return thus to the scenario of a precessing microquasar, where the two-sided jet suffers from variable Doppler boosting because the precession continuously changes the angle between jet and line of sight (Massi et al. 2004, 2012). Deriving the precessional period from the available radio images is not straightforward, because the images reflect the variation of the projected angle on the sky plane and, therefore, a combination of the ejection angle, inclination, and angle of the precession cone. The astrometry provides less biased results (Dhawan et al. 2006; Massi et al. 2012). The astrometry of the VLBA observations indicates that the peak in consecutive images describes a well-defined ellipse, six to seven times larger than the orbit, with a period of about 27–28 d (Massi et al. 2012). Assuming for the microquasar in LS I +61°303 the same core-shift effect typically observed in blazars, the peak of each image would then correspond to that part of the jet where the emission becomes optically thick at the observing frequency (Kovalev et al. 2008). Based on this assumption, Massi et al. (2012) interpreted the ellipse as the possible cross-section of the precession cone of the jet at the distance where the emission at 8.4 GHz becomes optically thick. The determined time span of 27–28 d to complete the ellipse is a first estimate of the precession period.

The most likely cause for precession of an accretion disk of a compact object is an assymetric supernova explosion of the progenitor. As a result the compact object could be tilted (Fragile et al. 2007). In this case either the accretion disk is coplanar with the compact object and, therefore, subject to the gravitational torque of the Be star or, instead, the accretion disk is coplanar with the orbit but tilted with respect to the compact object which induces, in the context of general relativity, Lense-Thirring precession if the compact object rotates (Massi & Zimmermann 2010). A deep investigation of these or other mechanisms of precession requires the knowledge of the precession parameters, such as the period of precession and the angle of the precession cone.

In this paper we present a timing analysis of 6.7 years of Green Bank Interferometer (GBI) radio data aimed at a more accurate determination of the precession period. Sections 2.1 and 2.2 present the determined period, called P2. The section illustrates the case that, while the aim of our research was reached, i.e., we obtained a more precise value of P2, we were presented with an additional unexpected result: the beating between the orbital period P1 and the precessional period P2 gives rise to a new period, Paverage = 2/(ν1 + ν2) modulated by 1/(ν1 − ν2) = 1667 d, i.e., the long-term modulation. Is the found Paverage the periodicity of the observed radio outburst? Indeed, whereas in the literature P1 is generally referred to as the period of the radio outbursts, it is also well known that there are differences between the observed and predicted (for P = P1) outburst times, and that these timing residuals have vs. time a sawtooth pattern (Sect. 2.3). In Sect. 2.4 we present two important results. First, we demonstrate mathematically that indeed the sawtooth function adjusts P1 to Paverage. Then we show that the GBI data folded with Paverage present an offset of 13 d at the minimum of the long-term modulation equal to the sawtooth function and equal to that predicted by the beat between P1 and P2. In the same section we also present the corresponding physical scenario. In Sect. 3 we discuss the implications of our results for the observed periodicity in the equivalent width of the Hα emission line in LS I +61°303. In Sect. 4 we present our conclusions.

thumbnail Fig. 1

Periodograms of 8.3 GHz (filled squares) and 2.2 GHz (circles) data. a) Output of Lomb-Scargle method. The Lomb-Scargle analysis gives on the y-axis the significance of the frequency. Two frequencies are found at ν1 = 0.03775 d-1 (P1 = 26.49 d) and at ν2 = 0.03715 d-1 (P2 = 26.92 d). b) Output of Lomb-Scargle method for data ≥4σ. c) Output of PDM. The most likely period yields the minimum dispersion and appears as a minimum in the PDM curve, i.e., specular to the maximum in the Lomb-Scargle plot. d) Output of PDM for data ≥4σ.

2. Data analysis and results

We analyzed here 6.7 years of the NASA/NRAO GBI LS I +61°303 database at 2.2 GHz and 8.3 GHz. The database covers three periods: 49 379.975–50 174.710 MJD, 50 410.044–51 664.879 MJD, and 51 798.333–51 823.441 MJD. The samples, flux densities at both frequencies with their corresponding errors, in each of the three time intervals are almost continuous with an average of eight observations per day. In order to search for possible periodicities we used the Lomb-Scargle method, which is very efficient on irregularly sampled data (Lomb 1976; Scargle 1982). We used the algorithms of the UK Starlink software package, PERIOD1. For both data sets at 2.2 GHz and 8.3 GHz (Fig. 1) we obtained the same result: two periods P1 = 26.49 ± 0.07 d (ν1 = 0.03775 ± 0.00010 d-1) and P2 = 26.92 ± 0.07 d (ν2 = 0.03715 ± 0.00010 d-1). The statistical significance of a period is calculated in PERIOD following the method of Fisher randomization as outlined in Linnel Nemec & Nemec (1985). The advantage of using a Monte-Carlo- or randomization-test is that it is distribution-free and that it is not constrained by any specific noise models (Poisson, Gaussian, etc.). The fundamental assumption is that if there is no periodic signal in the time series data, then the measured values are independent of their observation times and are likely to have occurred in any other order. One thousand randomized time series are formed and the periodograms calculated. The proportion of permutations that give a peak power higher than that of the original time series would then provide an estimate of p, the probability that, for a given frequency window there is no periodic component present in the data with this period. A derived period is defined as significant for p < 0.01, and a marginally significant one for 0.01 < p < 0.10 (Linnel Nemec & Nemec 1985). For both periods P1 = 26.49 ± 0.07 d (frequency window 0.0374−0.0379 d-1) and P2 = 26.92 ± 0.07 d (frequency window 0.0369−0.0374 d-1) and for both data sets at 8.3 GHz and 2.2 GHz we obtained 0.00 < p < 0.01.

2.1. Relative importance between the two periods and previous observations

thumbnail Fig. 2

Long-term modulation and period analysis. a) 8.3 GHz GBI radio data averaged over 3 d. b) Lomb-Scargle analysis results: two frequencies at ν1 = 0.03775 d-1 (P1 = 26.49 d) and ν2 = 0.03715 d-1 (P2 = 26.92 d). The small window shows the peak at 1/1667 d-1 present in all periodograms. c) Sum of two sinusoidal functions at 26.49 d and 26.92 d, with an amplitude ratio 1/0.7. d) Lomb-Scargle analysis results: two frequencies at 126.49d-1\hbox{$\frac{1}{26.49}~{\rm d^{-1}}$} and 126.92d-1\hbox{$\frac{1}{26.92}~{\rm d^{-1}}$}. The significance of the two frequencies becomes identical in the periodogram only for an amplitude ratio 1/1. e) Long-term modulation (1667 d) of a 26.49 d periodic outburst. f) Lomb-Scargle analysis results: one frequency at 126.49d-1\hbox{${1\over 26.49}~{\rm d}^{-1}$}. The small window to the left shows the peak at Paverage present in the periodogram of a simulation of long-term modulation (1667 d) of a 26.70 d periodic outburst. g) Sine wave of periodicity P1, modulated by a sine wave of periodicity 1667 d and corrected by a sawtooth function. h) Lomb-Scargle analysis results: two frequencies at 126.49d-1\hbox{$\frac{1}{26.49}~{\rm d^{-1}}$} and 126.92d-1\hbox{$\frac{1}{26.92}~{\rm d^{-1}}$} as in Figs. 2b and d.

Figure 1a shows the results of the Lomb-Scargle analysis for the data at 8.3 GHz and 2.2 GHz. There P1 is dominating over P2 for a factor of 1.8 at 2.2 GHz, and a factor of 1.5 at 8.3 GHz. Figure 1b shows the results of the Lomb-Scargle analysis, if only data with flux density ≥4σ are used. In this case the two periods have a more comparable significance, i.e., there is a factor of 1.4 at 2.2 GHz, and a factor of 1.2 at 8.3 GHz. We compare the Lomb-Scargle results with those obtained with the phase dispersion minimization (PDM) method (Stellingwerf 1978). The results of the PDM analysis on the whole data sets are shown in Fig. 1c; the results on data ≥4σ are shown in Fig. 1d. The results of the PDM analysis agree very well with those of Lomb-Scargle: there is a different significance of the two periods when data with low signal-to-noise-ratio (snr) are present in the analyzed data set, i.e., P1 dominates.

This could be the explanation why in the past the second period P2 was unseen. Taylor & Gregory (1982) found a period of 26.52 ± 0.04 d in their radio data set; in 1984 they wrote that part of the previous measurements were taken with the source in a weak state and repeated the analysis using new data and also including the old ones, obtaining the value of 26.496 ± 0.008 d (Taylor & Gregory 1984). In 1997 Ray et al. reported new observations and gave a period of 26.69 ± 0.02 d, i.e., coincident with our average Paverage=2ν1+ν2=20.03775+0.03715d=26.70± 0.05d\hbox{$P_{\rm average}= \frac{2}{\nu_1+\nu_2} = \frac{2}{0.03775 + 0.03715}~{\rm d}=26.70~\pm\,0.05~{\rm d}$}, as discussed in Sect. 2.2. These observations (Ray et al. 1997) are those of our first GBI interval, i.e., 49 379.975−50 174.710 MJD. In Fig. 2a one sees that these data sample only the interval of maximum activity. The difference between the results of Taylor & Gregory (1984) and of Ray et al. (1997) is puzzling for both groups. Ray et al. (1997) discuss how their best estimate of the period is significantly different (9σ) from the 26.496 ± 0.08 d value of Taylor & Gregory (1984). Gregory et al. (1999), faced with the difference between their value and the Ray et al. (1997) results, discuss how unlikely a sudden change in period would be. In the light of our present result we see that it is not a sudden change in period but the presence of two periods that in the Ray et al. data have comparable significance. As noted above, the value of Ray et al. (1997) 26.69 ± 0.02 d corresponds to our Paverage, as discussed in the next section.

2.2. Beating: long-term modulation and Paverage

The two frequencies ν1 = 0.03775 d-1 (P1 = 26.49 d) and ν2 = 0.03715 d-1 (P2 = 26.92 d) are only slightly different. This produces a beating, i.e., a new frequency is formed νaverage=ν1+ν22\hbox{$ \nu_{\rm average}={{\nu_1+\nu_2}\over 2}$}, modulated with νbeat = ν1 − ν2. For the sum of two sine functions the following identity holds sin(2πν1t)+sin(2πν2t)=2cos(2πν1−ν22t)sin(2πν1+ν22t),\begin{equation} \label{eq:beat1} \sin\left(2\pi \nu_1 t\right) + \sin\left(2\pi \nu_2 t\right) = 2\cos \left(2\pi{{\nu_1-\nu_2}\over 2} t\right)\sin\left(2\pi {{\nu_1+\nu_2}\over 2}t\right), \end{equation}(1)where the beat frequency (or frequency of the envelope) νbeat = ν1 − ν2 is twice the frequency of the cosine term. In our case, the term 1ν1−ν2=10.03775−0.03715d\hbox{${1\over \nu_1-\nu_2}=\frac{1}{0.03775 - 0.03715}~{\rm d}$} is equal to 1667 ± 393 d2.

Figure 2c shows the sum of two sine functions with different amplitudes fb(t)=sin(2πν1t)+asin(2πν2t)=2acos(2πν1−ν22t)sin(2πν1+ν22t)+(1−a)sin(2πν1t),\begin{eqnarray} f_{\rm b}(t) & = & \sin\left(2\pi \nu_1 t\right) + a\sin\left(2\pi \nu_2 t\right)\nonumber\\ & = & 2a\cos \left(2\pi{{\nu_1-\nu_2}\over 2} t\right)\sin\left(2\pi {{\nu_1+\nu_2}\over 2}t\right)\nonumber\\ \label{eq:beat2} & & + \left(1 -a\right)\sin\left(2\pi\nu_1t\right), \end{eqnarray}(2)with a = 0.7, whereas Fig. 2e shows the function fm(t)=(1+bsin(2πνmt))sin(2πν1t),\begin{equation} \label{eq:mod} f_{\rm m}(t) = (1+b\sin(2\pi\nu_{\rm m}t))\sin(2\pi\nu_1t), \end{equation}(3)with b = 0.7, ν1=126.49d-1\hbox{$\nu_1 = \frac{1}{26.49}~{\rm d^{-1}}$}, and νm=11667d-1\hbox{$\nu_{\rm m} = \frac{1}{1667}~{\rm d^{-1}}$}. As one can see, both Eqs. (2) and (3) are able to reproduce the long-term modulation; however, the periodograms are rather different. The periodogram of Eq. (2), shown in Fig. 2d agrees well with the periodogram of the GBI data of Fig. 2b. On the contrary, as one can see in Fig. 2f, in the periodogram of Eq. (3) only P1 = 26.49 d is present in the frequency range 0.036−0.039 d-1.

2.3. Sawtooth function

Gregory (1999) demonstrated the existence of a long-term modulation of the peak outburst flux. Gregory & Neish (2002) indicated that the modulation in radio properties may stem from periodic ejections of a shell (density enhancement) of gas in the equatorial disk of the Be star.

The long-term modulation is also present in the timing residuals of the outbursts, i.e., the difference between observed and predicted (for P = P1) outburst time (Gregory et al. 1999). The observational result is that timing residuals show a surprising sawtooth pattern, i.e., with a gradual rise from 0 to about 6 d, but then a rapid fall to a large negative value of about –7 d. The sawtooth function is shown in Figs. 2 and 8a of Gregory et al. (1999) and here in Fig. 3a. In detail the trend is as follows: observed and predicted outburst times coincide at the peak of the long-term modulation (i.e., at the radio maximum) resulting in a timing residual equal to τ = 0; the timing residual grows linearly with time and at the minimum of the long-term modulation reaches a maximum of about 6 d where it sharply switches to about –7 d. Also surprising is that the transition τ ≃ 6 d to τ ≃ − 7 d is not related to a strong change in flux density; it occurs at a time when the amplitude is low, i.e., at the minimum of the long-term modulation (Figs. 6a and 7b of Gregory et al. 1999). After this transition, τ starts to grow linearly with time reaching the value τ = 0 only after about 800 d at the new maximum of the long-term modulation.

We performed a test by using the sawtooth function to correct the model of Eq. (3) (P1 modulated by 1667 d) and to verify if the corrected model is able to reproduce the observed spectrum (i.e., P1 and P2). First we have to define the sawtooth function. The slope of the sawtooth pattern in Fig. 8a of Gregory et al. (1999), is about 0.008; we generate therefore the sawtooth function (Fig. 3a) with a period of 1667 d as3τ(t)=0.008 fmod (t,1667)\begin{equation} \tau(t) = 0.008\,\mathrm {fmod}\,(t, 1667) \end{equation}(4)and include it in Eq. (3), which becomes: fm(t)=(1+bsin(2πνmt))sin(2πν1(t−τ(t))).\begin{equation} \label{eq:mod+st} f_{\rm m}(t) = (1 + b\sin(2\pi\nu_{\rm m}t))\sin(2\pi\nu_1(t - \tau(t))). \end{equation}(5)To study the effect of the sawtooth function without any bias resulting from large holes in the sampling, we performed both simulations with original sampling and with regular sampling. We obtained the same results and here we show those with regular sampling.

The resulting periodogram of the simulated data (Fig. 2g) is shown in Fig. 2h. One sees that the model with the single periodicity P1, modulated by 1667 d, once corrected by the sawtooth function, is able to reproduce the results of our spectral analysis that is the two periods P1 and P2. In the next section we show that when one directly uses the two found periods P1 and P2, the sawtooth function is naturally explained.

2.4. Period of the observed outburst and Paverage

thumbnail Fig. 3

a) Sawtooth function based on timing residuals observed by Gregory et al. (1999). b) Sawtooth function based on the loss of synchronization of the sine functions f(P1) and f(P2) of Fig. 3c. c) Sine waves with periods P1 (black line) and P2 (red line) starting with synchronized peaks at t = 0. d) Sum of the sine waves. The first five bars are at a regular interval equal to Paverage, then the next bar follows at a distance of only 13.25 d (P1/2), followed by four bars again at a regular interval equal to Paverage as before. The two central bars are symmetric with respect to a peak of f(P2) (Fig. 3c) that is nearly equidistant from the preceeding peak and the delayed peak of f(P1) (see Appendix).

The sawtooth function results from the comparison between observed and predicted (for P = P1) outburst time. Here we will show first analytically (Eq. (6)) and then with the GBI observations that the observed outburst occurs at P = Paverage, i.e., at 1/νaverage of Eq. (1).

Analytically, adjusting P1 = 26.49 d (ν1 = 0.03775 d-1) by the given sawtooth function τ(t) to ν1(t−τ(t))=ν1t(1−0.008)=(0.03775×0.992)t=0.03745t\begin{equation} \nu_1(t - \tau(t)) = \nu_1 t(1 -0.008) = (0.03775\times 0.992)~t = 0.03745~t \end{equation}(6)gives as a result Paverage=10.03745d=26.70d\hbox{$P_{\rm average} = \frac{1}{0.03745}~{\rm d} = 26.70~{\rm d}$} as in Sect. 2.1, (with variations between 26.65−26.70 d for slopes between 0.006−0.008). This is an important result. It implies that the timing residuals between predicted (at P1) and observed outbursts are equal to the residuals between predicted (at P1) and Paverage.

thumbnail Fig. 4

a) Radio light curves (data averaged over 1 d) vs. Φorbit, with Φorbit related to (t−t0)P1\hbox{${(t-t_0)\over P_1}$}, with P1 = 26.49 d and t0 = JD2 443 366.775 (Gregory 2002), before (black) and after (green) (2 400 000.5 + 50 841)JD (see Appendix). b) Same as a) but for Φprecession, related to P2. c) Radio light curves vs. Φaverage for Paverage = 26.70 d. The small window shows the simulated data (Eq. (2)) identically folded. d) Radio light curve vs. Φaverage with t0 changed by Δt (Δt = 0 before (2 400 000.5 + 50 841)JD and Δt = 13.25 d after it, see Appendix). e) Sketch of the precessing jet in LS I +61°303 (out of scale).

We therefore ascertained that the observed outburst has periodicity Paverage. However, if we use Eq. (3) and we simply substitute P1 by Paverage, the periodogram fails to reproduce the observed periodogram with P1 and P2 and just shows Paverage (small window on the left in Fig. 2f). Is Paverage only an apparent periodicity, just produced by P1 and P2?

Let us fold the GBI data on the periods P1, P2, and Paverage (Figs. 4a–c). The data folded with the orbital period P1 show the broad cluster that is well known in the literature (e.g., Fig. 2c in Massi & Kaufman 2009), where flares above 100 mJy occur from about phase 0.4 to about phase 0.9. The broadening is due to the differences between the observed and predicted (by P1) outburst time that also causes the sawtooth pattern. Figure 4b, shows for the first time the data folded with the precessional period P2. The cluster of the large flares is also evident and superimposed to scattered smaller flares. If now we fold the data with Paverage one would expect, since it is the average of P1 and P2, a clustering rather similar to that in Figs. 4a and b. The result, shown in Fig. 4c, is completely different. First of all, there are two clusters. Second, each one of the two clusters is not as broad as those with P1 and P2, i.e., the clustering is better. Where does the double clustering come from? We used the color green for data after 50 841 MJD, located in the minimum. A harmonic might theoretically give rise to two possible clusters, but in this case green and black points had to be present in both clusters. On the contrary, all data before the minimum, i.e. the black points, cluster at one phase and all points after that, i.e., green points, cluster at another phase. The points before and after 50 841 MJD cluster separately with a shift of 0.5 in phase (or about 13 d).

We have also folded the simulated data of Eq. (2) with Paverage. In this case the dependency on P1 and P2 is very simple, just two sine functions. Nevertheless, the same kind of double clustering shown by the GBI data occurs for the simulated data of Eq. (2) as one can see in the box of Fig. 4c. The fact that a simple sum of sine functions in P1 and P2 produces the same jump, as the GBI data, if folded with Paverage, implies either that this simple mathematical form is true in LS I +61°303, or, more likely, that the jump is a property of the beating process and its value only depends on the two periods P1 and P2.

In mathematical terms, as shown in the Appendix, the beat of the two sine functions f(P1) and f(P2) has a phase reset when the delay of f(P2) with respect to f(P1) becomes larger than P2/2. Until that point f(P1) preceeds f(P2) (and Paverage) giving rise to a positive timing residual. After that point f(P2) (and Paverage) preceeds f(P1). This produces the jump from a large timing residual to a nearly equally large, but now negative, timing residual observed in the sawtooth function (see Fig. 3c and d and the derived sawtooth function in Fig. 3b, in impressing agreement with the observed one of Fig. 3a).

In physical terms the jump is illustrated in Fig. 4d. The radio maximum of the long-term modulation results when the ejection, periodical at P = P1, occurs at the smallest angle with respect to the line of sight and the Doppler boosting is largest (Kaufman Bernadó et al. 2002). Because of the precession P2 the angle of the ejection changes, and the radio minimum corresponds to an ejection occurring at the largest angle with respect to the line of sight. At this point the ejection has travelled half a precession cone and turns onto the other half of the precession cone, i.e., from I to II in Fig. 4d. This causes the phase jump in the data folded with Paverage and the jump in the sawtooth function of Fig. 3a.

Finally, in Fig. 4d, we present the data folded with Paverage, where we fold the data before 50 841MJD (black points) using the usual t0 = 43   366.275MJD (Gregory 2002), whereas for the data after 50 841MJD (green points) we use t0 + 13.25 d, correcting for the jump at the minimum. The better folding of the large flares with Paverage with respect to those with the two real periodicities P1 and P2 is the observational evidence for the result of Eq. (6), i.e., the periodicity of the radio outburst is Paverage.

3. Periodicities in the equivalent width of the Hα emission line

Our results have important implications for the short- and long- Hα variations for the Be star of the LS I +61°303 system. Zamanov et al. (1999) analyzed Hα spectra of LS I +61°303 and determined the same 26.5 d radio period, and in addition determined that the Hα emission line equivalent width (EW) varies over the same time scale as the long-term radio modulation. Moreover, Zamanov et al. (1999) tried to find evidence for long-term periodicities in other line parameters like the important B/R ratio.

Most Be/X-ray binaries show asymmetric split Hα profiles, the “blue” (B) or “violet” (V) peak and the “red” (R) peak. B/R (or V/R) variability refers then to the variation of the relative strength of the blue to the red peak. B/R variability cycles are common in Be stars forming canonical Be/X-ray binaries containing accreting X-ray pulsars. As stated above, Zamanov and collaborators tried to find evidence in LS I +61°303 for long-term periodicities in the B/R ratio, but they conclude that “unfortunately, the results of this search turned out to be negative” (Zamanov et al. 1999).

V/R variability is explained in terms of a nonaxisymmetrical equatorial disk in which a one-armed perturbation (a zone in the disk with higher density) propagates (Reig 2011). As a possible explanation for the long-term modulation, Gregory & Neish (2002) indicated that it may stem from periodic ejections of a shell of gas in the equatorial disk of the Be star. Gregory & Neish (2002) commenting the fact that there are no periodic variations in the Hα V/R ratio for LS I +61°303, nevertheless tried to test the one-armed density wave model predictions. Their conclusion is that the radio behavior “is at odds with the predictions of a one-armed density wave model. Thus, the one-armed density wave model does not agree quantitatively with the measurements for LS I +61°303” (Gregory & Neish 2002).

In other words, the origin of the long-term variation of the Hα emission line from the disk of the Be star does not seem to be related to structural disk variations (no B/R variations), whereas it is clearly related to a periodical change in the number of emitted photons (the EW(Hα)). A likely candidate as agent of these variations is the relativistic precessing jet that could well be able to produce EW(Hα) variations of the same timescales as those we observe in its radio emission.

4. Conclusions and discussion

Our timing analysis and simulations give the following results:

  • 1.

    The timing analysis of 6.7 years of GBI data ofLS I +61°303 resultsin two frequenciesν1 = 0.03775 d-1 (P1 = 26.49 ± 0.07 d) and ν2 = 0.03715 d-1 (P2 = 26.92 ± 0.07 d). The aim of our research, to obtain a better determination of the precessional period, indicated by the astrometry to be of 27–28 d, has therefore been reached. The period exists and it is sligthly above the orbital period.

  • 2.

    An additional and totally unexpected result of our research is that these two periodicities give rise to a Paverage=2ν1+ν2=20.03775+0.03715=26.70±0.05d\hbox{$P_{\rm average} = \frac{2}{\nu_1 + \nu_2} = \frac{2} {0.03775 +0.03715}= 26.70 \pm 0.05~{\rm d}$}, modulated with Pbeat=1ν1−ν2=10.03775−0.03715=1667±393d\hbox{$P_{\rm beat}={1\over\nu_1-\nu_2}={1\over {0.03775 -0.03715}}= 1667 \pm 393~{\rm d}$}. In other words, the long-term periodicity is equal to the beat of P1 and P2, and Pbeat modulates a new period, Paverage.

  • 3.

    We have shown that the sawtooth function derived in the past (Gregory et al. 1999) by comparing observed and predicted (for P = P1) outburst time compares P1 to Paverage. Our result therefore confirms the controversial value of 26.69 ± 0.02 d obtained in the past by Ray et al. (1997) for the periodicity of the outburst.

  • 4.

    Paverage is only an apparent periodicity, a result of the beat of P1 and P2. When the GBI data are folded with Paverage the data, instead of clustering at one specific phase, cluster at two phases separated by about 0.5 (13.25 d, i.e., P1/2) depending on whether the data are before or after the minimum of the long-term modulation. We have shown that the beat of P1 and P2 reproduces the same double clustering as the GBI data. We find that the time in the minimum of the long-term modulation when the phase jump occurs is the time when the relative delay between the two functions f(P1) and f(P2) have reached the maximum delay of P2/2. In a physical scenario this corresponds, as discussed below, to the point where the ejection has travelled half a precession cone and turns on the other half of the precession cone.

  • 5.

    Paverage and Pbeat are related to each other. They are both produced by the beat between P1 and P2, the two real periodicities. Pbeat, the long-term periodicity, is like Paverage, just an apparent periodicity.

  • 6.

    P1 and the long-term modulation alone, cannot reproduce the observed periodogram (i.e., P2). It would be possible to reproduce the observed periodogram only by adding a sawtooth function. In this scenario, the sawtooth function must be produced by a physical process continuously changing P1 to Paverage until a shift in orbital phase of +6 d/26.49 d is reached. Then another process would suddenly shift the outburst from +6 d/26.49 d to –7 d/26.49 d in orbital phase, but without changing the amplitude of the outburst (remaining in the minimum). Instead, we have shown that using directly the two periods P1 and P2 that we found in the periodogram of the GBI data, one naturally explains the two apparent periods Pbeat and Paverage, and the phase jump in Paverage.

  • 7.

    Implications of our results for variations at other wavelengths indicate the precessing relativistic jet as the agent responsible for the observed Hα variations, equal to the variations that we observe in radio emission of the jet.

In conclusion, LS I +61°303 seems to be one more case in astronomy of a “beat”, i.e., a phenomenon occurring when two physical processes create stable variations of nearly equal frequencies. The very small difference in frequency creates a long-term variation of period 1/(ν1 − ν2). The first astronomical case was that of a class of Cepheids, afterwards called beat Cepheids (Oosterhoff 1957).

The two periods in LS I +61°303 are the precession P2 of an ejection having orbital occurrence P1. Following the results of the radio spectral index analysis by Massi & Kaufman Bernadó (2009), the steady ejection of relativistic electrons associated to the compact object in LS I +61°303 increases until it terminates in a transient jet twice along the orbit. As quoted in the introduction, this agrees well with the Bondi & Hoyle (1944) accretion in an eccentric orbit that predicts two events in the system LS I +61°303: one around periastron, and the second shifted towards apastron (Taylor et al. 1992; Marti & Paredes 1995; Bosch-Ramon et al. 2006; Romero et al. 2007). However, the ejected relativistic electrons around periastron suffer strong inverse Compton losses because they are exposed to stellar photons, and only at the second accretion/ejection peak, occuring rather displaced from periastron, the relativistic electrons survive inverse Compton losses and produce the large observed radio outburst with period P1 = 26.49 d (Bosch-Ramon et al. 2006). This ejection periodically changes direction (Massi et al. 2012) with a period established in the present paper of 26.92 d. The radio maximum of the long-term modulation results when the ejection occurs at the smallest angle with respect to the line of sight and the Doppler boosting is largest (Kaufman Bernadó et al. 2002). The radio minimum results when the ejection occurs at large angles with respect to the line of sight. The point where the ejection has travelled half a precession cone and turns on the other half of the precession cone gives rise to the jump in phase for the apparent Paverage.

This research has brought up many questions in need of further investigation. Future work needs to be done to establish the radio vs. Hα relationship, the possible physical processes under the observed precession, and the Doppler boosting effects. Future observations are needed to estimate the angle of the precession cone; an estimate of this angle could result from comparing VLBA observations at the same orbital phase but interlapsed (1667/2)d, that is by about 31 orbital cycles. From the different position angle of the radio structures one could infer the aperture of the cone. Of course the two images interlapsed by 31 cycles should be done with high precision at the same orbital phase. Compare in Fig. 1 in Massi et al. (2012) image A (or image B) with image I (or J) taken only one cycle later than A (B), but at a sligtly differente phase ΔΦ = 0.016. The small phase difference already causes differences in the images. In Albert et al. (2008) it is discussed that two images taken 10 cycles apart and at the same orbital phase Φ = 0.62 are higly similar. Of course 10 cycles are still far from the requested 31 cycles. However, the slightly different position angle between the two images in Fig. 1 of those authors indicates that this comparison could be a good investigative tool to estimate the precession angle in future observations.

2

(1P1−1P2)-1=1658±382\hbox{$\left({1\over P_1} - {1\over P_2}\right)^{-1} = 1658\pm 382$} when using ≥6 digits.

3

Fmod is a function implemented in the math library of C.

Acknowledgments

We would like to thank the anonymous referee for helpful comments. We thank Lisa Zimmermann and Jürgen Neidhöfer for reading the manuscript and for the several interesting discussions. The Green Bank Interferometer is a facility of the National Science Foundation operated by the NRAO in support of NASA High Energy Astrophysics programs.

References

  1. Albert, J., Aliu, E., Anderhub, H., et al. 2008, ApJ, 684, 1351 [NASA ADS] [CrossRef] [Google Scholar]
  2. Bondi, H., & Hoyle, F. 1944, MNRAS, 104, 273 [NASA ADS] [CrossRef] [Google Scholar]
  3. Bosch-Ramon, V., Paredes, J. M., Romero, G. E., & Ribó, M. 2006, A&A, 459, L25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Cornwell, T., & Fomalont, E. B. 1999, in Synthesis Imaging in Radio Astronomy II, eds. G. B. Taylor, C. L. Carilli, & R. A. Perley, ASP Conf. Ser., 180, 187 [Google Scholar]
  5. Dhawan, V., Mioduszewski, A., & Rupen, M. 2006, in Proc. VI Microquasar Workshop, 52.1 [Google Scholar]
  6. Fender, R. P., Bell Burnell, S. J., & Waltman, E. B. 1997, Vist. Astron., 41, 3 [NASA ADS] [CrossRef] [Google Scholar]
  7. Fender, R. P., Belloni, T. M., & Gallo, E. 2004, MNRAS, 355, 1105 [NASA ADS] [CrossRef] [Google Scholar]
  8. Fragile, P. C., Blaes, O. M., Anninos, P., & Salmonson, J. D. 2007, ApJ, 668, 417 [NASA ADS] [CrossRef] [Google Scholar]
  9. Gregory, P. C. 1999, ApJ, 520, 361 [NASA ADS] [CrossRef] [Google Scholar]
  10. Gregory, P. C. 2002, ApJ, 575, 427 [NASA ADS] [CrossRef] [Google Scholar]
  11. Gregory, P. C., & Neish, C. 2002, ApJ, 580, 1133 [NASA ADS] [CrossRef] [Google Scholar]
  12. Gregory, P. C., Peracaula, M., & Taylor, A. R. 1999, ApJ, 520, 376 [NASA ADS] [CrossRef] [Google Scholar]
  13. Grundstrom, E. D., Caballero-Nieves, S. M., Gies, D. R., et al. 2007, ApJ, 656, 437 [NASA ADS] [CrossRef] [Google Scholar]
  14. Han, X., & Hjellming, R. M. 1992, ApJ, 400, 304 [NASA ADS] [CrossRef] [Google Scholar]
  15. Hannikainen, D. C., Wu, K., Stevens, J. A., et al. 2006, Chin. J. Astron. Astrophys. Suppl., 6, 010000 [Google Scholar]
  16. Kaufman Bernadó, M. M., Romero, G. E., & Mirabel, I. F. 2002, A&A, 385, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. Kovalev, Y. Y., Lobanov, A. P., Pushkarev, A. B., & Zensus, J. A. 2008, A&A, 483, 759 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Linnell Nemec, A. F., & Nemec, J. M. 1985, AJ, 90, 2317 [NASA ADS] [CrossRef] [Google Scholar]
  19. Lomb, N. R. 1976, Ap&SS, 39, 447 [NASA ADS] [CrossRef] [Google Scholar]
  20. Marti, J., & Paredes, J. M. 1995, A&A, 298, 151 [NASA ADS] [Google Scholar]
  21. Martí-Vidal, I., & Marcaide, J. M. 2008, A&A, 480, 289 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Massi, M. 2005, Habilitation Thesis [arXiv:astro-ph/0506731] [Google Scholar]
  23. Massi, M., & Aaron, S. 1999, A&AS, 136, 211 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Massi, M., & Kaufman Bernadó, M. 2009, ApJ, 702, 1179 [NASA ADS] [CrossRef] [Google Scholar]
  25. Massi, M., & Zimmermann, L. 2010, A&A, 515, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Massi, M., Ribó, M., Paredes, J. M., et al. 2004, A&A, 414, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Massi, M., Ros, E., & Zimmermann, L. 2012, A&A, 540, A14 [Google Scholar]
  28. Mirabel, I. F. 2012, Science, 335, 175 [NASA ADS] [CrossRef] [Google Scholar]
  29. Moldón, J., Ribó, M., & Paredes, J. M. 2012, A&A, 548, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Oosterhoff, P. T. 1957, Bull. Astr. Inst. Netherlands, 13, 320 [Google Scholar]
  31. Peracaula, M., Martí, J., & Paredes, J. M. 1997, A&A, 328, 283 [NASA ADS] [Google Scholar]
  32. Pooley, G. G., & Fender, R. P. 1997, MNRAS, 292, 925 [NASA ADS] [CrossRef] [Google Scholar]
  33. Porter, J. M., & Rivinius, T. 2003, PASP, 115, 1153 [NASA ADS] [CrossRef] [Google Scholar]
  34. Ray, P. S., Foster, R. S., Waltman, E. B., et al. 1997, ApJ, 491, 381 [NASA ADS] [CrossRef] [Google Scholar]
  35. Reig, P. 2011, Ap&SS, 332, 1 [NASA ADS] [CrossRef] [Google Scholar]
  36. Rodriguez, L. F., & Mirabel, I. F. 1997, ApJ, 474, L123 [NASA ADS] [CrossRef] [Google Scholar]
  37. Romero, G. E., Okazaki, A. T., Orellana, M., & Owocki, S. P. 2007, A&A, 474, 15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Scargle, J. D. 1982, ApJ, 263, 835 [NASA ADS] [CrossRef] [Google Scholar]
  39. Stellingwerf, R. F. 1978, ApJ, 224, 953, 2317 [NASA ADS] [CrossRef] [Google Scholar]
  40. Stewart, I. M., Fenech, D. M., & Muxlow, T. W. B. 2011, A&A, 535, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Taylor, A. R., & Gregory, P. C. 1982, ApJ, 255, 210 [NASA ADS] [CrossRef] [Google Scholar]
  42. Taylor, A. R., & Gregory, P. C. 1984, ApJ, 283, 273 [NASA ADS] [CrossRef] [Google Scholar]
  43. Taylor, A. R., Kenny, H. T., Spencer, R. E., & Tzioumis, A. 1992, ApJ, 395, 268 [NASA ADS] [CrossRef] [Google Scholar]
  44. Urry, C. M., & Padovani, P. 1995, PASP, 107, 803 [NASA ADS] [CrossRef] [Google Scholar]
  45. Valtaoja, E., Terasranta, H., Urpo, S., et al. 1992, A&A, 254, 71 [NASA ADS] [Google Scholar]
  46. Zamanov, R. K., Martí, J., Paredes, J. M., et al. 1999, A&A, 351, 543 [NASA ADS] [Google Scholar]

Appendix A: Beat and sawtooth function

Let us assume two sine functions f1 and f2 with period P1 = 26.49 d and P2 = 26.92 d. Let us start during the radio maximum, i.e., with zero timing residual between the peaks of f1 and f2. In Fig. 4e a zero timing residual corresponds to an ejection with the smallest angle with respect to the line of sight, that is with the strongest Doppler boosting. Since P1 is shorter than P2, each cycle f(P2) has a delay of P2 − P1 = 0.43 d. The period formed by the beat Paverage has a timing residual τ between its peak and that of f(P1) half of that delay, i.e., 0.43/2 = 0.21 d. This agrees well with τ = 0.008 × P1 = 0.21 d of Eq. (4).

f(Paverage) is always between f1 and f2 and its peaks have a regular distance of Paverage = 26.70 d (Fig. 3d). However, at the minimum the accumulated delay of f(P2) with respect to f(P1) is almost P2/2. That is, as shown in Fig. 3c, the peak of f(P2) is nearly equidistant from the two peaks of f(P1), the peak about 13 d before and the peak about 13 d after it. The beat with the first peak gives rise to a peak of f(Paverage) at about 827.8 d, whereas the beat with the second peak gives rise to a peak of f(Paverage) at 841.05 d. The difference between these two consecutive peaks of f(Paverage) is 13.25 d. This corresponds to about 0.5 in phase, when the data are plotted in phase for Paverage (Fig. 4c), and also corresponds to the jump of about 13 d in the sawtooth function observed by Gregory et al. (1999). After that point, f(P2) is preceeding f(P1). In the physical scenario of Fig. 4e this jump, or reset of the phase of Paverage, corresponds to the point where the ejection has travelled half a precession cone and turns onto the other half of the precession cone. Figure 3b shows τ vs. time, resulting from this simple analysis based on the sine functions of Fig. 3; the resulting slope of the function is indeed 0.008, as in Fig. 3a.

All Figures

thumbnail Fig. 1

Periodograms of 8.3 GHz (filled squares) and 2.2 GHz (circles) data. a) Output of Lomb-Scargle method. The Lomb-Scargle analysis gives on the y-axis the significance of the frequency. Two frequencies are found at ν1 = 0.03775 d-1 (P1 = 26.49 d) and at ν2 = 0.03715 d-1 (P2 = 26.92 d). b) Output of Lomb-Scargle method for data ≥4σ. c) Output of PDM. The most likely period yields the minimum dispersion and appears as a minimum in the PDM curve, i.e., specular to the maximum in the Lomb-Scargle plot. d) Output of PDM for data ≥4σ.
In the text
thumbnail Fig. 2

Long-term modulation and period analysis. a) 8.3 GHz GBI radio data averaged over 3 d. b) Lomb-Scargle analysis results: two frequencies at ν1 = 0.03775 d-1 (P1 = 26.49 d) and ν2 = 0.03715 d-1 (P2 = 26.92 d). The small window shows the peak at 1/1667 d-1 present in all periodograms. c) Sum of two sinusoidal functions at 26.49 d and 26.92 d, with an amplitude ratio 1/0.7. d) Lomb-Scargle analysis results: two frequencies at 126.49d-1\hbox{$\frac{1}{26.49}~{\rm d^{-1}}$} and 126.92d-1\hbox{$\frac{1}{26.92}~{\rm d^{-1}}$}. The significance of the two frequencies becomes identical in the periodogram only for an amplitude ratio 1/1. e) Long-term modulation (1667 d) of a 26.49 d periodic outburst. f) Lomb-Scargle analysis results: one frequency at 126.49d-1\hbox{${1\over 26.49}~{\rm d}^{-1}$}. The small window to the left shows the peak at Paverage present in the periodogram of a simulation of long-term modulation (1667 d) of a 26.70 d periodic outburst. g) Sine wave of periodicity P1, modulated by a sine wave of periodicity 1667 d and corrected by a sawtooth function. h) Lomb-Scargle analysis results: two frequencies at 126.49d-1\hbox{$\frac{1}{26.49}~{\rm d^{-1}}$} and 126.92d-1\hbox{$\frac{1}{26.92}~{\rm d^{-1}}$} as in Figs. 2b and d.
In the text
thumbnail Fig. 3

a) Sawtooth function based on timing residuals observed by Gregory et al. (1999). b) Sawtooth function based on the loss of synchronization of the sine functions f(P1) and f(P2) of Fig. 3c. c) Sine waves with periods P1 (black line) and P2 (red line) starting with synchronized peaks at t = 0. d) Sum of the sine waves. The first five bars are at a regular interval equal to Paverage, then the next bar follows at a distance of only 13.25 d (P1/2), followed by four bars again at a regular interval equal to Paverage as before. The two central bars are symmetric with respect to a peak of f(P2) (Fig. 3c) that is nearly equidistant from the preceeding peak and the delayed peak of f(P1) (see Appendix).
In the text
thumbnail Fig. 4

a) Radio light curves (data averaged over 1 d) vs. Φorbit, with Φorbit related to (t−t0)P1\hbox{${(t-t_0)\over P_1}$}, with P1 = 26.49 d and t0 = JD2 443 366.775 (Gregory 2002), before (black) and after (green) (2 400 000.5 + 50 841)JD (see Appendix). b) Same as a) but for Φprecession, related to P2. c) Radio light curves vs. Φaverage for Paverage = 26.70 d. The small window shows the simulated data (Eq. (2)) identically folded. d) Radio light curve vs. Φaverage with t0 changed by Δt (Δt = 0 before (2 400 000.5 + 50 841)JD and Δt = 13.25 d after it, see Appendix). e) Sketch of the precessing jet in LS I +61°303 (out of scale).
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.