© The Authors 2023

1. Introduction

The determination of the large-scale accretion geometry and orbital parameters of X-ray binaries is a fundametally important problem. Various techniques can be employed to examine the geometry of these systems, such as photometry, spectroscopy, imaging, and timing. A special place in this list belongs to polarimetry, which is known to be most sensitive to changes in geometry. The geometrical properties can be determined by tracking the changes in polarization degree (PD) and polarization angle (PA) as a function of the orbital phase. The stochastic variability on timescales comparable to the orbital period may significantly alter the average polarization profile. Dense coverage of an entire orbital cycle is needed to reliably determine the accretion geometry, shape, and orientation of the binary components.

The orbital parameters in binary systems are conventionally studied using optical and infrared polarimetry. For the low-mass X-ray binaries in outburst, emission in these wavelengths can be composed of several components: an (irradiated) accretion disk, a wind, a jet, and a hot accretion flow (Poutanen & Veledina 2014; Uttley & Casella 2014). Optical polarimetry has been used as a fine tool for distinguishing between them (Veledina et al. 2019; Kosenkov et al. 2020). In (near-)quiescence, optical polarimetry has helped constrain the role of the non-stellar components in total spectra (Kravtsov et al. 2022) and to determine the misalignment of the black hole (BH) and orbital spins (Poutanen et al. 2022). For high-mass X-ray binaries, emission in the infrared, optical, and ultraviolet bands is completely dominated by the donor star. The emission can be scattered by different large-scale components in the binary such as the accretion stream, disk, outflow or jet. The polarization signal in this case can reveal the location, orientation, and physical properties of the scattering component (Jones et al. 1994).

With the launch of the Imaging X-ray Polarimetry Explorer (IXPE; Weisskopf et al. 2022), the polarimetric field gained a second wind. It became possible to directly link the orientation of the large-scale binary components that are probed by the optical and infrared wavelengths to the innermost accretion geometry with the help of X-ray polarimetry. The prototypical BH X-ray binary Cygnus X-1 became the first target of these studies (Krawczynski et al. 2022). The week-long IXPE exposure has been accompanied by the global multiwavelength campaigns, allowing it to cover a large fraction of its 5.6 d orbital period.

Cyg X-1 is the first discovered BH X-ray binary and a well-studied system (Bowyer et al. 1965). It is a persistent source and a high-mass binary hosting a supergiant ∼40 M_⊙ donor star of spectral type O in a nearly circular orbit (eccentricity e ∼ 0.02) with the most massive Galactic BH M_BH = 21.2 ± 2.1 M_⊙ known to date (Miller-Jones et al. 2021). The donor is close to filling its Roche lobe, and the compact object accretes the matter through the focused stellar wind (Gies & Bolton 1986a). Accretion proceeds through the disk, whose emission is often seen in the X-rays (Gierliński et al. 1997, 1999; Zdziarski & Gierliński 2004), and a fraction of matter leaves the system in the form of the jet (Stirling et al. 2001; Fender et al. 2006; Miller-Jones et al. 2021).

Optical radiation is dominated by the light of the companion star and shows pronounced variations at an orbital period that is caused by the asymmetric shape of the donor (Kemp et al. 1983; Gies & Bolton 1986b; Brocksopp et al. 1999b; Orosz et al. 2011). Optical polarization measurements also show pronounced orbital variability (Nolt et al. 1975). The observed double-peak sinusoidal variations of the PD and PA are consistent with the scenario in which polarization arises from Thomson scattering of the donor star radiation by optically thin matter located within the binary (Brown et al. 1978; Milgrom 1978). This pattern is typical for binary systems and was observed in a variety of sources from classical (Piirola 1980; Berdyugin & Harries 1999) to gamma-ray binaries (Kravtsov et al. 2020). The synchronization with the orbital phase indicates that the source of the polarization is connected to the orbital motion of the BH around the companion star.

Polarization may originate from the accretion stream, its impact point on the accretion disk or the disk matter itself, or it might be related to the outflow/jet. Which component causes variations in the polarization in Cyg X-1 is unknown. The shapes of the PD and PA profiles have been used in several works to constrain orbital parameters such as inclination and eccentricity (Kemp et al. 1978, 1983; Karitskaya 1981; Dolan & Tapia 1989; Nagae et al. 2009).

The long orbital period of Cyg X-1 became an obstacle to tracking the polarimetric variations over a large fraction of a single cycle, and the average profile was obtained by including data from many cycles. This approach may lead to a substantial deviation of the obtained mean profile from the individual cycles, however, because the system is known to show substantial superorbital variability (Priedhorsky et al. 1983; Karitskaya et al. 2001; Poutanen et al. 2008; Zdziarski et al. 2011). This means that the scattering matter gradually rearranges within the binary over the superorbital period, leading to a systematic bias in the determination of the binary inclination from the mean orbital profile. On the other hand, the data obtained during a single orbital period at one telescope are unavoidably under-sampled (Dolan & Tapia 1989; Nagae et al. 2009). We performed multi-observatory polarimetric observations of Cyg X-1 that for the first time covered up to 30% of the orbit in one cycle. This was achieved by observing with nearly identical polarimeters from telescopes that were separated by ∼140° in longitude.

In this paper, we present the results of joint analysis of the new high-precision optical polarimetric observations of Cyg X-1, historical polarimetric data obtained in 1975–1987, and optical flux measurements, which allowed us to make a new attempt to qualitatively and quantitatively constrain the geometry of Cyg X-1. The paper is organized as follows. In Sect. 2 we describe the observations of the source and of the field stars to determine the contribution of the interstellar polarization. In Sect. 3 we present the main results of the study: the variations in the polarization on different timescales. We present a model for the observed orbital and superorbital variability of the polarization in Sect. 4. Finally, we summarize our findings in Sect. 5.

2. Data acquisition and analysis

We performed high-precision optical polarimetric observations of Cyg X-1 with the broad-band BVR polarimeters DIPol-2 (Piirola et al. 2014) and DIPol-UF (Piirola et al. 2021). DIPol-2 is mounted on the remotely controlled 60 cm Tohoku telescope (T60) at Haleakala Observatory, Hawaii. DIPol-UF is a visitor instrument installed at the 2.56 m Nordic Optical Telescope (NOT) at the Observatorio del Roque de los Muchachos (ORM) on La Palma, Spain. DIPol-2 and DIPol-UF are high-precision double-image CCD polarimeters, capable of measuring polarization simultaneously in three optical (BVR) bands. The polarization of the sky (even if it changes during observations) is optically eliminated by the design of the instruments. The instrumental polarization of both instruments is low (< 10⁻⁴) and is well calibrated by observing 15–20 unpolarized standard stars. The zero-point of the PA was determined by observing the highly polarized standards HD 204827 and HD 161056. Each measurement of the Stokes parameters (q_obs, u_obs) took about 20 s, and more than 200 individual measurements were obtained during the average observing night. A more detailed description of the methods and calibrations can be found in Piirola et al. (2020) and Kravtsov et al. (2022).

