© ESO, 2014

1. Introduction

A recent timing analysis of the GBI radio data of LS I +61°303 has revealed two frequencies: and (Massi & Jaron 2013). The period P₁ agrees with the value of (Gregory 2002) associated to the orbital period of the binary system formed by a compact object and a Be star (Grundstrom et al. 2007). Period P₂ agrees with the previous estimate of 27−28 days by radio astrometry of the precessional period of the radio structure (Massi et al. 2012).

The timing analysis results seem to confirm LS I +61°303 as the second case of a radio-emitting X-ray binary with precessional period (P₂) of its radio jet very close to the orbital (P₁) period of the binary system (Massi et al. 2012; Massi & Jaron 2013). Whereas the several known precessing X-ray binaries (Larwood 1998) have a precession period an order of magnitude longer than the orbital period (as predicted for a tidally forced precession by the companion star), the different case of GRO J1655-40 was discovered in 1995 (Hjellming & Rupen 1995). This black hole candidate, which has an orbital period of 2.601    ±    0.027 days (Bailyn et al. 1995), revealed a radio jet with a precessional period of 3.0   ±   0.2 days (Hjellming & Rupen 1995). LS I +61°303 is a high-mass X-ray binary, and GRO J1655-40 is a low-mass X-ray binary. The mechanism proposed in both cases to explain the precession is based on general relativity effects around the compact object. Lense-Thirring precession has been analysed for GRO J1655-40 by Martin et al. (2008) and for LS I +61°303 by Massi & Zimmermann (2010).

When two frequencies are only slightly different, a beating results, i.e. an interference producing a new frequency: their average, ${\mathit{\nu }}_{\mathrm{average}}\mathrm{=}\frac{{\mathit{\nu }}_{\mathrm{1}}\mathrm{+}{\mathit{\nu }}_{\mathrm{2}}}{\mathrm{2}}$, which gets modulated with the beat frequency, ν_beat = ν₁ − ν₂. As shown in Massi & Jaron (2013) for LS I +61°303, the term 1/ν_beat results and is compatible with the observed long-term modulation of the radio flux density of 1667    ±    8 days, attributed in the past to variations in the wind of the Be star (Gregory & Neish 2002). Evidence of beating in LS I +61°303 comes from the fact that the periodicity of the radio outbursts is indeed the one predicted by the beat theory, which is ${\mathit{P}}_{\mathrm{outburst}}\mathrm{=}{\mathit{P}}_{\mathrm{average}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{\nu }}_{\mathrm{average}}}\mathrm{=}\mathrm{26.70}\hspace{0.17em}\mathrm{\pm }\hspace{0.17em}\mathrm{0.05}$ days (Massi & Jaron 2013; Jaron & Massi 2014). Since the two frequencies ν_average and ν_beat are related, the long-term modulation should also be a result of the beating process between P₁ an P₂.

The aim of this paper is to prove, in the context of a physical scenario, how the mutual relationship between P₁ and P₂ can generate the long-term periodicity. We linked the two periodicities, P₁ and P₂, to two physical processes: the periodical (P₁) increase of relativistic electrons in a conical jet and the periodical (P₂) Doppler boosting of the emitted radiation by relativistic electrons because of jet precession. The hypothesis that the compact object in LS I +61°303, accreting material from the Be wind, undergoes a periodical (P₁) increase in the accretion at a particular orbital phase along an eccentric orbit was suggested and developed by several authors (Taylor et al. 1992; Marti & Paredes 1995; Bosch-Ramon et al. 2006; Romero et al. 2007). The hypothesis that a precessing jet, with an approaching jet having large excursions in its position angle, should give rise to appreciable Doppler boosting effects is supported by the morphology of Massi et al. (2004), Dhawan et al. (2006), and Massi et al. (2012) images showing extended radio structures changing from two-sided to one-sided morphologies at different position angles.

Our model describes a precessing (P₂) jet with a periodically (P₁) varying number of relativistic electrons. Section 2 describes our analytical model and how to calculate its flux density S_model. In Sect. 3, all parameters for calculating S_model are derived by fitting the model to 6.5 years of GBI radio flux density observations. In Sect. 4 we present our results and in Sect. 5 our conclusions.

2. Model definition

In this section we want to analytically determine the flux density of a precessing radio jet having a periodical variation in the number of its emitting particles. Section 2.1 describes the type of emitting source we are assuming for our model, which is the conical jet used for microquasars and AGNs after the seminal paper of Blandford & Königl (1979). In particular, we adopt Kaiser’s analytical treatment (Kaiser 2006) to express the changes of magnetic field and plasma properties along the jet axis. The flux density of our jet, which is S_model in Eq. (1), contains two contributions: the Doppler factors and the intrinsic synchrotron emission of the conical jet. The Doppler factors are determined in Sect. 2.2.

The most important parameter introduced to define the Doppler boosting term is the angle η between jet axis and line of sight. This is the angle that continuously varies because of precession. Section 2.3 deals with the second contribution of Eq. (1), i.e., the intrinsic synchrotron emission of the conical jet. The difference with respect to Kaiser’s model, where the jet is seen at η ~ 90°, is that here we are dealing with a precessing jet, i.e., with a continuously varying angle η. This implies the calculation of the optical depth through a stratified jet (Sect. 2.3 and Appendix). Also, in our case, the number of relativistic particles, N_rel, changes along the orbit of LS I +61°303, as described in Sect. 2.4.

2.1. The model

The general framework for LS I +61°303 emission in the radio energy range is that of synchrotron radiation emitted from a relativistic electron population of density N_rel moving in a magnetic field B. However, this simple scenario does not reproduce the radio spectra of LS I +61°303. A simple magnetized plasma containing particles with a power-law energy distribution produces a power-law spectrum (S ∝ ν^α) with spectral slope α < 0 above the critical frequency, while below this critical frequency self-absorption effects become important and the emission becomes optically thick with a spectral slope α = 2.5 (Longair 1994; Kaiser 2006). This spectrum of a uniform, self-absorbed synchrotron source is not the kind of radio spectrum observed in LS I +61°303. The optically thick emission observed in LS I +61°303 is much flatter with α in the range of 0.0 − 0.5 (Massi & Kaufman Bernadó 2009; Zimmermann 2013), similar to the values α ~ 0.0 − 0.6 observed in microquasars (Fender 2001). Optically thick synchrotron emission in microquasars and in radio loud AGN is modelled with conical jets where changes in magnetic field and plasma density along the flow give rise to the observed flat spectra (Blandford & Königl 1979; Kaiser 2006).

In the framework of this kind of jet model, described by unconstant physical quantities, we follow Kaiser’s model for adiabatic jets (Kaiser 2006), to which one can refer for more details. We use his notation when not otherwise specified. We therefore assume a conical jet, where all quantities are radially constant and change only along the jet axis x, with x expressed as x = x₀l, where x₀ is an arbitrary position along the x-axis defining the dimensionless coordinate l (see Fig. A.1 in Appendix). The cone radius is r(l) = r₀l^a₁, and we choose a₁ = 1 to have a cone with a constant opening angle ξ (in Kaiser’s notation tanξ = r₀/x₀).

