The First Detection of ¹³C¹⁷O in a Protoplanetary Disk: A Robust Tracer of Disk Gas Mass

The mass of a disk sets a limit on the material available for forming a planetary system and can influence the mode of giant planet formation. Most disk gas masses rely on observations of CO that are extrapolated to a total gas mass by assuming a constant $\mathrm{CO}$/${{\rm{H}}}_{2}$ abundance ratio in the disk. However, if the line emission is optically thick this will result in an underestimated disk mass (e.g., Bergin & Williams 2017). Spatially resolved line emission from both ${}^{13}{{\rm{C}}}^{16}{\rm{O}}$ and ${}^{12}{{\rm{C}}}^{18}{\rm{O}}$ are now commonly detected in protoplanetary disks and have been used to determine disk gas masses (e.g., Ansdell et al. 2016; Long et al. 2017). The second-most rare CO isotopologue, ${}^{13}{{\rm{C}}}^{18}{\rm{O}}$, has been detected in the TW Hya disk, and this emission has been proposed to be optically thin and traces the disk midplane, whereas the ${{\rm{C}}}^{18}{\rm{O}}$ emission is optically thick within the midplane snow line (Zhang et al. 2017). The underestimation of disk gas mass due to optically thick emission will be more significant in more massive gas-rich disks, i.e., those around Herbig Ae/Be stars versus those around T Tauri stars.

Low CO gas masses, with respect to the dust mass, have been consistently measured in disks, and complementary HD observations imply that this is because of the depletion of gas-phase CO in disks relative to that in the interstellar medium (ISM; Bergin et al. 2013; McClure et al. 2016). CO can be depleted from the gas phase via freezeout onto the icy grains in the cold midplane, and subsequent conversion to ${\mathrm{CO}}_{2}$ and more complex organic species, e.g., ${\mathrm{CH}}_{3}\mathrm{OH}$ (e.g., Bosman et al. 2018). Additionally, photodissociation via far-UV radiation destroys CO in the upper disk atmosphere, and isotope-selective photodissociation can enhance the various isotopologue ratios relative to ${}^{12}{{\rm{C}}}^{16}{\rm{O}}$ in the atmosphere (e.g., Miotello et al. 2014). These chemical effects make the conversion from CO gas mass to total gas mass a nontrivial task.

The protoplanetary disk around HD 163296 has been well characterized with the Atacama Large Millimeter/submillimeter Array (ALMA). Band 6 and 7 observations show rings in both the continuum and the CO gas emission (e.g., Isella et al. 2016; Notsu et al. 2019; see Figure 1(a)). There are four proposed ≈0.5–2 M_J planets in this disk inferred from the dust and gas rings, and deviations from Keplerian motion in the CO gas kinematics (Isella et al. 2016; Liu et al. 2018; Pinte et al. 2018; Teague et al. 2018). Recently, ≈5 au resolution observations of the continuum emission revealed an additional gap and ring in the inner disk as well as an azimuthal asymmetry in one of the previously detected rings (see Figure 1(a); Isella et al. 2018). Hence, the proposed planet-induced structures in the HD 163296 disk make it an excellent observational laboratory to study planet formation.

We present the first detection of ${}^{13}{{\rm{C}}}^{17}{\rm{O}}$ in a protoplanetary disk providing a strong constraint on the CO gas mass in the HD 163296 disk.