Cyg X-1 was observed for 27 nights between 2022 March 28 and July 28 at the T60 and for 6 nights, 2022 May 15–21, at the NOT (see Table 1). Taking advantage of the ∼140° difference in the longitude of the observatories, we covered 30% of the orbit of Cyg X-1 simultaneously with IXPE observations (red stripe in Fig. 1). We refer to this dataset as season 2 (or just S2) hereafter. The typical errors on the 30-min-averaged measurement of the Stokes parameters are σ_p ≈ 0.004% for the NOT and ≈0.01% for the T60 data.

Fig. 1. Observed normalized Stokes parameters for Cyg X-1 in the B band. The vertical dotted black lines limit the observational seasons S1 and S3. The two vertical dashed red lines show the start and end of the IXPE campaign on 2022 May 15–21 (season 2). The vertical solid purple lines show the start and end of the TESS observations.

We additionally used the historical observations carried out at the Pine Mountain Observatory (PMO), United States, in 1975–1987 (see Table 1 and Fig. 2). The reduced data were presented in parts in several papers (Kemp et al. 1978, 1979, 1983), but have never been published as a raw dataset. With the kind permission of the Pine Mountain Observatory staff, we are publishing these data (hereafter PMO data) in the public domain for the first time¹.

Fig. 2. Long-term variations in the normalized Stokes parameters q and u of Cyg X-1, measured in 1975–1987. The vertical dashed lines separate 11 observing seasons, which is roughly equal to one year of observations. The red crosses show season-averaged values.

The description of the observational techniques can be found in Kemp & Barbour (1981). Each PMO observation is a nightly average value with a typical integration time of several hours. We also used previously unpublished polarimetric observations of Cyg X-1, carried out in 2002 April–May with the TurPol polarimeter (Piirola 1973, 1988), installed on NOT, and with the 60 cm KVA telescope, ORM, La Palma.

In addition to the polarimetric data, we retrieved and analyzed the publicly available² Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014) optical light curves of Cyg X-1. We used calibrated data with a 2-min cadence (PDCSAP_FLUX), obtained in sectors 54 and 55 (2022 July – September; see Fig. 3).

Fig. 3. TESS optical light curve of Cyg X-1. The vertical dashed lines separate the consecutive orbital periods T1–T11.

3. Results

3.1. Average intrinsic polarization and its secular changes

The observed polarization of Cyg X-1 is the sum of the intrinsic polarization of the source and the interstellar (IS) polarization component that arises from the dichroism of the dust grains located between the observer and the target. The IS polarization was estimated and subtracted from the observed polarization. To find a reliable estimate of IS polarization in the source direction, we observed a sample of six field stars (see Fig. 4) that are close to Cyg X-1, as indicated by their Gaia parallaxes (see Fig. 5). We considered the wavelength dependence of the observed polarization (to exclude stars with intrinsic polarization) and took both the angular separation and proximity to the target into account. We conclude that the polarization of star Ref 2 from our sample can serve as the IS polarization estimate. Hereafter, we denote the normalized Stokes parameters of Ref 2 as (q_IS, u_IS) and subtract them from the observed values of the target (q_obs, u_obs) to obtain the Stokes parameters of the intrinsic polarization (q_int, u_int). These are translated into the intrinsic PD P and angle θ,

$$\begin{array}{c}\hfill P=\sqrt{q_{\text{int}}^2+u_{\text{int}}^2},\theta =\frac{1}{2}\arctan 2(u_{\text{int}},q_{\text{int}}).\end{array}$$(1)

Fig. 4. Finding chart and polarization properties of Cyg X-1 and the field stars. Top panel: polarization map for Cyg X-1 and the field stars. The length of the bars corresponds to the PD, and the direction corresponds to the PA (measured from north to east). Bottom panel: observed normalized Stokes parameters q and u for Cyg X-1 (stars) and field stars (circles). The blue, green, and red points with 1σ error bars correspond to the B, V, and R filters, respectively.

Fig. 5. Dependence of the PD (top panel) and PA (bottom panel) on the Gaia parallax for Cyg X-1 (red circle) and the field stars (black circles) in the R band. The horizontal error bars correspond to the errors on the Gaia parallaxes. The errors on the PD and PA are smaller than the symbol size.

The uncertainty on the PD is equal to the uncertainty of the individual Stokes parameters, and uncertainty on the PA in radians was estimated as σ_θ = σ_p/(2P) (Serkowski 1962; Kosenkov et al. 2017). The average observed and intrinsic polarization of Cyg X-1, as well as the interstellar polarization estimates, are listed in Table 2.

A large fraction of the observed polarization, about 4% out of a total 4.5–5%, has an IS origin (Kemp et al. 1979; Nagae et al. 2009). Subtracting the interstellar component from the observed polarization, we find an intrinsic PD of 0.8%±0.2% with a PA of 155° ±5° (or equivalently, −25°; see Table 2). This value is comparable to the characteristic optical PDs in other accreting BH X-ray binaries in outburst (Kosenkov et al. 2017; Veledina et al. 2019). The average intrinsic optical PA matches that measured in the X-rays within the errors (${\theta }_{\text{X}}=-22\stackrel{°}{.}2\pm 2\stackrel{°}{.}2$, Krawczynski et al. 2022).

The uniquely long history of polarimetric studies of Cyg X-1 allowed us to track the long-term evolution of the average PD almost 50 years back. We split the PMO V-band observations into 11 bins (PMO1–PMO11; see Fig. 2), each about a year long, and calculated the average values of the observed Stokes parameters within each bin. We plot them in the (q, u) diagram in Fig. 6 (colored crosses) along with our NOT+T60 2022 average polarization (red circle), KVA+NOT 2002 data (empty circle), and other published data (orange square and green diamond; Dolan & Tapia 1989; Nagae et al. 2009). We show the estimated value of the IS polarization with a blue square.

Fig. 6. Normalized Stokes parameters (qobs, uobs) of Cyg X-1 in the V band. The red circle with the error bar is the average polarization of Cyg X-1 in 2022. The crosses of different colors (from cool to warm) with error bars correspond to the average polarization in each of the 11 seasons of PMO observations in 1975–1987. The dashed blue and solid red lines show the directions of the average intrinsic polarization vector of Cyg X-1 for PMO and our observations, respectively. The other symbols correspond to the data obtained in other epochs, as described in the inset.

