Apparent motion of the circumstellar envelope of CQ Tau in scattered light

The study of spiral structures in protoplanetary disks is of great importance for understanding of processes in the disks, including planet formation. Bright spiral arms were detected in the disk of young star CQ Tau by Uyama et al. (2020) in H and L bands. The spiral arms are located inside the gap in millimeter size dust, recovered earlier using ALMA observations (Ubeira Gabellini et al. 2019). To explain the gap, Ubeira Gabellini et al. (2019) proposed the existence of a planet with the semimajor axis of 20 AU. We obtained multiepoch observations of a spiral feature in the circumstellar envelope of CQ Tau in Ic band using a novel technique of differential speckle polarimetry. The observations covering period from 2015 to 2021 allow us to estimate the pattern speed of spiral: $-0.2\pm1.1^{\circ}$/yr (68% credible interval, positive value indicates counter-clockwise rotation), assuming face-on orientation of the disk. This speed is significantly smaller than expected for a companion-induced spiral, if the perturbing body has the semimajor axis of 20 AU. We emphasize that the morphology of the spiral structure is likely to be strongly affected by shadows of a misaligned inner disk detected by Eisner et al. (2004).


INTRODUCTION
Scattered light observations of some protoplanetary disks at high angular resolution revealed spiral structures with characteristic sizes comparable to or several times larger than the Solar system size (Grady et al. 1999;Muto et al. 2012;Grady et al. 2013;Benisty et al. 2015;Muro-Arena et al. 2020). Models predict the spiral waves to emerge due to the influence of a companion (Juhász et al. 2015;Dong et al. 2016) or due to a gravitational instability of the disk . Perhaps the best tool to discriminate between those hypotheses is to measure the pattern speed of the spiral structure using observations at multiple epochs (Ren et al. 2020). In the case of a companion-induced spiral the pattern speed equals to the Keplerian speed of the companion. At the same time, a gravitational instability spiral wave moves at the local Keplerian speed at all distances from the star.
Given the origin of a spiral is known, the analysis of its observational parameters can be indispensable for constraining the conditions in the disk. It can help to determine the position and mass of the companion, which is still embedded in the disk and therefore may be difficult to detect by direct techniques (Zhu et al. 2015;Fung & Dong 2015). This is especially relevant as long as the observations of exoplanets in the process of formation are still rare (see PDS70b for exception (Keppler et al. 2018)). Gravitational instability induced waves can be used to estimate the mass of the disk, one of its most important parameters .
The interpretation of the spiral structures observations has to take into account the effects of uneven and variable illumination, caused by disk selfshadowing (Ren et al. 2020;Xie et al. 2021). These include global variations of brightness of the whole arms (e.g. MWC758 ) or narrow shadow lanes stretching in the radial direction (e.g. HD135344B, Stolker et al. (2016)), or both (e.g. HD100453, Benisty et al. (2017); Long et al. (2017)). Moreover, the shadows are expected to create pressure gradients strong enough to induce the spiral density waves, as was demonstrated by Montesinos et al. (2016). In this case their pattern speed is expected to be the same as in the case of spirals induced by the gravitational instability (Montesinos & Cuello 2018).
A prominent depression was found in the gas and mmsized dust distribution of the CQ Tau circumstellar disk by Ubeira Gabellini et al. (2019) using high angular resolution ALMA observations in 1.3 mm continuum and CO lines. According to Ubeira Gabellini et al. (2019), the inclination and the position angle of the major axis of the gas disk are 35 • and 55 • , respectively. The sizes of gas and large dust (1 µm−1 cm) cavities are ≈ 20 AU and ≈ 52 AU, respectively. Ubeira Gabellini et al. (2019) proposed that the cavity is cleared by a planet with a mass of 6 − 9M Jup and a semi-major axis of 20 AU. Taking into account the mass of the star 1.67M ⊙ the expected orbital period is 75 yr.
The existence of a planet is further supported by the spiral departures from Keplerian kinematics of the gas found by Wölfer et al. (2020). Uyama et al. (2020) reported an excess brightness of the spiral in L ′ -band and proposed that it may be explained by an additional thermal emission due to a disk-planet interaction. The possible connection between the spiral and the putative planet can be tested by measurement of spiral pattern speed, such analysis is favoured for CQ Tau by the short expected orbital period of the planet.
The SED of CQ Tau shows a significant NIR excess, which requires the presence of small dust particles in the cavity (Ubeira Gabellini et al. 2019). The inner disk with typical size of 0.45 AU around the star was detected in K-band using Palomar Testbed Interferometer (Eisner et al. 2004). These facts and UX Ori variability indicate that the inner part of the circumstellar disk may affect the illumination conditions of the outer disk, which should be taken into account during analysis.
In the current paper, we present observations of the circumstellar envelope of CQ Tau using differen-tial speckle polarimetry (DSP) at 48 epochs distributed over the period of more than 5 years in I c -band (effective wavelength is 0.806 µm). DSP is a novel technique based on analysis of series of short-exposure frames obtained through a dual-beam polarimeter without an adaptive optics correction. The spiral structure with the characteristic size of 0.1 − 0.15 ′′ , which corresponds to 15−25 AU at the distance of the source, was detected in the majority of cases; our observations allow us to trace the evolution of its morphology. The work is organized as follows. Sections 2 and 3 contain the description of the observational data and processing methods, respectively. In the section 4, we compare our data with those available in the literature and investigate the variability of the morphology of the spiral structure of CQ Tau. The conclusions are given in section 5.

