© ESO 2020

1 Introduction

One of the primary goals of observations of protoplanetary disks is to improve our understanding of the environment in which planets form. Measurements of the density and temperature distribution in these disks are of fundamental importance. Spatially resolved observations of rotational lines allow measurement of the gas temperature, but there are caveats. For instance, while it is relatively straightforward to measure the gas temperature in the surface layers of the disk, it is much harder to do that in the cold midplane region.

Protoplanetary disks are known to have a vertical temperature gradient: the surface layers are warmer than the midplane regions (e.g., D’Alessio et al. 1998). The temperature in these warm surface layers can be probed with optically thick molecular lines such as the J = 2−1 CO line, which is readily accessible with the Atacama Large Millimeter Array (ALMA). As long as the level populations of the molecule are in local thermodynamic equilibrium (LTE), and the CO line is optically thick, the brightness temperature of the CO emission at the local line center is a direct measurement of the gas temperature in that layer (Weaver et al. 2018). This is an application of the Eddington-Barbier rule for LTE radiative transfer: the observed intensity I_ν is roughly equal to the Planck function B_ν(T) at the depth in the photosphere where τ_ν ≃ 1. For highly optically thick lines, the τ_ν ≃ 1 surface of a protoplanetary disk at the local line center is located high above the midplane. The inferred gas temperature is therefore that of the warm surface layers of the disk, not that of the cold midplane.

However, if one is interested in the conditions under which dust growth and planet formation take place, the relevant temperature is that of the much cooler midplane. The high optical depth of a CO line typically shields the midplane regions from view. To probe these regions one could resort to optically thin molecular lines. However, for optically thin lines, the Eddington-Barbier rule cannot be applied, and the inference of the gas temperature from such optically thin lines becomes much more model dependent and insecure.

The problem of high optical depth of a line that shields the midplane can however be overcome for nonface-on disks, as was shown by Rosenfeld et al. (2013) and as we further demonstrate in this paper.

The solution is to make use of the morphology of the line emission channel maps of a disk in Keplerian rotation when seen at a moderate inclination. Such maps show one-sided “teacup-handle”-shaped emission at the location where the projected Kepler velocity is equal to the velocity-offset of the channel. As shown by radiative transfer models (e.g., Pavlyuchenkov et al. 2007), the vertical temperature structure tends to split this teacup handle into two: one from the upper warm surface layer, and one from the lower warm surface layer. The cold midplane does not show significant emission. This is in part due to the possibility that the CO gas is frozen out onto dust grains in these regions. But even if no freeze-out occurs, the lower temperature of the gas in the midplane regions simply causes the gas to emit less radiation. The split teacup handle thus takes the shape of the double wings of a dragonfly, albeit one-sided.

One of the protoplanetary disks for which this “dragonfly wings” geometry of CO rotational line channel maps is most clearly seen is HD 163296 (Rosenfeld et al. 2013; Pinte et al. 2018a; Isella et al. 2018). The CO 2–1 maps from the DSHARP campaign (Isella et al. 2018) have sufficient signal-to-noise ratio (S/N) that even the weak thermal emission from between the two dragonfly wings is measurable. This is the key to the method described in this paper.

In this paper we demonstrate that if the observations show nonzero brightness between the upper and lower “wings” of the dragonfly, this can be directly translated, using the Eddington-Barbier rule, into the midplane gas temperature. This measurement is direct: it does not require radiative transfer modeling to be interpreted. We apply this to the CO channel maps of the source HD 163296.

2 Testing the method using radiative transfer

2.1 Model setup

To demonstrate the principle, we set up a radiative transfer model with the RADMC-3D code¹ (Dullemond et al. 2012) in order to qualitatively reproduce the CO J = 2–1 channel maps of HD 163296 from the DSHARP campaign (Isella et al. 2018). We do not aim to make an exact replica of the disk. The main aim is to demonstrate the radiative transfer effects on the channel maps, so as to learn how to “read” them.

The parameters of the star are taken from Andrews et al. (2018): M_* = 2.04 M_⊙, L_* = 17 L_⊙, T_* = 9332 K, and the distance is 101 pc. The disk has an inclination with the line of sight of i = 46.7°, and a position angle of 133.3°. We assume the disk is axially symmetric. We take the disk gas surface density to be (1)

with r₀ = 400 au and Σ₀ = 1 g cm⁻². Our standard model has r_d = 150 au and η = 1, leading to a disk gas mass of about M_disk ≃ 0.039 M_⊙. This model is consistent with a typical Lynden-Bell & Pringle viscous disk model, with a power-law behavior close to the star and a weak exponential tail in the outer disk regions. However, in Sect. 4 we also show the results of η = 2, which more closely matches the outer regions, and hints at an outer truncation of the disk.

The full axisymmetric 3D gas density structure is constructed from Σ(r) by demanding hydrostatic equilibrium in cylindrical coordinates (r_cyl, z). The procedure for this is described in Appendix C. This also uniquely determines the azimuthal velocity v_ϕ (r_cyl, z). We choose a radial grid between r_cyl = 10 au and r_cyl = 800 au with 128 logarithmically spaced grid cells. The vertical extent of the z-grid is dependent on r_cyl: from the midplane upward we set up an equally spaced z-grid of 64, with the top of the grid proportional to the cylindrical radius: z_max (r_cyl) = 0.7 r_cyl. We assume mirror symmetry in the equatorial plane.

Once the model is set up in cylindrical coordinates, we map the model onto a spherical coordinate grid (r_spher, θ), because the radiative transfer tool we use, the RADMC-3D package, requires spherical coordinates. We choose the same logarithmic gridding for r_spher as we did for r_cyl. The θ-grid has 64 grid cells linearly spaced in the range π∕2 − 0.7 ≤ θ ≤ π∕2 (where θ = 0 is the north pole, θ = π∕2 is the equator, and θ = π is the south pole).

We choose a simple analytic prescription for the gas temperature, the inspiration for which is taken from the gas temperatures we actually measure with our method (see Sect. 3, Eqs. (4) and (5)). Alternatively, we could have computed the temperature structure self-consistently using Monte Carlo radiative transfer. However, the dust optical depth at large radii is low. This means that for the radii we are interested in here (r ≳ 50 au), the “warm surface, cool midplane” temperature structure is primarily dependent on the conical midplane shadow cast by the denser inner disk regions (r ≪ 50 au). Given that these inner regions are not the target of our study, a prescribed temperature profile is equally sufficient. The prescription is described in Appendix D.