The blue and red lines connect the PMO and NOT+T60 2022 data with the IS estimate, respectively. The length and direction (from the IS estimate toward the data points) of these lines correspond to the vectors of the average intrinsic polarization for different epochs. The vector directions match with a high accuracy (Δθ_int < 1°). This supports our choice of the reference star Ref 2 as an estimate of interstellar polarization because the alignment of the intrinsic polarization vectors is unlikely to be accidental. We note that the other historical values shown in Fig. 6 are substantially scattered in the q − u plane, despite their small error bars. This may be caused by the orbital variations: At least several different orbital periods must be averaged to obtain a robust estimate of the average polarization.

The average intrinsic polarization for our NOT+T60 2022 data differs significantly from the PMO data (|ΔP_int|≈0.4%), indicating secular changes in the PD while preserving a constant PA. The decrease in the intrinsic PD may be caused by the decrease in the scattered flux, tailored to the secular changes of the wind density, changes in the accretion disk size, and/or its spatial orientation. Figure 2 shows that the one-year polarization averages change slightly from one season to the next. The different spread of all data points shows that the amplitude of the variability also varies on yearly timescales.

3.2. Short-term variability of the orbital profiles

The significant variability in the Stokes (q_obs, u_obs) parameters with an amplitude of about 0.1%–0.2% is clearly visible in our 2022 observations (see Fig. 1). We performed a timing analysis of our BVR polarimetric data that revealed that the main variation period of the Stokes parameters has not changed since the 1970s and is equal to half of the orbital period within the errors. To study the possible changes in the average orbital profiles of the polarization over decades, we folded our data and the PMO polarization data with the orbital phase, adopting the period P_orb from the photometric ephemeris (Brocksopp et al. 1999a). To suppress the stochastic and instrumental noise, we split the data into 18 orbital phase bins, in which we calculated weighted average values of the Stokes parameters that were subsequently used to obtain the PD P and PA θ. The comparison of our light curve and the PMO polarization light curves is shown in Fig. 7. Except for the systematic offset between our and PMO data, the nature of which is discussed above, the shapes of the average PD and PA orbital profiles agree exceptionally well with each other.

Fig. 7. Orbital profiles of the observed and intrinsic polarization of Cyg X-1. Left panels: observed PD, P, (top panel) and PA, θ, (bottom panel) of Cyg X-1 in the V band, folded with the orbital period. The blue circles correspond to the data obtained in 2022. The red crosses correspond to the PMO data obtained in 1975–1987. Each point corresponds to the average value, calculated within a phase bin of width Δϕ = 1/18. The typical 1σ uncertainty is smaller than the symbol size. Right panels: same as the left panel, but showing the intrinsic polarization of Cyg X-1.

To determine how the shape of optical light curves changes from one period to another, we split TESS photometric data into 11 consecutive orbital periods T1–T11 (see Fig. 3). The shape of each individual profile is far from the double-sine wave that is expected in the case of ellipsoidal variations that are caused by the rotation of a tidally distorted star around the center of mass. A short-period variability is superimposed on the main double-sine curve, which leads to changes in the amplitudes and phases of the main maxima/minima. In Fig. 8 we show the orbital profiles of the intrinsic polarization in R band together with TESS photometric profiles. The optical polarization and optical flux both show double-sinusoidal orbital variations with minima in the conjunctions (phases 0 and 0.5) and maxima in the quadratures (phases 0.25 and 0.75).

Fig. 8. Orbital profiles of the polarization and flux of Cyg X-1 in the optical band. Top panel: intrinsic polarization of Cyg X-1 in the R filter, folded with the orbital period (seasons 1–3 are plotted). Each circle with the 1σ error bar shows the average polarization, calculated within a 30-min bin. Bottom panel: TESS magnitude of Cyg X-1, folded with the orbital period. Different colors (from cold to warm) correspond to different orbital periods T1–T11. The solid black line shows the average orbital profile.

Photometric variations arise from the nonspherical shape of the tidally distorted companion, whose visible area (and hence the flux) is largest around the quadratures. If the scattering of the donor star emission occurs in a region that is connected to the compact object, we expect the scattering angle to reach a maximum of 90° in the same phases, corresponding to the maximum PD in the case of Thomson scattering. In the conjunctions, the visible area of the supergiant approaches its minima, resulting in the lowest flux, while at the same time, the scattering angle reaches minimum or maximum, leading to a smaller PD. The short-term changes in both the flux and polarization, which are superimposed on the periodic variations, can originate from one or more mechanisms: the pulsations of the main star, spots on its surface, inhomogeneities of the wind, or eclipses of the bright parts of the disk by the infalling matter. We note that despite the correlation of flux and polarization, the polarization variability cannot be explained by variations in the unpolarized flux alone. For the intrinsic PD P = F_pol/F_tot ∼ 0.01, the unpolarized flux variations of about ΔF_tot ∼ 5 × 10⁻³F_tot give a negligibly small polarization variability:

$$\begin{array}{c}\hfill \mathrm{\Delta }P\approx P\frac{\mathrm{\Delta }{F}_{\mathrm{tot}}}{F_{\mathrm{tot}}}\sim 5\times {10}^{-5},\end{array}$$

while the observed one is at least factor 20 of larger.

With our exceptionally dense orbital coverage, we can compare the profiles of a single cycle with the average one, as given by the PMO data. In Fig. 9 we show the profile obtained from the 8-day-long monitoring of Cyg X-1 during season 2, overlaid on the average profile of the polarization variability in the V band. The figure shows that although the overall shapes of the polarization variability curves are roughly consistent with the patterns of 1975–1987, the amplitude of the variations is substantially higher in our season 2 data, where the harmonic content is also richer. These facts support the statement of Dolan & Tapia (1989) about the existence of nonorbital polarization variability and the importance to account for it when extracting orbital parameters. In the following sections, we describe the modeling of the polarization variability curves with different analytical models and discuss how the short-term variability affects the results.

Fig. 9. Variability in the intrinsic PD (upper panel) and PA (bottom panel) of Cyg X-1 in V band, measured in May 2022 (blue circles). Each circle corresponds to the average value, calculated within a 30-min bin. The 1σ errors are smaller than the symbol size. The red crosses correspond to the average binned polarization, measured by Kemp with colleagues during 1975–1987, with a constant shift in PD by ΔP = −0.4%.

3.3. Superorbital evolution of polarization profiles

In addition to the short-term variability, indications of long-term changes in the Cyg X-1 polarization profiles have been reported by several authors. Kemp et al. (1983) suggested long-term optical polarization variations at the superorbital period of 294 d, discovered in the X-rays (Priedhorsky et al. 1983). The authors discussed several models that could explain the variations, including the precession of the accretion disk and the obscuration of the scattering medium. Comparing the average optical polarization obtained between 1975 and 2006, Nagae et al. (2009) found secular variations in the average polarization component of Cyg X-1.

