A dust and gas cavity in the disc around CQ Tau revealed by ALMA

ABSTRACT

1 INTRODUCTION

Protoplanetary discs are the natural outcome of the star formation process (Shu, Adams & Lizano 1987). Material infalling from the pre-stellar core is channeled into the central star and, due to angular momentum conservation, eventually forms a disc. Since the earlier stages of this process, planetary systems form from the material present in the circumstellar discs. The structure and evolution of protostellar discs are important keys in understanding the planet formation process. The main scenarios currently considered for planet formation are the core accretion model (e.g. Pollack et al. 1996) and disc instability (e.g. Boss 1997). According to the former scenario, planets form through the sequential aggregation of the solid component present in protoplanetary discs and eventually form planetesimals (e.g. review of Testi et al. 2014). The embryo accretes through a balance between collision of planetesimals and fragmentation until it has obtained most of its mass (Birnstiel, Fang & Johansen 2016). After this, a rapid gas accretion can occur, which leads to the giant planet formation. The alternative scenario is related to the development of gravitational instabilities (GI) during the initial stages of the disc evolution. In the outer part of self-gravitating discs, if the disc–star mass ratio is higher than the disc aspect ratio (M_disc/M_⋆ ≳ H/R), the rapid growth of the density perturbations induced by GI may produce bound clumps, although the effectiveness of this model to produce Jupiter mass planets is controversial (e.g. review of Kratter & Lodato 2016). Independently of the exact planet formation process, it is clear that grain growth plays an important role in the evolution and dynamics of the disc.

Young massive planets are expected to imprint signatures on the dust and gas structure of their parent protoplanetary disc, such as cavities, gaps or asymmetries, which can be detectable in the scattered light and thermal emission of discs. Thanks to the new facilities, like the Atacama Large Millimeter/submillimeter Array (ALMA), SPHERE on the Very Large Telescope (VLT), and GPI at the Gemini Telescope, such signatures are now observed with relative ease (e.g. Garufi et al. 2013; van der Marel et al. 2013; ALMA Partnership et al. 2015; Benisty et al. 2015, 2017, 2018; Andrews et al. 2016; Isella et al. 2016; Pérez et al. 2016; Fedele et al. 2017; van der Plas et al. 2017; Pohl et al. 2017; Clarke et al. 2018; Fedele et al. 2018; Hendler et al. 2018; Long et al. 2018; Dipierro et al. 2018a; Liu et al. 2019). Of particular interest are those discs whose emission shows a dip at NIR wavelengths in the Spectral Energy Distribution (SED), indicating a decrease in the NIR grain opacity. Although there is no agreed explanation of what causes the dip in the NIR, possibly this can be related to the presence of a depletion of small dust grain in the inner disc regions. The flux emission at longer wavelengths resembles the typical emission coming from a dust-rich object in the outer regions (Strom et al. 1989; Calvet et al. 2005). These discs are referred to as transitional discs (TDs) (Espaillat et al. 2014). The term ‘transitional’ indicates that these objects may be in a transition phase from optically thick gas/dust-rich discs extending inward to the stellar surface to objects where the disc has been dispersed. Whether all discs go through such a phase is however still debated. These discs are anyway excellent candidates to test planet formation theories.

Quantifying the amount of dust and gas content inside the cavity of transitional discs is of fundamental importance in order to distinguish between different clearing mechanisms. The main processes invoked so far are (1) dispersal by (photoevaporative) winds (e.g. Clarke, Gendrin & Sotomayor 2001), (2) dynamical clearing by a (sub-)stellar companion (e.g. Papaloizou et al. 2007), and (3) MHD winds and dead zones (e.g. Flock et al. 2015; Pinilla et al. 2016). Grain growth was considered as another possible mechanism (e.g. Dullemond & Dominik 2005; Ciesla 2007; Birnstiel, Andrews & Ercolano 2012b). However, recently it was found not to alter the gas profile significantly (Bruderer 2013) and to have difficulties in explaining the large observed cavities at (sub-)mm wavelengths (Birnstiel & Andrews 2014). Photoevaporation, dynamical clearing, and dead zones are the most accredited scenarios. These processes are not mutually exclusive and can probably operate at the same time (Williams & Cieza 2011; Rosotti et al. 2013).

One of the main goals of current planet-formation studies is to find planets still embedded in protoplanetary discs, in order to catch the planet formation process as it happens. A potential detection of a planet in discs with cavities may put strong constraints on formation time-scales and eventually explain characteristics of observed exoplanets (Simbulan et al. 2017). The statistics of planet candidate detection in cavities of discs is still low, but the capabilities of new instruments are beginning to provide us with some new companion candidates (e.g. Quanz et al. 2013; Reggiani et al. 2018). Recently, Keppler et al. (2018) have detected and confirmed a point source within the gap of the transition disc around PDS 70. Considering the difficulties to obtain a secure direct detection, a complementary method consists of inferring the presence of protoplanets by looking at structures in discs, and comparing them with planet–disc hydrodynamical simulations to derive constraints on planetary masses (Jin et al. 2016; Rosotti et al. 2016; Dipierro et al. 2018b; Liu et al. 2018).

Before ALMA, the spatial distribution of gas and dust were considered to be similar, and the mm-continuum data were used to trace the gas content assuming a gas-to-dust ratio of ∼100 (similar to the interstellar ratio). Recent observations, on the contrary, reveal a discrepancy between the gas and dust disc sizes. The disc emission observed at short wavelengths (or in CO emission lines) is more extended than the disc emission observed at long wavelengths (e.g. Tazzari et al. 2016) and there is now evidence for the presence of gas inside the dust cavity (Andrews et al. 2012; van der Marel et al. 2013; Bruderer et al. 2014). Moreover, dust grains have probably grown into larger aggregates and their spectral index may even vary across the disc (Rodmann et al. 2006; Guilloteau et al. 2011; Pérez et al. 2012; Birnstiel, Klahr & Ercolano 2012a; Menu et al. 2014; Pérez et al. 2015; Tazzari et al. 2016; Liu et al. 2017). This leads to different distributions of micron- and millimetre-sized grains (Garufi et al. 2013; van der Marel et al. 2013; Feldt et al. 2017; Pohl et al. 2017). In order to overcome the uncertainties in the estimate of such profiles, a direct measurement of molecular gas is necessary to determine the spatial distribution of gas and how the gas-to-dust ratio varies in protoplanetary discs (Bergin et al. 2013; Miotello, Bruderer & van Dishoeck 2014; Williams & Best 2014; Miotello et al. 2017; Zhang et al. 2017).