OBSERVATIONS
The observations were made using SPeckle Polarimeter (SPP), which is a dual-beam polarimeter and speckle camera equipped with an Electron Multiplying CCD (EMCCD) sensor Andor iXon 897. The instrument is installed at the 2.5-m telescope of Caucasian Mountain Observatory of Sternberg Astronomical Institute of Lomonosov Moscow State University (Kornilov et al. 2014). Splitting of the orthogonally polarized beams is done by a Wollaston prism, both images formed by the latter are located side-by-side on the same detector, see Fig. 2. The instrument includes a rotating halfwave plate (HWP) for the modulation of the polarization state of incoming radiation. The atmospheric dispersion is compensated by a dedicated unit consisting of two rotatable direct vision Amici prisms. The plate scale of the instrument is 0.0206 ′′ /pix. The plate scale determination and calibration of the position angle of the instrument is performed by observations of wide binaries, as described by Safonov et al. (2017). The field of view of the instrument is a rectangle with the dimensions of 5 ′′ × 10 ′′ .
The SPP can operate in several bands of visible wavelength range, although in the present work we use the observations in the I c band only, at an effective wavelength λ = 806 nm. At the wavelengths longer than 460 nm, including the whole I c band, the detector oversamples the diffraction limited point spread function of the feeding telescope, which has an aperture diameter D = 2.5 m. The main detector of SPP has a negligible effective readout noise, smaller than 0.1 e − , and can run at high frame rates very close to 100% duty cycle thanks to the frame transfer technology. These features of the SPP, and other EMCCD-based speckle interferometers in general, ensure the reliable detection of speckle struc- ture of the images degraded by the atmospheric optical turbulence, which subsequently allows the computation of the Fourier spectra of these images. In our case the resulting Fourier transforms are used as an input for the differential speckle polarimetry (DSP) algorithm. The detailed description of the instrument design and calibration procedures is presented in paper (Safonov et al. 2017), which also discusses the implementation of the speckle interferometry and integral polarimetry with the SPP.
The observations for the present study were obtained in so called fast polarimetry mode. In this mode the detector obtains a series of short-exposure images at 33 frames per second, which gives exposure time of 30 ms for a single frame. At the same time the HWP continuously rotates at 300 • /sec. The rotation of the HWP serves two purposes. Firstly, it allows to measure both Stokes parameters corresponding to linear polarization. Secondly, it switches the images corresponding to orthogonal polarizations in the beams of the polarimeter each 5 frames, allowing for application of double difference technique (Bagnulo et al. 2009). The double difference discards most of the instrumental effects associated with the differences between two channels of the polarimeter.
In total we obtained 48 series of CQ Tau in the period from October 29th, 2015 to February, 2nd, 2021, all of them in I c band, they are listed in Table 2. The series duration varied from 2000 to 60000 frames (11000-12000 frames in most cases). The integral total flux polarimetry extracted from these observations is presented in Fig. 1 and Table 2. Some of observations were secured while the instrument was installed in the Cassegrain fo-cus, and some in the Nasmyth focus, as indicated in the observational log, Table 2. For the observations in the Nasmyth focus the polarimetric measurements were corrected for the instrumental polarization by the application of the inverse Mueller matrix of the telescope, see details in (Safonov et al. 2017). In some cases we estimated the magnitude by quasi-simultaneous observations of the standard star HIP25001.