The evolution of the strength of the magnetic field along the jet is parametrized as in Kaiser (2006) as B = B₀l^− a₂. Depending on the magnetic field geometry, a purely parallel magnetic field B ≡ B_∥ to the jet axis or a purely perpendicular magnetic field B ≡ B_⊥ to the jet axis imply different values for a₂. From Kaiser’s Table 1 we have for B ~ B_∥, a₂ = 2, while for B ~ B_⊥, a₂ = 1. As a matter of fact, the most likely configuration is a helical magnetic field configuration with both the parallel and the perpendicular components. Therefore, in the procedure described in Sect. 3 we have calculated the jet emission in both of these limiting cases.

The number density of the relativistic electrons has the usual power-law dependency on the electron energy: N_rel = κE^− p. As in Kaiser (2006), we represent the evolution of the electron distribution density along the jet by setting κ = κ₀l^− a₃, with a₃ = 2(2 + p)/3 and p = 1.8 (see Sect. 4.1). But in our case we also introduce a temporal dependence for electron density, setting κ₀ as function of time with periodicity P₁ (see Sect. 2.4).

In this framework we express the total flux density of the two jets, the approaching, and the receding jet, as ${\mathit{S}}_{\mathrm{model}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{=}{\mathit{S}}_{\mathrm{a}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathrm{\left(}{\mathit{\delta }}_{\mathrm{a}}\mathrm{\right(}\mathit{t}\mathrm{)}{\mathrm{)}}^{\mathit{k}\mathrm{-}\mathit{\alpha }}\mathrm{+}{\mathit{S}}_{\mathrm{r}}\mathrm{(}\mathit{t}\mathrm{\left)}\mathrm{\right(}{\mathit{\delta }}_{\mathrm{r}}\mathrm{\left(}\mathit{t}\mathrm{\right)}{\mathrm{)}}^{\mathit{k}\mathrm{-}\mathit{\alpha }}$(1)where δ_i is the Doppler boosting term, with i = a for the approaching jet and i = r, for the receding jet. The emission term, S_i(t), similarly represents the emission of each jet. Parameter k accounts for the ejecta properties, with k = 2 for a continuous jet (used here) and k = 3 for discrete condensations; α is the spectral index.

The quantities S_i and δ_i depend on the two physical periodicities P₁, the orbital period, and P₂, the precessional period, as described in the following.

2.2. The Doppler boosting effect and its dependency on P₂

Synchrotron radiation from relativistic electrons spiralling in the magnetic field of the approaching jet becomes enhanced by the Doppler effect (i.e., “boosted”), while the radiation from electrons of the receding jet becomes attenuated (i.e., “de-boosted”). Doppler boosting and de-boosting strongly depend on the angle between the jet and the observer’s line of sight (Urry & Padovani 1995; Mirabel & Rodríguez 1999). A precession of the jet implies a variation in this angle and therefore a variable Doppler boosting.

The Doppler boosting and de-boosting terms δ_i(t) are defined as $\begin{array}{ccc}{\mathit{\delta }}_{\mathrm{a}}& \mathrm{=}& \frac{\sqrt{\mathrm{1}\mathrm{-}{\mathit{\beta }}^{\mathrm{2}}}}{\mathrm{1}\mathrm{-}\mathit{\beta }\cos \mathit{\eta }}\\ {\mathit{\delta }}_{\mathrm{r}}& \mathrm{=}& \frac{\sqrt{\mathrm{1}\mathrm{-}{\mathit{\beta }}^{\mathrm{2}}}}{\mathrm{1}\mathrm{+}\mathit{\beta }\cos \mathit{\eta }}\end{array}$where β = v/c (with v the electron velocity component along the jet axis).

Angle η is the angle between observer’s line of sight and the jet axis. To follow its changes with time due to the jet precession we express η in terms of other quantities, throughout the spherical triangle of sides ζ,ψ, and η shown in Fig. 1. Angle ζ is the angle the jet precession axis makes with the direction of the line of sight. If the jet precession axis is tilted by an angle f with respect to the orbital axis of the binary system, then ζ = i + f, where i is the system inclination angle; ψ is the opening angle of the cone described by the jet during its precession period, P₂. In the spherical triangle of Fig. 1, the angle Ω, between sides ψ and ζ, changes with time because of precession and is defined as $\mathrm{\Omega }\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}\mathrm{+}\mathrm{\Delta }\mathrm{)}\mathit{/}{\mathit{P}}_{\mathrm{2}}\mathit{,}$(4)where t₀ = JD 2 443 366.775 (Gregory 2002, and references therein), and Δ is the offset to be determined in our model to synchronize P₂ to P₁.

From spherical geometry (Smart & Green 1977, Eq. (8)) it follows that the time (i.e. Ω) dependency of the angle η is $\cos \mathit{\eta }\mathrm{=}\cos \mathit{\zeta }\cos \mathit{\psi }\mathrm{+}\sin \mathit{\zeta }\sin \mathit{\psi }\cos \mathrm{\Omega }\mathit{.}$(5)Inserting the η expression into Eqs. (2) and (3), we get a temporal modulation of the Doppler boosting factor and we introduce in the observed flux density the dependency on the precession period P₂.

2.3. Intrinsic jet emission

In our precessing conical jet model, angle η, between the jet axis and the line of sight, must be clearly variable and can assume all values in the range 90° > η > ξ, with ξ the jet opening angle. This is an important difference with respect to Kaiser’s model, where the jet is seen almost perpendicular to its axis (i.e. η ~ 90°), or to Potter & Cotter (2012)’s model where the jet is observed at an angle smaller than the opening angle of the jet (i.e., η ≤ ξ).

Since in our model we allow time-dependent values for angle η, the plasma optical depth also will depend on time. In fact, the optical depth must be computed across the jet along a path that makes an angle η with the jet axis and intersects regions characterized by different values of physical quantities (since in our model all quantities vary along the jet axis, i.e., along x = x₀l) (Fig. A.1-bottom). To make the computation easier, we assume a pyramidal shape for the jet with a square basis of side 2r₀l and with a lateral surface facing the line of sight. In that way the integration path inside the jet does not depend on other parameters related to the curved jet surface. The perpendicular direction to this surface makes an angle [90° − (η − ξ)] with the line of sight for the case of the approaching jet and [90° − (η + ξ)] for the receding jet (see Fig. A.1-top in the Appendix).

With these assumptions the approaching jet emission can be written as ${\mathit{S}}_{\mathrm{a}}\mathrm{=}{\mathrm{\int }}_{\mathrm{1}}^{\mathit{L}}{\mathit{r}}_{\mathrm{0}}{\mathit{x}}_{\mathrm{0}}{\mathit{I}}_{{\mathit{\nu }}_{\mathrm{a}}}\mathrm{\left(}\mathit{\eta ,l}\mathrm{\right)}\frac{\sin \mathrm{(}\mathit{\eta }\mathrm{-}\mathit{\xi }\mathrm{)}}{{\mathit{D}}^{\mathrm{2}}}\mathit{l}\mathrm{d}\mathit{l}$(6)where D is the distance of LS I +61°303, and the integral is from l = 1 to l = L, with x₀L being the jet length and L operationally determined as the limiting value above which there are no more contributions to the jet flux. For the receding jet, we have the same expression with a change of sign in the angular factor that takes the projection of the emitting surface in the direction perpendicular to the line of sight into account. Clearly the difference between the emission of the approaching and of the receding jet becomes negligible for cases where ξ ≪ η.

