© ESO 2020

1. Introduction

Gamma-ray binaries constitute a subclass of high-mass 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 best-studied γ-ray binaries and was observed in the last few decades over the whole range of the electromagnetic spectrum, from the radio to the very high-energy γ-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₁ = 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 Lomb-Scargle timing analysis of 37 years of radio data resulted in the detection of the second period, P₂ = 26.935 ± 0.013 d (Massi & Torricelli-Ciamponi 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₂ = 26.926 ± 0.005 d (Wu et al. 2018). The beat of these two periods $P_{\text{beat}}={(P_1^{-1}-P_2^{-1})}^{-1}\approx 1660$ d is very close to the period of the long-term superorbital variability P_sup ≈ 1700 d observed in the radio (Massi & Jaron 2013), Hα emission (Zamanov et al. 2013), X-rays (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 two-year 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 co-aligned with the orbit of the compact companion (Nagae et al. 2009).

In this paper, we present the results of our BVR high-precision 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 broad-band BVR polarimeter Dipol-2 (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 R-passbands. The total integration time with a typical ten-second 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⁻⁵ for all passbands and is measured with the accuracy of a few parts per million (10⁻⁶). 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.

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%, ${\theta }_{\mathrm{is}}=126\stackrel{°}{.}5\pm 3\stackrel{°}{.}7$.

The values of the average intrinsic polarization of LS I +61° 303 in the B, V, and R-bands 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 wavelength-dependent PD, which is an indication of the important role of bound-free 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 Lomb-Scargle method¹ of spectral analysis (Lomb 1976; Scargle 1982), using the ASTROPY package (Astropy Collaboration 2018). The Lomb-Scargle 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⁻⁵, 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⁴ 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⁴ 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: Lomb-Scargle periodogram for normalized Stokes parameter u in the B band. Bottom panel: Lomb-Scargle 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 Porb = 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 small-amplitude variations in the observed polarization of LS I +61° 303. Moreover, a continuous change of polarization implies co-planar 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 short-term changes in the distribution and density of the light scattering material in the Be decretion disk.

Fig. 6. Orbital variability of intrinsic PD Pint 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 qint (left column) and uint (right column), in B, V, and R bands (top, middle and bottom panels, respectively) with the orbital phase folded with the period Porb = 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 co-rotating 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 co-rotating light scattering envelope can be represented through a Fourier series of the orbital longitude λ = 2πϕ (where ϕ is a phase of the orbital period):

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill q_{\mathrm{int}}& =q_0+q_{1}cos\lambda +q_{2}sin\lambda +q_{3}cos2\lambda +q_{4}sin2\lambda ,\hfill \\ \hfill u_{\mathrm{int}}& =u_0+u_{1}cos\lambda +u_{2}sin\lambda +u_{3}cos2\lambda +u_{4}sin2\lambda .\hfill \end{array}\end{array}$$(1)

We fit the intrinsic Stokes parameters q_int and u_int with these functions, using q₀–q₄ and u₀–u₄ as free parameters. The best fit parameters together with the χ² 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 Lomb-Scargle periodogram close to the frequency corresponding to one half of the orbital period. The parameters

$$\begin{array}{c}\hfill A_q=\sqrt{\frac{q_4^2+q_3^2}{q_2^2+q_1^2}},A_u=\sqrt{\frac{u_4^2+u_3^2}{u_2^2+u_1^2}},\end{array}$$(2)

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):

$$\begin{array}{cc}& {\left(\frac{1-cosi}{1+cosi}\right)}^4=\frac{{(u_3+q_4)}^2+{(u_4-q_3)}^2}{{(u_4+q_3)}^2+{(u_3-q_4)}^2}\hfill \end{array}$$(3)

$$\begin{array}{cc}& tan2\mathrm{\Omega }=\frac{A+B}{C+D},\hfill \end{array}$$(4)

where

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill A& =\frac{u_4-q_3}{{(1-cosi)}^2},B=\frac{u_4+q_3}{{(1+cosi)}^2},\hfill \\ \hfill C& =\frac{q_4-u_3}{{(1+cosi)}^2},D=\frac{u_3+q_4}{{(1-cosi)}^2}·\hfill \end{array}\end{array}$$(5)

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 straight-line 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₀ (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 pseudo-vector 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 pseudo-vectors (q₀, u₀).

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₀, u₀) and a typical scattering fraction f₀ 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 χ²/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 χ²/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 co-planar (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 χ² 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 dot-dashed lines correspond to per-passband B, V, and R distributions, respectively. Lower-triangle 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)² 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 χ² 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 χ² 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 χ² value of the fit is much smaller than in the case of constrained Ω (50 versus 67) and the best-fit Ω ≈ 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 co-planar 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.

Fig. 10. Same as Fig. 8, but for free Ω.

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 Lomb-Scargle 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 Lomb-Scargle 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⁻⁵, 10⁻¹⁰, and 10⁻¹⁴, 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: Lomb-Scargle periodogram for brightness variability of LS I +61° 303 in the R band. Bottom panel: Lomb-Scargle 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 Porb = 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 (sinusoidal-like 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 Paredes-Fortuny 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 Paredes-Fortuny 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

$$\begin{array}{c}\hfill m(\varphi )=-Acos[2\pi (\varphi -{\varphi }_0)]+m_0,\end{array}$$(6)

where ϕ₀ (phase of the peak), A (amplitude), and m₀ (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 χ² ≈ 1).

Fig. 14. Variability of LS I +61° 303 optical brightness in V-band 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.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 Paredes-Fortuny 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 small-amplitude 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 Paredes-Fortuny 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 Paredes-Fortuny 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.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 X-ray and non-X-ray 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)³. 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 Be-type binary with spectral lines closely resembling so-called 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 late-type secondary, which were detected and studied with high-resolution 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 high-mass 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 phase-locked) 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, X-ray, 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 high-energy emission (see Sect. 4, as well as Zamanov et al. 2014; Paredes-Fortuny 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 high-precision 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 co-planar 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 high-temperature 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 co-planar 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, X-ray, and TeV light curves. Our results thus open a new avenue to model the broad-band 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.

Acknowledgments

We acknowledge support from the Magnus Ehrnrooth foundation, the Finnish National Agency for Education (EDUFI) Fellowship (VK), the Russian Science Foundation grant 20-12-00364 (VK, SST and JP), the Academy of Finland grant 333112, and the ERC Advanced Grant HotMol ERC-2011-AdG-291659 (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/2-1 and from the bwHPC project of the state of Baden-Württemberg. The Dipol-2 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 SAW-2011-KIS-7. We are grateful to the Institute for Astronomy, University of Hawaii, for the allocated observing time.

References

Appendix A: Polarization from Thomson scattering by an orbiting cloud

We consider an orbital plane with the rotation axis $\widehat{\mathit{n}}$ inclined by angle i to the line of sight. Let us choose the polarization basis $({\widehat{\mathit{e}}}_1,{\widehat{\mathit{e}}}_2)$ such that vector ${\widehat{\mathit{e}}}_1$ is along the projection of the vector $\widehat{\mathit{n}}$ on the plane of the sky, ${\widehat{\mathit{e}}}_2$ lies at the interception of the plane of the sky, and the orbital plane is perpendicular to ${\widehat{\mathit{e}}}_1$. This basis can be supplemented by ${\widehat{\mathit{e}}}_3$, which coincides with the direction to the observer $\widehat{\mathit{o}}$ to form a right-handed coordinate system. In this system, $\widehat{\mathit{n}}=(sini,0,cosi)$. 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 ($\widehat{\mathit{n}}$, $\widehat{\mathit{o}}$), 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

$$\begin{array}{c}\hfill r(\lambda )=\frac{a(1-e^2)}{1+ecos(\lambda -{\lambda }_{\mathrm{p}})},\end{array}$$(A.1)

Fig. A.1. Geometry of orbit.

where a is the orbit semi-major axis and λ_p is the longitude of the periastron. The unit vector towards the cloud is

$$\begin{array}{c}\hfill \widehat{\mathit{r}}=(-cosicos\lambda ,-sin\lambda ,sinicos\lambda ),\end{array}$$(A.2)

and the scattering angle Θ is given by

$$\begin{array}{c}\hfill \mu =cos\mathrm{\Theta }=\widehat{\mathit{r}}·\widehat{\mathit{o}}=sinicos\lambda .\end{array}$$(A.3)

The total flux observed at Earth from the system F_tot consists of the direct radiation from the Be star F_* = L_*/(4πD²) and the scattered flux F_sc = L_sc(μ)/(4πD²), 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

$$\begin{array}{c}\hfill L_{\mathrm{sc}}(\mu )=L_{\ast }{f}_{\mathrm{sc}}l(\mu ),\end{array}$$(A.4)

where

$$\begin{array}{c}\hfill l(\mu )=\frac{3}{8}(1+{\mu }^2)\end{array}$$(A.5)

is the Thomson scattering indicatrix,

$$\begin{array}{c}\hfill f_{\mathrm{sc}}=\frac{N_{\mathrm{e}}{\sigma }_{\mathrm{T}}}{4\pi r^2(\lambda )}\end{array}$$(A.6)

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 cross-section. We can scale f_sc to some typical value f₀ as

$$\begin{array}{c}\hfill f_{\mathrm{sc}}=f_0{\left[\frac{a(1-e^2)}{r(\lambda )}\right]}^2=f_0{[1+ecos(\lambda -{\lambda }_{\mathrm{p}})]}^2.\end{array}$$(A.7)

The polarization degree of scattered light is

$$\begin{array}{c}\hfill P_{\mathrm{sc}}=\frac{1-{\mu }^2}{1+{\mu }^2}=\frac{{sin}^2\mathrm{\Theta }}{1+{cos}^2\mathrm{\Theta }},\end{array}$$(A.8)

but the observed one is diluted by the star. The polarized flux for the assumed Thomson scattering indicatrix is

$$\begin{array}{c}\hfill F_{\mathrm{sc}}{P}_{\mathrm{sc}}=F_{\ast }{f}_{\mathrm{sc}}l(\mu )\frac{1-{\mu }^2}{1+{\mu }^2}=F_{\ast }{f}_{\mathrm{sc}}\frac{3}{8}(1-{\mu }^2).\end{array}$$(A.9)

The polarization degree of the total flux is P = F_scP_sc/F_tot, but if we assume that scattered radiation contributed very little to the total flux F_sc ≪ F_*, we get

$$\begin{array}{c}\hfill P=\frac{3}{8}{f}_{\mathrm{sc}}(1-{\mu }^2)=\frac{3}{8}{f}_0{sin}^2\mathrm{\Theta }{[1+ecos(\lambda -{\lambda }_{\mathrm{p}})]}^2.\end{array}$$(A.10)

The polarization (pseudo-)vector $\widehat{\mathit{p}}$, defined by the direction of dominant oscillations of electromagnetic waves, is perpendicular to the scattering plane, so

$$\begin{array}{c}\hfill \widehat{\mathit{p}}=\frac{\widehat{\mathit{o}}\times \widehat{\mathit{r}}}{|\widehat{\mathit{o}}\times \widehat{\mathit{r}}|}=\frac{1}{sin\mathrm{\Theta }}(sin\lambda ,-cosicos\lambda ,0),\end{array}$$(A.11)

and $sin\mathrm{\Theta }=\sqrt{1-{\mu }^2}=\sqrt{1-{sin}^{2}i{cos}^2\lambda }$. The position angle χ of $\widehat{\mathit{p}}$ is given by the following formulae:

$$\begin{array}{cc}\hfill cos\chi & ={\widehat{\mathit{e}}}_1·\widehat{\mathit{p}}=\frac{sin\lambda }{sin\mathrm{\Theta }},\hfill \end{array}$$(A.12)

$$\begin{array}{cc}\hfill sin\chi & ={\widehat{\mathit{e}}}_2·\widehat{\mathit{p}}=-\frac{cosicos\lambda }{sin\mathrm{\Theta }}·\hfill \end{array}$$(A.13)

Because the normalized Stokes parameters are defined as q = Pcos(2χ) and u = Psin(2χ), we need

$$\begin{array}{cc}\hfill cos(2\chi )& =\frac{{sin}^{2}i-(1+{cos}^{2}i)cos(2\lambda )}{2{sin}^2\mathrm{\Theta }},\hfill \end{array}$$(A.14)

$$\begin{array}{cc}\hfill sin(2\chi )& =-\frac{cosisin(2\lambda )}{{sin}^2\mathrm{\Theta }}·\hfill \end{array}$$(A.15)

Combining that with expression (A.10) for polarization degree P, we get

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill q& =\frac{3f_0}{16}[{sin}^{2}i-(1+{cos}^{2}i)cos2\lambda ]{[1+ecos(\lambda -{\lambda }_{\mathrm{p}})]}^2,\hfill \\ \hfill u& =\frac{3f_0}{8}[-cosisin2\lambda ]{[1+ecos(\lambda -{\lambda }_{\mathrm{p}})]}^2.\hfill \end{array}\end{array}$$(A.16)

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Ω:

$$\begin{array}{cc}\hfill q_{\mathrm{obs}}& =qcos(2\mathrm{\Omega })-usin(2\mathrm{\Omega }),\hfill \end{array}$$(A.17)

$$\begin{array}{cc}\hfill u_{\mathrm{obs}}& =qsin(2\mathrm{\Omega })+ucos(2\mathrm{\Omega }).\hfill \end{array}$$(A.18)

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

$$\begin{array}{c}\hfill tan\left(\frac{E}{2}\right)=\sqrt{\frac{1-e}{1+e}}tan\left(\frac{\lambda -{\lambda }_{\mathrm{p}}}{2}\right)\end{array}$$(A.19)

and the mean anomaly

$$\begin{array}{c}\hfill M=E-esinE,\end{array}$$(A.20)

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:

$$\begin{array}{c}\hfill {\varphi }_{\mathrm{orb}}=M/(2\pi )+{\varphi }_{\mathrm{p}}.\end{array}$$(A.21)

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₀. 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₀ and u₀, which describe permanent polarization and are not related to the orbital motion of the scattering cloud.

All Tables

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: Lomb-Scargle periodogram for normalized Stokes parameter u in the B band. Bottom panel: Lomb-Scargle 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 Porb = 26.496 d.

In the text

Fig. 6. Orbital variability of intrinsic PD Pint 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 qint (left column) and uint (right column), in B, V, and R bands (top, middle and bottom panels, respectively) with the orbital phase folded with the period Porb = 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 dot-dashed lines correspond to per-passband B, V, and R distributions, respectively. Lower-triangle 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: Lomb-Scargle periodogram for brightness variability of LS I +61° 303 in the R band. Bottom panel: Lomb-Scargle 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 Porb = 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 V-band 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 Paredes-Fortuny et al. (2015).
In the text

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