Differential polarimetric visibility estimation
All observations were processed by the DSP method. For each frame processing started with the subtraction of bias and background. Then the images F L , F R corresponding to the beams of the polarimeter were extracted and their Fourier transforms F L,i (f ) and F R,i (f ) were computed. Here L and R indices mean left and right beam of the polarimeter. i is the frame number in series, f is the vector of spatial frequency. The distortion of Wollaston prism and atmospheric dispersion compensator prisms was corrected in the Fourier space. These steps are illustrated in Fig. 2 The resulting spectra were combined in the following way (we omit the dependence on spatial frequency f for brevity): Here θ i is the position angle of the HWP for the frame i. * stands for the complex conjugation. The parenthe-  . Scheme illustrating the first steps of the differential speckle polarimetry processing using a single frame as an example. The first row displays the raw frame with two images corresponding to orthogonal polarizations. In the second row the extracted subframes are presented. Third row shows the Fourier transforms of the subframes. First and third panel contain the amplitude, second and fourth contain the argument of the complex Fourier transform.
ses * indicate averaging over frames in series. These equations represent the demodulation of signal, initially modulated by rotation of the HWP. They decompose the signal into harmonics h. As long as the detector obtains ≈ 40 frames per revolution of the HWP, we calculated values (1) and (2) up to N h = 20. N e is the average number of photons in a single frame. The general idea of equations (1,2) is to average the cross-spectra of images corresponding to Stokes Q and I and subsequently normalize them by the average power spectrum of Stokes I image. Previously in (Safonov et al. 2019) (hereinafter S19) we demonstrated that R ch and R sh corresponding to h = 4 can be considered as an estimators of differential polarimetric vis-ibility (DPV) R Q and R U : R ch and R sh for h = 0 characterize the optics after the HWP. Other harmonics (h = 0 and h = 4) are expected to be zero and thus can be used for the noise estimation (S19, appendix B). The factors R Q,ins (f ) and R U,ins (f ) characterize the effect of instrumental polarization, they are discussed in below. The DPV is defined by equations: Here I, Q, U are the Fourier transforms of the Stokes parameters distributions in the object. One can see that DPV is the ratio of visibilities of the object in orthogonal polarizations and can be defined for two Stokes parameters describing linear polarizations Q and U . The measurement of modulus of R Q and R U was demonstrated by Norris et al. (2012) using adaptive optics assisted aperture masking at NaCo/VLT. Unlike that experiment, in SPP the full pupil of the 2.5-m telescope is employed, which provides larger throughput of the system. Apart from this, no adaptive optics correction is used in our case. The signal in the Fourier spectrum of the image at high frequencies is lower than it would be with such a correction. On the other hand, the optical system is simpler and more suitable for the precision polarimetry. Note that DSP allows one to obtain not only modulus, but also the argument of R Q and R U .
The quantities in equations (1-5) are defined in the instrument reference system. The formulae for their conversion in horizontal and equatorial reference frames are given in S19, appendix B. The correction for instrumental terms R Q,ins and R U,ins is done in horizontal reference frame, as long as these factors are associated with the telescope, which has an alt-az mount. This operation is described in Appendix A. Then R Q and R U are converted to equatorial reference frame. For observations secured in the Cassegrain focus the instrumental terms R Q,ins and R U,ins are closer than 10 −4 to unity, much smaller than the typical level of signal for CQ Tau. Therefore the correction for the instrumental polarization effect was not applied for this focus.
As an example of measurement of R Q and R U , we provide the data for CQ Tau obtained in the Cassegrain focus at UT 2015-10-29 00:51 in Fig. 3. Here the DPV can be traced to frequencies 0.7 − 0.8f c , where f c is Figure 3. Columns 1-4: the measurements of RQ and RU for CQ Tau, unpolarized star HIP38325, polarization standard HD19820, from the top to the bottom. Columns 1-2 contain the absolute part of signal, the columns 3-4 contain its argument. The axes are shown in spatial frequency, normalized by cut-off frequency λ/D, where D = 2.5 m is the aperture diameter. The rightmost column presents the image of the envelope in the polarized intensity relative to the total flux of the star, reconstructed by the algorithm presented in section 3.2. Dashes indicate the orientation of polarization. Please note different brightness scales, the unit is polarized flux per pixel divided by the total non-polarized flux from the star. Grey circle in the upper right image indicates typical angular resolution. In all images the north is up, the east is left.
the cut-off frequency D/λ. D = 2.5 m is the telescope aperture diameter, λ = 0.806 µm is the wavelength of the observation. Hence the cut-off frequency f c = 3.1 × 10 6 RAD −1 in our case.
In order to characterize the quality of data we estimated the noise of R Q and R U using method presented in S19. Then we average the noise in an annular region of Fourier space limited by 0.3f c < |f | < 0.4f c . For the series obtained at 2015-10-29 00:51 for CQ Tau it is 6.5 × 10 −3 RAD, in terms of RMS. This quantity computed for all series is presented in seventh column of Table 2.
The Fig. 3 also contains measurement for the unpolarized star HIP38325 secured at UT 2015-10-30 02:40 and the polarization standard star HD19820 secured at UT 2019-12-12 20:16. In both cases the I c band was used. Both stars have no circumstellar dusty envelopes.
One can see that R Q and R U for CQ Tau significantly deviate from the constant level. This behaviour is different from the unpolarized star and polarization standard star. The structure of |R Q | and |R U | for CQ Tau exhibits butterfly pattern, which indicates scattering envelope, as in the case of W Hya, reported by Norris et al. (2012). Interestingly, there is significant structure in the argument of R Q and R U as well, favouring an asymmetry of the polarized flux distribution.