In our model we focus on the CO 2–1 line. In principle we may need to worry about dust extinction, especially close to the star (inwardof about 10 au) or in the known dust rings. In this model we assume however that the regions we probe are sufficiently dust-depleted that the dust extinction does not play a role.

The most uncertain parameter of our model is the ¹²C¹⁶O abundance. We use an abundance of 10⁻⁴ relative to H₂. We do not account for freeze-out or photodissociation of CO. For the case of HD 163296, the star is a Herbig star, the disk midplane temperature is likely to be above 20 K (if only marginally) throughout the disk, meaning that freezeout is probably avoided. For T Tauri stars this may not be the case. Photodissociation is in principle important for the CO abundance in the surface layers. It can be included in a simplified way using a simple column density threshold along the lines of Rosenfeld et al. (2013). But since the total CO column density of the warm surface layers is not our real concern in this paper, we ignore a treatment photodissociation. The CO molecular data are taken from the LAMDA database in Leiden². We focus on the J = 2–1 line. With RADMC-3D we can now compute channel maps of the CO 2–1 line. We assume that the level populations of the CO molecules are in LTE. Microturbulence is not included.

2.2 Resulting channel maps

In Fig. 1 the channel map at δv = 1.76 km s⁻¹ is shown, where δv is defined as the line-of-sight velocity relative to the systemic velocity (5.8 km s⁻¹) of the system. This image shows the familiar Keplerian teacup handle shape: a teardrop-shaped ring of emission from the locations where the projected Kepler velocity is redshifted by 1.76 km s⁻¹. In fact, the teacup handle has a bright front and back layer, separated by a dim region in between. These bright front and back rings, which look like the twin wings of a dragonfly, come from the CO emission from the warm surface layers (Pavlyuchenkov et al. 2007). In between the two dragonfly wings the emission is much weaker because the temperature in the disk midplane region is much lower than in the surface layers.

If the CO emission were optically thin, then both the front-side and the back-side dragonfly wings would be equally bright, and they would both be closed teardrop-shaped rings. We would not be able to distinguish between front side and back side. The dust opacity is negligible in this model at this wavelength and distance from the star. For the real disk around HD 163296 the dust extinction plays a role, but only near the known dust rings at 67 and 100 au (Isella et al. 2018).

However, the CO line emission is highly optically thick (see Fig. 2). The space between the two dragonfly wings is therefore not “hollow” but is highly opaque, and emits optically thick line emission at the midplane temperature.

Therefore, the back-side CO emission is extincted by the cooler CO gas at the midplane. Only a narrow rim of the bright CO emission from the back side is seen, and at some point this backside dragonfly wing disappears entirely behind the optically thick cool CO emission. Only the front-side bright wing is complete.

This geometry is a well-known outcome of radiative transfer models of line emission channel maps of protoplanetary disks (e.g., Pavlyuchenkov et al. 2007; Rosenfeld et al. 2013; Pinte et al. 2018a). To better understand the geometry of what we see in Fig. 1, it is useful to realize that the teacup-handle-shaped object (from lower to upper surface, including the opaque medium between the dragonfly wings) is highly optically thick, while the rest of the disk is transparent at this wavelength. One may therefore regard this teacup handle-shaped object as a truly 3D object. The surface of this 3D object is the τ ~ 1 surface, but this surface is not only at the top and the bottom (the disk surface), but also at the inner and outer edge of the teacup handle. The locations of the inner and outer edge depend on the velocity channel one is looking at, and the intrinsic line width. At any given channel we are looking at the surface of the teacup handle. The brightness is, with the Eddington-Barbier rule, the Planck function at the temperature of the gas at this surface. The dark lanes between the bright dragonfly wings represents the Planck function at the lower temperature of the midplane. By looking at the inner edge and the outer edge of the teacup handle we can thus directly measure the gas temperatureinside the disk, as a function of radius. Different channels allow us to obtain several independent “cuts” through the disk, each yielding the temperature as a function of radius, though in a slightly different range in r.

Fig. 1 CO 2–1 channel map of the model of the disk of HD 163296 at δv = 1.76 km s−1. The color scale stretches from Iν = 0 to 1.2 × 10−12 erg cm−2 s−1 Hz−1 ster−1. The same color scale is used for all subsequent channel map plots.

Fig. 2 Optical depth of the CO 2–1 channel map shown in Fig. 1. The color scale is logarithmic, spanning from τ = 10−2 to τ = 10+5.

Fig. 3 Same as Fig. 1, but now with two paths overplotted that represent equal-radial-velocity curves in the midplane. The inner one is chosen to probe the inner edge of the 3D teacup handle, while the outer one is chosen to probe the outer edge.

2.3 Extracting intensity along iso-velocity curves

Let us test the Eddington-Barbier temperature extraction with our model. To associate the projected Kepler velocity with the location on the channel map, we need to compute the iso-velocity curves corresponding to the channel: the curves for which the projected Kepler velocity along the line-of-sight equals δv.

The azimuthal velocity v_ϕ(r) is the Kepler velocity, that is, . The iso-velocity curve is defined by the equation (2)

which gives a relation between r and ϕ in the plane of the disk. This equation is also valid for the case where v_ϕ(r) deviates from Kepler (see Appendix C). Further, if we know how v_ϕ (r, z) varies with height z above the midplane, we can also use it to compute isovelocity curves for line emission from above and below the midplane. For every pair r, ϕ (and possibly z) we can now project this location onto the sky (see e.g., the appendix of Isella et al. 2018).

Figure 3 shows the same teacup handle as before, but now we overplot the midplane isovelocity curves for two different velocities. The inner curve is for a velocity that is 0.35 km s⁻¹ higher than δv = 1.5 km s⁻¹, while the outer curve is for a velocity that is 0.20 km s⁻¹ lower than δv = 1.5 km s⁻¹. These roughly probe the inner and outer surfaces of the teacup handle, respectively.

