Free Access
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/0004-6361/202038745
Published online 20 November 2020

© 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 P1 = 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, P2 = 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 P2 = 26.926 ± 0.005 d (Wu et al. 2018). The beat of these two periods d is very close to the period of the long-term superorbital variability Psup ≈ 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−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) Pobs and the PA θobs.

Table 1.

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.

thumbnail 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, Pint 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.

thumbnail 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.

thumbnail 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 Pis and θis of the IS polarization in all passbands and the average values Pint 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: Pis = 2.20 ± 0.18%, .

Table 2.

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 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 method1 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 Porb = 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 PV = 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 104 light curves, varying the observed data points around the mean values in accordance with the observational errors. In the second case, we also simulated 104 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 Ppol = 13.244 ± 0.012 d, which differs by less than 1σ from Porb/2. Thus, for the very first time, the orbital period of LS I +61° 303 was obtained directly from the polarimetric measurements.

thumbnail 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).

thumbnail 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 Pint and the PA θint of the intrinsic polarization (Fig. 6) and for the normalized Stokes parameters qint and uint (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.

thumbnail 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).

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

(1)

We fit the intrinsic Stokes parameters qint and uint with these functions, using q0q4 and u0u4 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 Lomb-Scargle periodogram close to the frequency corresponding to one half of the orbital period. The parameters

(2)

Table 3.

Fourier coefficients and their errors of the best fits to the data in BVR bands with Eq. (1).

giving the ratio of the amplitudes of the second to the first harmonic, attain the values Aq = 0.9, 1.0, 3.3 and Au = 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):

(3)

(4)

where

(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 f0 (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 (q0, u0).

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 (q0, u0) and a typical scattering fraction f0 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 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 χ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.

thumbnail 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.

thumbnail 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 Ω.

Table 4.

Best fit parameter estimates for the model described in Sect. 3.3 and Appendix A.

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 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.

thumbnail 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.

thumbnail 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 Porb = 26.4960 ± 0.0028 d. For the BVR bands, the determined periods are PB = 26.56 ± 0.05, PV = 26.58 ± 0.04, and PR = 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).

thumbnail 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 Porb 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.

thumbnail 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 Porb and fit with the function

(6)

Table 5.

Parameters of the best fit to the optical photometry of LS I +61° 303 in the V band with function (6) for eleven observing seasons.

where ϕ0 (phase of the peak), A (amplitude), and m0 (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).

thumbnail 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 ≈ 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.

thumbnail 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 = 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)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 PPol = 13.244 d, which is close to half of the orbital period Porb = 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.


2

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 (e2/2)cos[2(λ − λp)].

3