Reconstruction of the image of the envelope in polarized intensity
If the object can be represented as a sum of a point-like star and an extended envelope, the Fourier transforms of Stokes parameters distributions will have the following form: Here I ⋆ is the Fourier transform of the image of the star, a value equal to the flux from the star I ⋆ for all frequencies. q ⋆ , u ⋆ are the dimensionless Stokes parameter of the star. I e (f ), Q e (f ), U e (f ) -Fourier transforms of the envelope image in Stokes parameters I, Q, U , respectively.
In the following we take into account that the polarization of the star is small q ⋆ ≪ 1, u ⋆ ≪ 1, and the envelope is faint in comparison with the star Q e ≪ I ⋆ , U e ≪ I ⋆ . Substituting (6-8) into (5) and keeping only first order terms, we obtain for the Fourier transforms of the envelope flux distributions: Here the quantities on the left side of equations are normalized by the total stellar flux and additionally include the contribution of polarized direct stellar radiation: As one can see from Fig. 3, R Q and R U , and subsequently Q ′ e (f ) and U ′ e (f ) have acceptable signal-tonoise ratio up to frequencies 0.7-0.8f c . Therefore we multiplied them by low-pass filter in order to suppress the wings of point spread function in the image space and reduce the effect of the noise. For this spatial filter we used the optical transfer function (OTF) of a circular aperture: Here f s = 1/α s is the sampling frequency, α s is the angular scale of the detector. Therefore the image produced by applying inverse Fourier transform to them will have the pixel size of α s . We artificially increased the region in which Q ′ e (f )G(f ) and U ′ e (f )G(f ) are defined two times by padding it with zeros. This operation did not add or remove any information from the data, however the pixel size in resulting image decreased to α s /2, making the appearance of the object in the image space smoother.
Filtered and padded Q ′ e (f ) and U ′ e (f ) were transformed into the image space by applying inverse FFT. This reconstruction process results in the estimation of polarized intensity in the circumstellar environment Q ′ (α) and U ′ (α), where α is a two-dimensional vector of spatial coordinate relative to the point-like unpolarized star. The quantities Q ′ (α) and U ′ (α) are then used for the computation of the polarized intensity I p (α) and the angle of polarization χ(α): Recall that the polarized intensity is the product of the total intensity and the fraction of polarization. The examples of I p (α) and χ(α) are presented in Fig. 3 in the rightmost column. As long as spatial filter is non-zero up to frequencies 0.7f c ≡ 2.17 × 10 6 RAD −1 , the reconstructed image of the circumstellar envelope has typical resolution of ≈ 0.1 ′′ . For the unpolarized star HIP38325 the image does not contain any significant structure. On the other hand, in the image of CQ Tau envelope, a spiral arm in 0.12 ′′ to the North from the star is present. There is also a less prominent feature in 0.1 ′′ to the South from the star. Both features demonstrate an azimuthal polarization pattern which indicates that they are parts of a reflection nebula associated with the star.
Note that the estimations Q ′ (α) and U ′ (α) contain the contribution from the polarized direct radiation of the star. This effect is prominent in the case of polarization standard HD19820. In the bottom right panel of Fig. 3 one can see a bright source of polarized radiation coinciding with the star. In case of CQ Tau there is such a source as well, which is discussed in section 4.2.

