© ESO 2020

1 Introduction

With modern telescopes such as the Atacama Large Millimeter Array (ALMA) at millimetre-wavelengths or the SPHERE (Spectro-Polarimetric High-contrast Exoplanet REsearch) ins- trument at the Very Large Telescope (VLT/SPHERE) at optical wavelengths, it has become possible to produce spatially resolved images of planet-forming disks surrounding young low and intermediate mass stars down to spatial scales of a few astronomical units. These observations have revealed structures such as spiral arms, vortices, and, most prominently, gaps and rings in dust continuum observations. Large high-spatial-resolution ALMA surveys indicate that dust gaps and rings are common in observed planet-forming disks (Andrews et al. 2018; Long et al. 2018). Compared to ALMA, rings and gaps are not as frequently observed in scattered light images, which might be caused by observational biases. VLT/SPHERE scattered light observations show gaps and rings mostly in bright and extended disks (Garufi et al. 2018), furthermore shallow gaps are likely to be harder to detect in the optical compared to the millimetre regime due to optical depth effects. However, considering those factors, Garufi et al. (2018) concluded that scattered light observations also indicate that rings and gaps in disks are a normality.

The most common explanation of the observed dust gaps are forming planets that carve gaps and cause a pile up of larger dust particles in pressure bumps, which are observed as mostly azimuthally symmetric ring structures (e.g. Pinilla et al. 2012; Dipierro et al. 2015; Zhang et al. 2018). However, other interpretations are also possible such as enhanced pebble growth at molecular ice lines (e.g. Pinilla et al. 2017; Zhang et al. 2015; Owen 2020), dust evolution and inefficient fragmentation of dust particles (gaps in scattered light, Birnstiel et al. 2015), or magnetized disks (e.g. Flock et al. 2015; Béthune et al. 2017; Riols et al. 2020). Larger sample studies of continuum observations do not exclude any of those scenarios, however, the ice-line scenario is unlikely to be a universal mechanism as most observed gaps and rings are not associated with the expected locations of molecular ice lines (e.g. Huang et al. 2018; Long et al. 2018; van der Marel et al. 2019). The various gap formation theories differ especially in their predictions for the gas. For example in the planet scenario a dust gap is associated with a gas gap, but this is not the case for the ice line scenario or for dust gaps produced at the edge of dead zones in magnetized disks. A single planet might also produce multiple gaps in the gas and dust as shown for example by Zhang et al. (2018). However, the models of Ziampras et al. (2020) indicate that radiative effects can suppress multiple gap formation (i.e. each gap requires a planet). Accurate gas and dust surface density measurements of planet-forming disks that allow us to estimate the gap depths and shapes are required to provide constraints for the various gap formation scenarios.

In a different manner to the dust, gas observations do not provide a clear picture of the presence of gaps and rings in the gas disk. One reason for this is that spectral line observations are more expensive in terms of observing time compared to continuum observations. Furthermore the interpretation of molecular line data is more complex due to thermo-chemical and kinematic effects. For example ring structures in molecular line observations can be purely caused by chemical effects such as the freeze-out of molecules, photo-desorption, or changing dust properties (e.g. Öberg et al. 2015 molecule: DCO⁺; Cleeves 2016 CO; Bergin et al. 2016 C₂H; Salinas et al. 2017 DCO⁺, N₂ D⁺; Cazzoletti et al. 2018 CN; Qi et al. 2019 N₂ H⁺). This complexity allows for different interpretations of the same gas observations as for example discussed by van der Marel et al. (2018) for the case of HD 163296. So far indications of the presence of gas gaps that might have been produced by planets exist for only a few disks: e.g. HL Tau (Yen et al. 2016); HD 163296 (Isella et al. 2016; Teague et al. 2018); AS 209 (Favre et al. 2019); TW Hya (Teague et al. 2017).

Not included in the above list are transitional disks with their large (r ≳ 20 au), strongly dust depleted (often no dust is detected at millimetre wavelengths) inner cavities. Observations of transition disks clearly show strong gas depletion (by factors of ≳100) within their dust cavities. However, extensive studies and thermo-chemical models have shown that the gas cavity is usually smaller than the dust cavity and also that the gas is not as strongly depleted as the dust (e.g. Bruderer et al. 2014; Carmona et al. 2014; van der Marel et al. 2016; Drabek-Maunder et al. 2016; Dong et al. 2017; Boehler et al. 2018; Ubeira Gabellini et al. 2019). Although the physical conditions in the cavities are different compared to gaps in full disks (e.g. cavities are also optically thin in the dust at infrared wavelengths), the thermo-chemical models used for transitional disk studies have already addressed many aspects of dust and gas gap modelling (e.g. determining gas surface density profiles) and we use here a very similar modelling approach.

In this work we focus on the disk around the Herbig Ae/Be star HD 163296, which is of particular interest for studying gas observations as three generations of ALMA CO molecular line and continuum observations exist (e.g. de Gregorio-Monsalvo et al. 2013; Flaherty et al. 2015; Isella et al. 2016) that were all combined to produce the DSHARP (Disk Substructure At High Angular Resolution Project) dataset of HD 163296 for the dust at 1.25 mm and the ¹²CO J = 2−1 line (Isella et al. 2018). Furthermore, HD 163296 is likely to be the only disk where indirect evidence for forming planets was reported using different observational methods. Besides the indications for planets from dust gaps in millimetre observations (Isella et al. 2016, 2018), scattered light images also show gaps and ring structures (Grady et al. 2000; Muro-Arena et al. 2018). The strongest constraints for ongoing planet formation in HD 163296 come from the detection of kinematic signatures in CO spectral line observations, which are most likely caused by planets. Such signatures were found for the two dust gaps at r ≈ 86 and r ~ 145 au (Teague et al. 2018, 2019; Pinte et al. 2020) but also outside the observed millimetre continuum disk at r ~ 250 au (Pinte et al. 2018a). However, it is not yet clear if all of the dust gaps in HD 163296 host a planet. The masses of the potential planets derived via different methods can also vary by an order of magnitude (e.g. Pinte et al. 2020).

In this work we use HD 163296 as an example to study the impact of dust gaps and possible gas gaps on high spatial resolution spectral line observations. We focus on the DSHARP dataset as it provides the highest spatial resolution for the dust and gas emission. We model this data with a self-consistent gas and dust model for HD 163296 using the radiation thermo-chemical disk code PRODIMO (PROtoplanetary DIsk MOdel). The main purpose of this paper is not to present a detailed best fit model for HD 163296, but to study the impact of thermo-chemical processes on observables such as radial intensity profiles and channel maps, with the aim of quantifying the importance of thermo-chemical processes for accurate measurements of the gas surface density, kinematic signatures, and pressure gradients.

We first describe our modelling approach and our fiducial models for the disk of HD 163296 (Sect. 2). Using these models we investigate the origin of radial features in the observed CO line emission (Sect. 3) and discuss the impact of thermo-chemical processes on measurements of the rotation velocity and pressure gradients (Sect. 4). In Sect. 5 we present synthetic channel maps, discuss the presence of gas gaps, compare our results to previous modelling, and discuss possible model improvements. Our conclusions are presented in Sect. 6.

2 Methods

2.1 Radiation thermo-chemical modelling

To model the disk of HD 163296, we used the radiation thermo-chemical disk code PRODIMO¹ (PROtoplanetary DIsk MOdel; Woitke et al. 2009, 2016; Kamp et al. 2010). PRODIMO consistently solves for the dust radiative transfer, gas thermal balance, and chemistry for a given static two-dimensional dust and gas density structure. Furthermore, PRODIMO provides modules to produce synthetic observables such as spectral lines (Woitke et al. 2011), spectral energy distributions (SED; Thi et al. 2011), and images. For the chemistry we use a chemical network with 235 gas and ice phase species and 2844 chemical reactions including freeze-out of atoms, molecules and PAHs (polycyclic aromatic hydrocarbons), photo-, thermal-, and cosmic-ray desorption of ices, X-ray chemistry, H₂ formation on grains, and excited H₂ chemistry. This chemical network is described in detail in Kamp et al. (2017). The chemistry is solved consistently with the heating and cooling balance for the gas temperature. Important heating and cooling processes are photo-electric heating, PAH-heating, chemical heating by exothermic reactions, X-ray heating, viscous heating, thermal-accommodation on dust grains, and line cooling by atomic and molecular species. Further details on the implementation of the heating and cooling balance and the various processes can be found in Woitke et al. (2009, 2011) and Aresu et al. (2011).

2.2 HD 163296 model