In Fig. 4 the intensity profiles along these paths are shown. The intensity is represented in the form of brightness temperature T_b in units of Kelvin. The measured intensity I_ν and the brightness temperature T_b are related via the Planck function according to (3)

The Eddington-Barbier rule simply states that the measured brightness temperature T_b is about equal to the actual gas temperature T_gas at the location along the line of sight where the disk becomes optically thick. If this location is in the disk midplane, then we write T_gas as T_mid. Overplotted in that figure is the actual midplane temperature. The paths start near the star, and then follow the paths in a clockwise direction.

One can see that as expected the curves of T_b and T_mid match very well in the region where the paths probe the inner or outer edge of the teacup handle, which represent the disk midplane. Indeed, the Eddington-Barbier rule works well, and we can therefore use it to probe the midplane temperaturedirectly. In contrast, where the path passes along the bright emission from the front surface layer, we instead get a much brighter intensity because, again according to Eddington-Barbier, we probe the surface layers.

In reality we measure the channel maps as an average over a certain channel width w, that is averaged between δv − w∕2 and δv + w∕2, instead of at the precise velocity δv. We investigate this in Appendix A. It turns out that the effect on the measured midplane temperature is relatively small for channel widths up to a few times the intrinsic line width.

One should also bear beam smearing and noise in mind, both of which limit the feasibility of the method close to the star. Finally, the effect of turbulent broadening may have an effect similar to the channel-averaging: it will widen the surface “wings” a bit more than the thermal broadening does (see Appendix B). In fact, this “wing-widening” is what makes it theoreticallypossible to measure the microturbulence in the disk. For HD 163296 however, only an upper limit could be determined (Flaherty et al. 2015).

Nevertheless, with sufficient resolution and sensitivity, the method we present here constitutes a powerful tool to directly measure the midplane temperature of a protoplanetary disk without having to make any model assumptions.

Fig. 4 Brightness temperature (using the full Planck function) along the two paths shown in Fig. 3. Left: inner path. Right: outer path. Overplotted is the actual midplane temperature of the model. In the regions where the path passes through the part of the image where the darker midplane emission is visible, the match between the synthetic observation and the true gas temperature is nearly perfect, as expected. The deviation between the two curves elsewhere is due to the fact that the path passes through the much brighter surface layer emission.

3 Analysis of the observed CO channel maps of HD 163296

Armed with the knowledge of Sect. 2 we now turn to the CO J = 2–1 channel maps of the disk around HD 163296 of the DSHARP campaign (Andrews et al. 2018; Isella et al. 2018), with the goal of directly measuring the midplane temperature of this disk. The source HD 163296 is one of the best-studied protoplanetary disks (e.g., Grady et al. 2000; Tilling et al. 2012; Rosenfeld et al. 2013; Isella et al. 2016, 2018; Pinte et al. 2018a; Muro-Arena et al. 2018; Liu et al. 2018; Booth et al. 2019; Teague et al. 2019).

Figure 5 shows the channel maps, the top eight of which show the blueshifted channels relative to the systemic velocity, while the bottom eight show the redshifted channels. We note that the extent of the maps (~ 650 au), and therefore of the CO emission, is much larger than the radius of the outermost massive dust ring at 100 au (Isella et al. 2018; Huang et al. 2018).

In using the Eddington-Barbier rule, we assume that the CO lines of this disk are optically thick and in LTE throughout the vertical structure of the disk, all the way out to about 400 au. The models of Sect. 2 show that without any freeze-out, the optical depth of the CO lines at line center can easily be of the order of 1000 or more. Partial freeze-out of CO is therefore possible without violating our assumption of the line being optically thick. For instance, Qi et al. (2015) find that CO freeze-out reduces the gas-phase CO abundance by a factor of about ten, which would in our case easily keep the optical depth of the CO well above unity. But freeze-out is nonetheless an important caveat, and is discussed in more detail in Sect. 5.

We also assume that the dust is optically thin everywhere. This is not correct around 100 and 75 au, the locations of known dust rings, but the method of extraction of the midplane gas temperature described in Sect. 2 requires the teacup handle to be well-resolved, which limits us to regions beyond 100 au anyway.

From the channel maps of Fig. 5 we see that the teacup handle of the channels between δv = −0.8 km s⁻¹ and δv = +0.8 km s⁻¹ are not closed. The CO emission appears to end somewhere between 400 and 500 au, cutting off the outer parts of the teacup handle. We omit these channels from our analysis from here onward, even though our technique can to a certain extent also be applied to these channels.

The remaining channels all show very clear teacup handle shapes. All of these images show the expected optical depth effects discussed in Sect. 2. The bright emission is from the warm disk surface facing us. The emission from the warm disk surface on the other side is mostly extincted, with the exception of a narrow lining, also seen in the radiative transfer models. For the positive velocity channels (and mirrored for the negative velocity channels), this narrow lining is only visible on the southwest side for the inner boundary of the teacup handle and to the northeast side for the outer boundary, exactly as predicted by the radiative transfer models. The disappearance of the lining on the northeast side for the inner boundary and southwest side for the outer boundary of the teacup handle is due to the CO self-absorption. It is therefore clear that there is cooler gas between the two warm surface layers, which is rich enough in CO to be optically thick in the line. The simplest explanation is that the CO at the disk midplane is not frozen out entirely, or not at all. The midplane is thus optically thick in the CO line, meaning that the Eddington-Barbier method for measuring the midplane temperature should work. A slightly more complex scenario could be that the CO is frozen out near the midplane, but not in the layers immediately below the warm surface layers (see e.g., the illustration in Fig. 3 of Pinte et al. 2018b). In that case, our method would not measure the midplane temperature, but only the temperature slightly above the midplane. We discuss this scenario in Sect. 5.

Fig. 5 Observed channel maps of CO J = 2–1 of HD 163296. The channel velocity δv is relative to the system velocity. The horizontal and vertical coordinates are in au.

3.1 Midplane temperature

To probe the midplane temperature we measure the intensity I_ν at locations in the channel maps that we estimate to be the inner- and outer edges of the teacup handle at the midplane. The regions at which we probe I_ν for this purpose are indicated in the channel maps in Fig. 6.

These regions (two per panel) were constructed using isovelocity curves. Given that a protoplanetary disk does not rotate at exactly a Keplerian rotation rate (see Appendix C), we slightly offset the isovelocity curves in velocity so that they nicely probe the inner and outer edges of the teacup handle, as shown in Fig. 6.