We find signatures of the long-term variability in our 2022 polarimetric data. Figure 10 shows the change in the polarization profiles with superorbital phase, separated roughly by a month (see Fig. 1). The changes in the average values of Stokes parameters along with the changes in the amplitude and profiles of the orbital variations are significant. In Fig. 11 we show the superorbital profile of the V-band polarization of Cyg X-1. The average values of the PD for seasons S1 – S3 (empty red circles), folded with the superorbital period (P_sup = 294 d, JD₀ = 2 440 000), are consistent with the same part of the superorbital profile as observed in PMO data.

Fig. 10. Orbital variations in the observed normalized Stokes parameters q (left) and u (right) of Cyg X-1 in the B band for different seasons (seasons 1, 2, and 3 from top to bottom). The horizontal dashed lines show the weighted average values of the corresponding parameters.

Fig. 11. Intrinsic PD of Cyg X-1 in the V band for PMO (blue circles) and S1–S3 (empty red circles) data, folded with the superorbital period. The red points are shifted by a constant ΔP = −0.4% in the vertical direction to take secular changes in the PD into account.

4. Modeling

To explain the behavior of the polarization at different timescales, we considered several possibilities for the geometry of scattering matter. We started with the generic model for polarization production in binary systems, in which polarization arises from the Thomson scattering of the companion star radiation by a cloud of an optically thin matter near the compact object (Brown et al. 1978). The key assumption of the model is the corotation of the scattering material with the secondary, in our case, with the compact object. The PD in this case peaks at orbital phases at which the scattering angle is 90°. For a circular orbit, the Stokes parameters of the linear polarization vary as a sine-like wave at twice the orbital frequency. In the case of eccentric orbit and/or asymmetry of the distribution of the light-scattering material about the orbital plane, the profiles become skewed and can be described by adding the first harmonic of the orbital period. Alternatively, the appearance of the first harmonic can be related to the presence of an optically thick scattering material. Below, we study the harmonic content of the polarization profiles and consider different possibilities for the geometry of the scattering matter.

4.1. Fourier method

The polarization profiles corresponding to the case of the optically thin corotating scatterer in a circular orbit can be decomposed into Fourier series of the orbital longitude λ = 2πϕ (where ϕ is the orbital phase),

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

We employed Bayesian inference implemented as the Markov chain Monte Carlo (MCMC; Goodman & Weare 2010) ensemble sampler in the emcee package (Foreman-Mackey et al. 2013) in PYTHON to fit the orbital profiles of the Stokes parameters observed during season S2 with Eq. (2). The best-fit curves are shown in Fig. 12. Following the approach described in Drissen et al. (1986) and Kravtsov et al. (2020), we used the obtained Fourier coefficients to derive the inclination i of the binary,

$$\begin{array}{c}\hfill {\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},\end{array}$$(3)

Fig. 12. Variability in the observed Stokes parameters of Cyg X-1, obtained during season 2 in B, V, and R bands (top, middle, and bottom panels, respectively). Each circle corresponds to the average value, calculated within a 30-min bin. The 1σ errors are smaller than the symbol size. The solid black lines correspond to the best fit with the Fourier series given by Eq. (2). The dashed red lines in the middle panels show the best fit of PMO historical V -band data with the same model, shifted vertically to overlap our data.

and the position angle Ω of the orbital axis on the sky,

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

where

$$\begin{array}{c}\hfill \begin{array}{cc}& 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)

By fitting the orbital polarization profiles obtained in S2 with Eq. (2), we obtained formal values of the inclination i = 125° ±5° (i > 90° indicates the clockwise apparent motion of the compact object on the sky) and the position angle Ω = 129° ±5° of the orbital axis on the sky. However, the formal errors on the estimated orbital parameters obtained from the error propagation are underestimated and hence do not correspond to their actual confidence intervals, which are determined primarily by the internal properties of the model (2) and the amplitude of the stochastic variability in the data. The inclination estimates corresponding to the best-fit Fourier coefficients are always biased toward higher values (Aspin et al. 1981; Simmons et al. 1982; Wolinski & Dolan 1994). The confidence intervals on the orbital parameters for different signal-to-noise ratios can be obtained using Fig. 4 of Wolinski & Dolan (1994): 1σ and 2σ confidence intervals were calculated for four levels of data quality given by γ = 0.5N(A/σ_p)², where σ_p is the standard deviation of noise in the data, A is the amplitude of the polarimetric variability, and N is the number of observations. Our value of γ_obs = 0.5 × 100 × (6.7)² ≈ 2200 lies between their grid points (γ = 120 000 and γ = 300). To calculate the confidence intervals on the inclination that we obtained for our S2 data, we therefore performed our own Monte Carlo simulations following the procedure described in Wolinski & Dolan (1994): We modeled the Stokes parameters for different values of i ranging from 90° to 180° using the standard Brown et al. (1978) model (Eq. (6) in that paper). Then, we simulated the Gaussian noise in q and u by adding the fluctuations of the variance ${\sigma }_{\text{p}}^2=0.5N{A}^2/{\gamma }_{\text{obs}}$. The Fourier model (Eq. (2)) was then fit to the simulated data using the MCMC approach, and the inclination i′ was calculated using Eq. (3). In Fig. 13 we show the inclination estimates i′ (black points) with the 1σ and 2σ confidence intervals (solid blue and red lines) as a function of the true input inclination i.

Fig. 13. Estimated 1σ and 2σ confidence intervals on the true inclination i for given estimate i′ (solid blue and red lines). The vertical dashed black line corresponds to the best-fit inclination of Cyg X-1. The vertical dotted blue lines correspond to the 1σ error on the best-fit inclination.

The inclination i ≈ 125° that we derived from the best-fit Fourier coefficients of Cyg X-1 in V and R bands (shown as the dashed vertical black line in Fig. 13) is close to the so-called critical angle $i_{\text{crit}}^{\prime }$. Above this angle, the 1σ confidence interval on the orbital inclination extends to i = 180° (Wolinski & Dolan 1994). This means that using high-precision polarimetry, we can only place a lower limit of 180° > i > 120° on the inclination value of the Cyg X-1 orbit. We note that previous polarimetrically derived inclination values (Kemp et al. 1978; Dolan & Tapia 1989; Nagae et al. 2009) are most likely overestimated because they were obtained by modeling the data with larger error bars, for which the critical angle is expected to be smaller than $i_{\text{crit}}^{\prime }$ ∼ 130°. Our lower limit on the inclination i > 120° and the clockwise direction of the orbital motion on the sky are consistent with the value i = 153° ±1° from Miller-Jones et al. (2021).

