© ESO, 2013

1. Introduction

In 1995, Hjellming & Rupen (1995) discovered GRO J1655-40, the first case of a microquasar having a precessing radio jet with a precessional period P_prec = 3 ± 0.2d (Hjellming & Rupen 1995) rather close to the orbital period P_orbit = 2.601 ± 0.027d (Bailyn et al. 1995). The source LS I + 61°303 is a probable second case with a precessional period P₂ = 26.92 ± 0.07d of the radio jet (Massi & Jaron 2013) very close to the orbital period P₁ = 26.4960 ± 0.0028d (Gregory 2002; Grundstrom et al. 2007). Fast variations of the radio structure in LS I + 61°303 were the results from MERLIN observations (Massi et al. 2004) and have been confirmed by VLBA observations by Dhawan et al. (2006), there interpreted as due to a cometary tail of a pulsar. Radio spectral index analysis (Massi & Kaufman Bernadó 2009), 3D simulations (Romero et al. 2007), and VLBA images (Massi et al. 2012) indicated, however, a microquasar scenario and the radio astrometry of 27−28 d suggested a slightly different value of the precession period from the orbital one (Massi et al. 2012). Finally, timing analysis resulted in P₂ = 26.92 ± 0.07d (Massi & Jaron 2013). Lense-Thirring precession has been invoked as a possible mechanism in both objects, for GRO J1655-40 by Martin et al. (2008), and for LS I + 61°303 by Massi & Zimmermann (2010).

In addition to these two periods, i.e., the orbital one P₁ and the precessional one P₂, the radio properties of LS I + 61°303 feature two additional periods. Ray et al. (1997) determined a P = 26.69 ± 0.02d for the periodical large radio outburst. Gregory (2002) determined a P = 1667 ± 8d for a long-term periodicity modulating the radio outbursts. Massi & Jaron (2013) show that the intrinsic periodicities of the source are the orbital period P₁ = 1/ν₁ and the precessional period P₂ = 1/ν₂ giving rise to a beating which results in a P_average = 2/(ν₁ + ν₂) = P_outburst = 26.70 ± 0.05d, consistent with the period found by Ray et al. (1997), modulated by a P_long = 1/(ν₁ − ν₂) = 1667 ± 393d, consistent with the long-term periodicity determined by Gregory (2002).

The importance of LS I + 61°303 has grown after the detection of emission in the TeV band (Albert et al. 2006). Nowadays, the source is observed at all wavelengths of the electromagnetic spectrum. In addition LS I + 61°303 is now included in the search for neutrino sources. While the high energy emission could be due to Inverse Compton, i.e., of leptonic origin, there could also be a hadronic component. The source LS I + 61°303 and other gamma-ray binaries could therefore play a role as accelerators of Galactic cosmic rays (Romero et al. 2005; Abbasi et al. 2012). What sets LS I + 61°303 apart from other gamma-ray binaries is its periodicity; however, this has not yet been fully exploited because of the uncertainty in predictions of its outburst. In fact, P_average, that is the periodicity of the outbursts, cannot be used directly for predictions because of its periodical reset at the minimum of the long-term modulation inherent to the beating process (Massi & Jaron 2013). On the other hand, using P₁ to predict the occurrence of the radio outbursts implies timing residuals between predicted and observed outbursts up to ~7 d and is affected by a systematic trend with a sawtooth pattern (Gregory et al. 1999, see their Fig. 2). In this paper we aim to determine the best procedure for predicting the periodical radio outburst. We determine the reset affecting P_average and compare the procedures using P₁ and P_average.

2. Timing residuals

In this section we analyze the timing residuals defined as the difference between observed and predicted outburst times, for a given period. We use the Green Bank Interferometer (GBI) radio data at 8.3 GHz for this analysis. They cover the interval from 49379.975 MJD to 51823.441 MJD. During the minimum of the long-term modulation the outburst loses its characteristic shape and becomes a broad low curve without a definite peak (Paredes et al. 1990). We find this phase, defined by Paredes et al. (1990) as “quiet phase”, also in the GBI data, in the interval 50500 MJD to 51000 MJD, where no maxima occur and where consequently large timing residuals would affect our tests for P₁ and P_average. We therefore remove this interval that biases our analysis. This interval in terms of the long-term modulation phase Θ, $\mathrm{\Theta }\mathrm{=}\frac{\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}}{{\mathit{P}}_{\mathrm{long}}}\mathrm{-}\mathrm{int}\left(\frac{\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}}{{\mathit{P}}_{\mathrm{long}}}\right)\mathit{,}$(1)with t₀ = 43366.275MJD and P_long = 1667d, corresponds to Θ = 0.28−0.58, which is therefore the interval where no clear radio outbursts occur.

