The TW Hya Rosetta Stone Project I: Radial and vertical distributions of DCN and DCO⁺

Figure 1 shows integrated emission (moment-0) maps and radial profiles of the DCN and DCO⁺ J=2–1, 3–2, and 4–3 rotational line emission toward TW Hya. The maps are constructed by integrating emission across all channels that show signal above 3$\sigma$ (Table 2). We do not detect any molecular line hyperfine structure, i.e., it is either too weak or unresolved. The radial profiles are derived from the moment-0 maps by azimuthally integrating the inclination-corrected maps ($i=7^{\circ}$ and P.A.=$155^{\circ}$) in narrow rings. These values are slightly different from those used by (Huang et al., 2018), and were selected because they provided the best visual fit to the data when overplotting a Keplerian model on top of channel maps (see Appendix A). Displayed uncertainties are extracted using the rms per beam in the maps and taking into account the number of independent beams in each ring.

All DCO⁺ and DCN moment-0 maps show qualitatively similar central depressions or holes in the line emission. The central holes are not identical for DCO⁺ and DCN, however. Figure 2 shows an overlay of the higher resolution 4–3 lines, which indicate that the DCN hole has a somewhat smaller radius. On closer inspection, it seems like DCN and DCO⁺ share a radial component that peaks at $\sim$0$\farcs$6. Interior to this, DCN has a shoulder at $\sim$0$\farcs$4. Exterior to 0$\farcs$6 the emission pattern clearly differs for DCN and DCO⁺. The DCN emission decreases quickly with radius, followed by a low-intensity halo stretching out to larger radii. DCO⁺ displays a broad second component that encompasses the DCN halo. This second DCO⁺ component peaks at $\sim$1″. None of these features correspond to previously noted dust rings, but we note that the DCN shoulder at $\sim$0$\farcs$4 or $\sim$25 au coincides with a dust gap , and a break in the ¹²CO radial emission profile (Huang et al., 2018). Furthermore, the DCN peak nearly coincides with a second dust gap at 41 au ( $\sim$0$\farcs$69). The precise dust and gas properties of these gaps are unknown, and it is therefore unclear if the dust gaps are causing an increased DCN emission. One possible causal connection is that dust gaps likely present increased UV penetration, and therefore increased dissociation of CO in carbon atoms, which could fuel nitrile production. If this is the case, we would expect future observations to reveal excess emission of other nitriles at the same locations.

The emission morphologies of each molecule are consistent across the different transitions. This consistency indicates that all lines of each molecule likely originate in the same disk regions, and that observed emission substructures trace changes in column density, and not only changes in excitation conditions across the disk. We cannot rule out that some of the missing DCO⁺ and DCN emission towards the disk center is due to continuum opacity, since Huang et al. (2018) finds that the continuum $<$20 au may be, in part, optically thick. Continuum opacity cannot be the whole explanation for these central cavities, however, since other optically thin lines, including ¹³C¹⁸O, present centrally peaked emission (Zhang et al., 2017). This suggests a chemical cause for the central cavities, for example the turn-off of the dominant deuterium fractionation pathway, and this is discussed in §5.

Table 2 lists the integrated emission across the entire DCO⁺ disk (out to 2$\farcs$5), and within a radius of 1$\farcs$7 or 100 au (where most DCO⁺ and DCN emission originates). The listed uncertainties are based on the measured rms per beam in line emission-free moment-0 maps, produced from line-free channels in the relevant spw, multiplied by the square root of the number of beams within $r=$2$\farcs$5 and 1$\farcs$7, respectively. We suspect that the DCN 2–1 line emission is underestimated due to a combination of unaccounted emission from hyperfine line emission, and incomplete flux recovery of this low SNR line.

Figure 3 shows the extracted spectra of the DCN and DCO⁺ lines using Keplerian masks (Pegues et al., 2020) to enhance the SNR. The 4–3 and 3–2 lines show Keplerian profiles, while the resolution of the 2–1 lines are too poor to resolve the characteristic double-peak structure. The DCN lines are broader than the DCO⁺ lines, which has two origins: unresolved hyperfine structure, and more emission emitting at smaller disk radii. We also inspected spectra from individual pixels and saw no evidence for substantial non-Keplerian motion or line self-absorption.

