© L. Wölfer et al. 2021

1 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 substructures 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 that could also account for the observations, such as the magnetorotational instability (MRI; e.g. Flock et al. 2015, 2017; Riols & Lesur 2019), zonal flows (e.g. Uribe et al. 2015), the compositional baroclinic instability (Klahr & Bodenheimer 2004) and gravitational instability (Kratter & Lodato 2016).

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). As 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 therefore of paramount importance to directly access and characterise 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 (1)

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 characterise 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. (2018a, 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. (2020a) to measure the gas-dust coupling as well as the width of gas pressure bumps. Pinte et al. (2018, 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 reported in a few discs of the ALMA DSHARP large program (Andrews et al. 2018; Pinte et al. 2020), but more data are needed to confirm the robustness of such claims. Similarly, a possible signature for an embedded planet in the HD 100546 disc was 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. 2015, 2018).

One type of disc, the so-called transition disc, is of particular interest, because examples 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 us 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, characterisation of the present dust and gas cavity and possible formation mechanisms, we analyse 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 the method used to extract and model the gas kinematics. Our results are discussed in Sect. 5 and summarised in Sect. 6.

2 Observations

2.1 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 2018) (RA = 05^h35^m58.47^s, Dec = +24°44′54.09″; J2000). The intermediate-mass star (1.67 M_⊙; Garcia 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⁻⁸ –10⁻⁷ M_⊙ yr⁻¹ (Donehew & Brittain 2011; Mendigutía et al. 2012).

The disc around CQ Tau represents one of the first discs whose millimetre (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 itsdust 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, 2003; Chapillon et al. 2008). The average dust opacity coefficient was constrained by Banzatti et al. (2011) using VLA (1.3–3.6 cm), PdBI (2.7–1.3 mm), and SMA (0.87 mm) observations, probing significant grain growth in the disc with up to centimetre (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 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 ¹³ CO and C¹⁸ 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 Fig. 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.

2.2 Data reduction

We present 1.3 mm ALMA observations of the CQ Tau system in band 6, combining datasets from cycles 2, 4, and 5 (2013.1.00498.S, PI: L. Pérez; 2016.A.00026.S, 2017.1.01404.S., PI: L. Testi), 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 configurations, but they share a similar spectral setup. The longest baseline from the first project extended to 1091 m, while for the latest two it was increased to 3700 m and 8500 m respectively, thus enhancing the spatial resolution. In all three observationsthe ALMA correlator was configured to observe the 1.3 mm dust continuum emission, as well as the molecular lines ¹² CO J = 2−1, ¹³ CO J = 2−1, and C¹⁸ O J = 2−1.

After applying ALMA standard pipeline calibration, we followed a similar processing as for the DSHARP data calibration (Andrews et al. 2018), using CASA 5.4.1. We started by flagging the channels located at ± 25 km s⁻¹ from each spectral line, and averaged the remaining channels to 125 MHz width channels, which were combined with the data from the continuum spectral windows. As a next step, we aligned the dust continuum emission and checkedby comparison the flux calibration of each individual execution, to ensure they all have the same flux.

To enhance the signal-to-noise ratio (S/N), self-calibration was performed in two stages. First, we self-calibrated 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 was then combined with the extended baseline datasets from the observationsobtained in Cycles 4 and 5 (datasets #2 and 3), and four phase calibrations with solution intervals of 900, 360, 150, and 90 s, as well as one amplitude calibration with a solution interval of 360 s.

All the dust continuum emission calibration steps, including the centroid shifting and self-calibration tables, were then applied to the molecular line emission channels. The continuum emission was subtracted using the uvcontsub task, and image cubes were 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 was calculated with the package keplerian_mask¹, 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⁻¹. These values were chosen after some initial fits for the gas kinematics, explained later in Sect. 4.2. We further chose an inner and outer radius of 0 and 2′′ respectively and convolved with the beam rescaled by 1.5 times its size. Some important characteristics of the data are given in Table 1.

3 Observational results

3.1 Integrated intensity and 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 ¹² 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 we 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 S/N. This results in a peak S/N of 40 for the ¹² CO, 20 for the ¹³CO, and 14 for the C¹⁸ O velocity-integrated intensity as well as 29for the ¹²CO, 19 for the ¹³CO, and 18 for the C¹⁸O peak intensity.

While the optically thick ¹²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 ¹³CO and C¹⁸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 in the peak intensity of all isotopologues and integrated intensity of C¹⁸ O by roughly 35–60% (with respect to the peak) in the northwest and southeast parts of the disc along the minor axis (symmetric pair of dips), which becomes especially prominent in the C¹⁸ O data and co-locates with a possible under-brightness in near-infrared (NIR) scattered light (see Fig. 4 of Uyama et al. 2020). Similar to the NIR under-brightness, the drop in peak intensity appears to be more pronounced (relative to the peak ~7–12% deeper dip) on the southeast side of the disc (compare also Sect. 3.2).

Such symmetric features are often linked to the presence of a misaligned inner disc, casting a shadow over the outer disc (e.g. 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 should be exercised, because 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 used the DALI (Dust And LInes; Bruderer et al. 2012; Bruderer 2013) model presented by Ubeira Gabellini et al. (2019) and convolved the spectral image cubes with the beam of our observation. The resulting peak intensity map is shown in the left panel of Fig. C.1 for ¹² 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 northwest dip) appears slightly brighter in the model map. The disc vertical structure is likely playing a role here with the northwest side being the far side of the disc, 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 ¹²CO shows that the disc is also brighter, and thus likely warmer, on the northeastern side. While ¹³ CO does not show any strong east–west asymmetry, the disc is clearly brighter on the southwestern side of C¹⁸ O, which instead of a higher temperature may trace a small over-density due to the lower optical depth of the line.

Fig. 1 ALMA observations of the velocity-integrated intensity (left panels) and peak brightness temperature (right panels) of the 12 CO (top panels),13CO (middle panels), and C18O (bottom panels) J = 2−1 transition. The conversion from flux to brightness temperature was performed with the Planck law. The contours of the continuum are overplotted on top of the 12 CO integrated intensity at 20, 100, and 130 σ (1 σ = 11 μJy beam−1). The synthesized beam is shown in the bottom left corner of each panel.

3.2 Radial and azimuthal cuts

Figure 2 presents the normalised 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 a position angle of 235° and inclination of 35° (compare results of Sect. 4.2.2) were used.

As discussed above we indeed find an opposite east–west asymmetry in the curves for ¹² CO and C¹⁸ O, while the intensity of ¹³CO is relatively symmetric about the y-(minor) axis. The intensity of the left (east) peak compared to the right (west) peak is roughly 10% higher for the ¹²CO and 25% lower for theC¹⁸O peak intensity. This matches the brightness and temperature differences seen in the maps of Fig. 1, where the east side of the disc is brighter in ¹² CO but fainter in C¹⁸O.

In addition, the profiles underline that the under-brightness appears stronger on the southeast side 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 roughly 35–45% at the dips with the southeast dip being (relative to the peak) 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 on 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 from the brightness asymmetries and from the symmetric pair of dips and are therefore unlikely to account for the variations.

Figure 3 displays the normalised 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 the 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¹⁸ O between ~65 and 85 au of about 3% (relative to the peak value) in addition to the intensity drop at the ~ 20 au cavity in ¹³ CO and C¹⁸ O. A corresponding slight dip is present in the ¹³CO peak intensity profile. Being more optically thin, C¹⁸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 because 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).

Fig. 2 Azimuthal variation of the peak intensity for an annulus of 20–40 au (0.′′ 12-0.′′25), normalised to the peak value and shown for the three CO isotoplogues. The uncertainties are given as 1 σ.

Fig. 3 Azimuthally averaged radial intensity profiles of the continuum (red), 12CO (blue), 13CO (orange), and C18O (green) data. The profiles are derived from the corresponding peak intensity maps and normalised to the peak value.

Fig. 4 Residuals of the 12CO brightness temperature after subtraction of an azimuthally averaged radial profile. Two spirals SpT1 and SpT2 are spanningan azimuth of ~100° and >180° respectively between ~10 and 180 au.

4 Analysis

4.1 Temperature structure

As the ¹²CO emission is optically thick and in local thermodynamic equilibrium (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 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 and 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).

Fig. 5 Rotation map of the three CO isotopologues and the corresponding uncertainties of 12 CO (top right panel) calculated with bettermoments. Regions below 4 σ (12 CO) and 3.5 σ (13 CO, C18 O) are masked out. The maximum and minimum velocities along the red- and blueshifted major axes are overlaid with grey lines on the 12CO rotation pattern.

4.2 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 ¹² CO and 3.5 σ for ¹³ CO and C¹⁸ O, is shown in Fig. 5. Besides the velocity field of all three isotopologues, the corresponding error map of ¹² CO is included in the top right panel. These statistical uncertainties are calculated by linearising 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⁻¹. In the central regions, the uncertainties increase due to beam smearing.

The isovelocities of the outer disc in the ¹²CO velocity map match those of a Keplerian flat disc model (e.g. Rosenfeld et al. 2013), whereas significant distortions can be seen in the inner disc (up to ~0.′′5) with the kinematics in the centre being slightly twisted and the blue- and redshifted parts bending in opposite directions. In the top left panel of Fig. 5, the maximum and minimum velocities along the red- and blueshifted major axes are overlaid, emphasising 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 because 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 therefore focus on a razor-thin disc Keplerian model.

The rotation pattern of ¹³CO shows a similar twisting and bending to the ¹²CO emission while no substructure can be discerned in the less bright C¹⁸ O.

4.2.1 Analysis of the gas rotation velocity

To analyse the gas kinematics of the disc around CQ Tau and characterise the apparent deviations from Keplerian velocity, we fit a Keplerian profile (2)

with (r, ϕ) being the deprojected cylindical coordinates, i the inclination of the disc and v_LSR the systemic velocity to the rotation map of ¹²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₀, y₀), i, and the disc position angle PA are used. The latter is measured between the north and the redshifted semi-major axis in an easterly direction. As a first step, the starting positions of the free fit parameters are optimised 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 parameterisation 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 summarised in Table D.1.

Fig. 6 Best-fit Keplerian rotation model (left panel) and the residuals calculated by subtracting the model from the observed rotation velocity (right panel), shown for 12 CO. The residuals inside a radius of 1.5 times the FWHM of the beam are masked out. In the residuals a spiral Spv is spanning more than one azimuth between ~40 and 180 au.

4.2.2 Razor-thin disc model

For the razor-thin disc models we fixed the object’s distance and inclination to 162 pc and 35° (Ubeira Gabellini et al. 2019) respectively, and fitted for the disc centre (x₀, y₀ ∈{−0.′′5, 0.′′5}), systemic velocity (v_LSR ∈{−5 km s⁻¹, 20 km s⁻¹}), stellar mass (M_*∈{0.1 M_⊙, 5 M_⊙}), and disc position angle (PA ∈{−360°, 360°}).

In the first two runs (runs 1, 2 in Table D.1), we attempted to fit the entire disc, choosing an outer radius of 1′′ (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 setups result in very similar fit parameters, yet returning a slightly smaller stellar masswhen 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), or annuli of size 0.′′ 2 (runs 7–10). Overall we find that v_LSR slightly increases towards the outer disc, while PA is relatively constant. The largest scatter is found for the stellar mass M_*, which ranges 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 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 a better approach than convolving the model map as it is done in eddy. However, generating channel maps, as opposed to a simple rotation map, 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 approximately two times the 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 the context of the scatter that is found in the stellar mass. One possibility is that the uncertainties in the velocity centroid are underestimated. We therefore performed an additional run of model 2 with the velocity errors increased by a factor of ten. This returns very similar fit parameters with the uncertainties also being a factor of ten larger, albeit 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.

Figure 6 shows the results from run 2, 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⁻¹, 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 from beam smearing effects (Teague et al. 2016, 2018c). We therefore masked out the residuals inside a radius of 1.5 times the FWHM of the beam. They clearly reveal the non-Keplerian structure of the rotation velocity, showing significant spiral features (Sp_v) 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.

As twisted kinematics are sometimes linked to the presence of a misaligned inner disc, we performed several runs, adding the parameterisation of a potential warp in the (flat disc) model. Nevertheless, 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 cannot be fitted by our simple model. For the same reason, we were not able to obtain any constraints on the emission surface.

4.3 Analysis of the spiral structure

Three significant spirals are observed in the residuals presented in Sects. 4.1 and 4.2: two in the brightness temperature and one in the rotation velocity of ¹² 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 southeast 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 blueshifted side, the spiral may be driving such a large velocity perturbation that the redshifted side is bending the other way, therefore resembling a warp. The geometry of the disc agrees with the southeast 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 (3)

as well as a spiral, (4)

where r represents the radius and ϕ the polar angle of the spiral. As Eq. (4) is similar to the equation of a logarithmic spiral, we refer to it as such in the following. The resulting parameters are presented in Table 2 and the corresponding spirals are shown in 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 a radius of 1.5 times 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 residuals with white and black dashed lines.

Both functions provide a good and very similar parameterisation for the small temperature spiral Sp_T1, with the overall linear 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 that is obvious in the polar-deprojected plot and is thus missing the inner part of the spiral. The logarithmic parameterisation 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 parameterisations, with the overall nature being very close to linear. As a comparison, the Sp_T2 spiral is alsoincluded 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ϕⁿ. The NIR spiral matches the anchoring point of the temperature spirals and follows the course of Sp_T2 at the inner edge, suggesting they are connected.

In addition to the ¹²CO data, we search for features in the residuals of ¹³CO and C¹⁸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 logarithmic spirals Sp_T2 and Sp_v observed in ¹²CO are overlaid for comparison. Even though no clear spiral structure can be found in the brightness temperature of either ¹³ CO or C¹⁸ O, the launching point of the spiral is clearly visible in the ¹³CO residuals and some indication of a spiral matching the logarithmic parameterisation of Sp_T2 is present. Similarly, the velocity residuals of ¹³CO also show signs of a spiral similar to that observed in ¹²CO, yet slightly more tightly wound, while no substructure can be distinguished in C¹⁸ O.

Using (5)

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 the 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 parameterisation 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 decreases 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 parameterisation). This is expected, because both functions provide a similar parameterisation for these spirals. The largest difference is seen for Sp_T2, where the Archimedean spiral was not able to sufficiently reproduce the inner parts of the spiral.