The intensity I_ν, emerging from the jet surface, i.e., at τ = 0, in the direction η, is deduced from the radiative transport equation when considering that the plasma inside the jet is stratified along the jet axis direction, i.e., along a direction that makes an angle η with respect to the line of sight ${\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{\eta ,l}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{\tau }}_{\mathrm{end}}\mathrm{\left(}\mathit{l}\mathrm{\right)}}\frac{{\mathit{J}}_{\mathit{\nu }}}{{\mathit{\chi }}_{\mathit{\nu }}}{\mathrm{e}}^{\mathrm{-}{\mathit{\tau }}^{\mathrm{\prime }}\mathit{/}\cos \mathit{\eta }}\mathrm{d}\left[\frac{{\mathit{\tau }}^{\mathrm{\prime }}}{\cos \mathit{\eta }}\right]\mathit{,}$(7)where τ is the optical depth defined as $\mathit{\tau }\mathrm{\left(}\mathit{l}\mathrm{\right)}\mathrm{=}\mathrm{-}{\mathrm{\int }}_{\mathit{l}}^{\mathrm{\infty }}{\mathit{\chi }}_{\mathit{\nu }}\mathrm{d}\mathit{x}\mathrm{=}{\mathit{\tau }}_{\mathrm{0}}{\mathit{l}}^{\mathrm{1}\mathrm{-}{\mathit{a}}_{\mathrm{3}}\mathrm{-}\mathrm{(}\mathit{p}\mathrm{+}\mathrm{2}\mathrm{)}{\mathit{a}}_{\mathrm{2}}\mathit{/}\mathrm{2}}\mathit{.}$(8)The upper integration limit τ_end(l) represents the jet optical depth at the specific value of l and takes into account that, at each distance, x₀l, from the beginning of the jet, the integration path inside the jet involves different regions of the jet itself (see Appendix). The optical depth inside the jet is also a function of the angle η (see expressions (A.5) and (A.6)) inducing higher values of the optical depth when the jet is pointing toward the observer, i.e. for low η values (but in the limit η ≥ ξ). The consequences of this dependence will be analysed further in a subsequent paper.

The exponent for l in expression (8), when taking Kaiser’s Table 1 into account for all possible a_i values, is always negative, so the optical depth is at its maximum at the beginning of the jet: ${\mathit{\tau }}_{\mathrm{0}}\mathrm{=}\mathit{\tau }\mathrm{(}\mathit{l}\mathrm{=}\mathrm{1}\mathrm{)}\mathrm{=}\frac{{\mathit{\chi }}_{\mathrm{0}}{\mathit{x}}_{\mathrm{0}}}{\mathrm{-}\mathrm{[}\mathrm{1}\mathrm{-}{\mathit{a}}_{\mathrm{3}}\mathrm{-}\mathrm{(}\mathit{p}\mathrm{+}\mathrm{2}\mathrm{)}{\mathit{a}}_{\mathrm{2}}\mathit{/}\mathrm{2}\mathrm{]}}{\mathit{\nu }}^{\mathrm{-}\mathrm{(}\mathit{p}\mathrm{+}\mathrm{4}\mathrm{)}\mathit{/}\mathrm{2}}\mathrm{·}$(9)Emission and absorption quantities in the above expressions are given, as in Longair (1994) and following the Kaiser (2006) parametrization, by $\begin{array}{ccc}{\mathit{J}}_{\mathit{\nu }}& \mathrm{=}& {\mathit{J}}_{\mathrm{0}}{\mathit{\nu }}^{\mathrm{-}\mathrm{(}\mathit{p}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{/}\mathrm{2}}{\mathit{l}}^{\mathrm{-}{\mathit{a}}_{\mathrm{3}}\mathrm{-}\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{a}}_{\mathrm{2}}\mathit{/}\mathrm{2}}\frac{\mathrm{W}}{{\mathrm{m}}^{\mathrm{3}}\mathrm{Hz}}\\ {\mathit{\chi }}_{\mathit{\nu }}& \mathrm{=}& {\mathit{\chi }}_{\mathrm{0}}{\mathit{\nu }}^{\mathrm{-}\mathrm{(}\mathit{p}\mathrm{+}\mathrm{4}\mathrm{)}\mathit{/}\mathrm{2}}{\mathit{l}}^{\mathrm{-}{\mathit{a}}_{\mathrm{3}}\mathrm{-}\mathrm{(}\mathit{p}\mathrm{+}\mathrm{2}\mathrm{)}{\mathit{a}}_{\mathrm{2}}\mathit{/}\mathrm{2}}{\mathrm{m}}^{-1}\mathit{,}\end{array}$where ${\mathit{J}}_{\mathrm{0}}\mathrm{=}\mathrm{2.3}\mathrm{\times }{\mathrm{10}}^{-25}\mathrm{(}\mathrm{1.3}\mathrm{\times }{\mathrm{10}}^{\mathrm{37}}{\mathrm{)}}^{\mathrm{(}\mathit{p}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{/}\mathrm{2}}\mathit{a}\mathrm{(}\mathit{p}\mathrm{)}{\mathit{B}}_{\mathrm{0}}^{\mathrm{(}\mathit{p}\mathrm{+}\mathrm{1}\mathrm{)}\mathit{/}\mathrm{2}}\mathit{\kappa }0$(12)and ${\mathit{\chi }}_{\mathrm{0}}\mathrm{=}\mathrm{3.4}\mathrm{\times }{\mathrm{10}}^{-9}\mathrm{(}\mathrm{3.5}\mathrm{\times }{\mathrm{10}}^{\mathrm{18}}{\mathrm{)}}^{\mathit{p}}\mathit{b}\mathrm{(}\mathit{p}\mathrm{)}{\mathit{B}}_{\mathrm{0}}^{\mathrm{(}\mathit{p}\mathrm{+}\mathrm{2}\mathrm{)}\mathit{/}\mathrm{2}}{\mathit{\kappa }}_{\mathrm{0}}\mathit{,}$(13)with the constants a(p) and b(p) as in Longair (1994) and κ₀ as determined in the next section.

2.4. Periodical increase in relativistic particles along the orbit

To take into account the observed orbital modulation in the radio flux we make the working hypothesis that the relativistic electron input in the jet changes with time. This hypothesis is plausible since as discussed below, the accretion is variable along an eccentric orbit. Our type of representation therefore includes the classical hypothesis that the presence of relativistic electrons in the jet is related to matter accretion on the compact object but makes no hypothesis for how they have been accelerated.

For LS I +61°303, with e = 0.72 ± 0.15 (Casares et al. 2005), the Bondi & Hoyle (1944) accretion theory in an eccentric orbit predicts two maxima (Taylor et al. 1992; Marti & Paredes 1995; Bosch-Ramon et al. 2006; Romero et al. 2007). One maximum is around periastron at Φ_periastron = 0.23 ± 0.02 (Casares et al. 2005), with the orbital phase Φ defined as the fractional part of $\frac{\mathrm{(}\mathit{t}\mathrm{-}{\mathit{t}}_{\mathrm{0}}\mathrm{)}}{{\mathit{P}}_{\mathrm{1}}}$. However, near periastron the ejected relativistic electrons suffer severe inverse Compton losses because of the proximity to the Be star; a high energy outburst is expected but no or negligible radio emission (Bosch-Ramon et al. 2006). For the second accretion peak, the compact object star is much farther from the Be star, and inverse Compton losses are lower. The relativistic electrons can propagate out of the orbital plane, and we observe a radio outburst (Bosch-Ramon et al. 2006).

In modelling the radio emission we assume a periodical function for relativistic electron density: ${\mathit{\kappa }}_{\mathrm{0}}\mathrm{=}\mathit{A}\mathit{f}\mathrm{(}\mathrm{\Phi }\mathrm{-}{\mathrm{\Phi }}_{\mathrm{r}}\mathrm{)}$(14)where f(Φ − Φ_r) is equal to 1 for Φ = Φ_r and expresses (times A) the periodical (P₁) increase in relativistic electrons due to Bondi’s accretion in an eccentric orbit. In this first stage (focussed on Doppler boosting effects and on data periodicity reproduction), we do not include relativistic electron energetic losses but instead allow a different slope for the onset and for the decay so as to take into account that electrons take some time to lose their energy. The top of Fig. 4 shows the resulting κ₀ function for the fit solution illustrated in Figs. 2−7.

3. Model fitting

Table 1 reports the parameters, described in the previous section, needed to calculate, by Eq. (1), the flux density of our model, S_model. The values for the parameters were found by minimizing the χ² between S_model, computed at a given t_i, and the simultaneous observed flux density with the Green Bank interferometer at 8 GHz. That is, we found those values of the parameters that give the lowest value of χ²: ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{i}}\frac{\mathrm{\left(}{\mathit{S}}_{\mathrm{GBI}}\mathrm{\right(}{\mathit{t}}_{\mathit{i}}\mathrm{)}\mathrm{-}{\mathit{S}}_{\mathrm{model}}\mathrm{(}{\mathit{t}}_{\mathit{i}}\mathit{,}\mathrm{parameters}\mathrm{)}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}}$(15)where σ is the error of GBI data at 8 GHz.