Given that each isovelocity curve should disappear behind the CO self-absorption for more than half of its length (because geometrically, we can only see part of the inner and outer edge of the teacup handle if it is optically thick), we limit the curves accordingly. The exact choice of the limited domain along the isovelocity curves is done by eye. Given that the dim emission between the dragonfly wings becomes unresolved when coming too close to the star, we set a minimum radius for each curve, which we also determine by eye. This procedure yields two arcs for each panel along which to determine the brightness temperature: one for the inner edge of the teacup handle and one for the outer edge. For every point along these arcs, the corresponding radius is known through the Keplerian velocity model. Finally, to obtain a S/N we expand each arc in the vertical direction (with respect to the disk geometry) by averaging over 10 arcs with aspect ratios (z∕r) between −0.05 and +0.05. The complete procedure described above defines not just two arcs but two areas in the channel maps, as shown in Fig. 6.

Each arcyields an intensity function I_ν(r), which is converted into a brightness temperature T_b(r) using the Planck function, (Eq. (3)). According to the Eddington-Barbier rule, this brightness temperature is the actual gas temperature at the inner and outer edge of the teacup handle. Each arc is thus an independent measurement of T_mid (r), but along varying segments of r. The results are plotted in Fig. 7.

The midplane gas temperature as a function of radius, as seen in this figure, is almost constant, not departing significantly from 20 K; it is a bit higher at 100 au (about 25 K) and a bit lower at 400 au (about 18 K), but otherwise has a surprisingly shallow slope. The noise is mostly around ± 5 K, though inward of 200 au it increases to about ±10 K. The temperature estimates from the outer edges appear to be a bit less noisy than those from the inner edge. This is presumably because the outer edge is more visible than the inner edge, and there is thus less danger of accidently picking up flux from the bright surface layers.

Using the least-squares method we fit a simple parameterized model to the data: (4)

which is shown in green in Fig. 7.

Another way to analyze these T_mid data is by binning them in radial intervals with the width of the average deprojected beam , where b_maj and b_min are the FWHM beam sizes along the major and minor axes of the beam, respectively, and i = 46.7° is the inclination of the disk (see Appendix H of Dullemond et al. 2018). For each radial bin we decipher which of the tracks shown in Fig. 6 contributes to this bin. Each track, when binned, thus yields an independent measurement of the temperature in that bin. It is not easy to estimate the error on each of these measurements. Instead, we compute the mean μ and the variance σ² of the measurements within the same bin, and estimate the error as Err where N is the number of tracks contributing to that radial bin. The results are shown in Fig. 8 where a “wiggly” curve can be seen.

It is not clear whether these wiggles are real or just an artefact of the analysis method. However, one dip stands out strongly: the dip near 100 au. This is the result of the dust ring. It is a combination of the extinction effect described inIsella et al. (2018) and the continuum-subtraction performed on the CO data. Overplotted with a dotted line is the power law fit obtained by Flaherty et al. (2015) for their high-resolution CO 3–2 data. Their fit is , which they find using a forward radiative transfer modeling approach in combination with an MCMC fitting procedure. Their finding is in good agreement with our model-independent temperature measurement.

Fig. 6 Observed channel maps of CO J = 2–1 of HD 163296. The channel velocity δv is relative to the system velocity. Overplotted at the areas where the midplane temperature is measured. In each panel the top areaprobes the outer edge of the teacup handle, while the bottom area probes the inner edge.

Fig. 7 Midplane temperature measurements within the various areas shown in Fig. 6. The color blue is for the area on the inside of each teacup handle, while the orange color is for the one on the outside. The green thick line is the best-fit parameterized model.

Fig. 8 Midplane temperature measurements. As Fig. 7, but now binned in radial bins.

Fig. 9 As Fig. 5, but now with the predicted locations of the CO surface emission for Keplerian disk (dotted line; corresponding to η = 1) and subkeplerian disk (solid line; corresponding to η = 2).

3.2 Surface temperature

In principle we can apply a very similar procedure to probe the surface temperature. For each channel we just need to know the exact shape of thefront-side surface of the teacup handle, and then sample the intensity along that curve. However, given that this bright ring is rather narrow, a model of its shape must be rather accurate to avoid accidently missing the bright intensity. The shapes of these isovelocity curves on the surface of the disk depend on the flaring/nonflaring geometry of the surface as well as on the possible deviation from Keplerian rotation (see Appendix C). While these parameters can be extracted from the observations, the procedure will be rather sensitive to the proper determination of these parameters.

Earlier studies determined the surface temperature profile in a different and more robust way. Isella et al. (2018) derived the following surface temperature profile from the same dataset: (5)

4 Subkepler rotation due to disk outer edge

Serendipitously, our analysis of the channel maps reveals tentative evidence for subkepler rotation in the very outer regions of the disk around HD 163296. Such subkepler rotation can be produced by a rather sharp exponential cut-off of the outer disk (a disk outer edge), leading to a strong negative pressure gradient d ln p∕d lnr ≪ 0, which then causes the subkepler rotation. This seems to be consistent with the fading away of the CO emission beyond 500 au as seen in the low-velocity channel maps and in zero-moment maps.

The most direct way of spotting the evidence for the subkepler rotation is by overplotting the expected shape of the surface wings of theteacup handle, assuming a disk without an outer edge (i.e., our standard model given by Eq. (1) with η = 1). In that casethe rotation velocity of the gas is very close to Keplerian. This expected shape is seen as the blue dotted lines in Fig. 9. One can see that the blue dotted lines predict the location of the surface emission rather well for the high-velocity channels, but for the lower-velocity channels (δv = ±1.12 km s⁻¹ and δv = ±1.44 km s⁻¹) the predicted surface wings are much larger than what the data would suggest.

If we now modify our model to η = 2, leading to a much stronger decay of the surface density in the outer regions (see Fig. 10), the pressure gradient in the outer disk regions becomes substantial, leading to significant deviation from Keplerian rotation (see Fig. 11). The corresponding predicted emission regions are shown with the blue solid lines in Fig. 9. One can see that now the low-velocity channels are much more closely matched.