We begin our characterization of DCO⁺ and DCN using disk-averaged rotational diagrams. We use the DCN and DCO⁺ integrated fluxes within 100 au, which encompasses the main emission features. We exclude the DCN 2–1 emission from the rotational diagram analysis because it is likely underestimated (see discussion above) – if it is included, the derived temperature is above 100 K and the fit to the other lines is poor. The main flux uncertainty is a 10% absolute flux calibration uncertainty, which is added in quadrature to the rms-based integrated emission errors when calculating the rotational diagram. The molecular line data was all taken from CDMS and is listed in Table 2. We used the following partition functions (also from CDMS) calculated at [0, 9.375, 18.75, 37.5, 75, 150, 225, 300] K: [0, 5.769, 11.1866, 22.0293, 43.7220, 87.1365, 130.5570, 173.9803] and [0, 17.2240, 33.3906, 65.7550, 130.5095, 262.2213, 409.8604, 586.3727] for DCO⁺ and DCN, respectively.

To calculate the rotational diagrams, we follow the MCMC procedure outlined in Loomis et al. (2018a), using the emcee package (Foreman-Mackey et al., 2013). Figure 4 shows rotational diagrams corresponding to random draws from the posterior probability distributions of the excitation temperatures and column densities, with the optical depth corrected values of $N_{\rm u}/g_{\rm u}$ plotted against $E_{\rm u}$. The disk-averaged column densities of both molecules are similar: 2.6 and 3.9$\times 10^{12}$ cm^-2, respectively for DCN and DCO⁺. DCO⁺ appears colder than DCN (33 vs 66 K), but we note that the fit to the DCN data is very uncertain since we had to remove the 2–1 line, and we cannot exclude that the actual DCN excitation temperature is $<$40 K.

We next use DCO⁺ 4–3, 3–2, and 2–1 radial emission profiles (Fig. 5, top panel) to carry out the same rotational diagram procedure in radial bins across the disk. The radial profiles were generated from moment-0 maps with a resolution of 0$\farcs$5, and hence appear smoother than the profiles in Fig. 1. Similar to the disk-averaged rotational analysis, the absolute calibration uncertainty was added in quadrature to the rms-based uncertainty shown in Fig. 5. The middle and bottom panels of Fig. 5 show the resulting radial temperature and column density profiles. DCO⁺ appears to display a rapid decrease in excitation temperature within the inner gap, i.e. out to $\sim$25 au, but the emission levels within the gap are low and the uncertainties too high to claim a certain trend. Between 25–70 au, the best-fit DCO⁺ excitation temperatures are 37–39.5 K, suggesting that in the bulk of the disk, DCO⁺ originates in a layer that is more elevated than the CO snow surface, which is expected at $\sim$20–25 K. Beyond 70 au, which coincides with the outer edge of the pebble disk (Andrews et al., 2012, 2016; Huang et al., 2018), the DCO⁺ temperature begins to drop and reaches $<$30 K at 90 au. At even larger radii there is a possible temperature inversion, but the SNR is too low to determine that this is real. We note that the lower DCO⁺ excitation temperature in the outer disk has a disproportional impact on the disk-averaged excitation temperature due to the relatively larger emission area of the 70–120 au portion of the disk compared to the inner 70 au, which explains why the disk-averaged DCO⁺ excitation temperature is lower than the ‘characteristic’ excitation temperature.

The rotational diagram results have two important caveats. First the emission surfaces of the different DCO⁺ transitions are not necessarily the same, and this may skew the temperature upwards. We explore below whether the DCO⁺ emission could be consistent with a colder origin. Second, the rotational diagram analysis, as well as the parametric models below assume LTE, and if some of the DCO⁺ is originating from very elevated disk layers this assumption may not hold. While we cannot exclude non-LTE excitation, we note that it typically results in sub-thermally excited lines and therefore an under-prediction of the kinetic temperature, and should therefore not impact the main result here that DCO⁺ appears to originate from a warmer disk layer than expected.