We used the function minimization software MINUIT of the CERN Program Library. In particular, we first set an initial model using the program SCAN, then we optimize the model with the program SIMPLEX to perform the minimization using the simplex method of Nelder & Mead (1965). The best solution is for B ≡ B_∥, ζ = 46, ψ = 40, Δ =  −922.5, β = 0.28, A = 3.84, Φ_r = 0.6, L = 14, and ξ = 1.5. All figures (Figs. 2–7) have then been made with this model. However, the non-linearity of the problem makes the solution depending on the initial model, and we also found other solutions with only slight worse fits. This uncertainty is considered in the parameter errors (shown in Table 1) that give the dispersion with respect to the mean value found for each parameter in the different solutions. The solutions for the different magnetic configurations, i.e. cases with B ≡ B_∥ or B ≡ B_⊥ are within the given errors. The solutions in Table 1 therefore apply to both cases of B.

4. Results and discussion

We compare here model data and observations, both given as a function of time and of the phases, as the orbital phase Φ and the long-term modulation phase Θ. Model data have been generated with the parameter values as specified in Fig. 2 caption.

4.1. Long-term modulation and Doppler boosting effects

Figure 2a shows GBI data (averaged over 2 days) and model data, i.e. S_model, vs. time. The model reproduces the data at this “large scale”. Now we analyse the plots in detail. In Fig. 2a, two bars are given around the minimum of the long-term modulation around 50 900 MJD and two other bars at a higher flux density phase of LS I +61°303 around 51 375 MJD. In terms of the long-term modulation phase Θ, the two time intervals correspond to Θ = 0.52 and Θ = 0.83, respectively.

One can see the zoom of the outbursts at Θ = 0.83 in Fig. 2b and how the model also reproduces on a “fine” scale the periodical large outburst occurring at that Θ phase. The model data are folded as a function of the orbital phase Φ in Fig. 3a. Along with S_model of Eq. (1), we also plot the intrinsic emission of the approaching jet S_a (blue points) and that of the receding jet S_r (red points). We note that the intrinsic emission of the two jets is different. This difference between the approaching and the receding jet, which is due to optical depth effects and is important for low values of η, will be analysed in detail in a subsequent paper. We now discuss the effects of Doppler boosting. Figures 3e and f show (black points) the Doppler factors δ^k − α for k = 2 and α(t) calculated from GBI observations at 2.2 GHz and 8.3 GHz. As a comparison, the violet curve in the figures corresponds to Doppler factors with α = (1 − p)/2 and p equal to the value p = 1.8 used in our model for the electron energy distribution (see Sect. 2.1). As can be seen in Fig. 3c, the minimum η value, which determines the maximum Doppler boosting, occurs at Φ = 0.56, which is rather close to the orbital phase Φ = 0.6 of the peak of the intrinsic emission S_a (Fig. 3a, blue points). This coincidence of strong intrinsic emission and maximum Doppler boosting gives rise to the large outbursts observed at this Θ phase in LS I +61°303.

We likewise examine the relationship between Doppler effects and the minimum of the long-term modulation. One can see the zoom at about 50 900 MJD in Fig. 2d. In the two consecutive outbursts, the model reproduces the periodicity and the lower amplitude of the outbursts. The model data folded in function of the orbital phase are given in Fig. 3b. As shown in Fig. 3d, the angle η between the jet and the line of sight has its minimum value at about Φ = 0.25. This implies that at this Θ phase the maximum Doppler boosting (Fig. 3 f) occurs in an unfavorable situation, i.e., when there is little emission to amplify (low values of S_a and S_b in Fig. 3b).

4.2. Amplitude modulation and orbital shift of the outburst occurrence

During the maximum of the long-term modulation, the peak of the outburst occurs at orbital phase Φ ~ 0.6 (Paredes et al. 1990). Afterwards, the orbital phase of the peak of the outburst changes, a phenomenon analysed in the past by Paredes et al. (1990, with references there) in terms of orbital phase shift, and by Gregory et al. (1999) in terms of timing residuals. In this section we show how our model of a precessing (P₂) jet with periodically (P₁) varying ejection reproduces i) the observed gradual shift of outburst occurrence from Φ = 0.6 to longer orbital phases (Φ ~ 0.9) towards the minimum; ii) the so-called “quiet phase emission with no distictive maximum” discussed in Paredes et al. (1990); and iii) the observed re-appearing of outbursts having a distinctive maximum, even if still low, at orbital phase Φ ~ 0.5.

Figures 3c and d discussed in the previous section show the variations of angle η formed by the jet with the line of sight, vs. orbital phase for the two extreme situations: Θ = 0.83 and Θ = 0.52, which is close to the maximum and the minimum of the long-term modulation. We now consider the curves of η for several different Θ values (Fig. 4 top) and the related outbursts (Fig. 4 bottom).

At Θ = 0.86, the minimum of η that corresponds to the maximum Doppler boosting occurs at Φ = 0.6. In the same Fig. 4 (top) the adopted relativistic electron density is drawn, as a function of P₁, with its maximum at Φ = Φ_r = 0.6. Clearly at Θ = 0.86 the maximum Doppler boosting and the maximum ejection occur at the same orbital phase; that is, the electrons are ejected at the minimum η angle with respect to the line of sight. After 200 days, at Θ = 0.98 the maximum Doppler boosting occurs at Φ = 0.7 (red curve at Θ = 0.98 in Fig. 4 top), and it matches the decay phase of the ejection. The result is the red dotted curve of Fig. 4 (bottom): an outburst slightly lower than that at Θ = 0.86 and shifted at Φ = 0.7. At Θ = 0.1 we can appreciate the evolution of the same phenomenon: the maximum Doppler boosting (MDB; green curve in Fig. 4 top) is further shifted in phase and this causes a shifted and lower outburst (green curve in Fig. 4 bottom). In the interval starting after Θ = 0.265 (with the peak of the outburst at Φ ≃ 0.85) until Θ ~ 0.45, the outburst loses its typical shape, which is now a rather broad low outburst, and it becomes very difficult to define a peak (see the curve at Θ = 0.45 in Fig. 4 bottom). This is the so-called “quiet phase emission, with no distinctive maximum” noticed by Paredes et al. (1990) during the minimum of the long-term modulation. Finally, the outburst starts to resume its shape at Θ = 0.52. As one can see in Fig. 4, there is some boosting at the onset of the new ejection peaking at Φ = 0.6. The result is an outburst with peak at Φ ≃ 0.5. Our model, therefore, is able to reproduce the orbital phase shift of Fig. 5 of Paredes et al. (1990).