2.1. Prediction using the orbital period P₁

Gregory et al. (1999) showed (see their Figs. 2 and 8) that when radio outbursts are predicted to occur with P₁ then there are residuals of −5.5 to + 7 days, and that these residuals follow a sawtooth trend. In this section we analytically determine the sawtooth function and remove it from the data to finally obtain noise-limited residuals. These residuals are compared to residuals based on P_average in the next section.

We divided the GBI data into bins of size P = P₁. A direct search for the data point with the greatest flux within each bin gives a residual of $\mathrm{\Delta }{\mathit{t}}_{\mathit{i}}\mathrm{=}{\mathit{t}}_{\max \mathit{,}\mathit{i}}\mathrm{-}{\mathit{t}}_{\mathit{i}}\mathit{.}$(2)One sees that the resulting timing residuals in Fig. 1a are not dominated by noise; on the contrary, a systematic pattern is present. The shape, as pointed out by Gregory et al. (1999), is a sawtooth pattern with a gradual rise and a rapid fall, and periodicity P_long equal to that present in the peak flux. This long-term modulation was first calculated to be about 1584 d (Gregory et al. 1999) and later established to be 1667 ± 8 d (Gregory 2002). The analytical shape for the sawtooth function has been determined in Massi & Jaron (2013) to be¹$\mathit{\tau }\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathit{m}\mathrm{·}\mathrm{fmod}\hspace{0.17em}\mathrm{\left(}\mathit{t,}\mathrm{1667}\mathrm{\right)}\mathrm{+}\mathit{c}\mathit{.}$(3)We remove τ(t) from the timinig residuals shown in Fig. 1a and compute the standard deviation σ for the noise-limited residuals. We minimize σ by varying m. Figure 1b shows the difference between the data in Fig. 1a and the sawtooth function. Figure 1b, with σ = 1.87d, will therefore be our term of quantitative comparison when using P_average in the next section.

Fig. 1 a) Timing residuals for the radio outbursts expected to occur with orbital periodicity P1. b) Difference between the residuals from Fig. 1a and the sawtooth function discussed in Sect. 2.1.

Fig. 2 Top: standard deviation σ as a function of the periodicity Paverage of the radio outbursts. Bottom: timing residuals for the assumption that the outbursts occur with a periodicity of Paverage = 26.704d and that there is a phase jump during the long-term minimum.

2.2. Predictions with P_average

To predict the outburst occurrence in LS I + 61°303 one generally determines for a given epoch that time t for which Φ(P₁) = 0.6, where $\mathrm{\Phi }\mathrm{=}\frac{\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}}{{\mathit{P}}_{\mathrm{1}}}\mathrm{-}\mathrm{int}\left(\frac{\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}}{{\mathit{P}}_{\mathrm{1}}}\right)\mathit{,}$(4)with t₀ = 43366.275MJD. This method follows from the analysis of the orbital shift by Paredes et al. (1990). During the maximum of the long-term periodicity of 1667 d the outburst occurs at Φ(P₁) = 0.6. This, in terms of the timing residuals in Fig. 1a, corresponds to a delay of zero; afterwards the phase of the outbursts increases until a maximum of about Φ(P₁) = 0.9 (Paredes et al. 1990). In terms of the timing residuals in Fig. 1a, this corresponds to a timing residual of (0.9−0.6)P₁ ≃ 8d. The jump in the sawtooth function from maximum positive delay to maximum negative delay occurs when the peak of the outburst, after the nearly quiescent period at the minimum, appears again but at an earlier phase than Φ = 0.6, i.e., at phase Φ = 0.4 (Paredes et al. 1990). In terms of timing residuals, this corresponds to (0.4−0.6)P₁ = −5d.