Fig. 7 Deprojected and rotated (top panels) as well as polar-deprojected (bottom panels) TB residuals of 12CO with overlaid Archimedean and logarithmic spirals. The blue dashed lines show the fit of the spiral observed in the NIR by Uyama et al. (2020) extrapolated to large azimuthal angles. The solid blue line highlights the spiral region used by Uyama et al. (2020) to obtain the fit.

5 Discussion

Both the gas temperature and the rotation velocity of ¹² 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 S/N of 2–6 suggests they are 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 of 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 S/N.

Although no significant spiral structure is observed in the ¹³ CO and C¹⁸ O lines, indications for such features are present in both the brightness temperature and rotation velocity residuals of ¹³ CO and may be made accessible with deeper high-resolution observations.

All three spirals observed in ¹²CO are well described by a modified logarithmic spiral, and two (i.e. not SP_T2) are well described by an Archimedian (linear) spiral with radially decreasing pitch angles (Fig. 10). The angles found for the logarithmic parameterisation 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).

As the optically thick ¹²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 would be expected for the midplane and the surface layers because the midplane would 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). It is worth considering, nevertheless, that 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, therefore 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 that could be used to further distinguish possible launching scenarios, but these 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 (Huang et al. 2018), G17.64+0.16 (Maud et al. 2019), MWC 758 (Boehler et al. 2018; Dong et al. 2018), and HD100453 (Rosotti et al. 2020b). 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 ¹³ CO possibly suggest a spiral slightly more tightly wound than the corresponding spiral in ¹² CO. A more tightly wound spiral would show smaller pitch angles, which would support the findings of Juhász & Rosotti (2018).