One of the methods currently used to probe the spatial distribution of gas is to detect the line emission of the most abundant molecular species, i.e. CO and its less abundant isotopologues. The dominant constituent of the gaseous disc, H₂, lacking a permanent electric dipole moment, is difficult to detect, whereas the second most abundant molecule, CO, can be a possible alternative to probe the gas content. ¹²CO is the most abundant isotopologue, but its low-J rotational transitions become optically thick at low column density. For this reason it is a poor tracer of the gas content, however, its emission lines are a good probe to constrain the temperature at the disc surface. Less abundant isotopologues, such as ¹³CO and, in particular, C¹⁸O have less optically thick lines, which can be used to trace the gas down to the midplane (van Zadelhoff et al. 2001; Dartois, Dutrey & Guilloteau 2003). Importantly, optically thin tracers can give a more accurate measure of the gas column density.

The main processes regulating CO abundances, freeze-out, and isotope-selective photodissociation, need to be taken into account when interpreting CO isotopologues lines (Miotello et al. 2014, 2016). Recent ALMA disc surveys (Ansdell et al. 2016; Pascucci et al. 2016; Long et al. 2017; Miotello et al. 2017) found very low CO-based gas masses, confirming the discrepancy between CO and hydrogen deuteride mass estimates in brighter discs (e.g. Bergin et al. 2013; McClure et al. 2016). Different processes sequestering elemental carbon may be at play in protoplanetary discs, affecting the CO emission significantly (Bruderer et al. 2012; Bergin et al. 2013; Favre et al. 2013; Miotello et al. 2017; Molyarova et al. 2017; Yu et al. 2017; see Section 4.3).

CQ Tau is an ideal candidate to perform a comparative analysis between observations and simulations. While it is not a transitional disc by definition (since it has NIR excess), CQ Tau is known to have a cavity at mm wavelengths in the dust continuum (Tripathi et al. 2017; Pinilla et al. 2018). In this paper, we report ALMA observations of the system in both continuum and CO integrated intensity at unprecedented angular resolution and model the gas and dust emission to determine the disc density structure. This structure is then compared to global three-dimensional hydrodynamical simulations of a disc hosting a planet. The paper is organized as follows. In Section 2 we describe the properties of the target and the ALMA observations. In Section 3 we present the method used for the modeling. In Section 4 we described the modeling results from the physical–chemical code DALI (Dust And LInes; Bruderer et al. 2012). In Section 5 we studied the possibility of the cavity to be cleared by an embedded planet through a set of hydrodynamical simulations, which allowed us to derive the planet mass and location. Moreover, we discuss the results from the hydrosimulations with the physical–chemical code DALI and compare them with other possible mechanism to open a cavity (Section 6). Finally we conclude our work with Section 7.

2 OBSERVATIONS

2.1 Target properties

CQ Tau is a variable star of the UX Ori class with a spectral type F2 (see Table 1). It is a young (∼10 Myr), nearby (162 pc, Gaia Collaboration 2016), intermediate-mass (1.67 M_⊙) star, surrounded by a massive circumstellar disc (Natta, Grinin & Mannings 2000). The continuum emission at mm-wavelength was observed with different instruments: OVRO interferometer (Mannings & Sargent 1997); the Plateau de Bure interferometer at 2.9 and 1.2 mm (Natta et al. 2000); VLA at 7 mm and 3.6 cm (Testi et al. 2001). The analysis of the SED at mm-wavelengths allows us to infer that dust has grown up to grains larger than the typical ISM size (Testi et al. 2001, 2003; Chapillon et al. 2008). The disc was already observed with the the IRAM array by Chapillon et al. (2008), who detect a weak disc emission in ¹²CO. Mendigutía et al. (2012) gave an upper limit for the accretion rate of CQ Tau equal to |$\log (\dot{M}_{\rm acc})\lt -8.3~\mathrm{M}_{\odot }$| yr⁻¹, whereas Donehew & Brittain (2011) measured |$\log (\dot{M}_{\rm acc})\sim -7~\mathrm{M}_{\odot }$| yr⁻¹. This suggests that the CQ Tau disc has a high accretion rate.

Banzatti et al. (2011) presented intermediate resolution observations taken with VLA (1.3–3.6 cm), PdBI (2.7–1.3 mm), and SMA (0.87 mm). They estimate the dust opacity index to be β ∼ 0.6 ± 0.1, assuming that the opacity follows a power-law function of the wavelength λ: κ_λ ∝ λ^−β. Combining these new data with previous measurements, the authors were able to probe the significant grain growth occurring in the disc, with the largest grains growing up to ∼cm size. Trotta et al. (2013) found also an indication of grain growth variation with radius, where larger grains were found in the inner disc, compared to the outer disc. In a recent work by Tripathi et al. (2017), the authors detected a cavity in the continuum intensity profile of CQ Tau in SMA archival data. Pinilla et al. (2018) presented CQ Tau observations with ALMA finding the presence of a dust cavity in the continuum. The disc was also observed in scattered light (PDI) by Benisty et al. (in preparation).

2.2 Data

2.2.1 Observational strategy and data reduction

We present Band 6 ALMA data of CQ Tau (RA = 05^h35^m58.46712^s; Dec. = + 24^○44′54.0864″) from three separate programs executed during Cycles 2, 4, and 5 (2013.1.00498.S, PI: L. Pérez, 2016. A.00026.S, 2017.1.01404.S, PI: L. Testi). The ALMA correlator was configured to observe simultaneously CO(2–1), ¹³CO(2–1), and C¹⁸O(2–1), as well as the nearby continuum at 1.3 mm. The array configurations used for Cycles 2, 4, and 5 are, respectively, C34.6, C40.7, and C43.8.

Data calibration was performed by ALMA following the standard Quality Assessment procedures for each of the three programmes. After downloading the data sets, we reapplied the calibration tables and extracted the CQ Tau data from each of the four data sets (the Cycle 5 observations included two separate executions). We corrected the data for the known proper motion of the star and performed one combined iteration of phase only self-calibration. The combined data sets were then imaged to produce the continuum and integrated line intensity images used in this paper. For the purpose of the analysis presented in this paper, we imaged the data using a Briggs weighting with a robust parameter of + 0.5 and a Gaussian restoring beam of 0.″15, corresponding to a spatial resolution of |${\sim }24\,$| au. This choice was made to maximize the sensitivity to the line emission and for the purpose of deriving average radial intensity profiles of the disc emission. A detailed analysis of the disc at full angular and spectral resolution will be presented in a forthcoming paper. The observation set used in this paper is shown in Fig. 1, the noise levels achieved are 30 μJy beam⁻¹ in the continuum, and 5 mJy beam⁻¹ km s⁻¹ in the line-integrated intensity maps. We expect the intensity scale of these observations to be accurate within ∼ 10 per cent, as described in the ALMA Technical Handbook.¹

Figure 1.

ALMA observation of the continuum at 1.3 mm on the upper left panel and zero-moment maps of ¹²CO (upper right), ¹³CO (lower left), and C¹⁸O (lower right) J = 2-1. On top on the ¹²CO image we overplot the contours of the continuum, respectively, at 20, 40, 80, 160, and 320 σ. The white circle in the lower right corner indicates the size of the beam: 0.″15, corresponding to |${\sim } 24\,$| au.