Here we show that the procedure above described to predict the outburst occurrence should be replaced by the determination whether Φ(P_average) = 0.7 or Φ(P_average) = 0.2 (as in Eq. (4) with P_average instead of P₁), where the ambiguity between 0.7 or 0.2 is determined as shown in the next section. We demonstrate by the residuals how effective this procedure is.

We compute the timing residuals for P_average following the procedure of Sect. 2.1 for P = P_average; the residuals are shown in the lower panel of Fig. 2. The value in Massi & Jaron (2013) for P_average is 26.70 ± 0.05 d. Therefore, we compute the value of P_average minimizing the standard deviation σ of the residuals for a number of trial periodicities equally distributed in P = 26.65−26.75d, i.e., within the error bar given in Massi & Jaron (2013). In the plot, the minimum σ shows a plateau with center at 26.704d and a half width of 0.04 d, that define therefore our improved value for P_average and its error.

The residuals in Fig. 2 Bottom have σ = 1.60d, even slightly better than the “residuals of the residuals” with σ = 1.87d shown in Fig. 1b, obtained by the complex removal of the sawtooth function from the residuals in Fig. 1a.

2.3. Determination of the periodical reset

Figure 3 shows that the folding of the data with the orbital and precessional periods P₁ and P₂ results in one clustering in both cases, whereas the folding with P_average, which is simply their average and therefore one might expect an intermediate clustering, results in a double clustering. In particular, all data before the minimum of the long-term modulation P_long (black in Fig. 3) cluster at one phase, and all points after the minimum (green) cluster at the other phase. This is what a beat of P₁ and P₂ predicts: at the minimum of P_long the function in P₂ has arrived at a delay with respect to the function in P₁ equal to P₁/2, then one cycle later it preceeds the function in P₁ instead of lagging behind it as before the minimum. This gives rise to a jump in the function in P_average, being f(P_average) always between f(P₁) and f(P₂) (Massi & Jaron 2013).

From our tests we find that the phase Φ_outburst changes during every minimum of the long-term modulation such that ${\mathrm{\Phi }}_{\mathrm{outburst}}\mathrm{=}\{\begin{array}{c}\\ \mathrm{0.7}& \mathrm{for}\mathit{K}\mathrm{even}\mathit{,}\\ \mathrm{0.2}& \mathrm{for}\mathit{K}\mathrm{odd}\mathit{.}\end{array}$(5)Here, the integer number K counts the long-term modulation cycles for t ≥ 49174MJD $\mathit{K}\mathrm{=}\mathrm{int}\left(\frac{\mathit{t}\mathrm{-}\mathrm{49174}\mathrm{MJD}}{{\mathit{P}}_{\mathrm{long}}}\right)\mathit{,}$(6)where 49174MJD = (50841−1667)MJD is coincident with the time of the phase jump related to the GBI data (Massi & Jaron 2013).

For the purpose of illustrating the predicted times of outburst we use the function $\mathit{f}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}\mathit{A}{\left[\frac{\mathrm{1}}{\mathrm{2}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\cos \left(\frac{\mathrm{2}\mathit{\pi }}{\mathit{P}}\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}\mathrm{)}\mathrm{-}\mathrm{2}\mathit{\pi }{\mathrm{\Phi }}_{\mathrm{outburst}}\right)\right]}^{\mathit{n}}\mathit{,}$(7)where A is a scaling factor which we set to 150 mJy to roughly match the average amplitude of the outbursts, P is the periodicity of the outbursts, t₀ = 43366.275MJD as usual, Φ_outburst is the phase of the outbursts, and n is an exponent to sharpen the peaks; we chose n = 8. In Fig. 4 two intervals of the GBI 8.3 GHz data are plotted together with the predicted outbursts. The continuous curve is the prediction made by Eq. (7) with P = P_average = 26.704d and Φ_outburst computed according to Eq. (5). The dashed curve is the result inserting P₁ and Φ_outburst = 0.6. There are times when the prediction with the orbital period P₁ is quite accurate, as shown in Fig. 4b, but in Fig. 4a, e.g., the deviations from the predictions are as large as several days. The prediction with P = 26.704d and setting Φ_outburst according to the number K is in good agreement with the observations in both cases.