HD 163296 was observed with ALMA in Band 7 during Cycle 3 on 2016 September 16 (2015.1.01259.S, PI: S. Notsu). See Notsu et al. (2019) for details on the calibration and self-calibration of the data. The spectral windows have a resolution of 1953.125 kHz. There are 14 ${}^{13}{{\rm{C}}}^{17}{\rm{O}}$ J = 3–2 hyperfine structure lines that lie between 321.851 and 321.852 GHz. All lines lie within a frequency range less than the spectral resolution of the data; hence, we are observing the blending of all of the lines. The ${}^{13}{{\rm{C}}}^{17}{\rm{O}}$ molecular data we use is from Klapper et al. (2003) and was accessed via the Cologne Database for Molecular Spectroscopy (CDMS; Müller et al. 2005). We detected ${}^{13}{{\rm{C}}}^{17}{\rm{O}}$ initially via a matched filter analysis⁸ (Loomis et al. 2018) using a Keplerian mask assuming a disk position angle of 132° and an inclination of 42° (e.g., Isella et al. 2016). The resulting signal-to-noise ratio (S/N) is ≈3.5. The filter response is shown in Figure 1(b) with the black line marking the ${}^{13}{{\rm{C}}}^{17}{\rm{O}}$ J = 3–2 transition after correction for the source velocity (5.8 km s⁻¹).

Zoom In Zoom Out Reset image size

The line imaging was conducted using CLEAN with CASA version 4.6.0. The native spectral resolution of the data is 1.8 km s⁻¹; however, in order to optimize the S/N the final images were generated with a 3 km s⁻¹ channel width and a uv taper of 05 resulting in a synthesized beam of 087 × 051 (100°). Figures 1(c) and (d) present the ${}^{13}{{\rm{C}}}^{17}{\rm{O}}$ integrated intensity map and the intensity-weighted velocity map, respectively. The integrated intensity map was made using channels ±6 km s⁻¹ about the source velocity. The peak integrated intensity is 0.55 Jy beam⁻¹ km s⁻¹ with an rms noise level of 0.08 Jy beam⁻¹ km s⁻¹ (S/N = 7) that was extracted from the spatial region beyond the detected line emission. The intensity-weighted velocity map was made in the same manner but also with a 3σ clip.

We also use archival data to benchmark our modeling, including the ${}^{12}{{\rm{C}}}^{16}{\rm{O}}$, ${}^{13}{{\rm{C}}}^{16}{\rm{O}}$, ${}^{12}{{\rm{C}}}^{18}{\rm{O}}$ J = 2–1 transitions reported in Isella et al. (2016), and the ${}^{12}{{\rm{C}}}^{16}{\rm{O}}$ J = 3–2 ALMA Science Verification data.⁹ All integrated intensity maps were deprojected and azimuthally averaged and are shown in Figures 3(a)–(e). The errors are the standard deviation of intensity of the pixels in each bin divided by the number of beams per annulus (e.g., Carney et al. 2018). In addition, Figure 3(e) shows the ${}^{12}{{\rm{C}}}^{17}{\rm{O}}$ J = 3–2 total integrated intensity value with its associated errors (Qi et al. 2011). All data plotted assume a source distance of 122 pc (van den Ancker et al. 1998). Although GAIA DR2 puts this source at 101.5 pc (Gaia Collaboration et al. 2018), in order to compare to previous analyses we use the previous value. We discuss the impact of the revised distance in Section 4.

Previous observations of the HD 163296 disk with the SMA and ALMA have detected multiple $\mathrm{CO}$ isotopologues: ${}^{12}{{\rm{C}}}^{16}{\rm{O}}$, ${}^{13}{{\rm{C}}}^{16}{\rm{O}}$, ${}^{12}{{\rm{C}}}^{18}{\rm{O}}$, and ${}^{12}{{\rm{C}}}^{17}{\rm{O}}$ (Qi et al. 2011; Isella et al. 2016). The models that were used to reproduce the line emission in Qi et al. (2011) recover the following global isotope ratios:

$$\begin{eqnarray*}\begin{array}{rcl}n{(}^{12}{{\rm{C}}}^{16}{\rm{O}})/n{(}^{13}{{\rm{C}}}^{16}{\rm{O}}) & = & 67\,\pm \,8,\\ n{(}^{12}{{\rm{C}}}^{16}{\rm{O}})/n{(}^{12}{{\rm{C}}}^{18}{\rm{O}}) & = & 444\,\pm \,88,\\ n{(}^{12}{{\rm{C}}}^{18}{\rm{O}})/n{(}^{12}{{\rm{C}}}^{17}{\rm{O}}) & = & 3.8\,\pm \,1.7,\end{array}\end{eqnarray*}$$