Open in new tabDownload slide

The line emission from the disc shows a clear rotational pattern dominated by the Keplerian motion around the central star. To compute the integrated intensity maps analyzed in this paper, we first compensate for the Keplerian velocity pattern, then we integrated the emission in each spatial pixel using an optimized range of velocities. The procedure we followed is very similar to the one presented in Yen et al. (2018) and Ansdell et al. (2018). The parameters used for the Keplerian motion subtraction compensation are stellar mass M_⋆ = 1.67 M_⊙, distance from the Sun d=162 pc, position angle of the disc major axis pa=55 deg, and inclination i = 35 deg, which optimize the subtraction of the Keplerian pattern for our data set and are consistent with the values derived by Chapillon et al. (2008) (noting the different definition of the position angle, complementary to ours).

2.2.2 Observational results

The map of the continuum emission (upper left of Fig. 1) presents a compact disc with a clear cavity. The integrated intensity ¹²CO map (upper right of Fig. 1) reveals the presence of gas inside the cavity seen in the continuum. The emission is centrally concentrated and the majority of the intensity comes from the inner part of the disc. Considering that we expect ¹²CO emission to be optically thick, the increase in emission in the inner regions of the cavity is likely due to higher temperatures, and not directly related to the gas surface density (this will be better quantified with the models described in the following section). The ¹³CO integrated intensity map (lower left panel of Fig. 1) shows a decrease of emission in the inner regions of the cavity. The C¹⁸O integrated intensity map (lower right panel of Fig. 1) shows a cavity at the same location as in ¹³CO. The continuum and CO isotoplogues radial intensity profiles, shown in Fig. 2, are obtained with a radial cut on the major axis. The difference between the profiles of the three CO isotopologues can most likely be explained by the different optical depths (e.g. Fedele et al. 2017). The error bars shown in Fig. 2 are computed as the dispersion of the measured points inside the bin.

Figure 2.

Observational profiles done with a radial cut on the major axis of the continuum (blue), ¹²CO (orange), ¹³CO (grey), and C¹⁸O (red). Each profile is normalized to the peak value. The lines are derived from the moment zero maps.

Open in new tabDownload slide

Focusing on the outer disc, the continuum extends radially up to ∼0|${^{\prime\prime}_{.}}$|5 in radius, while the ¹²CO emission extends up to ∼1″. As it was already shown for other discs (Panić et al. 2009; Andrews et al. 2012; Rosenfeld et al. 2013; Canovas et al. 2016; Ansdell et al. 2018), the gas disc extent on average seems to be a factor of 2 larger in radius with respect to the mm dust discs. The emission shows azimuthal asymmetries, in both the continuum and lines. For the purpose of this study, we limit our analysis to the general properties. A follow-up paper will be focused on them in detail.

3 METHODS

The aim of this work is to obtain a physical–chemical disc model that can simultaneously reproduce the dust and gas emission properties, with a particular focus on the disc cavity. In addition to the thermal continuum and CO isotopologues maps previously described, we additionally recover information from the SED, with data taken from previous studies, as described in the Appendix (Table A1).

3.1 DALI

In order to determine the gas and dust temperature and molecular abundances at each position in the disc, we use the physical–chemical code DALI (Bruderer et al. 2012; Bruderer 2013) that self-consistently computes the physical, thermal, and chemical disc structure. Given a density structure and a stellar spectrum as inputs, the continuum radiative transfer is solved using a Monte Carlo method to calculate the dust temperature T_dust and local continuum radiation field from UV to mm wavelengths. The chemical composition of the gas is then obtained by a chemical network calculation of all the cells. The chemical abundances enter a non-LTE excitation calculation of the main atoms and molecules. The gas temperature T_gas is then obtained from a balance between heating and cooling processes. Since both the chemistry and the molecular excitation depend on T_gas and vice versa, the problem is solved iteratively. The code is 2D in R-z. Finally, spectral image cubes are created with a ray tracer fixing the inclination of the disc in the sky at 33^○, assuming perfect azimuthal symmetry. Note that in the DALI modelling the hydrostatic equilibrium is not solved.

3.2 Gas and dust profiles

Figure 3.

Surface density profiles of the best representative model for the three components considered: gas (red line), small grains (green line), and large grains (blue line). The dashed lines show the profiles of the surface density of the disc if the cavity was not present.

Open in new tabDownload slide

Two different dust populations are considered into DALI, not only large (1 |${\mu m}$|–1 cm) but also small (0.005–1 |${\mu m}$|⁠) grains. While the radial distribution of large grains is described by equation (2), small grains follow the gas distribution presented in equation (1), with a different surface density normalization (Σ₀), taken equal to Σ_{0, small}. The size distribution of dust grains is n(a) ∝ a^−q, where q is the grain size distribution index, taken as q = 3 for large grains and as q = 3.5 for small grains, and the maximum size of grains considered is |$a_{\rm max} = 1\,$|cm. We choose a value of q = 3 for large grains after performing a test with a value of 3.5. In order to reproduce the SED profile at mm size grains, indeed, we needed to give more weight to the mm-size grains (see Fig. 4).

Figure 4.

SED profile. The black dots are the observations as described in Table A1 in the Appendix. The blue dashed line is our best representative model (see Table 2).

Open in new tabDownload slide

Both continuum and line ALMA data show the presence of a cavity in the inner regions of the disc as shown in Fig. 1. The size of the cavity and the amount of depletion are not the same for the dust and gas components. We call R_cav the gas cavity radius, while the dust cavity radius is set by the location of the Gaussian ring peak, R_dust. Within R_cav, the gas (and the small dust) surface density is lowered by a factor of δ_gas, whereas the large dust depletion is described by the Gaussian profile of equation (2), given a width of σ (see Fig. 3).

For the stellar photosphere we considered a blackbody with |$L_{*}=10\, \mathrm{L}_{\odot }$| and |$T_{*}=6900\,$| K (Table 1). We do not consider accretion on to the star, since the accretion luminosity contributes to less than |$10{{\ \rm per\ cent}}$| to the total luminosity of the star (Meeus et al. 2012) and such amount is mainly just important for the inner disc heating. The inner radius for both gas and dust in our models is set to 0.2 au, following the definition of the dust sublimation radius: R_sub = 0.07(L_*/L_⊙)^1/2 au (assuming a dust sublimation temperature of 1500 K; Dullemond & Dominik 2005).

3.3 Determining the best representative model