In contrast to the inclination, the value Ω obtained from Eq. (4) is an unbiased estimate of the true position angle of the projection of the orbital axis on the sky. We used the same MCMC approach as for the inclination to calculate the confidence interval on this angle. Our value Ω = 129° ±10° (or Ω = 129° −180° = − 51° ±10° because of the ±180° ambiguity) is consistent within 3σ with those that were determined by the direction of the intrinsic polarization θ_int ≈ −25° and the position angle of the jet on the sky Ω_jet ≈ −26° (Miller-Jones et al. 2021).

We emphasize that this low accuracy is not a result of the polarization measurement errors (which are smaller than 0.01% for the whole set of our data). Figure 12 shows a remarkable intrinsic scatter of the S2 data points around the fit curves that is especially noticeable for the Stokes q parameter. This aperiodic noise is explained by an additional suborbital variability component that appears on timescales shorter than one orbital cycle. Thus, the key assumption on corotation of the light-scattering material over (at least) a few consecutive orbital cycles does not hold for the Cyg X-1 binary system. Therefore, the traditional Fourier fit up to second harmonics made on polarization data cannot provide meaningful estimates of the orbital inclination, regardless of data quality, quantity, and sampling frequency.

4.2. Eccentric model

While for a circular orbit, theory predicts a smooth change in polarization with the dominant second harmonic of the orbital period, the eccentricity of the orbit shifts all the changes in the polarization toward the periastron. The polarization depends on the scattering angle, which changes according to the orbital motion of the scattering cloud. In the case of an eccentric orbit, this angle changes with different rates in different parts of the orbit, resulting in unequal distances between consecutive maxima/minima of the orbital Stokes parameters curves (this effect was observed for the binary with e ≈ 0.4 in Berdyugin & Tarasov 1998). Therefore, the orbital curves of Stokes parameters can be used for an independent estimation of the orbital eccentricity. We adopted the Thomson scattering model from the appendix of Kravtsov et al. (2020) to describe the orbital changes in the polarization of Cyg X-1. By fitting this model to the V-band season 2 data, we were able to place 3σ upper limit on the eccentricity of Cyg X-1 orbit. The eccecntricity is e < 0.08.

Figure 14 shows the (q, u) plane of the average Stokes parameters of Cyg X-1 obtained during season 2, together with the best-fits with the Fourier series (Eq. (2)) and the model of the Thomson scattering by a cloud on an eccentric orbit (Kravtsov et al. 2020). The latter model (which is a special case of the general BME model for symmetrically distributed matter in the orbital plane) cannot reproduce the pretzel-like shape of the trace left by orbital variations of Cyg X-1 on the (q, u)-plane: The additional source of asymmetry is needed. To explain a similar pattern, Kemp et al. (1978) proposed a model in which the scattering region is eclipsed by the secondary body for half the orbit. This model requires a high (90° > i > 130°) orbital inclination, which contradicts the latest results (including this article).

Fig. 14. Average Stokes parameters of Cyg X-1 obtained during season 2 in the R band. The color coding and size of the circle correspond to the orbital phase and 1σ uncertainty, respectively. The solid blue curve corresponds to the best fit with the Fourier series given by Eq. (2). The dashed orange line corresponds to the to the best-fit model of a scattering cloud on an eccentric orbit from the appendix of Kravtsov et al. (2020).

4.3. Polarization by Thomson scattering off a precessing accretion disk

In this section, we present a model of polarization from Thomson scattering by a tilted precessing accretion disk, which can naturally explain the asymmetric pattern of the polarization variability observed in Cyg X-1 without requiring a highly inclined or eccentric orbit. We considered the following geometry: The orbit with an eccentricity e is inclined by an angle i to the line of sight $\widehat{\mathit{o}}$ (Fig. 15). The accretion disk, surrounded by a cloud of electrons, rotates around the optical companion together with the compact object. The disk axis ${\widehat{\mathit{n}}}_{\mathrm{d}}$ is inclined by an angle β to the orbital axis $\widehat{\mathit{n}}$ (see Fig. 16). The axis of the disk can precess about the orbital axis with the period T_sup.

Fig. 15. Geometry of the system. Optical emission of the companion star (blue circle at the origin) is scattered by the tilted accretion disk (orange) and the optically thin material (indicated by the black circle in the center of the disk).

Fig. 16. Geometry of the precessing disk. The disk axis $\widehat{{\mathit{n}}_{\mathbf{d}}}$ is inclined by the angle β with respect to the orbital axis $\widehat{\mathit{n}}$ and rotates around it.

To describe the orbital motion, we introduced the coordinate system ($\widehat{\mathit{x}}$, $\widehat{\mathit{y}}$, $\widehat{\mathit{z}}$), in which the $\widehat{\mathit{z}}$-axis is directed along the orbital axis $\widehat{\mathit{n}}$, the vector $\widehat{\mathit{x}}$ lies in the orbital plane and its projection on the sky is directed to the south, and the vector $\widehat{\mathit{y}}$ forms the right-handed basis. In this basis, $\widehat{\mathit{o}}=(sini,0,cosi)$, $\widehat{\mathit{n}}=(0,0,1)$, and ${\widehat{\mathit{n}}}_{\mathrm{d}}=(sin\beta cos\gamma ,sin\beta sin\gamma ,cos\beta )$. The angle γ is the azimuth of the projection of the disk axis onto the orbital plane measured from $\widehat{\mathit{x}}$ to $\widehat{\mathit{y}}$. To describe the polarization, we used the polarization basis (${\widehat{\mathit{e}}}_1$, ${\widehat{\mathit{e}}}_2$), in which the vector ${\widehat{\mathit{e}}}_1=(-cosi,0,sini)$ lies along the projection of the vector $\widehat{\mathit{n}}$ on the plane of the sky, and ${\widehat{\mathit{e}}}_2=(0,-1,0)$ is perpendicular to ${\widehat{\mathit{e}}}_1$ and lies in the plane of the orbit.

The distance between the compact object and the optical companion varies with the orbital longitude λ, measured from $\widehat{\mathit{x}}$ to $\widehat{\mathit{y}}$ . It can be expressed as

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

where a is the semimajor axis of the orbit, and λ_p is the longitude of the periastron. The unit vector pointing toward the compact object is

$$\begin{array}{c}\hfill \widehat{\mathit{r}}=(cos\lambda ,sin\lambda ,0),\end{array}$$(7)

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}$$(8)