For this paper we focus on the DSHARP ¹²CO J = 2−1 spectral line and the 1.25 mm continuum observations. As a starting point we used an existing model from the DIANA² project (DIsc ANAlysis) that matches the SED and line fluxes from near-infrared to millimetre wavelengths within a factor of about two (see Woitke et al. 2019; Dionatos et al. 2019 for details). The model includes a separate zone for the inner disk (r ≲ 2.5 au) producing a jump in the gas and dust surface density. Otherwise this model has a smooth density distribution in the gas and dust. The HD 163296 DIANA model was also used as a starting point by Muro-Arena et al. (2018) to model VLT/SPHERE scattered light images and the pre-DSHARP continuum images (Isella et al. 2016) that already show clear signatures of dust gaps and rings in both the thermal and scattered light dust emission.

We updated the DIANA HD 163296 model considering the new Gaia distance (d = 101 pc), adapted stellar parameters, and the disk mass of ~0.2 M_⊙ as derived by Booth et al. (2019) from ¹³C¹⁷O observations. This disk mass is about a factor of three lower than in the original DIANA model (Woitke et al. 2019) and a factor of about three higher than in the older HD 163296 PRODIMO model of Tilling et al. (2012).

The model presented here is in reasonable agreement with the spatially resolved DSHARP data (neglecting the gaps) and still matches the SED and line fluxes to a similar quality as the original DIANA model. All relevant parameters of this initial model are listed in Table A.1 and further details are provided in Appendix A. We use this model as a starting point for our models with dust and gas gaps and as a reference for a disk with a smooth density distribution. How we include gaps and rings in the model is described in Sect. 2.4.

2.3 Observational dataset and synthetic observables

For this work we focus on the HD 163296 DSHARP³ dataset (Andrews et al. 2018) that is described in detail in Isella et al. (2018). The beam size of the dust continuum observations at 1.25 mm is 0.″ 038 × 0.″048 (~ 4.3 au) and the rms noise is 0.023 mJy beam⁻¹. The ¹²CO J = 2−1 line observations have a beam size of 0.″ 104 × 0.″095 (~ 10 au), a channel width of 0.32 km s⁻¹ (but the spectral resolution is 0.64 km s⁻¹ due to Hanning smoothing), and a rms noise per channel of 0.84 mJy beam⁻¹ (Isella et al. 2018).

From this dataset we produced azimuthally averaged radial intensity profiles for the dust and the gas emission following the approach of Huang et al. (2018). We used elliptical apertures, with the minor axis given by the disk inclination i = 46.7° and the orientation given by the position angle of the disk PA = 133.3°. The width of the aperture is about one third of the beam width. For the error bars we used three times the root mean square of the data divided by the square root of independent beams within the elliptical aperture (Isella et al. 2016). This procedure gives a nearly identical dust radial profile as derived by Huang et al. (2018), with the difference that we neglect the azimuthally asymmetric structure observed in one of the gaps (i.e. we did not remove it as was done in Huang et al. 2018). For the gas we apply the same procedure on the integrated intensity (moment 0) image of the ¹²CO J = 2−1 line. The resulting radial intensity profiles will be discussed in Sect. 3.

To produce synthetic observables such as continuum images and spectral line cubes, we use the line radiative transfer module of PRODIMO (Woitke et al. 2011). The molecular data for the line excitation calculation for the CO molecule are from the LAMDA database (Leiden Atomic and Molecular Database; Jankowski & Szalewicz 2005; Yang et al. 2010; Schöier et al. 2005). For C¹⁸ O we assume a fixed ¹⁶O∕¹⁸O isotopologueratio of 498.7 (Scott et al. 2006). The modelled synthetic images are convolved with a 2D kernel using the beam properties of the observations. For the spectral lines the synthetic channel maps are additionally re-gridded to the velocity spacing of the observations. This was done with the CASA software package (Common Astronomy Software Applications, McMullin et al. 2007) using the image functions convolve2d and regrid, respectively. For a comparison to the observations we produce synthetic radial intensity profiles in the same manner as was done for the observations.

Fig. 1 Gas (Σg, solid lines) and dust (Σd, dashed lines) surface densities for the inner 180 au of our fiducial disk models. Top panel: initial SMOOTH model without any gaps or rings. The grey star symbol indicates the gas surface density of 44.4 g cm−2 as measured by Booth et al. (2019) using 13C17O J = 3−2 observationswith a spatial resolution of ≈70 au. Bottom panel: models using the same fitted dust surface density profile but varying gas surface density profiles (see Sect. 2.4). All models are described in Table 1.

2.4 Fiducial models: surface densities

For the presentation and discussion of our results we focus on four fiducial models called SMOOTH, DUSTGAPS, GASGAPS, and GASGAPS S. The only differences between these models are the input gas (Σ_g) and dust (Σ_d) radial surface density profiles as shown in Fig. 1. The model configurations are summarized in Table 1. The SMOOTH model is the starting model with a smooth gas and dust radial surface density profile.

To introduce the dust gaps in our model we follow the approach used by Pinte et al. (2016) for HL Tau and Muro-Arena et al. (2018) for HD 163296. We fit the radial dust intensity profile by adapting the dust surface density via an iterative procedure but keep all other dust properties (i.e. grain size distribution) fixed. To measure the depth of the resulting dust gaps we use the approach of Huang et al. (2018); the gap depth is then given by $${\Delta }_{\text{gap}}=\frac{{\Sigma }_{\text{d}}(r_{\text{ring}})}{{\Sigma }_{\text{d}}(r_{\text{gap}})}.$$(1)

We find gap depths of roughly Δ_gap ~ 1700, 50 and 3 for the dust gaps at r = 48 (dust gap one, DG1), 85 (DG2), and 150 au (DG3), respectively. These numbers are roughly consistent with the results from Isella et al. (2016). They found depletion factors of >100, 70, and 6 for their narrow dust gap model, using rectangular gaps to fit the pre-DSHARP data (six times lower spatial resolution). We note that Isella et al. (2016) measured the gap depth relative to a smooth dust density profile. The presented dust disk model provides us with an accurate model for the dust emission, which is crucial for the interpretation of the continuum subtracted line data. In this DUSTGAPS model only Σ_d is adapted whereas Σ_g is the same as in the SMOOTH model.

To produce gas gaps we apply a parameterized approach using parameters such as gap location, gap width, and depletion factor. The shape of the gaps are a combination of Gaussian wings and a flat bottom (for details see Oberg et al. 2020). The gap locations and widths are chosen to be similar to the dust gaps. The depths of the gaps are given relative to Σ_g from the SMOOTH model. For the GASGAPS model we use gas depletion factors of 1000, 20, and 5 for DG1 to DG3 (i.e. similar to the dust). To study the impact of shallow gas gaps we introduce the model GASGAPS S where we chose a constant gas depletion factor of three for all three gas gaps. The resulting Σ_g profile for the gaps is comparable to the hydrodynamic modelling results of Teague et al. (2018). However, we do not use their profile because we want to be consistent with our gas gap modelling approach. We note that the Σ_g profiles around DG3 are very similar for both models with gas gaps (see Fig. 1). We emphasize that we did not fit the gas data to produce the gas surface density profiles, they are merely used to show the impact of gas gaps on observables.

The resulting gas and dust surface density profiles are shown in Fig. 1. For these models we self-consistently calculate the two-dimensional dust and gas temperature structure, the chemical abundances, and the synthetic observables. In Appendix A we discuss further details of the models and show the resulting two-dimensional density structure, temperature structure, and far-UV radiation field for each model. The CO chemistry within the gaps is discussed in Sect. 4 where the CO disk structure is also shown (Fig. 5).

3 Origin of the radial features in the CO line emission

In Fig. 2 we show the observed azimuthally averaged radial intensity profiles for the dust and the ¹²CO J = 2−1 line in comparison to the SMOOTH model. Even though the SMOOTH model has no radial features, by construction the modelled ¹²CO J = 2−1 radial profile is in good agreement with the observations. For r > 30 au the maximum deviation is about a factor of two but for most radii it is significantly smaller. The largest discrepancies (up to a factor of 3.2) are actually in the inner 30 au where the optically thick dust disk makes the interpretation of the data challenging (e.g. Weaver et al. 2018). However, for our study we focus on the regions around the three prominent dust gaps indicated as DG1, DG2, and DG3 in Fig. 2.

Similarly to earlier lower spatial resolution data (Isella et al. 2016; van der Marel et al. 2018), radial features are visible in the azimuthally averaged ¹²CO J = 2−1 radial intensity profile. However, the DSHARP data now clearly shows that the three weak peaks in the ¹²CO J = 2−1 radial intensity profile coincide with the location of the three prominent dust gaps. In this section we investigate the physical origin of these radial features.