We start our investigation by exploring the parameters listed in Table 2, taken from the values reported in literature (e.g. Testi et al. 2003; Chapillon et al. 2008; Banzatti et al. 2011; Trotta et al. 2013). We do not carry out any χ² minimization analysis as the complexity of the DALI disc models does not allow us to cover a wide set of parameter values. Nevertheless, we explore the parameter space starting from a sparse grid of models, which is then refined around the values obtained for the best representative models. The sub-mm continuum emission helps us to fix the parameters related to the radial distribution of large grains in the disc, the size of the dust cavity, and the total dust mass. The temperature can vary, changing the vertical structure of the disc and the position or scale height of the cavity wall, which determines the amount of direct irradiation by the star that the disc can receive. Finally, the mid-IR part of the SED and the CO line radial profiles are sensitive to both density structure and temperature of the disc. In Table 2 the range of values used for each parameter is listed, where the best representative value is shown in bold face.

3.3.1 SED and radial profile of the continuum emission

Our first step is to look for a model that reproduces the SED profile between 0.36 |${\mu m}$| and 3.55 cm. The optical part of the SED (from 0.36 to 3.4 |${\mu m}$|⁠) is affected by extinction and it is corrected considering the extinction law of Cardelli, Clayton & Mathis (1989), with R_v = 3.1. In order to match the stellar spectrum, an extinction of |$A_{\rm v} = 1.9\,$|mag is required, consistent with previous works (e.g. Garcia Lopez et al. 2006). In our models the vertical structure (h_c, ψ) is the main factor that can introduce a shadowing and lower the amount of exposed disc surface accessible to stellar light (Woitke et al. 2016). These are key parameters for the final temperature of the disc and to reproduce the observed NIR and FIR excess.

In order to describe the mm part of the SED, we vary R_dust, σ, and Σ_{0, dust}, until we match simultaneously the radial profile of the ALMA continuum observations and the longer wavelengths points in the SED.

3.3.2 CO emission radial profiles

The following step is to reproduce the gas profiles of the three CO isotopologues. The CO isotopologues radial profiles are more radially extended than the continuum emission (as it is possible to see from the normalized profiles in Fig. 2). As many previous studies suggest, a single surface density profile for both gas and dust is not able to reproduce the continuum emission and the CO integrated intensity maps simultaneously, indicating that optical depth effects are not sufficient to explain the different radial extent (e.g. Facchini et al. 2017) as it is shown in the DALI modelling described in Section 4.

4 MODEL RESULTS

4.1 SED best representative model

In order to reproduce the SED profile (Fig. 4), we follow the approach described in Section 3.3.1. The two main observational constraints to set the best representative model for the SED are the near-infrared (NIR) and far-infrared (FIR) excesses that give information about the vertical structure of the disc, and the total amount of small dust present in the disc. We vary the values of h_c and ψ, which describe how much the disc is vertically extended and flared. In particular, the parameter that mostly affects the SED is h_c (Woitke et al. 2016). Our best representative values for these parameters are h_c = 0.125 and ψ = 0.05. The range of values under which our models match the observed SED are 0.1 < h_c < 0.15 and 0.05 < ψ < 0.1. We note that a variation in the stellar luminosity value can affect the choice of such parameters. Longer wavelengths are affected mainly by the total amount of mass present in the disc. Once the correct amount of total mass is found, we fix this value and vary the other parameters in order to find a good agreement with the ALMA continuum profile. As it is possible to see in Fig. 4, the observations are well described by our best representative model. Note that model and data diverge at the wavelength of 3.57 cm. This may be explained by free–free emission, which we do not model, and affects the integrated flux at cm wavelengths. We note also that the SED shows a NIR excess together with a resolved cavity in both dust and gas seen in the ALMA data. The small grains are the main responsible for the emission at such wavelengths and a small amount in the cavity is enough to reproduce the NIR emission.

4.2 Radial profiles

In order to compare our models with our ALMA observations, we convolve the model results with a Gaussian beam of of 0.″15, which simulates the resolution of the observations, both for the continuum and for the CO isotopologues. We perform a radial cut along the major axis of the continuum and moment zero maps, which has a better resolution compared to an azimuthally averaged profile, in order to compare with our 2D models.

4.2.1 A Gaussian ring for large grains

Inspired by the work of Pinilla et al. (2018), we have employed a Gaussian functional profile, as previously discussed (see equation 2). As best representative parameters, we find a centre peak position of |$R_{\rm dust}=53\,$| au and a Gaussian width of |$\sigma =13\,$| au. Our results are consistent with the Pinilla et al. (2018) work where the best model had a Gaussian peak radius of |${\sim } \, 46\,$| au and, respectively, an inner and outer width of |${\sim } \, 11\,$| and |${\sim } \, 17\,$| au. The peak of the Gaussian is |$\Sigma _{\rm 0, dust}=0.6\, {\rm g\, cm^{-2}}$|⁠. The simulated continuum emission profile is plotted in Fig. 5 (blue line) along with the data points (black line). The data are always matched by our model within the error bars, and from now on we will refer to this as the best representative model for the large dust component.

Figure 5.

Dust continuum taken at 1.3 mm with ALMA. The radial cut on the major axis of observation (black line) is shown along with our best representative model (blue line): Gaussian profile.

Open in new tabDownload slide

4.2.2 Gas radial profile

We use the CO isotopologues observations as a proxy for the gas distribution. In Fig. 6 we plot the ¹²CO, ¹³CO, and C¹⁸O observed radial profiles, along with a selection of our DALI models. In the disc inner regions, where CO lines are more optically thick, the gas temperature plays a big role and strongly depends on the small grains distribution. We implement a relatively high surface density normalization of small grains Σ_{0, small} = 0.0375 g cm⁻². Setting a higher column density of small grains, while keeping fixed the amount of gas, shifts the dust τ = 1 layer upwards and, consequently, a more significant amount of UV photons is absorbed there. The location of the CO isotopologues τ = 1 layer does not change significantly but the gas temperature at the emitting layer is reduced as the UV flux is reduced. Moreover, a smaller amount of UV flux can produce a decrease of the dust temperature and consequently, due to thermal balance, can contribute to a lower gas temperature. Accordingly, simulated integrated intensities are lower. A smaller Σ_{0, small} would lead to higher integrated intensities than what is shown in Fig. 6, with a worst match with the observational data. All CO isotopologues model results shown in Fig. 6 are obtained assuming a gas profile (equation 1) described by the following parameters: γ_gas = 0.3, Σ_{0, gas} = 2.5 g cm⁻², and a cut-off radius of R₀ = 56 au. The models always overpredict the emission in the outer disc (¹³CO in particular). Our focus is however on the cavity, as we are interested in constraining which may be the clearing mechanism. Future studies will attempt to give a better description of the CQ Tau outer disc.

Figure 6.

ALMA observations radial cut on the major axis (black line) for gas in the three different isotopologues of CO (top: ¹²CO, middle: ¹³CO, bottom: C¹⁸O). Along with our best representative model (red line), on the right we show additional models varying R_cav, while on the left we vary δ_gas. Moreover, we present an additional model (orange) that still describes well the ¹²CO profile and it is useful for the model results (see below).

