A highly non-Keplerian protoplanetary disc: Spiral structure in the gas disc of CQ Tau

In the past years, high angular resolution observations have revealed that circumstellar discs appear in a variety of shapes with diverse substructures being ubiquitous. This has given rise to the question of whether these substructures are triggered by planet-disc interactions. Besides direct imaging, one of the most promising methods to distinguish between different disc shaping mechanisms is to study the kinematics of the gas disc. In particular, the deviations of the rotation profile from Keplerian velocity can be used to probe perturbations in the gas pressure profile that may be caused by embedded planets. In this paper we aim to analyze the gas brightness temperature and kinematics of the transitional disc around the CQ Tau star in order to resolve and characterize substructure in the gas, caused by possible perturbers. For our analysis we use spatially resolved ALMA observations of 12CO, 13CO and C18O (J=2-1). We further extract robust line centroids for each channel map and fit a number of Keplerian disc models to the velocity field. The gas kinematics of the CQ Tau disc present non-Keplerian features, showing bent and twisted iso-velocity curves in 12CO and 13CO. Significant spiral structures are detected between 10-180 au in both the brightness temperature and the rotation velocity of 12CO after subtraction of an azimuthally symmetric model, which may be tracing planet-disc interactions with an embedded planet or low-mass companion. We identify three spirals, two in the brightness temperature and one in the velocity residuals, spanning a large azimuth and radial extent. The brightness temperature spirals are morphologically connected to spirals observed in NIR scattered light in the same disc, indicating a common origin. Together with the observed large dust and gas cavity, the spirals support the hypothesis of a massive embedded companion in the CQ Tau disc.