4.3. Orbital and outburst periodicities

Figure 2b, where the fits of two consecutive large outbursts are shown, has already proved the ability of the model also reproducing the shorter orbital periodicity, P₁ = 26.49 d. Here the comparison is extended to the whole data set, by folding, as shown in Fig. 5a, all 6.7 yr of GBI data and model data with P₁. However, P₁ refers to the orbital periodical ejection of relativistic electrons, the real periodicity of the outburst is P_average = 26.70 d (Sect. 1). Because P_average is the average of P₁ and P₂, one would expect, when folding the data with P_average, a clustering similar to or somehow intermediate to what is obtained using P₁ and P₂ (shown in Fig. 5b). In contrast, Fig. 5c reveals a double clustering. In Fig. 5c, we use a different colour for model data and observations before or after the minimum of the long-term modulation. The GBI points before (red points) and after (green points) 50 841 MJD cluster separately with a shift of 0.5 in phase. This effect is discussed in Massi & Jaron (2013) and Jaron & Massi (2014). What we want to point out here is that our model indeed reveals the same double clustering as the GBI observations: model data before (black points) and after (blue points) 50 841 MJD cluster separately with a shift of 0.5 in phase as well.

4.4. Spectral index

The top panel of Fig. 6 shows the observed spectral index, from GBI data at ν = 2.2 GHz and 8.3 GHz, folded with the 1667 d long-term modulation. The GBI data indicate that at some Θ phases, the emission in LS I +61°303 develops an optically thick spectrum (α ≥ 0). This is the flat or inverted spectrum typical of jets in microquasars as discussed in Sect. 2.1. Our model of a precessing conical jet with a periodical increase in intrinsic emission must therefore be able to reproduce this spectral characteristic of LS I +61°303 emission.

To make the comparison, we calculated the observed flux, S_model, at 2.2 GHz just by replacing this frequency in equations shown in Sect. 2 and using the parameters obtained by fitting the model to GBI data at 8.3 GHz, i.e. the set of parameters used to generate all figures. Indeed the model spectral index shown at the bottom of Fig. 6 is consistent with the observed one (Fig. 6 top). In particular, we note that the emission becomes optically thick at Θs where the maximum emission of relativistic electrons occurs at the smallest observing angles. In this case, in fact, high values of the optical depth are obtained and the emission becomes optically thick (see Sect. 2.3). A deep analysis of the evolution of the opacity in the jet vs. Φ and vs. Θ is the argument of a paper in preparation. Whereas the model presented here consists of only a conical jet, the following paper will include the transient jet as well. In microquasars, the so-called transient jet always occurs after an optically thick phase (i.e., after the inverted/flat spectrum phase) and is associated to a large optically thin outburst (Fender et al. 2004; Massi 2011). In LS I +61°303, the large optically thin outburst occurs about five days after the optical thick outburst and dominates at 2.2 GHz (Massi & Kaufman Bernadó 2009). The optically thin outburst suddenly reverses α from positive values to the low value α ≃  − 0.5, and this corresponds, as discussed in Massi & Kaufman Bernadó (2009, Sect. 5.1), to a new population of energetic electrons with an index power law p = 2, which is steeper than the one associated to the conical jet, modelled here, with p = 1.8. It is rather interesting that existing analyses of the origin of high energy emission in LS I +61°303 invoke not only two populations of electrons, but also populations with almost the same index as results from radio analysis (p = 1.8 for the conical jet and p = 2 for the transient jet). Zabalza et al. (2011) exclude GeV and TeV emission in LS I +61°303 being produced from the same parent particle population and determine an index ≤1.8 for Fermi data and 2.1 for MAGIC data.

4.5. VLBA maps vs. model

Hints of a fast precessing jet in LS I +61°303 arose when two consecutive MERLIN images revealed a surprising variation of 60° in position angle in only one day (Massi et al. 2004). A set of VLBA observations, A-J, interlapsed three days and covering an interval of 30 days were performed by Dhawan et al. (2006). The VLBA images, about 5 mas in size (therefore more than a factor 10 the orbital size, 2 × a = 0.45 mas) were showing highly varying features. With our model we can determine the variation in the ejection angle during the VLBA runs and compare these variations with the observed variations in position angle of the maps. The changes of η are due to precession, but we put them as usual as a function of the orbit, that is of Φ(P₁), and relate them to changes in VLBA maps also given as a function of Φ(P₁) in Dhawan et al. (2006) and Massi et al. (2012). The trend of η vs. Φ(P₁) is shown in Fig. 7 together with the high dynamic range VLBA images by Massi et al. (2012).

In Fig. 7 we see, as the plot of η vs. Φ reveals, that run H was performed at the minimum η, and D was performed nearly at the maximum η. These two extreme situations of η for the two runs, H and D, are sketched in Fig. 7 (right corner). This important information implies that two runs at similar η, but one performed before and the other after run D, refer to two jets pointing in opposite directions with respect to the axis of the precession cone (Fig. 7, left corner). If the model is correct, the position angle of the associated radio structures should reflect this different orientation. Indeed, the structure for run E points towards the SW, whereas runs C and J show structures pointing to SE. Similarly, the structure at run B points E, whereas that for run F points W. Finally, for A, I, H, G, our model results in a low η angle, and this corresponds to a jet pointing closest to the earth. One can see that the related VLBA structures develop indeed a N-S feature. In run A (G), a transition of the geometry between B (F) and H, I is even traceable.

4.6. Derived jet parameters

The model fitting of Sect. 3 determines the values of two normalization factors present in our model, i.e. the normalization factor for the jet emission $\mathit{Q}\mathrm{=}\frac{\mathrm{2}{\mathit{x}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{r}}_{\mathrm{0}}{\mathit{J}}_{\mathrm{0}}}{\mathrm{-}\mathrm{(}\mathrm{2}\mathrm{-}\mathrm{2}{\mathit{a}}_{\mathrm{3}}\mathrm{-}\mathrm{(}\mathit{p}\mathrm{+}\mathrm{2}\mathrm{\left)}{\mathit{a}}_{\mathrm{2}}\mathrm{\right)}{\mathit{D}}^{\mathrm{2}}}\mathit{,}$(16)and the normalization factor of the jet optical depth $\mathit{q}\mathrm{=}\mathrm{(}\mathrm{2.2}\mathrm{\times }{\mathrm{10}}^{\mathrm{9}}{\mathrm{)}}^{\mathrm{-}\mathrm{(}\mathit{p}\mathrm{+}\mathrm{4}\mathrm{)}\mathit{/}\mathrm{2}}\frac{\mathrm{2}{\mathit{\chi }}_{\mathrm{0}}{\mathit{r}}_{\mathrm{0}}}{\mathrm{-}{\mathit{\kappa }}_{\mathrm{0}}\mathrm{(}\mathrm{2}\mathrm{-}\mathrm{2}{\mathit{a}}_{\mathrm{3}}\mathrm{-}\mathrm{(}\mathit{p}\mathrm{+}\mathrm{2}\mathrm{\left)}{\mathit{a}}_{\mathrm{2}}\mathrm{\right)}}$(17)where r₀ = x₀  tan  ξ.