The alternative model of massive X-ray binaries suggested by Harmanec (1985), which assumes that X-ray 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 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

  1. Abadi, M., Agarwal, A., Barham, P., et al. 2015, TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems, Software available from tensorflow.org [Google Scholar]
  2. Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, ApJ, 701, L123 [NASA ADS] [CrossRef] [Google Scholar]
  3. Ackermann, M., Ajello, M., Ballet, J., et al. 2013, ApJ, 773, L35 [NASA ADS] [CrossRef] [Google Scholar]
  4. Albert, J., Aliu, E., Anderhub, H., et al. 2006, Science, 312, 1771 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  5. Aragona, C., McSwain, M. V., Grundstrom, E. D., et al. 2009, ApJ, 698, 514 [NASA ADS] [CrossRef] [Google Scholar]
  6. Archambault, S., Archer, A., Aune, T., et al. 2016, ApJ, 817, L7 [NASA ADS] [CrossRef] [Google Scholar]
  7. Aspin, C., Simmons, J. F. L., & Brown, J. C. 1981, MNRAS, 194, 283 [NASA ADS] [CrossRef] [Google Scholar]
  8. Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
  9. Baluev, R. V. 2008, MNRAS, 385, 1279 [NASA ADS] [CrossRef] [Google Scholar]
  10. Berdyugin, A. V., & Tarasov, A. E. 1998, Astron. Lett., 24, 111 [Google Scholar]
  11. Brown, J. C., McLean, I. S., & Emslie, A. G. 1978, A&A, 68, 415 [NASA ADS] [Google Scholar]
  12. Brown, J. C., Aspin, C., Simmons, J. F. L., & McLean, I. S. 1982, MNRAS, 198, 787 [NASA ADS] [CrossRef] [Google Scholar]
  13. Brown, A. G. A., Vallenari, A., Prusti, T., et al. 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Casares, J., Ribas, I., Paredes, J. M., Martí, J., & Allende Prieto, C. 2005, MNRAS, 360, 1105 [NASA ADS] [CrossRef] [Google Scholar]
  15. Chandrasekhar, S., & Breen, F. H. 1947, ApJ, 105, 435 [NASA ADS] [CrossRef] [Google Scholar]
  16. Chernyakova, M., & Malyshev, D. 2020, Conf. Proc. Multifrequency Behaviour of High Energy Cosmic Sources - XIII, submitted [arXiv:2006.03615] [Google Scholar]
  17. Chernyakova, M., Neronov, A., Molkov, S., et al. 2012, ApJ, 747, L29 [NASA ADS] [CrossRef] [Google Scholar]
  18. Chernyakova, M., Babyk, I., Malyshev, D., et al. 2017, MNRAS, 470, 1718 [NASA ADS] [CrossRef] [Google Scholar]
  19. Drissen, L., Lamontagne, R., Moffat, A. F. J., Bastien, P., & Seguin, M. 1986, ApJ, 304, 188 [NASA ADS] [CrossRef] [Google Scholar]
  20. Duane, S., Kennedy, A., Pendleton, B. J., & Roweth, D. 1987, Phys. Lett. B, 195, 216 [NASA ADS] [CrossRef] [Google Scholar]
  21. Golding, N. 2019, J. Open Sour. Softw., 4, 1601 [CrossRef] [Google Scholar]
  22. Gregory, P. C. 2002, ApJ, 575, 427 [NASA ADS] [CrossRef] [Google Scholar]
  23. Grundstrom, E. D., Caballero-Nieves, S. M., Gies, D. R., et al. 2007, ApJ, 656, 437 [NASA ADS] [CrossRef] [Google Scholar]
  24. Harmanec, P. 1985, Bull. Astron. Inst. Czechoslov., 36, 327 [Google Scholar]
  25. Harmanec, P., Koubský, P., Nemravová, J. A., et al. 2015, A&A, 573, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Harrison, F. A., Ray, P. S., Leahy, D. A., Waltman, E. B., & Pooley, G. G. 2000, ApJ, 528, 454 [NASA ADS] [CrossRef] [Google Scholar]
  27. Heiles, C. 2000, AJ, 119, 923 [NASA ADS] [CrossRef] [Google Scholar]
  28. Hutchings, J. B., & Crampton, D. 1981, PASP, 93, 486 [NASA ADS] [CrossRef] [Google Scholar]
  29. Kosenkov, I. A., Berdyugin, A. V., Piirola, V., et al. 2017, MNRAS, 468, 4362 [NASA ADS] [CrossRef] [Google Scholar]
  30. Koubský, P., Harmanec, P., Brož, M., et al. 2019, A&A, 629, A105 [EDP Sciences] [Google Scholar]
  31. Leahy, D. A. 2001, Int. Cosm. Ray Conf., 6, 2524 [Google Scholar]
  32. Li, J., Torres, D. F., & Zhang, S. 2014, ApJ, 785, L19 [NASA ADS] [CrossRef] [Google Scholar]
  33. Lomb, N. R. 1976, Ap&SS, 39, 447 [Google Scholar]
  34. Manset, N., & Bastien, P. 2000, AJ, 120, 413 [NASA ADS] [CrossRef] [Google Scholar]
  35. Massi, M., & Jaron, F. 2013, A&A, 554, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Massi, M., & Torricelli-Ciamponi, G. 2016, A&A, 585, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Massi, M., Ros, E., & Zimmermann, L. 2012, A&A, 540, A142 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Massi, M., Jaron, F., & Hovatta, T. 2015, A&A, 575, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Mendelson, H., & Mazeh, T. 1989, MNRAS, 239, 733 [NASA ADS] [Google Scholar]
  40. Nagae, O., Kawabata, K. S., Fukazawa, Y., et al. 2006, PASJ, 58, 1015 [CrossRef] [Google Scholar]
  41. Nagae, O., Kawabata, K. S., Fukazawa, Y., et al. 2009, AJ, 137, 3509 [CrossRef] [Google Scholar]
  42. Paredes, J. M., Marziani, P., Marti, J., et al. 1994, A&A, 288, 519 [NASA ADS] [Google Scholar]
  43. Paredes-Fortuny, X., Ribó, M., Bosch-Ramon, V., et al. 2015, A&A, 575, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Piirola, V., Berdyugin, A., & Berdyugina, S. 2014, in Ground-based and Airborne Instrumentation for Astronomy V, Proc. SPIE, 9147, 91478I [Google Scholar]
  45. Piirola, V., Berdyugin, A., Frisch, P. C., et al. 2020, A&A, 635, A46 [CrossRef] [EDP Sciences] [Google Scholar]
  46. 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]
  47. Quirrenbach, A., Bjorkman, K. S., Bjorkman, J. E., et al. 1997, ApJ, 479, 477 [NASA ADS] [CrossRef] [Google Scholar]
  48. R Core Team 2019, R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria [Google Scholar]
  49. Scargle, J. D. 1982, ApJ, 263, 835 [Google Scholar]
  50. 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]
  51. Simmons, J. F. L., & Boyle, C. B. 1984, A&A, 134, 368 [NASA ADS] [Google Scholar]
  52. Simmons, J. F. L., Aspin, C., & Brown, J. C. 1982, MNRAS, 198, 45 [NASA ADS] [CrossRef] [Google Scholar]
  53. Sobolev, V. V. 1949, Uch. Zap. Leningrad Univ., 16 [Google Scholar]
  54. Stefl, S., Harmanec, P., Horn, J., et al. 1990, Bull. Astron. Inst. Czechoslov., 41, 29 [Google Scholar]
  55. Tarasov, A. E., Berdyugina, S. V., & Berdyugin, A. V. 1998, Astron. Lett., 24, 316 [Google Scholar]
  56. Taylor, A. R., & Gregory, P. C. 1982, ApJ, 255, 210 [NASA ADS] [CrossRef] [Google Scholar]
  57. Taylor, A. R., Kenny, H. T., Spencer, R. E., & Tzioumis, A. 1992, ApJ, 395, 268 [NASA ADS] [CrossRef] [Google Scholar]
  58. Wolinski, K. G., & Dolan, J. F. 1994, MNRAS, 267, 5 [NASA ADS] [CrossRef] [Google Scholar]
  59. Wood, K., Bjorkman, J. E., Whitney, B., & Code, A. 1996, ApJ, 461, 847 [NASA ADS] [CrossRef] [Google Scholar]
  60. Wu, Y. W., Torricelli-Ciamponi, G., Massi, M., et al. 2018, MNRAS, 474, 4245 [NASA ADS] [CrossRef] [Google Scholar]
  61. Zaitseva, G. V., & Borisov, G. V. 2003, Astron. Lett., 29, 188 [NASA ADS] [CrossRef] [Google Scholar]
  62. Zamanov, R. K., Martí, J., Paredes, J. M., et al. 1999, A&A, 351, 543 [NASA ADS] [Google Scholar]
  63. Zamanov, R., Stoyanov, K., Martí, J., et al. 2013, A&A, 559, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. 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 right-handed 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

