# External or internal companion exciting the spiral arms in CQ Tau?

Iain Hammond, Valentin Christiaens, Daniel J. Price, Maria Giulia

Ubeira-Gabellini, Jennifer Baird, Josh Calcino, Myriam Benisty, Giuseppe

Lodato, Leonardo Testi, Christophe Pinte, et al.

► To cite this version:

Iain Hammond, Valentin Christiaens, Daniel J. Price, Maria Giulia Ubeira-Gabellini, Jennifer Baird, et al.. External or internal companion exciting the spiral arms in CQ Tau?. Monthly Notices of the Royal Astronomical Society, 2022, 515, pp.6109-6121. 10.1093/mnras/stac2119. insu-03860276

HAL Id: insu-03860276

https://insu.hal.science/insu-03860276

Submitted on 17 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# External or internal companion exciting the spiral arms in CQ Tau?

Iain Hammond<sup>1,★</sup> Valentin Christiaens<sup>1,2</sup> Daniel J. Price<sup>1</sup> Maria Giulia Ubeira-Gabellini,<sup>3</sup>

Jennifer Baird,<sup>1</sup> Josh Calcino,<sup>1,4</sup> Myriam Benisty,<sup>5,6</sup> Giuseppe Lodato,<sup>1,3</sup> Leonardo Testi,<sup>7,8</sup>

Christophe Pinte<sup>1,5</sup> Claudia Toci<sup>1,3</sup> and Davide Fedele<sup>8,9</sup>

Accepted 2022 July 25. Received 2022 July 19; in original form 2022 February 21

## ABSTRACT

We present new high-contrast images in near-infrared wavelengths ($\lambda_c$ = 1.04, 1.24, 1.62, 2.18, and 3.78 $\mu$m) of the young variable star CQ Tau, aiming to constrain the presence of companions in the protoplanetary disc. We reached a $K_s$-band contrast of 14 mag with SPHERE/IRDIS at separations greater than 0.4' from the star. Our mass sensitivity curve rules out giant planets above 4 $M_{Jup}$ immediately outside the spiral arms at ~60 au and above 2–3 $M_{Jup}$ beyond 100 au to 5$\sigma$ confidence assuming hot-start models. We do, however, detect four spiral arms, a double-arc and evidence for shadows in scattered light cast by a misaligned inner disc. Our observations may be explained by an unseen close-in companion on an inclined and eccentric orbit. Such a hypothesis would also account for the disc CO cavity and disturbed kinematics.

**Key words:** techniques: image processing – planet–disc interactions – protoplanetary discs.

## 1 INTRODUCTION

The advent of dedicated high-contrast facilities, such as the near-infrared (NIR) Spectro-Polarimetric High contrast imager for Exoplanets REsearch (SPHERE; Beuzit et al. 2008) instrument and the Atacama Large Millimeter Array (ALMA), triggered a revolution in imaging of protoplanetary discs. The enhanced angular resolution and sensitivity has led to the discovery of substructures in protoplanetary discs, such as annular gaps (e.g. Andrews et al. 2018; Long et al. 2018), inner cavities (e.g. Canovas et al. 2015; Villenave et al. 2019), crescent-shape asymmetries (e.g. Casassus et al. 2013), spiral arms (e.g. Muto et al. 2012; Rosotti et al. 2020), and shadows (e.g. Keppler et al. 2020), which are interpreted as possible signposts for the presence of embedded companions.

Hydrodynamical simulations of perturbing bodies in protoplanetary discs have had some success in recreating observed substructures. A few examples of discs where these models have been applied are HD 142527 (e.g. Avenhaus et al. 2014; Price et al. 2018), MWC 758 (e.g. Benisty et al. 2015; Dong et al. 2015; Reggiani et al. 2018; Calcino et al. 2020), IRS 48 (e.g. van der Marel et al. 2013, 2016; Calcino et al. 2019), HD 100453 (e.g. Benisty et al. 2017; Gonzalez et al. 2020; Rosotti et al. 2020), HD 169142 (e.g.

Bertrang et al. 2018; Toci et al. 2020b), HD 100546 (e.g. Follette et al. 2017; Fedele et al. 2021), and PDS 70 (e.g. Keppler et al. 2018; Toci et al. 2020a). In some cases, different companions have been proposed to explain the same substructures (e.g. for MWC 758 from Calcino et al. 2020 and Baruteau et al. 2019).

Spiral arms have long been hypothesized to arise due to planet–disc interactions (Goldreich & Tremaine 1979, 1980; Ogilvie & Lubow 2002). The amplitude of the spirals depends on the mass of the planet with respect to the so-called thermal mass (the mass for which the Hill radius equals the local gas scale height) and properties of the disc such as the temperature profile. In general, low-mass planets (e.g. masses much lower than the thermal mass) are expected to produce a single inner and outer wake. When the planetary mass increases to approximately the thermal mass or higher, multiple spirals are generated (Bae & Zhu 2018), both interior and exterior to the planet. Perturbations from such a massive planet external to the spiral arms initially seemed to be an exciting and robust explanation for the apparent *m* = 2 spirals detected in several protoplanetary discs (most notably MWC 758 and HD 135344B; Garufi et al. 2013; Benisty et al. 2015). However, no convincing detection has emerged, despite mass limits of the order of 3 $M_{Jup}$ exterior to the spirals assuming hot-start formation scenario for the cases of HD 135344B and MWC 758 (Maire et al. 2017; Reggiani et al. 2018; Boccaletti et al. 2021).

Radial velocity and direct imaging surveys have found super-Jupiter mass planets with non-negligible eccentricities (Bitsch,

<footnote n="1"><sup>1</sup>School of Physics and Astronomy, Monash University, Vic 3800, Australia</footnote>

<footnote n="2"><sup>2</sup>Space sciences, Technologies & Astrophysics Research (STAR) Institute, Université de Liège, Allée du Six Août 19c, B-4000 Sart Tilman, Belgium</footnote>

<footnote n="3"><sup>3</sup>Dipartimento di Fisica, Università Degli Studi di Milano, Via Celoria, 16, Milano, I-20133, Italy</footnote>

<footnote n="4"><sup>4</sup>Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA</footnote>

<footnote n="5"><sup>5</sup>Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="6"><sup>6</sup>Unidad Mixta Internacional Franco-Chilena de Astronomía (CNRS, UMI 3386), Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile</footnote>

<footnote n="7"><sup>7</sup>European Southern Observatory (ESO), Karl-Schwarzschild-Str 2, D-85748 Garching, Germany</footnote>

<footnote n="8"><sup>8</sup>INAF, Osservatorio Astrofisico di Arcetri, Largo Enrico Fermi 5, I-50125 Firenze, Italy</footnote>

<footnote n="9"><sup>9</sup>INAF, Osservatorio Astrofisico di Torino, Via Osservatorio 20, I-10025 Pino Torinese, Italy</footnote>

<footnote>*E-mail: iain.hammond@monash.edu</footnote>
<table><caption>**Table 1.** Stellar and disc properties of CQ Tau.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Reference</th></tr></thead><tbody><tr><td>Distance (pc)</td><td>149.37<sup>+1.35</sup><sub>-1.33</sub></td><td>1</td></tr><tr><td>Age (Myr)</td><td>8.9<sup>+2.8</sup><sub>-2.5</sub></td><td>2</td></tr><tr><td>Spectral type</td><td>F2</td><td>3, 4</td></tr><tr><td>T<sub>eff</sub> (K)</td><td>6850 ± 100</td><td>2, 5</td></tr><tr><td>Av (mag)</td><td>0.41</td><td>2</td></tr><tr><td>Log(L) (L<sub>☉</sub>)</td><td>0.87</td><td>2</td></tr><tr><td>M (M<sub>☉</sub>)</td><td>1.47</td><td>2</td></tr></tbody></table>

References: (1) Gaia Collaboration (2021), (2) Vioque et al. (2018), (3) Herbig (1960), (4) Natta et al. (2001), and (5) Dodin & Suslina (2021).

Trifonov & Izidoro 2020). If these eccentric orbits were present during the protoplanetary disc phase, the resulting spiral structure would be drastically different to the spiral morphology expected from perturbers on circular orbits (Zhu & Zhang 2022). Calcino et al. (2020) showed that a super-Jupiter-mass planet on an eccentric orbit can produce multiple spiral arms outside of a central cavity. Several studies (e.g. Muley, Fung & van der Marel 2019; Baruteau et al. 2021) have found massive planets become eccentric as they begin to carve deep gaps which deplete the co-orbital Lindblad torques that tend to quickly damp any growth in eccentricity (Goldreich & Tremaine 1980; Artymowicz 1993). In this context, the peculiar spiral morphology of CQ Tau (Uyama et al. 2020) makes it a promising laboratory to better understand the origin of spiral arms in protoplanetary discs.

CQ Tau is a young (~9 Myr, Vioque et al. 2018) variable star of intermediate mass (1.47 M<sub>☉</sub>, Vioque et al. 2018), located near the Taurus star-forming region (RA = 05h:35m:58.47s, Dec. = +24°44:54.09; J2000) at a distance of 149 pc (Gaia Collaboration 2021). It is spectral-type F2 (see Table 1) and shows UX Ori-type irregular photometric variability, possibly indicative of the presence of a companion (Section 4.3). CQ Tau is surrounded by a massive circumstellar disc, historically observed in continuum emission with multiple instruments (OVRO interferometer; Mannings & Sargent 1997; PdBI, Natta, Grinin & Mannings 2000; VLA, Testi et al. 2001). Its spectral energy distribution (SED) shows infrared and mm-wavelength emission, allowing us to infer the presence of grains larger than those in the interstellar medium (Testi et al. 2001, 2003; Chapillon et al. 2008), up to cm-sized (Banzatti et al. 2011). There is evidence of grain growth variation with larger grains present in the inner disc (Trotta et al. 2013). Moreover, CQ Tau has a relatively high-accretion rate of the order of $10^{-8}-10^{-7}$ M<sub>☉</sub>yr<sup>-1</sup> (Donehew & Brittain 2011; Mendigutía et al. 2012; Guzmán-Díaz et al. 2021).