Taking expressions (12) and (13) into account and assuming p = 1.8 and D =   2.3 Kpc (Gregory et al. 1979), we can express the constant quantities Q and q in terms of the jet physical parameters $\begin{array}{ccc}\mathit{Q}& \mathrm{=}& \frac{{\mathit{x}}_{\mathrm{0}}^{\mathrm{3}}\mathit{tg\xi }{\mathit{B}}_{\mathrm{0}}^{\mathrm{1.4}}\mathit{A}}{\mathrm{-}\mathrm{(}\mathrm{2}\mathrm{-}\mathrm{2}{\mathit{a}}_{\mathrm{3}}\mathrm{-}\mathrm{(}\mathit{p}\mathrm{+}\mathrm{2}\mathrm{\left)}{\mathit{a}}_{\mathrm{2}}\mathrm{\right)}}\mathrm{4.4}\mathrm{\times }{\mathrm{10}}^{-50}\\ \mathit{q}& \mathrm{=}& \frac{{\mathit{x}}_{\mathrm{0}}\mathit{tg\xi }{\mathit{B}}_{\mathrm{0}}^{\mathrm{1.9}}}{\mathrm{-}\mathrm{(}\mathrm{2}\mathrm{-}\mathrm{2}{\mathit{a}}_{\mathrm{3}}\mathrm{-}\mathrm{(}\mathit{p}\mathrm{+}\mathrm{2}\mathrm{\left)}{\mathit{a}}_{\mathrm{2}}\mathrm{\right)}}\mathrm{3.8}\mathrm{\times }{\mathrm{10}}^{-3}\end{array}$and derive, for the solution of Sect. 3, the values for the magnetic field intensity at the jet basis, B₀, and that for the normalization distance x₀:

${\mathit{B}}_{\mathrm{0}}\mathrm{\simeq }\mathrm{9}\mathrm{G}{\mathit{x}}_{\mathrm{0}}\mathrm{\simeq }\mathrm{1.6}\mathrm{\times }{\mathrm{10}}^{\mathrm{13}}\mathrm{cm}\mathit{.}$An independent estimate for the magnetic field strength is the value of the equipartition B_eq ≃ 0.7 G derived from VLBI data (Massi et al. 1993) for a source size d ≃ 5.5    ×    10¹³ cm, assuming equipartition between high energy particle energy and magnetic energy. To compare the result of B₀ ≃ 9 G for our model to B_eq, we compute the jet mean magnetic field $\begin{array}{c}{\mathrm{˜}}_{}\\ \mathit{B}\end{array}$ over its length L. When considering the magnetic field functional form (Sect. 2.1), the topology resulting in Sect. 3 and the value of L (Sect. 2.3), from the expression $\begin{array}{c}\mathrm{˜}\\ \mathit{B}\end{array}\mathrm{=}\frac{{\mathit{B}}_{\mathrm{0}}}{\mathit{L}}{\mathrm{\int }}_{\mathrm{1}}^{\mathit{L}}{\mathit{l}}^{\mathrm{-}{\mathit{a}}_{\mathrm{2}}}\mathrm{d}\mathit{l,}$(20)we obtain $\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{B}\end{array}\mathrm{\simeq }\mathrm{0.6}$G, which is consistent with the value of B_eq derived from the VLBI observation.

From the parameters deduced above, the electron energy content can also be estimated. Since synchrotron emission is centred on a peak spectral frequency ν (Ginzburg & Syrovatski 1965) such as $\mathit{\nu }\mathrm{=}\mathrm{1.8}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{\mathrm{6}}{\mathit{\gamma }}^{\mathrm{2}}\mathit{B,}$(21)the electron energy, E = γmc², can be derived for different values of the magnetic field strength along the jet and for the two frequencies, ν = 2.2    ×    10⁹  Hz ÷ 8.3    ×    10⁹ Hz, considered in this analysis. At the jet basis (l = 1), B = B_o, and γ ≃ 11 ÷ 22 results. For the jet end (l = L), we must take both the decreased magnetic field intensity and relativistic electrons adiabatic losses into account (synchrotron losses are negligible) so that an electron with γ(l = 1) at the jet basis arrives at the jet end with γ(L) = γ(1)L^− 2/3 (Kaiser 2006). To reproduce the observed emission for ν = 2.2    ×    10⁹  Hz ÷ 8.3    ×    10⁹ Hz at the jet end, we need $\mathit{\gamma }\mathrm{\left(}\mathit{L}{\mathrm{)}}^{\mathrm{2}}\mathit{B}\mathrm{\right(}\mathit{L}\mathrm{)}\mathrm{=}\mathit{\gamma }\mathrm{(}\mathrm{1}{\mathrm{)}}^{\mathrm{2}}{\mathit{B}}_{\mathit{o}}{\mathit{L}}^{\mathrm{-}\mathrm{2}\mathit{/}\mathrm{3}\mathrm{-}\mathrm{2}}\mathrm{=}\mathit{\nu }\mathit{/}\mathrm{1.8}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{\mathrm{6}}\mathit{,}$(22)i.e., for the solution of Sect.3, at least electrons with γ(1) ≃ 371 ÷ 741.

In conclusion, relativistic electrons with energies in the range 10  mc² ≲ E ≲ 740  mc² are needed to reproduce the radio emission analysed in this paper. Since the mean over the orbital period of the adopted angular expression (the one shown in Fig. 4 top, see Sect. 2.4) is 0.58  A, the luminosity injected in the two jets in the form of relativistic electrons results in ${\mathit{L}}_{\mathrm{rel}}\mathrm{=}\mathrm{2}\mathrm{\left(}\mathit{\pi }{\mathit{r}}_{\mathit{o}}^{\mathrm{2}}\mathit{c}\mathrm{\right)}\mathrm{0.58}\mathit{A}{\mathrm{\int }}_{{\mathit{E}}_{\min }}^{{\mathit{E}}_{\mathrm{Max}}}{\mathit{E}}^{\mathrm{1}\mathrm{-}\mathit{p}}\mathrm{d}\mathit{E}\mathrm{\simeq }\mathrm{5}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{\mathrm{34}}\mathrm{erg}\mathit{/}\mathrm{s}\mathit{.}$(23)This value is only an indicative lower limit for the luminosity in the injected relativistic electron distribution since it is strictly related to the specific part of the radio spectrum analysed in this paper. In particular the presence of synchrotron emission at higher frequencies will imply higher values for electron energies and, consequently, higher values of L_rel. Indeed, jets from X-ray binaries are supposed to radiate from radio to optical wavelengths. The jet origin for the emission is supported by the correlation of the emission in the different bands (Russell et al. 2013, and reference therein). In this respect, it is worth noting that in LS I +61°303 the timing analysis of optical data by Zaitseva & Borisov (2003) results in a period equal to P₁. This result confirms what has previously been reported by Mendelson & Mazeh (1989, 1994) and Paredes et al. (1994). Particularly intriguing is that the light curve folded with this period shows that the brightness increases in the orbital phase interval where the radio outburst occurs. Finally, as already pointed out in Sect. 4.4, the present analysis holds true for the conical jet with p = 1.8 considered here. Shocks in a transient jet may produce even more accelerated particles with a different energy power-law distribution.