where $n{(}^{{\rm{X}}}{{\rm{C}}}^{{\rm{Y}}}{\rm{O}})$ is the number density of the molecule. This is consistent with the carbon and oxygen isotope ratios observed in the ISM (Wilson 1999).

We make a first estimate of the column density of gas traced by the ${}^{13}{{\rm{C}}}^{17}{\rm{O}}$ emission under the assumption of optically thin emission in local thermodynamic equilibrium. Following Carney et al. (2019, their Equation (1), with molecular data obtained from CDMS; Müller et al. 2005), the average column density for the ${}^{13}{{\rm{C}}}^{17}{\rm{O}}$ within 50 au, assuming an excitation temperature of 50 K, is 7.1 × 10¹⁵ cm⁻². This is equivalent to an n_H column density of 2.65 × 10²⁵ cm⁻² (44.4 g cm⁻²) at 50 au. In comparison, the corresponding value for the ${}^{12}{{\rm{C}}}^{18}{\rm{O}}$ is 1.7 × 10¹⁶ cm⁻², resulting in an $n{(}^{12}{{\rm{C}}}^{18}{\rm{O}})$/$n{(}^{13}{{\rm{C}}}^{17}{\rm{O}}$) ratio of 2.5. Under the assumption that both the lines are optically thin (and taking the previously derived isotopic ratios), this value is a factor of 100 too small. Therefore, the ${}^{12}{{\rm{C}}}^{18}{\rm{O}}$ line emission is optically thick and the resulting gas mass derived from this tracer will be underestimated.

To quantify this more robustly, we utilize an existing disk model that has been shown to fit emission lines from multiple CO isotopologues (${}^{12}{{\rm{C}}}^{16}{\rm{O}}$, ${}^{13}{{\rm{C}}}^{16}{\rm{O}}$, ${}^{12}{{\rm{C}}}^{18}{\rm{O}}$, and ${}^{12}{{\rm{C}}}^{17}{\rm{O}}$) to model our ${}^{13}{{\rm{C}}}^{17}{\rm{O}}$ detection (Qi et al. 2011). The density (hydrogen nuclei density, n_H) and temperature of the disk are shown in Figures 2(a) and (b). The CO abundance distribution, shown in Figure 2(c), was determined by setting n(CO) to a constant fractional abundance of 6.0 × 10⁻⁵ with respect to H₂ in the molecular layer following Qi et al. (2011). This assumes that 41% of volatile carbon is in the form of CO (Graedel et al. 1982). This abundance was reduced by a factor of 10⁻⁴ in the midplane where T_gas ≤ 19 K and by a factor of 10⁻⁸ in the atmosphere where the vertically integrated hydrogen column density, σ(n_H), from the disk surface is <1.256 × 10²¹ cm⁻². The depleted value in the midplane is consistent with the CO abundances derived from chemical models including nonthermal desorption (Walsh et al. 2010). The photodissociation and freezeout boundaries are shown in white contours overlaid on Figures 2(a) and (b).

Zoom In Zoom Out Reset image size

The first model that we test, Model 1, uses a constant ${}^{13}{{\rm{C}}}^{17}{\rm{O}}$ fractional abundance of 5.39 × 10⁻¹⁰ relative to H₂. This assumes isotope ratios that are consistent with the observations and modeling from Qi et al. (2011). Model 1 has a total disk mass of 0.089 M_⊙. Using the CDMS data for ${}^{13}{{\rm{C}}}^{17}{\rm{O}}$ we generated a LAMDA-like file in order to model the J = 3–2 hyperfine components in the Line Modeling Engine (LIME¹⁰ ; Brinch & Hogerheijde 2010). Synthetic images cubes were computed assuming the appropriate position angle and inclination of the source, and the resulting images were smoothed with a Gaussian beam to the spatial resolution of the observations using the CASA task, imsmooth. The generated integrated intensity map was then deprojected and azimuthally averaged. The radial profiles from Model 1 (orange) are shown alongside the observations in Figures 3(a)–(e).