Fig. 2 Azimuthally averaged radial intensity profiles for the 12CO J = 2−1 line emission (top panel) and the 1.25 mm dust continuum emission (bottom panel). The grey solid lines with error bars show the observations. The blue solid lines show the SMOOTH model. Bottom left of each panel: we show the full-width at half-maximum (FWHM) of the beam for the gas and continuum observations, respectively. The vertical dotted grey lines indicate the location of the three prominent dust gaps.

3.1 Impact of temperature

The radial gas and dust intensity profiles for the DUSTGAPS model are shown in Fig. 3. Our dust disk model nicely fits the observed radial dust intensity profile, but also improves the match for the ¹²CO J = 2−1 radial profile compared to the SMOOTH model. The modelled profile shows similar peaks at the location of the dust gaps as the observations. We emphasize that the gas density structure in the DUSTGAPS model is identical to the SMOOTH model.

The most likely explanation for the peaks in the line emission of the DUSTGAPS model is actually a temperature increase within the dust gaps (see Appendix A). This was already postulated by the van der Marel et al. (2018) modelling of pre-DSHARP ALMA data of HD 163296. They argue that even if the gas is depleted within dust gaps, the increase in temperature will lead to an increase in the line emission at the location of dust gaps.

To test this temperature scenario, we constructed a model called DUSTGAPS FREEZE, which includes thedust gaps but whose temperature and chemical structure is identical to the SMOOTH model. As seen in Fig. 3, the DUSTGAPS FREEZE model still shows peaks in the ¹²CO J = 2−1 radial profile at the location of the dust gaps, although there is no temperature increase within the dust gaps. For DG2 and DG3 the radial line intensity profiles of the DUSTGAPS and DUSTGAPS FREEZE model are nearly identical, whereas for DG1 the peak is significantly weaker in the DUSTGAPS FREEZE model. This means that for DG1 thetemperature increase within the gap partly explains the increase in the line emission. However, for DG2 and DG3 this is not the case and temperature and chemical changes are not the cause of the peaks in the ¹²CO J = 2−1 radial profile.

From the bottom panel of Fig. 3 we see that the dust radial profile of the DUSTGAPS and DUSTGAPS FREEZE model are not identical. This is expected because in the DUSTGAPS FREEZE model the dust temperature structure is not consistent with the dust density structure. This implies that a modelling approach with a prescribed temperature structure will lead to inaccurate results, and self-consistent dust radiative transfer modelling is required to derive accurate dust surface density profiles. However, the changes in the dust radial profile in the DUSTGAP FREEZE model are too weak to have a significant impact on the peaks in the ¹²CO J = 2−1 radial profile. Thisis especially true for the outermost gap.

Isella et al. (2018) have shown that in the channel maps of the ¹²CO J = 2−1 line, clear signatures of dust extinction are visible (see also Fig. B.1). As the dust rings are optically thick (or close to optically thick, Dullemond et al. 2018), they absorb the line emission from the back side of the disk, as this emission has to pass through the midplane of the disk and the dust rings. In the next section we explore the impact of the dust absorption process on the radial intensity profiles of spectral lines.

Fig. 3 Same as Fig. 2 but for the DUSTGAPS model and the DUSTGAPS FREEZE model. For the latter the gas and dust temperature and chemical structure is the same as in the SMOOTH model. Top panel: the SMOOTH model (blue solid line) is also shown for reference.

3.2 Contribution of the back side of the disk

The observational data clearly shows the emission from the front and back side of the disk, as the disk is inclined and the ¹²CO J = 2−1 emission layer is high up in the disk. Furthermore, the line emission from the back side of the disk shows radial features that are likely caused by the dust rings (see Isella et al. 2018 and Fig. B.1). We note that those features in the back side emission cannot be identified in the pre-DSHARP data used in Isella et al. (2016) and van der Marel et al. (2018) due to lower spatial resolution and lower sensitivity. The high quality and spatial resolution of the DSHARP line data allows us to separate the emission from the two sides of the disk. By applying Keplerian masking, wecan separate the back and front side emission and produce two sets of azimuthally averaged radial ¹²CO J = 2−1 intensity profiles, one for the front and one for the back side of the disk (see Appendix B for details).

The observed ¹²CO J = 2−1 radial profiles for the back and front side of the disk are shown in Fig. 4. For DG2 and DG3 the radial features in the line emission are clearly seen in the profile for the back side of the disk, but almost vanish in the profile for the front side. For DG1 the situation is much less clear. The reason is most likely the non-perfect Keplerian masking, which becomes more difficult in the inner disk as the front and back side layers are not as well resolved due to the lower height of the emission layers.

To verify this scenario we took our DUSTGAPS model and separated the front and back half of the disk emission during the line radiative transfer step to achieve an exact separation of the two disk halves. For the resulting two line cubes we then applied the same procedure as for the full line cube to produce the radial line intensity profiles. As seen in Fig. 4, the modelled profiles show the same behaviour as the observational data. However, for the back side of the emission the radial features are much more pronounced than in the data.

Both the observational data and the model indicate that the peaks in the full ¹²CO J = 2−1 radial profile are actually a consequence of the presence of dust gaps and rings. As the dust gaps are optically thin, the line emission from the back side of the disk is unaffected, whereas the rings absorb a significant fraction of the emission. This produces the radial features in the back side emission, which are strong enough to show up in the full line radial profile. However, in the model the absorption of the line emission from the back side of the disk is too strong, as is clearly seen in Fig. 4. This indicates that our dust model overestimates the extinction properties of the dust and the optical depth in the rings, despite the fact that the model nicely fits the observed thermal dust emission.

For our model the unsettled dust opacity properties at λ = 1.25 mm for the chosen dust size distribution are ${\kappa }_{\text{d}}^{\text{abs}}=1.24\text{g}{\text{cm}}^{\text{−2}}$ (absorption) and ${\kappa }_{\text{d}}^{\text{sca}}=7.56\text{g}{\text{cm}}^{\text{−2}}$ (scattering). Using the parameter ${\varepsilon }_{\text{d}}:={\kappa }_{\text{d}}^{\text{abs}}/({\kappa }_{\text{d}}^{\text{abs}}+{\kappa }_{\text{d}}^{\text{sca}})$ like Isella et al. (2018), we get ɛ_d = 0.14 for our dust model. Isella et al. (2018) derived values in the range of ɛ_d ~ 0.4−0.6 for the different dust rings, using the continuum and line data combined with dust temperature estimates from various models. This indicates that our model has a too high scattering opacity at millimetre wavelengths and therefore overestimates the total dust optical depth ${\tau }_{\text{d}}^{\text{ext}}$ within the rings. However, our dust model fits the data because scattering can reduce the observed dust emission as has been shown by Zhu et al. (2019), Liu (2019) and Ueda et al. (2020). Consequently, the model overestimates the absorption of the back side ¹²CO J = 2−1 line emission, resulting in too weak line emission at the location of the dust rings (see Fig. 4).

So far the impact of the dust disk on the back side of the gas emission has only been directly observed for HD 163296 (Isella et al. 2018). However, it is likely that such an effect occurs in all disks with prominent dust gaps and rings as long as the dust rings are not optically thin. Furthermore, we can only directly see these effects in inclined disks, as for disks with low inclination the back side of the disk is hidden behind the front side of the disk. Even for the case of an optically thin line, where we observe the sum of the emission from the front and back half of the disk, the dust absorption effect will still produce features in radial intensity profiles and in channel maps. This introduces some degeneracies for the interpretation of radial gas intensity profiles as radial features can be produced by varying temperatures or the presence of optically thick dust rings, and in most cases a combination of those two.

We show that a dust model that fits the millimetre dust image is not necessarily good enough to interpret line observations. To fix our model, it is most likely that we require a more sophisticated dust model that allows for radial variations of the dust properties (i.e. dust size distribution). This scenario is supported by the fact that our model underestimates the line emission from the back side of the disk at the location of the dust rings by factors of five to ten, but is in much better agreement (within a factor of two) at the location of the dust gaps (see Fig. 4).

As already mentioned by Isella et al. (2018), gas lines can provide crucial constraints for the dust properties and a model as presented here is ideal to study this in more detail. However, this is beyond the scope of this paper, but will be explored in future work where we also include an improved dust model.

Fig. 4 Radial intensity profiles for the 12CO J = 2−1 line shown separately for the front and back side of the disk. The solid and dashed black lines show the profiles for the front and back side of the disk for the observations. The solid and dashed red lines are the front and back side profiles derived from the DUSTGAPS model. The grey solid and dashed lines with error bars show the observed full (i.e. both sides) radialprofiles for the gas and dust, respectively, for reference. For easier comparison we multiplied the dust radial profile by a factor of 30 and the profiles for the back side of the disk by a factor of two.