4.7. The low/hard X-ray state

The varying morphology of VLBA images, plus the lack of accretion disk features (e.g., black body radiation) in the X-ray spectrum of LS I +61°303, have led some authors to favour non-accretion models (e.g., Dubus 2013). VLBA images were discussed in Sect. 4.5, and in this section we analyse the X-ray state of LS I +61°303. In fact, an accreting object presents rather different X-ray states (McClintock & Remillard 2006). While the requested black body component is the characteristic of the thermal state (McClintock & Remillard 2006), radio and X-ray observations of LS I +61°303 point towards a very low low/hard X-ray state.

As shown by Corbel et al. (2013), in the early phase of jet formation, the low density of particles in the jet produces an optically thin synchrotron power-law spectrum; as the density of the jet plasma increases, this leads to a transition to higher optical depth and to a flat/inverted radio spectrum (α ≥ 0). When sources having these self-absorbed compact radio jets are observed in X-rays, they reveal a low/hard X-ray state (Gallo et al. 2003; Fender et al. 2004). Concerning LS I +61°303, it is known that its persistent radio emission develops a flat/inverted radio spectrum at some orbital phases (exactly where an increased accretion is predicted) as shown by GBI observations at 2.2 GHz and 8.3 GHz (Massi & Kaufman Bernadó 2009). In addition, a flat spectrum has been recently observed by the Effelsberg 100-m telescope at four frequencies, from 4.85 GHz to 14.6 GHz (Fig. 4.13, Zimmermann 2013).

Modelling of the low/hard state, especially at very low luminosities, favours an accretion disk truncated at a large radius (Fender 2001; Gardner & Done 2013, and references therein). The origin of the X-ray power-law emission with a typical photon index Γ ~ 1.7 (Remillard & McClintock 2006) (for LS I +61°303 in the range 1.5−1.8, Sidoli et al. 2006; Chernyakova et al. 2006; Zimmermann & Massi 2012) is still controversial, if it is due to a Comptonizing corona or/and to the jet (Markoff et al. 2001). Nevertheless, it is well established that this X-ray power-law emission that characterizes the low/hard state is correlated with the radio emission from the self-absorbed jet (${\mathit{L}}_{\mathit{R}}\mathrm{\propto }{\mathit{L}}_{\mathrm{X}}^{\mathrm{0.7}}$, Gallo et al. 2003). The luminosity values for LS I +61°303 are L_X = (1 − 6) × 10³⁴ ergs/s and L_R = (1 − 17) × 10³¹ erg/s (Combi et al. 2004). From Eqs. (3) and (5) in Merloni et al. (2003), we obtain ${\mathit{\xi }}_{\mathit{R}\mathrm{X}}\mathrm{=}\frac{\log {\mathit{L}}_{\mathit{R}}\mathrm{-}{\mathit{\xi }}_{\mathrm{RM}}\log \mathit{M}\mathrm{-}{\mathit{b}}_{\mathrm{R}}}{\log {\mathit{L}}_{\mathrm{X}}}$(24)where ξ_RM = 0.69 − 0.89 and b_R = 3.26 − 11.38 (Merloni et al. 2003). For a mass of 3 M_⊙ (where the influence of the mass difference between a neutron star of 1.4 M_⊙ and a stellar mass black hole of 3 − 15  M_⊙ is negligible in the result), we obtain ξ_RX   = 0.55 − 0.85. This value is consistent with the solutions ξ_RX = 0.49 − 0.71 by Merloni et al. (2003) and with the value of ~0.7 by Gallo et al. (2003) for accreting black holes in low/hard state, and it is clearly below the value >1.4 given by Migliari & Fender (2006) for accreting neutron stars in the low/hard state. Even though the nature of the compact object in LS I +61°303 is still unknown (Casares et al. 2005), the agreement with the low/hard state of black hole X-ray binaries would indeed favour the black hole scenario.

The luminosity for accreting black holes in the low/hard state is about L_X ~ 10³⁶  erg/s that may drop to L_X = 10^30.5 − 10^33.5 erg/s (McClintock & Remillard 2006) at the lowest phase, which is called the quiescent state. With its luminosity L_X = (1 − 6) × 10³⁴ erg/s, LS I +61°303 is clearly in a very low low/hard X-ray state. That is a similar case to MWC 656, the recently discovered Be-black hole system, which is in the quiescent state (Casares et al. 2014). During the low-hard state the liberated energy of the accretion is thought to be converted into magnetic energy that powers the relativistic jet typically observed in this state (Smith 2012, and references therein). In a system in quiescence, Gallo et al. (2006) demonstrate that the outflow kinetic power accounts for a sizable fraction of the accretion energy budget.

5. Conclusions

In this paper we have developed a model that reproduces the observed long-term modulation of the radio emission of the binary system LS I +61°303. The hypothesis is that the radio emission is due to a jet precessing with a period of . The jet is anchored on the accretion disk of a compact object orbiting, with a period of , the Be star. Variable accretion around the eccentric orbit increases the relativistic electron density in the jet at a particular orbital phase. We fitted our theoretically derived flux density, S_model(t) (see Eq. (1)) to 6.5 yr GBI radio observations and determined the parameters of Table 1. We found that the precessing (P₂) jet model with intrinsic periodical variations (P₁) of emitting particles reproduces all seven observational facts:

• 1.

  the long-term modulation (Gregory 2002, andreferences there), see Figs. 2aand 4 bottom;

• 2.

  the orbital shift of the peak of the radio outburst (Paredes et al. 1990), see Fig. 4 bottom;

• 3.

  the periodical outburst with orbital occurrence (Gregory 2002, and references there), see Figs. 2b and c and Fig. 5a;

• 4.

  the double clustering for data folded with $\frac{{\mathit{P}}_{\mathrm{1}}\mathrm{+}{\mathit{P}}_{\mathrm{2}}}{\mathrm{2}}$ (Massi & Jaron 2013), see Fig. 5c;

• 5.

  the spectral index evolution (Massi & Kaufman Bernadó 2009), see Fig. 6;

• 6.

  a variation of the ejection angle suggestive of the rapid rotation in position angle of VLBA maps (Dhawan et al. 2006; Massi et al. 2012), see Fig. 7;

• 7.

  the strength of the magnetic field (Massi et al. 1993), see Sect. 4.6.

The maximum of the long-term modulation occurs when the maximum ejection of relativistic electrons takes place and at the same time the approaching jet is forming the smallest possible angle with the line of sight (i.e., the Doppler boosting is at its maximum). After a cycle of duration P₂, when the approaching jet again forms the smallest angle η, the ejection is now already on its decay phase (being P₁ < P₂) and its Doppler boosting gives rise to an outburst with a lower amplitude than that of one cycle before. As we quantitatively show in Sect. 4.2, the peak of the outburst gets lower and lower at each cycle until the most unfavorable situation occurs, i.e., when the jet is pointing towards us while the main ejection is at its minimum. After this point, it happens that when the approaching jet again forms the smallest η angle, there is the beginning of a new ejection of particles and the observed flux density slowly rises again. The period of 1667   ±   8 days (Gregory 2002) of the long-term modulation is then just the number of P₁ cycles needed to compensate for the difference P₂ − P₁ and synchronize P₁ again to P₂, which is just [P₂/(P₂ − P₁)] P₁ (Massi & Jaron 2013).