The observed flux F_tot = F_* + F_sc is the sum of the flux produced by the optical companion F_* and the scattered flux F_sc. We assumed that the latter is produced by Thomson scattering (in an optically thin regime) of stellar radiation by the accretion disk and the surrounding cloud of electrons. In this case, 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 )=\frac{3}{8}(1+{\mu }^2)L_{\ast }(f_{\mathrm{sc}}^{\mathrm{cloud}}+f_{\mathrm{sc}}^{\mathrm{disk}}),\end{array}$$(9)

where l(μ) = 3(1 + μ²)/8 is the Thomson scattering indicatrix, and $f_{\mathrm{sc}}^{\mathrm{cloud}}$ and $f_{\mathrm{sc}}^{\mathrm{disk}}$ are fractions of radiation scattered by the cloud and the disk, respectively. In both cases, this fraction is proportional to the total number of free electrons N_e in a cloud/disk and drops with the distance as 1/r²(λ),

$$\begin{array}{cc}& f_{\mathrm{sc}}^{\mathrm{cloud}}=\frac{N_{\mathrm{e}}^{\mathrm{cloud}}{\sigma }_{\mathrm{T}}}{4\pi r^2(\lambda )},\hfill \end{array}$$(10)

$$\begin{array}{cc}& f_{\mathrm{sc}}^{\mathrm{disk}}=\frac{N_{\mathrm{e}}^{\mathrm{disk}}{\sigma }_{\mathrm{T}}}{4\pi r^2(\lambda )}cos\mathrm{\Psi },\hfill \end{array}$$(11)

where $cos\mathrm{\Psi }=-{\widehat{\mathit{n}}}_{\mathrm{d}}·\widehat{\mathit{r}}=-sin\beta cos\gamma cos\lambda -sin\beta sin\gamma sin\lambda =-sin\beta cos(\gamma -\lambda )$. The cos Ψ term is proportional to the effective area of the disk intercepting the stellar radiation, which depends on the position of the disk in the orbit λ and the orientation of its axis ${\widehat{\mathit{n}}}_{\mathrm{d}}$, defined by two angles: the inclination β of the disk, and its azimuth γ. The latter angle can change with time due to precession as γ = ±2πφ_sup + γ₀, where φ_sup is the precession phase, and γ₀ is the angle γ at zero phase. The sign determines the direction (counterclockwise or clockwise) of the precession.

We scaled $f_{\mathrm{sc}}^{\mathrm{cloud}}$ and $f_{\mathrm{sc}}^{\mathrm{disk}}$ to the typical values $f_0^{\text{cloud}}$ and $f_0^{\text{disk}}$ as

$$\begin{array}{cc}& f_{\mathrm{sc}}^{\mathrm{cloud}}=f_0^{\mathrm{cloud}}{\left[\frac{a(1-e^2)}{r(\lambda )}\right]}^2=f_0^{\mathrm{cloud}}{[1+ecos(\lambda -{\lambda }_{\mathrm{p}})]}^2,\hfill \end{array}$$(12)

$$\begin{array}{cc}& f_{\mathrm{sc}}^{\mathrm{disk}}=f_0^{\mathrm{disk}}{\left[\frac{a(1-e^2)}{r(\lambda )}\right]}^{2}cos\mathrm{\Psi }=f_0^{\mathrm{disk}}{[1+ecos(\lambda -{\lambda }_{\mathrm{p}})]}^{2}cos\mathrm{\Psi }.\hfill \end{array}$$(13)

The PD of scattered radiation in Thomson regime can be expressed in terms of the scattering angle Θ as

$$\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}$$(14)

The observer measures the PD P = F_scP_sc/F_tot of the total flux F_tot, most of which is unpolarized and produced by the optical companion star. The polarized flux of the scattered radiation 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}$$(15)

Therefore, the total PD

$$\begin{array}{c}\hfill P=\frac{F_{\mathrm{sc}}{P}_{\mathrm{sc}}}{F_{\mathrm{tot}}}\approx \frac{3}{8}[f_{\mathrm{sc}}^{\mathrm{cloud}}+f_{\mathrm{sc}}^{\mathrm{disk}}](1-{\mu }^2),\end{array}$$(16)

where we assumed that F_sc ≪ F_* and substituted $f_{\mathrm{sc}}=f_{\mathrm{sc}}^{\mathrm{cloud}}+f_{\mathrm{sc}}^{\mathrm{disk}}$.

The normalized Stokes parameters of linear polarization are defined as q = Pcos(2χ) and u = Psin(2χ), where χ is the position angle of the polarization pseudo-vector $\widehat{\mathit{p}}$ in the polarization basis (${\widehat{\mathit{e}}}_1$, ${\widehat{\mathit{e}}}_2$),

$$\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 }}(-cosisin\lambda ,cosicos\lambda ,sinisin\lambda ),\end{array}$$(17)

where $sin\mathrm{\Theta }=\sqrt{1-{\mu }^2}=\sqrt{1-{sin}^{2}i{cos}^2\lambda }$. Thus, the expressions for the Stokes parameters can be written as

$$\begin{array}{c}\hfill \begin{array}{cc}& q=\frac{3}{8}[f_{\mathrm{sc}}^{\mathrm{cloud}}+f_{\mathrm{sc}}^{\mathrm{disk}}](1-{\mu }^2)cos(2\chi ),\hfill \\ \hfill & u=\frac{3}{8}[f_{\mathrm{sc}}^{\mathrm{cloud}}+f_{\mathrm{sc}}^{\mathrm{disk}}](1-{\mu }^2)sin(2\chi ),\hfill \end{array}\end{array}$$(18)

where the polarization angle χ is defined by the expressions

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

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

The explicit expressions for cos(2χ) and sin(2χ) are

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

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

Combining Eqs. (18) with (21) and (22), we obtain

$$\begin{array}{c}\hfill \begin{array}{cc}& q=\frac{3}{16}[{sin}^{2}i-(1+{cos}^{2}i)cos2\lambda ][f_{\mathrm{sc}}^{\mathrm{cloud}}+f_{\mathrm{sc}}^{\mathrm{disk}}],\hfill \\ \hfill & u=-\frac{3}{8}cosisin2\lambda [f_{\mathrm{sc}}^{\mathrm{cloud}}+f_{\mathrm{sc}}^{\mathrm{disk}}].\hfill \end{array}\end{array}$$(23)

The $f_{\mathrm{sc}}^{\mathrm{disk}}$ term depends on cos Ψ, reflecting the difference in the amount of scattered radiation for different orientations of the disk axis ${\widehat{\mathit{n}}}_{\mathrm{d}}$ relative to the source of the light. In our model, the disk is not transparent: it has two sides (top and bottom), only one of which is illuminated at any given time. The top of the disk is illuminated when cos Ψ > 0, the bottom of the disk is bright when cos Ψ < 0. The top of the disk is visible for the observer when $cos\mathrm{\Sigma }={\widehat{\mathit{n}}}_{\mathrm{d}}·\widehat{\mathit{o}}>0,$ and the bottom of the disk is visible when cos Σ < 0. Therefore, the disk is illuminated and visible only when the product cos Ψ cos Σ is positive, or $f_{\mathrm{sc}}^{\mathrm{disk}}=0$ when cos Ψ cos Σ ≤ 0.