Fig. 5 CO number density and the 12CO J = 2−1 emitting layer (hatched area) for the four fiducial models. The upper and lower bound of the emitting layer at each radius correspond to 15 and 85 percent of the total flux, integrated from top to bottom. The white solid contour shows Tdust = 25 K. To determine the emission layer we assume that we are looking face-on onto the disk.

4 Impact of thermo-chemical processes on the measured rotation velocity

Teague et al. (2018) use ALMA CO isotopologue observations to map the velocity field in the disk of HD 163296. They found deviations from the expected rotation velocity at the position of the observed dust gaps, and concluded that such deviations are likely caused by planets. Planets open up gaps in the gas and therefore affect the radial pressure gradients and consequently the measured rotation velocity.

If radial pressure gradients and the vertical extent of the disk are considered, the rotation velocity of a Keplerian disk is given by (e.g. Weidenschilling 1977; Rosenfeld et al. 2013) $$\frac{v_{\text{rot}}^2}{r}=\frac{G{M}_{*}}{{(r^2+z^2)}^{3/2}}+\frac{1}{\rho }\frac{\partial P}{\partial r},$$(2)

where v_rot is the rotation velocity, G the gravitational constant, M_* the mass of the star, r, z are the radial and vertical spatial coordinates, ρ the gas density, and P the local gas pressure. We neglect here the impact of self-gravity since Rosenfeld et al. (2013) found that this term is of less importance. The second term of Eq. (2) was used by Teague et al. (2018) to infer radial density gradients and profiles. However, they neglected the impact of temperature on the pressure gradient as their observations, in particular changes in the line width across the gaps, suggest that the temperature is not the dominating factor.

With the model presented here, for the first time we can study the impact of the temperature on v_rot in a self-consistent manner. To do this we use a similar presentation as Teague et al. (2018) but measure v_rot directly in the model instead of the observables. In this way we avoid problems arising from observational biases and limitations (e.g. spatial resolution).

The line emitting region of the ¹²CO J = 2−1 line in the models is shown in Fig. 5. At each radial grid point we integrate the line emission vertically from the top towards the disk midplane and define the top and bottom border of the region where the flux reaches 15 and 85% of the total flux, respectively. Within this line emitting layer we calculate several vertically averaged quantities as a function of radius, which are shown in Fig. 6. The quantity ⟨δz⟩ is the change in the height of the line emitting region relative to the SMOOTH model; ⟨d∕g⟩ is the dust to gas mass ratio; ⟨n_⟨H⟩⟩ is the total hydrogen number density ($n_{\langle \text{H}\rangle }=n_{\text{H}}\text{}+\text{}2n_{{\text{H}}_{\text{2}}}$); and ⟨T_d⟩ and ⟨T_g⟩ are the average dust and gas temperatures, respectively. Furthermore, we calculate the average rotation velocity ⟨v_rot⟩ by using Eq. (2) and vertically average v_rot in the same manner as the other quantities. For the presentation in Fig. 6 we use ⟨δv_rot⟩, which we define, following Teague et al. (2018), as the deviation from the expected Keplerian velocity v_Kep (i.e. neglecting the disk height and pressure gradients) $$\langle \delta v_{\text{rot}}\rangle =\frac{\langle v_{\text{rot}}\rangle -v_{\text{Kep}}}{v_{\text{Kep}}}\times 100[\%].$$(3)

We also applied the above-described procedure for the C¹⁸O J = 2−1 line, which is shown in Fig. E.1. However, as the results are qualitatively very similar to the ¹²CO J = 2−1 line, here we only discuss the results for the ¹²CO J = 2−1 line.

Fig. 6 Physical properties in the line emitting region of the 12CO J = 2−1 line forthe four fiducial models (see legend). Top row: gas and dust surface densities Σ (all models with gaps have the same dust surface density profile). The other panels show quantities averaged over the line emitting region (from top to bottom): shift in the height of the line emitting region ⟨δz⟩; dust to gas mass ratio ⟨d∕g⟩; total hydrogen number density ⟨n⟨H⟩⟩; dust temperature ⟨Td⟩; gas temperature ⟨Tg⟩; and the change in the rotational velocity ⟨δvrot⟩ (see Eq. (3)). The vertical dotted lines indicate the locations of the dust gaps DG1, DG2, and DG2.

4.1 Line emitting layer

The height of the line emitting layer changes by up to ⟨δz⟩≈−2 au (see Figs. 5 and 6), when dust gaps are present (i.e SMOOTH versus DUSTGAPS model). The reason is the reduced optical depth within the dust gaps, which increases the temperature but also causes additional photo-dissociation of CO in the upper layers. The impact of the latter effect is however limited as CO is efficiently self-shielding within the gaps (see also Facchini et al. 2018). The change in temperature also affects the population levels of the CO molecule. In regions where the temperature increases, higher levels are populated, the ¹²CO J = 2−1 line intensity decreases in those layers, and the line emitting layer of ¹²CO J = 2−1 moves deeper into the disks with cooler temperatures and higher densities.

In the GASGAPS model the line emitting layer moves even deeper into the disk as now the gas density also drops within the dust gaps. The emitting layer within DG1 almost reaches the disk midplane. For this gap, Σ_g was reduced by a factor of 1000 with respect to the SMOOTH model. In the model with shallow gas gaps (GASGAPS S) the line emitting layer also moves deeper into the disk at the location of the dust gaps with respect to the DUSTGAPS model.

4.2 Rotation velocity

For ⟨δv_rot⟩ (bottom panel of Fig. 6) we first discuss some general properties and then focus on the behaviour of ⟨δv_rot⟩ around DG2 and DG3. Finally we discuss DG1, which shows significantly different physical properties from the other two gaps.

As seen from Eq. (2), v_rot is sensitive to changes in the height of the emission layer and changes in the pressure gradient, which is affected by temperature and density. For the SMOOTH model ⟨δv_rot⟩ < 0 at all radii because the emission layer is high above the midplane (Fig. 5) and the global radial surface density gradient used in our model is negative. The contribution of the latter to ⟨δv_rot⟩ is about three percentage points at r = 180 au, and becomes negligible for r < 100 au (see also Pinte et al. 2018b). Self-gravity would enhance v_rot and therefore would slightly shift ⟨δv_rot⟩ towards zero. However, as already noted we neglect the self-gravity term here because it is of less importance (see Rosenfeld et al. 2013)and would not significantly affect our conclusions.

Around the location of DG2 and DG3 we see a distinct pattern in v_rot in the models with gas gaps. This pattern is caused by the radial change of the pressure gradient, and is qualitatively similar to what was derived by Teague et al. (2018) from the pre-DSHARP data of HD 163296. Across a gap (in radial direction) the density first decreases, reaches its minimum, and then increases again, which results in a negative pressure gradient (⟨δv_rot⟩ < 0), followed by a positive one (⟨δv_rot⟩ > 0) (Teague et al. 2018).

In the DUSTGAPS model we see an inverse pattern at DG2 compared to the models with a gas gap. As there is no gas depletion within the dust gap, this pattern is solely caused by the temperature gradient within DG2. As the gas temperature increases in the gap, due to dust depletion, the sign of the ⟨δv_rot⟩ pattern is switched compared to the case of a gas gap. With respect to the SMOOTH model, ⟨δv_rot⟩ is at most ± 3% in the DUSTGAPS model around DG2.

For DG1 the gas temperature changes by about a factor of two within the dust gap, which is significantly higher than for DG2 and DG3. Due to the strong dust depletion that results in ⟨d∕g⟩ < 10⁻⁴ (depending on the gas depletion), the cooling of the gas via thermal accommodation on grains becomes less efficient. This is even stronger if the gas is also depleted as the cooling rate is proportional to the dust and gas number density (Woitke et al. 2009). Furthermore, photo-electric heating via PAHs⁴ is more efficient inside the gap due to the enhanced far-UV field (see e.g. Fig. A.2). For DG1 in the GASGAPS model viscous heating also kicks in, as the line emitting layer is close to the midplane, resulting in even higher gas temperatures.

As a consequence of the gas temperature increase in DG1, the ⟨δv_rot⟩ pattern is dominated by the temperature gradient in the GASGAPS S model. Only in the GASGAPS model, where the gas is depleted by a factor of 1000 in DG1, do we see the expected ⟨δv_rot⟩ pattern for a gas gap but with weaker amplitudes (i.e. the gas gap appears less deep). As already noted by Teague et al. (2018) and van der Marel et al. (2018), a temperature increase within the gap will influence the derived gas surface density profile and also the ⟨δv_rot⟩ pattern. This will lead to an underestimation of the mass of the forming planet that has opened up the gap. However, our models indicate that in the case of a deep dust gap, such as DG1 in HD 163296, the ⟨δv_rot⟩ pattern can even be completely dominated by the temperature gradient and hide the presence of a gas gap.