Zoom In Zoom Out Reset image size

Model 1 underpredicts the ${}^{13}{{\rm{C}}}^{17}{\rm{O}}$ peak emission in the integrated intensity map by a factor of 2.5, yet provides a reasonable fit to the other lines (within a factor of two). The higher spatial resolution observations are effected by dust opacity within ≈50 au (see Isella et al. 2016); therefore, we focus on reproducing the data beyond 50 au.

To better fit the ${}^{13}{{\rm{C}}}^{17}{\rm{O}}$ observations we globally increase the gas mass of Model 1. This was done by initially multiplying n_H by a factor of 1.5 and then increasing this factor in steps of 0.5 until the best by-eye fit of 3.5 was found. The results for Model 2 are shown in Figure 3 (purple). Model 2 provides a better fit to all of the lines. This model assumes a smooth radial gas density structure contrary to the most recent observations. However, our work is focused on reproducing the global disk mass rather than the underlying small-scale gas surface density variations. Model 2 has a total disk mass of 0.31 M_⊙.

We note that a similar fit can be obtained using a different CO snow-line location at 90 au as determined in Qi et al. (2015). This requires a corresponding increase in gas mass (×3.5) within the snow line, and we obtain a similar ${}^{12}{{\rm{C}}}^{18}{\rm{O}}$ column density profile as in Qi et al. (2015) beyond the snow line. Both of these models use the same underlying physical structure but have different CO snow-line locations and levels of CO depletion beyond the snow line. The Qi et al. (2011) model has simpler assumptions regarding the freezeout of CO, consistent with other work (e.g., Williams & Best 2014), and was found to be a slightly better fit to the observations.

We also generated optical depth maps for a face-on disk to recover the maximum value of τ for each transition. We then radially averaged these maps and plotted the resulting optical depth of the ${}^{12}{{\rm{C}}}^{18}{\rm{O}}$ J = 2–1 and ${}^{13}{{\rm{C}}}^{17}{\rm{O}}$ J = 3–2 transitions for both models in Figure 2(d). It can be seen that ${}^{12}{{\rm{C}}}^{18}{\rm{O}}$ is optically thick within the CO snow line (155 au) in both models, whereas the ${}^{13}{{\rm{C}}}^{17}{\rm{O}}$ remains optically thin across the full radial extent of the disk.

4.1. Comparison to Other Mass Estimates

4.2. The Impact of CO Chemistry on the Disk Mass

4.3. Constraints on the Location of the CO Snow Line

4.4. Is the Disk Gravitationally Stable?

The potential exoplanet population currently probed with ALMA, via the ringed depletion of continuum emission, is gas giant planets on wide orbits. In the case of HD 163296 this would imply a multiple giant planet system and indeed, the presence of such a system has already been proposed (Isella et al. 2016; Liu et al. 2018; Pinte et al. 2018; Teague et al. 2018). The formation of massive planets on wide orbits can in some cases be achieved by core accretion, but a more economical route might involve the gravitational fragmentation of the outer regions of the disk (Boss 2011). Our new, higher disk mass estimate prompts us to investigate whether such processes may have occurred (or be occurring) in the HD 163296 disk.

The stability of a disk against fragmentation can be quantified via the Toomre Q parameter (Toomre 1964):

$$\begin{eqnarray*}&&Q=\displaystyle \frac{{c}_{s}\kappa }{\pi G{\rm{\Sigma }}},\end{eqnarray*}$$