To confirm this, we need to compute the synthetic channel maps from the RADMC-3D radiative transfer model. In Fig. 12 one can see the model predictions for the case of η = 2 (all other parameters are kept the same as in the standard model). In our view, the match to the data is reasonably good, even for the low-velocity channels, confirming that the disk with a strong exponential cut-off (η = 2) is consistent with the data.

To verify that the normal disk without an outer edge (η = 1) is not consistent with the data, Fig. 13 compares the predictions of channel δv = +1.12 km s⁻¹ between three models: the original Lynden-Bell & Pringle model with only a very weak cut-off of η = 1 (we call this “no edge”); a model with an outer edge (η = 2) for the CO, but not for the H₂+He gas (thus still close to Keplerian rotation); and a model with an η = 2 outer edge for both the CO and the H₂+He gas (the model shown in Fig. 12), and hence strongly subkeplerian rotation.

This comparison shows that the model without an outer edge (only a weak η = 1 tapering according toLynden-Bell & Pringle; left panel of Fig. 13) does not match the data at all. There is far too much CO emission at large radii (r ≳ 500 au). This is, in itself, not surprising. From the zero-moment maps it was already known before that the CO emission suddenly drops beyond r ≳ 500 au (e.g., Isella et al. 2018). The middle panel of Fig. 13 is meant to verify whether or not an η = 2 exponential cut-off only in the CO abundance but not in the pressure-providing H₂ +He gas could explain the lack of emission at large radii. In this case the pressure gradient remains the same as in the standard model, and thus the disk shows near-Keplerian rotation. As one can see in the figure, the reduced CO abundance at large radii indeed suppresses the emission at large radii, consistent with the data. However, the teacup handle is not closed, which is inconsistent with the data. Only when we apply the η = 2 cut-off to the CO and the H₂+He (right panel of Fig. 13) do we see that the teacup handle is closed, as in the observed data.

The evidence for this subkeplerian rotation in the outer disk regions is admittedly inconclusive. In the outer regions of the disk (r ≳ 200 au), from where this evidence is derived, the S/N becomes suboptimal, meaning that we cannot exclude other effects that could explain the shape of the low-velocity channel maps.

However, if the evidence holds up, then this would indicate that the outer disk has an outer cut-off that is sharper than predicted from viscous disk spreading, that is, the disk has an outer edge. There are many possible reasons for this to be the case: for example, the viscous disk model may not be a good description of the evolution of this protoplanetary disk; a process such as stellar-FUV-driven outside-in photoevaporation by the central star (Gorti et al. 2009) may be truncating the disk; there may be an unseen companion that is massive enough to truncate the disk from the outside, but not massive enough to be detected through direct imaging; or finally it could also be that the turbulent viscosity of the disk is so low that the disk has not yet had time to viscously expand.

Fig. 10 Comparison of the surface density profile of the standard model of this paper (η = 1, blue curve) and that of the disk with an outer edge (η = 2, orange curve). For clarity, the axes are logarithmic.

Fig. 11 Comparison of the midplane azimuthal gas velocity profile of the standard model of this paper (η = 1, blue curve) and that of the disk with an outer edge (η = 2, orange curve). Overplotted with a dotted line is the midplane Kepler velocity. In contrast to Fig. 10, the axes are here linear.

5 Effect of freeze-out of CO

The method of measuring the midplane gas temperature described in Sects. 2 and 3 relies on the ¹²CO J = 2 − 1 line being optically thick even at the cold midplane. If CO freeze-out onto dust grains does not play a role, then this condition is easily met. However, there is evidence that CO freeze-out is taking place in the disk around HD 163296. The most compelling evidence is the detection of ring-shaped N₂ H⁺ J = 3 − 2 emission (Qi et al. 2015). For any disk in which the midplane temperature drops below ~ 20 K in the outer regions, CO freeze-out is likely to occur. This could drastically lower the optical depth of the gas-phase CO in the J = 2 − 1 line in the midplane regions.

In fact, for HD 163296 our measurements of Sect. 3 yield midplane temperatures that are not far from the temperature at which CO molecules tend to freeze out onto dust grains. It is therefore not unlikely that (a) part of the gaseous CO is frozen out (reducing its gas-phase abundance), and consequently that (b) the temperature we measure is in fact the temperature of the “CO snow surface”.

The CO snow surface is an axisymmetric curved surface. At the midplane, where the temperature at a given radius is the lowest, the location of this snow surface (defined at the midplane as the snow line) can be relatively far in, down to perhaps even 75 au (Qi et al. 2015, when correcting for the new Gaia distance). However, the gas is warmer at the surface, and the CO stays in the gas phase much farther out. It might therefore be that, while the surface emission remains unaffected by CO freeze-out, the midplane CO emission and extinction may be strongly suppressed; in the most extreme case it is even absent altogether. In that case, our method of measuring the gas temperature by employing Eddington-Barbier will instead measure the gas temperature at the bottom of the nonfrozen-out surface layer on the far side (see e.g., the illustration in Fig. 3 of Pinte et al. 2018b). Since that is the freeze-out temperature by definition, we would be measuring the freeze-out temperature as a function of radius instead of the midplane temperature. This by itself is not uninteresting, because it allows us to calibrate the CO freeze-out temperature under different conditions. This temperature is dependent on gas pressure (e.g., Davis 2005), and so it is not surprising that we find a shallow, but gradual decrease of the temperature with radius (cf. Fig. 8).

If we wish to measure the midplane temperature, and not the freeze-out temperature, the following question is important: how strongly must the CO freeze-out deplete the gas-phase CO to start affecting the reliability of our method? For our model, the ¹²CO J = 2–1 line is very optically thick (of theorder of τ ~ 10³ or higher, see Fig. 2). A factor five or ten reduction in the gas-phase CO abundance, such as suggested by Qi et al. (2015), would still leave the midplane optically thick at line-center, as required for our method. But what if the depletion is more extreme?

To test this, we adjust our model in the following way. We set the CO freeze-out temperature to T_freezeout = 23 K, so that the midplane freeze-out radius lies at about 95 au. We now mimic the freeze-out by reducing the CO abundance by a factor f_CO everywhere where T <T_freezeout. We experiment with several values of f_CO. In Fig. 14 we show the resulting channel maps for δv = 1.44 km s⁻¹, for three values of f_CO: f_CO = 1 (left), f_CO = 10⁻² (middle), and f_CO = 10⁻⁴ (right).