We note that it is possible that we over-predict the temperature increase in DG1 in our models. This can be seen in the radial line intensity profiles where the peak at DG1 is too pronounced (see Sect. 3 and also Fig. 8). The measurements of the ¹²CO J = 2−1 brightness temperature by Isella et al. (2018) and Teague et al. (2019) seem to show a peak at the location of DG1 but not as strong as in our model. A reason for a too high temperature in the models could be our dust model (see Sect. 3). We make the dust gaps by removing dust grains of all sizes, including the small ones (i.e. < 1 μm). This leads to less efficient gas cooling via thermal accommodation on grains and to a strong reduction of the optical depth at far-UV wavelengths and consequently to more efficient gas heating. As already noted in Sect. 3, an improved dust model is required (i.e. radial variation of grain size distribution) to improve our fit to the data and to derive a more accurate temperature structure (see also Appendix C for a comparison to the models of Facchini et al. 2018). Nevertheless, our models indicate that in deep dust gaps such as DG1 the impact of the temperature on ⟨δv_rot⟩ is significant, contrary to the two shallow dust gaps.

In Fig. 7 we present again ⟨δv_rot⟩ for the GASGAPS S model, but now showing the individual components that determine the ⟨δv_rot⟩ profile. We note that for all three gas gaps in the GASGAPS S model the gas is depleted by a factor of three with respect to the SMOOTH model (see Sect. 2.4). Figure 7 clearly shows that only for DG3, the ⟨δv_rot⟩ profile is fully determined by the density gradient. For DG2 the temperature gradient has some limited impact on ⟨δv_rot⟩, but for DG1 the ⟨δv_rot⟩ profile is mostly shaped by the temperature gradient. This figure demonstrates the complexity of ⟨δv_rot⟩ measurements but also that each dust gap is quite unique with respect to the ⟨δv_rot⟩ profile.

Our analysis in this section shows the complex interplay of dust and gas gaps and thermo-chemical processes that all have an impact on the measured rotational velocity. Nevertheless, our analysis confirms the interpretation of Teague et al. (2018) that the measured deviations δv_rot for DG2 and DG3 in HD 163296 are most likely caused by gas gaps. For DG1 the quality of the data used by Teague et al. (2018) did not allow for a firm conclusion. However, our dust model for HD 163296 shows that the properties of DG1 are quite different to DG2 and DG3, in particular it is significantly deeper (see also Isella et al. 2016). DG1 is also the only dust gap in HD 163296 that was detected in scattered light observations with VLT/SPHERE (Muro-Arena et al. 2018). The extreme dust depletion in DG1 makes the interpretation of gas line observationsaround DG1 more complex as our modelling shows. To fully understand the nature and origin of DG1, molecular line observations with a similar spatial resolution as the DSHARP continuum observations would be beneficial.

Fig. 7 Deviations ⟨δvrot⟩ from the Keplerian velocity for the GASGAPS S model. Here we show the individual contributions to ⟨δvrot ⟩ from the height of the 12CO J = 2−1 emission layer (grey solid line), the density gradient (brown solid line), and the temperature gradient (pink solid line). The green solid line shows ⟨δvrot⟩ including allthe components (same as in the bottom panel of Fig. 6).

5 Potential gas gaps

Our analysis in Sect. 4 has shown that both dust and gas gaps influence the ¹²CO J = 2−1 line excitation and emitting region and the measured rotation velocity. In this section we discuss to what extent it is possible to infer the presence of gas gaps from the radial intensity profiles and channel maps of our models. Furthermore we also discuss our results with respect to previous gas gap modelling of HD 163296 from the literature.

We note that we did not attempt to precisely fit the ¹²CO J = 2−1 radial intensity profiles or the channel maps by modulating the gas density. First of all the ¹²CO J = 2−1 is highly optically thick and therefore mostly insensitive to the gas density, and secondly as discussed in Sect. 3.2 the limitations of our dust model do not allow for accurate fitting of the gas lines. Nevertheless, our models and synthetic observables are still instructive to devise a strategy for interpreting real observations.

5.1 Radial intensity profiles

In Fig. 8 we show the ¹²CO J = 2−1 and C¹⁸O J = 2−1 radial intensity profiles for our models with gaps. For none of the models is a clear signature of a gas gap, namely a drop in intensity at the position of the dust gaps, visible. This is simply because of the high disk mass M_disk ~ 0.2 M_⊙, which makes even the C¹⁸O J = 2−1 line optically thick in HD 163296 (see also Booth et al. 2019). Although in the GASGAPS model C¹⁸ O J = 2−1 becomes optically thin in the centre of DG1 (measured along the z-axis of the disk), the inclination of the disk and the limited spatial resolution prevents a direct view into the gap. If we view the model face-on we indeed see, at least in the GASGAPS model, a decrease in the line intensity across the dust gaps, although for ¹²CO J = 2−1 we only see thisin the two outermost gaps (see Fig. D.1). This shows the limitation of azimuthally averaged radial intensity profiles for measuring gas gap properties even if the disk inclination is considered by de-projecting the data.

Fig. 8 Azimuthally averaged radial intensity profiles for the 12CO J = 2−1 (coloured solid lines) and C18O J = 2−1 (dashed lines) spectral lines. The models shown are indicated in the legend. The grey solid line with error bars are the 12CO J = 2−1 line observations.

5.2 Channel maps

In Fig. 9 we show five selected velocity channels of the DSHARP ¹²CO J = 2−1 line observationsin comparison to our model results. The modelled line cube was convolved with the same beam size as the observations, however we did not include the noise in the synthetic observables to show the differences between the models more clearly, and also because of the high quality of the data. The goal here is not an exact comparison to the observations but rather to investigate what signatures of dust and/or gas gaps become visible in channel maps.

The channel maps for the SMOOTH model show the perfectly smooth butterfly pattern as expected for an azimuthally symmetric disk with a smooth density and temperature structure. In the DUSTGAPS model we see clearly the impact of the dust rings on the CO emission from the back side of the disk, which is also apparent in the observations (Isella et al. 2018). However, the front side of the disk (the bright emission) also shows clear differences compared to the SMOOTH model. At the location of the dust gaps the emission is brighter, which is caused by the higher temperaturewithin the gaps (see also Sect. 3). It is also apparent that the shape of the butterfly is no longer completely smooth, especially close to the position of the dust gaps. As there are no changes in the gas density compared to the SMOOTH model, those changes in the observed velocity can only be caused by the change of height of the emission layer and the temperature gradient within the dust gaps (see Sect. 4).

In the models with gas gaps the structure in the channel maps becomes more pronounced as expected from changes in the pressure gradient. Compared to the data it seems that the distortions in the channel maps are too pronounced in the GASGAPS model. However, both the DUSTGAPS and GASGAPS S model appear to be roughly consistent with the data. Such distortions are also visible in the observed channel maps, but the spectral resolution is likely to be too low to unambiguously confirm their presence. Compared to the ¹²CO J = 2−1 radial profiles, the channel maps of the models show the differences in the physical model much more clearly (gas gaps or not etc.). Considering the channel maps and the δv_rot profiles (Sect. 4), our model with shallow gas gaps (GASGAPS S) is the one most consistent with the data.

5.3 Comparison to other models

Several models are published in the literature that attempted to model the pre-DSHARP dust and CO isotopologue data published in Isella et al. (2016). To our knowledge the DSHARP dust and ¹²CO J = 2−1 data have not yet been modelled in detail with a dust radiative transfer or thermo-chemical code. The spatial resolution of this pre-DSHARP data is about six times lower for the continuum and two times lower for the ¹²CO J = 2−1 line compared to the DSHARP data.

The most similar approach to our modelling was done by van der Marel et al. (2018) with the radiation thermo-chemical disk code DALI (Dust And LInes, e.g. Bruderer et al. 2012; Bruderer 2013). Similarly to our results, they found that the radial line intensity profiles can actually show peaks at the position of the dust gaps even if gas gaps are present. However, they concluded that this is a temperature effect, whereas our analysis in Sect. 3.2 shows that at least for the gaps DG2 and DG3 the peaks are likely to be a result of the line emission from the back side (see Sect. 3.2). Another difference is that in their models with deep gas gaps (depletion factors of 10 to 100) the line intensity at the position of the dust gaps also drops. They use rectangular gaps for both the gas and the dust (meaning the density depletion factor is constant across the gaps) while we use Gaussian-shaped gap profiles. Furthermore, there are differences in the underlying global disk structure and in the dust model (e.g. treatment of settling, dust opacities), which makes a detailed comparison difficult. However, both studies find the same physical effect of an increased radiation field and temperature in dust gaps. van der Marel et al. (2018) concluded that with the data available at that time, it was not possible to confidently infer the presence of gas gaps in HD 163296 and that the dust rings and gaps could also be caused by icelines, which do not require any gas gaps. In that case one would expect a ⟨δv_rot⟩ pattern similar to our DUSTGAPS model (see Fig. 6), but this pattern is not consistent with the observations of Teague et al. (2018), at least not for DG2 and DG3.