Open in new tabDownload slide

We explored different values for the size and depletion of the cavity (respectively, R_cav and δ_gas) to estimate the best representative model and to gain insights about their degeneracies (Fig. 6). Observations agreed with a gas depletion factor of 10⁻³ <δ_gas < 10⁻¹. The top left panel of Fig. 6 shows that the ¹²CO profile cannot be described by a model with a depletion factor of 10⁻³, otherwise the ¹²CO emission would also show a depression at small radii. Furthermore, ¹³CO and C¹⁸O profiles suggest that the depletion factor has to be smaller than δ_gas ≲ 10⁻¹ (models shown in blue and yellow do not match the data). Finally, the cavity size is constrained to be larger than 15 au and smaller than 25 au. Our best representative model is shown with a red line in Fig. 6 and assumes R_cav = 20 au and δ_gas = 10⁻². The gas has a smaller cavity radius (⁠|$R_{\rm cav}=20\,$| au) than the large dust grains (⁠|$R_{\rm dust}=53\,$| au). The depletion in gas within the cavity is found to be similar to the one of the large grains (see Fig. 3).

4.3 Dust-to-gas ratio

Fig. 7 shows the cumulative function of mass as a function of the radial distance from the star of both the gas (red) and the dust (blue, small, and large grains together) for our representative model. The dust-to-gas ratio is not uniform across the disc, and is close to 0.1. We compute the total mass for the three different components of our best representative model: M_gas = 3, 2 × 10⁻³M_⊙; |$M_{\rm small \, dust} = 4.9 \times 10^{-5} \mathrm{M}_{\odot }$|⁠; |$M_{\rm large \, dust} = 2.3 \times 10^{-4} \mathrm{M}_{\odot }$|⁠. The global dust-to-gas ratio for this model of CQ Tau is ∼0.09. This value is higher than the typical ISM values of ∼0.01, but it is consistent with recent findings from ALMA surveys of protoplanetary discs in different star-forming regions (Ansdell et al. 2016, 2017; Pascucci et al. 2016; Long et al. 2017; Miotello et al. 2017). The CO isotopologues emission lines are found to be generally faint compared with the disc continuum emission, resulting in low CO-based gas masses and high CO-based dust-to-gas mass ratios. This may be interpreted as rapid disc physical evolution via loss of gas, or alternatively, as chemical evolution due to carbon depletion processes. The latter may occur if carbon is sequestered from CO and it is either locked into large icy bodies in the mid-plane or converted into more complex and less volatile molecules (Aikawa et al. 1996; Bergin et al. 2014; Du, Bergin & Hogerheijde 2015; Eistrup, Walsh & van Dishoeck 2016; Kama et al. 2016; Yu et al. 2016). Current CO and continuum ALMA observations alone cannot distinguish between the two scenarios, but high dust-to-gas mass ratios hint at disc evolution, either physical or chemical. Independent gas mass measurements of discs have been obtained thanks to the detection on hydrogen deuteride with the PACS instrument on the Herschel Space Telescope in few discs (Bergin et al. 2013; McClure et al. 2016). When compared with CO isotopologues observations, such data show that volatile carbon needs to be depleted by different factors depending on the source (Favre et al. 2013; McClure et al. 2016). Unfortunately, no facility is available today to reproduce this kind of observations. Alternatives to constrain the level of volatile carbon and oxygen depletion in discs come from the observations of atomic recombination lines, such as [CI] (Kama et al. 2016) or from the detection of hydrocarbon lines (Bergin et al. 2016; Cleeves et al. 2018).

Figure 7.

Radial profile of the cumulative function of mass of our best representative model (as shown in Fig. 3). The red line shows the gas behaviour, while the blue line shows the dust component considering both the small and large grains; the dashed line shows the cumulative mass function without considering the presence of a cavity for both gas (red) and dust (blue). The global dust-to-gas ratio is ∼0.09.

Open in new tabDownload slide

4.4 A further test for large dust grains: self-similar density profile

In order to model the large dust grains density profile, we also attempted to use a self-similar density profile as done for the gas (see equation 1). In this case, we used the same cut-off radius as for the gas (R₀ = 56 au), but with a different γ_dust and surface density normalization of large dust grains Σ_{0, dust}, allowed to vary independently. The best parameters for large dust grains are found to simultaneously reproduce both the observations of the ALMA continuum radial profile and the SED. The power-law index of the surface density distribution is γ_dust = −0.7 and the surface density normalization is Σ_{0, dust} = 0.85 g cm⁻².

To model the inner region of the disc, we consider that the size of the cavity and the amount of depletion can be different between the dust and the gas. We choose as dust cavity radius R_{cav, dust} and as dust cavity depletion δ_dust that are allowed to be different from the gas ones. The dust cavity has a radius of 28 au with a depletion of a factor δ_dust = 10⁻². The proximity of the cavity edge with the cut-off radius implies that the dust is concentrated in a narrow ring.

The self-similar density profile reproduces the observed SED (it overlaps the Gaussian ring SED, blue curve in Fig. 4). The resulting continuum emission profile is plotted in grey in Fig. 8 along with the data points (black line). From this plot it is possible to see that the self-similar density profile still describes well the large grains behaviour.

Figure 8.

Dust continuum taken at 1.3 mm with ALMA. The radial cut on the major axis of observation (black line) is shown along with our best model obtained assuming a self-similar density profile of large dust grains that still can nicely describe our data (see Section 4.4). The dust cavity radius is at 28 au and the depletion factor is of δ_dustcav = 10⁻².

Open in new tabDownload slide

5 IS THE CAVITY INDUCED BY AN EMBEDDED PLANET?

The cavity visible in the CQ Tau disc can be explained by several mechanisms: dispersal by (photoevaporative) winds, dynamical clearing by a (sub-)stellar companion, MHD winds, and dead zones. In this section, we focus on the study of the dynamical clearing induced by a companion using both analytical considerations (Section 5.1) and hydrodynamical simulations (Section 5.2). We note that the profiles used for the fitting procedure with DALI are based on simplified physical models that qualitatively represent the characteristic features of the real profile. As a consequence, performing hydrodynamical simulations becomes of paramount importance in order to verify that the DALI models can be qualitatively reproduced in hydrodynamical models. This mechanism, indeed, leads to a hole/cavity empty of large grains that are filtered out, while the gas and small dust grains can still continue to be accreted by the central star (Rice et al. 2006; Zhu et al. 2012).

5.1 Analytical considerations

We explore the possibility that a massive planet located in the cavity region at a separation R_p ≈ 15−25 au is responsible for the gas and dust density structure we observe.

5.2 Numerical simulations

We perform a set of 3D numerical simulations of a planet and a gas + dust accretion disc using the code phantom (Price et al. 2018) in order to find the mass of the planet producing the density structure expected from the DALI modeling of the observations.