To compare the calculations with the observed Stokes parameters (q_obs, u_obs), the orientation of the orbit on the sky must be taken into account: The projection of the orbital axis on the sky makes an angle Ω with direction to the north (we note that Ω defined in this way differs by π/2 from the longitude of the ascending node that is commonly used instead). The observed Stokes parameters (q_obs, u_obs) can be obtained by rotating the vector (q, u) by an angle 2Ω,

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

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

The modeled Stokes parameters are functions of the orbital longitude λ and need to be computed as functions of the orbital phase ϕ. While for the circular (or nearly circular) orbit, λ can be calculated as λ = 2π(ϕ + ϕ_p)+λ_p, where ϕ_p is the phase of the periastron, for the eccentric orbit, we need to solve Kepler’s equation: From the true anomaly of the orbit λ − λ_p, we can find the eccentric anomaly E,

$$\begin{array}{c}\hfill tan\left(\frac{E}{2}\right)=\sqrt{\frac{1-e}{1+e}}tan\left(\frac{\lambda -{\lambda }_{\text{p}}}{2}\right),\end{array}$$(26)

and then the mean anomaly M,

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

which can be converted into the orbital phase as

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

Thus, the free parameters of the model are the inclination i, the eccentricity e, the longitude of periastron λ_p, the position angle of the projection of the orbit axis on the sky Ω, the phase of periastron ϕ_p, the inclination of the disk β and its initial position angle γ₀ in the orbital plane, the period of precession T_sup (which can be set to be infinite for the nonprecessing case), and the scattering fractions $f_0^{\text{disk}}$ and $f_0^{\text{cloud}}$. In order to fit the data, there is a need for additional constant Stokes parameters q₀ and u₀, which describe the average polarization.

We fit the described model to S2 data by adopting the orbital parameters of Cyg X-1: eccentricity e = 0.02, inclination i = 153° (Miller-Jones et al. 2021), and inclination of the disk β = 20° (Ibragimov et al. 2007). The position angle of the orbital axis was set to be Ω = −26° to match the position angle of the jet. The contributions from the disk and the cloud were assumed to be $f_0^{\text{disk}}=0.75$ and $f_0^{\text{cloud}}=0.25$. Because the S2 data cover only one full orbital cycle, the precession period T_sup was set to be much longer than the orbital period so that possible precession of the disk was not taken into account.

The solid blue and dashed red lines in Figs. 17 and 18 show the fits of the model described above to S2 data with and without scattering off the tilted accretion disk, respectively. Although the reduced χ² of the fits does not differ dramatically (χ²[disk + cloud] = 1.01 versus χ²[cloud] = 1.23), the asymmetric pretzel-like trace of the polarization on the (q, u)-plane cannot be reproduced by the scattering cloud alone. Additionally, the model with the tilted accretion disk predicts the changes in the shape of the orbital polarization profiles with the precession phase: If the superorbital variability observed from the radio to the X-rays is related to the disk precession, the pretzel will make a complete turn around its center once per superorbital period T_sup. To detect this effect, a significant part of the superorbital period must be covered with continuous high-precision optical polarimetric observations.

Fig. 17. Same as Fig. 14, but showing the model curves at the (q, u)-plane calculated with (solid blue line) and without (dashed red line) scattering by the accretion disk.

Fig. 18. Stokes parameters of Cyg X-1, obtained during season 2 (light crosses) together with the best-fit models with (solid blue line) and without (dashed red line) scattering by the tilted accretion disk.

5. Summary

We presented new high-precision polarimetric observations of the BH X-ray binary Cyg X-1. Combining them with the 12-year-long PMO observations performed in 1975–1987, we were able to study the polarization behavior at the timescales ranging from hours to decades. The interstellar polarization, which dominates the observed optical polarization P_obs ∼ 4.5%, was accurately measured and subtracted from the data, allowing us to determine the intrinsic polarization of Cyg X-1 as P_int = 0.82%±0.15% with a PA θ_int = 155° ±5°. The alignment of the X-ray and optical PA, as well as the stability of this angle during the secular PD change, indirectly support our estimate of intrinsic polarization of Cyg X-1. Around-the-clock monitoring of the polarization with two telescopes located in different hemispheres allowed us to track the evolution of the polarization within one orbital cycle with a temporal resolution that is unrivalled so far. The intrinsic polarization of Cyg X-1 shows the orbital variations with two pronounced peaks in the quadratures and two minima in conjugations, most probably produced by Thomson scattering of the companion star radiation by matter that is gravitationally bound to the black hole. The amplitudes of the two consecutive polarization minima measured within one orbital cycle differ significantly, which implies asymmetry of the scattering matter in the orbital plane. We suggest that a tilted accretion disk could be the source of this asymmetry. We find that a misalignment of β ≳ 15° can reproduce the orbital polarization variations. This is in line with the recent finding of a significant misalignment between the orbital and jet axes (Zdziarski et al. 2023). Our modeling of orbital variations in the Stokes parameters allowed us to constrain the eccentricity e < 0.08 and inclination of the orbit i > 120°.

In addition to the orbital variations, we found a significant change (ΔP_int ≈ −0.4%) in the average intrinsic PD of Cyg X-1 on a timescales of several decades while preserving the constant intrinsic PA. The decrease in the PD indicates the change in the fraction of scattered radiation, which in turn depends on the amount of scattering material and its effective scattering cross section. This may reflect secular changes in the size/shape of the accretion disk and/or changes in its spatial orientation. We note that the asymmetry of the (q, u)-plane trace of the polarization can be purely artificial. It may result from a complex superposition of periodic and nonperiodic variations and the orbital phase sampling. A long-term high-precision monitoring program with good orbital and superorbital coverage is needed to exclude this.

Analyzing high-precision TESS photometric data, we found stochastic variations of the flux on timescales shorter than the orbital period. Together with the stochastic variability found in the optical polarization, this suggests that one or several additional components are at play: pulsations of the optical companion, spots on its surface, wind clumpiness, eclipses of the bright part of the accretion disk by the infalling matter, and precession of the accretion disk.

Acknowledgments