Isella et al. (2016) also modelled the CO line data to derive gas depletion factors within the dust gaps. They used a model with a parameterized two-dimensional temperature structure and a simplified chemistry (e.g. fixed CO abundance and freeze-out at a certain temperature) and varied the CO surface densities for each line individually to fit the observational data. In their model the temperature does not change due to the presence of dust gaps or gas gaps, and also the dust and gas temperature are equal. For their narrow gap model (which is more consistent with the DSHARP data), they found gas depletion factors of 0–2.5, 30–70, and 3–6 for the gaps DG1 to DG3. For DG1 and DG3 this is consistent with our GASGAPS S model, but the deep gas gap at the position of DG2 seems to be inconsistent with our results, in particular from the channel maps. Liu et al. (2018) used the same approach as Isella et al. (2016) for the temperature, chemical, and radiative transfer modelling but derived the gas surface density via hydrodynamic simulations including planets. In their best model the gas gap at the position of DG2 is rather shallow with a gas depletion factor of only a few and for DG1 they found a gas depletion factor of about ten. However, their model is also still in reasonable agreement with the pre-DSHARP gas observations. In summary those models without thermo-chemistry seem to indicate the presence of gas gaps but the derived gas depletion factors can be quite different. The model of Liu et al. (2018) is roughly consistent with our model with shallow gas gaps (DUSTGAPS S, gas depletion factor of three in all gaps), and this is also our model that matches the DSHARP data best.

Fig. 9 Dust continuum images and channel maps for the model series in comparison to the observations. Each row shows, from left to right, the observations, the SMOOTH, the DUSTGAPS, the GASGAPS, and the GASGAPS S model. Toprow: the dust continuum (note the logarithmic scale) is shown for comparison; five remaining rows: 12CO J = 2−1 emission for selected velocity channels. The velocity is given in the top left corner of each row. The dashed ellipses in the channel maps indicate the location of the dust gaps for reference. The filled white ellipses in the bottom left corner of each panel indicate the beam size of the observations.

5.4 How to improve

Deriving ga surface density profiles from observations is of particular importance to constrain hydrodynamic models and the origin of thegaps (e.g. van der Marel et al. 2018). If the gaps are produced by forming planets, it is also possible to derive planet masses froma comparison to hydrodynamic models. As an example we summarize here planet mass estimates that were derived for the DG2 in HD 163296. We have chosen DG2 because for that gap at least five different methods were used in the literature to estimate the planet mass, including direct imaging. Zhang et al. (2018) estimated a planet mass in the range of M_p ~ 0.07−0.6 M_J using the constraints from the DSHARP dust observations; using the pre-DSHARP dust and gas data Isella et al. (2016) and Liu et al. (2018) found M_p ~ 0.3 M_J and M_p ~ 0.46 M_J, respectively; Teague et al. (2018) derived M_p ~ 1 M_J by measuring pressure gradients and Pinte et al. (2020) reported (tentatively) M_p ~ 1−3 M_J using the measured local deviations from Keplerian velocities in the ¹²CO J = 2−1 DSHARP data. Taking the extreme values of the reported mass results in a range of M_p = 0.07−3 M_J, nearly two orders of magnitude. These mass estimates are consistent with planet-mass upper limits derived from high-contrast imaging with VLT/SPHERE (M_p = 3−7 M_J; Mesa et al. 2019) and Keck/NIRC2 (M_p = 4.5−6.5M_J; Guidi et al. 2018).

As noted by Pinte et al. (2020), planet mass measurements using constraints from gas observations tend to give higher planet masses than methods that only use constraints from dust observations. The methods using only gas observations neglect the impact of thermo-chemical processes. However, this most likely leads to an underestimation of the planet mass (see Sect. 4). Hydrodynamic models that use constraints from the dust observations suffer from a limited knowledge of the dust properties (e.g. sizes), the turbulence, or the disk scale-height. For example, the estimated planet mass can increase by about a factor of two for a turbulence alpha of α = 10⁻³ compared to a case with α = 10⁻⁴ (e.g. Zhang et al. 2018). The upper limits derived from direct imaging depend on the system age and the planet and/or sub-stellar evolution model used (e.g Mesa et al. 2019).

For HD 163296 the available datasets are of exceptional quality but nevertheless inferring the presence of planets and in particular estimating their mass is still uncertain. As shown in this work, the interpretation of such data is very complex and the interplay of the dust and gas and thermo-chemical effects has to be considered. However, to improve this situation multi-wavelength data for the gas and the dust are also required. Multiple CO line transitions or observations of other molecules with a similar or better quality as the DSHARP data would provide much stronger constraints on the disk geometry (e.g. scale-height), temperature structure, and chemistry. For the dust, multi-wavelength observations and polarisation measurements can provide constraints on the radial variations of dust properties and also the disk geometry (e.g. Muro-Arena et al. 2018; Ohashi & Kataoka 2019). Only combined modelling of the dust and gas will allow us to unlock the full potential of such high-quality observations.

For the models as presented here, more sophisticated dust modelling considering radial variations in the dust properties is required (for first attempts see e.g. Akimkin et al. 2013; Facchini et al. 2017; Greenwood et al. 2019). Such models will allow us to consider both dust and gas observations to constrain the radial dust properties and, by using various molecular tracers, to eventually derive robust gas surface density profiles. Such information is crucial for hydrodynamic simulations and would provide the required input to improve the implementation of thermo-chemical processes in full hydrodynamic disk simulations (for first attempts see e.g. Vorobyov et al. 2020; Gressel et al. 2020).

6 Conclusions

We have presented a self-consistent thermo-chemical model to interpret the high spatial resolution DSHARP dust and gas data for the disk around the Herbig Ae/Be star HD 163296. We used the radiation thermo-chemical disk code PRODIMO (PROtoplanetary DIsk MOdel) to investigate the importance of the dust component and thermo-chemical processes for the interpretation of high spatial resolution gas observations. We studied the impact of dust and gas gaps on observables such as radial intensity profiles and channel maps, and the role of thermo-chemical processes for measuring pressure gradients in disks. Our main findings are:

• 1.

  Azimuthally averaged radial ¹²CO J = 2−1 intensity profiles derived from the DSHARP data clearly show peaks at the location of the three prominent dust gaps (DG1, DG2, and DG3) of HD 163296 (Sect. 3). The main contribution of those intensity peaks does not come from the increase of temperature within the dust gaps but from the emission of the back side of the disk. The line emission from the back side of the disk is partly absorbed by the dust rings but can pass through the dust gaps causing a bumpy radial ¹²CO J = 2−1 intensityprofile. A possible gas temperature increase within the dust gap can further enhance those peaks. This effect is not unique to HD 163296 but will affect any millimetre-line observations including optically thin lines as long as the dustrings are at least marginally optically thick. A proper dust model is therefore crucial for interpreting high quality millimetre-line observations.

• 2.

  We analysed the impact of thermo-chemical effects on pressure gradients and the measured rotational velocities v_rot (Sect. 4). We find that for deep dust gaps, such as DG1 at r = 48 au in HD 163296, temperature gradients within the gap can have a significant impact on the measured radial v_rot profile. For some cases the temperature gradient can even dominate the profile and hide the presence of any gas density gradient and gas gap. For the two other prominent dust gaps of HD 163296 we find that the impact of temperature gradients is negligible and the measured v_rot profile provides a good estimate of the density gradient and the presence of gas gaps. This is consistent with the findings of Teague et al. (2018) and confirms their interpretation that those gaps are most likely carved by planets. However, the nature and origin of the deep dust gap DG1 (dust depletion factor of 100 to 1700) remains unclear, as the interpretation of the gas data for this gap is more complex. Higher spatial resolution gas data is required to draw a firm conclusion on the origin of that gap (e.g. one or multiple planets, dead zone, ice line).