It is difficult to determine whether the spirals observed in the temperature and the velocity trace the same underlying perturbation. Even though the spirals Sp_T2 and Sp_v do not fullyoverlap, they appear to align in the outer parts of the disc between ~130 and 180 au and from 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. However, we note that the pitch angles are only a rough estimate because 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, andtherefore there may be a physical mechanism behind these differences. The authors suggest that layers with different thermal properties are traced, arguing that close to the disc’s 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 in surface density and thus in a higher CO opacity. This will move the τ = 1 layer to a higher altitude, where the temperature is 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 coupled with 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. 2020b). 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. As 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). As this results in a rapid increase in 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, and consequently become more tightly wound at larger distances, resulting in small pitch angles. It is therefore 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. (2018, 2019) are detected and the two sides of the disc cannot be resolved due to spectral and spatial resolution limitations. As several studies (Ubeira Gabellini et al. 2019; Uyama et al. 2020), 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.

Fig. 8 Deprojected and rotated (top panels) as well as polar-deprojected (bottom panels) velocity residuals of 12 CO with overlaid Archimedean and logarithmic spirals. The spiral SpT2 observed in TB is overlaid as a grey dashed line. The residuals inside a radius of 1.5 times the FWHM of the beam are masked out.

Fig. 9 Deprojected and rotated TB (top) and velocity (bottom) residuals of 13CO (left) and C18O (right). The logarithmic spirals SpT2 and Spv observed in 12CO are overlaid.The residuals inside a radius of 1.5 times the FWHM of the beam are masked out in the velocity residuals.