Fig. 3 Left: 8.3 GHz GBI radio data vs. Φorbit, where Φorbit(= P1) = 26.49d. Middle: 8.3 GHz GBI radio data vs Φprecession, where Φprecession(= P2) = 26.92d. Right: 8.3 GHz GBI radio data vs. Φaverage for ${\mathit{P}}_{\mathrm{average}}\mathrm{=}\frac{\mathrm{2}}{{\mathit{\nu }}_{\mathrm{1}}\mathrm{+}{\mathit{\nu }}_{\mathrm{2}}}\mathrm{=}\mathrm{26.704}\mathrm{d}$. Data before the minimum at 50841 MJD are black circles, data after 50841 MJD are green squares.

Fig. 4 Illustration of the predicted outbursts, GBI 8.3 GHz data in red, prediction with P = P1 as a dashed line, and prediction with P = Paverage with phase jump as a continuous line, left: time interval with a clear difference between the two predictions, right: time interval where the prediction with P1 is also in good agreement with the observation.

3. Conclusions

The usual procedure used until now to predict the outburst occurrence is based on the orbital period P₁ = 26.4960 ± 0.0028d. Recently, Massi & Jaron (2013) determined that the periodicity of the radio outburst is P_average = 26.70 ± 0.05d. However, the direct use of P_average is not possible because of the phase jump during the minimum of the long-term modulation. In this paper we determined a straightforward way of calculating the phase term and therefore of predicting the radio outbursts using P_average. We obtained the following results:

• 1.

  We showed that the standard deviation of the timing residualsbecomes minimal and equal to σ = 1.6d for a periodicity of P_average = 26.704 ± 0.004d when takinginto account the phase jump of 0.5 during the minimum of thelong-term modulation. This value of P_average is our new improvedestimation of the periodicity of the observed radio outbursts.

• 2.

  The quantity which resets the phase Φ_outburst is the integer number of long-term cycles K as introduced in Eq. (6). The phase Φ_outburst is either 0.2 or 0.7 depending on whether K is odd or even, respectively, as described in Eq. (5).

• 3.

  We confirm the results of Paredes et al. (1990) that during the minimum of the long-term modulation not only is the outburst low, but it even becomes a broad curve without a definite peak. We established that this occurs in the interval Θ = 0.28−0.58.

Acknowledgments

We thank James Urquhart, Lisa Zimmermann, and Jürgen Neidhöfer for reading the manuscript and for useful discussions and the anonymous referee for constructive comments which improved the paper. 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

All Figures

	Fig. 1 a) Timing residuals for the radio outbursts expected to occur with orbital periodicity P1. b) Difference between the residuals from Fig. 1a and the sawtooth function discussed in Sect. 2.1.
In the text

	Fig. 2 Top: standard deviation σ as a function of the periodicity Paverage of the radio outbursts. Bottom: timing residuals for the assumption that the outbursts occur with a periodicity of Paverage = 26.704d and that there is a phase jump during the long-term minimum.
In the text

Fig. 3 Left: 8.3 GHz GBI radio data vs. Φorbit, where Φorbit(= P1) = 26.49d. Middle: 8.3 GHz GBI radio data vs Φprecession, where Φprecession(= P2) = 26.92d. Right: 8.3 GHz GBI radio data vs. Φaverage for ${\mathit{P}}_{\mathrm{average}}\mathrm{=}\frac{\mathrm{2}}{{\mathit{\nu }}_{\mathrm{1}}\mathrm{+}{\mathit{\nu }}_{\mathrm{2}}}\mathrm{=}\mathrm{26.704}\mathrm{d}$. Data before the minimum at 50841 MJD are black circles, data after 50841 MJD are green squares.

In the text

	Fig. 4 Illustration of the predicted outbursts, GBI 8.3 GHz data in red, prediction with P = P1 as a dashed line, and prediction with P = Paverage with phase jump as a continuous line, left: time interval with a clear difference between the two predictions, right: time interval where the prediction with P1 is also in good agreement with the observation.
In the text