• 3.

  We showed that azimuthally averaged radial line intensity profiles are difficult to interpret and are not sufficient to derive accurate gas surface density profiles for HD 163296. Both of our gas gap models with shallow gas gaps (gas depletion factor of a few) and deep gas gaps (depletion factor >10) are roughly consistent with the observed ¹²CO J = 2−1 radial intensity profile. However, a comparison of the models with the observed channel maps indicates that deep gas gaps are inconsistent with the data (Sect. 5). Also the model with only dust gaps produces distortions in the channel maps that seem to be consistent with observations. However, the dust gap only model would produce a distinctly different profile for the radial v_rot profile. This shows that considering the velocity information of gas line observations is crucial for the interpretation of dust and gas gaps in disks.

The existing high spatial resolution molecular line observations, which are mostly optically thick lines, and limitations in the models do not allow for the derivation of accurate gas surface density profiles. The large ALMA programme “The Chemistry of Planet Formation”, will provide high spatial resolution data for various molecules, and the modelling will be further improved (e.g. in our case radially varying dust properties). Analysis of such data with models as presented here will provide more stringent constraints for hydrodynamic simulations and gap formation scenarios. It will also improve the mass estimates for forming planets and might tell us if there is a planet in each observed dust gap. Such information will be crucial to relate the forming planet population to the observed exoplanet population (see e.g. Lodato et al. 2019; Ndugu et al. 2019).

Acknowledgements

We want to thank the anonymous referee for a very constructive report that improved the paper. C.R., G.M.-A., and C.G. acknowledge funding from the Netherlands Organisation for Scientific Research (NWO) TOP-1 grant as part of the research programme “Herbig Ae/Be stars, Rosetta stones for understanding the formation of planetary systems”, project number 614.001.552. This research has made use of NASA’s Astrophysics Data System. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013), matplotlib (Hunter 2007) and scipy (Virtanen et al. 2020). We would like to thank the Center for Information Technology of the University of Groningen for their support and for providing access to the Peregrine high performance computing cluster. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2013.1.00366.S, ADS/JAO.ALMA#2013.1.00601.S, and ADS/JAO.ALMA#2016.1.00484.L. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST 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.

Appendix A: Disk model

In Table A.1 we list the stellar and disk parameters for our model of HD 163296. These parameters are fixed for all models presented in the paper as we only vary the underlying gas and dust surface densities. However, we note that the gas to dust ratio changes as a function of radius and height depending on the presence of dust gaps and rings and because of dust settling. The used dust composition is a result of the HD 163296 DIANA model (Woitke et al. 2019). To improve the match to the DSHARP observations, which were not included in the DIANA model, we optimized a few parameters by hand (e.g. a_pow, α_set, or R_tap). However, the new model is still very similar to the DIANA model, and the agreement with the observational data used for the DIANA model is equally good.

The used settling prescription (Dubrulle settling; Dubrulle et al. 1995) depends on the local density (see Woitke et al. 2016 for details) and therefore on the gas and dust depletion within the gaps. As a consequence the dust settling is slightly stronger in the models with gas gaps. However, this effect is not strong enough to affect our dust density profile significantly and the fit to the data is equally good for all models with dust gaps.

For the vertical density structure we stick to the parameterized approach used for the DIANA model. In this approach the scale height of the disk is given by $H(r)=H(100\text{au})\times {(r/100\text{au})}^{\beta }$. We therefore do not self-consistently calculate the hydrostatic equilibrium structure. However, this parameterized approach is a good approximation for the vertical hydrostatic equilibrium and is commonly used for thermo-chemical and radiative transfer modelling ofdisks (e.g. Muro-Arena et al. 2018; van der Marel et al. 2018). Furthermore, the parameterized approach makes the model calculations faster by factors of 10 to 100. Further details on the differences between parameterized and hydrostatic disk models are discussed in Woitke et al. (2016).

In Figs. A.1–A.4 we show the two-dimensional disk density structure, temperature structure, and the radiation field for the SMOOTH, DUSTGAPS, GASGAPS, and GASGAPS S models. Also shown are the input surface densities forthe gas and the dust for reference. A comparison of Figs. A.1 and A.2 clearly shows the impact of the dust gaps on the radiation field and temperature structure. The visual extinction A_V decreases within dust gaps and the first dust gap becomes even optically thin (A_V < 1). As a consequence the dust and gas temperature are always higher within the dust gaps with respect to the regions around the gaps (the rings).

Fig. A.1 Two-dimensional disk structure for the SMOOTH model. From top left to bottom right: surface densities, gas density ρgas, dust density (ρdust), far-UV radiation field (χ), gas temperature (Tgas), and dust temperature (Tdust). The 2D contour plots only show the region from r = 20 to r = 180 au, to clearly show the physical properties within the prominent dust gaps (e.g. Fig. A.2). The contour lines in the plots for the gas and dust density correspond to the values indicated in the colour bar. The red contour line in the radiation field plot (first column bottom row) indicates a visual extinction of unity.

Fig. A.2 Same as Fig. A.1 but for the DUSTGAPS model.

Fig. A.3 Same as Fig. A.1 but for the GASGAPS model.

Fig. A.4 Same as Fig. A.1 but for the GASGAPS S model.

Appendix B: Separating the back and front side

We used Keplerian masking to produce ¹²CO J = 2−1 radial intensity profiles for the front and back side of the disk (see Fig. 4). To construct the Keplerian masks we applied the CLEAN_mask routine from the python package imgcube⁵ developed by R. Teague. As input for the routine we used the same values for the stellar mass, inclination, and position angle as for the models. The height of the emission layer was matched by eye by comparing the generated masks with the observational data. We slightly adapted the imgcube routine CLEAN_mask to produce separate masks for the front and back side of the disk. To produce a mask for the visible part of the back side, we simply removed the areas of the back side mask where it overlaps with the front side mask. In Fig. B.1 we show the two masks for one channel of the ¹²CO J = 2−1 line cube. Using those masks we produced two line cubes one for the front side and one for the back side of the disk. For these line cubes we then applied the same procedure as for the full cube (see Sect. 2.3) to produce radial intensity profiles.

Fig. B.1 One velocity channel of the 12CO J = 2−1 line cube showing the Keplerian masks for the front (white solid line) and back side (red dashed line) of the disk. The colour scale shows the line intensity in Jy beam−1.

Appendix C: Comparison to Facchini et al. (2018)

As noted in the main text, our model does not allow for radially varying dust properties (i.e. dust size distribution). This could have an impact on the physical properties within the gap. For example, if the small dust grains (≲ 1 μm) are not strongly depleted within a dust gap detected at millimetre wavelengths, they could efficiently shield the deeper regions of the gap from far ultraviolet (FUV) radiation. It is out of the scope of this paper to study this in detail but we want to compare our results concerning the physical and chemical conditions within the gaps to the thermo-chemical model of Facchini et al. (2018) that includes radially varying dust properties.

Facchini et al. (2018) modelled a disk around a solar-like star (M_* = 1 M_⊙, L_* = 3 L_⊙) with a gap at r ~ 20 au. The gas density structure was given by hydrodynamic modelling including a planet (with different masses) that produced the gap in the gas density distribution. To determine the radial dust density structure, they used a 1D dust evolution code that considers a radially varying grain size distribution; for the vertical structure they included dust settling. This 2D dust and gas density structure is used as input for the radiation thermo-chemical disk code DALI (e.g. Bruderer et al. 2012; Bruderer 2013) to determine the disk radiation field, gas and dust temperature, and chemical abundances. We want to emphasize that the Facchini et al. (2018) model is not representative for HD 163296 and their modelled gap is also much closer to the star (r = 20 au) than the three prominent gaps for HD 163296, which are all located at r ≳ 50 au.

Similar to our results, they find that due to the presence of the dust gaps the FUV radiation can penetrate deeper into the disk resulting in a stronger FUV radiation field compared to the surrounding disk, and they also find that the dust temperature increases inside the gap. However, for the gas temperature, their model shows a decrease within the dust gap in the deeper layers of the disk. This is different to our models where the gas temperature always increases within dust gaps (see Figs. A.2–A.4). van der Marel et al. (2018) also find that the gas temperature increases within the gaps in their model for HD 163296. They also used DALI, but with a simpler dust model than Facchini et al. (2018). The reason for the difference in the gas temperature might be that for DG1 in our model, where we have the strongest temperature increase, the FUV radiation field within the gap is stronger by a factor of ≳ 10 than in the model of Facchini et al. (2018). This could be caused by the different treatments of the small dust grain population, which can shield FUV radiation efficiently. However, we also note that in our model the stellar UV excess in the range of 6− 13.6 eV is about three orders of magnitude higher than in the model of Facchini et al. (2018), simply because of the different stellar properties used in the models.