Introduction
Due to angular momentum conservation, circumstellar discs are the natural outcome of the star formation process when infalling material from a molecular cloud core is channeled towards the newly formed central star. These accretion discs composed of gas and dust represent the nurseries of planetary systems. They evolve and ultimately disperse while giving birth to various objects with the evolutionary processes significantly influencing the ongoing planet formation. At the same time, the forming planets will backreact on the disc and affect its evolution and structure, resulting in a highly coupled and complex problem.
High angular resolution observations indeed show that circumstellar discs are commonly marked by a variety of substruc-tures in the gas and especially the dust, such as gaps, rings or even cavities, as well as spiral arms and azimuthal asymmetries (e.g. van der Marel et al. 2013;Casassus 2016;Andrews et al. 2018;Cazzoletti et al. 2018;Long et al. 2018;Andrews 2020). Such substructures might be caused by embedded (proto-) planets (e.g. Lin & Papaloizou 1986;Zhang et al. 2018;Lodato et al. 2019), suggesting that planet formation occurs already in early evolutionary stages. However, there exist other mechanisms such as the magnetorotational instability (MRI) (e.g. Flock et al. 2015Flock et al. , 2017Riols & Lesur 2019), zonal flows (e.g. Uribe et al. 2015), the compositional baroclinic instability (Klahr & Bodenheimer 2004) or gravitational instability (Kratter & Lodato 2016) that could also account for the observations. One way to distinguish between the different scenarios and to understand possible planet-disc interactions is to directly image a young planet in its environment (e.g. Keppler et al. 2018;Wagner et al. 2018). Since this technique is only feasible for very few, massive objects, that are not affected by dust extinction (Sanchis et al. 2020), another promising method is to look for perturbations that are induced in the velocity field of the rotating gas. In this context, studying the gas component can help to access the different dynamical processes that are shaping the disc and reveal a number of previously undetected substructures. The density structure of dust grains that are typically probed by ALMA observations is determined by the gas dynamics. It is thus of paramount importance to directly access and characterize the gas kinematics to distinguish between various scenarios.
Different and complementary image analysis techniques to probe disc kinematics are being developed. For a geometrically thick disc around a single star that is both in radial and vertical hydrostatic equilibrium, the gas rotation velocity v rot is given by with r being the cylindrical radius, M * the mass of the star, ρ gas the gas density and ∂P/∂r representing the radial pressure gradient. Identifying deviations from Keplerian rotation can therefore be used to probe the local pressure gradient and to characterize the shape of the perturbation. Additional deviations may arise for a massive disc, due to its gravitational field. This technique has recently been used by Teague et al. (2018aTeague et al. ( , 2019a to constrain the gas surface density profile of the HD 163296 disc, leading to the kinematical detection of two embedded Jupiter-mass planets as well as significant meridional flows. In addition, Teague et al. (2018b) report a vertical dependence on the pressure maxima, studying the gas kinematics of AS 209. The deviations from Keplerian rotation are further used by Rosotti et al. (2020b) to measure the gas-dust coupling as well as the width of gas pressure bumps. Pinte et al. (2018Pinte et al. ( , 2019 detect 'kink-' features in the iso-velocity contours of HD 163296 and HD 97048 data respectively, consistent with a Jupiter-mass planet (∼ 2 M Jup ). Tentative detections of such azimuthally located features have also been found in a few discs of the ALMA DSHARP large program Pinte et al. 2020), but more data are needed to confirm the robustness of such claim. Similarly, a possible signature for an embedded planet in the HD 100546 disc is presented by Casassus & Pérez (2019) who reveal a Doppler-flip in the residual kinematical structure after subtracting a Keplerian best-fit model, as expected from a planet-disc interaction model (e.g. Pérez et al. 2015Pérez et al. , 2018. One type of discs, the so-called transition discs, are of particular interest, since they show evidence for dust (and gas) depleted inner regions (e.g. Strom et al. 1989;Ercolano & Pascucci 2017). Sometimes treated as being in a transition phase from an optically thick disc to disc dispersal, transition discs may enable to probe various mechanisms that play a role during disc evolution and represent excellent candidates to catch planet formation in action. Detecting a planet in discs with cavities may link planet formation with fully formed planetary systems and put constraints on the formation processes and timescales.
In this work we study the transitional disc around CQ Tau. Following Ubeira Gabellini et al. (2019), who focused on the radial profiles, characterization of the present dust and gas cavity and possible formation mechanisms, we analyze the gas component of the CQ Tau disc both in terms of its velocity and temperature structure, finding significant spiral structures. The paper is organised as follows: In Sect. 2 we describe the observations and data reduction, whereas the observational results are presented in Sect. 3. A description and analysis of the spiral structure is shown in Sect. 4 alongside with the method used to extract and model the gas kinematics. Our results are discussed in Sect. 5 and summarised in Sect. 6.

Target
The variable star CQ Tau (UX Ori class) is a YSO of spectral type F2 located in the Taurus star forming region at a distance of ∼ 162 pc (Gaia Collaboration et al. 2018 Lopez et al. 2006;Ubeira Gabellini et al. 2019) has an estimated age of ∼ 10 Myr and is surrounded by a massive circumstellar disc (Natta et al. 2000), which is found to have a high accretion rate of the order of 10 −8 − 10 −7 M yr −1 (Donehew & Brittain 2011;Mendigutía et al. 2012).
The disc around CQ Tau represents one of the first discs whose mm continuum was observed with different instruments (e.g. OVRO interferometer (Mannings & Sargent 1997); PdBI (Natta et al. 2000); VLA (Testi et al. 2001)) in order to constrain its dust properties. An analysis of the spectral slope at mm-wavelengths reveals that dust grains have grown to larger sizes than the typical ISM size (Testi et al. 2001(Testi et al. , 2003Chapillon et al. 2008). The average dust opacity coefficient is constrained by Banzatti et al. (2011) using VLA (1.3-3.6 mm), PdBI (2.7-1.3 mm) and SMA (0.87 mm) observations, probing significant grain growth in the disc with up to cm-sized grains. Trotta et al. (2013) further find that larger grains are present in the inner disc with respect to the outer disc, indicating a variation of grain growth with radius.
Subsequent high resolution gas and dust observations have revealed the CQ Tau disc to be a transition disc with an inner cavity. Tripathi et al. (2017) detect a gap in the 880 µm continuum emission of new and archival SMA data and Pinilla et al. (2018) report a dust cavity of ∼ 46 au in ALMA observations of the mm continuum, fitting the intensity profile in the visibility plane. Ubeira Gabellini et al. (2019) present recent ALMA observations, confirming a large cavity of 53 au radius (peak of the Gaussian dust ring) in the 1.3 mm continuum as well as a smaller gas cavity of 20 au in the 13 CO and C 18 O emission, fitting the surface density profiles. The authors performed 3D hydrodynamical simulations which suggest a hidden planet of several M J located at ∼ 20 au as a possible cause for the observed gas and dust depleted regions. Even though such a planet could not be detected in combined Keck/NIRC2 and Subaru/AO188+HiCIAO observations of CQ Tau (Uyama et al. 2020), due to a lack of contrast (compare their Figure 3), the data reveal the presence of a small spiral seen in small dust grains on scales of 30-60 au, that might be induced by a companion candidate.

Data reduction
We present 1.3 mm ALMA observations of the CQ Tau system in band 6, combining datasets from cycle 2, 4 and 5 (2013.1.00498.S, PI: 2016.A.00026.S, 2017, previously presented at a lower angular resolution in Ubeira Gabellini et al. (2019), and detailed in Table A.1. These three projects have different antenna configuration, but they share a similar spectral setup. The longest baseline from the first project extends to 1091 m, while for the latest two it is increased to 3700 m and 8500 m respectively, thus enhancing the spatial resolution. In all three observations the ALMA correlator was configured to observe the 1.3 mm dust continuum emission, as well as the molecular lines 12 CO J = 2 − 1, 13 CO J = 2 − 1, and C 18 O J = 2 − 1. After applying ALMA standard pipeline calibration, we follow a similar processing as for the DSHARP data calibration , using CASA 5.4.1. We start by flagging the channels located at ±25 km s −1 from each spectral line, and average the remaining channels to 125 MHz width channels, which are combined with the data from the continuum spectral windows. As a next step, we align the dust continuum emission and check by comparison the flux calibration of each individual execution, to ensure they all have the same flux.
To enhance the signal-to-noise ratio, self-calibration is performed in two stages. First, we self-calibrate the shorter baseline dataset, corresponding to the Cycle 2 observations (dataset #1 in Table A.1), by applying three steps of phase-only calibration using solution intervals of 300, 120 and 30 s, and one step of amplitude calibration using the whole observation time range as a solution interval. This self-calibrated short baseline dataset is then combined with the extended baseline datasets from the observations obtained in Cycles 4 and 5 (datasets #2 and 3), and four phase calibrations with solution intervals of 900, 360, 150, 90 s, as well as one amplitude calibration with solution interval of 360 s.
All the dust continuum emission calibration steps, including the centroid shifting and self-calibration tables, are then applied to the molecular line emission channels. The continuum emission is subtracted using the uvcontsub task, and image cubes are generated for each isotopologue (compare channel maps in Appendix B) using a robust parameter of 0.6 and Keplerian masking. This value was found to give the best trade off between spatial resolution and sensitivity. The Keplerian mask is calculated with the package keplerian_mask 1 , using an inclination of 35 • , position angle of 235 • , distance of 162 pc, stellar mass of 1.54 M and systemic velocity of 6.17 km s −1 . These values were chosen after some initial fits for the gas kinematics, explained later in Sect. 4.2. We further choose an inner and outer radius of 0 and 2 respectively and convolve with the beam rescaled by 1.5 times its size. Some important characteristics of the data are given in Table 1.

Integrated and peak brightness temperature maps
In Fig. 1 we show the velocity integrated intensity maps (left panels) along with the peak brightness temperature maps (right panels) for the three different CO isotopologues. The underlying channel maps are presented in Appendix B. The 12 CO integrated intensity is overlaid by the contours of the 1.3 mm continuum. To compute the integrated as well as the peak intensity maps 1 https://github.com/richteague/keplerian_mask we have used the bettermoments code described in Teague & Foreman-Mackey (2018) and then converted from flux density units to units of Kelvin with the Planck law. In addition, the Keplerian mask described in Sect. 2.2 is applied in the velocity integrated intensity maps to enhance the signal-to-noise. This results in a peak SNR of 40 for the 12 CO, 20 for the 13 CO and 14 for the C 18 O velocity integrated intensity as well as 29 for the 12 CO, 19 for the 13 CO and 18 for the C 18 O peak intensity.
While the optically thick 12 CO data do not trace any cavities in the gas distribution, a significant gas cavity is observed in the inner disc region as seen in the optically thinner 13 CO and C 18 O emission, shown also by Ubeira Gabellini et al. (2019) at a lower angular resolution. In addition to Ubeira Gabellini et al.
(2019) we clearly note a drop of the peak intensity of all isotopologues and integrated intensity of C 18 O by roughly 35-60 % (with respect to the peak) in the north-west and south-east part of the disc along the minor axis (symmetric pair of dips), that becomes especially prominent in the C 18 O data and co-locates with a possible under-brightness in NIR scattered light (see Fig.  4 of Uyama et al. 2020). Similar to the NIR under-brightness, the drop of peak intensity appears to be more pronounced (∼ 7-12 % deeper dip) in the south-east side of the disc (compare also Sect. 3.2).
In the past such symmetric features have been linked to the presence of a misaligned inner disc, casting a shadow over the outer disc (Marino et al. 2015;Facchini et al. 2018;Casassus & Pérez 2019). However, the co-location of this under-brightness with the minor axis of the disc suggests caution, since beam dilution can lead to artificial azimuthal features due to the low compact emission of the line central channels. To test for this effect we have used the DALI (Dust And LInes Bruderer et al. 2012;Bruderer 2013)  1 of Appendix C for 12 CO. Two clear dips appear along the minor axis, with the results being similar for the more optically thin lines. Thus, beam dilution can mostly account for the strong under-brightness seen in the data. Similar to the data, the upper side of the disc (including the north-west dip) appears slightly brighter in the model map. The disc vertical structure is likely playing a role here with the north-west side being the far side of the disk, thus associated to a larger projected emission area and a less severe beam dilution. In agreement with that, the upper side of the disc is also found to be brighter in some of the channels (compare Appendix B).
The brightness temperature map of 12 CO shows that the disk is also brighter, and thus likely warmer, in the north-eastern side. While 13 CO does not show any strong east-west asymmetry, the disc is clearly brighter in the south-western side of C 18 O, which instead of a higher temperature may trace a small over-density due to the lower optical depth of the line.  Figure 2 presents the normalized azimuthal variations of the peak intensity (before T B conversion) for the three CO isotopologues. Each profile is shown for the ring of strongest intensity between 0.12 − 0.25 (∼ 20-40 au) with the uncertainty being the rms computed on the peak intensity map. To rotate and deproject the maps an inclination and position angle of 35 • and 235 • respectively (compare results of Sect. 4.2.2) were used.

Radial and azimuthal cuts
As discussed above we indeed find an opposite east-west asymmetry in the curves for 12 CO and C 18 O, while the intensity of 13 CO is relatively symmetric along the azimuth. The intensity of the left (east) peak compared to the right (west) peak is roughly roughly 10 % higher for the 12 CO and 25 % lower for the C 18 O peak intensity. This matches the brightness and temperature differences seen in the maps of Fig. 1  In addition, the profiles underline that the under-brightness appears stronger in the south-east of the disc. For comparison, the azimuthal profiles derived from the DALI model are shown in the right panel of Fig. C.1 for all three isotopologues. In contrast to the data no strong east-west asymmetry is seen in the peaks. The peak intensity further drops by 35-45 % at the dips with the south-east dip being about 2-6 % deeper compared to the northwest dip. The asymmetry found in the under-brightness is thus more pronounced in the data. Together with the under-brightness seen in NIR, where beam dilution cannot be invoked to explain the latter, this supports the assumption that an additional shadowing may occur in the south side of the disc. While the continuum shows two clumps (top left panel of Fig. 1) they are present at a different location than the brightness asymmetries as well as the symmetric pair of dips and are thus unlikely to account for the variations. Figure 3 displays the normalized radial intensity profiles for the three CO isotopologues as well as the 1.3 mm continuum. The curves are obtained from the azimuthally averaged intensity per annuli of size 0.02 (∼ 3.2 au) from the peak intensity maps (before T B conversion), again using an inclination and position angle of 35 • and 235 • respectively. The error bars are calculated as the standard deviation per annulus divided by the square root of number of independent beams in the annulus. Compared with the profiles shown by Ubeira Gabellini et al. (2019) for the integrated intensity, we notice a drop of the peak intensity in the radial profile of C 18 O between ∼ 65-85 au of about 3 % in addition to the intensity drop at the ∼ 20 au cavity in 13 CO and C 18 O. A corresponding slight dip is present in the 13 CO peak intensity profile. Being more optically thin, C 18 O is mostly tracing the column density. Therefore the observed feature may be indicative of a depleted region around 75 au, possibly carved by an unseen companion. Continuum absorption is unlikely to account for the dip since the peak of the continuum flux lies around 52 au rather than 75 au. Another explanation may be the enhancement of emission around ∼ 90 au at the edge of the continuum, rather than a dip, potentially caused by an enhanced desorption of CO ices by increased UV or a temperature inversion (e.g. Cleeves et al. 2016;Facchini et al. 2017

Temperature structure
Since the 12 CO emission is optically thick and in LTE at these low rotational transitions, the brightness temperature (top right panel of Fig. 1) can be used as a probe of the gas kinetic temperature. In this context the gas temperatures up to 75 K that we observe are as they would be expected in the upper disc layers (Bruderer et al. 2014). To uncover small perturbations in this temperature structure, similar to Teague et al. (2019b) we subtract an azimuthally averaged radial T B profile similar to the one shown in the bottom panel of Fig. 3. This leaves significant spiral structure in the resulting residuals as shown in Fig. 4. Two clear spirals are observed, a smaller one (Sp T1 ) spanning an azimuth of ∼ 100 • between ∼ 10−180 au and a larger spiral (Sp T2 ) covering more than half an azimuth at a similar radial extent. Both spirals have the same orientation. The small spiral seen in the NIR by Uyama et al. (2020) at radii of ∼ 30-60 au co-locates with the anchoring point of the large spiral observed here (compare Fig. 7 in Sect. 4.3).

Velocity structure
The bettermoments package can also be used to generate the observed rotation velocity of the gas (similar to a moment 1 map). The code fits a quadratic model to the brightest pixel as well as the two neighbouring pixels to find the centroid of the line in pixel coordinates. Compared to other methods this approach is more robust to noise or errors in the line shape, allowing a precision that is greater than the velocity resolution. The resulting velocity structure of the disc, masking regions below 4 σ for 12 CO and 3.5 σ for 13 CO and C 18 O, is shown in Fig. 5. Besides the velocity field of all three isotopologues the corresponding error map of 12 CO is included in the top right panel. These statistical uncertainties are calculated by linearizing and propagating the uncertainty from the fluxes to the centroid estimate. For most regions the achieved precision is well below the channel width of 500 m s −1 . In the central regions the uncertainties increase due to beam smearing.
The isovelocities of the outer disc in the 12 CO velocity map match those of a Keplerian flat disc model (e.g. Rosenfeld et al. 2013), whereas significant distortions can be noticed in the in- ner disc (up to ∼ 0.5 ) with the kinematics in the centre being slightly twisted and the blue-and red-shifted parts bending in opposite directions. In the top left panel of Fig. 5 the maximum and minimum velocities along the red-and blue-shifted minor and major axis are overlaid, emphasizing the non-Keplerian term present in the inner disc. For the case of a razor-thin Keplerian disc a dipole morphology would be expected that is symmetric about the semi-major axis. The isovelocities further hint towards a rather non-or slightly elevated or flared emission surface since the lobes of the rotation pattern are overall not distinctively bent away (in one direction) from the disc major axis, although this may be resulting from the perturbing spiral structure. We still attempted to fit for the emission surface, however none of the fits converged. In the following we will thus focus on a razor-thin disc Keplerian model. The rotation pattern of 13 CO shows a similar twisting and bending as the 12 CO emission while no substructure can be discerned in the less bright C 18 O.

Analysis of the gas rotation velocity
To analyze the gas kinematics of the disc around CQ Tau and characterize the apparent deviations from Keplerian velocity we fit a Keplerian profile with (r, φ) being the deprojected cylindical coordinates, i the inclination of the disc and v LSR the systemic velocity to the rotation map of 12 CO shown in Fig. 5 using the eddy code (Teague 2019). The associated uncertainties are included in the fit. In order to deproject the sky-plane coordinates (x, y) into the midplane cylindrical coordinates (r, φ), the disc centre (x 0 , y 0 ), i and the disc position angle PA are used. The latter is measured between the north and the red-shifted semi-major axis in an easterly direction. As a first step, the starting positions of the free fit parameters are optimized with scipy.optimize with their posterior distributions estimated using the MCMC sampler. In this context we used 200 walkers, 5000 steps to burn in and 5000 additional steps to sample the posterior distribution function. For all of our models we assumed flat priors that were allowed to vary over a wide range. The uncertainties of the posterior distributions represent the 16th to 84th percentiles about the median value.
In addition to the razor-thin disc model, a parameterization for the emission surface as well as a warped structure can be included in the model. The results of our modelling are reported in the following and summarized in Table D.1 of Appendix D.

Razor-thin disc model
For the razor-thin disc models we fixed the object's distance and inclination to 162 pc and 35 In the first two runs we attempted to fit the entire disc, choosing an outer radius of 1.0 (162 au) to exclude possible noise at the disc's edge. For the second run we further set an inner boundary of 0.25 (40 au), which corresponds to the inner edge of the Gaussian ring of the dust continuum, as obtained by Ubeira Gabellini et al. (2019). Both set-ups result in very similar fit parameters, yet returning a slightly larger stellar mass when the disc centre is excluded. In addition, we tried to fit specific regions of the disc, including only the inner disc (run 3,4), outer disc (run 5,6) as well as annuli of size 0.2 (run 7-10). Overall we find that v LSR is slightly increasing towards the outer disc, while PA is relatively constant. The largest scatter is found for the stellar mass M * , that is ranging from 1.47 to 1.65 M , driven by the model trying to account for the non-Keplerian structure in the rotation map.
All thin disc models rapidly converged with a Gaussian posterior distribution function (PDF), resulting in relatively similar residuals when the model is subtracted from the velocity data. We tried both convolving the models with the beam of the observation and not using the convolution, with both approaches returning comparable results. In this context we note that convolving channel maps prior to collapsing them into the rotation map would be the better approach than convolving the model map as it is done in eddy. Generating channels maps, as opposed to a simple rotation map, however requires far more model assumptions which is why we choose not to do it. The effects are negligible outside the disc centre (i.e. outside ∼2 beam FWHM) and thus do not significantly affect our results.
The posterior distributions presented in Table D.1 show very small and likely underestimated uncertainties, especially in context of the scatter that is found in the stellar mass. One possibility is that the uncertainties in the velocity centroid are underestimated. We thus performed an additional run of model 2 with the velocity errors increased by a factor of 10. This returns very similar fit parameters with the uncertainties also being a factor 10 larger, however still too small to account for the observed scatter. Therefore the small uncertainties cannot only be explained by underestimated velocity errors but result from systematic uncertainties in our model.
In Fig. 6 the results from run 2 are presented, including the best-fit model (left panel) and the corresponding residuals (right panel) after subtraction from the data. The model corresponds to a position angle PA = 235 • , systemic velocity v LSR = 6173 m s −1 and a stellar mass of M * = 1.57 M . While the residuals are less than about 5 % of the velocity data outside of ∼ 0.2 , they (partly) grow to more than 50 % in the very inner disc (< 0.1 ) which suffers strongly of beam smearing effects (Teague et al.   2018c, 2016). We thus masked out the residuals inside of 1.5 the FWHM of the beam. They clearly reveal the non-Keplerian structure of the rotation velocity, showing significant spiral features that cover more than one azimuth at radii of 40-180 au, with the same orientation and a similar location as the spirals observed in the gas temperature. Since twisted kinematics are sometimes linked to the presence of a misaligned inner disc, we performed several runs, adding the parameterization of a potential warp in the (flat disc) model. Yet none of these models converged. Besides being limited by the spatial resolution of the data (0.121 × 0.098 ) the kinematics are strongly dominated by the large spiral structure in the outer disc. Thus a small feature such as a warp in the very inner disc regions can not be fit by our simple model. For the same reason we were not able to obtain any constraints on the emission surface.

Analysis of the spiral structure
Three significant spirals are observed in the residuals presented in Sect. 4.1 and Sect. 4.2, two in the brightness temperature and one in the rotation velocity of 12 CO. All spirals show the same orientation, suggesting a counter-clockwise rotation of the disc if the spirals are trailing, and they cover a large azimuth and radial extent. A counter-clockwise rotation implies that the south-east side of the disc is closer to the observer, which should be visible in the rotation map as a bending of the high velocity components towards the north. While this is true for the blue-shifted side, the spiral may be driving such a large velocity perturbation, that the red-shifted side is bending the other way, then resembling the presence of a warp. The geometry of the disc agrees with the SE dip in the peak brightness temperature map being dimmer due to beam dilution (see Sect. 3.1).
We reproduce the spiral morphology with different functional forms, in particular an Archimedean (linear) spiral as well as a spiral r = a + be kφ (4) where r represents the radius and φ the polar angle of the spiral. Since Eq. 4 is similar to the equation of a logarithmic spiral, we will refer to it as such in the following. The resulting parameters are presented in Table 2 Fig. 7 for the brightness temperature and in Fig. 8 for the rotation velocity. In both cases we plot the deprojected and rotated maps, using the inclination and position angle found in Sect. 4.2.2 (top plots). Again the inner disc regions are masked out inside of 1.5 the FWHM of the beam for the velocity residuals. Additionally the polar-deprojected maps are shown (bottom plots). The spiral fits are overlaid as Sp T1 and Sp T2 in the temperature and as Sp v in the velocity with white and black dashed lines. Both functions provide a good and very similar parameterization for the small temperature spiral Sp T1 , with the overall lin-ear nature becoming clear when looking at the polar-deprojected map. For the second, large temperature spiral Sp T2 , the linear spiral is not able to account for the curvature obvious in the polardeprojected plot and is thus missing the inner part of the disc. The logarithmic parameterization on the other hand is able to better represent the anchoring point of the spiral. The velocity spiral Sp v is again well represented by both functions. No large differences can be seen between the two parameterizations, with the overall nature being very close to linear. As a comparison the Sp T2 spiral is also included in the velocity residuals. While the two spirals are co-located in the outer regions, they deviate from each other towards the inner disc for both spiral functions.
The spiral found in the NIR is overlaid on the brightness temperature residuals in Fig. 7. Here the solid line represents the location of the NIR spiral, while the dashed line shows the extrapolation of the fit performed by Uyama et al. (2020) for this spiral with a function r = a + bφ n . The NIR spiral matches the anchoring point of the temperature spirals and follows the course of Sp T2 at the inner edge, suggesting them to be connected.
In addition to the 12 CO data we search for features in the residuals of 13 CO and C 18 O, presented in Fig. 9 for T B (top panels) and the velocity (bottom panels). The residuals are calculated again by subtracting an azimuthally symmetric model, using the best-fit model from run 2 for the velocity. The loga- rithmic spirals Sp T2 and Sp v observed in 12 CO are overlaid for comparison. Even though no clear spiral structure can be found in the brightness temperature of either 13 CO or C 18 O, the launching point of the spiral is clearly visible in the 13 CO residuals and some indication of a spiral matching the logarithmic parameterization of Sp T2 is present. Similarly, the velocity residuals of 13 CO also show indication for a spiral similar to that observed in 12 CO, yet slightly more tightly wound, while no substructure can be distinguished in C 18 O. Using we calculate the pitch angle β for all three spirals. The resulting angles are shown in Fig. 10 for the linear and logarithmic fit.
Here the solid lines are shown at the location where the spirals are present, while the dashed lines represent the extrapolation towards smaller radii. We further included the constant pitch angle found by Uyama et al. (2020), who fitted a logarithmic spiral of r = be kφ to a spiral feature observed in the NIR. As mentioned above this spiral is located at the anchoring point of spiral Sp T2 , yet spanning a smaller radial extent of ∼ 30-60 au and azimuth of roughly 50 • . The pitch angle of the NIR spiral, which is given as ∼ 34 • , seems consistent with the pitch angle found from the linear parameterization for Spiral Sp T2 , however in these radial regions the linear fit failed to reproduce the spiral morphology (shaded region in Fig. 10). Here the logarithmic spiral provided a better approximation with the resulting pitch angles lying 20-40 • above that of the NIR spiral. The pitch angle is decreasing faster with radius for the logarithmic spirals, with the angles being larger until ∼ 126 au (Sp T1 ), 147 au (Sp T1 ) and 106 au (Sp v ) compared to the linear spiral, yet overall comparable for Sp T1 and Sp v (keeping in mind an uncertainty due to the by-eye parameterization). This is expected, since both functions provide a similar parameterization for these spirals. The largest difference is seen for Sp T2 , where the Archimedean spiral could not sufficiently reproduce the inner parts of the spiral.

Discussion
Both the gas temperature and the rotation velocity of 12 CO show significant spiral structure over the bulk of the disc when an azimuthally symmetric model is subtracted from the observations (Figs. 7,8). Together with a similar feature found in the NIR (Uyama et al. 2020) the extent of the structures over a large azimuth (> 180 • temperature, > 360 • velocity) and radius (10-180 au temperature, 40-180 au velocity) with a SNR of 2-6 suggests them to be real features. Both the large temperature spiral SP T2 and the NIR spiral appear to follow a similar course, strongly suggesting a link between the two. The location of the spirals matches those of the prominent distortions and bendings occurring in the velocity field and are possibly the cause for the latter. Higher spatial and spectral resolution will be necessary to resolve the very inner disc regions and confirm the presence of the spirals at a higher SNR.
Even though no significant spiral structure is observed in the 13 CO and C 18 O lines, indications for such features are present in both the brightness temperature and rotation velocity residuals of 13 CO and may be made accessible with deeper high resolution observations.
All three spirals observed in 12 CO are well described by a modified logarithmic spiral as well as two (except SP T2 ) by an Archimedian (linear) spiral with radially decreasing pitch angles (Fig. 10). The angles found for the logarithmic parameterization mostly lie above those of a linear one. Overall all three spirals are loosely wound and consequently show relatively large pitch angles. These characteristics may be explained through Lindblad-resonance driven spiral wakes of a massive embedded companion (e.g. Ogilvie & Lubow 2002;Rafikov 2002;Bae & Zhu 2018a,b).
Since the optically thick 12 CO is tracing higher disc layers, significantly smaller pitch angles, and thus more tightly wound spirals, are expected at the midplane if the disc is passively heated with a positive vertical temperature gradient (compare Juhász & Rosotti 2018). In this context, to distinguish between different spiral launching scenarios it is crucial to use observations spanning the full vertical extent. If the spirals were for example caused by gravitational instability, similar pitch angles are expected for the midplane and the surface layers since the midplane will be heated by shocks in that case. Given the small number of spirals, gravitational instability however seems a rather unlikely cause for the spiral structure we observe (e.g. Cossins et al. 2009;Hall et al. 2020). On the other hand, several spiral arms could possibly appear as only two or three arms due to resolution effects (Dipierro et al. 2014) and the disc around CQ Tau happens to be relatively massive, thus gravitational instability cannot be ruled out at this point.
So far no clear spiral structures have been found in the dust continuum (or optically thin lines) of CQ Tau to further distinguish possible launching scenarios but could be made accessible with higher spectral and spatial resolution data. Spiral arms have been observed with ALMA in the continuum of several discs, including for example Elias 2-27 (Pérez et al. 2016), IM Lup and WaOph 6 , G17.64+0.16 (Maud et al. 2019), MWC 758 (Boehler et al. 2018;Dong et al. 2018) and HD100453 (Rosotti et al. 2020a). For the latter, counterparts to the observed NIR spirals (Wagner et al. 2015;Benisty et al. 2017) were not only found in the dust continuum but also the CO emission, enabling the authors to study the thermal structure of the disc and link the spirals to a known binary companion. The velocity residuals of 13 CO possibly suggest a spiral slightly more tightly wound than the corresponding spiral in 12 CO. A more tightly wound spiral would show smaller pitch angles, which would support Juhász & Rosotti (2018).
It is difficult to determine whether the spirals observed in the temperature and the velocity are tracing the same underlying perturbation. Even though the spirals Sp T2 and Sp v do not fully overlap, they appear to align in the outer parts of the disc between ∼ 130-180 au and 0 to -180 • , hinting towards the same formation mechanism. On the other hand, the calculated pitch angles for the according spirals differ by several degrees. We note however, that the pitch angles are only a rough estimate, since no actual fit was performed and that the actual pitch angles may lie much closer for Sp T2 and Sp v .
Similarly Teague et al. (2019b) observe temperature and velocity spirals in TW Hydra, that appear to align but do not fully overlap, thus there may be a physical mechanism behind these differences. The authors suggest that layers with different ther-mal properties are traced, arguing that close to the discs surface spirals in the velocity should be more pronounced due to efficient cooling, while the heat produced by spirals would be more efficiently trapped closer to the midplane. In the case of a companion, the spiral density waves created by either a planet or binary companion will lead to an increase of surface density and thus in a higher CO opacity. This will move the τ = 1 layer to a higher altitude, where the temperature in generally higher, resulting in the observed spiral substructure in the gas temperature (Phuong et al. 2020a,b).
Even though it is impossible to fully disentangle all three velocity components (v r , v z , v φ ), the same orientation in the residuals as well as the full azimuthal coverage hint towards a vertical perturbation (compare Appendix B in Teague et al. 2019b). This is consistent with the (potentially) companion launched spirals in HD 100453 (Rosotti et al. 2020a). As shown by Pinte et al. (2019) an embedded planet will cause perturbations in all three velocity directions, with their strength decreasing with height above the midplane for radial and rotational motions, whereas they increase with height for vertical motions. Gas flowing towards the midplane and falling into the observed cavity could be another explanation for the vertical motions. Since the vertically moving material will receive more stellar light it may appear brighter regardless of the underlying surface density.
In case the spirals are indeed launched by an embedded planet outside of its orbit, they are expected to converge towards the planet location (Juhász et al. 2015;Bae & Zhu 2018a,b). Since this results in a rapid increase of the pitch angle towards the planet, a possible companion is expected to be located inside of ∼ 25 au in our case. This is consistent with the findings of Ubeira Gabellini et al. (2019), who propose an unseen planet of 6 − 9 M Jup to explain the deep dust and gas cavity. Given such a location and mass, a companion is not expected to be noticeably affected by extinction in any band (Sanchis et al. 2020). We note, that dynamically launched spirals tend to open up only close to the planet consequently becoming more tightly wound at larger distances, resulting in small pitch angles. It is thus puzzling, that the spirals we observe are still very open far from the possible companion.
In the channel maps presented in Figs. B.1, B.2 and B.3 no clear kinks similar to Pinte et al. (2018Pinte et al. ( , 2019 are detected and the two sides of the disc cannot be resolved due to spectral and spatial resolution limitations. Since several studies (Uyama et al. 2020;Ubeira Gabellini et al. 2019) including ours on the other hand indicate the presence of a massive companion of CQ Tau at < 25 au, imprints on the iso-velocities are expected and may be made accessible with higher resolution or sensitivity.

Summary
In this work we presented high angular resolution ALMA observations of 12 CO, 13 CO and C 18 O J = 2 − 1 data of the disc around CQ Tau and used the 12 CO data to analyze the gas temperature and kinematics of the disc. The main results of this analysis are summarised in the following.
The morphology of the significant spiral structure observed in the brightness temperature and rotation velocity of 12 CO together with the number of spirals and large pitch angles supports a dynamical launching scenario, for example an embedded planet or binary rather than gravitational instabilities. Such a companion is expected to be relatively massive and to be located inside of ∼ 25 au, which is in agreement with Ubeira Gabellini et al. (2019). Further multi-line observations at a higher velocity resolution as well as 3D modelling are required to further distinguish the different mechanisms.
In addition to the gas cavity an intensity drop can be seen in the north-west and south-east side of the disc in all three isotopologues. Postprocessing the DALI model presented by Ubeira Gabellini et al. (2019) revealed the under-brightness to be caused by beam dilution, rather than by a temperature or line width effect. Since the dip in the south-east appears to be more pronounced, and is co-located with the under-brightness in scattered light, some additional shadowing may still occur, potentially caused by the spiral itself or misaligned regions in the disc. Figure 11 presents an illustration of the possible morphology of the disc around CQ Tau, taking into account the results from Ubeira Gabellini et al. (2019), Uyama et al. (2020) and this work: a dust and gas cavity are present in the disc at ∼ 50 and ∼ 20 au respectively which could be explained by a massive companion inside of 20 au. Spiral structure is found in both T B , v rot and the NIR, further supporting the hypothesis of an unseen companion. The inner disc regions, including their position angle and inclination, remain unresolved.
Altogether it appears that the disc around CQ Tau is far from a Keplerian disc. Thus a more detailed non-Keplerian model is required to describe the gas rotation and could be addressed in a future work. Especially the possibility of a massive companion, either a binary star or very massive planet, needs to be further explored. To construct such a model higher angular and spectral resolution ALMA data are essential. Furthermore, NIR interferometry observations at milliarcsec resolution with VLTI-Gravity (Gravity Collaboration et al. 2017) could help to constrain the inclination and position angle of the innermost disc, providing information on the presence of misaligned regions. Combining dust and gas observations of different molecular lines with high SNR is needed to constrain the vertical temperature profile and to further analyse the clearly observed spiral structure of the gas disc.   [arcsec] [arcsec]