The first set of models focuses on the inner cavity seen in both DCO⁺ and DCN emission. The first abundance model (A1) has a constant abundance exterior to 25 au, and five orders of magnitude lower abundance in the inner hole. The abundance exterior to 25 au is estimated by eye to fit the 4–3 lines. Figures 7 (first panel), and 8 (first panel) show that this prescription does not provide a good fit to either DCO⁺ or DCN emission, but for different reasons. In the case of DCO⁺, the inner part of the emission profile is well fit by Model A, while the outer disk is not, indicative of that the DCO⁺ abundance is higher in the outer than inner disk. In the case of DCN, this model produces too little emission in the hole, and too much emission at large radii, beyond 60 au.

We first address the DCO⁺ discrepancy by exploring whether changing the abundance profile exterior to the hole from constant to an increasing power law provides a better fit (Model A2). The second panel of Fig. 7 shows that a DCO⁺ abundance power law with a power law index of 0.75 does indeed provide a reasonable fit to the DCO⁺ emission profile when imaged at this resolution. Note, however, that our higher resolution data shows a DCO⁺ double-peak, rather than a broad single ring, which could not be explained by such a continuous power law. We also note that the model that fits the 4–3 data predicts a 2–1 emission level that is close to the one observed. This shows that the our observations do not rule out the presence of DCO⁺ close to the disk midplane, as long as there is a substantial amount of DCO⁺ in the warm, upper disk layers.

The DCN A1 discrepancy indicates that the DCN abundance is lower in the outer than the inner disk if DCN is emitting from all disk layers. One possible explanation is that DCN is only emitting from warmer gas, of which there is a limited amount in the outer disk. To test this, we apply a 20 K temperature boundary to the A1 model (Model A3), i.e. a constant abundance exterior to 25 au wherever the temperature is $>$20 K, and a five orders of magnitude lower abundance everywhere else. In our disk model, the midplane drops below 20 K at 44 au. The second panel of Fig. 8 shows that this model fits the 4–3 data quite well in the outer disk, but naturally does not fix the underabundance noted towards the disk center, which could be addressed either by making the inner cavity $\sim$5 au smaller or by implementing a smaller DCN depletion factor (not shown). The 2–1 data is always overpredicted, indicative of that the DCN is even warmer.

In a second set of of models ,we apply temperature boundaries, rather than radial cut-offs, with the aim of exploring whether the observed emission profiles can be explained by DCN and DCO⁺ temperature-dependent formation alone. The first temperature models (B1a and B1b) assume constant abundances within lower and upper temperature bounds across the disk. We explored several different boundaries, and found that the inner radial profile requires a maximum temperature cut-off at 25–30 K; 27 K provides the best fit for DCO⁺ (B1a), and 30 K for DCN (B1b). We set the temperature minimum cut-off of 20 K, based on model predictions for CO freeze-out. The B1 model provides a good fit to the DCN 4–3 data, but over-predicts the 2–1 line emission. The 2–1 emission is also over-predicted for DCO⁺ in the inner disk by 20–30%. This is qualitatively consistent with the results from the rotational diagram analysis, which indicated that a substantial amount of DCO⁺ originates from temperatures above 30 K in the inner 70 au of the disk.

In addition to the mismatch between the relative levels of 2–1 and 4–3 emission, the B1a model under-predicts the DCO⁺ emission in the outer disk. We explore whether this mismatch can be explained by a second DCO⁺ component around the dust edge, where UV photons may penetrate deeper into the disk, bringing excess cold CO into the gas-phase. We modify B1a such that exterior to 60 au (the pebble disk edge is around 70 au), DCO⁺ is present $<$27 K, while interior to 60 au, the B1a model boundaries of 20$<$T$<$27 K still apply. This results in a more correct shape of the radial profile and may also explain the presence of a double-peaked DCO⁺ radial profile. We emphasize that neither the B1a/b or B2 models correctly predict the 4–3/2–1 line ratios in the inner disks for DCO⁺ and DCN. Since the B1 and B2 temperature boundaries simulate a cold emission layer, this mismatch suggests that neither molecule is mainly present in the cold midplane.