Our model, meant in the first instance to explain the long-term modulation affecting the radio emission of LS I +61°303, can become a powerful tool for understanding the physical processes at work in this gamma-ray binary. For this purpose the model should be extended to include internal shocks in the conical flow, i.e., the transient jet. As discussed in Sect. 4.4. the transient jet observed in LS I +61°303 is associated to a second population of relativistic electrons, with p = 2.0, which is therefore different from the population with p = 1.8, forming the conical jet analysed in this paper. It is quite interesting that two similar populations are invoked to explain emission in the TeV and GeV band (Zabalza et al. 2011). The important implications, if confirmed, are that the GeV emission is related to inverse Compton of the electrons of the conical jet and that the TeV emission is related to electrons of the transient jet.

This model, with conical jet and transient jet can be tested by fitting it to the radio spectral index data analysed in Massi & Kaufman Bernadó (2009). The theory of accretion in an eccentric orbit (Sect. 2.4) also predicts an ejection around periastron passage, along with the ejection at orbital phase Φ = 0.6 considered here (i.e., towards apoastron). By fitting the model to the radio spectral index data, one could verify if each of the two populations of electrons is generated twice along the orbit, i.e., that the transient jet occurs at periastron as well.

Indeed, at different epochs, emission towards apoastron (Albert et al. 2009) and periastron (Acciari et al. 2011) was detected by Cherenkov telescopes. Different Doppler boosting of the two ejections along the eccentric orbit could explain the observed variations. Maximum Doppler boosting occurring between the two ejections could explain the broad gamma-ray emission without any clear orbital modulation observed in Fermi-LAT data (Hadasch et al. 2012; Ackermann et al. 2013).

Finally, considering that one contribution to the gamma-ray emission consists of upscattered stellar photons, the role of Doppler boosting in gamma-ray emission is probably even more important than in the radio band. In fact, the amplification due to the Doppler factor for external inverse Compton is higher than for synchrotron emission (Kaufman Bernadó et al. 2002, and reference therein).

Acknowledgments

We thank Brenda Burton, Lars Fuhrmann, Frédéric Jaron, Jürgen Neidhofer, and Lisa Zimmermann for interesting discussions and the anonymous referee for the careful reading of our manuscript and the useful comments made for its improvement. The Very Long Baseline Array is operated by the USA National Radio Astronomy Observatory, which is a facility of the USA National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Green Bank Interferometer is a facility of the National Science Foundation operated by the NRAO in support of NASA High Energy Astrophysics programmes. MINUIT is the Function Minimization and Error Analysis CERN Program Library. We thank Dirk Muders for technical support.

References

Appendix A: Optical depth through a stratified jet

The integration upper limit in the intensity expression in Eq. (7), τ_end(l), represents the jet’s maximum optical depth. Figure A.1 shows a cut of the jet structure where the line of sight makes an angle η with the jet axis, l. When looking at the figure it is easy to understand how photons coming out from the jet side facing the observer at generic distance x₀l originate along the integration path inside the jet, i.e., along the segment AB for the case of the approaching jet shown in the figure.

If we call l_A the value of the axial distance corresponding to point A and l_B that corresponding to point B, the following relationships describe the geometrical projection of the optical path AB on the jet axis and on the jet radius, respectively: $\begin{array}{ccc}\mathit{AB}\cos \mathit{\eta }& \mathrm{=}& {\mathit{x}}_{\mathrm{0}}\mathrm{(}{\mathit{l}}_{\mathit{A}}\mathrm{-}{\mathit{l}}_{\mathit{B}}\mathrm{)}\\ \mathit{AB}\sin \mathit{\eta }& \mathrm{=}& \mathit{r}\mathrm{\left(}{\mathit{l}}_{\mathit{A}}\mathrm{\right)}\mathrm{+}\mathit{r}\mathrm{\left(}{\mathit{l}}_{\mathit{B}}\mathrm{\right)}\mathrm{=}{\mathit{r}}_{\mathrm{0}}\mathrm{(}{\mathit{l}}_{\mathit{A}}\mathrm{+}{\mathit{l}}_{\mathit{B}}\mathrm{)}\mathit{.}\end{array}$The jet region included between l_A and l_B, i.e., l_A ≥ l ≥ l_B, is the one that contributes to the emission coming out of the jet at l = l_A. It is evident from the figure that the emitting region is different for different values of l_A. To quantify the depth of the involved region, the value of l_B can easily be derived from the above expressions and results in ${\mathit{l}}_{\mathit{B}}\mathrm{=}{\mathit{l}}_{\mathit{A}}\frac{\tan \mathit{\eta }\mathrm{-}\tan \mathit{\xi }}{\tan \mathit{\eta }\mathrm{+}\tan \mathit{\xi }}\mathrm{·}$(A.3)All the above relationships hold in the limit that the angle η, between the jet axis and the line of sight, is larger than the jet opening angle ξ = r₀/x_o; i.e., tan  η > tan  ξ.

Once l_B is known following Eq. (8), the optical depth at the distance l_A is ${\mathit{\tau }}_{\mathit{A}}\mathrm{-}{\mathit{\tau }}_{\mathit{B}}\mathrm{=}\mathrm{-}{\mathrm{\int }}_{{\mathit{l}}_{\mathit{B}}}^{{\mathit{l}}_{\mathit{A}}}{\mathit{\chi }}_{\mathit{\nu }}\mathrm{d}\mathit{x}\mathrm{=}{\mathit{\tau }}_{\mathrm{0}}{\mathit{l}}_{\mathit{A}}^{\mathit{C}}\mathrm{[}\mathrm{1}\mathrm{-}\mathrm{(}{\mathit{l}}_{\mathit{B}}\mathit{/}{\mathit{l}}_{\mathit{A}}{\mathrm{)}}^{\mathit{C}}\mathrm{]}$(A.4)where C = 1 − a₃ − (p + 2)a₂/2 and τ_A = 0 assuming no emission or absorbtion between the observer and the jet.

Since l_A represents our generic distance along the jet axis, we can define for the approaching jet the optical depth τ_end ≡ τ_B as $\mathrm{\left[}{\mathit{\tau }}_{\mathrm{end}}{\mathrm{]}}_{\mathrm{app}}\mathrm{\right(}\mathit{l}\mathrm{)}\mathrm{=}{\mathit{\tau }}_{\mathrm{0}}{\mathit{l}}^{\mathit{C}}\left[{\left(\frac{\tan \mathit{\eta }\mathrm{-}\tan \mathit{\xi }}{\tan \mathit{\eta }\mathrm{+}\tan \mathit{\xi }}\right)}^{\mathit{C}}\mathrm{-}\mathrm{1}\right]\mathrm{·}$(A.5)For the receding jet the same procedure can be performed by paying attention to some different signs (note that also for the receding jet the positive direction for l is defined going away from the core), to arrive at $\mathrm{\left[}{\mathit{\tau }}_{\mathrm{end}}{\mathrm{]}}_{\mathrm{rec}}\mathrm{\right(}\mathit{l}\mathrm{)}\mathrm{=}{\mathit{\tau }}_{\mathrm{0}}{\mathit{l}}^{\mathit{C}}\left[\mathrm{1}\mathrm{-}{\left(\frac{\tan \mathit{\eta }\mathrm{+}\tan \mathit{\xi }}{\tan \mathit{\eta }\mathrm{-}\tan \mathit{\xi }}\right)}^{\mathit{C}}\right]\mathrm{·}$(A.6)

All Tables

All Figures