(A.1)

thumbnail 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

(A.2)

and the scattering angle Θ is given by

(A.3)

The total flux observed at Earth from the system Ftot consists of the direct radiation from the Be star F* = L*/(4πD2) and the scattered flux Fsc = Lsc(μ)/(4πD2), where L* and Lsc 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

(A.4)

where

(A.5)

is the Thomson scattering indicatrix,

(A.6)

is the fraction of scattered photons that is related to the total number of free electrons Ne in a cloud and the distance between the Be star and the scattering cloud, and σT is the Thomson cross-section. We can scale fsc to some typical value f0 as

(A.7)

The polarization degree of scattered light is

(A.8)

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

(A.9)

The polarization degree of the total flux is P = FscPsc/Ftot, but if we assume that scattered radiation contributed very little to the total flux Fsc ≪ F*, we get

(A.10)

The polarization (pseudo-)vector , defined by the direction of dominant oscillations of electromagnetic waves, is perpendicular to the scattering plane, so

(A.11)

and . The position angle χ of is given by the following formulae:

(A.12)

(A.13)

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

(A.14)

(A.15)

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

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

(A.17)

(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

(A.19)

and the mean anomaly

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

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

All Tables

Table 1.

Log of polarimetric observations of LS I +61° 303.

Table 2.

Observed PD and PA of LS I +61° 303, interstellar polarization and average intrinsic polarization.

Table 3.

Fourier coefficients and their errors of the best fits to the data in BVR bands with Eq. (1).

Table 4.

Best fit parameter estimates for the model described in Sect. 3.3 and Appendix A.

Table 5.

Parameters of the best fit to the optical photometry of LS I +61° 303 in the V band with function (6) for eleven observing seasons.

All Figures

thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail Fig. 10.

Same as Fig. 8, but for free Ω.

In the text
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail Fig. A.1.

Geometry of orbit.

In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.