We use a locally isothermal equation of state: the gas temperature is a radial power-law that produces a disc thickness H = h_c(R/R₀)^ψR (where ψ = 0.05; bold parameter in Table 2). The planet and the star are modelled by two sink particles (Bate, Bonnell & Price 1995). The planet is initially on a circular orbit and is allowed to accrete and migrate. The sink particles exert the standard gravitational force on gas and dust particles. The motion of the sinks is computed step by step from the force they exert on each other and from the back-reaction of the disc. Particles are considered accreted when their distance from the sink is below the sink radius (see that next section), and they are gravitationally bound to it.

The profile of ‘large’ dust grains in the DALI modeling describes the properties of a vast population of different grain sizes (1μm ≲ a_large ≲ 1cm). In our simulations however, we do not follow a dust population with a grain size distribution, but we only model individual sizes.

To describe the dynamics of dust grains, we adopt the one fluid algorithm developed in Laibe & Price (2014), Price & Laibe (2015), and Ballabio et al. (2018) for large grains (for which we considered two different sizes: a_large = 0.1mm and a_large = 0.5mm). All grain sizes we used are characterized by a Stokes Number St < 1 throughout the entire disc, so that the dust dynamics is correctly described by the one fluid SPH algorithm (Ballabio et al. 2018).

We performed some simulations using also an additional family of dust grains of size |$a_{\rm small}=1\, \mu m$|⁠, with the aim to investigate the dynamics of small dust particles that are part of the modeling with DALI. We notice that, as expected and previously discussed in Section 4.2, very small dust grains (⁠|$a_{\rm small}=1\, {\rm \mu m}$|⁠) reproduce closely the dynamics of the gas, being tightly coupled to it. As a consequence, we choose not to investigate the dynamics of |$a_{\rm small}=1 \, {\rm \mu m}$| grains since the information about their distribution can be inferred by merely rescaling the gas profile with the appropriate dust-to-gas ratio.

Angular momentum transfer, allowing gas accretion, is provided by SPH artificial viscosity. We choose the artificial viscosity parameter α_AV = 0.2 in order to have an equivalent Shakura & Sunyaev (1973) turbulent parameter α_SS ≈ 0.005−0.01 throughout the entire disc. Our choice of α_SS is dictated by the minimum viscosity required in an SPH numerical simulation in order to properly provide shock resolution.

Each simulation is evolved for |$t=60\, t_{\rm orb}$|⁠, which corresponds to ∼4500 yr. We run all our simulations with N_part = 7.5 × 10⁵ particles. We also run two simulations using a larger number of SPH particles [N_part = (1.5–2.2) × 10⁶] in order to check the reliability of our results at lower resolution.

5.2.1 Initial conditions

The initial conditions of our simulations consist of a planet with mass ranging M_p = 3–9M_J at a separation from the central star of |$R_{\rm p}=20\, {\rm au}$|⁠. The choice of the planet mass and its separation from the star have been made starting from the considerations in equation (5).

We run simulations using three different couples (M_p, h_c) for which the gap opening criterion in equation (5) is satisfied. We note that the disc aspect-ratio h is poorly constrained from observations considering that it is strongly affected by the estimate of the luminosity of the star (as previously discussed).

The planet is initially completely embedded in the disc, that is assumed to have an unperturbed density profile (no cavity is present at the beginning of the simulation neither in the gas nor in the dust density profiles). The gas initial surface density profile uses the parameters obtained from the power-law DALI density model (bold parameters in Table 2, i.e. dashed red curve in Fig. 3, respectively) without depletion. The large dust grains initial density profile uses for simplicity the self-similar DALI model introduced in Section 4.4, but using a tapering radius larger than R₀ (⁠|$R_{\rm taper}=70\, {\rm au}$|⁠): the steep profile at large radii observed in the large dust grain surface density profile from the DALI modeling is the natural outcome of large dust grains radial drift; as the simulation starts, the large grains drift inward and the dust density profile steepens. This demonstrates that the narrow dust ring observed with ALMA is most likely due to radial drift of large dust grains.

It should be noted that thicker discs are more efficient in transporting the angular momentum throughout the disc. As a consequence, as can be seen in equation (4), for a fixed mass of the planet, the gap opening in the gas becomes progressively less efficient as the disc aspect-ratio grows. This, in fact, introduces a degeneracy for the couples of parameters (M_p, h), that is confirmed by the results of our numerical simulations.

In Table 3 we report a summary of the simulations that better reproduce the profiles from the DALI modeling. We note that our high-resolution runs (simulations 7 and 8 in Table 3) show overall good agreement with lower resolution runs except in the very inner region of the disc (see the end of Section 5.2.2).

5.2.2 Results

In all numerical simulations reported in Table 3 the planet rapidly (⁠|$t\approx 15\, t_{\rm orb}$|⁠) carves a gap in the initially unperturbed disc (both in the gas and in the dust), as expected from the gap opening criterion presented in equation (5).

Fig. 9 shows the results of our numerical simulations, from left to right with M_p = 3, 6, and 9 M_J, respectively, for the various species. The first and second rows show the azimuthally averaged density profiles for gas (orange solid lines), for 0.5-mm-sized dust (green solid lines), and for 0.1-mm-sized dust (blue solid lines). For a direct comparison, the same plots also show the surface density profiles obtained from the DALI modeling. The red dashed line is the reference model for the gas (boldface parameters in Table 2), which has a gas depletion factor of δ_gas = 10⁻², while the orange dotted line is the same model but with a smaller gas depletion in the cavity, δ_gas = 10⁻¹, and larger cavity size, |$R_{\rm cav}=25\, {\rm au}$|⁠. These gas profiles produce the red and orange emission profiles in Fig. 6, respectively. The blue dashed line shows the reference Gaussian ring density profile for the ‘large dust grains’ obtained from DALI, which produces the blue curve in the continuum emission profile in Fig. 5; the grey dotted curve is the best-matching power-law density profile for the large dust grains obtained with DALI and discussed in Section 4.4, the dust continuum emission associated to it is the grey curve in Fig. 8. The third and fourth rows of Fig. 9 show the surface density colour plots of the corresponding simulations in the upper panels (Refs 1, 2, and 3 in Table 3).

Figure 9.