Fig. 10 Pitch angle of the three spirals observed in the TB and rotation velocity residuals of 12CO, shown for the linear and logarithmic spiral. The constant pitch angle of the NIR spiral found by Uyama et al. (2020) is included as a reference. The shaded region represents the disc regions where the linear fit fails to reproduce the morphology of spiral SpT2.

6 Summary

In this work, we present high-angular-resolution ALMA observations of ¹² CO, ¹³ CO, and C¹⁸ O J = 2−1 data of the disc around CQ Tau and use the ¹²CO data to analyse the gas temperature and kinematics of the disc. The main results of this analysis are summarised as follows.

The morphology of the significant spiral structure observed in the brightness temperature and rotation velocity of ¹² 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 and 3D modelling are required to further distinguish the different mechanisms.

In addition to the gas cavity, an intensity drop can be seen on the northwest and southeast sides 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. As the dip in the southeast 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 T_B, v_rot, and 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. Therefore, a more detailed non-Keplerian model is required to describe the gas rotation and could be addressed in a future work. The possibility of a massive companion in particular, 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 2017) could help to constrainthe 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 S/N is needed to constrain the vertical temperature profile and to further analyse the clearly observed spiral structure of the gas disc.

Fig. 11 Illustration of the possible morphology of the disc around CQ Tau.