High-resolution observations tracing the mm-dust (e.g. Tripathi et al. 2017; Pinilla et al. 2018; Ubeira Gabellini et al. 2019) and gas distribution (e.g. Ubeira Gabellini et al. 2019; van der Marel et al. 2021; Wölfer et al. 2021) shows the presence of an inner cavity. Ubeira Gabellini et al. (2019) analysed the dust and gas surface density profiles, finding a Gaussian dust ring at a radius of 53 au in the 1.3mm continuum (also Francis & van der Marel 2020) and a less prominent gas cavity of 20–25 au in the <sup>13</sup>CO and C<sup>18</sup>O emission. Based on the hydrodynamical simulations, the authors suggested it to be caused by an undetected planet of mass 6–9 M<sub>Jup</sub> located at ~20 au. Recently, Uyama et al. (2020) presented L'-band Keck/NIRC2 and H-band polarimetric Subaru/AO188 + HiCIAO observations of CQ Tau. In the small grain distribution they detected a spiral arm with a high pitch angle (34 ± 2°), possibly induced by an unseen companion candidate, as well as a potential secondary spiral. No companion was detected, despite mass sensitivity limits

of ~7 M<sub>Jup</sub> immediately beyond the prominent NIR spiral (~60 au; Uyama et al. 2020). This is in tension with hydrodynamical models (Sanchis et al. 2020), which predict a 1–2 M<sub>Jup</sub> planet should be readily detectable at radii ≥5–10 au, with a contrast of 14–16 mag, respectively, due to negligible disc extinction in L'-band. Finally, Wölfer et al. (2021) analysed the gas brightness temperature and kinematics of the disc, finding non-Keplerian gas substructures, such as three spiral structures spanning radii from ~10 to 180 au and connected with the spirals observed in NIR scattered light, and a nearly symmetric temperature dip approximately co-spatial with the semiminor axis of the disc.

In this paper, we first present new NIR images and constraints on possible companions in the outer disc of CQ Tau. We present our data reduction and post-processing in Section 2, our spiral arm characterization and detection limits in Section 3, and discuss the indicators of a companion in the disc of CQ Tau in Section 4.

## 2 OBSERVATIONS AND DATA REDUCTION

The observations used in this work are summarized in Table 2. All are new to this work, with the exception of program 098.C-0760 (2017 October 6), which is presented in appendix E of Bohn et al. (2022). Together, they represent all currently available ESO high-contrast imaging data for CQ Tau. Section 2.1 presents our PDI reduction, Section 2.2 our IRDIFS reduction, and Section 2.3 our NaCo reduction and new pipeline.

### 2.1 VLT/SPHERE polarimetric data sets

Two polarimetric data sets were captured with the Infrared Differential Imaging Spectrometer (IRDIS; Dohlen et al. 2008) instrument in dual-polarization mode (de Boer et al. 2020; van Holstein et al. 2020b) on 2017 October 6 and 2018 February 18 with the SPHERE instrument at the Very Large Telescope (VLT). The exposure time was set to 50 s for individual coronagraphic (apodized Lyot coronagraph N_ALC_YJH_S, diameter ~185 mas) images and 0.8375 s for non-coronagraphic images providing a total of 10 and 18 polarimetric cycles, respectively. The J-band filter ($\lambda_c$ = 1.24 µm) was used for both observations.

Both data sets were reduced using the IRDIS Data reduction for Accurate Polarimetry pipeline (IRDAP v1.3.3; van Holstein et al. 2020a, b). This automated pipeline incorporates image calibration (dark subtraction, flat-fielding, and correction of bad pixels), pre-processing and post-processing for polarimetric data. The position of the central star was interpolated using the two centre calibration files; one taken at the beginning of the observation and other at the end. For each polarimetric cycle, four cubes each corresponding to linear polarization directions of $Q^+$, $Q^-$, $U^+$, and $U^-$ were obtained. Using these frames, the clean Stokes $Q$ and $U$ were determined using the standard double-difference method, corrected for instrument polarization with the Mueller matrix model and converted to polar coordinates ($Q_\phi$ and $U_\phi$). By fitting a 2D Gaussian to the unsaturated flux frames of the coronagraphic data set, we measured a full width at half-maximum (FWHM) of 3.61 pixels after excluding one bad unsaturated cube.

We explored polarized phase function effects by producing two estimated total intensity images with the coronagraphic $Q_\phi$ image. The first using an in-house processing routine assuming typical grain sizes smaller than the wavelength of the observation (1.24 µm). For independent confirmation, we used the total intensity and phase function routines in DISKMAP (Stolker et al. 2016) using
<table><caption>**Table 2.** Summary of the archival VLT high-contrast imaging observations on CQ Tau used in this work.</caption><thead><tr><th>Date</th><th>Strategy</th><th>Program</th><th>Instrument</th><th>Filter</th><th>Plate scale<sup>(a)</sup><br>(mas px<sup>-1</sup>)</th><th>Coronagraph</th><th>T<sub>int</sub><sup>(b)</sup><br>(min)</th><th>&lt;β&gt;<sup>(c)</sup><br>(arcsec)</th><th>ΔPA<sup>(d)</sup><br>°</th></tr></thead><tbody><tr><td>2016-12-21</td><td>ADI</td><td>298.C-5014</td><td>IRDIS</td><td>Ks</td><td>12.265 ± 0.009</td><td>N.ALC_YJH_S</td><td>127</td><td>0.69</td><td>44.4</td></tr><tr><td>2016-12-21</td><td>A/SDI</td><td>298.C-5014</td><td>IFS</td><td>YJH</td><td>7.46 ± 0.02</td><td>N.ALC_YJH_S</td><td>128</td><td>0.69</td><td>44.4</td></tr><tr><td>2017-10-06</td><td>PDI</td><td>098.C-0760</td><td>IRDIS</td><td>J</td><td>12.263 ± 0.009</td><td>N.ALC_YJH_S</td><td>33</td><td>0.41</td><td>–</td></tr><tr><td>2018-02-18</td><td>PDI</td><td>098.C-0760</td><td>IRDIS</td><td>J</td><td>12.238 ± 0.009</td><td>–</td><td>20</td><td>0.65</td><td>–</td></tr><tr><td>2018-11-28</td><td>ADI</td><td>1101.C-0092</td><td>NACO</td><td>L'</td><td>27.208 ± 0.088</td><td>AGPM</td><td>157</td><td>0.55</td><td>67.4</td></tr></tbody></table>

(a)Maire et al. (2016) for IRDIS and IFS, Launhardt et al. (2020) for NaCo. (b)Total integration time calculated after bad frame removal.

(c)Average seeing at λ = 500 nm achieved during the sequence, returned by ESO's DIMM station. (d)Cumulative parallactic angle variation.

our r<fnref n="2" />-scaled image mapped with the scattering surface described in Ubeira Gabellini et al. (2019).

## 2.2 VLT/SPHERE IRDIFS data set

CQ Tau was observed by SPHERE in IRDIFS-EXT mode on 2016 December 21 (Program ID 298.C-5014(A), PI: L. Testi). This mode allowed the simultaneous acquisition of Integral Field Spectroscopy (IFS; Claudi et al. 2008) and IRDIS data, in the YJH ($\lambda_c$ = 1.04, 1.24, and 1.62 µm) and Ks ($\lambda_c$ = 2.18 µm) bands, respectively. Observations were made with a stabilized pupil to enable the application of angular differential imaging-based (ADI, Marois et al. 2006a) post-processing algorithms to model and subtract the emission from the star. For both data sets, the N.ALC_YJH_S coronagraph (inner working angle ~99 mas) was used and the total integration time was 128 min (Table 2). IRDIS images had an individual exposure time of 16 s, while IFS images were 32 s each. For both IFS and IRDIS, images containing satellite spots (Marois et al. 2006b; Langlois et al. 2013; labelled CENTER) were acquired every ~31 min for a total of five CENTER cubes.

### 2.2.1 IFS

The IFS data are acquired with a spectral resolution of R ~ 30 and span a field of view of 1".73 × 1".73. We performed basic data calibration using the ESO SPHERE pipeline<fnref n="1" /> (v0.36), involving dark subtraction, detector flat-fielding, anamorphism correction (Maire et al. 2016), wavelength calibration, and spectral cube construction.

Bad pixel correction and satellite spot-based centring of the images utilized routines implemented in the Vortex Image Processing<fnref n="2" /> package (VIP; Gomez Gonzalez et al. 2017). In order to increase the signal to noise (S/N) of the four satellite spots used for centring, we applied a high-pass filter to our CENTER images where the image itself was subtracted with a median low-pass filtered version of the image. We used a 2D Moffat function to fit the satellite spot positions and to derive the centroid of the star for each CENTER image. For the OBJECT images, the stellar position behind the coronagraph was interpolated using the CENTER images interspersed throughout the sequence.

We post-processed our calibrated, centred cubes using principal component analysis (PCA) functions implemented in VIP (Amara & Quanz 2012; Gomez Gonzalez et al. 2017). We used two flavours of PCA, leveraging either the spectral diversity alone with spectral differential imaging (SDI) or both the spectral and angular diversity present in the data (ASDI; Ubeira-Gabellini et al. 2020).