Comparison with data from literature
The fact that the northern spiral is brighter and farther from the star than the southern one suggests that it lies on the far edge of the flared disk, which has more favorable conditions for scattering the stellar radiation. This assumption is supported by the kinematics data of Wölfer et al. (2020), which show that the disk rotates counter-clockwise and the northwest edge is farther.
Two spiral arms were detected in the disk of CQ Tau earlier by Uyama et al. (2020) in H band using Subaru/HiCIAO. The approximation of the spiral arms from their work is compared to our image in Fig. 4. As one can see, the brighter spiral found by Uyama et al. (2020) is slightly more distant from the star and at approximately the same position angle as in our images. Larger distance from the star is probably explained by the longer wavelength in case of HiCIAO observations: 1.6 µm. Indeed, at the wavelength of our observation 0.8 µm, the optical depth rises faster along the line of sight connecting the surface of the disk and the star. Thus one can expect that the region dominating in scattering light is closer to the star. The fainter arm from (Uyama et al. 2020) at 0.3 ′′ from the star turned out to be inaccessible to us. DSP is less sensitive to fainter outer envelopes than adaptive optics assisted coronagraphs like HiCIAO due to less efficient separation of the radiation of bright central source and faint circumstellar material.
The comparison with ALMA data in 1.3 mm continuum by Ubeira Gabellini et al. (2019) is presented in the same Fig. 4. In millimeter waves the envelope appears as a relatively large ring. The structures found by Uyama et al. (2020) and in present work are located inside this ring. Presumably the spiral in scattered light lies on the inner side of dust torus facing the star. We note that ALMA traces mm-sized dust particles, mean-while visible and NIR radiation is scattered by much smaller particles. Spirals in scattered light inside a mm-wave dust ring was reported for EM * SR 21 by (Muro-Arena et al. 2020).
Palomar Testbed Interferometer observations of Eisner et al. (2004) in K band demonstrated the existence of an inner disk around CQ Tau with the characteristic size of 2.7 mas (gaussian fit), which corresponds to 0.45 AU at the distance of the source. The major axis of the inner disk has P.A.=104 ± 6 • , the inclination of the disk is 48 ± 5 • . The disk with these parameters is indicated in Fig. 4 by the red line. One can see that the orientations of the inner disk and the larger disk found by mm-wave observations are very different.
If the inner disk is geometrically thick, its inclination can be even larger. The fact that CQ Tau is a UX Ori variable (Berdyugin et al. 1990b;Natta et al. 2000;Shakhovskoj et al. 2005) favours a larger inclination, the variable stars of this type display irregular brightness minimums induced by passages of dust clouds across the line of sight (Grinin 1988). The minimums are accompanied by the rise of polarization. The UX Ori variability is typical for young stellar objects with highly inclined disks. On the basis of high resolution ESPaDOnS spectra Dodin & Suslina (2021) found recently that the clouds, occasionally eclipsing CQ Tau, has structures comparable in size with the stellar diameter.