Acknowledgements

We thank the referee Simon Casassus for his constructive comments and suggestions that helped to improve our work. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00498.S, ADS/JAO.ALMA#2016.A.000026.S, and ADS/JAO.ALMA#2017.1.01404.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan) and KASI (Republic ofKorea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, auI/NRAO and NAOJ. We acknowledge the support from the DFG Research Unit "Planet Formation Witnesses and Probes: Transition Discs" (FOR 2634/1, ER 685/8-1 and ER 685/11-1). N.K. acknowledges support provided by the Alexander von Humboldt Foundation in the framework of the Sofja Kovalevskaja Award endowed by the Federal Ministry of Education and Research. GR acknowledges support from the Netherlands Organisation for Scientific Research (NWO, program number 016.Veni.192.233). This work was partly supported by the Italian Ministero dell’Istruzione, Università e Ricerca (MIUR) through the grant Progetti Premiali 2012 iALMA (CUP C52I13000140001), by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) ref. No. FOR 2634/1 TE 1024/1-1, and by the DFG cluster of excellence Origin and Structure of the Universe (www.universe-cluster.de), and by the EU Horizon 2020 research and innovation programme, Marie Sklodowska-Curie grant agreement NO 823 823 (Dustbusters RISE project).

Appendix A ALMA observing log

Appendix B Channel maps

Fig. B.1 12CO J = 2− 1 line imaged with a channel width of Δv = 0.5 km s−1 and a beam of 0.′′121 × 0.′′098.

Fig. B.2 13CO J = 2− 1 line imaged with a channel width of Δv = 0.7 km s−1 and a beam of 0.′′128 × 0.′′103.

Fig. B.3 C18O J = 2− 1 line imaged with a channel width of Δv = 1.0 km s−1 and a beam of 0.′′129 × 0.′′103.

Appendix C Beam dilution

Fig. C.1 Left: peak intensity map from postprocessing the DALI model presented by Ubeira Gabellini et al. (2019), shown for 12CO. Right: azimuthal variations of the peak intensity of 12CO, 13 CO, and C18 O for an annulus of 20–40 au (0.′′12-0.′′25), derived from the DALI model and normalised to the peak value.

Appendix D Posterior distributions of the modelling

References

All Tables