In the inner regions of the TW Hya disk, both DCN and DCO⁺ emission (and column densities) increase rapidly with increasing radius starting around 20–25 au, corresponding to midplane temperatures of 27–30 K. In the absence of multi-line observations, this observation would have been in line with expectations, since if DCO⁺ forms from reactions between the cold gas tracer H₂D⁺ and CO, DCO⁺ should be most abundant in the midplane just interior to the CO snowline (e.g. Mathews et al., 2013; Aikawa et al., 2018). However, the inferred warm DCO⁺ excitation temperature at 25 au shows that DCO⁺ (and DCN) cannot be primarily emitting from the midplane. Instead the measured excitation temperature of $\sim$40 K at $25-60$ au places a substantial amount of the DCO⁺ in an elevated disk layer. This discovery could point to a relatively inefficient low-temperature H₂D⁺ fractionation chemistry, and an efficient deuterium enrichment through reactions with, e.g., CH₂D⁺. The latter is expected to proceed at higher temperatures than the H₂D⁺ chemistry and may therefore be consistent with a 40 K excitation temperature (Wootten, 1987; Parise et al., 2009; Favre et al., 2015; Roueff et al., 2015).

This proposed scenario presents a new puzzle, however, which is that we do not observe much DCN and DCO⁺ in the innermost midplane region. If we are observing a deuterium fractionation chemistry that is active at $\sim$40 K, we would naively expect abundant DCO⁺ and DCN in the inner disk midplane regions that correspond to this temperature, roughly $\sim$10-20 au. Exterior to the DCO⁺ and DCN cavity, a possible solution to the observed emission pattern is that there are two DCO⁺ and DCN production zones, cold and warm, respectively, which when observed towards a face-on disk masquerade as the luke-warm emission layer we observe. Indeed the A1 model in §4 shows that this is a possibility from an excitation point of view. Whether this scenario is chemically plausible is less clear, and we still need an explanation why DCO⁺ in TW Hya appears much warmer than in e.g. the HD 162936 disk (Flaherty et al., 2017). In the inner disk the lack of DCO⁺ and DCN may in part be explained by continuum blocking out some fraction of the molecular emission, making the hole seem deeper than it really is. However, as discussed previously, we do not think that it is likely that continuum opacity alone is responsible for the central cavity. Instead we suggest that the lack of DCO⁺ and DCN in the inner disk, and the elevated temperature and therefore elevated location of DCO⁺ exterior to 25 au, together point towards a relatively inefficient DCO⁺ and DCN production in the disk midplane at all disk radii.

One possible explanation for the lower than expected DCO⁺ midplane abundances is that CO is depleted in the disk far beyond the CO freeze-out zone due to e.g. chemical processing and diffusion (e.g. Meijerink et al., 2009; Xu et al., 2017; Schwarz et al., 2018). If CO depletion begins at 40 K instead of 25 K, this would explain why DCO⁺ production is low both in the inner disk midplane, and close to the CO snow surface in the outer disk. This explanation is supported by observations that show that the TW Hya disk is very CO-depleted throughout the disk molecular layers including in the inner disk (Favre et al., 2013; Kama et al., 2016; Schwarz et al., 2016; Zhang et al., 2019). CO depletion could also diminish DCN formation in the same locations, since CO depletion from the gas would result in a depletion of the overall carbon reservoir that controls DCN and HCN production. However, there is both observational (Hily-Blant et al., 2010) and theoretical (Long et al. subm.) evidence that HCN (and by extension DCN) is not very sensitive to CO depletion, and this idea should therefore be considered highly speculative. DCN may be present at elevated disk layers simply because it mainly forms through the CH₂D⁺ pathway.

If CO depletion controls where in a disk there is an active deuterium fractionation chemistry, the TW Hya results may be far from universal. CO depletion through either chemistry or diffusive flows is expected to become more severe with disk age, and TW Hya has an unusually old disk. By contrast, we may expect to find colder and more mid-plane oriented DCO⁺ in disks that are less depleted in CO. Interestingly in the disk around Herbig Ae star HD 163296, which has been shown to be much less depleted in CO than TW Hya (Zhang et al., 2019), the DCO⁺ excitation temperature is low ($<$20 K) (Flaherty et al., 2017). By contrast, most DCO⁺ in the disk around Herbig Ae star HD 169142, appears to be warm (Carney et al., 2018). Additional observations towards a sample of young and old T Tauri and Herbig Ae disks would be key to resolve if the DCO⁺ (and DCN) chemistry migrates to elevated disk layers over time, and if there are systematic differences between T Tauri and Herbig Ae disks.