Is the spiral induced by a planet?
The pattern speed of a spiral induced by a planet is defined by the orbital motion of that planet. Therefore if the spiral detected by Uyama et al. (2020) and confirmed by our observations is indeed induced by the planet proposed by Ubeira Gabellini et al. (2019), its pattern speed should be ≈ 4.8 • /yr, which corresponds to orbital period of 75 yr. The direction of motion should be counter-clockwise.
To test this hypothesis we considered the following model of object variability. Firstly, as long as CQ Tau is a UX Ori variable, we allowed for a change in brightness and total polarization of an unresolved source. Secondly, we assumed that the envelope can rotate as a whole. For rotation modeling we did not deproject the disk because, according to Ubeira Gabellini et al. Using these assumptions we approximated each individual measurement of DPV with a reference one. For the reference observation we took the data obtained at 2018-11-25 23:49 which has relatively low noise level and lies roughly in the middle of considered time period. Note that this approximation employs the DPV measurements directly, not the images reconstructed using For each observation we obtained a set of parameters q ′′ ⋆ , u ′′ ⋆ , ζ, η by minimization of χ 2 statistic. The first two of these parameters characterize the change in total polarization of the unresolved source, the third defines the change in its brightness. We interested most in parameter η, which is the rotation angle between the reference observation and a given one. The estimations of η are presented in Fig. 5, as a function of time.
In the same figure their weighted linear approximation by η = η 0 + ωt is presented. η values exhibit significant spread, which is probably due to variable shadowing of the envelope. However the overall trend is absent. The parameters of the linear fit are: η 0 = −0.4 ± 1.6 • , ω = −0.2 ± 1.1 • /yr (68% credible intervals). Recall that the pattern speed expected for a planet-induced spiral is ω e = 4.8 • /yr (direction counter-clockwise), the corresponding trend is indicated by the red line in figure. Therefore we reject the hypothesis that the spiral in scattered light detected in our study is caused by a planet with the orbital period of 75 yr at the level of 4.3σ.