By comparing the left and middle panels of Fig. 14, one can see that a factor of 100 reduction of gas-phase CO near the midplane still leaves the midplane region optically thick. However, an optically thin layer is formed between the midplane and the warm surface layer. One can see through this layer, and the slightly warmer temperature of the bottom of this layer is revealed. In the right panel, with a factor of 10 000 reduction of gas-phaseCO, the entire region between the warm surface layers is optically thin. The image looks like an empty “oyster shell”. In this case, one is looking at the bottom of the warm surface layer, as suggested by Pinte et al. (2018b).

The effect of the freeze-out is perhaps somewhat subtle in the δv = +1.44 km s⁻¹ channel shown in Fig. 14. It becomes more pronounced in the δv = 1.12 km s⁻¹ channel, as shown in Fig. 15, because lower-velocity channels emphasize the very outer regions. In this figure we also compare the standard disk model (η = 1) to the cut-off disk model (η = 2). The η = 1 models clearly produce more emission at large radii than is seen in the observations, that is, they are inconsistent. Also, the geometry for the frozen-out cases does not seem to match the observations. The η = 2 model matches the lack of emission at very large radii. The freeze-out cases (bottom-middle and bottom-right panels) produce a geometry reminiscent of a “tulip flower”. This is caused by the detachment of the top and bottom surface layers, because the gas in between is no longer optically thick. While a more detailed comparison to the observations (using an ALMA simulator) would be necessary to rule this geometry out, a simple visual comparison to the lower-left panel of Figs. 6 and 9 reveals no clear evidence of this tulip flower geometry.

These results depend not only on the CO depletion factor f_CO, but also on the disk mass. Booth et al. (2019) recently found evidence that the gas mass of HD 163296 is M_disk = 0.21 M_⊙, which is about five times more massive than in the model of this paper. Given that our model assumes LTE, increasing the number density of CO everywhere by that same factor of five would require a five times stronger depletion to obtain the same results for the disk’s optical depth in ¹²CO J = 2 − 1 as in Fig. 14.

Astrochemical models such as those of Aikawa & Nomura (2006) and Walsh et al. (2010) predict gas phase CO depletion factors as strong as f_CO ~ 10⁻⁴. In that case, if we probe beyond the CO snowline radius, our method likely measures the freeze-out temperature, not the midplane temperature (at least for the cold outer regions).

On the other hand, in many protoplanetary disks, including HD 163296, the outer regions of the disk are likely to have a strongly reduced dust content, as a result of radial drift of dust grains (Birnstiel & Andrews 2014). Indeed, for HD 163296 the dust thermal emission is confined to the inner 200 au, with most of the emission coming from within 120 au. With less dust surface area to freeze on to,the timescale for freeze-out increases. Turbulent mixing could subsequently increase the abundance of gas-phase CO in the deep interiors of the disk (Xu et al. 2017). Whether or not this effect is enough to make the midplane optically thick in the ¹²CO J =2–1 line is unclear and requires further modeling.

In conclusion: the possibility that CO freeze-out could render the J = 2 − 1 line optically thin between the two warm surface layers of the disk could jeopardize the efficacy of our method for measuring the midplane gas temperature. For disks with even higher inclination it will be possible to directly check if the gas between the warm surface layers is optically thick or thin: in the case of an optically thin gas layer, we should see the back-side and front-side emission detach from each other, since we will then be able to look through this otherwise optically thick region.

Fig. 12 As in Fig. 6, but now the synthetic channel maps obtained from our RADMC-3D model. Shown here is the case of η = 2 (see Eq. (1)), i.e., the model with an outer edge and, as a result, strong subkeplerian rotation.

Fig. 13 Comparing the model predictions for the δv = +1.12 km s−1 channel. Left: original standard model without outer edge (or more precisely: with only the standard Lynden-Bell and Pringle cutoff of η = 1). Middle: same model as left for the H2+He gas, but a η = 2 edge for the CO gas. Right: both H2+He and CO gas outer edge with η = 2, i.e. with the resulting strong deviations from Keplerian rotation.

Fig. 14 Comparison of model predictions for the δv = +1.44 km s−1 channel for different levels of CO freeze-out. Left: no freeze-out, middle: freeze-out by a factor of hundred, right: freeze-out by a factor of ten-thousand.

Fig. 15 As in Fig. 14, but now for a lower velocity channel (δv = +1.12 km s−1), thus emphasizing the outer disk regions. Compared are the standard model without the outer disk cut-off (η = 1) in the top row, and the model with a steeper disk cut-off (η = 2) and resulting strong subKepler rotation in the bottom row. We note that this is the same velocity channel as used in Fig. 13, and the leftmost and rightmost panels of that figure are the same as the top-left and bottom-left panels here, respectively.

6 Conclusions

Here, we analyze the CO channel maps of HD 163296 in more detail than in our earlier studies. Our findings can be summarized as follows.

• 1.

  By employing the Eddington-Barbier rule to line channel maps of sufficiently inclined protoplanetary disks, it is possible to directly measure the midplane gas temperature of the disk from the channel maps without having to do model fitting. The only model assumption is that the line is very optically thick, especially near the midplane. Strong freeze-out could invalidate this assumption. The method requires the vertical structure of the CO emission to be well spatially resolved, and therefore the method lends itself best to the outer regions of the disk.

• 2.

  For HD 163296 we find the midplane temperature to be almost flat around 20 K. There is a gentle slope, from ~25 K at 100 au down to ~18 K at 400 au.

• 3.

  If the observations do not show emission between the dragon fly wings, our method can provide an upper limit on the midplane temperature, if the S/N and beam smearing issues are properly accounted for.

• 4.

  The fact that the temperature we find for HD 163296 is close to the CO freeze-out temperature could mean that we in fact measure the temperature right at the vertical CO freeze-out surface, as suggested for IM Lupi in Pinte et al. (2018b). If the CO freeze-out factor is only moderate (a factor of five or ten), as in Qi et al. (2015), the midplane would still be optically thick. However, chemical models tend to predict much stronger freeze-out, which would yield an optically thin midplane. Our analysis does not provide conclusive evidence for either scenario.