This paper is based on observations made with the Nordic Optical Telescope, owned in collaboration by the University of Turku and Aarhus University, and operated jointly by Aarhus University, the University of Turku, and the University of Oslo, representing Denmark, Finland, and Norway, the University of Iceland and Stockholm University at the Observatorio del Roque de los Muchachos, La Palma, Spain, of the Instituto de Astrofisica de Canarias. The DIPol-2 and DIPol-UF polarimeters were built in cooperation between the University of Turku, Finland, and the Leibniz-Institut für Sonnenphysik, Germany. We are grateful to the Institute for Astronomy, University of Hawaii, for the allocated observing time. VK acknowledges support from the Vilho, Yrjö and Kalle Väisälä Foundation and Suomen Kulttuurirahasto. AAZ has been supported by the Polish National Science Center under the grant 2019/35/B/ST9/03944.

References

All Tables

All Figures

	Fig. 1. Observed normalized Stokes parameters for Cyg X-1 in the B band. The vertical dotted black lines limit the observational seasons S1 and S3. The two vertical dashed red lines show the start and end of the IXPE campaign on 2022 May 15–21 (season 2). The vertical solid purple lines show the start and end of the TESS observations.
In the text

	Fig. 2. Long-term variations in the normalized Stokes parameters q and u of Cyg X-1, measured in 1975–1987. The vertical dashed lines separate 11 observing seasons, which is roughly equal to one year of observations. The red crosses show season-averaged values.
In the text

	Fig. 3. TESS optical light curve of Cyg X-1. The vertical dashed lines separate the consecutive orbital periods T1–T11.
In the text

Fig. 4. Finding chart and polarization properties of Cyg X-1 and the field stars. Top panel: polarization map for Cyg X-1 and the field stars. The length of the bars corresponds to the PD, and the direction corresponds to the PA (measured from north to east). Bottom panel: observed normalized Stokes parameters q and u for Cyg X-1 (stars) and field stars (circles). The blue, green, and red points with 1σ error bars correspond to the B, V, and R filters, respectively.

In the text

	Fig. 5. Dependence of the PD (top panel) and PA (bottom panel) on the Gaia parallax for Cyg X-1 (red circle) and the field stars (black circles) in the R band. The horizontal error bars correspond to the errors on the Gaia parallaxes. The errors on the PD and PA are smaller than the symbol size.
In the text

Fig. 6. Normalized Stokes parameters (qobs, uobs) of Cyg X-1 in the V band. The red circle with the error bar is the average polarization of Cyg X-1 in 2022. The crosses of different colors (from cool to warm) with error bars correspond to the average polarization in each of the 11 seasons of PMO observations in 1975–1987. The dashed blue and solid red lines show the directions of the average intrinsic polarization vector of Cyg X-1 for PMO and our observations, respectively. The other symbols correspond to the data obtained in other epochs, as described in the inset.

In the text

Fig. 7. Orbital profiles of the observed and intrinsic polarization of Cyg X-1. Left panels: observed PD, P, (top panel) and PA, θ, (bottom panel) of Cyg X-1 in the V band, folded with the orbital period. The blue circles correspond to the data obtained in 2022. The red crosses correspond to the PMO data obtained in 1975–1987. Each point corresponds to the average value, calculated within a phase bin of width Δϕ = 1/18. The typical 1σ uncertainty is smaller than the symbol size. Right panels: same as the left panel, but showing the intrinsic polarization of Cyg X-1.

In the text

Fig. 8. Orbital profiles of the polarization and flux of Cyg X-1 in the optical band. Top panel: intrinsic polarization of Cyg X-1 in the R filter, folded with the orbital period (seasons 1–3 are plotted). Each circle with the 1σ error bar shows the average polarization, calculated within a 30-min bin. Bottom panel: TESS magnitude of Cyg X-1, folded with the orbital period. Different colors (from cold to warm) correspond to different orbital periods T1–T11. The solid black line shows the average orbital profile.

In the text

Fig. 9. Variability in the intrinsic PD (upper panel) and PA (bottom panel) of Cyg X-1 in V band, measured in May 2022 (blue circles). Each circle corresponds to the average value, calculated within a 30-min bin. The 1σ errors are smaller than the symbol size. The red crosses correspond to the average binned polarization, measured by Kemp with colleagues during 1975–1987, with a constant shift in PD by ΔP = −0.4%.

In the text

	Fig. 10. Orbital variations in the observed normalized Stokes parameters q (left) and u (right) of Cyg X-1 in the B band for different seasons (seasons 1, 2, and 3 from top to bottom). The horizontal dashed lines show the weighted average values of the corresponding parameters.
In the text

	Fig. 11. Intrinsic PD of Cyg X-1 in the V band for PMO (blue circles) and S1–S3 (empty red circles) data, folded with the superorbital period. The red points are shifted by a constant ΔP = −0.4% in the vertical direction to take secular changes in the PD into account.
In the text

Fig. 12. Variability in the observed Stokes parameters of Cyg X-1, obtained during season 2 in B, V, and R bands (top, middle, and bottom panels, respectively). Each circle corresponds to the average value, calculated within a 30-min bin. The 1σ errors are smaller than the symbol size. The solid black lines correspond to the best fit with the Fourier series given by Eq. (2). The dashed red lines in the middle panels show the best fit of PMO historical V -band data with the same model, shifted vertically to overlap our data.

In the text

	Fig. 13. Estimated 1σ and 2σ confidence intervals on the true inclination i for given estimate i′ (solid blue and red lines). The vertical dashed black line corresponds to the best-fit inclination of Cyg X-1. The vertical dotted blue lines correspond to the 1σ error on the best-fit inclination.
In the text

Fig. 14. Average Stokes parameters of Cyg X-1 obtained during season 2 in the R band. The color coding and size of the circle correspond to the orbital phase and 1σ uncertainty, respectively. The solid blue curve corresponds to the best fit with the Fourier series given by Eq. (2). The dashed orange line corresponds to the to the best-fit model of a scattering cloud on an eccentric orbit from the appendix of Kravtsov et al. (2020).

In the text

	Fig. 15. Geometry of the system. Optical emission of the companion star (blue circle at the origin) is scattered by the tilted accretion disk (orange) and the optically thin material (indicated by the black circle in the center of the disk).
In the text

	Fig. 16. Geometry of the precessing disk. The disk axis $\widehat{{\mathit{n}}_{\mathbf{d}}}$ is inclined by the angle β with respect to the orbital axis $\widehat{\mathit{n}}$ and rotates around it.
In the text

	Fig. 17. Same as Fig. 14, but showing the model curves at the (q, u)-plane calculated with (solid blue line) and without (dashed red line) scattering by the accretion disk.
In the text

	Fig. 18. Stokes parameters of Cyg X-1, obtained during season 2 (light crosses) together with the best-fit models with (solid blue line) and without (dashed red line) scattering by the tilted accretion disk.
In the text