The spiral as an effect of shadowing by the inner disk
The increased optical depth for the lines of sight passing near the equator of the inner disk should lead to the shadows on the outer disk. Such a mechanism was proposed to explain the appearance of disks in scattered light for a number of objects: HD100453 (Benisty et al. 2017), DoAr 44 (Casassus et al. 2018), RXJ1604.3-2130(Pinilla et al. 2018. Given the small inclination of the outer disk, the shadow position should roughly coincide with the line of intersection of the inner and outer disks equator planes. In our case, the brighter northern spiral arm indicates that the northern edge of the inner disk is farther from us. Taking into account the inclinations and position angles of the disks, the expected position angles of the shadow are 140 • , 320 • , as indicated in Fig. 4. Mutual configuration of the disks is illustrated in Fig. 6. The angle between the axes of rotation of the inner and outer disks is 40 • . The orientations of the shadows match the positions of the gaps in the spiral structure. Therefore the spirals detected by Uyama et al. (2020) and in the present work may be actually rings with shadows cast by the inner disk, rather than physical structures in scattering dust.

CONCLUSIONS
The structures found in the protoplanetary disks in scattered light show variability on timescales from days to years. Both physical motion (Ren et al. 2018(Ren et al. , 2020 and illumination changes (Stolker et al. 2016;Benisty et al. 2017;Casassus et al. 2018;Pinilla et al. 2018) are observed. In each specific case it is important to differentiate these two basic explanations. The analysis of variations in the morphology of envelopes is greatly facilitated by multiepoch observations made by a homogeneous technique.
In the current paper, we present observations of the circumstellar envelope of a young star CQ Tau by differential speckle polarimetry in the I c band at the 2.5m telescope of the Caucasian Mountain Observatory of Lomonosov Moscow State University, covering the period from 2015 to 2021. We detect a spiral arm at the characteristic distances 0.1−0.15 ′′ north of the star. The structure is located at the same position angle as the spiral arm discovered previously by Uyama et al. (2020) at a slightly larger stellocentric distance. We also report the presence of a southern, closer to the star, arm of the spiral.
To explain the large gap in the gaseous and dust disk of CQ Tau Ubeira Gabellini et al. (2019) proposed the existence of a planet on the orbit with the semimajor axis of 20 AU. This possible planet could potentially generate structures appearing as spiral arms in scattered light, as was proposed for SAO 206462 (Xie et al. 2021), MWC 758 (Ren et al. 2018), AB Aur (Boccaletti et al. 2020). The pattern speed of the respective spiral in this case is dictated by the orbital motion of the planet: 4.8 • /yr. On basis of observations covering the period of more than 5 yr we estimated the spiral pattern speed of CQ Tau to be −0.2 ± 1.1 • /yr (68% confidence interval). We conclude that the spiral in scattered light is not caused by that planet. However it is not excluded that the spiral is generated by a planet further away from the star, at distances of 40-50 AU.
The relative stability of the spiral structure may indicate that its morphology is largely determined by shadowing of the outer disk by the misaligned inner one, like in cases of HD100453 (Benisty et al. 2017), DoAr 44 (Casassus et al. 2018), RXJ1604.3-2130(Pinilla et al. 2018. The gaps in the spiral structure seen in scattered light at position angles of ≈ 160 • and ≈ 340 • correspond to the expected positions of the shadows from the inner disk found earlier by Eisner et al. (2004). A slow precession of the inner disk is still possible. In the subsequent analysis and interpretation of the spiral structure of CQ Tau, it is important to consider the possibility of its non-uniform illumination.

ACKNOWLEDGMENTS
We are grateful to Sergei Lamzin for the useful comments on manuscript. Constructive comments by anonymous referee allowed us to improve the analysis and presentation. We acknowledge the financial support of the Russian Science Foundation Public Monitoring Committee 20-72-10011. Scientific equipment used in this study was bought partially by the funds of the M. V. Lomonosov Moscow State University Program of Development.

A. CORRECTION OF INSTRUMENTAL POLARIZATION
The instrumental polarization induced by a mirror can be modelled using the formalism of Stokes vectors and Mueller matrices. The Stokes vector of outgoing radiation is the product of Stokes vector of incoming radiation and the Mueller matrix of the mirror. For the integral polarimetry the Mueller matrix of the mirror is averaged over the beam. Thus, mirrors acting in a symmetrical configuration, e.g. the primary and the secondary, do not modify the polarization of incoming radiation, except for a negligible depolarization. When the SPP is mounted in the Cassegrain focus the whole system up to beamsplitter possesses axial symmetry, thus ensuring low level of the instrumental polarization. Measurements of unpolarized stars demonstrated the instrumental polarization lower than 10 −4 (Safonov et al. 2017).
In the Nasmyth focus an oblique reflection by the diagonal M3 mirror is added, axial symmetry is broken, and the instrumental polarization becomes quite substantial. It varies from 2 to 4% over the wavelength range of SPP. In (Safonov et al. 2017) we constructed a model of instrumental polarization induced by M3 mirror of 2.5-m telescope and provided a method for its correction.
While the instrumental polarization is definitely important for measuring the integral polarization, for diffraction limited imaging it adds a new level of complexity. The fact that the reflection geometry is different for different points in the pupil lead to emergence of so called differential polarization aberrations (Breckinridge et al. 2015). The first order term of these aberrations manifests itself as a phase ramp in R Q,ins and R U,ins , S19. The latter paper contains  extensive analysis of instrumental terms R Q,ins and R U,ins and proposes the following model for them: Here f x , f y are the components of the two-dimensional frequency vector. q ins , u ins is the polarization of an intrinsically unpolarized star detected by the instrument. s ⋆ q,ins , t ⋆ q,ins , s ⋆ u,ins , t ⋆ u,ins are the components of so called instrumental polaroastrometric signal (IPS). Equations (A1,A2) assume that R Q,ins , R U,ins are defined in the horizontal reference system.
The quantities q ins , u ins , s ⋆ q,ins , t ⋆ q,ins , s ⋆ u,ins , t ⋆ u,ins were determined in S19 from measurements of point-like unpolarized stars, see their values for I c band in Table 1. In the Nasmyth focus the IPS is large and should be corrected through the division of observed R c4 and R s4 by (A1) and (A2), respectively. This operation is performed in the horizontal reference system. Fig. 7 demonstrates the same quantities as in the fifth column of Fig. 3, but obtained in the Nasmyth focus and corrected for the instrumental polarization. Again for comparison we give the observations for a point-like unpolarized star and a point-like polarization standard, HIP8370 and HD25443, respectively. One can see that the spiral feature is persistent for the CQ Tau. It is not observed neither for the unpolarized star, nor for the polarization standard.
The image of the polarization standard HD25443 exhibits a ring in polarized light, which corresponds to the first bright ring of the Airy function. Note that this ring has a constant polarization angle equal to the integral polarization angle of the star. In this respect it differs from the azimuthal polarization pattern, as in the case of CQ Tau, which is a genuine feature of a scattering reflection nebula.

B. APPROXIMATION OF ONE OBSERVATION BY ANOTHER
The section 4.2 is based on the approximation of a given observation k by a reference one j. The differential polarization visibility for the reference observation can be written as: according to equations (6,7). Recall that q ⋆ is the polarization of the star, Q e is the Fourier transform of envelope image in Stokes parameter Q, I ⋆ is the stellar flux. If the envelope rotated by η (positive for counter-clockwise rotation), then the DPV will be: The rotated version of Q e,r will depend on both Q e and U e due to transformation of Stokes parameters in the rotated reference system. For the equations describing rotation of R Q and R U see appendix B of S19. Suppose that the stellar flux increases by ζ and stellar polarization changes to q ′ ⋆ . The DPV will change in the following way: Expressing Q e from (B4) and substituting it in (B5), we obtain formula describing the transformation of the measurement: Similar relation can be formulated for Stokes U : The second terms in these equations will be considered as new parameters: The comparison of measurements is performed by computation of the residual of the following form: where Here R k Q , R k U are the DPV measurements for the observation k, R j * Q,r , R j * U,r are the same, but for observation j, modified according to equations (B6) and (B7). σ R,k , σ R,j * are the errors of estimation of DPV, the method of their determination can be found in Appendix B of S19.
For each observation we determine four parameters: q ′′ ⋆ , u ′′ ⋆ , ζ, η by minimization of the total residual: Here the summation takes place for the frequencies smaller than 0.7f c , where f c is the cut-off frequency. Note that the total residual defined by (B10) is technically a χ 2 statistic. The posterior probability of parameters was sampled using the affine-invariant ensemble Markov Chain Monte Carlo method (Goodman & Weare 2010). 50 walkers and 2 × 10 4 iterations in total were used. The example of the resulting posterior for observation 2019-10-26 22:10 approximated by 2018-11-25 23:49 is displayed in the left part Fig. 8. The best values for parameters are: q ′′ ⋆ = −0.06 ± 0.07%, u ′′ ⋆ = −1.32 ± 0.07%, ζ = 0.67 ± 0.05, η = 8.8 ± 2.7 (values after the plus-minus sign are 68% confidence intervals). The components of the corresponding residual (B9) are presented in the right part of Fig. 8. The absence of significant departures indicates that the fitting of one observation by modified version of another has satisfactory quality. reference observation we took the data obtained on 2020-10-14 23:00. Then the procedure described in section 4.2 was repeated for this reference. The resulting rotation angles η 1 and η 2 for two references are directly compared in the left panel of the Fig. 9. The CDF of η 1 − η 2 normalized by its uncertainty is presented in the right panel of the same figure. These normalized deviations are expected to follow shifted standard normal distribution. The shift emerges because the envelope may have rotated between moments 2018-11-25 23:49 and 2020-10-14 23:00. One can see that the distribution of deviations indeed follows the normal distribution for the majority of data points. However the wings of the distribution are larger than expected.
The dependence of rotation angles η, determined using the alternative reference observation, on time were fitted with the linear law, as it was done in section 4.2. The parameters of the fit are η 0 = 1.4 ± 1.7 • , ω = −0.7 ± 0.9 • /yr (68% credible intervals). They coincide within errors with the values obtained using the observation for 2018-11-25 23:49 as a reference. We conclude that the choice of the reference observation has no major effect on the results of section 4.2. 1.21 ± 0.04 159.9 ± 1.0 9.6 ± 0.1 0.0137 32 -54 0.90 Table 2 continued Note-Focal stations: C -Cassegrain, N -Nasmyth. p and χ are the fraction and angle of polarization, respectively. Ic band magnitude estimated by quasi-simultaneous observations of HIP25001. At some epochs these observation were not conducted. σ is uncertainty of R averaged over region 0.3fc < |f | < 0.4fc, see Section 2. z is zenith angle, h⊙ is the Sun altitude, β is seeing measured by MASS-DIMM instrument (Kornilov et al. 2014). The exposure of a single frame was 30 ms for all observations. All observations were conducted in Ic band, which effective wavelength is 0.806 nm.