• 5.

  We find tentative evidence for strong subkepler motion in the very outer regions of the disk around HD 163296. This is consistent with a steep negative pressure gradient there, suggesting that the outer edge seen in the CO emission is not just the removal (photodissociation or freeze-out) of CO, but the outer edge of the gaseous disk itself.

Acknowledgements

C.P.D. thanks Rich Teague for useful discussions, and for a preview to the results of his inference of v_ϕ from the same channel map dataset. C.P.D. acknowledges work presented by Alex Fontana in his Masters thesis that partly inspired the work presented in this paper. C.P.D. and N.D. acknowledge funding by the DFG under grant number DU 414/20-1, which is part of SPP 1992 “Diversity of Extrasolar Planets”. S.A. acknowledges support from the National Aeronautics and Space Administration under Grant No. 17-XRP17_2-002 issued through the Exoplanets Research Program. We thank the anonymous referee for an important and useful suggestion that has helped us improve the paper.

Appendix A Effect of the channel width

The high-resolution CO ALMA observations from the DSHARP campaign have a channel width of 0.32 km s⁻¹. This is not negligible compared to the intrinsic line width. In fact, when one overplots the channel maps of two neighboring channels, one notices that the teacup handles are hardly overlapping, meaning that we are clearly on the border of unresolving the line in frequency space. Instead of the DSHARP data, we could use the earlier data from Isella et al. (2016), since these data have higher velocity resolution, but they have lower spatial resolution and less sensitivity.

In this appendix we use the RADMC-3D model to verify the size of the deviation from the true midplane temperature when applying the Eddington-Barbier rule to a channel with width 0.32 km s⁻¹. To this end we repeat the experiment we did in Sect. 2, but now we make 32 “subchannel” images between δv = 1.76 ± 0.16 km s⁻¹, i.e., in total over a channel width of 0.32 km s⁻¹. We then take the mean of all 32 subchannel images and use that to extract the midplane temperature.

Figure A.1 shows the effect of the band-averaging. The shape of the averaged teacup handle (panel on the right) is slightly less clear than in the original single-frequency maps, but overall the shape remains the same. The resulting Eddington-Barbier midplane temperature extractions are shown in Fig. A.2. While the match between the extraction (blue curve) and the real midplane temperature (orange curve) is less good than for the single-frequency case (Fig. 4), the error is only about 1 to 2 Kelvin. This is well within the noise of the data of HD 163296.

Fig. A.1 CO2–1 channel map as in Fig. 1, but now for a finite-width channel with a width of 0.32 km s−1 (see Appendix A). Left: zero-width channel map at δv = 1.76 − 0.16 km s−1. Middle: zero-width channel map at δv = 1.76 + 0.16 km s−1. These are the two extreme sides of the 0.32 km s−1 band. Left: low-velocity extreme. Middle: high-velocity extreme. In between, another 30 intermediate “subchannel” maps are computed. The mean of all these subchannel maps is shown to the right.

Fig. A.2 As Fig. 4, but now for the band-averaged channel map model of Fig. A.1-right. Left: inner path. Right: outer path. Overplotted is the actual midplane temperature of the model.

Appendix B Effect of unresolved turbulence

In the standard models of this paper we set the unresolved turbulence to zero. The only line broadening is due to the thermal broadening of the CO molecules themselves. Here we demonstrate the effect of possible turbulent broadening on the shape of the teacup handle. We start from the subkepler (η = 2) model of Sect. 4. We then add turbulence at three different magnitudes compared to the local sound speed: v_turb = 0.1 c_s, v_turb = 0.3 c_s, and v_turb = 1 c_s. The results are shown in Fig. B.1.

Fig. B.1 Demonstration of the effect of unresolved turbulence in the disk on the δv = 1.12 km s−1 channel map. Only for very strong turbulence (vturb = cs) do the effects become so large that they are easily seen in the maps. This is why it is so challenging to measure the strength of turbulence from the channel maps, as studied by Flaherty et al. (2015).

Appendix C Two-dimensional disk structure and orbital velocity

Due to the radial gas pressure gradient, the protoplanetary disk does not orbit with exact Kepler velocity. The vertical gas pressure gradient also affects the orbital velocity of the gas, albeit in an indirect way. In cylindrical coordinates it is possible to set up an axially symmetric gas disk structure in which the gravity from the star, the gas pressure gradients in the disk, and the centrifugal force are in exact balance (Rosenfeld et al. 2013). Let us define r_cyl as the cylindrical radial coordinate, while r_spher is the spherical radial coordinate, and z is the vertical coordinate measured from the midplane of the disk. They are related via . The gravitational force per unit mass of gas is (C.1)

The centrifugal force is (C.2)

and the pressure gradient force is (C.3)

where is the gas pressure, with c_s the isothermal sound speed. Force balance is reached when f_g + f_c + f_p = 0. For the e_z component this leads to the equation (C.4)

The right-hand-side does not depend on p or ρ, but does depend on z. If the temperature structure, and thus c_s(r_cyl, z), is a given function of r_cyl and z, then Eq. (C.4) can be directly numerically integrated from the midplane upward, assuming we know the midplane pressure p(r_cyl, 0). This yields the vertical pressure profile p(r_cyl, z). By dividing this by the known , we obtain the density profile ρ(r_cyl, z). Typically we do not know p(r_cyl, 0) a priori, butinstead we prescribe the radial surface density profile (C.5)

Numerically we start the vertical integration of Eq. (C.4) with a trial value of p(r_cyl, 0), integrate Eq. (C.4), and rescale ρ(r_cyl, z) to match Eq. (C.5). This immediately gives us the 2D density structure.

Next we demand that also the radial forces be in balance. This yields the equation (C.6)

Although this looks like a differential equation, it is not, because we already know the 2D structure of p(r_cyl, z), and so we can numerically evaluate the third term of the equation, as well as the second one. This immediately leads to the value of v_ϕ (r_cyl, z). The resulting solution is in perfect 2D force balance. It can happen that Eq. (C.6) is not solvable, i.e., it yields imaginary v_ϕ velocity. In that case the radial pressure gradient force is apparently larger than the gravitational force, in which case no azimuthal velocity can be found that leads to equilibrium. In practice, this part of the disk will likely escape through a thermal wind.