All Figures

Fig. 1 ALMA observations of the velocity-integrated intensity (left panels) and peak brightness temperature (right panels) of the 12 CO (top panels),13CO (middle panels), and C18O (bottom panels) J = 2−1 transition. The conversion from flux to brightness temperature was performed with the Planck law. The contours of the continuum are overplotted on top of the 12 CO integrated intensity at 20, 100, and 130 σ (1 σ = 11 μJy beam−1). The synthesized beam is shown in the bottom left corner of each panel.

In the text

	Fig. 2 Azimuthal variation of the peak intensity for an annulus of 20–40 au (0.′′ 12-0.′′25), normalised to the peak value and shown for the three CO isotoplogues. The uncertainties are given as 1 σ.
In the text

	Fig. 3 Azimuthally averaged radial intensity profiles of the continuum (red), 12CO (blue), 13CO (orange), and C18O (green) data. The profiles are derived from the corresponding peak intensity maps and normalised to the peak value.
In the text

	Fig. 4 Residuals of the 12CO brightness temperature after subtraction of an azimuthally averaged radial profile. Two spirals SpT1 and SpT2 are spanningan azimuth of ~100° and >180° respectively between ~10 and 180 au.
In the text

	Fig. 5 Rotation map of the three CO isotopologues and the corresponding uncertainties of 12 CO (top right panel) calculated with bettermoments. Regions below 4 σ (12 CO) and 3.5 σ (13 CO, C18 O) are masked out. The maximum and minimum velocities along the red- and blueshifted major axes are overlaid with grey lines on the 12CO rotation pattern.
In the text

	Fig. 6 Best-fit Keplerian rotation model (left panel) and the residuals calculated by subtracting the model from the observed rotation velocity (right panel), shown for 12 CO. The residuals inside a radius of 1.5 times the FWHM of the beam are masked out. In the residuals a spiral Spv is spanning more than one azimuth between ~40 and 180 au.
In the text

	Fig. 7 Deprojected and rotated (top panels) as well as polar-deprojected (bottom panels) TB residuals of 12CO with overlaid Archimedean and logarithmic spirals. The blue dashed lines show the fit of the spiral observed in the NIR by Uyama et al. (2020) extrapolated to large azimuthal angles. The solid blue line highlights the spiral region used by Uyama et al. (2020) to obtain the fit.
In the text

	Fig. 8 Deprojected and rotated (top panels) as well as polar-deprojected (bottom panels) velocity residuals of 12 CO with overlaid Archimedean and logarithmic spirals. The spiral SpT2 observed in TB is overlaid as a grey dashed line. The residuals inside a radius of 1.5 times the FWHM of the beam are masked out.
In the text

	Fig. 9 Deprojected and rotated TB (top) and velocity (bottom) residuals of 13CO (left) and C18O (right). The logarithmic spirals SpT2 and Spv observed in 12CO are overlaid.The residuals inside a radius of 1.5 times the FWHM of the beam are masked out in the velocity residuals.
In the text

	Fig. 10 Pitch angle of the three spirals observed in the TB and rotation velocity residuals of 12CO, shown for the linear and logarithmic spiral. The constant pitch angle of the NIR spiral found by Uyama et al. (2020) is included as a reference. The shaded region represents the disc regions where the linear fit fails to reproduce the morphology of spiral SpT2.
In the text

	Fig. 11 Illustration of the possible morphology of the disc around CQ Tau.
In the text

	Fig. B.1 12CO J = 2− 1 line imaged with a channel width of Δv = 0.5 km s−1 and a beam of 0.′′121 × 0.′′098.
In the text

	Fig. B.2 13CO J = 2− 1 line imaged with a channel width of Δv = 0.7 km s−1 and a beam of 0.′′128 × 0.′′103.
In the text

	Fig. B.3 C18O J = 2− 1 line imaged with a channel width of Δv = 1.0 km s−1 and a beam of 0.′′129 × 0.′′103.
In the text

	Fig. C.1 Left: peak intensity map from postprocessing the DALI model presented by Ubeira Gabellini et al. (2019), shown for 12CO. Right: azimuthal variations of the peak intensity of 12CO, 13 CO, and C18 O for an annulus of 20–40 au (0.′′12-0.′′25), derived from the DALI model and normalised to the peak value.
In the text