### 2.2.2 IRDIS

IRDIS data consisted of 160 science cubes. We used an in-house reduction pipeline (Christiaens et al. 2021) to reduce the data set, after adaptation to both support satellite spot centring and handle classical imaging data. The pipeline uses ESO's Common Pipeline library ESOREX (v3.13.5) recipes for calibration, and VIP functions for pre- and post-processing. The breakdown of the pipeline for this observation is as follows: (i) dark current subtraction, (ii) flat-fielding, (iii) PCA-based sky subtraction, (iv) bad pixel correction, (v) centring via satellite spots, (vi) anamorphism correction, (vii) finer centering using a background star (see Appendix B), (viii) bad frame rejection (two frames), (ix) final PSF creation, (x) cropping, and (xi) post-processing with median-ADI, PCA-ADI in full-frame and annular PCA-ADI. We measured a FWHM of 5.24 pixels.

## 2.3 VLT/NACO ADI data set

CQ Tau was observed with VLT/NaCo (Lenzen et al. 2003; Rousset et al. 2003) on the night of 2018 November 28 in L'-band ($\lambda_c$ = 3.78 µm) as part of the Imaging Survey for Planets around Young stars (Launhardt et al. 2020) program. An Annular Groove Phase Mask (AGPM; Mawet et al. 2005, 2013) was used to suppress stellar signal. Science frames were captured in pupil-stabilized mode, enabling the application of classical ADI and PCA-ADI post-processing algorithms. The integration time was 0.35 s per frame.

Since the ESO pipeline does not support coronographic imaging, we developed a new NaCo data reduction pipeline<fnref n="3" /> utilizing functions from VIP, and used it to calibrate, pre- and post-process the data.

First, the systematically bad lower-left quadrant, defined by every eighth column returning the saturated detector count, was corrected in each cube by detecting pixel columns that featured the saturated detector count and were replaced with the median of the neighbouring 5 × 5 unaffected pixels. Cubes were automatically classified based on header keywords. For the purpose of ADI, the angle between north and the vertical axis for each science image was calculated using the starting and ending parallactic and pupil position angles. Subsequently, the rotation required to align each frame with true north was found using linear interpolation. Each derotation angle was corrected for the inherent offset from true north for NaCo (0.572 ± 0.178°; Launhardt et al. 2020).

We used a new PCA dark subtraction technique, adapted from the PCA-based sky subtraction routine implemented in VIP, which aims to model and subtract the dark current (Appendix A). The flat frames obtained at different airmasses were sorted based on their median pixel value and used to create a master flat frame.

<footnote n="1">¹http://www.eso.org/sci/software/pipelines/sphere/</footnote>

<footnote n="2">²https://github.com/vortex-exoplanet/VIP</footnote>

<footnote n="3">³https://github.com/IainHammond/NACO_pipeine</footnote>
All science, sky, and unsaturated frames were then divided by the master flat frame. ‘Hot’ and ‘cold’ pixels were detected using the master flat and replaced with the median value of neighbouring 5 × 5 unaffected pixels. Transient bad pixels (e.g. due to cosmic rays) and NaN values were also identified and sigma-filtered in 5 × 5 pixel-boxes. NaCo’s inconsistent detector integration time in the first few frames of each cube was accounted for by measuring the flux in an annulus from 3 arcsec to the edge of the frame. The first six frames systematically featured a flux $2\sigma$ above the median of their respective cube and were hence removed, and the remaining 94 rescaled to the median flux. A low-pass filter was applied to each unsaturated non-coronagraphic frame to determine the approximate location of the star on either of the three good detector quadrants. For each unsaturated frame, images where the star was dithered on other quadrants of the detector were used for sky subtraction to isolate the stellar PSF. A 2D Gaussian was fit to a subframe containing the median-combined PSFs, resulting in a FWHM of 4.19 pixels.

Each science cube was paired with a sky cube captured closest in time and aligned with respect to each other using dust grains present in both cubes. We applied VIP’s PCA-based sky subtraction algorithm with one principal component to remove the background flux, using the science cube minus the median frame of the sky cube as input and the sky cube minus the median frame of the sky cube as the PCA library.

We performed a two-step centring algorithm on the science cubes to correct for the uncertainty in location of the star behind the AGPM centre. Each science cube was first duplicated and median-combined into a single frame to enhance the speckle pattern around the star. A high-pass filter removed large-scale variations greater than the FWHM. By looping over each median-combined frame, we: (i) performed a global centring by fitting a negative 2D Gaussian to the area around the AGPM centre and (ii) aligned frames with respect to each other by maximizing the cross-correlation of the speckle pattern expressed in logarithmic intensity. This iterative process was repeated five times to ensure the error in the centring is small. Final centring shifts were applied to individual (not median) frames of each science cube. Any frame with pixel shifts discrepant by more than 0.4 pixel with respect to the median pixel shifts along x and y was rejected. Additionally, we calculated the Pearson correlation coefficient between individual frames to the master median-combined frame in a 31 × 31 pixel subframe around the centre, rejecting frames with a correlation factor $1\sigma$ below the median correlation factor.

The final master cube contained 28 010 frames, temporally subsampled to 2801 frames by taking the median of every 10 frames. We post-processed the cube using VIP implementations of median-ADI and PCA-ADI algorithms (both in full-frame and concentric annuli). We explored 1–20 principal components, with a 1 × FWHM inner mask and 1–3 × FWHM rotation threshold for PCA-ADI in concentric annuli.

# 3 RESULTS AND ANALYSIS

Fig. 1 shows our final six images of CQ Tau with VLT/NaCo, SPHERE/IRDIS, and SPHERE/IFS. Our ADI images are shown with one principal component to minimize geometric biases introduced by ADI (Milli et al. 2012), while the ASDI image is shown with three principal components to better emphasize each spiral feature.

## 3.1 Spiral arms and shadows

We detect spiral arms in all images. These are the most conspicuous in the $Q_\phi$ image, which is used for labelling and tracing of the spirals (Fig. 2). We detect a known prominent spiral arm (we have labelled MNRAS 515, 6109–6121 (2022))

S1) to the north-east in all data sets and is robust through all tested principal components for our ADI observations. We observe a radial extent of $0.2' - 0.4'$ for S1, consistent with Uyama et al. (2020). The NaCo data set features a kink in the L'-band spiral emission of the prominent spiral, however, we attribute this to the co-spatial, systematically bad horizontal pixels on the detector (Fig. A1, left) or a geometric bias of ADI and is unlikely to be a naturally occurring phenomena. A smaller tightly wound spiral (labelled S2) is found also to the north-east. The non-coronagraphic $Q_\phi$ image best shows the presence of sub-mm grains down to several au separation, in addition to a spiral in the north (labelled S3) and emission to the south (S4) recently reported in Safonov et al. (2022). There is a double-arc structure present to the SW (labelled A1 and A2) in our IRDIS ADI, PDI, and IFS A/SDI images. It is unclear whether this emission has been detected in NaCo data. Features S1 and A1 appear significantly more connected (as if a prolongation of S1) in our IFS data compared to $Q_\phi$, with the exception of an azimuthally extended discontinuity in the SE in the IFS data. However, this may be a post-processing artefact.

To characterize the spiral arms, we applied a spiral fitting routine to the $Q_\phi$ images. The scattering surface of the disc, derived through both the SED (Ubeira Gabellini et al. 2019) and CO rotation profiles (Bohn et al. 2022), was mapped using DISKMAP’s power-law function. The image was deprojected and $r^2$-scaled to account for radially decreasing stlocentric flux, knowing the Gaia distance to the system and the inclination ($35^\circ$) and PA ($235^\circ$) of the outer disc. We applied a low-pass filter (Gaussian kernel equal to half the FWHM) followed by a high-pass Laplacian filter to enhance disc features and improve the automatic detection of local radial maxima – the latter accomplished using the same technique as in Casassus et al. (2021). Fig. 2 shows the traces detected using this method as cyan crosses. We fit a logarithmic spiral to the trace for S1–S4 and A1–A2 using scipy optimize, which uses a Nelder-Mead minimization, to find the best-fitting parameters of the log spiral equation ($r = ae^{b\theta}$) with the pitch angle defined as $\tan^{-1}(1/b)$. As we have two functions for mapping the disc surface, contrasting pitch angles were measured for each spiral. For this reason, our final pitch angle uncertainty corresponds to the difference between the pitch angles of the two disc shapes and the two $Q_\phi$ images (S4 could only be traced in the non-coronagraphic image). The surface derived from the CO rotation profiles is close to geometrically flat in the inner $\lesssim 50$ au and is unlikely to occur in nature, however, it is meaningful to explore all our available functions. Our final fits are presented in Table 3 and visualized in Fig. 2. Fig. 3 shows the pitch angle variation found using consecutive arcs subtending $30^\circ$. Our errors are dominated by the uncertainty in the mapping of the disc surface.

We find pitch angles of $21.2^\circ$, $4.9^\circ$, $11.6^\circ$, and $25.2^\circ$ for S1–S4. Features A1 and A2 are more consistent with circular arcs. While S2 has a shallow pitch angle, the evolution of the trace is more supportive of a spiral (Fig. 3). S1 and S2 are opening radially outward whereas the opposite is seen for S3 and S4. Our deprojected $Q_\phi$ images and IFS image suggest that S1 and S3 may connect through A1 (Section 4.1).