where c_s is the sound speed of the gas, κ is the epicyclic frequency (equal to the angular velocity Ω in a Keplerian disk), and Σ is the surface density of the gas. Toomre Q values of 1 or less imply that the disk is susceptible to fragmentation, but simulations have shown that disks with Q ≲ 1.7 begin to undergo instabilities in the form of nonaxisymmetric spirals (Durisen et al. 2007). We calculate Q across the disk (Figure 4, orange) accounting for the lower mass due to the new source distance, assuming a g/d of 260, and the midplane temperature structure of our model (Figure 2(b)). We find that the minimum value of Q is ≈6 at ≈110 au, suggesting that the disk is currently gravitationally stable (in agreement with recent work from Powell et al. 2019).

Zoom In Zoom Out Reset image size

The relatively large age of the HD 163296 system brings into question its stability earlier in its lifetime. The determination of previous disk masses is complicated by processes including episodic accretion (e.g., Mendigutía et al. 2013) and the decrease in accretion rate with time (e.g., Venuti et al. 2014). The magnitude of these effects are still under debate (see Hartmann et al. 2016 for a review), so we therefore assume that all of the accreted mass once resided in the disk, and that the accretion rate ($\dot{M}$) has been constant over the disk lifetime. HD 163296 has an estimated age of ${6.03}_{-0.14}^{+0.43}$ Myr and $\mathrm{log}\dot{M}=-{6.81}_{-0.15}^{+0.16}$ M_⊙ yr⁻¹ (C. Wichittanakom et al. 2019, in preparation); thus, under these assumptions, we estimate the disk mass at 0.1 Myr to be ${1.13}_{-0.28}^{+0.51}$ M_⊙.

The resulting minimum Toomre Q values for this star–disk configuration¹¹ would be in the range of 1.3–0.7 (Figure 4, purple), placing regions of the disk from ∼50–220 au in the regime of instability. Such behavior early in the disk lifetime has implications for the trapping and growth of dust (Rice et al. 2004) and the chemical composition of the disk (Evans et al. 2015). This previous unstable state could also be the source of the four massive planets currently proposed to reside in the disk around HD 163296.

We have presented the first detection of ${}^{13}{{\rm{C}}}^{17}{\rm{O}}$ in a protoplanetary disk showcasing the power of this optically thin isotopologue as a tracer of disk gas mass. This work provides robust evidence that disks are more massive than previously assumed (see also Zhu et al. 2019). Future observations of this tracer in more sources may help to address the discrepancy between the masses of disks and the observed exoplanet population (Manara et al. 2018).

We thank an anonymous referee for constructive comments that improved the clarity of several sections of the paper. We thank Andrea Isella for the ${}^{12}\mathrm{CO}$, ${}^{13}\mathrm{CO}$, and ${{\rm{C}}}^{18}{\rm{O}}$ J = 2–1 data that was vital in our model comparison. We also thank Chumpon Wichittanakom for providing us with updated stellar parameters for HD 163296 prepublication of their paper. A.B. acknowledges the studentship funded by the Science and Technology Facilities Council of the United Kingdom (STFC). C.W. acknowledges funds from the University of Leeds. J.D.I. and C.W. acknowledge support from the STFC under ST/R000549/1. S.N. is grateful for support from JSPS (Japan Society for the Promotion of Science) Overseas Research Fellowships and the ALMA Japan Research Grant of NAOJ Chile Observatory, NAOJ-ALMA-211. This work is supported by MEXT/JSPS KAKENHI grant Nos. 16J06887, 17K05399, 19K03910. Part of ALMA Data analysis was carried out on the Multi-wavelength Data Analysis System operated by the Astronomy Data Center (ADC), National Astronomical Observatory of Japan. This paper makes use of the following ALMA data: 2011.0.00010.SV, 2013.1.00601.S, and 2015.1.01259.S. ALMA is a partnership of European Southern Observatory (ESO) (representing its member states), National Science Foundation (USA), and National Institutes of Natural Sciences (Japan), together with National Research Council (Canada), National Science Council and Academia Sinica Institute of Astronomy and Astrophysics (Taiwan), and Korea Astronomy and Space Science Institute (Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, Associated Universities, Inc/National Radio Astronomy Observatory (NRAO), and National Astronomical Observatory of Japan.