Issue 
A&A
Volume 643, November 2020



Article Number  A170  
Number of page(s)  15  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202038745  
Published online  20 November 2020 
Orbital variability of the optical linear polarization of the γray binary LS I +61° 303 and new constraints on the orbital parameters
^{1}
Department of Physics and Astronomy, University of Turku, 20014 Turku, Finland
email: vakrau@utu.fi
^{2}
Space Research Institute of the Russian Academy of Sciences, Profsoyuznaya Str. 84/32, Moscow 117997, Russia
^{3}
Department of Astrophysics, St. Petersburg State University, Universitetskiy Pr. 28, Peterhof, 198504 St. Petersburg, Russia
^{4}
School of Physical Sciences and CfAR, Dublin City University, Dublin 9, Ireland
^{5}
Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, Dublin 2, Ireland
^{6}
Institut für Astronomie und Astrophysik Tübingen, Universität Tübingen, Sand 1, 72076 Tübingen, Germany
^{7}
Graduate School of Sciences, Tohoku University, Aobaku, 9808578 Sendai, Japan
^{8}
LeibnizInstitut für Sonnenphysik, Schöneckstr. 6, 79104 Freiburg, Germany
^{9}
Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 968221897, USA
^{10}
Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden
Received:
25
June
2020
Accepted:
14
September
2020
We studied the variability of the linear polarization and brightness of the γray binary LS I +61° 303. Highprecision BVR photopolarimetric observations were carried out with the Dipol2 polarimeter on the 2.2 m remotely controlled UH88 telescope at Mauna Kea Observatory and the 60 cm Tohoku telescope at Haleakala bservatory (Hawaii) over 140 nights in 2016−2019. We also determined the degree and angle of the interstellar polarization toward LS I +61° 303 using two out of four nearby field stars that have Gaia’s parallaxes. After subtracting the interstellar polarization, we determined the position angle of the intrinsic polarization θ ≃ 11°, which can either be associated with the projection of the Be star’s decretion disk axis on the plane of sky, or can differ from it by 90°. Using the LombScargle method, we performed timing analyses and period searches of our polarimetric and photometric data. We found statistically significant periodic variability of the normalized Stokes parameters q and u in all passbands. The most significant period of variability, P_{Pol} = 13.244 ± 0.012 d, is equal to one half of the orbital period P_{orb} = 26.496 d. The fits of the polarization variability curves with Fourier series show a dominant contribution from the second harmonic which is typical for binary systems with circular orbits and nearly symmetric distribution of light scattering material with respect to the orbital plane. The continuous change of polarization with the orbital phase implies coplanarity of the orbit of the compact object and the Be star’s decretion disk. Using a model of Thomson scattering by a cloud that orbits the Be star, we obtained constraints on the orbital parameters, including a small eccentricity e < 0.2 and periastron phase of ϕ_{p} ≈ 0.6, which coincides with the peaks in the radio, Xray, and TeV emission. These constraints are independent of the assumption about the orientation of the decretion disk plane on the sky. We also extensively discuss the apparent inconsistency with the previous measurements of the orbital parameters from radial velocities. By folding the photometry data acquired during a threeyear time span with the orbital period, we found a linear phase shift of the moments of the brightness maximum, confirming the possible existence of superorbital variability.
Key words: binaries: general / gamma rays: stars / polarization / stars: emissionline, Be / stars: individual: LS I +61 303
© ESO 2020
1. Introduction
Gammaray binaries constitute a subclass of highmass binary systems with the emission peaking in the GeV band (see recent review by Chernyakova & Malyshev 2020). In fact, LS I +61° 303 is one of the beststudied γray binaries and was observed in the last few decades over the whole range of the electromagnetic spectrum, from the radio to the very highenergy γrays (Taylor & Gregory 1982; Paredes et al. 1994; Zamanov et al. 1999; Harrison et al. 2000; Abdo et al. 2009; Albert et al. 2006). This binary system consists of a B0 Ve star with a circumstellar disk (Casares et al. 2005) and a compact companion star orbiting the primary star on an apparently eccentric (e ≥ 0.5) orbit. The nature of the compact object, a black hole, or a neutron star is still unknown.
Using a Bayesian analysis of 20 years of radio data, Gregory (2002) determined the orbital period P_{1} = 26.4960 ± 0.0028 d. The compact object moving around the Be star and interacting with circumstellar matter produces the orbital variability seen in all parts of the spectrum (Taylor et al. 1992; Mendelson & Mazeh 1989; Leahy 2001; Grundstrom et al. 2007). A LombScargle timing analysis of 37 years of radio data resulted in the detection of the second period, P_{2} = 26.935 ± 0.013 d (Massi & TorricelliCiamponi 2016), which is consistent with a previously determined period of morphological changes in the radio structure, mapped in the VLBI images (Massi et al. 2012). Recently, the VLBA astrometry increased the accuracy of the period to P_{2} = 26.926 ± 0.005 d (Wu et al. 2018). The beat of these two periods d is very close to the period of the longterm superorbital variability P_{sup} ≈ 1700 d observed in the radio (Massi & Jaron 2013), Hα emission (Zamanov et al. 2013), Xrays (Chernyakova et al. 2012; Li et al. 2014), and γrays (Ackermann et al. 2013).
From the analysis of the optical spectropolarimetric data, Nagae et al. (2006, 2009) showed that the polarization position angle (PA) was stable over the twoyear period at θ ≃ 25°. According to these studies, the linear polarization in LS I +61° 303 arises as a result of Thomson scattering in the Be star’s decretion disk, and this disk is most likely coaligned with the orbit of the compact companion (Nagae et al. 2009).
In this paper, we present the results of our BVR highprecision photopolarimetric observational campaign of LS I +61° 303, which allowed us to determine the orbital period directly from the variability of the Stokes parameters of linear polarization for the first time and provide new constarints on the orbital parameters. In addition, we discuss our photometry data obtained for this object in terms of both orbital and superorbital variability.
2. Photopolarimetric observations
2.1. Observed polarization
The observations of LS I +61° 303 were performed with the broadband BVR polarimeter Dipol2 (Piirola et al. 2014) mounted on the remotely controlled 2.2 m UH88 telescope at Mauna Kea Observatory and the 60 cm Tohoku telescope (T60) at Haleakala Observatory, Hawaii. The object was observed during 140 nights from 2016 September 15 to 2019 September 1 (MJD 57646–58728). Every night, 24 to 48 measurements of the normalized Stokes parameters q and u were made simultaneously in the B, V, and Rpassbands. The total integration time with a typical tensecond exposure was 25−50 min. A summary of our observations is given in Table 1. The observational errors of q and u were computed as the standard errors of the weighted mean values and are in the range of 0.01−0.03%. In each of our five observing runs, the value of instrumental polarization was determined from observations of 15−20 nearby bright, unpolarized stars. The value of instrumental polarization is ≤5 × 10^{−5} for all passbands and is measured with the accuracy of a few parts per million (10^{−6}). To determine the PA zero point, the highly polarized standard stars HD 204827 and HD 25443 were observed. Detailed descriptions of the calibration and observation procedures can be found in Kosenkov et al. (2017) and Piirola et al. (2020). Table 2 shows the observed average values of the polarization degree (PD) P_{obs} and the PA θ_{obs}.
Log of polarimetric observations of LS I +61° 303.
We also extracted the fluxes of LS I +61° 303 and the closest field star (star 4, see Fig. 1) in the B, V, and R passbands from our polarimetry images acquired with the T60 telescope, and we used them for the measurements of relative brightness variations of LS I +61° 303.
Fig. 1. Polarization map of LS I +61° 303 (red circle at the origin) and field stars (black circles) in the R band. The length of the bars corresponds to the degree of linear polarization P, and the direction corresponds to the PA θ (measured from the north to the east). 
2.2. Interstellar and intrinsic polarization
In order to obtain the degree and direction of the intrinsic linear polarization, P_{int} and θ_{int}, it is necessary to estimate the parameters of the interstellar (IS) polarization. For this purpose, we observed the polarization of four of the nearest field stars, located at the angular distance of < 7′ from LS I +61° 303 (Fig. 1) and with known parallaxes. Figure 2 shows the dependence of the observed PD P and PA θ on parallax for LS I +61° 303 and the field stars (parallaxes are taken from the Gaia Data Release 2; Prusti et al. 2016; Brown et al. 2018). The directions of the IS polarization for the field stars lie in the 110° −125° interval. Assuming that the field stars are intrinsically unpolarized, we see that the degree of IS polarization decreases with distance (i.e., it grows with parallax) from ≈4% to ≈2% in the closest vicinity of the LS I +61° 303. This unusual behavior may result from the depolarization effect in several dust clouds with nearly orthogonal orientations of the Galactic magnetic field, located along the line of sight toward LS I +61° 303.
Fig. 2. Dependence of observed PD P and PA θ on parallax (Gaia DR2) for LS I +61° 303 (red) and field stars (black) in the B band. The numbering of the field stars is the same as in Fig. 1. The error bars correspond to the 1σ errors. 
For a more detailed study of IS polarization, we analyzed the data from Heiles’ stellar polarization catalog (Heiles 2000) in a 10° ×10° region around LS I +61° 303. The dependence of PD on distance for 232 stars in this area of the sky is shown in Fig. 3. We see that there is a significant scatter of the IS PD for distant stars. Moreover, the PD of LS I +61° 303 is smaller than that of all field stars at similar distances. It means that LS I +61° 303 has a significant intrinsic polarization of which the direction does not match that of the IS polarization. The large scatter in the degree of IS polarization at the distances d > 2 kpc may be linked to the complex structure of the IS medium inside the Heart Nebula (IC 1805) that is located in close proximity to LS I +61° 303. Thus, a careful approach, taking into account both proximity to the line of sight and proximity in distance, must be used for estimating the IS polarization for this object.
Fig. 3. Dependence of linear PD P (Heiles 2000) on distance (evaluated from the inverse of the Gaia DR2 parallaxes) for the stars in the 10° ×10° area around LS I +61° 303 in the R band. The red star indicates our value of the observed polarization of LS I +61° 303, the orange crosses correspond to the four nearby field stars shown in Fig. 1. The vertical bar corresponds to the 1σ error in polarization. 
We chose the average polarization of two nearby field stars (#1 and #2), which are close in distance to LS I +61° 303, as the best estimate for the IS polarization in the direction of the binary. Table 2 shows the estimated values P_{is} and θ_{is} of the IS polarization in all passbands and the average values P_{int} and θ_{int} of the intrinsic polarization for LS I +61° 303, obtained after subtracting IS polarization. The values of IS polarization derived by us are in good agreement with those obtained by Nagae et al. (2006) from the polarization in Hα emission line: P_{is} = 2.20 ± 0.18%, .
Observed PD and PA of LS I +61° 303, interstellar polarization and average intrinsic polarization.
The values of the average intrinsic polarization of LS I +61° 303 in the B, V, and Rbands are the same within the errors. This suggests that Thomson scattering on free electrons in the disk around the Be star is likely the polarization mechanism responsible for the constant component of polarization in LS I +61° 303. This conclusion is in agreement with the results obtained by Nagae et al. (2006, 2009). It is interesting to note that the direction of the average intrinsic polarization in all passbands differs from the value of 25° derived by Nagae et al. (2006). The PA of the intrinsic optical polarization of LS I +61° 303, obtained by us as the average for the BVR bands, is ≃11°. This difference, however, is most likely a result of uncertainty in the determination of the IS polarization and does not imply physical changes of the disk orientation with time. As we mentioned above, we used the average polarization of two (close in distance) field stars as an estimate of the IS polarization, whereas Nagae et al. (2009) used polarization in Hα and the empirical Serkowski law (Serkowski 1973), which is not applicable when the structure of the IS medium is complex (i.e., consisting of multiple clouds of different properties), as is demonstrated by a peculiar distance dependence of the PD for the nearby field stars shown in Fig. 3.
The intrinsic polarization PA allows us to put constraints on the orientation of the Be circumstellar disk on the sky. As was shown by Quirrenbach et al. (1997) for four Be stars, the polarization vector is parallel to the projection of the rotation axis of the circumstellar disk. However, all their stars showed wavelengthdependent PD, which is an indication of the important role of boundfree hydrogen absorption in the envelope of a Be star. On the other hand, our intrinsic PD does not depend on the wavelength, implying that the electron scattering dominates. In this case, if the disk is optically thick, the polarization vector of radiation escaping from the disk may be perpendicular to the PA of the rotational axis (Chandrasekhar & Breen 1947; Sobolev 1949; Wood et al. 1996). In the following, we consider both possibilities.
3. Orbital variability
3.1. Polarization variability
The variability curves for the Stokes parameters q and u in the V band, obtained over three years of observations of LS I +61° 303, are shown in Fig. 4. Our observations have revealed a small but significant variability with the amplitude of 0.2−0.3% in all passbands. To study this variability, we performed a timing analysis of our BVR polarimetric data. For this purpose, we applied the LombScargle method^{1} of spectral analysis (Lomb 1976; Scargle 1982), using the ASTROPY package (Astropy Collaboration 2018). The LombScargle periodograms for the normalized Stokes parameter u in the BVR passbands are shown in Fig. 5. We see that the period of the highest peak in each band is close to one half of the orbital period P_{orb} = 26.4960 ± 0.0028 d, as is expected in binary systems. Despite significant nonperiodic scatter, this peak is clearly present in all passbands and both Stokes q and u. For example, the period in the V band in Stokes u is P_{V} = 13.25 ± 0.06 d with a false alarm probability ∼10^{−5}, which was independently estimated using an analytical approximation of Baluev (2008) and a bootstrap method from the ASTROPY package. The error on the found period was estimated in two different ways. In the first case, we simulated 10^{4} light curves, varying the observed data points around the mean values in accordance with the observational errors. In the second case, we also simulated 10^{4} light curves assuming sinusoidal variations with the fixed period P = 13.25 d and the amplitude, time stamps, and errors as in the observed data. For both cases, we computed the periodograms of the light curves and estimated an error on the found period as the half width of the distribution of the highest peak positions. The value of the error in the first case was from three to ten times larger than in the second case, depending on the passband, due to the presence of the superorbital variability and other possible periodicities in the data. For the final error estimate, we chose the largest of the two. The mean for the three passband periods in Stokes u is P_{pol} = 13.244 ± 0.012 d, which differs by less than 1σ from P_{orb}/2. Thus, for the very first time, the orbital period of LS I +61° 303 was obtained directly from the polarimetric measurements.
Fig. 4. Observed Stokes parameters q (top panel) and u (bottom panel) with 1σ errors for LS I +61° 303 in the V band, measured from 2016 September 15 to 2019 September 1 (MJD 57646–58728). 
Fig. 5. Top panel: LombScargle periodogram for normalized Stokes parameter u in the B band. Bottom panel: LombScargle periodogram zoom for the normalized Stokes parameter u in the B, V, and R bands (solid blue, dashed green, and dotted red curves, respectively). The horizontal dashed line corresponds to false alarm probability (FAP) = 1%. The vertical red line corresponds to one half of the orbital period P_{orb} = 26.496 d. 
Folding the data with the orbital period using the ephemeris given in Aragona et al. (2009), we obtained the phase curves for the PD P_{int} and the PA θ_{int} of the intrinsic polarization (Fig. 6) and for the normalized Stokes parameters q_{int} and u_{int} (Fig. 7). The synchronization of the variable component of linear polarization in LS I +61° 303 with the phase of the orbital period clearly indicates that it arises due to the orbital motion of the compact component. In the binary system with the compact companion and Be star, the extended hot region with enhanced electron density can be formed around or near a neutron star or a black hole in a process of their interaction with the disk material. Thomson scattering of stellar photons on such a structure orbiting Be star can explain the regular smallamplitude variations in the observed polarization of LS I +61° 303. Moreover, a continuous change of polarization implies coplanar orientation of the orbit and the Be star decretion disk. This conclusion directly follows on from the observed polarization variability and does not require any assumptions on the nature of compact companion (see Nagae et al. 2009). There is also a strong nonperiodic component, which arises due to the long and shortterm changes in the distribution and density of the light scattering material in the Be decretion disk.
Fig. 6. Orbital variability of intrinsic PD P_{int} and PA θ_{int} of LS I +61° 303 in V band. The filled blue squares with 1σ errors correspond to the average values of the individual observations (represented by gray circles with error bars) and the standard errors of the mean calculated within the phase bins of width Δϕ = 0.091. The vertical dashed lines correspond to the phases of periastron (ϕ = 0.275) and apastron (ϕ = 0.775) as derived by Aragona et al. (2009). 
Fig. 7. Variability of intrinsic Stokes parameters q_{int} (left column) and u_{int} (right column), in B, V, and R bands (top, middle and bottom panels, respectively) with the orbital phase folded with the period P_{orb} = 26.496 d. The filled squares with 1σ errors correspond to the average values of the individual observations (gray circles with error bars) and the standard errors of the mean calculated within phase bins of width Δϕ = 0.091. The vertical dashed lines correspond to the phases of periastron (ϕ = 0.275) and apastron (ϕ = 0.775) as derived by Aragona et al. (2009). The red dashed lines show the best fit with Fourier series given by Eq. (1). The black solid and blue dotted lines correspond to the best fit model of a scattering cloud on an eccentric orbit from Sect. 3.3, for an Ω constrained close to θ_{int} and with a free Ω, respectively. 
The dominant second harmonic in the variability of the Stokes parameters is typical for a binary system with symmetric density distribution of the scattering matter about the orbital plane and a circular orbit (Brown et al. 1978). The polarization arising from light scattering on gaseous structures like clouds, streams, and so on, peaks at orbital phases where the scattering angle is close to 90°. For a corotating envelope and a circular binary orbit, this gives two prominent symmetric peaks in polarization separated by the phase interval ≃0.5. In the case of a (weak) asymmetry of the distribution of the light scattering material about the orbital plane, the first harmonic appears. This is apparently observed in LS I +61° 303. As we see in Figs. 6 and 7, there are two nearly symmetric peaks near the phases 0.5 and 1.0. They are most pronounced in the variations of the Stokes u and the PA θ.
3.2. Modeling polarization with circular orbit
According to the common approach (Brown et al. 1978), the phase curves of the Stokes parameters in binary systems with a circular orbit and corotating light scattering envelope can be represented through a Fourier series of the orbital longitude λ = 2πϕ (where ϕ is a phase of the orbital period):
We fit the intrinsic Stokes parameters q_{int} and u_{int} with these functions, using q_{0}–q_{4} and u_{0}–u_{4} as free parameters. The best fit parameters together with the χ^{2} values and the degrees of freedom (d.o.f.) are given in Table 3. The best fit curves are shown in Fig. 7 (red dashed lines). We see that variabilities of Stokes parameters q and u are dominated by the second harmonic of the orbital period. This is expected, because of the position of the highest peak in the LombScargle periodogram close to the frequency corresponding to one half of the orbital period. The parameters
giving the ratio of the amplitudes of the second to the first harmonic, attain the values A_{q} = 0.9, 1.0, 3.3 and A_{u} = 2.5, 5.5, 5.4 for the B, V, and R bands, respectively.
It is possible to derive the inclination i of the binary orbit and the position angle Ω of the projection of the orbital axis (see Eqs. (A.17) and (A.18)) from the best fit Fourier coefficients in Eq. (1) to the observed (or intrinsic) Stokes parameters q and u (Drissen et al. 1986):
where
However, there is a bias in the polarimetrically derived orbit inclination because of the unavoidable noise (finite accuracy) in the polarization data (Aspin et al. 1981; Simmons et al. 1982; Wolinski & Dolan 1994). The inclination of orbit i derived from the best fit Fourier coefficients is always biased toward a higher value. For noisy data, the inclination approaches 90° with wide confidence intervals extending to very low values (see Wolinski & Dolan 1994). Similar bias can also be induced by stochastic noise, arising from an intrinsic nonperiodic component of the polarization variability (Manset & Bastien 2000). As the errors on q and u increase, a straight line becomes an acceptable fit to the (q, u) light curves. Because a straightline fit to the (q, u) data yields a 90° inclination for the system (Brown et al. 1978), the derived value of i is biased toward 90° when a high noise level is present in the data.
Due to the strong nonperiodic component in the polarization variability of LS I +61° 303, the inclination of the orbit derived from the Fourier coefficients for all passbands is close to 90°, with 2σ confidence intervals extending down to values of i ≤ 20°. Thus, assuming a (nearly) circular orbit in LS I +61° 303, we can only say that our polarization data are consistent with orbit inclination i in the range ≈20° −90°. The position angle Ω, derived from Eq. (4), is about 30° for all passbands. As was shown by Wolinski & Dolan (1994), there is no bias in the value of Ω derived from polarimetry. This is due to the fact that the value of Ω is determined by the positions of the polarization maxima and minima, which are not affected by the noise in the same way as the amplitude and shape of the polarization variability curve.
3.3. Modeling polarization with an eccentric orbit
The latest value of eccentricity e derived for LS I +61° 303 from the radial velocity variations in the optical spectral lines is 0.54 ± 0.03 (Aragona et al. 2009). For such an eccentric orbit, one should expect fast changes in polarization occurring near the periastron passage (Brown et al. 1982; Simmons & Boyle 1984). According to the orbital geometry shown in Fig. 3 of Aragona et al. (2009), we should expect one sharp peak just before, and another one soon after, the periastron passage, which occurs at phase 0.275. These two peaks must be followed by a gradual change of polarization during the remaining part of the orbital motion toward apastron. Such a behavior of intrinsic polarization has been detected in the interacting binary HD 187399, which has an orbit eccentricity of e ≃ 0.4 (Berdyugin & Tarasov 1998), but this is not observed in LS I +61° 303.
Using the approach proposed by Simmons & Boyle (1984), we tried to model the observed shape of the polarization variability curve in LS I +61° 303. We considered a simple scenario of a small cloud of free electrons orbiting the central illuminating source and scattering its radiation (see Appendix A for the model). Under this assumption, the modeled Stokes parameters depend on the eccentricity e, the orbit inclination i, as well as on the longitude of the periastron λ_{p}, and typical scattering fraction f_{0} (see Eq. (A.16)). Polarization, produced by such a system, depends on time through the orbital longitude λ, which can be connected to the observed orbital phase ϕ through an eccentric anomaly (Eq. (A.19)) and a Kepler equation (Eq. (A.20)). The orbital phase is then corrected for the phase of the periastron ϕ_{p} (corresponding to λ = λ_{p}) using Eq. (A.21). The model Stokes pseudovector should also be rotated to account for the position angle Ω of the projection of the orbital axis (see Eqs. (A.17) and (A.18)) relative to the north. The constant levels of intrinsic polarization corresponding to the emission from the decretion disk is represented by the pseudovectors (q_{0}, u_{0}).
We employed Bayesian inference to fit this model to the intrinsic polarization light curves of LS I +61° 303 in three filters simultaneously, allowing different levels of constant polarization (q_{0}, u_{0}) and a typical scattering fraction f_{0} in each band. We applied a Hamiltonian Monte Carlo algorithm (Duane et al. 1987) implemented in the GRETA package (Golding 2019) for R (R Core Team 2019), which utilizes the TENSORFLOW backend (see Abadi et al. 2015). First, we took the eccentricity e = 0.54 and the argument of the periastron ω = λ_{p} ± 90° = 40° from Aragona et al. (2009), but other parameters were allowed to vary in the broadest possible intervals. We were not able to fit the observed variability of the Stokes parameters q and u for the whole range of reasonable inclinations (i.e., i = 20° −90°) resulting in χ^{2}/d.o.f. > 180/54. Then, we fixed only eccentricity e ≃ 0.5 and tried to fit the curves by varying other orbital parameters. We could only fit variations of one Stokes parameter reasonably well, with the resulting χ^{2}/d.o.f. = 177/53. Thus, we were not able to find a set of orbital parameters providing a sufficiently good fit for both Stokes q and u for e ≥ 0.5. This can be simply explained by the influence of eccentricity on the shape of the polarization curve. We can see that in our data (Fig. 7) the distance between the neighboring maxima is Δϕ ≈ 0.4, while in the case of the eccentricity e ≈ 0.5, it should be Δϕ ≈ 0.15.
Finally, we allowed all parameters, including eccentricity e, to vary. At the first step, we assumed that the projection of the decretion disk axis on the sky coincides with the PA of the average intrinsic polarization. If the orbit of the compact object and the decretion disk are nearly coplanar (see Sect. 3.1), we expect that Ω ≈ θ_{int}. We thus limited the prior distribution of Ω to follow a Gaussian with the peak at 11° (θ_{int} from Table 2) with an arbitrarily chosen standard deviation of 10° and truncation interval of [ − 10° ;40° ]. The best fit parameters together with the values of χ^{2} and the d.o.f. are given in the upper part of Table 4, and a comparison of the model with the data is shown in Fig. 7. The posterior distributions are shown in Fig. 8. A corresponding orbital geometry is demonstrated in the left panel of Fig. 9.
Fig. 8. Posterior distributions for parameters of model described in Sect. 3.3 and Appendix A for Ω = θ_{int} ± 10°. Diagonal panels: distributions of model parameters. The blue solid, green dashed, and red dotdashed lines correspond to perpassband B, V, and R distributions, respectively. Lowertriangle panels: joint posterior distributions of two parameters. The green, blue, and violet contours correspond to 0.68, 0.95, and 0.997 probability levels of parameters i, Ω, λ_{p}, and ϕ_{p}. The shades of blue, green, and red correspond to the same probability levels for B, V, and R filters, respectively. 
Fig. 9. Relative orbit of a compact object in LS I +61° 303 around Be star (yellow circle at the origin), which lies at the ellipse focus. The orbital parameters are taken from Table 4. The red dashed line is the major axis of the orbit. The black dots on the ellipse are spaced by the Δϕ = 0.1. Because of a 180° degeneracy in λ_{p}, the orbit can be also rotated by 180°. Left panel: orbit for constrained orientation relative to the decretion disk plane with Ω = θ_{int} ± 10°. Right panel: orbit assuming free Ω. 
We see that our eccentric orbit fit gives an upper limit on e of 0.15, with the best fit value e ≃ 0.06. The high value of inclination i obtained from the eccentric orbit fit is most likely a result of a similar bias as for the circular orbit model. The best fit value of Ω ≈ 28° coincides with the value of 30° obtained from the fits with the circular orbit model in Sect. 3.2. The best fit phase of the periastron ϕ_{p} ≈ 0.62 differs significantly from the usually assumed value of 0.275. We note that some of the parameters are degenerate: for example, λ_{p} and Ω can be substituted by λ_{p} ± 180° and Ω ± 180° without affecting the Stokes parameters. Furthermore, at low eccentricities (such as the best fit values that we obtained), λ_{p} and ϕ_{p} are correlated (see Eqs. (A.19)–(A.21)). This effect can be observed in the joint posterior distribution of λ_{p} and ϕ_{p} shown in Fig. 8.
As it follows from the formulae given in the appendix of the paper by Simmons & Boyle (1984)^{2} a large orbit eccentricity should lead to the appearance of a noticeable third harmonic of the orbital period in the variations of q and u. However, in the observed polarization variability of LS I +61° 303, the contribution from this harmonic is insignificant. A good fit can only be obtained if we assume a nearly circular orbit (producing a strong second harmonic). Thus, the behavior of the variable component of linear polarization in LS I +61° 303 is in an apparent contradiction with a high (e ≥ 0.5) value of the eccentricity of the binary orbit.
We note that the χ^{2} value of the fit with constrained Ω seems too large (67 for 52 d.o.f.). Therefore, at the next step, we relaxed the constraints on Ω, choosing a very wide prior interval. The best fit parameters together with the values of χ^{2} and the d.o.f. are shown in the lower part of Table 4. The posterior distribution is shown in Fig. 10. The comparison of the model to the data is shown in Fig. 7 with the blue dotted line. A sketch of the possible orbit is pictured in the right panel of Fig. 9. We note that the χ^{2} value of the fit is much smaller than in the case of constrained Ω (50 versus 67) and the bestfit Ω ≈ 120°, which differs by 90° from the value obtained in the first case. Such a value for Ω can be interpreted in two different ways. This might indicate that the orbit of the compact object is nearly perpendicular to the disk plane, which is not likely (see Sect. 3.1). A more probable interpretation is that the projection of the decretion disk axis makes a 90° angle to the average intrinsic PA, and the orbit of the compact object is nearly coplanar to the disk (i.e., Ω ≈ θ_{int} + 90°). In this geometry, the polarization vector of radiation scattered in a cloud is predominantly perpendicular to the PA of the average polarization associated with the disk. This explains why, at the phases where the scattering angle is about 90° and maximum polarization of the scattered radiation is expected (i.e., at ϕ ≈ 0.2 and 0.7), the minima are observed in both q and u. Despite a very different geometry, the best fit eccentricity is identical to the case of Ω ≈ θ_{int}. Also the best fit phase of the periastron ϕ_{p} is still very close to 0.6. Thus, we conclude that the obtained values for the eccentricity and the periastron phase are robust and do not depend on our assumptions on the disk and compact object orbit orientations. We discuss the implications of the obtained orbital parameters in Sect. 5.
4. Optical photometry
In addition to the optical polarimetry, we also studied the brightness variations of LS I +61° 303 in the B, V, and R passbands by measuring the ratio of the fluxes from the binary and the nearest field star #4. For the images taken with the UH88 telescope, the small size of the field of view places star #4 too close to the edge of the image field, preventing us from using these images for reliable relative photometry. Thus, only the images taken with the T60 telescope have been used. An example of the light curve in the V band, obtained from three years of observations, is shown in Fig. 11.
Fig. 11. Variation of relative amplitude of LS I +61° 303 around the mean in the V passband as a function of superorbital phase. The shaded blue bands correspond to the eleven observing seasons (S1–S11). 
As in the case of the Stokes parameters, we used a LombScargle timing analysis to find the period of the variability of the brightness in LS I +61° 303. The frequency of the highest peak on the LombScargle periodograms (Fig. 12) is close, but not exactly equal to the orbital period P_{orb} = 26.4960 ± 0.0028 d. For the BVR bands, the determined periods are P_{B} = 26.56 ± 0.05, P_{V} = 26.58 ± 0.04, and P_{R} = 26.60 ± 0.03 d, with the false alarm probabilities ∼10^{−5}, 10^{−10}, and 10^{−14}, respectively. The errors of the periods were estimated via the randomization of the data within 6σ interval and recalculation of the periodogram (number of recalculations N = 1000).
Fig. 12. Top panel: LombScargle periodogram for brightness variability of LS I +61° 303 in the R band. Bottom panel: LombScargle periodogram zoom for brightness variability in the BVR bands (dotted blue, dashed green, and solid red lines, respectively). The horizontal dotted line corresponds to the FAP = 1%. The vertical red dotted line marks the orbital period P_{orb} = 26.496 d. 
The light curves of LS I +61° 303, folded with the orbital period P_{orb} are shown in Fig. 13. The shape of these curves (sinusoidallike wave with an amplitude of about 0.03 mag, minimum near the periastron and maximum in the apastron) is in a good agreement with the previous results obtained by Mendelson & Mazeh (1989), Zaitseva & Borisov (2003), and ParedesFortuny et al. (2015). Like in the case of polarization variability, there is a noticeable scatter around the light curves in all passbands, which indicates the presence of a nonperiodic component.
Fig. 13. Variability of LS I +61° 303 brightness in BVR bands (top, middle, and bottom panels, respectively) with the orbital period. The filled squares with 1σ errors correspond to the average values of the individual observations (gray crosses) and the standard errors of the mean calculated within the phase bin of width Δϕ = 0.091. 
In order to find a possible superorbital phase shift of the maximum of brightness, which is discussed for example in ParedesFortuny et al. (2015), we divided our data into eleven subsets (see Table 5 and Fig. 11). We then folded them with the orbital period P_{orb} and fit with the function
where ϕ_{0} (phase of the peak), A (amplitude), and m_{0} (vertical offset) are free parameters. The data and the best fit curves are shown in Fig. 14. The values of the best fit parameters are given in Table 5. The errors on the parameters were calculated as a square root of the diagonal elements of the covariance matrix of the fit. The measurement errors are smaller than the intrinsic scatter in the photometric data. For a more accurate estimation of errors of the fit parameters, we used the standard deviation of the data instead of the measurement errors for the fitting (in that case, the reduced χ^{2} ≈ 1).
Fig. 14. Variability of LS I +61° 303 optical brightness in Vband for different parts of the data (seasons S1–S11 from top to bottom). The black solid lines correspond to the best fit of the data with function (6). The vertical dashed lines give the phases of the best fit maxima ϕ_{0}, and the corresponding ±1σ, ±2σ, and ±3σ confidence intervals are shown with varying shades of blue. The vertical scale [0.1, −0.1] is the same in all panels. 
We see from Fig. 14 that there is a significant phase shift by Δϕ_{0} ≈ 0.3 seen between observing seasons S1 and S11. The quality of fits for the first two seasons, S1 and S2, are quite low, but the linear shift of the phase of the maximum of brightness is apparent even without them. Thus, our new photometry data are in agreement with the results obtained by ParedesFortuny et al. (2015). The relatively short time baseline (three years) does not allow us to conclude whether this shift is periodic (superorbital) or linear. We have also studied a dependence of the amplitude of brightness variability over time (see bottom panel of Fig. 15). As is seen, there is no apparent trend here, unlike the one seen in the phase of the maximum. We must emphasize that, due to the presence of the nonperiodic brightness fluctuations, studying possible periodic (superorbital) events is complicated. These fluctuations introduce significant noise, which may completely diminish smallamplitude gradual changes occurring over long time scales.
Fig. 15. Dependence of brightness maximum orbital phase ϕ_{0} (top panel) and the amplitude A (bottom panel) of the sinusoidal fits on the superorbital phase. The solid blue lines correspond to the linear fit of the data, while the ±1σ and ±3σ confidence intervals are shown in dark and light blue. The red squares show the parameters from Table 3 of ParedesFortuny et al. (2015). 
It is worth pointing out that the observed light curve does not necessarily imply a high eccentricity of the binary orbit. On the contrary, if the brightness variations result from disruption of the disk by the compact object at the periastron passage (see ParedesFortuny et al. 2015), one would expect light curve asymmetry, that is, a fast dip near the previously assumed periastron phase ϕ_{p} ≃ 0.275 with more gradual growth toward apastron at ϕ_{a} ≃ 0.775. The observed light curves in all passbands, apart from being noisy, do not show such asymmetry. The small and nearly equal amplitude in all passbands and the shape of light curve variability, can be similarly (if not better) explained by the occultation of the Be star or a bright emission area, by a dense gaseous cloud. The brightness variability can be interpreted in terms of viewing geometry, and this explanation does not require an eccentric orbit and periodic disk disruption. We note that the phase of the maximum brightness, ϕ_{0} = 0.78 ± 0.15, differs significantly from the previously assumed periastron phase ϕ_{p} = 0.275, it is, however, much closer to our estimate of the periastron phase of ≈0.6. We discuss this fact in Sect. 5.2.
5. Orbital parameters of LS I +61° 303
5.1. Eccentricity
The estimates of the orbit eccentricity e derived for LS I +61° 303 by different methods, vary across a wide range from 0.3 to 0.8 (see Grundstrom et al. 2007, and references therein). The most reliable and direct way to reconstruct the binary orbit is to measure radial velocity (RV) variations of the stellar atmospheric spectral lines. Several efforts have been made to find orbital parameters for LS I +61° 303 from the RV variations (Hutchings & Crampton 1981; Casares et al. 2005; Grundstrom et al. 2007; Aragona et al. 2009). The latest solutions obtained from the RV curves for the He I and He II lines give the value of eccentricity from ≃0.54 (Aragona et al. 2009) to ≃0.72 (Casares et al. 2005).
The common feature seen in all RV curves for the LS I +61° 303 is a high degree of scatter of individual measurements around the fitting curve, which is particularly pronounced in the orbital phases after the previously assumed periastron phase around 0.3 and the 0.8−0.9 phase range. Another feature is the apparent presence of a “secondary bump” near the phase 0.7−0.8 (Grundstrom et al. 2007; Aragona et al. 2009).
We want to emphasize that such a shape of the RV curve is not exceptional, but rather typical for many Be components in Xray and nonXray binary stars. The formal solution of such a curve often results in substantial eccentricity for the Be star orbit. An explanation of this phenomenon was proposed by Harmanec (1985, see his Figs. 1 and 3)^{3}. The optical component in most of such systems is embedded in a dense disk. The spectral lines, seen in the optical and UV spectra and usually associated with this star, show complex profiles, often with emission wings. From that point of view, LS I +61° 303 is a typical Betype binary with spectral lines closely resembling socalled shell lines, which are identified in spectra of many interacting binary stars (see binary spectra published in Casares et al. 2005; Grundstrom et al. 2007; Aragona et al. 2009). Even the lines that do not show prominent emission are not purely atmospheric, but may be formed in the different parts of the optically thick disk or gaseous shell around the star. The RVs, measured for these lines (apart from the orbital motion) are affected by the complex kinematics of the gas motion in the different parts of the disk. Thus, they may not adequately reproduce the orbital motion of the Be star, and the solution of RV curves may result in spurious (and different for the different lines) values of the eccentricity of the orbit (Harmanec 1985; see also more recent papers by Harmanec et al. 2015; Koubský et al. 2019).
The complications with RV curves are well illustrated by the behavior of the shell lines in the Be binary star KX And (Stefl et al. 1990). The formal solutions of the RV curves obtained for the six different groups of lines formed in the vicinity of the primary give an orbit eccentricity value ranging from 0.22 to 0.64 (see Table 5 in Stefl et al. 1990). In contrast, the RV variations of the spectral lines of the latetype secondary, which were detected and studied with highresolution spectroscopy by Tarasov et al. (1998), revealed zero eccentricity and proved circular orbit in KX And.
We believe that the high value of eccentricity (e ≥ 0.5) obtained from the RV curves for the LS I +61° 303 (and perhaps for the LS 5039, see Aragona et al. 2009) is spurious. Of course, there is a morphological difference between the Be systems with normal (nondegenerate) secondary components and highmass X/γray binary systems with black hole/neutron star companions. However, the remarkable resemblance of the RV curves of the Be component in systems like KX And and many massive X/γray binaries, including LS I +61° 303 (see Fig. 3 in Harmanec 1985), cannot be ignored. Thus, the real orbit in LS I +61° 303 might be only slightly eccentric and close to circular. This is also suggested by the behavior of the periodic (orbitally phaselocked) variable polarization that we observed in the system.
5.2. Phase of the periastron
In addition to the high eccentricity, the RV measurements resulted in the value for the orbital phase of the periastron ϕ_{p} = 0.22−0.23 (Casares et al. 2005), ϕ_{p} = 0.3 (Grundstrom et al. 2007) and ϕ_{p} = 0.275 (Aragona et al. 2009). However, the peak on the light curves does not occur on this orbital phase at most wavelengths. For example, the peaks in the radio, Xray, and TeV light curves occur around the phase ϕ ≈ 0.6 (Massi et al. 2015; Archambault et al. 2016; Chernyakova et al. 2017). The peak of the optical brightness occurs with some delay relative to the peak in the highenergy emission (see Sect. 4, as well as Zamanov et al. 2014; ParedesFortuny et al. 2015). It is difficult to understand the nature of the phase delay between the passage of the periastron by a compact object and the maximum of radiation, especially in the case of high eccentricity. Our value for the orbital phase of periastron ϕ_{p} ≈ 0.6 (see Table 4 and Fig. 9 for orbital parameters and a possible geometry), obtained from the modeling of the optical polarization, may explain all observational facts mentioned above as an interaction of the compact object, moving around the Be star, with the densest part of the circumstellar disk, near the periastron. A new set of orbital parameters should now be considered when modeling light curves in different energy bands.
6. Summary
Our new highprecision BVR measurements of the linear polarization of the γray binary LS I +61° 303 revealed periodic orbital variability in all passbands. The timing analysis of the Stokes parameters yielded the first ever detection of a polarimetric period P_{Pol} = 13.244 d, which is close to half of the orbital period P_{orb} = 26.496 d. The continuous change of polarization with the orbital phase implies that the variations arise from the orbital motion of the compact star whose orbit is coplanar with the Be star’s decretion disk. The mechanism producing orbital variation of the polarization is most likely Thomson scattering of the stellar light in the hightemperature region around the compact object. The orbital variability curve is dominated by the second harmonic, which is typical for binaries with close to circular orbit and nearly symmetric distribution of the light scattering material with respect to the orbital plane. This implies that the high eccentricity of the binary orbit in LS I +61° 303 derived from the solution of the RV curves may not be real and that the true orbit is close to circular.
After the determination and subtraction of the interstellar polarization component, we obtained the PA of the constant component of polarization associated with the Be decretion disk at θ_{int} ≃ 11°. Although this value differs from the previously determined ≃25° (Nagae et al. 2006, 2009), we believe that the difference is most likely due to the uncertainties in the determinations of the interstellar polarization component.
We considered two cases: when the position angle of the disk normal nearly coincides with θ_{int}, and when it differs by 90° from it. Using a model of a scattering cloud at an elliptical orbit nearly coplanar with the disk, we modeled orbital variations of Stokes parameters and constrained the eccentricity at e < 0.15 and the phase of the periastron at ϕ_{p} ≈ 0.6. Both constraints are independent of the assumption on the disk orientation. The longitude of the periastron was found to be λ_{p} ≈ 146° or 225° for the two cases. The obtained value for the periastron phase differs significantly from the previously assumed phase, which is based on the RV measurements, but is very close to the peaks of the radio, Xray, and TeV light curves. Our results thus open a new avenue to model the broadband emission from the enigmatic γray binary LS I +61° 303.
We also found photometric orbital variability of LS I + 61° 303 in B, V, and R filters with amplitudes Δm ≈ 0.1 mag. The phase shift of the brightness maximum between the data sets acquired over the period of three years can be approximated with a simple linear model.
The paper by Simmons & Boyle (1984) contains an error in Eq. (A.2): inside the last square brackets, the “+” sign is missing between the terms 2ecos(λ − λ_{p}) and (e^{2}/2)cos[2(λ − λ_{p})].
The alternative model of massive Xray binaries suggested by Harmanec (1985), which assumes that Xray components in massive systems are not neutron stars or black holes, is outdated. However, his interpretation of the RV curves of spectral lines formed in the vicinity of the visible component in Be binary systems is correct at least in some cases.
Acknowledgments
We acknowledge support from the Magnus Ehrnrooth foundation, the Finnish National Agency for Education (EDUFI) Fellowship (VK), the Russian Science Foundation grant 201200364 (VK, SST and JP), the Academy of Finland grant 333112, and the ERC Advanced Grant HotMol ERC2011AdG291659 (SVB, AVB). We thank the Academy of Finland (projects 317552, 331951) and the German Academic Exchange Service (DAAD, projects 57405000 and 57525212) for travel grants. DM acknowledges support from the DFG grant MA 7807/21 and from the bwHPC project of the state of BadenWürttemberg. The Dipol2 polarimeter was built in cooperation by the University of Turku, Finland, and the Leibniz Institut für Sonnenphysik, Germany, with support from the Leibniz Association grant SAW2011KIS7. We are grateful to the Institute for Astronomy, University of Hawaii, for the allocated observing time.
References
 Abadi, M., Agarwal, A., Barham, P., et al. 2015, TensorFlow: LargeScale Machine Learning on Heterogeneous Systems, Software available from tensorflow.org [Google Scholar]
 Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, ApJ, 701, L123 [NASA ADS] [CrossRef] [Google Scholar]
 Ackermann, M., Ajello, M., Ballet, J., et al. 2013, ApJ, 773, L35 [NASA ADS] [CrossRef] [Google Scholar]
 Albert, J., Aliu, E., Anderhub, H., et al. 2006, Science, 312, 1771 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Aragona, C., McSwain, M. V., Grundstrom, E. D., et al. 2009, ApJ, 698, 514 [NASA ADS] [CrossRef] [Google Scholar]
 Archambault, S., Archer, A., Aune, T., et al. 2016, ApJ, 817, L7 [NASA ADS] [CrossRef] [Google Scholar]
 Aspin, C., Simmons, J. F. L., & Brown, J. C. 1981, MNRAS, 194, 283 [NASA ADS] [CrossRef] [Google Scholar]
 Astropy Collaboration (PriceWhelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
 Baluev, R. V. 2008, MNRAS, 385, 1279 [NASA ADS] [CrossRef] [Google Scholar]
 Berdyugin, A. V., & Tarasov, A. E. 1998, Astron. Lett., 24, 111 [Google Scholar]
 Brown, J. C., McLean, I. S., & Emslie, A. G. 1978, A&A, 68, 415 [NASA ADS] [Google Scholar]
 Brown, J. C., Aspin, C., Simmons, J. F. L., & McLean, I. S. 1982, MNRAS, 198, 787 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, A. G. A., Vallenari, A., Prusti, T., et al. 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Casares, J., Ribas, I., Paredes, J. M., Martí, J., & Allende Prieto, C. 2005, MNRAS, 360, 1105 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S., & Breen, F. H. 1947, ApJ, 105, 435 [NASA ADS] [CrossRef] [Google Scholar]
 Chernyakova, M., & Malyshev, D. 2020, Conf. Proc. Multifrequency Behaviour of High Energy Cosmic Sources  XIII, submitted [arXiv:2006.03615] [Google Scholar]
 Chernyakova, M., Neronov, A., Molkov, S., et al. 2012, ApJ, 747, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Chernyakova, M., Babyk, I., Malyshev, D., et al. 2017, MNRAS, 470, 1718 [NASA ADS] [CrossRef] [Google Scholar]
 Drissen, L., Lamontagne, R., Moffat, A. F. J., Bastien, P., & Seguin, M. 1986, ApJ, 304, 188 [NASA ADS] [CrossRef] [Google Scholar]
 Duane, S., Kennedy, A., Pendleton, B. J., & Roweth, D. 1987, Phys. Lett. B, 195, 216 [NASA ADS] [CrossRef] [Google Scholar]
 Golding, N. 2019, J. Open Sour. Softw., 4, 1601 [CrossRef] [Google Scholar]
 Gregory, P. C. 2002, ApJ, 575, 427 [NASA ADS] [CrossRef] [Google Scholar]
 Grundstrom, E. D., CaballeroNieves, S. M., Gies, D. R., et al. 2007, ApJ, 656, 437 [NASA ADS] [CrossRef] [Google Scholar]
 Harmanec, P. 1985, Bull. Astron. Inst. Czechoslov., 36, 327 [Google Scholar]
 Harmanec, P., Koubský, P., Nemravová, J. A., et al. 2015, A&A, 573, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Harrison, F. A., Ray, P. S., Leahy, D. A., Waltman, E. B., & Pooley, G. G. 2000, ApJ, 528, 454 [NASA ADS] [CrossRef] [Google Scholar]
 Heiles, C. 2000, AJ, 119, 923 [NASA ADS] [CrossRef] [Google Scholar]
 Hutchings, J. B., & Crampton, D. 1981, PASP, 93, 486 [NASA ADS] [CrossRef] [Google Scholar]
 Kosenkov, I. A., Berdyugin, A. V., Piirola, V., et al. 2017, MNRAS, 468, 4362 [NASA ADS] [CrossRef] [Google Scholar]
 Koubský, P., Harmanec, P., Brož, M., et al. 2019, A&A, 629, A105 [EDP Sciences] [Google Scholar]
 Leahy, D. A. 2001, Int. Cosm. Ray Conf., 6, 2524 [Google Scholar]
 Li, J., Torres, D. F., & Zhang, S. 2014, ApJ, 785, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Lomb, N. R. 1976, Ap&SS, 39, 447 [Google Scholar]
 Manset, N., & Bastien, P. 2000, AJ, 120, 413 [NASA ADS] [CrossRef] [Google Scholar]
 Massi, M., & Jaron, F. 2013, A&A, 554, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Massi, M., & TorricelliCiamponi, G. 2016, A&A, 585, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Massi, M., Ros, E., & Zimmermann, L. 2012, A&A, 540, A142 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Massi, M., Jaron, F., & Hovatta, T. 2015, A&A, 575, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mendelson, H., & Mazeh, T. 1989, MNRAS, 239, 733 [NASA ADS] [Google Scholar]
 Nagae, O., Kawabata, K. S., Fukazawa, Y., et al. 2006, PASJ, 58, 1015 [CrossRef] [Google Scholar]
 Nagae, O., Kawabata, K. S., Fukazawa, Y., et al. 2009, AJ, 137, 3509 [CrossRef] [Google Scholar]
 Paredes, J. M., Marziani, P., Marti, J., et al. 1994, A&A, 288, 519 [NASA ADS] [Google Scholar]
 ParedesFortuny, X., Ribó, M., BoschRamon, V., et al. 2015, A&A, 575, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Piirola, V., Berdyugin, A., & Berdyugina, S. 2014, in Groundbased and Airborne Instrumentation for Astronomy V, Proc. SPIE, 9147, 91478I [Google Scholar]
 Piirola, V., Berdyugin, A., Frisch, P. C., et al. 2020, A&A, 635, A46 [CrossRef] [EDP Sciences] [Google Scholar]
 Prusti, T., de Bruijne, J. H. J., Brown, A. G. A., et al. 2016, A&A, 595, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Quirrenbach, A., Bjorkman, K. S., Bjorkman, J. E., et al. 1997, ApJ, 479, 477 [NASA ADS] [CrossRef] [Google Scholar]
 R Core Team 2019, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria [Google Scholar]
 Scargle, J. D. 1982, ApJ, 263, 835 [Google Scholar]
 Serkowski, K. 1973, in Interstellar Dust and Related Topics, eds. J. M. Greenberg, & H. C. van de Hulst, IAU Symp., 52, 145 [CrossRef] [Google Scholar]
 Simmons, J. F. L., & Boyle, C. B. 1984, A&A, 134, 368 [NASA ADS] [Google Scholar]
 Simmons, J. F. L., Aspin, C., & Brown, J. C. 1982, MNRAS, 198, 45 [NASA ADS] [CrossRef] [Google Scholar]
 Sobolev, V. V. 1949, Uch. Zap. Leningrad Univ., 16 [Google Scholar]
 Stefl, S., Harmanec, P., Horn, J., et al. 1990, Bull. Astron. Inst. Czechoslov., 41, 29 [Google Scholar]
 Tarasov, A. E., Berdyugina, S. V., & Berdyugin, A. V. 1998, Astron. Lett., 24, 316 [Google Scholar]
 Taylor, A. R., & Gregory, P. C. 1982, ApJ, 255, 210 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, A. R., Kenny, H. T., Spencer, R. E., & Tzioumis, A. 1992, ApJ, 395, 268 [NASA ADS] [CrossRef] [Google Scholar]
 Wolinski, K. G., & Dolan, J. F. 1994, MNRAS, 267, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Wood, K., Bjorkman, J. E., Whitney, B., & Code, A. 1996, ApJ, 461, 847 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, Y. W., TorricelliCiamponi, G., Massi, M., et al. 2018, MNRAS, 474, 4245 [NASA ADS] [CrossRef] [Google Scholar]
 Zaitseva, G. V., & Borisov, G. V. 2003, Astron. Lett., 29, 188 [NASA ADS] [CrossRef] [Google Scholar]
 Zamanov, R. K., Martí, J., Paredes, J. M., et al. 1999, A&A, 351, 543 [NASA ADS] [Google Scholar]
 Zamanov, R., Stoyanov, K., Martí, J., et al. 2013, A&A, 559, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zamanov, R., Martí, J., Stoyanov, K., Borissova, A., & Tomov, N. A. 2014, A&A, 561, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Polarization from Thomson scattering by an orbiting cloud
We consider an orbital plane with the rotation axis inclined by angle i to the line of sight. Let us choose the polarization basis such that vector is along the projection of the vector on the plane of the sky, lies at the interception of the plane of the sky, and the orbital plane is perpendicular to . This basis can be supplemented by , which coincides with the direction to the observer to form a righthanded coordinate system. In this system, . The center of the coordinates is on the source of light (i.e., the Be star).
Let us introduce longitude λ in the orbital plane measured from the southern part of the line formed by the intersection of the orbital plane with the plane containing (, ), as in Simmons & Boyle (1984; see our Fig. A.1 and their Fig. 1). The orbital distance to the scattering cloud varies with longitude as
Fig. A.1. Geometry of orbit. 
where a is the orbit semimajor axis and λ_{p} is the longitude of the periastron. The unit vector towards the cloud is
and the scattering angle Θ is given by
The total flux observed at Earth from the system F_{tot} consists of the direct radiation from the Be star F_{*} = L_{*}/(4πD^{2}) and the scattered flux F_{sc} = L_{sc}(μ)/(4πD^{2}), where L_{*} and L_{sc} are corresponding luminosities and D is the distance. Assuming that Thomson scattering (in an optically thin regime) is responsible for scattered flux, the angular distribution of scattered luminosity can be represented as
where
is the Thomson scattering indicatrix,
is the fraction of scattered photons that is related to the total number of free electrons N_{e} in a cloud and the distance between the Be star and the scattering cloud, and σ_{T} is the Thomson crosssection. We can scale f_{sc} to some typical value f_{0} as
The polarization degree of scattered light is
but the observed one is diluted by the star. The polarized flux for the assumed Thomson scattering indicatrix is
The polarization degree of the total flux is P = F_{sc}P_{sc}/F_{tot}, but if we assume that scattered radiation contributed very little to the total flux F_{sc} ≪ F_{*}, we get
The polarization (pseudo)vector , defined by the direction of dominant oscillations of electromagnetic waves, is perpendicular to the scattering plane, so
and . The position angle χ of is given by the following formulae:
Because the normalized Stokes parameters are defined as q = Pcos(2χ) and u = Psin(2χ), we need
Combining that with expression (A.10) for polarization degree P, we get
The observed Stokes parameters also depend on the position angle Ω of the projection of the orbital axis, which are obtained by rotating vector (q, u) by angle 2Ω:
The obtained Stokes parameters are functions of the longitude λ and need to be computed as functions of the orbital phase. From the true anomaly of the orbit λ − λ_{p}, we find the eccentric anomaly
and the mean anomaly
which can be converted to the orbital phase measured in the interval [0,1] using an additional free parameter, the phase of the periastron ϕ_{p}:
Thus, the model free parameters are the inclination i, eccentricity e, latitude of the periastron λ_{p}, the PA of the projection of the orbit axis on the sky Ω, the phase of the periastron ϕ_{p}, and the scattering fraction f_{0}. We note that Ω defined this way differs by π/2 from the PA of the ascending node usually used to describe binary star orbits. In order to fit the data, there is a need for additional constant Stokes parameters q_{0} and u_{0}, which describe permanent polarization and are not related to the orbital motion of the scattering cloud.
All Tables
Observed PD and PA of LS I +61° 303, interstellar polarization and average intrinsic polarization.
All Figures
Fig. 1. Polarization map of LS I +61° 303 (red circle at the origin) and field stars (black circles) in the R band. The length of the bars corresponds to the degree of linear polarization P, and the direction corresponds to the PA θ (measured from the north to the east). 

In the text 
Fig. 2. Dependence of observed PD P and PA θ on parallax (Gaia DR2) for LS I +61° 303 (red) and field stars (black) in the B band. The numbering of the field stars is the same as in Fig. 1. The error bars correspond to the 1σ errors. 

In the text 
Fig. 3. Dependence of linear PD P (Heiles 2000) on distance (evaluated from the inverse of the Gaia DR2 parallaxes) for the stars in the 10° ×10° area around LS I +61° 303 in the R band. The red star indicates our value of the observed polarization of LS I +61° 303, the orange crosses correspond to the four nearby field stars shown in Fig. 1. The vertical bar corresponds to the 1σ error in polarization. 

In the text 
Fig. 4. Observed Stokes parameters q (top panel) and u (bottom panel) with 1σ errors for LS I +61° 303 in the V band, measured from 2016 September 15 to 2019 September 1 (MJD 57646–58728). 

In the text 
Fig. 5. Top panel: LombScargle periodogram for normalized Stokes parameter u in the B band. Bottom panel: LombScargle periodogram zoom for the normalized Stokes parameter u in the B, V, and R bands (solid blue, dashed green, and dotted red curves, respectively). The horizontal dashed line corresponds to false alarm probability (FAP) = 1%. The vertical red line corresponds to one half of the orbital period P_{orb} = 26.496 d. 

In the text 
Fig. 6. Orbital variability of intrinsic PD P_{int} and PA θ_{int} of LS I +61° 303 in V band. The filled blue squares with 1σ errors correspond to the average values of the individual observations (represented by gray circles with error bars) and the standard errors of the mean calculated within the phase bins of width Δϕ = 0.091. The vertical dashed lines correspond to the phases of periastron (ϕ = 0.275) and apastron (ϕ = 0.775) as derived by Aragona et al. (2009). 

In the text 
Fig. 7. Variability of intrinsic Stokes parameters q_{int} (left column) and u_{int} (right column), in B, V, and R bands (top, middle and bottom panels, respectively) with the orbital phase folded with the period P_{orb} = 26.496 d. The filled squares with 1σ errors correspond to the average values of the individual observations (gray circles with error bars) and the standard errors of the mean calculated within phase bins of width Δϕ = 0.091. The vertical dashed lines correspond to the phases of periastron (ϕ = 0.275) and apastron (ϕ = 0.775) as derived by Aragona et al. (2009). The red dashed lines show the best fit with Fourier series given by Eq. (1). The black solid and blue dotted lines correspond to the best fit model of a scattering cloud on an eccentric orbit from Sect. 3.3, for an Ω constrained close to θ_{int} and with a free Ω, respectively. 

In the text 
Fig. 8. Posterior distributions for parameters of model described in Sect. 3.3 and Appendix A for Ω = θ_{int} ± 10°. Diagonal panels: distributions of model parameters. The blue solid, green dashed, and red dotdashed lines correspond to perpassband B, V, and R distributions, respectively. Lowertriangle panels: joint posterior distributions of two parameters. The green, blue, and violet contours correspond to 0.68, 0.95, and 0.997 probability levels of parameters i, Ω, λ_{p}, and ϕ_{p}. The shades of blue, green, and red correspond to the same probability levels for B, V, and R filters, respectively. 

In the text 
Fig. 9. Relative orbit of a compact object in LS I +61° 303 around Be star (yellow circle at the origin), which lies at the ellipse focus. The orbital parameters are taken from Table 4. The red dashed line is the major axis of the orbit. The black dots on the ellipse are spaced by the Δϕ = 0.1. Because of a 180° degeneracy in λ_{p}, the orbit can be also rotated by 180°. Left panel: orbit for constrained orientation relative to the decretion disk plane with Ω = θ_{int} ± 10°. Right panel: orbit assuming free Ω. 

In the text 
Fig. 10. Same as Fig. 8, but for free Ω. 

In the text 
Fig. 11. Variation of relative amplitude of LS I +61° 303 around the mean in the V passband as a function of superorbital phase. The shaded blue bands correspond to the eleven observing seasons (S1–S11). 

In the text 
Fig. 12. Top panel: LombScargle periodogram for brightness variability of LS I +61° 303 in the R band. Bottom panel: LombScargle periodogram zoom for brightness variability in the BVR bands (dotted blue, dashed green, and solid red lines, respectively). The horizontal dotted line corresponds to the FAP = 1%. The vertical red dotted line marks the orbital period P_{orb} = 26.496 d. 

In the text 
Fig. 13. Variability of LS I +61° 303 brightness in BVR bands (top, middle, and bottom panels, respectively) with the orbital period. The filled squares with 1σ errors correspond to the average values of the individual observations (gray crosses) and the standard errors of the mean calculated within the phase bin of width Δϕ = 0.091. 

In the text 
Fig. 14. Variability of LS I +61° 303 optical brightness in Vband for different parts of the data (seasons S1–S11 from top to bottom). The black solid lines correspond to the best fit of the data with function (6). The vertical dashed lines give the phases of the best fit maxima ϕ_{0}, and the corresponding ±1σ, ±2σ, and ±3σ confidence intervals are shown with varying shades of blue. The vertical scale [0.1, −0.1] is the same in all panels. 

In the text 
Fig. 15. Dependence of brightness maximum orbital phase ϕ_{0} (top panel) and the amplitude A (bottom panel) of the sinusoidal fits on the superorbital phase. The solid blue lines correspond to the linear fit of the data, while the ±1σ and ±3σ confidence intervals are shown in dark and light blue. The red squares show the parameters from Table 3 of ParedesFortuny et al. (2015). 

In the text 
Fig. A.1. Geometry of orbit. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.