The PDI and SDI data sets suggest a tentative shadow on the outer disc at two different position angles. A dimming is present at $\mathrm{PA}\approx140^\circ-160^\circ$ most clearly visible in $Q_\phi$, also reported in $^{12}$CO, $^{13}$CO, and $C^{18}O$ $J=2-1$ peak brightness temperature (Wölfer et al. 2021).

## 3.2 Companion detection limits

Apart from the point source at $\sim 2.2'$ separation (327 au, or $\sim 250$ au from the edge of the dust-continuum ring) that we determined to be
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Our reduced and annotated multiwavelength near-IR scattered light observations of CQ Tau. From left to right, top to bottom: VLT/NaCo ADI (L-band, NPC = 1), SPHERE/IRDIS ADI (Ks-band, NPC = 1), SPHERE/IRDIS PDI (J-band, r²-scaling), SPHERE/IFS SDI (J-band, NPC = 1), and ASDI (YJH-bands, NPC = 3). The regions attenuated by a coronagraph are masked, with a slightly larger mask on the IFS data to cover the bright signal at the edge of the coronagraphic mask. ADI data sets are displayed with one principal component to preserve features and limit geometric bias. All images use a linear colour scale and arbitrary cuts to best present all features. The disc semimajor and semiminor axis is overlaid on the $Q_\phi$ image. Grey patches represent the FWHM of the respective non-coronagraphic PSF.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Annotated versions of our coronagraphic (left) and non-coronagraphic (right) $Q_{\phi}$ images from Fig. 1, after deprojection, $r^2$-scaling, and low-pass and high-pass filtering. Best-fitting logarithmic spirals for features S1–S4 and A1–A2 are presented in blue, along with their traces (cyan crosses). Flux is shown in arbitrary units on a linear scale with limits chosen to best show all the features in the images.</figcaption></figure>

<table><caption>**Table 3.** Summary of logarithmic clock-wise spiral traces.</caption><thead><tr><td>Feature</td><td>Pitch angle (°)</td><td>Root PA (°)</td><td>ΔPA (°)</td></tr></thead><tbody><tr><td>S1</td><td>21.2 ± 3.1</td><td>140</td><td>100</td></tr><tr><td>S2</td><td>4.9 ± 3.8</td><td>70</td><td>50</td></tr><tr><td>S3</td><td>11.6 ± 1.4</td><td>58</td><td>120</td></tr><tr><td>S4</td><td>25.2 ± 2.1</td><td>270</td><td>90</td></tr><tr><td>A1</td><td>3.2 ± 3.8</td><td>260</td><td>60</td></tr><tr><td>A2</td><td>-0.3 ± 1.9</td><td>258</td><td>30</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Pitch angle of features S1–S4 and A1–A2. Each feature was traced using consecutive arcs subtending 30° to show the evolution of the pitch angle over the entire trace. Feature A2 subtended approximately 30° and thus only has one point. Uncertainties represent the difference between pitch angles using our two models for the disc scattering surface (Ubeira Gabellini et al. 2019; Bohn et al. 2022) and two $Q_{\phi}$ images.</figcaption></figure>

a background object (Appendix B), our final post-processed images did not contain any reliable point sources in $K_s$ or $L'$-band over a wide range of tested principal components (1–20). Fig. 4 shows the $5\sigma$ contrast curve, and Fig. 5 the mass sensitivity limits for our two ADI

data sets without subtracting the S1 spiral. The $3\sigma$ sensitivity at 0.3 and 0.6 was first determined using VIP's contrast curve function. The purpose of this step was to infer the approximate flux recovered of synthetic planets with S/N of 10 using full-frame PCA-ADI to 20 principal components. Knowing the approximate flux, we created six ADI cubes each containing a synthetic companion at 0.3 and 0.6, separated in PA by 60° forming a 'spiral' of uniquely positioned synthetic companions in each cube. The cubes were post-processed with full-frame and annular PCA-ADI techniques to 1–20 principal components, the latter using $3 \times$ FWHM pixel annuli and a rotation threshold of $1–3 \times$ FWHM for the inner and outer annuli. We compute the S/N for the synthetic companions in each cube to determine the optimal number of principal components, at each separation, ignoring negative side lobes extending azimuthally from the companions and accounting for small sample statistics as in Mawet et al. (2014). Our contrast at each separation was then the optimal contrast achieved with our range of principal components. The final $5\sigma$ contrast curve (Fig. 4) presents the optimal contrast obtained with our post-processing techniques and provides deeper contrast than previously reported (Asensio-Torres et al. 2021) in the noise-limited regime.

We computed the expected absolute magnitude of 0.5–104 $M_{Jup}$ giant planets and brown dwarfs using hot-start COND evolution models (Baraffe et al. 2003), with a linear interpolation between the 5 and 10 Myr magnitudes to 8.9 Myr. We use an apparent magnitude of 6.173 (2MASS, Cutri et al. 2003) and 5.043 (AllWISE, Cutri et al. 2021) in $K_s$ and $L'$-band, respectively, and the updated Gaia distance (Table 1) for CQ Tau to convert our contrast curve into absolute magnitude limits. Based on the analysis of Sanchis et al. (2020), we can exclude a significant effect of disc extinction at these distances from the star and these values of the potential planet mass. Our $K_s$-band mass sensitivity curve rules out giant planets above 4 $M_{Jup}$ outside the spiral region and 2–3 $M_{Jup}$ beyond 100 au to $5\sigma$ confidence (Fig. 5). The end of the NIR spiral region is characterized by a sharp drop in the contrast curve occurring at ~60 au. At that location, we enter the noise-limited regime and our contrast is no
<figure><img src="image_5.png" /><figcaption>**Figure 4.** 5σ contrast (in magnitudes) with respect to the star computed for our NaCo and IRDIS ADI data sets. The contrast for each radial separation was selected using the optimal contrast achieved by varying the number of principal components in Section 3.2. The ALC2 (IRDIS) and AGPM (NaCo) coronagraphs hide the inner 14 and 16 au, respectively.</figcaption></figure>