A second possible explanation for the inferred low levels of DCO⁺ in the TW Hya disk midplane is that the disk midplane regions are chemically quenched due to a lack of ionization. Close to the CO snow surface, there may be too little ionizing radiation to drive a H₂D⁺- or CH₂D⁺-mediated chemistry, and the measured excitation temperature of DCO⁺ may reflect the coldest disk layer where an ion-molecule mediated deuterium fractionation chemistry is efficient. TW Hya has been inferred to have a low level of ionization throughout most of the disk (Cleeves et al., 2015), which supports this scenario. If this is the primary explanation for the lack of DCO⁺ and DCN in the midplane, we would expect to see a decreasing DCO⁺ temperature with radius, since ionization should increase in the outer, more tenuous disk regions. Indeed such a decrease is detected exterior to 60 au, but we note that there are also other possible reasons for this decrease in DCO⁺ temperature, including release of cold CO into the gas-phase through photodesorption (Öberg et al., 2015b). Additional high-resolution DCO⁺ and DCN observations towards disks with estimated ionization levels are needed to test this hypothesis.

Finally, we note that it is an open question whether we should also expect DCO⁺ and DCN in the warm disk atmosphere in the inner disk. Favre et al. (2015) predicts that DCO⁺ should form abundantly in this disk region, while Aikawa et al. (2018) have come to a different conclusion. More theoretical work is needed to resolve this, but in the meantime we note that in the case of TW Hya, there is no evidence that the inner disk atmosphere is an important source of deuterated molecules.

The above discussion is relevant for both DCN and DCO⁺. We now proceed with exploring reasons for observed differences between the two molecules. First, based on emission profiles and toy models, the DCN cavity is somewhat smaller and/or less empty than the DCO⁺ one. This suggests that there is at least one warm deuterium fractionation pathway that mainly affects DCN. There is tentative evidence for this warmer formation channel is important for DCN throughout the disk, since the disk-averaged DCN excitation temperature appears higher than that of DCO⁺. This, however, needs to be revisited with deeper DCN 2–1 observations.

In the outer disk of TW Hya, the DCN and DCO⁺ radial profiles also diverge. While DCN presents a halo exterior to the pebble disk emission, DCO⁺ has a much more substantial second emission component close to the edge of the pebble disk. Similar differences have been seen in other disks, most notably towards IM Lup and HD 163296 (Huang et al., 2017; Öberg et al., 2015b; Salinas et al., 2017). The origin of a DCO⁺ peak at the pebble disk edge is likely a result of increased penetration of UV radiation in the less shielded outer disk regions, which results in CO sublimation due to either a thermal inversion (Cleeves, 2016), or enhanced CO ice photodesorption (Öberg et al., 2015b; Huang et al., 2016; Aikawa et al., 2018). Cold, CO-rich gas constitute an ideal environment for DCO⁺ formation through the H₂D⁺ channel. DCN clearly requires something in addition to this to form efficiently, though the presence of the DCN halo suggests that a small amount of DCN also forms under these conditions. One possible avenue to test whether the DCO⁺ and DCN in the outer disk originates with H₂D⁺ would be to add observations of N₂D⁺. N₂D⁺ only forms through reactions with H₂D⁺ and could therefore be used to map out where this pathway is active (see e.g. Pagani et al., 2007; Salinas et al., 2017; Aikawa et al., 2018; Caselli et al., 2019). An important complication, is that N₂D⁺ is only expected where there is substantial CO freeze-out and not seeing N₂D⁺ can therefore not be used to rule out the H₂D⁺ pathway.

In summary, there is evidence for active deuterium fractionation chemistry in the TW Hya disk. However, much of the observed emission from DCO⁺ and DCN appears to originate well above the CO snow surface, and some process, perhaps CO depletion or low disk midplane ionization, may be limiting the efficiency of the cold pathway in the lower disk layers and in the inner disk midplane.