Profiles (top row) and colour plots (central and bottom row) of the surface density of both gas and dust. Different columns show the results of simulations for different planet masses M_p and aspect ratios h_c. The first row and the second row show the azimuthally averaged surface density of gas and dust, respectively. In particular from our simulations of gas (orange solid curve), large dust grains |$a_{\rm large}=0.1 \, {\rm mm}$| (blue solid curve, from left to right simulations 1, 2, 3, respectively) and |$a_{\rm large}=0.5 \, {\rm mm}$| (green solid curve, from left to right simulations 4, 5, and 6, respectively) after |$t\approx 60\, t_{\rm orb}$| of evolution, to be compared to the surface density profiles obtained from the DALI modeling for the gas (red dashed and orange dotted curves, producing the CO isotopologues emission in Fig. 6 with the same colours) and for the large dust grains (blue dashed and grey dotted curves producing the continuum emission profiles in Fig. 5 with the same colours). Dashed curves represent the best denisty models from DALI, while dotted ones are those DALI models that best fit the simulations. The third row shows the surface density colour-plot of the gas (from simulations 1, 2, and 3, as indicated in the left top corner of the images). The fourth row shows the surface density colour plot of the large dust grains |$a_{\rm large}=0.1 \, {\rm mm}$| (from simulations 1, 2, and 3, as indicated in the left top corner of the images). The colour scale is logarithmic; the large white dot represents the central star and the small one shows the planet position.

Open in new tabDownload slide

In Table 3 we report also two convergence test runs using a larger number of particles and the same parameter choice as our fiducial model (N_part = 1.5 × 10⁶ and N_part = 2.2 × 10⁶ simulations 7 and 8, respectively). These high-resolution runs are in good agreement with the lower resolution one throughout the entire disc for both gas and dust, apart from the inner disc region where at low resolution the lower density produces an increase of SPH artificial viscosity, and thus a spurious fast evacuation of the innermost part of the cavity.

We finally note that the mass of the planet does not change significantly during our simulations. In particular, in all cases, we observe a growth of the planet mass at the end of the simulation of a fraction of M_J. The correspondent growth rate at the end of the simulation is of |$\dot{M}_{\rm p}\sim 5\times 10^{-5} \, {\rm M_{\rm J}\, yr^{-1}}$|⁠, consistent with the expected theoretical prediction for the planet accretion rate (see equation 15 of D’Angelo & Lubow 2008).

5.2.3 Radiative transfer modelling

We have processed our hydrodinamical model results, using a radiative transfer tool (RADMC3D; Dullemond et al. 2012) for the dust component, in order to test if our SPH model is qualitatively consistent with the observations. The images obtained were, then, convolved with a Gaussian beam of 0.″15 as done in DALI and compared with the observational data. In particular, we found the best agreement with the data in the model with the planet at R_p = 20 au, M_comp = 6 M_Jup and the size of large grains of 0.1 mm (see Fig. 10). We show the observations at 1.3 mm (black line) along with the synthetic dust flux profiles convolved based on radiative transfer calculation. The plots show the model presented in Section 5.2.2 in green considering also the inner disc component; the model in yellow, instead, does not take into account the presence of inner disc and uses a cavity lowered of the same factor from the planetary radius up to the sublimation radius. Our simulation (both yellow and green lines) well describe the outer disc, which reaches a semistable state in few dynamical time-scales, whereas the inner disc is not well reproduced by the hydrodynamical simulations. Its flux is, indeed, higher than the one of observations, leading to a disc where the carved cavity is not well seen (see Fig. 10). We expect, however, the inner disc within the planetary orbit to drain with time on a longer time-scale than the outer disc and this might affect the dust profile in the inner regions. With a lower inner disc emission, we would be able to reproduce qualitatively the depth and separation of the planet (yellow line in Fig. 10), main goal of hydromodelling we provided. A more detailed analysis on the inner disc will be addressed in future works.

Figure 10.

Radial profiles of the dust component of the hydrodinamical simulations after radiative transfer modelling. The black dots are our data in the continuum (1.3 mm); the green and yellow lines are, respectively, the model with and without the presence of the inner disc. The best model here listed is using R_p = 20 au, M_comp = 6 M_Jup. The size of large grains is taken as 0.1 mm.

Open in new tabDownload slide

6 DISCUSSION

We here investigate if a planet could be responsible for the formation of a cavity in the gas and dust profile similar to the one found through the chemical–physical code DALI. We note the following:

• A planet mass of the order of 3M_J does not reproduce the strong gas depletion in the inner disc inferred from the thermo-chemical modeling with DALI. On the other hand, a larger planet mass (M_] ≈ 6 − 9M_J) does a much better job, although the 6M_J case requires a slightly thinner disc (h = 0.1) than that obtained from the DALI modeling.

• Overall, we find that both our gas and dust density profiles qualitatively reproduce the ones inferred from thermo-chemical modeling of the observed fluxes. However, we note two issues. First, concerning the gas, the best representative DALI models predict a gas depletion factor of the order of 10⁻², which is not achieved in our simulations, that only show a gas depletion of 10⁻¹ (see red and orange curves in CO isotopologues DALI model Fig. 6 and simulations Fig. 9). Such a high gas depletion is essentially required by the observed gap in ¹³CO. A higher gas depletion might be achieved using smaller values of α_SS and h_c, or much larger values of the planet mass. The latter possibility is unlikely at this planetary location, since it would imply that the planet mass falls within the brown dwarf mass regime. Direct imaging observations, indeed, ruled out the presence of a companion with mass higher than 10 M_Jup beyond 20 au (Benisty et al., in preparation).

• Concerning the dust, our hydrodynamical model involving one planet reproduces well the ring feature at ≈50 au, but not the inner region of the cavity, where an inner dust disc close to the star still remains. While a semistationary state has been reached outside the planet orbit (i.e. we do not expect the gap edge location, depletion at the planet location and outer density profile to change significantly at longer time-scales), we note that the dust in the inner disc decouples from the gas and is not effectively accreted on to the central star along with the gas. This implies that further depletion of the dust in the inner cavity might occur at later times.

• Good agreement between the outcome of our simulations and the DALI gas density profile can be obtained with the model having |$\delta _{\rm gas}=10^{-1},\, R_{\rm cav}=25$|⁠; in particular, we note that the location of the cavity edge for the gas in this model (‘orange’ CO isotopologues DALI model in Fig. 6) is perfectly consistent with the density structure provided by a planet placed at |$R_{\rm p}=20\, {\rm au}$| in our simulations. This DALI model, despite not being the best representative for the gas, reproduces fairly well the fluxes for ¹²CO and C¹⁸O, but it fails in reproducing the ¹³CO.

• The dust profiles from the numerical simulations appear to best reproduce the DALI self-similar model (Section 4.4) with δ_dust = 10⁻² and R_cav = 28 au (grey model in Figs 8 and 9). The hydrodynamic density profiles, in particular those with a_large = 0.1 mm, reproduce the location of the cavity edge and depth of the large dust grains density distribution of this model (see blue solid curves compared to grey dotted ones in Fig. 9). This does not occur for larger grains (⁠|$a_{\rm large}=0.5\, {\rm mm}$|⁠) that instead produce a deeper gap with respect to that in the modeling, as shown in the second row of Fig. 9.

