Long-term periodicity in LSI+61303 as beat frequency between orbital and precessional rate
aa r X i v : . [ a s t r o - ph . H E ] M a r Astronomy&Astrophysicsmanuscript no. massijaron c (cid:13)
ESO 2017December 13, 2017
Long-term periodicity in LS I +61 ◦
303 as beat frequency betweenorbital and precessional rate
M. Massi and F. Jaron
Max-Planck-Institut f¨ur Radioastronomie, Auf dem H¨ugel 69, D-53121 Bonn, Germanye-mail: mmassi, fjaron, @mpifr-bonn.mpg.de
Received 2012;
ABSTRACT
Context.
In the binary system LS I + ◦
303 the peak flux density of the radio outburst, which is related to the orbital period of26 . ± . ± −
28 d.
Aims.
How close is the precessional period of the radio jet to the orbital period? Any periodicity in the radio emission should berevealed by timing analysis. The aim of this work is to establish the accurate value of the precessional period.
Methods.
We analyzed 6.7 years of the Green Bank Interferometer database at 2.2 GHz and 8.3 GHz with the Lomb-Scargle andphase dispersion minimization (PDM) methods and performed simulations.
Results.
The periodograms show two periodicities, P = . ± .
07 d ( ν = . − ) and P = . ± .
07 d ( ν = . − ).Whereas radio outbursts have been known to have nearly orbital occurrence P with timing residuals exhibiting a puzzling sawtoothpattern, we probe in this paper that they are actually periodical outbursts and that their period is P average = ν + ν = . ± .
05 d.The period P average as well as the long-term modulation P beat = ν − ν = ±
393 d result from the beat of the two close periods, theorbital P and the precessional P periods. Conclusions.
The precessional period, indicated by the astrometry to be of 27–28 d, is P = .
92 d. The system LS I + ◦
303 seemsto be one more case in astronomy of beat, i.e., a phenomenon occurring when two physical processes create stable variations of nearlyequal frequencies. The very small di ff erence in frequency creates a long-term variation of period 1 / ( ν − ν ). The long-term modulationof 1667 d results from the beat of the two close orbital and precessional rates. Key words.
Radio continuum: stars - X-rays: binaries - X-rays: individual (LS I + ◦
303 ) - Gamma-rays: stars
1. Introduction
The TeV emitting source LS I + ◦
303 has radio characteris-tics that make it unique not only among the small number ofgamma-ray emitting systems given in the X-ray binary class,a class of binary systems where a neutron star or a black holeis orbiting around a normal star, but also among the largergroup of the radio emitting X-ray binary systems (Fender et al.1997; Massi 2005; Mirabel 2012). The peak flux density ofthe radio outburst, which is related to the orbital period of26 . ± . ± ff erent spectral characteristics. There is afirst outburst with a flat / inverted spectrum and a second opti-cally thin outburst associated with di ff erent conditions, as in-dicated by its high amplitude, the spectral index, and the H α emission line measurements (Massi & Kaufman Bernad´o 2009;Grundstrom et al. 2007). The complex spectral sequence foundin LS I + ◦
303 finds a natural explanation in the disk-jet cou-pling model for microquasars: first, there is a continuous outflowwith a flat or inverted spectrum, then an event triggers a shockin this slow optically thick outflow (Fender et al. 2004), and thegrowing shock creates the optically thin outburst (Valtaoja et al.1992; Hannikainen et al. 2006). One of the characteristics thatmake LS I + ◦
303 unique among the other radio emitting X-ray binary systems is that this spectral evolution, between in-verted and optically thin spectra, may occur twice during the or-bital period (Massi & Kaufman Bernad´o 2009). This agrees well with the Bondi & Hoyle (1944) accretion in an eccentric orbit(as in LS I + ◦ + ◦
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 + ◦ +
105 (Pooley & Fender 1997;Rodriguez & Mirabel 1997, 20–40 min).The radio morphology of the system also shows unique char-acteristics. The resolved extended structure changes position an-gle (i.e., angle of the projection of the jet onto the plane ofthe sky, measured from north through east) with the surpris-ing large variation of 60 ◦ in only one day (Massi et al. 2004;Dhawan et al. 2006). Moreover, the jet is sometimes one-sidedand at other times two-sided. Because of both of these vari-ations in position angle and morphology, the hypothesis thatLS I + ◦
303 might be a precessing microquasar was broughtforth (Massi et al. 2004). The one-sidedness of jets is usually at-tributed to relativistic bulk motion along a relatively small angleto the line of sight, which leads to Doppler boosting of the jet anddeboosting of the counterjet emission (Urry & Padovani 1995).A variation of that angle due to precession would cause vari-able Doppler deboosting of the counter jet, making it appear for
1. Massi and F. Jaron: Orbital and precessional periodicities in LS I + ◦ larger angles (double-sided jet) and disappear for smaller valuesof the angle to the line of sight (one-sided jet, blazar like).In 2006, VLBA observations by Dhawan et al. (2006) mea-sured the same large rotation of 60 ◦ / day in their images asMassi et al. (2004). Some of the VLBA images, showing again aone-sided structure were, however, interpreted by Dhawan et al.(2006) as a cometary tail pointed away from the companion Bestar, that is in favor of a pulsar model rather than of the pre-cessing microquasar model. The reanalysis of this VLBA dataset and the resulting higher dynamic range of the self-calibratedmaps has actually revealed a double-sided structure in severalimages (Massi et al. 2012). Before we illustrate how the re-sults in Massi et al. (2012) brought us to the present investiga-tion on the precessional period, let us consider two importantpoints in the two following paragraphs, the first concerning self-calibration and the second the pulsar model.Self-calibration of interferometric data is a well-establishedtechnique (Cornwell & Fomalont 1999). It may fail at SNR < σ are present in 6 out the 12 images. Concerning the e ff ectson an image of strong variations of flux density during observa-tions (Stewart et al. 2011), the source LS I + ◦
303 shows a lowradio flux density at all orbital phases, apart from the maximumof the long-term modulation. In these epochs a large outburstlasting few days occurs around apastron. This means that dur-ing the maximum of the long-term modulation one might expecta reduction of the dynamic range of the produced maps aroundapastron, i.e., weak features will be lost. This is not the case forthe VLBA observations of Dhawan et al. (2006) performed to-ward the minimum of the long-term modulation.We are therefore faced in LS I + ◦
303 with a changingstructure from a double-sided structure to a one-sided struc-ture. Are such variations possible also in the pulsar scenario?Simulations (Mold´on et al. 2012) show that the emission fromthe cometary tail of a pulsar for a particular orientation and incli-nation of the orbit, after almost one orbital cycle and at a partic-ular orbital phase may look like a double-sided nebula at a fixedposition angle. In these conditions, the fading and expanding lastpart of the cometary tail may appear detached from the bright-est part of the cometary tail which is closer to the orbit. Thispossibility is clearly ruled out for LS I + ◦ ff erent orbital phases and at di ff erent position angles(Massi et al. 2012).We return thus to the scenario of a precessing microquasar,where the two-sided jet su ff ers from variable Doppler boost-ing because the precession continuously changes the angle be-tween jet and line of sight (Massi et al. 2004, 2012). Derivingthe precessional period from the available radio images is notstraightforward, because the images reflect the variation of theprojected angle on the sky plane and, therefore, a combinationof the ejection angle, inclination, and angle of the precessioncone. The astrometry provides less biased results (Dhawan et al.2006; Massi et al. 2012). The astrometry of the VLBA obser-vations indicates that the peak in consecutive images describesa well-defined ellipse, six to seven times larger than the orbit,with a period of about 27-28 d (Massi et al. 2012). Assuming forthe microquasar in LS I + ◦
303 the same core-shift e ff ect typ-ically observed in blazars, the peak of each image would thencorrespond 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 el-lipse as the possible cross-section of the precession cone of thejet at the distance where the emission at 8.4 GHz becomes op-tically thick. The determined time span of 27-28 d to completethe ellipse is a first estimate of the precession period.The most likely cause for precession of an accretion diskof a compact object is an assymetric supernova explosion ofthe progenitor. As a result the compact object could be tilted(Fragile et al. 2007). In this case either the accretion disk iscoplanar with the compact object and, therefore, subject to thegravitational torque of the Be star or, instead, the accretiondisk is coplanar with the orbit but tilted with respect to thecompact object which induces, in the context of general rela-tivity, Lense-Thirring precession if the compact object rotates(Massi & Zimmermann 2010). A deep investigation of these orother mechanisms of precession requires the knowledge of theprecession parameters, such as the period of precession and theangle of the precession cone.In this paper we present a timing analysis of 6.7 years ofGreen Bank Interferometer (GBI) radio data aimed at a more ac-curate determination of the precession period. Sections 2.1 and2.2 present the determined period, called P . The section illus-trates the case that, while the aim of our research was reached,i.e., we obtained a more precise value of P , we were presentedwith an additional unexpected result: the beating between the or-bital period P and the precessional period P gives rise to a newperiod, P average = / ( ν + ν ) modulated by 1 / ( ν − ν ) = P average the period-icity of the observed radio outburst? Indeed, whereas in the lit-erature P is generally referred to as the period of the radio out-bursts, it is also well known that there are di ff erences betweenthe observed and predicted (for P = P ) outburst times, and thatthese timing residuals have vs time a sawtooth pattern (Sect.2.3). In Sect. 2.4 we present two important results. First, wedemonstrate mathematically that indeed the sawtooth functionadjusts P to P average . Then we show that the GBI data foldedwith P average present an o ff set of 13 d at the minimum of thelong-term modulation equal to the sawtooth function and equalto that predicted by the beat between P and P . In the samesection we also present the corresponding physical scenario. InSect. 3 we discuss the implications of our results for the ob-served periodicity in the equivalent width of the H α emissionline in LS I + ◦
2. Data analysis and results
We analyzed here 6.7 years of the NASA / NRAO GBILS I + ◦
303 database at 2.2 GHz and 8.3 GHz. Thedatabase covers three periods: 49379.975 − − − ffi cient on irreg-ularly sampled data (Lomb 1976; Scargle 1982). We usedthe algorithms of the UK Starlink software package, PERIOD( ). For both data sets at2.2 GHz and 8.3 GHz (Fig. 1) we obtained the same result:two periods P = . ± .
07 d ( ν = . ± . − )and P = . ± .
07 d ( ν = . ± . − ). Thestatistical significance of a period is calculated in PERIODfollowing the method of Fisher randomization as outlined
2. Massi and F. Jaron: Orbital and precessional periodicities in LS I + ◦ Power
Frequency (d -1 )a 0 600 1200 1800 Power
Frequency (d -1 )b 0.6 0.7 0.8 0.9 1 Rel. Variance
Frequency (d -1 )c 0.7 0.8 0.9 1 Rel. Variance
Frequency (d -1 )d Fig. 1.
Periodograms of 8.3 GHz (filled squares) and 2.2 GHz (circles) data. a: output of Lomb-Scargle method. The Lomb-Scargleanalysis gives on the y-axis the significance of the frequency. Two frequencies are found at ν = . − ( P = .
49 d) and at ν = . − ( P = .
92 d). b: output of Lomb-Scargle method for data ≥ σ . c: output of PDM. The most likely period yieldsthe 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 ≥ σ .in Linnel Nemec & Nemec (1985). The advantage of using aMonte-Carlo- or randomization-test is that it is distribution-free and that it is not constrained by any specific noise mod-els (Poisson, Gaussian, etc.). The fundamental assumption isthat if there is no periodic signal in the time series data, thenthe measured values are independent of their observation timesand are likely to have occurred in any other order. One thou-sand randomized time series are formed and the periodogramscalculated. The proportion of permutations that give a peakpower higher than that of the original time series would thenprovide an estimate of p , the probability that, for a given fre-quency window there is no periodic component present in thedata with this period. A derived period is defined as signifi-cant for p < .
01, and a marginally significant one for 0 . < p < .
10 (Linnel Nemec & Nemec 1985). For both periods P = . ± .
07 d (frequency window 0.0374 − − )and P = . ± .
07 d (frequency window 0.0369 − − )and for both data sets at 8.3 GHz and 2.2 GHz we obtained0 . < p < . Figure 1 a shows the results of the Lomb-Scargle analysis for thedata at 8.3 GHz and 2.2 GHz. There P is dominating over P for a factor of 1.8 at 2.2 GHz, and a factor of 1.5 at 8.3 GHz.Figure 1 b shows the results of the Lomb-Scargle analysis, ifonly data with flux density ≥ σ are used. In this case the twoperiods have a more comparable significance, i.e., there is a fac- tor of 1.4 at 2.2 GHz, and a factor of 1.2 at 8.3 GHz. We comparethe Lomb-Scargle results with those obtained with the phase dis-persion minimization (PDM) method (Stellingwerf 1978). Theresults of the PDM analysis on the whole data sets are shownin Fig. 1 c; the results on data ≥ σ are shown in Fig. 1 d. Theresults of the PDM analysis agree very well with those of Lomb-Scargle: there is a di ff erent significance of the two periods whendata with low signal-to-noise-ratio (snr) are present in the ana-lyzed data set, i.e., P dominates.This could be the explanation why in the past the secondperiod P was unseen. Taylor & Gregory (1982) found a pe-riod of 26 . ± .
04 d in their radio data set; in 1984 theywrote that part of the previous measurements were taken withthe source in a weak state and repeated the analysis usingnew data and also including the old ones, obtaining the valueof 26 . ± .
008 d (Taylor & Gregory 1984). In 1997 Ray etal. reported new observations and gave a period of 26 . ± .
02 d, i.e., coincident with our average P average = ν + ν = . + . d = . ± .
05 d, as discussed in Sect. 2.2. Theseobservations (Ray et al. 1997) are those of our first GBI inter-val, i.e., 49379.975 − σ ) from the 26 . ± .
08 d value of Taylor & Gregory(1984). Gregory et al. (1999), faced with the di ff erence betweentheir value and the Ray et al. (1997) results, discuss how unlikely
3. Massi and F. Jaron: Orbital and precessional periodicities in LS I + ◦ a sudden change in period would be. In the light of our presentresult we see that it is not a sudden change in period but thepresence of two periods that in the Ray et al. data have compara-ble significance. As noted above, the value of Ray et al. (1997)26 . ± .
02 d corresponds to our P average , as discussed in thenext section. P average The two frequencies ν = . − ( P = .
49 d) and ν = . − ( P = .
92 d) are only slightly di ff erent. This pro-duces a beating, i.e., a new frequency is formed ν average = ν + ν ,modulated with ν beat = ν − ν . For the sum of two sine functionsthe following identity holdssin (2 πν t ) + sin (2 πν t ) = (cid:18) π ν − ν t (cid:19) sin (cid:18) π ν + ν t (cid:19) , (1)where the beat frequency (or frequency of the envelope) ν beat = ν − ν is twice the frequency of the cosine term. In our case, theterm ν − ν = . − . d is equal to 1667 ±
393 d. Figure 2 c shows the sum of two sine functions with di ff erentamplitudes f b ( t ) = sin (2 πν t ) + a sin (2 πν t ) = a cos (cid:18) π ν − ν t (cid:19) sin (cid:18) π ν + ν t (cid:19) + (1 − a ) sin (2 πν t ) , (2)with a = .
7, whereas Fig. 2 e shows the function f m ( t ) = (1 + b sin(2 πν m t )) sin(2 πν t ) , (3)with b = . ν = . d − , and ν m = d − . As one can see,both Eqs. 2 and 3 are able to reproduce the long-term modula-tion; however, the periodograms are rather di ff erent. The peri-odogram of Eq. 2, shown in Fig. 2 d agrees well with the peri-odogram of the GBI data of Fig. 2 b. On the contrary, as one cansee in Fig. 2 f, in the periodogram of Eq. 3 only P = .
49 d ispresent in the frequency range 0.036 − − . Gregory (1999) demonstrated the existence of a long-term mod-ulation of the peak outburst flux. Gregory & Neish (2002) in-dicated that the modulation in radio properties may stem fromperiodic ejections of a shell (density enhancement) of gas in theequatorial disk of the Be star.The long-term modulation is also present in the timing resid-uals of the outbursts, i.e., the di ff erence between observed andpredicted (for P = P ) outburst time (Gregory et al. 1999). Theobservational result is that timing residuals show a surprisingsawtooth 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. Thesaw tooth function is shown in Figs. 2 and 8 a of Gregory et al.(1999) and here in Fig. 3 a. In detail the trend is as follows: ob-served and predicted outburst times coincide at the peak of thelong-term modulation (i.e., at the radio maximum) resulting ina timing residual equal to τ =
0; the timing residual grows lin-early with time and at the minimum of the long-term modula-tion reaches a maximum of about 6 d where it sharply switchesto about -7 d. Also surprising is that the transition τ ≃ τ ≃ − ( P − P ) − = ±
382 when using ≥ occurs at a time when the amplitude is low, i.e., at the mini-mum of the long-term modulation (Figs. 6 a and 7 b of Gregoryet al. 1999). After this transition, τ starts to grow linearly withtime reaching the value τ = P modulated by 1667 d) and to verify ifthe corrected model is able to reproduce the observed spectrum(i.e., P and P ). First we have to define the sawtooth function.The slope of the sawtooth pattern in Fig. 8 a of Gregory et al.(1999), is about 0.008; we generate therefore the sawtooth func-tion (Fig. 3 a) with a period of 1667 d as τ ( t ) = .
008 fmod ( t , f m ( t ) = (1 + b sin(2 πν m t )) sin(2 πν ( t − τ ( t ))) . (5)To study the e ff ect of the sawtooth function without any biasresulting from large holes in the sampling, we performed bothsimulations with original sampling and with regular sampling.We obtained the same results and here we show those with reg-ular sampling.The resulting periodogram of the simulated data (Fig. 2 g) isshown in Fig. 2 h. One sees that the model with the single peri-odicity P , modulated by 1667 d, once corrected by the sawtoothfunction, is able to reproduce the results of our spectral analysisthat is the two periods P and P . In the next section we showthat when one directly uses the two found periods P and P , thesawtooth function is naturally explained. P average The sawtooth function results from the comparison between ob-served and predicted (for P = P ) outburst time. Here we willshow first analytically (Eq. 6) and then with the GBI observa-tions that the observed outburst occurs at P = P average , i.e., at1 / ν average of Eq. 1.Analytically, adjusting P = .
49 d ( ν = . − ) bythe given sawtooth function τ ( t ) to ν ( t − τ ( t )) = ν t (1 − . = (0 . × . t = . t (6)gives as a result P average = . d = .
70 d as in Sect.2.1, (with variations between 26.65 − − P ) and observed outburstsare equal to the residuals between predicted (at P ) and P average .We therefore ascertained that the observed outburst has peri-odicity P average . However, if we use Eq. 3 and we simply sub-stitute P by P average , the periodogram fails to reproduce theobserved periodogram with P and P and just shows P average (small window on the left in Fig. 2 f). Is P average only an apparentperiodicity, just produced by P and P ?Let us fold the GBI data on the periods P , P , and P average (Figures 4 a-c). The data folded with the orbital period P showthe broad cluster that is well known in the literature (e.g., Fig. 2 cin Massi & Kaufman 2009), where flares above 100 mJy occurfrom about phase 0.4 to about phase 0.9. The broadening is dueto the di ff erences between the observed and predicted (by P )outburst time that also causes the sawtooth pattern. Figure 4 b,shows for the first time the data folded with the precessional fmod is a function implemented in the math library of C.4. Massi and F. Jaron: Orbital and precessional periodicities in LS I + ◦
50 100 150 200 250 300 49600 50400 51200
Flux density (mJy)
Time (MJD)
GBI Observations (8 GHz) a 50 100 150 200 250 300 49600 50400 51200
Flux density (mJy)
Time (MJD) sin(1/26.49) + 0.7 sin(1/26.92) c 50 100 150 200 250 300 49600 50400 51200
Flux density (mJy)
Time (MJD) sin(1/26.49) [1 + 0.7 sin(1/1667)] e 50 100 150 200 250 300
Flux density (mJy)
Time (MJD) model e + sawtooth function g 0 200 400 600 800 1000 1200 1400 1600 0.036 0.038 0.04
Power
Frequency (d -1 )b Power
Frequency (d -1 )d Power
Frequency (d -1 )f Power
Frequency (d -1 )h 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 ν = . − ( P = .
49 d) and ν = . − ( P = .
92 d). The small window shows the peak at1 / − present in all periodograms. c: Sum of two sinusoidal functions at 26.49 d and 26.92 d, with an amplitude ratio 1 / . d − and . d − . The significance of the two frequencies becomes identicalin the periodogram only for an amplitude ratio 1 /
1. e: Long-term modulation (1667 d) of a 26.49 d periodic outburst. f: Lomb-Scargleanalysis results: one frequency at . d − . The small window to the left shows the peak at P average present in the periodogram ofa simulation of long-term modulation (1667 d) of a 26.70 d periodic outburst. g: Sine wave of periodicity P , modulated by a sinewave of periodicity 1667 d and corrected by a sawtooth function. h: Lomb-Scargle analysis results: two frequencies at . d − and . d − as in Figs. 2 b and d.period P . The cluster of the large flares is also evident and su-perimposed to scattered smaller flares. If now we fold the datawith P average one would expect, since it is the average of P and P , a clustering rather similar to that in Fig. 4 a and b. The re-sult, shown in Fig. 4 c, is completely di ff erent. First of all, thereare two clusters. Second, each one of the two clusters is not asbroad as those with P and P , i.e., the clustering is better. Wheredoes the double clustering come from? We used the color greenfor data after 50841 MJD, located in the minimum. A harmonicmight theoretically give rise to two possible clusters, but in thiscase green and black points had to be present in both clusters. Onthe 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 50841 MJDcluster separately with a shift of 0.5 in phase (or about 13 d).We have also folded the simulated data of Eq. 2 with P average .In this case the dependency on P and P is very simple, just twosine functions. Nevertheless, the same kind of double clusteringshown by the GBI data occurs for the simulated data of Eq. 2as one can see in the box of Fig. 4 c. The fact that a simple sumof sine functions in P and P produces the same jump, as theGBI data, if folded with P average , implies either that this simplemathematical form is true in LS I + ◦ P and P .In mathematical terms, as shown in the Appendix, the beat ofthe two sine functions f ( P ) and f ( P ) has a phase reset whenthe delay of f ( P ) with respect to f ( P ) becomes larger than P /
2. Until that point f ( P ) preceeds f ( P ) (and P average ) giv-ing rise to a positive timing residual. After that point f ( P ) (and P average ) preceeds f ( P ). This produces the jump from a largetiming residual to a nearly equally large, but now negative, tim-ing residual observed in the sawtooth function (see Fig. 3 c andd and the derived saw tooth function in Fig. 3 b, in impressingagreement with the observed one of Fig. 3 a). In physical terms the jump is illustrated in Fig. 4 d. The radiomaximum of the long-term modulation results when the ejec-tion, periodical at P = P , occurs at the smallest angle withrespect to the line of sight and the Doppler boosting is largest(Kaufman Bernad´o et al. 2002). Because of the precession P the angle of the ejection changes, and the radio minimum corre-sponds to an ejection occurring at the largest angle with respectto the line of sight. At this point the ejection has travelled half aprecession cone and turns onto the other half of the precessioncone, i.e., from I to II in Fig. 4 d. This causes the phase jumpin the data folded with P average and the jump in the saw toothfunction of Fig. 3 a.Finally, in Fig. 4 d, we present the data folded with P average ,where we fold the data before 50841MJD (black points) usingthe usual t = . t + .
25 d, correct-ing for the jump at the minimum. The better folding of the largeflares with P average with respect to those with the two real peri-odicities P and P is the observational evidence for the result ofEq. 6, i.e., the periodicity of the radio outburst is P average .
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 + ◦
303 system.Zamanov et al. (1999) analyzed H α spectra of LS I + ◦
303 anddetermined the same 26.5 d radio period, and in addition deter-mined that the H α emission line equivalent width (EW) variesover 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 / Rratio.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
5. Massi and F. Jaron: Orbital and precessional periodicities in LS I + ◦ -8-6-4-2 0 2 4 6 8 Timing residuals (d)
Time (MJD)a -8-6-4-2 0 2 4 6 8
0 500 1000 1500
Timing residuals (d)
Time (50000 MJD +)b -0.8-0.4 0 0.4 0.8
720 760 800 840 880 920 960 A m p li t ude Timec-0.4-0.2 0 0.2 0.4
720 760 800 840 880 920 960 A m p li t ude Timed
Fig. 3. a: Sawtooth function based on timing residuals observedby Gregory et al. (1999). b: Sawtooth function based on the lossof synchronization of the sine functions f ( P ) and f ( P ) ofFig. 3 c. c: Sine waves with periods P (black line) and P (redline) starting with synchronized peaks at t =
0. d: Sum of thesine waves. The first five bars are at a regular interval equal to P average , then the next bar follows at a distance of only 13.25 d( P / P average as before. The two central bars are symmetric withrespect to a peak of f ( P ) (Fig. 3 c) that is nearly equidistantfrom the preceeding peak and the delayed peak of f ( P ) (seeAppendix).strength of the blue to the red peak. B / R variability cycles arecommon in Be stars forming canonical Be / X-ray binaries con-taining accreting X-ray pulsars. As stated above, Zamanov andcollaborators tried to find evidence in LS I + ◦
303 for long-term periodicities in the B / R ratio, but they conclude that “un-fortunately, the results of this search turned out to be negative”(Zamanov et al. 1999).V / R variability is explained in terms of a nonaxisymmetricalequatorial disk in which a one-armed perturbation (a zone in thedisk with higher density) propagates (Reig 2011). As a possi-ble explanation for the long-term modulation, Gregory & Neish(2002) indicated that it may stem from periodic ejections of ashell of gas in the equatorial disk of the Be star. Gregory & Neish(2002) commenting the fact that there are no periodic variationsin the H α V / R ratio for LS I + ◦ + ◦ α emission line from the disk of the Be star does not seemto be related to structural disk variations (no B / R variations),whereas it is clearly related to a periodical change in the numberof emitted photons (the EW(H α )). A likely candidate as agent ofthese variations is the relativistic precessing jet that could wellbe able to produce EW(H α ) variations of the same timescales asthose we observe in its radio emission.
0 0.2 0.4 0.6 0.8 1
Flux density (mJy) Φ orbit(P ) a 200
0 0.2 0.4 0.6 0.8 1
Flux density (mJy) Φ precession(P ) b 200
0 0.2 0.4 0.6 0.8 1
Flux density (mJy) Φ average c 200
0 0.2 0.4 0.6 0.8 1
Flux density (mJy) Φ average c 200
0 0.2 0.4 0.6 0.8 1
Flux density (mJy) Φ average + offset d Fig. 4. a: Radio light curves (data averaged over 1 d) vs Φ orbit , with Φ orbit related to ( t − t ) P , with P = t = JD2443366.775 (Gregory 2002), before (black) and after(green) (2400000.5 + Φ precession , related to P . c: Radio light curves vs Φ average for P average = .
70 d. The small window shows thesimulated data (Eq. 2) identically folded. d: Radio light curvevs Φ average with t changed by ∆ t ( ∆ t = + ∆ t = .
25 d after it, see Appendix). e: Sketch ofthe precessing jet in LS I + ◦
303 (out of scale).
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 + ◦
303 results in two frequencies ν = . − ( P = . ± .
07 d) and ν = . − ( P = . ± .
07 d). The aim of our research, to obtain a better determi-nation of the precessional period, indicated by the astrometryto be of 27–28 d, has therefore been reached. The period ex-ists and it is sligthly above the orbital period.2. An additional and totally unexpected result of our researchis that these two periodicities give rise to a P average = ν + ν = . + . = . ± .
05 d, modulated with P beat =
6. Massi and F. Jaron: Orbital and precessional periodicities in LS I + ◦ ν − ν = . − . = ±
393 d. In other words, thelong-term periodicity is equal to the beat of P and P , and P beat modulates a new period, P average .3. We have shown that the sawtooth function derived in the past(Gregory et al. 1999) by comparing observed and predicted(for P = P ) outburst time compares P to P average . Our resulttherefore confirms the controversial value of 26 . ± . P average is only an apparent periodicity, a result of the beatof P and P . When the GBI data are folded with P average the data, instead of clustering at one specific phase, clusterat two phases separated by about 0.5 (13.25 d, i.e., P /
2) de-pending on whether the data are before or after the minimumof the long-term modulation. We have shown that the beat of P and P reproduces the same double clustering as the GBIdata. We find that the time in the minimum of the long-termmodulation when the phase jump occurs is the time whenthe relative delay between the two functions f ( P ) and f ( P )have reached the maximum delay of P /
2. In a physical sce-nario this corresponds, as discussed below, to the point wherethe ejection has travelled half a precession cone and turns onthe other half of the precession cone.5. P average and P beat are related to each other. They are both pro-duced by the beat between P and P , the two real periodic-ities. P beat , the long-term periodicity, is like P average , just anapparent periodicity.6. P and the long-term modulation alone, cannot reproducethe observed periodogram (i.e., P ). It would be possible toreproduce the observed periodogram only by adding a saw-tooth function. In this scenario, the sawtooth function mustbe produced by a physical process continuously changing P to P average until a shift in orbital phase of + / + / / P and P that we found in the periodogramof the GBI data, one naturally explains the two apparent pe-riods P beat and P average , and the phase jump in P average .7. Implications of our results for variations at other wave-lengths indicate the precessing relativistic jet as the agentresponsible for the observed H α variations, equal to the vari-ations that we observe in radio emission of the jet.In conclusion, LS I + ◦
303 seems to be one more case inastronomy of a “beat”, i.e., a phenomenon occurring when twophysical processes create stable variations of nearly equal fre-quencies. The very small di ff erence in frequency creates a long-term variation of period 1 / ( ν − ν ). The first astronomical casewas that of a class of Cepheids, afterwards called beat Cepheids(Oosterho ff + ◦
303 are the precession P of anejection having orbital occurrence P . Following the results ofthe radio spectral index analysis by Massi & Kaufman Bernad´o(2009), the steady ejection of relativistic electrons associated tothe compact object in LS I + ◦
303 increases until it termi-nates in a transient jet twice along the orbit. As quoted in theintroduction, this agrees well with the Bondi & Hoyle (1944)accretion in an eccentric orbit that predicts two events in thesystem LS I + ◦ ff er strong inverse Compton losses because they are exposed to stellar pho-tons, and only at the second accretion / ejection peak, occuringrather displaced from periastron, the relativistic electrons surviveinverse Compton losses and produce the large observed radiooutburst with period P = .
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 ra-dio maximum of the long-term modulation results when the ejec-tion occurs at the smallest angle with respect to the line of sightand the Doppler boosting is largest (Kaufman Bernad´o et al.2002). The radio minimum results when the ejection occurs atlarge angles with respect to the line of sight. The point wherethe ejection has travelled half a precession cone and turns on theother half of the precession cone gives rise to the jump in phasefor the apparent P average .This research has brought up many questions in need of fur-ther investigation. Future work needs to be done to establish theradio vs H α relationship, the possible physical processes un-der the observed precession, and the Doppler boosting e ff ects.Future observations are needed to estimate the angle of the pre-cession cone; an estimate of this angle could result from compar-ing VLBA observations at the same orbital phase but interlapsed(1667 / ff erentposition angle of the radio structures one could infer the apertureof the cone. Of course the two images interlapsed by 31 cyclesshould 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), butat a sligtly di ff erente phase ∆Φ = . ff erences in the images. In Albert et al.(2008) it is discussed that two images taken 10 cycles apart andat the same orbital phase Φ = .
62 are higly similar. Of course10 cycles are still far from the requested 31 cycles. However,the slightly di ff erent position angle between the two images inFig. 1 of those authors indicates that this comparison could be agood investigative tool to estimate the precession angle in futureobservations. Acknowledgements.
We would like to thank the anonymous referee for help-ful comments. We thank Lisa Zimmermann and J¨urgen Neidh¨ofer for read-ing the manuscript and for the several interesting discussions. The Green BankInterferometer is a facility of the National Science Foundation operated by theNRAO in support of NASA High Energy Astrophysics programs.
References
Albert, J., Aliu, E., Anderhub, H., et al. 2008, ApJ, 684, 1351Bondi, H., & Hoyle, F. 1944, MNRAS, 104, 273Bosch-Ramon, V., Paredes, J. M., Romero, G. E., & Rib´o, M. 2006, A&A, 459,L25Cornwell, T. & Fomalont, E. B. 1999, Synthesis Imaging in Radio AstronomyII, 180Dhawan, V., Mioduszewski, A., & Rupen, M. 2006, Proceedings of the VIMicroquasar Workshop, p. 52.1Fender, R. P., Bell Burnell, S. J., & Waltman, E. B. 1997, Vistas in Astronomy,41, 3Fender, R. P., Belloni, T. M., & Gallo, E. 2004, MNRAS, 355, 1105Fragile, P. C., Blaes, O. M., Anninos, P., & Salmonson, J. D. 2007, ApJ, 668,417Gregory, P. C. 1999, ApJ, 520, 361Gregory, P. C. 2002, ApJ, 575, 427Gregory, P. C., & Neish, C. 2002, ApJ, 580, 1133Gregory, P. C., Peracaula, M., & Taylor, A. R. 1999, ApJ, 520, 376Grundstrom, E. D., Caballero-Nieves, S. M., Gies, D. R., et al. 2007, ApJ, 656,437Han, X., & Hjellming, R. M. 1992, ApJ, 400, 304Hannikainen, D. C., Wu, K., Stevens, J. A., Vilhu, O., Rodriguez, J.,Hjalmarsdotter, L., & Hunstead, R. W. 2006, Chinese Journal of Astronomyand Astrophysics Supplement, 6, 010000
7. Massi and F. Jaron: Orbital and precessional periodicities in LS I + ◦ Kaufman Bernad´o, M. M., Romero, G. E., & Mirabel, I. F. 2002, A&A, 385, L10Kovalev, Y. Y., Lobanov, A. P., Pushkarev, A. B., & Zensus, J. A. 2008, A&A,483, 759Linnell Nemec, A. F., & Nemec, J. M. 1985, AJ, 90, 2317Lomb, N. R. 1976, Ap&SS, 39, 447Marti, J., & Paredes, J. M. 1995, A&A, 298, 151Mart´ı-Vidal, I. & J. M. Marcaide, J. M. 2008, A&A, 480, 289Massi, M., & Aaron, S. 1999, A&AS, 136, 211Massi, M., Rib´o, M, Paredes, J. M., et al. 2004, A&A, 414, L1Massi, M. 2005, arXiv:astro-ph / ff , P. T. 1957, Bulletin of the Astr. Inst. of the Netherlands, 13, 320Peracaula, M., Mart´ı, J., & Paredes, J. M. 1997, A&A, 328, 283Pooley, G. G., & Fender, R. P. 1997, MNRAS, 292, 925Porter, J. M., & Rivinius, T. 2003, PASP, 115, 1153Ray, P. S., Foster, R. S., Waltman, E. B., et. al, 1997, ApJ, 491, 381Reig, P. 2011, Ap&SS, 332, 1Rodriguez, L. F., & Mirabel, I. F. 1997, ApJ, 474, L123Romero, G. E., Okazaki, A. T., Orellana, M., & Owocki, S. P. 2007, A&A, 474,15Scargle, J. D. 1982, ApJ, 263, 835Stellingwerf, R.F. 1978, ApJ, 224, 953 2317Stewart, I. M., Fenech, D. M., & Muxlow, T. W. B. 2011, A&A, 535, A81Taylor, A. R., & Gregory, P. C. 1982, ApJ, 255, 210Taylor, A. R., & Gregory, P. C. 1984, ApJ, 283, 273Urry, C. M., & Padovani, P. 1995, PASP, 107, 803Valtaoja, E., Terasranta, H., Urpo, S., Nesterov, N. S., Lainela, M., & Valtonen,M. 1992, A&A, 254, 71Zamanov, R. K., Mart´ı, J., Paredes, J. M., et al. 1999, A&A, 351, 543Taylor, A. R., Kenny, H. T., Spencer, R. E., & Tzioumis, A. 1992, ApJ, 395, 268 Appendix A: Beat and sawtooth function
Let us assume two sine functions f and f with period P = .
49 d and P = .
92 d. Let us start during the radio maximum, i.e., with zero timing residualbetween the peaks of f and f . In Fig. 4 e a zero timing residual corresponds toan ejection with the smallest angle with respect to the line of sight, that is withthe strongest Doppler boosting. Since P is shorter than P , each cycle f ( P )has a delay of P − P = .
43 d. The period formed by the beat P average hasa timing residual τ between its peak and that of f ( P ) half of that delay, i.e.,0 . / = .
21 d. This agrees well with τ = . × P = .
21 d of Eq. 4. f ( P average ) is always between f and f and its peaks have a regular distanceof P average = .
70 d (Fig. 3 d). However, at the minimum the accumulated de-lay of f ( P ) with respect to f ( P ) is almost P /
2. That is, as shown in Fig. 3 c,the peak of f ( P ) is nearly equidistant from the two peaks of f ( P ), the peakabout 13 d before and the peak about 13 d after it. The beat with the first peakgives rise to a peak of f ( P average ) at about 827.8 d, whereas the beat with thesecond peak gives rise to a peak of f ( P average ) at 841.05 d. The di ff erence be-tween these two consecutive peaks of f ( P average ) is 13.25 d. This corresponds toabout 0.5 in phase, when the data are plotted in phase for P average (Fig. 4 c), andalso corresponds to the jump of about 13 d in the sawtooth function observed byGregory et al. (1999). After that point, f ( P ) is preceeding f ( P ). In the physicalscenario of Fig. 4 e this jump, or reset of the phase of P average , corresponds to thepoint where the ejection has travelled half a precession cone and turns onto theother half of the precession cone. Figure 3 b shows τ vs time, resulting from thissimple analysis based on the sine functions of Fig. 3; the resulting slope of thefunction is indeed 0.008, as in Fig. 3 a.vs time, resulting from thissimple analysis based on the sine functions of Fig. 3; the resulting slope of thefunction is indeed 0.008, as in Fig. 3 a.