We note that even without substantial pressure support, Eq. (C.6) shows that the balance between centrifugal and gravitational force leads to a reduction of |v_ϕ | with increasing distance from the midplane. The surface layers of a disk therefore tend to rotate slower than the midplane (see e.g., Rosenfeld et al. 2013).

Appendix D Simple temperature prescription

To keep the radiative transfer model of Sect. 2 as simple as possible, we decided to prescribe the 2D temperature structure of the disk rather than compute it. We use the measured midplane and surface temperature radial profiles (Eqs. (4) and (5)) as the basis of this prescription. We use cylindrical coordinates, and interpret the r in Eqs. (4) and (5) as the cylindrical radius r_cyl. We parameterize the z-dependence of the temperature T(r_cyl, z) as (D.1)

where (D.2)

This introduces two parameters: the dimensionless height of the surface layer h_s and the dimensionless width of the transition from midplane to surface temperature w_s. Equation (D.2) creates a smooth two-layer disk model, with T = T_mid below the surface layer and T = T_mid above.

References

All Figures

	Fig. 1 CO 2–1 channel map of the model of the disk of HD 163296 at δv = 1.76 km s−1. The color scale stretches from Iν = 0 to 1.2 × 10−12 erg cm−2 s−1 Hz−1 ster−1. The same color scale is used for all subsequent channel map plots.
In the text

	Fig. 2 Optical depth of the CO 2–1 channel map shown in Fig. 1. The color scale is logarithmic, spanning from τ = 10−2 to τ = 10+5.
In the text

	Fig. 3 Same as Fig. 1, but now with two paths overplotted that represent equal-radial-velocity curves in the midplane. The inner one is chosen to probe the inner edge of the 3D teacup handle, while the outer one is chosen to probe the outer edge.
In the text

Fig. 4 Brightness temperature (using the full Planck function) along the two paths shown in Fig. 3. Left: inner path. Right: outer path. Overplotted is the actual midplane temperature of the model. In the regions where the path passes through the part of the image where the darker midplane emission is visible, the match between the synthetic observation and the true gas temperature is nearly perfect, as expected. The deviation between the two curves elsewhere is due to the fact that the path passes through the much brighter surface layer emission.

In the text

	Fig. 5 Observed channel maps of CO J = 2–1 of HD 163296. The channel velocity δv is relative to the system velocity. The horizontal and vertical coordinates are in au.
In the text

	Fig. 6 Observed channel maps of CO J = 2–1 of HD 163296. The channel velocity δv is relative to the system velocity. Overplotted at the areas where the midplane temperature is measured. In each panel the top areaprobes the outer edge of the teacup handle, while the bottom area probes the inner edge.
In the text

	Fig. 7 Midplane temperature measurements within the various areas shown in Fig. 6. The color blue is for the area on the inside of each teacup handle, while the orange color is for the one on the outside. The green thick line is the best-fit parameterized model.
In the text

	Fig. 8 Midplane temperature measurements. As Fig. 7, but now binned in radial bins.
In the text

	Fig. 9 As Fig. 5, but now with the predicted locations of the CO surface emission for Keplerian disk (dotted line; corresponding to η = 1) and subkeplerian disk (solid line; corresponding to η = 2).
In the text

	Fig. 10 Comparison of the surface density profile of the standard model of this paper (η = 1, blue curve) and that of the disk with an outer edge (η = 2, orange curve). For clarity, the axes are logarithmic.
In the text

	Fig. 11 Comparison of the midplane azimuthal gas velocity profile of the standard model of this paper (η = 1, blue curve) and that of the disk with an outer edge (η = 2, orange curve). Overplotted with a dotted line is the midplane Kepler velocity. In contrast to Fig. 10, the axes are here linear.
In the text

	Fig. 12 As in Fig. 6, but now the synthetic channel maps obtained from our RADMC-3D model. Shown here is the case of η = 2 (see Eq. (1)), i.e., the model with an outer edge and, as a result, strong subkeplerian rotation.
In the text

Fig. 13 Comparing the model predictions for the δv = +1.12 km s−1 channel. Left: original standard model without outer edge (or more precisely: with only the standard Lynden-Bell and Pringle cutoff of η = 1). Middle: same model as left for the H2+He gas, but a η = 2 edge for the CO gas. Right: both H2+He and CO gas outer edge with η = 2, i.e. with the resulting strong deviations from Keplerian rotation.

In the text

	Fig. 14 Comparison of model predictions for the δv = +1.44 km s−1 channel for different levels of CO freeze-out. Left: no freeze-out, middle: freeze-out by a factor of hundred, right: freeze-out by a factor of ten-thousand.
In the text

Fig. 15 As in Fig. 14, but now for a lower velocity channel (δv = +1.12 km s−1), thus emphasizing the outer disk regions. Compared are the standard model without the outer disk cut-off (η = 1) in the top row, and the model with a steeper disk cut-off (η = 2) and resulting strong subKepler rotation in the bottom row. We note that this is the same velocity channel as used in Fig. 13, and the leftmost and rightmost panels of that figure are the same as the top-left and bottom-left panels here, respectively.

In the text

Fig. A.1 CO2–1 channel map as in Fig. 1, but now for a finite-width channel with a width of 0.32 km s−1 (see Appendix A). Left: zero-width channel map at δv = 1.76 − 0.16 km s−1. Middle: zero-width channel map at δv = 1.76 + 0.16 km s−1. These are the two extreme sides of the 0.32 km s−1 band. Left: low-velocity extreme. Middle: high-velocity extreme. In between, another 30 intermediate “subchannel” maps are computed. The mean of all these subchannel maps is shown to the right.

In the text

	Fig. A.2 As Fig. 4, but now for the band-averaged channel map model of Fig. A.1-right. Left: inner path. Right: outer path. Overplotted is the actual midplane temperature of the model.
In the text

	Fig. B.1 Demonstration of the effect of unresolved turbulence in the disk on the δv = 1.12 km s−1 channel map. Only for very strong turbulence (vturb = cs) do the effects become so large that they are easily seen in the maps. This is why it is so challenging to measure the strength of turbulence from the channel maps, as studied by Flaherty et al. (2015).
In the text