DCO⁺ and DCN have been observed and characterized towards TW Hya in a number of earlier studies. For reference, we extracted disk average column densities of DCO⁺ and DCN of 3.9 and 2.6$\times 10^{12}$ cm^-2 respectively, and a DCO⁺ peak column density of $\sim$7$\times 10^{12}$ cm^-2. Both are substantially higher compared to values derived from single dish observations of 3 and $<$0.4$\times 10^{11}$ cm^-2 for DCO⁺ and DCN, respectively (van Dishoeck et al., 2003; Thi et al., 2004). This difference can likely be explained by beam dilution in the single dish observations. The large difference in DCN and DCO⁺ column densities inferred from single dish observations is probably a beam dilution effect as well, since we find DCN to be more compact than DCO⁺.

DCO⁺ and DCN have also been marginally resolved by Qi et al. (2008) and Öberg et al. (2012), and these observations were used to derived radial column density profiles. Qi et al. (2008) found a peak column density of $\sim$4$\times 10^{12}$ cm^-2, close to our measurement. By contrast the estimates of the DCN disk averaged column density in Qi et al. (2008) and Öberg et al. (2012), are an order of magnitude lower than we find here. Some of this may be explained by beam averaging, since the synthesized beam in Öberg et al. (2012) was larger than the resolved DCN emitting region. The remaining difference can probably be accounted for by different disk temperature structure and DCN emitting layer assumptions.

DCO⁺ column densities and abundances have also been estimated towards a handful of other disks and the results are remarkably similar to those we find towards TW Hya; Teague et al. (2015) and Qi et al. (2015) found a DCO⁺ column densities towards DM Tau and HD 163296 of $\sim 10^{12}$ cm^-2, and Carney et al. (2018) and Salinas et al. (2018) found DCO⁺ abundances with respect to hydrogen towards HD 169142 and HD 163296 of $0.9-1.5\times 10^{-12}$ and $2-6\times 10^{-12}$, respectively. The consistent DCO⁺ column densities and abundances towards this sample of four disks is difficult to interpret, since the DCO⁺ emitting layer appears quite different in e.g. TW Hya and HD 163296. Finally, (Salinas et al., 2017) also estimated the DCN abundance in the HD 163296 disk, and found $\sim$10^-12 per hydrogen nuclei, which is again consistent with TW Hya.

Model predictions of DCN and DCO⁺ column density profiles, abundances, and emitting layers go back to the early 2000s. In a majority of models, DCO⁺ column densities across disks are $\sim 10^{12}$ cm^-2, in good agreement with the TW Hya findings (Aikawa & Herbst, 2001; Willacy, 2007; Aikawa et al., 2018), but as discussed below this may be a coincidence since the DCO⁺ distributions in TW Hya and in fiducial model disks appear quite different. A notable exception is Favre et al. (2015), who predicted substantially higher column densities due to efficient warm DCO⁺ formation. In contrast to our findings, DCN is predicted to be at least an order of magnitude less abundant than DCO⁺ in most models (Aikawa & Herbst, 2001; Aikawa et al., 2002; Willacy, 2007; Favre et al., 2015). The one exception is one model in Willacy (2007), which predicts similarly high DCN and DCO⁺ column densities. This model includes efficient ice photodesorption, which both increases the overall gas-phase carbon reservoir, and desorbs some of the DCN that forms through grain surface chemistry in their model. We speculate that the difference between models and observations with regard to the relative DCO⁺ and DCN abundances may be due to a high C/O ratio in the TW Hya disk, which would enhance both HCN production and the importance of the CH₂D⁺ fractionation pathway.

Models also predict shapes of radial profiles. In most models, DCO⁺ displays a prominent inner hole, while DCN does not (Aikawa & Herbst, 2001; Willacy, 2007; Willacy & Woods, 2009; Aikawa et al., 2018). This results in different DCO⁺ and DCN radial profiles across the disk, in contrast to what is observed in both TW Hya and HD 163296, where DCN and DCO⁺ appear to coincide at intermediate disk radii. This mismatch between models and observations suggests that the relative contributions of cold and warm deuterium fractionation pathways to DCO⁺ and DCN remain to be fully worked out in disks.