longer biased by the bright spirals. This curve does not rule out a brown-dwarf companion interior to the sub-mm ring, or a stellar companion in the gas cavity. The discrepancy between the contrast achieved with NaCo (*L'*-band) and SPHERE (*Ks*-band) is primarily due to the higher sky background in *L'*-band.

Uncertainty in the age of CQ Tau (Table 1) will influence our companion magnitudes (brighter at younger ages). For this reason, we re-computed our mass sensitivity using our constraints on the age of the system (6.3–11.7 Myrs), shown as the shaded region in Fig. 5. Since temporal subsampling of our NaCo data can reduce the S/N of a companion at large (>0.5) separation from the host star (Gomez Gonzalez et al. 2017), we also computed the *L'*-band mass sensitivity without binning with 20 principal components using both PCA-ADI flavours. No substantial change in contrast was detected, suggesting a median binning factor of 10 frames was not substantial enough to smear a candidate companion’s signal nor decorrelate the PSF. One can also infer that a single tested principal component will not provide optimal contrast, but instead sampling over a range of components will give the deepest contrast possible at each radial separation (as conducted in Fig. 5).

# 4 DISCUSSION

## 4.1 Where is the spiral-inducing companion?

Our non-detection of a wide orbit companion inducing the spiral arms in CQ Tau is consistent with the cases of HD 135344B and MWC 758, where a ~3 $M_{Jup}$ companion exterior to the spirals should have been detected to 5σ confidence (Maire et al. 2017; Reggiani et al. 2018; Boccaletti et al. 2021). Taking into account the multi-epoch astrometry of the point-source detected in Fig. B1, we rule out the possibility that this source is bound to the system and thus is not dynamically influencing the disc. While our *L'*-band images did not reach equivalent contrast to our *Ks*-band data, the Keck/NIRCC2 ADI data presented in Uyama et al. (2020) set companion mass limits of ~5 $M_{Jup}$ outside the prominent spiral to 0.9. The minimum mass required (on the order of several Jupiter masses) to excite such high-amplitude spiral arms, based on simulations (Dong et al. 2015; Juhász et al. 2015; Zhu et al. 2015; Dong & Fung 2017), is likely to be detectable in our IRDIS and NaCo data sets assuming

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Mass sensitivity curve for *Ks*-band and *L'*-band full frame PCA-ADI reductions determined with the hot-start COND evolution models (Baraffe et al. 2003) at 8.9 Myr. Shading around each curve corresponds to the age uncertainty. The dashed line represents an extrapolation of the models beyond 104 $M_{Jup}$. Red line gives the S1 spiral extent. Purple shading indicates the gas cavity. Vertical dotted line indicates the centre of the 1.3 mm dust continuum ring (Ubeira Gabellini et al. 2019, Wölfer et al. 2021).</figcaption></figure>

a ‘hot start’ (Fig. 5). Furthermore, CQ Tau does not feature a two-arm spiral pattern separated by ~180° expected for high-mass external companion (Fung & Dong 2015), although we recognize this conclusion could be influenced by illumination effects (Section 4.2). However, the recovered pitch angle for S3 is within expectations for inner spiral arms excited by a several Jupiter mass perturber (~10–15°; Dong et al. 2015). The pitch angle of the S1 feature in particular (34 ± 2° in Uyama et al. 2020, 21.2 ± 3.1° in this work due to our deprojection of the disc, Section 3.1) would be more favourably explained by an external perturber, however, can also form if an inner companion is on an eccentric orbit.

The mass sensitivity estimates are dependent upon a well-constrained age, distance, and apparent magnitude of the system. If the system is older than the 8.9 Myr used to compute the COND model magnitudes in Fig. 5, our 5σ detection limit will not reach as low mass and a companion may not be detected, and vice versa. We can, however, use the NICEST extinction map (Juvela & Montillaud 2016) based off the Two Micron All-Sky Survey to make the observation that CQ Tau is not spatially associated with the Taurus star-forming region, as defined in Lombardi, Lada & Alves (2010). Therefore, CQ Tau may be older than the 1–3 Myr Taurus cloud complex, in agreement with the aforementioned age range. A further consideration is the choice of post-formation luminosity of a companion, which can lead to significant magnitude discrepancies at the same age (Spiegel & Burrows 2012). High-contrast imaging is less supportive of cold-start models associated with core-accretion (Wang et al. 2018; Marleau, Mordasini & Kuiper 2019b) and favour warm- or hot-start initial conditions (although detection bias favours brighter, hot-start planets). The use of the hottest, most luminous models, AMES-DUSTY (Allard et al. 2001) would make our *Ks*-band images sensitive to lower-mass planets. On the other hand, Bern EXoplanet cooling tracks with COND atmospheres and warm-start initial conditions (Marleau et al. 2019a; Vigan et al. 2021 and recently applied to SPHERE data in Asensio-Torres et al. 2021) indicate a planet with mass on the order of ~8 $M_{Jup}$ may not be detected in our *Ks*-band images immediately outside the S1 spiral. CQ Tau is a relatively old disc allowing time for a slower planet formation, which may be at play. We therefore cannot entirely rule out an external perturber with our observations.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** The result of applying PCA-ADI to our coronagraphic $Q_{\phi}$ image with one principal component to test the effect of PCA-ADI on azimuthally extended signals. Flux is shown in arbitrary units on a linear scale with limits chosen to best show all the features in the image. Extended disc signals (i.e. spirals) have been strongly filtered.</figcaption></figure>

Searching for a planet or companion interior to the spirals is a logical next step. While Ubeira Gabellini et al. (2019) successfully fit the depletion of dust and gas tracers due to a 6–9 M<sub>Jup</sub> planet at 20 au with hydrodynamical models, spirals and non-axisymmetric features were not discussed. It is therefore natural to question whether an inner perturber could explain such features.

Our ADI data in Fig. 1 features several point-like sources and filtered spiral signals. However, PCA-ADI introduces geometric biases to azimuthally extended disc signals, transforming them into point-like sources. We therefore tested the application of PCA-ADI on our PDI data set. A fake ADI cube was built from rotating the $Q_{\phi}$ image using the same derotation angles as the IRDIS ADI observation, and subsequently post-processed with PCA-ADI to one principal component. The image, which we have labelled PDI PCA-ADI (Fig. 6), features the prominent single spiral arm S1 present in all our data sets, filtered versions of S2, S3, S4, A1, and potentially A2. The azimuthally extended emission from S2 and S3 has been most strongly filtered, such that the residual signals are analogous to those observed in our ADI data. Therefore, the distortion (or absence entirely) of these structures in ADI images is rather ADI filtering.

Our data sets suggest a single large-scale spiral arm (S1) in CQ Tau, accompanied by the existence of three smaller spirals (S2, S3, and S4) and a double-arc (A1 and A2). Fig. 7 shows the S1 spiral arm connects into a large-scale velocity perturbation in the $^{12}$CO data, suggesting that this is part of a broader spiral structure (Wölfer et al. 2021). A single large-scale or far-reaching spiral is not unusual, as seen in polarized light images of HD 100546, V1247 Ori, and HD 142527 (Ohta et al. 2016; de Boer et al. 2021; Avenhaus et al. 2014, respectively). In the case of HD 142527, all disc features have been qualitatively reproduced by the simulated interaction of the disc and an inclined, eccentric binary companion (Price et al. 2018). The presence of a planetary-mass companion or close-in binary remains ambiguous in the case of HD 100546 and V1247 Ori, however, candidates have been proposed (Kraus et al. 2017; Sissa et al. 2018; Fedele et al. 2021, Norfolk et al. 2022). Similarly, the binary system HD 34700 features a spiral of greater radial extent and pitch angle relative to a series of smaller spirals present in scattered light images (Monnier et al. 2019). Our mass limits support an analogous scenario,

<figure><img src="image_8.png" /><figcaption>**Figure 7.** $Q_{\phi}$ $r^2$-scaled contours overlaid on $^{12}$CO J = 2–1 rotation map, $v_0$, produced with the quadratic method in BETTERMOMENTS (Teague & Foreman-Mackey 2018). Regions with velocities not between 0 to 12 km s$^{-1}$ are masked. The arrow highlights the non-Keplerian velocity signature observed in the vicinity of the NIR spiral. CO data originally presented in Wölfer et al. (2021; data availability within).</figcaption></figure>

although CQ Tau features a smaller dust cavity compared to HD 142527 (~90 au) and lacks the extreme dust asymmetry of V1247 Ori, HD 142527, and HD 34700 (Pérez et al. 2014; Ohta et al. 2016; Benac et al. 2020, respectively).

Multi-epoch $I_c$-band ($\lambda_{\text{eff}} = 0.8$ µm) differential speckle polarimetry (DSP) observations by Safonov et al. (2022) spanning 2015–2021 detected a northern and southern spiral spanning $0^{\circ}.1-0^{\circ}.15$ from the star (which could account for some of the stellar variability, Dodin & Suslina 2021). The morphology recovered in our non-coronagraphic image at the same PA and separation is consistent with the DSP image, providing an independent confirmation this feature is real. While the authors suggested the northern spiral was S2, the feature is rather spatially associated with S3. A similar spiral to the south (PA $\approx 200^{\circ}$) is associated with emission best seen in our non-coronagraphic $Q_{\phi}$ image and is partially visible in our coronagraphic data sets (Fig. 1). Safonov et al. (2022) estimated the pattern speed of S3 to be $-0.2 \pm 1.1^{\circ}$ yr$^{-1}$, compared to $\sim 4.8^{\circ}$ yr$^{-1}$ if induced by a companion at 20 au in Ubeira Gabellini et al. (2019). If caused by a closer-in companion at shorter separation, one would expect a higher magnitude pattern speed of the spiral corresponding to the companion’s Keplerian velocity. The authors suggest S3 could instead be tracing the edge of an outer dust disc, however, our measured pitch angle is less supportive of this. This analysis would benefit from independent confirmation with multi-epoch SPHERE/IRDIS and Subaru/AO188 + HiCIAO polarimetric observations as only 135 d passed between our PDI data sets.

Fig. 2 and Table 3 show small (<4°) pitch angles for A1 and A2. Our images and spiral traces (Fig. 3) suggest that A1 could be connected to S1 and S3, forming part of the same spiral arm that has been visually segmented due to the illumination of the disc. If so, the trace of this feature would represent a pitch angle reversal which can only be explained with an inner companion on an eccentric orbit also predicted for MWC 758 (fig. 3 of Calcino et al. 2020). Such a perturber is favourable given our detection limits. In the case of A2, the arc-like appearance is reminiscent of HD 100453 (Benisty et al. 2017) and MWC 758 (Benisty et al. 2015; Reggiani et al. 2018). In both systems, this was attributed to the outer edge of the scattering
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Estimated total intensity of the r²-scaled coronagraphic $Q_\phi$ image after correcting for the degree of polarization with DISKMAP. Flux is displayed on a linear scale with an intensity cut from the 0.1–100 percentiles.</figcaption></figure>

surface on the bottom disc and located near the disc semiminor axis. However, in MWC 758, the arc is coincident with a dust trap (see fig. 1, Dong et al. 2018), which could suggest scattered light is probing the edge of a large dust concentration. We speculate the A2 feature may also be probing the edge of a large dust collection in CQ Tau considering its location along the disc semimajor (not minor) axis, near-zero pitch angle and proximity to the southern 1.3 mm continuum peak.

## 4.2 Are the shadows real? If so, what is causing them?

A drop in flux observed both towards the SE and west (this is most clear in $Q_\phi$ and when enhanced in Fig. 2) suggests shadowing on to the outer disc from the proposed inner disc ($r \approx 0.2$ au, $i \approx 48^\circ$, PA $\approx 106^\circ$ in Eisner et al. 2004 and illustrated in Safonov et al. 2022, $i \approx 29^\circ$, PA $\approx 140^\circ$ in Bohn et al. 2022). A drop in temperature co-located with the reduction in brightness was detected in CO isotopologues (Wölfer et al. 2021) to the south-east and north-west. The authors cite beam dilution as a possible mechanism due to its proximity to the disc semiminor axis, however, the re-detection of the southeast shadow in scattered light suggests this may be gas cooling in the presence of the shadow (e.g. HD 142527, Casassus et al. 2015, 2019). The large-scale scattered-light surface brightness asymmetry between the shadows closely resembles HD 143006, where two dark lanes give way to under brightness in half of the disc (Benisty et al. 2018). A warped disc model due to an inclined binary companion successfully reproduced the shadow, with a 10 $M_{Jup}$ circumbinary planet producing the dust cavity (Ballabio et al. 2021).

To test whether the shadows are real or a polarization phase function artefact, we show the estimated total intensity image in Fig. 8. As the shadows and large-scale brightness asymmetry still appear in the total intensity image, this implies the dimming is unlikely to be caused by the scattering properties of the dust grains in the disc.

Our observations are currently most supportive of a fairly axially symmetric shadow expected for a misaligned inner disc (Whitney et al. 2013), due to the presence of a massive companion in the cavity (Owen & Lai 2017). A recent study of the inner disc geometry using VLTI/GRAVITY (Bohn et al. 2022) found its inclination of

$29 \pm 3^\circ$ was a small deviation from that of the outer disc ($32 \pm 1^\circ$) determined with $^{12}$CO and $^{13}$CO line data, however, the position angle of the inner and outer disc major axis varied significantly ($140_{-9}^{+7}^\circ$ compared to $234 \pm 1^\circ$, respectively). The expected shadows cast on to the disc agree well with the darker regions and large-scale asymmetry. The discrepancy between the inner disc orientation presented in Eisner et al. (2004) could be due to its precession over the two epochs (e.g. Bate et al. 2000). Multi-epoch observations would further constrain the presence and movement of the inner disc. As the proposed inner disc precesses, the shadows should appear to ‘rock’ back and forth in azimuth over time (Nealon, Price & Pinte 2020). Irregularly changing patterns could indicate a variety of formation mechanisms (e.g. HD 135344B, Stolker et al. 2017) such as super-Keplerian transient shadows projected by an inner companion with a circumplanetary disc (e.g. GG Tau, Itoh et al. 2014; Keppler et al. 2020) or shadowing from the spiral itself (e.g. HD 34700, Montesinos et al. 2016).

## 4.3 What is causing the UX Ori variability?

UX Orionis stars are characterized by photometric variations of several magnitudes in the optical that coincide with increased extinction and polarization, hypothesized as being due to the disturbance of circumstellar dust by a companion or binary (Grinin et al. 1994, 2008). Shakhovskoj, Grinin & Rostopchina (2005) reported a 7800 d (~21.35 yr) cycle in the light-curve data for CQ Tau, with a secondary harmonic of 3900 d (~10.67 yr). When subtracted from the periodogram, a 1020 d (~2.79 yr) cycle is revealed (a cycle ratio close to that of UX Ori itself). American Association of Variable Star Observers data consisting of over 11 000 observations since 1946 presents the same approximate 10-year period. If the 21.35-yr cycle were to be caused by a companion orbiting the star with a combined mass of $1.5 M_\odot$, this equates to an orbital separation of 8.8 au. Interestingly, when following the theoretical evolution of co-planar binaries in protoplanetary discs, the semimajor axis of the cavity reaches $a_{cav} \approx 3.5a_{bin}$ (Artymowicz & Lubow 1994; Hirsh et al. 2020; Ragusa et al. 2020). In the case of CQ Tau with a gas cavity upper limit of 25 au, we also arrive at an orbital separation of ~8 au. However, models by Demidova, Sotnikova & Grinin (2010a) suggest the shorter cycle in UX Ori stars should correspond to a binary orbit (~2.2 au using the shorter cycle for CQ Tau) and its material streams, while the secondary cycle is caused by spiral waves (Sotnikova & Grinin 2007). Spectra from the Echelle SpectroPolarimetric Device for the Observation of Stars strongly supports the presence of a transient inhomogeneous dust cloud, comparable to the stellar radius, with crossing times corresponding to a distance of at least 0.9 au (Dodin & Suslina 2021).

Gaia DR3 finds a poor fit to the single-star model (Lindegren et al. 2021) indicated by the astrometric excess noise of 0.45 mas with a significance of 274, and dimensionless re-normalized unit weight error (RUWE) $\simeq 3.7$. The RUWE (ideally $= 1$ for a good fit to the single-star model), which is sensitive to unresolved binaries, is tentatively indicative of an unresolved object. The $15 \pm 2$ km s$^{-1}$ radial velocity of the star (in the Kinematic Local Standard of Rest frame, LSRK) is compatible with an M-dwarf companion on a ~2 au orbit as predicted by Demidova et al. (2010a), although measurements will be biased by accretion streams and the precession of, and obscuration by, the inclined inner disc. The presence of a circumstellar disc can also result in inflated RUWE values, prompting a higher cutoff for binaries in the presence of disc material. Fitton, Tofflemire & Kraus (2022) proposed such a cutoff at RUWE $= 2.5$, however, CQ Tau exceeds this value.
For comparison, these astrometric errors are larger than for PDS 70 (*d* = 112 pc) – a transition disc with two directly imaged protoplanets (e.g. Keppler et al. 2018; Benisty et al. 2021) – and smaller than those found for the stellar binary system GG Tau A (projected separation ~38 au, RUWE ≃ 6.4) at a comparable distance as CQ Tau. Most interestingly, the astrometric errors for CoKu Tau/4 (*d* = ~150 pc, astrometric excess noise of 0.8 mas with a significance of 293, RUWE ≃ 5.9), a near-equal-mass binary separated by 8 au (Ireland & Kraus 2008), are similar to those for CQ Tau. A binary companion would account for both the UX Ori-type variability and Gaia measurements, however, it is unclear at this time whether its separation would be 2 or 8 au.

## 5 CONCLUSION

We combined all available ESO high-contrast imaging data of CQ Tau to investigate the hypothesis of a companion dynamically influencing the disc. In addition, we developed a new, open-source NaCo data-reduction pipeline for reducing any coronagraphic NaCo data set with modern techniques, and implemented a PCA-based dark subtraction algorithm. Our key conclusions are as follows:

(i) We confirm the detection of an extended (0′.2 – 0′.4, Fig. 1) prominent spiral arm (Uyama et al. 2020), bright in Y–L'–band wavelengths co-spatial with the dust-continuum peak and large-scale velocity perturbation (Fig. 7), in addition to detecting three smaller spirals. Logarithmic spiral fitting resulted in a pitch angle of 21.2 ± 3.1°, 4.9 ± 3.8°, 11.6 ± 1.4°, and 25.2 ± 2.1°, respectively. Two arcs are detected, with one segment likely joining a larger spiral arm that appears discontinuous due to the illumination of the disc, while another is co-incident with the dust continuum (Section 4.1).

(ii) Companion detection limits appear to rule out a hot-start companion above 4 M<sub>Jup</sub> outside the spiral region in Ks-band and 2–3 M<sub>Jup</sub> beyond 100 au in 5σ confidence (Fig. 5). We detected NIR sources between ~ 0′.15 and ~ 0′.4 in our ADI data sets, however, we can not distinguish them from filtered disc features or residual speckles (Fig. 6).

(iii) Multi-epoch astrometry strongly suggests the point source at ~ 2′.2 separation is not bound to the system (nor a stellar flyby candidate), and CQ Tau should not be considered a binary based off this object (Appendix B).

(iv) Disc shadowing and illumination asymmetry is detected in our PDI data (Fig. 1). Two dark lanes are close to axially symmetric and likely caused by the misaligned inner disc best compared to HD 143006. Multi-epoch observations would constrain their change in azimuth.

(v) Our results (points i–iv), considered together with ALMA, GRAVITY, and DSP data (Ubeira Gabellini et al. 2019; Wölfer et al. 2021; Bohn et al. 2022; Safonov et al. 2022) support the hypothesis of an inner companion on an inclined orbit forming a misaligned inner disc (Section 4). Such a companion could be on an eccentric orbit driving the spiral pattern (Calcino et al. 2020; Debras, Baruteau & Donati 2021) and high Gaia RUWE, although dedicated models are required to re-produce the large pitch angle of the prominent spiral.

We suggest new SPH simulations that extend the work of Demidova et al. (2010a), Demidova, Grinin & Sotnikova (2010b), Demidova & Shevchenko (2015), and Calcino et al. (2020) to replicate the large-scale spiral structure. Upcoming Gaia astrometry and binary orbital fits could constrain the orbit of an unseen companion within the gas cavity. The narrow separation of the proposed companion makes it favourable for observations with the Keck Planet Imager and Characterizer in Vortex Fiber Nulling mode. Follow-up observations

with upcoming infrared high-contrast instruments on the Extremely Large Telescope could provide the contrast and resolution required to fully comprehend the mechanisms driving the substructures in CQ Tau.

## ACKNOWLEDGEMENTS

We thank the referee for their useful suggestions. IH acknowledges the support of a Research Training Program scholarship from the Australian government. IH, VC, DP, and CP acknowledge funding from the Australian Research Council via DP180104235 and FT170100040. VC also acknowledges funding from the Belgian F.R.S.-FNRS for financial support through a postdoctoral researcher fellowship. CT, GL, LT, and MGUG acknowledge funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement no. 823823 (RISE DUSTBUSTERS project). LT acknowledges support from the Italian Ministero dell'Istruzione, Università e Ricerca through the grant Progetti Premiali 2012-iALMA (CUP C52I13000140001), from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Ref no. 325594231FOR 2634/2 TE 1024/2-1, from the DFG Cluster of Excellence Origins (www.origins-cluster.de), and from the European Research Council (ERC) via the ERC Synergy Grant ECOGAL (grant 855130). MB acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant PROTOPLANETS no. 101002188). DF acknowledges the support of the Italian National Institute of Astrophysics (INAF) through the INAF Mainstream projects ARIEL and the 'Astrochemical Link between Circumstellar Disks and Planets', 'Proplanetary Disks Seen through the Eyes of New-generation Instruments' and by the PRIN-INAF 2019 Planetary Systems At Early Ages (PLATEA). This work has made use of the Multi-modal Australian ScienceS Imaging and Visualisation Environment (MASSIVE; www.massive.org.au). This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration (NASA). This work has made use of data from the European Space Agency (ESA) mission Gaia. Our data reduction and analysis made use of the PYTHON programming language and the ASTROPY (Astropy Collaboration et al. 2013, 2018), NUMPY (Harris et al. 2020), MATPLOTLIB (Hunter 2007), PHOTUTILS (Bradley et al. 2022), SCIKIT-IMAGE (van der Walt et al. 2014), SCIKIT-LEARN (Pedregosa et al. 2011), and SCIPY (Virtanen et al. 2020) packages.

## DATA AVAILABILITY

The reduced images will be made available on the Strasbourg astronomical Data Center (CDS) and through Monash University's *Bridges* research repository via doi:10.26180/20372814. Data reduction algorithms used in this article are available within the article.

## REFERENCES

Allard F., Hauschildt P. H., Alexander D. R., Tamanai A., Schweitzer A., 2001, ApJ, 556, 357

Amara A., Quanz S. P., 2012, MNRAS, 427, 948

Andrews S. M. et al., 2018, ApJ, 869, L41

Artymowicz P., 1993, ApJ, 419, 155

Artymowicz P., Lubow S. H., 1994, ApJ, 421, 651

Asensio-Torres R. et al., 2021, A&A, 652, A101
Astropy Collaboration, 2013, A&A, 558, A33

Astropy Collaboration, 2018, AJ, 156, 123

Avenhaus H., Quanz S. P., Schmid H. M., Meyer M. R., Garufi A., Wolf S., Dominik C., 2014, ApJ, 781, 87

Bae J., Zhu Z., 2018, ApJ, 859, 118

Ballabio G., Nealon R., Alexander R. D., Cuello N., Pinte C., Price D. J., 2021, MNRAS, 504, 888

Banzatti A., Testi L., Isella A., Natta A., Neri R., Wilner D. J., 2011, A&A, 525, A12

Baraffe I., Chabrier G., Barman T. S., Allard F., Hauschildt P. H., 2003, A&A, 402, 701

Baruteau C. et al., 2019, MNRAS, 486, 304

Baruteau C., Wafflard-Fernandez G., Le Gal R., Debras F., Carmona A., Fuente A., Rivière-Marichalar P., 2021, MNRAS, 505, 359

Bate M. R., Bonnell I. A., Clarke C. J., Lubow S. H., Ogilvie G. I., Pringle J. E., Tout C. A., 2000, MNRAS, 317, 773

Benac P. et al., 2020, ApJ, 905, 120

Benisty M. et al., 2015, A&A, 578, L6

Benisty M. et al., 2017, A&A, 597, A42

Benisty M. et al., 2018, A&A, 619, A171

Benisty M. et al., 2021, ApJ, 916, L2

Bertrang G. H.-M., Avenhaus H., Casassus S., Montesinos M., Kirschslager F., Perez S., Cieza L., Wolf S., 2018, MNRAS, 474, 5105

Beuzit J.-L. et al., 2008, in McLean I. S., Casali M. M., ed., Proc. SPIE Conf. Ser. Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II. SPIE, Bellingham, p. 701418

Bitsch B., Trifonov T., Izidoro A., 2020, A&A, 643, A66

Boccaletti A. et al., 2021, A&A, 652, L8

Bohn A. J. et al., 2022, A&A, 658, A183

Bradley L. et al., 2022, Zenodo, astropy/photutils: v1.5.0

Calcino J., Christiaens V., Price D. J., Pinte C., Davis T. M., van der Marel N., Cuello N., 2020, MNRAS, 498, 639

Calcino J., Price D. J., Pinte C., van der Marel N., Ragusa E., Dipierro G., Cuello N., Christiaens V., 2019, MNRAS, 490, 2579

Canovas H. et al., 2015, ApJ, 805, 21

Casassus S. et al., 2013, Nature, 493, 191

Casassus S. et al., 2015, ApJ, 812, 126

Casassus S. et al., 2021, MNRAS, 507, 3789

Casassus S., Pérez S., Osses A., Marino S., 2019, MNRAS, 486, L58

Chapillon E., Guilloteau S., Dutrey A., Piétu V., 2008, A&A, 488, 565

Christiaens V. et al., 2021, MNRAS, 502, 6117

Claudi R. U. et al., 2008, in McLean I. S., Casali M. M., ed., Proc. SPIE Conf. Ser. Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II. SPIE, Bellingham, p. 70143E

Cutri R. M. et al., 2003, VizieR Online Data Catalog, II/246

Cutri R. M. et al., 2021, VizieR Online Data Catalog, II/328

de Boer J. et al., 2020, A&A, 633, A63

de Boer J. et al., 2021, A&A, 649, A25

Debras F., Baruteau C., Donati J.-F., 2021, MNRAS, 500, 1621

Demidova T. V., Grinin V. P., Sotnikova N. Y., 2010b, Astron. Lett., 36, 498

Demidova T. V., Shevchenko I. I., 2015, ApJ, 805, 38

Demidova T. V., Sotnikova N. Y., Grinin V. P., 2010a, Astron. Lett., 36, 422

Dodin A. V., Suslina E. A., 2021, MNRAS, 503, 5704

Dohlen K. et al., 2008, in McLean I. S., Casali M. M., ed., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II. SPIE, Bellingham, p. 70143L

Donehew B., Brittain S., 2011, AJ, 141, 46

Dong R. et al., 2018, ApJ, 860, 124

Dong R., Fung J., 2017, ApJ, 835, 38

Dong R., Zhu Z., Rafikov R. R., Stone J. M., 2015, ApJ, 809, L5

Ducati J. R., Bevilacqua C. M., Rembold S. B., Ribeiro D., 2001, ApJ, 558, 309

Eisner J. A., Lane B. F., Hillenbrand L. A., Akeson R. L., Sargent A. I., 2004, ApJ, 613, 1049

Fedele D., Toci C., Maud L., Lodato G., 2021, A&A, 651, A90

Fitton S., Tofflemire B. M., Kraus A. L., 2022, Res. Notes Am. Astron. Soc., 6, 18

Follette K. B. et al., 2017, AJ, 153, 264

Francis L., van der Marel N., 2020, ApJ, 892, 111

Fung J., Dong R., 2015, ApJ, 815, L21

Gaia Collaboration, 2021, A&A, 649, A1

Garufi A. et al., 2013, A&A, 560, A105

Goldreich P., Tremaine S., 1979, ApJ, 233, 857

Goldreich P., Tremaine S., 1980, ApJ, 241, 425

Gomez Gonzalez C. A. et al., 2017, AJ, 154, 7

Gonzalez J.-F. et al., 2020, MNRAS, 499, 3837

Grinin V. P., Barsunova O. Y., Shugarov S. Y., Kroll P., Sergeev S. G., 2008, Astrophysics, 51, 1

Grinin V. P., The P. S., de Winter D., Giampapa M., Rostopchina A. N., Tambovseva L. V., van den Ancker M. E., 1994, A&A, 292, 165

Guzmán-Díaz J. et al., 2021, A&A, 650, A182

Harris C. R. et al., 2020, Nature, 585, 357

Herbig G. H., 1960, ApJ, 131, 632

Hirsh K., Price D. J., Gonzalez J.-F., Ubeira-Gabellini M. G., Ragusa E., 2020, MNRAS, 498, 2936

Hunter J. D., 2007, Comput Sci. Eng., 9, 90

Hunziker S., Quanz S. P., Amara A., Meyer M. R., 2018, A&A, 611, A23

Ireland M. J., Kraus A. L., 2008, ApJ, 678, L59

Itoh Y. et al., 2014, Res. Astron. Astrophys., 14, 1438

Juhász A., Benisty M., Pohl A., Dullemond C. P., Dominik C., Paardekooper S. J., 2015, MNRAS, 451, 1147

Juvela M., Montillaud J., 2016, A&A, 585, A38

Keppler M. et al., 2018, A&A, 617, A44

Keppler M. et al., 2020, A&A, 639, A62

Kraus S. et al., 2017, ApJ, 848, L11

Langlois M. et al., 2013, in Esposito S., Fini L., ed., Proceedings of the Third AO4ELT Conference, Firenze, p. 63

Launhardt R. et al., 2020, A&A, 635, A162

Lenzen R. et al., 2003, in Iye M., Moorwood A. F. M., ed., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes. SPIE, Bellingham, p. 944

Lindegren L. et al., 2021, A&A, 649, A2

Lombardi M., Lada C. J., Alves J., 2010, A&A, 512, A67

Long F. et al., 2018, ApJ, 869, 17

Maire A. L. et al., 2017, A&A, 601, A134

Maire A.-L. et al., 2016, in Evans C. J., Simard L., Takami H., ed., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI. SPIE, Bellingham, p. 990834

Mannings V., Sargent A. I., 1997, ApJ, 490, 792

Marleau G.-D., Coleman G. A. L., Leleu A., Mordasini C., 2019a, A&A, 624, A20

Marleau G.-D., Mordasini C., Kuiper R., 2019b, ApJ, 881, 144

Marois C., Lafrenière D., Doyon R., Macintosh B., Nadeau D., 2006a, ApJ, 641, 556

Marois C., Lafrenière D., Macintosh B., Doyon R., 2006b, ApJ, 647, 612

Mawet D. et al., 2013, The Messenger, 152, 8

Mawet D. et al., 2014, ApJ, 792, 97

Mawet D., Riaud P., Absil O., Surdej J., 2005, ApJ, 633, 1191

Mendigutía I., Mora A., Montesinos B., Eiroa C., Meeus G., Merín B., Oudmaijer R. D., 2012, A&A, 543, A59

Milli J., Mouillet D., Lagrange A.-M., Boccaletti A., Mawet D., Chauvin G., Bonnefoy M., 2012, A&A, 545, A111

Monnier J. D. et al., 2019, ApJ, 872, 122

Montesinos M., Perez S., Casassus S., Marino S., Cuadra J., Christiaens V., 2016, ApJ, 823, L8

Muley D., Fung J., van der Marel N., 2019, ApJ, 879, L2

Muto T. et al., 2012, ApJ, 748, L22

Natta A., Grinin V., Mannings V., 2000, in Mannings V., Boss A. P., Russell S. S., ed., Protostars and Planets IV. University of Arizona Press, p. 559

Natta A., Prusti T., Neri R., Wooden D., Grinin V. P., Mannings V., 2001, A&A, 371, 186

Nealon R., Price D. J., Pinte C., 2020, MNRAS, 493, L143

Norfolk B. J. et al., 2022, preprint (arXiv:2208.02542)
Ogilvie G. I., Lubow S. H., 2002, *MNRAS*, 330, 950

Ohta Y. et al., 2016, *PASJ*, 68, 53

Owen J. E., Lai D., 2017, *MNRAS*, 469, 2834

Pecaut M. J., Mamajek E. E., 2013, *ApJS*, 208, 9

Pedregosa F. et al., 2011, J. Mach. Learn. Res., 12, 2825

Pérez L. M., Isella A., Carpenter J. M., Chandler C. J., 2014, *ApJ*, 783, L13

Pinilla P. et al., 2018, *ApJ*, 859, 32

Price D. J. et al., 2018, *MNRAS*, 477, 1270

Ragusa E., Alexander R., Calcino J., Hirsh K., Price D. J., 2020, *MNRAS*, 499, 3362

Reggiani M. et al., 2018, *A&A*, 611, A74

Rosotti G. P. et al., 2020, *MNRAS*, 491, 1335

Rousset G. et al., 2003, in Wizinowich P. L., Bonaccini D., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 4839, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, SPIE, Bellingham, p. 140

Safonov B. S., Strakhov I. A., Goliguzova M. V., Voziakova O. V., 2022, *AJ*, 163, 31

Sanchis E., Picogna G., Ercolano B., Testi L., Rosotti G., 2020, *MNRAS*, 492, 3440

Shakhovskoj D. N., Grinin V. P., Rostopchina A. N., 2005, *Astrophysics*, 48, 135

Sissa E. et al., 2018, *A&A*, 619, A160

Sotnikova N. Y., Grinin V. P., 2007, *Astron. Lett.*, 33, 594

Spiegel D. S., Burrows A., 2012, *ApJ*, 745, 174

Stolker T. et al., 2017, *ApJ*, 849, 143

Stolker T., Dominik C., Min M., Garufi A., Mulders G. D., Avenhaus H., 2016, *A&A*, 596, A70

Teague R., Foreman-Mackey D., 2018, Res. Notes Am. Astron. Soc., 2, 173

Testi L., Natta A., Shepherd D. S., Wilner D. J., 2001, *ApJ*, 554, 1087

Testi L., Natta A., Shepherd D. S., Wilner D. J., 2003, *A&A*, 403, 323

Toci C., Lodato G., Christiaens V., Fedele D., Pinte C., Price D. J., Testi L., 2020a, *MNRAS*, 499, 2015

Toci C., Lodato G., Fedele D., Testi L., Pinte C., 2020b, *ApJ*, 888, L4

Tripathi A., Andrews S. M., Birnstiel T., Wilner D. J., 2017, *ApJ*, 845, 44

Trotta F., Testi L., Natta A., Isella A., Ricci L., 2013, *A&A*, 558, A64

Ubeira Gabellini M. G. et al., 2019, *MNRAS*, 486, 4638

Ubeira-Gabellini M. G., Christiaens V., Lodato G., Ancker M. v. d., Fedele D., Manara C. F., Price D. J., 2020, *ApJ*, 890, L8

Uyama T. et al., 2020, *AJ*, 159, 118

van der Marel N. et al., 2013, *Science*, 340, 1199

van der Marel N. et al., 2021, *AJ*, 161, 33

van der Marel N., van Dishoeck E. F., Bruderer S., Andrews S. M., Pontoppidan K. M., Herczeg G. J., van Kempen T., Miotello A., 2016, *A&A*, 585, A58

van der Walt S. et al., 2014, *Peer J.*, 2, e453

van Holstein R. G. et al., 2020a, Astrophysics Source Code Library, record ascl:2004.015

van Holstein R. G. et al., 2020b, *A&A*, 633, A64

Vigan A. et al., 2021, *A&A*, 651, A72

Villenave M. et al., 2019, *A&A*, 624, A7

Vioque M., Oudmaijer R. D., Baines D., Mendigutía I., Pérez-Martínez R., 2018, *A&A*, 620, A128

Virtanen P. et al., 2020, *Nature Methods*, 17, 261

Wang J. J. et al., 2018, *AJ*, 156, 192

Whitney B. A., Robitaille T. P., Bjorkman J. E., Dong R., Wolff M. J., Wood K., Honor J., 2013, *ApJS*, 207, 30

Wölfer L. et al., 2021, *A&A*, 648, A19

Zhu Z., Dong R., Stone J. M., Rafikov R. R., 2015, *ApJ*, 813, 88

Zhu Z., Zhang R. M., 2022, *MNRAS*, 510, 3986

<figure><img src="image_10.png" /><figcaption>**Figure A1.** Median and PCA-based dark subtraction of a NaCo science cube captured in the middle of the CQ Tau observation. Both are produced using the NaCo pipeline and are displayed on a linear scale with an intensity cut spanning 0.1–99.9 percentiles. The systematically bad quadrant has already been corrected as described in Section 2.3.</figcaption></figure>

model and subtract the varying (both spatially and temporally) bias level, dark current, and thermal radiation from the instrument captured in the dark frames. In contrast to PCA-based sky subtraction (Hunziker et al. 2018), the aim of PCA in this context is to best correct the read-out grid pattern and systematic erroneous horizontal bars on the Aladdin 2 detector (Fig. A1, left) common in the dark and science frames. All dark frames are first cropped to the size of the science frames and saved into a cube for use as the PCA library. The systematic bad quadrant, AGPM centre and its shadow cast on the detector were masked to prevent potentially erroneous pixels and the stellar PSF being included in the PCA library.

We begin an iterative loop over the science cubes to best correct the systematic horizontal bars: (i) add the difference of the median science pixel value and the median dark pixel value to the science cube; (ii) provide this science cube to VIP's *cube_subtract_sky_pca* routine with the median dark cube as the PCA library to perform PCA dark subtraction with one principal component; (iii) restore the science image by adding the difference calculated in (i) to the PCA dark subtracted image; (iv) the standard deviation of pixel values in a subframe around the spurious horizontal bars is calculated; (v) loop over steps (i–iv) and search for the optimal offset [added in step (i)] that minimizes the standard deviation [measured in step (iv)] using *scipy.minimize* (Nelder–Mead solver) that minimizes the standard deviation metric; (vi) end the minimization after 100 iterations or when the function converges on a minimum. We repeat this process for the sky and flat cubes.

Fig. A1 presents a significant reduction of background bias in the corrected science frames when using PCA-based dark subtraction compared to a median subtraction.

## APPENDIX B: BACKGROUND SOURCE

We resolved a point source at ~2''2 (327 au) separation from CQ Tau in our ADI and PDI data sets. Fig. B1 (right) shows the relative astrometry for the point source with respect to the star in the IRDIS ADI (Fig. B1, left), NaCo ADI, and both IRDIS PDI data sets. The predicted track for a static background object is reversed through time, starting with the 2018 NaCo measurement.

In order to provide maximum time baseline, we considered archival NASA/ESO *Hubble* Space Telescope (HST) observations from 2000 February 02 (Program ID 8216) captured by the Wide-Field and Planetary Camera 2 (PC2) instrument with the F675W ($\lambda_c$ = 0.6717 µm) and F814W ($\lambda_c$ = 0.8012 µm) filters. A second HST observation captured by the NICMOS instrument (camera 2) on 2005 January 08 with the F110W ($\lambda_c$ = 1.025 µm) filter was also

## APPENDIX A: PCA DARK SUBTRACTION WITH NACO

We adapted the PCA-based sky subtraction routine in VIP (*cube_subtract_sky_pca*, Gomez Gonzalez et al. 2017) to
<figure><img src="image_11.png" /><figcaption>**Figure B1.** Left: point-source detection in the 2016 IRDIS *Ks*-band data set (north is up, east is to the left). Brightness is shown on a log scale spanning 0–99.9 percentile in intensity. Right: empty points correspond to the expected astrometry for a background object at each epoch, filled points are the measured astrometry. The astrometry is consistent with that of a background object.</figcaption></figure>

considered. We conservatively assumed the uncertainty on each HST astrometric point to be the pixel scale of the corresponding instrument (0.0455 pixel⁻¹ for PC2, 0.075 pixel⁻¹ for NIC2) and assumed the star was at the intersection of the main diffraction artefacts in the PC2 dataset. The source was easily detected by eye in all data sets and its location was inferred by fitting a 2D Gaussian to a subframe encompassing the point source. In Fig. B1 (right), we compare the detected motion of the source with the expected trajectory for a background star, considering the proper motion of CQ Tau from Gaia (RA = 2.98 ± 0.06 mas yr⁻¹, Dec. = −26.36 ± 0.04 mas yr⁻¹; Gaia Collaboration 2021).

To explore whether this object previously interacted with the disc around CQ Tau, we calculated its *J*–*Ks* colour and inferred the

spectral-type using existing pre-main sequence (Pecaut & Mamajek 2013) and main sequence (Ducati et al. 2001) colours to be consistent with an M-type star. Knowing the contrast, apparent magnitude of CQ Tau in *Ks*-band (Section 3.2) and the expected absolute magnitude for a young M-type star, we infer its distance to be ≳2000 pc (compared to 149 pc for CQ Tau). As the astrometry of the point source is consistent with a background star with a small non-zero proper motion itself, we therefore rule out the possibility that this is a bound object or stellar flyby.

This paper has been typeset from a *TEX/LATEX* file prepared by the author.