According to Facchini et al. (2018) the decrease in gas temperature within the gap has a significant impact on the CO line emission, as most of the CO lines remain optically thick across the gap (similar to our model). Consequently, they find that their modelled gas gap is clearly visible in their synthetic CO line images, contrary to our models and the models of van der Marel et al. (2018). If the gas temperature decreased within the gaps of HD 163296, it would become even more difficult to reproduce the peaks in the observed radial ¹²CO J = 2−1 intensity profile. The only way to produce peaks would be via the absorption of line emission from the back side of the disk by the dust (see Sect. 3.2). However, it is unclear how strong this effect is in the model of Facchini et al. (2018) as they only present synthetic observables for the disk viewed face-on and also do not show channel maps.

This comparison indicates again that the dust model is crucial for determining the gas properties. However, a more systematic approach of thermo-chemical dust and gas gap modelling in disks considering, for example, different disk masses, gaps at different radii, different dust models, and stellar properties is necessary to improve our understanding of high spatial resolution gas line observations and the physical properties of the gas within dust gaps.

Appendix D: Radial profiles for the disk seen face-on

In Fig. D.1 we show the same models as in Fig. 8 but the disk is now seen nearly face-on (i = 5°). If the disk is seen face-on, emission from the back side of the disk will not be visible as both CO lines are optically thick. This is especially apparent for DG3 where the peak in the radial profile vanished in the DUSTGAPS model. For the other dust gaps the peak is still there but is less pronounced, as we now see mostly the impact of the temperature increase within the dust gap. For DG3 the temperature change is not significant.

For the models with gas gaps the emission of the ¹²CO J = 2−1 line decreases in DG2 and DG3 for the face-on case, as we now have a more direct view into the dust gaps. Furthermore, at DG2 and DG3 the gastemperature is lower in the line-emitting region (see Fig. 6) in the models with gas gaps compared to the DUSTGAPS model, hence the weaker emission. For DG1 the temperature increase is much stronger and we always see a peak in the ¹²CO J = 2−1 radial profile at the location of DG1. For the C¹⁸O J = 2−1 line the situation is similar except for DG1 in the GASGAPS model. There the gas depletion within the gap is so strong that C¹⁸ O J = 2−1 becomes optically thin, and the line emission drops significantly.

Fig. D.1 Same as Fig. 8 but the inclination for the models is i = 5°. The dotted vertical lines indicate the location of the three prominent dust gaps. The grey solid line with error bars shows the 12CO J = 2−1 observationsfor reference.

Appendix E: Radial properties for the C¹⁸O J = 2−1 line-emitting layer

In Fig. E.1 we show the same averaged quantities as a function of radius for the C¹⁸ O J = 2−1 line emitting layer as are shown in Fig. 6 for the ¹²CO J = 2−1 line (see Sect. 4 for a detailed description). Although the absolute numbers are different compared to the ¹²CO J = 2−1 line, the general picture is very similar. The ⟨δv_rot⟩ profile shows the same pattern as for the ¹²CO J = 2−1 line only with weaker amplitudes. This is expected as in our models the C¹⁸ O J = 2−1 line is optically thick. Only for DG1 in the GASGAPS model, with a gas depletion factor of 1000, does the C¹⁸ O J = 2−1 line become optically thin.

Fig. E.1 Same as Fig. 6 but for the C18O J = 2−1 line.

References

All Tables

All Figures

Fig. 1 Gas (Σg, solid lines) and dust (Σd, dashed lines) surface densities for the inner 180 au of our fiducial disk models. Top panel: initial SMOOTH model without any gaps or rings. The grey star symbol indicates the gas surface density of 44.4 g cm−2 as measured by Booth et al. (2019) using 13C17O J = 3−2 observationswith a spatial resolution of ≈70 au. Bottom panel: models using the same fitted dust surface density profile but varying gas surface density profiles (see Sect. 2.4). All models are described in Table 1.

In the text

Fig. 2 Azimuthally averaged radial intensity profiles for the 12CO J = 2−1 line emission (top panel) and the 1.25 mm dust continuum emission (bottom panel). The grey solid lines with error bars show the observations. The blue solid lines show the SMOOTH model. Bottom left of each panel: we show the full-width at half-maximum (FWHM) of the beam for the gas and continuum observations, respectively. The vertical dotted grey lines indicate the location of the three prominent dust gaps.

In the text

	Fig. 3 Same as Fig. 2 but for the DUSTGAPS model and the DUSTGAPS FREEZE model. For the latter the gas and dust temperature and chemical structure is the same as in the SMOOTH model. Top panel: the SMOOTH model (blue solid line) is also shown for reference.
In the text

Fig. 4 Radial intensity profiles for the 12CO J = 2−1 line shown separately for the front and back side of the disk. The solid and dashed black lines show the profiles for the front and back side of the disk for the observations. The solid and dashed red lines are the front and back side profiles derived from the DUSTGAPS model. The grey solid and dashed lines with error bars show the observed full (i.e. both sides) radialprofiles for the gas and dust, respectively, for reference. For easier comparison we multiplied the dust radial profile by a factor of 30 and the profiles for the back side of the disk by a factor of two.

In the text

	Fig. 5 CO number density and the 12CO J = 2−1 emitting layer (hatched area) for the four fiducial models. The upper and lower bound of the emitting layer at each radius correspond to 15 and 85 percent of the total flux, integrated from top to bottom. The white solid contour shows Tdust = 25 K. To determine the emission layer we assume that we are looking face-on onto the disk.
In the text

Fig. 6 Physical properties in the line emitting region of the 12CO J = 2−1 line forthe four fiducial models (see legend). Top row: gas and dust surface densities Σ (all models with gaps have the same dust surface density profile). The other panels show quantities averaged over the line emitting region (from top to bottom): shift in the height of the line emitting region ⟨δz⟩; dust to gas mass ratio ⟨d∕g⟩; total hydrogen number density ⟨n⟨H⟩⟩; dust temperature ⟨Td⟩; gas temperature ⟨Tg⟩; and the change in the rotational velocity ⟨δvrot⟩ (see Eq. (3)). The vertical dotted lines indicate the locations of the dust gaps DG1, DG2, and DG2.

In the text

	Fig. 7 Deviations ⟨δvrot⟩ from the Keplerian velocity for the GASGAPS S model. Here we show the individual contributions to ⟨δvrot ⟩ from the height of the 12CO J = 2−1 emission layer (grey solid line), the density gradient (brown solid line), and the temperature gradient (pink solid line). The green solid line shows ⟨δvrot⟩ including allthe components (same as in the bottom panel of Fig. 6).
In the text

	Fig. 8 Azimuthally averaged radial intensity profiles for the 12CO J = 2−1 (coloured solid lines) and C18O J = 2−1 (dashed lines) spectral lines. The models shown are indicated in the legend. The grey solid line with error bars are the 12CO J = 2−1 line observations.
In the text

Fig. 9 Dust continuum images and channel maps for the model series in comparison to the observations. Each row shows, from left to right, the observations, the SMOOTH, the DUSTGAPS, the GASGAPS, and the GASGAPS S model. Toprow: the dust continuum (note the logarithmic scale) is shown for comparison; five remaining rows: 12CO J = 2−1 emission for selected velocity channels. The velocity is given in the top left corner of each row. The dashed ellipses in the channel maps indicate the location of the dust gaps for reference. The filled white ellipses in the bottom left corner of each panel indicate the beam size of the observations.

In the text

Fig. A.1 Two-dimensional disk structure for the SMOOTH model. From top left to bottom right: surface densities, gas density ρgas, dust density (ρdust), far-UV radiation field (χ), gas temperature (Tgas), and dust temperature (Tdust). The 2D contour plots only show the region from r = 20 to r = 180 au, to clearly show the physical properties within the prominent dust gaps (e.g. Fig. A.2). The contour lines in the plots for the gas and dust density correspond to the values indicated in the colour bar. The red contour line in the radiation field plot (first column bottom row) indicates a visual extinction of unity.

In the text

	Fig. A.2 Same as Fig. A.1 but for the DUSTGAPS model.
In the text

	Fig. A.3 Same as Fig. A.1 but for the GASGAPS model.
In the text

	Fig. A.4 Same as Fig. A.1 but for the GASGAPS S model.
In the text

	Fig. B.1 One velocity channel of the 12CO J = 2−1 line cube showing the Keplerian masks for the front (white solid line) and back side (red dashed line) of the disk. The colour scale shows the line intensity in Jy beam−1.
In the text

	Fig. D.1 Same as Fig. 8 but the inclination for the models is i = 5°. The dotted vertical lines indicate the location of the three prominent dust gaps. The grey solid line with error bars shows the 12CO J = 2−1 observationsfor reference.
In the text

	Fig. E.1 Same as Fig. 6 but for the C18O J = 2−1 line.
In the text