• Our hydrodynamical simulations are not able to reproduce the Gaussian density profile of the large dust grains that provides the best match with the continuum flux (blue curves in Figs 5 and 9). We note that, in general, any planet model would predict a sharp disc truncation beyond the planet location instead of a shallow Gaussian depletion. As we show in Section 5.2.3, a sharp density feature is smoothed when synthetic dust flux profile have been created and then convolved. The results found can qualitatively describe the observations as well as the Gaussian profile if the inner disc present in the simulations is neglected.

• We find that the best match between our hydrodynamical simulations and DALI density profiles is provided by grains with a_large = 0.1 mm with respect to those with a_large = 0.5 mm; this suggests that the dynamics of the entire large dust grain population can be best described by the smaller (≈0.1 mm) sizes, rather than the larger ones (≈0.5 mm).

  A more accurate treatment would require a multigrain approach (Dipierro et al. 2018b; Hutchison, Price & Laibe 2018), in order to account for the different dynamics of different grain sizes characterized by different Stokes numbers.

6.1 Comparison with other possible scenarios

We here briefly discuss our results to explain the cavity together with other possible mechanisms.

Clearing mechanism. We find that a planet |$M_{\rm p}=9\, M_{\rm J}$| located at |$R_{\rm p}=20\, {\rm au}$| from the central star qualitatively (but not perfectly) reproduces the DALI modeling of gas and the large dust density distribution. The best match is obtained with |$0.1\, {\rm mm}$| grains, implying that large dust grains behave qualitatively as fluid composed by grains with that size. In particular, good agreement is obtained at large radii, where the simulations clearly show that the planet might produce a ring-like feature in the dust distribution. However, we were not able to reproduce the dust depletion, expected from the DALI modeling, in the inner cavity using a one-planet model, and at the end of our simulations we are still left with an inner disc of dust around the star.

Photoevaporation. Models of photoevaporation predict small cavities (≲10 au) and low accretion rates. The debate about what is the driving mechanism of photoevaporation between FUV, EUV, X-ray, and what are the exact values of mass-loss rate is still open. In order for the photoevaporation to start the accretion rate have to drop below the photoevaporation rate. As soon as it starts, a cavity free from gas and dust is generated, with the very large grains rapidly migrating inward (Alexander & Armitage 2007). CQ Tau was considered by Donehew & Brittain (2011) and Mendigutía et al. (2012) to have a high accretion rate (⁠|$\dot{M}_{\rm acc}\lesssim 10^{-7}~\mathrm{M}_{\odot }$| yr⁻¹), which can be considered high for disc dispersal to be the main mechanism acting in the formation of the CQ Tau cavity (Owen 2016). A combination between the clearing mechanism and photoevaporation can, instead, be considered as a possible option (e.g. Williams & Cieza 2011; Rosotti et al. 2013; Rosotti, Ercolano & Owen 2015) to be further investigated.

Dead zones. The suppression of magnetorotational instability inside the disc can create regions of low disc ionization, the so-called ‘dead zones’ (e.g. Flock, Henning & Klahr 2012). In these regions the rate of gas flow decreases with an accumulation of gas in the outer edge. At the same time the pressure maxima is able to trap large particles with a consequent ring-like feature in both gas and dust components (Flock et al. 2015). The ring is in general located at about the same radius (Pinilla et al. 2016) in both the components. However, when dead zones and MHD wind are present at the same time, a difference in the inner radii position can be observed. The presence of dead zones can increase the turbulence present in the disc and generate asymmetric vortices in the disc (Ruge et al. 2016) at (sub)mm wavelengths. Such structures can be also created by disc–planet interaction if the α viscosity is low. We cannot exclude this mechanism by our observations. Higher resolution data are need to better constrain the profile of the ring-like disc, mainly its symmetry and variability (Pinilla et al. 2016). However, dead zones can explain ring-like structures, but not cavities such as the one clearly present in the CQ Tau system.

A combination between the mechanism described above can be also a possible solution, which, however, we will not address in this paper.

7 CONCLUSIONS

Discs where it is possible to clearly distinguish a cavity are of particular interest in the study of protoplanetary disc evolution and of the planet formation process. In particular, we have studied the disc around the CQ Tau pre-main-sequence star. We have made use of spatially resolved ALMA observations of the dust continuum and of three CO isotopologues, ¹²CO, ¹³CO, and C¹⁸O, together with the observed SED of CQ Tau. We employed the chemical–physical code DALI to model the CQ Tau disc and to derive the surface density profile of the dust (small and large grains separately) and of the gas. Finally, we ran 3D SPH simulations to test whether an embedded planet could be able to create the observed disc structure. The main results from our analysis are summarized here:

• In the CQ Tau disc the gas radial extent is a factor of 2 broader than that one of the large dust grains;

• CQ Tau shows a clear cavity in both gas and dust, which can be described by different functional forms. Assuming that the mm-sized grains are distributed in a Gaussian ring, the Gaussian peak position is located at a distance of 53 au from the central star and the Gaussian width is |$\sigma =13\,$|au. The cavity is present also in the gas component. The cavity radius and level of depletion are degenerate, but they can be constrained with reasonable uncertainties: the radius is smaller than 25 au and larger than 15 au and the depletion has to be between 10⁻¹ and 10⁻³ with respect to the profile of a full self-similar disc.

• The computed dust-to-gas ratio is not radially constant throughout the disc. Moreover, the global dust-to-gas ratio is found to be ∼0.09, higher than the typical values assumed (∼0.01) and most likely because of carbon depletion.

• We performed a set of SPH simulations in order to investigate whether the presence of a massive planet might be responsible for the radial structure of the CQ Tau disc. We find that a massive planet with a mass of at least |$M_{\rm p}=6 - 9\, M_{\rm J}$| located at |$R_{\rm p}=20\, {\rm au}$| from the central star can reproduce well the ring feature at ≈50 au but not the dust distribution in the inner cavity region of the best representative density model predicted using DALI. Also, our simulations do not reproduce the strong gas depletion within the cavity that appears from the best-fitting DALI modeling, although we caution the reader of possible degeneracies in both the thermochemical modeling (see above) and the hydrodynamical modeling (for which a lower viscosity in the disc might reconcile the results).

ACKNOWLEDGEMENTS

We thank the anonymous referee for providing insightful comments. We are thankful to Antonella Natta, Antonio Garufi, Mario van den Ancker, and Carlo F. Manara for useful discussions. 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 of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, auI/NRAO and NAOJ. This work was partly supported by the Italian Ministero dellÍstruzione, Università e Ricerca through the grant Progetti Premiali 2012 – iALMA (CUP C52I13000140001), by the Deutsche Forschungs-gemeinschaft (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). AM and SF acknowledge an ESO Fellowship. GD and ER acknowledge financial support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 681601). MT has been supported by the DISCSIM project, grant agreement 341137 funded by the European Research Council under ERC-2013-ADG and by the UK Science and Technology research Council (STFC). GL, MGUG, AM, SF, LT, LP, and MT have received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 823823 (RISE DUSTBUSTERS project).

Footnotes

REFERENCES