© ESO 2018

1 Introduction

DCO⁺ is a good tracer of thedeuterium fractionation and ionization fraction of low temperature environments (Millar et al. 1989; Favre et al. 2015). Detections of DCO⁺ and other simple deuterated molecules towards protoplanetary disks are present in only a handful of T Tauri disks: TW Hya (van Dishoeck et al. 2003), DM Tau (Guilloteau et al. 2006; Teague et al. 2015), AS 209, IM Lup V4046 Sgr, and LkCa 15 (Öberg et al. 2010, 2011, 2015; Huang et al. 2017), and in the disks around the Herbig Ae stars MWC 480 (Huang et al. 2017) and HD 163296 (Qi et al. 2008, 2015; Mathews et al. 2013; Yen et al. 2016; Salinas et al. 2017). High angular resolution observations of some of these disks have revealed a surprisingly complex radial structure. The chemistry involved in the gas-phase formation of DCO⁺ is thought to be well understood and has been previously studied in the interstellar medium (ISM) and in disks (Turner 2001; Willacy 2007; Roueff et al. 2013; Favre et al. 2015). In this paper, we attempt to determine whether our current understanding of the DCO⁺ chemistry is sufficient to reproduce the complex radial structure seen in protoplanetary disks, particularly in the disk surrounding the Ae Herbig star HD 163296.

DCO⁺ forms in thegas phase through two different regimes: a low temperature deuteration (hereafter cold deuteration, CD) and a warm deuteration (hereafter warm deuteration, WD) channel. The CD and WD regimes are gradually enhanced at temperatures lower than ~30 and ~80 K, respectively (Millar et al. 1989; Albertsson et al. 2013). This enhancement is a direct consequence of the lower zero-point vibrational energy of simple deuterated molecules in comparison to their non-deuterated counterparts. The origin of highly deuterated species and of DCO⁺ in the solar nebula can be attributed to in situ synthesis in the primordial disk or to being inherited from the ISM. Deuterium is injected intothe chemistry via ion-molecule reactions, and is kept in it by the endothermic nature of the inverse reaction. For an effective enhancement of the deuterium fractionation ratio, the environment must be cold (typically tens of degrees Kelvin) and ionized.In the dense shielded ISM, the ionization environment is mainly dominated by galactic cosmic rays (CRs). In the cold and dense outer regions of protoplanetary disks where deuterium enrichment takes place, the ionization source comes from galactic cosmic rays (CRs) and to a lesser extent from radiation, both from the host star and from external sources in low-density regions above the midplane where the ultraviolet (UV) radiation can penetrate.

DCO⁺ was first proposed as a good CO snowline tracer in protoplanetary disks (Mathews et al. 2013). Recent chemical models (including the WD channel) conducted by Favre et al. (2015) have proposed DCO⁺ as a good tracer of the ionization degree in the inner regions of planet-forming disks rather than the temperature structure and the CO snowline. The WD channel formation channel involves ionized simple hydrocarbons, unlike the CD involving H₂ D⁺, which is highly reactive with CO. If the CD dominates the formation of DCO⁺ in disks, then it can be used as an indirect tracer of the CO snowline because its parent molecule, H₂ D⁺, is readily destroyed by CO in the gas phase. The models of Favre et al. (2015) show that the WD channel enhances the column density of DCO⁺ by a factor of 5 in the warm regions of a T Tauri-like disk, where CO is still in the gas phase, and is responsible for the bulk of the abundance. Observations of significant emission of DCO⁺ in the inner parts of protoplanetary disks have been already reported (Qi et al. 2015; Huang et al. 2017). In particular, Salinas et al. (2017) have seen DCO⁺ extending from ~50 to ~300 AU in the disk surrounding the Herbig Ae star HD 163296.

HD 163296 has a massive (0.089 M_⊙) inclined (44°) disk and its gaseous content, probed by ¹²CO emission, extends at least to ~500 AU (Qi et al. 2013; Mathews et al. 2013). These attributes and its proximity (122 pc, van den Ancker et al. 1998) make it an excellent candidate to study both the radial and vertical distribution of deuterated species such as DCO⁺. The disk has been observed in the millimeter regime revealing a distribution of mm-sized grains that extend up to ~230 AU with multiple rings and gaps (Isella et al. 2016).

The DCO⁺ radial distribution of the HD 163296 disk has been characterized by Salinas et al. (2017) using ALMA observations. The data are consistent with three regimes of different constant abundances defined by one inner radius at 50 AU and two breaks at 120 and 245 AU. They found that the first two regimes correlate well with the expected WD and CD channels traced by DCN and N₂ D⁺. The third regime correlates with the extent of the mm-sized dust grains which hints at a local decrease in UV opacity allowing photodesorption of CO and consequent DCO⁺ formation as proposed for the disk around IM Lup (Öberg et al. 2015). An interesting alternative to explain this excess DCO⁺ emission is releasing CO through thermal desorption (Cleeves 2016). Dust grain evolution models by Facchini et al. (2017) have predicted a thermal inversion of the dust temperature as a direct consequence of radial drift and settling in the disk surrounding HD 163296 given low turbulence values.

Our main goal is to implement a simple chemical network for the CD channel and a parametrized WD channel to reproduce both the location and the amount of the observed DCO⁺ in HD 163296. Specifically, we aim to constrain the relative contribution of the WD channel after taking into account the CD channel. Although several other effects can produce the observed ringlike feature at 245 AU and drop off at larger radii of the DCO⁺ emission, this paper only focuses on the chemical conditions that can lead to such structures. A fully self-consisting modeling including an accurate determination of the temperature profile, dust properties, and gas density profile is needed to distinguish between the possible origins of the DCO⁺ morphology.

In Sect. 2 we briefly describe the data and explain our modeling strategy. Section 3 contains theresults obtained from our different modeling approaches. Section 4 discusses the validity of these models and provides an interpretation of the model parameters. Finally, in Sect. 5 we summarize our findings and conclusions.

2 Methods

2.1 Previous observation of HD 163296

This study uses data of DCO⁺ J = 3–2 at 216.112 GHz in the disk surrounding the Herbig Ae star HD 163296 (α₂₀₀₀ = 17^h56^m51^s.21, δ₂₀₀₀ = − 21°57′22. ′′0) obtained by the Atacama Large (sub)Millimeter Array (ALMA) in Band 6 as a part of Cycle 2 on 2014 July 27–29 (project 2013.1.01268.S). The spectral resolution is 0.7 MHz corresponding to 0.085 km s⁻¹ with respect to the rest frequency of the line. The uv-coverage of the data ranges from 20 to 630 kλ.

The data were reduced as described by Salinas et al. (2017) and Carney et al. (2017). The DCO⁺ J = 3–2 line was continuum subtracted in the visibility plane using a first-order polynomial fit and imaged using a Briggs weighting of 0.5. The resulting synthesized beam has dimensions of 0.′′53 × 0.′′42 and is shown along with the obtained channel maps in Appendix A. Salinas et al. (2017) use a Keplerian mask to obtain an integrated intensity map. A radial emission profile can be constructed by taking the average value of concentric ellipsoids, centered at the star, in the integrated intensity map. The resulting DCO⁺ radial emission profile can be seen in Fig. 1.

The observational analysis of Salinas et al. (2017) modeled DCO⁺ as a three-region radial structure. They argue that these radial regions correspond roughly to the WD channel, the CD channel, and a third unidentified regime at larger radii (≳240 AU). Their analysis treated each region independently and assumed that the bulk of the DCO⁺ emission comes from the midplane. These assumptions do not impact their ability to constrain the regimes radially, but do not allow us to constrain the contribution of the WD channel in the HD 163296 disk. Salinas et al. (2017) only derived vertically averaged DCO⁺ abundances, but here we carry out full 2D modeling (radius and height) of the chemistry. Because deuteration, and the formation of DCO⁺, is a temperature dependent process a vertical treatment of its abundance is needed in disks, due to their strong vertical temperature gradient. For our analysis, we assume that the WD channel remains active where the CD channel operates, but its contribution to the DCO⁺ emission in the outer regions of the disk is negligible. If this is the case, reproducing the emission at large radii alone can provide a constraint on the amount of DCO⁺ produced bythe CD channel.

Fig. 1 Observed DCO+ radial profile (black continuous line). These profiles are obtained by averaging the values of concentric ellipsoids of the integrated intensity maps. The error bars of the observed DCO+ radial profilecorresponds to 3σ, where σ is the standard deviation of the values contained in the ellipsoid divided by the square root of the number of beams. The 50 AU bar corresponds to the semiminor axis of the synthesized beam and serves as a measure of the spatial resolution. The dashed lines show radial profiles of the emission from models with two different evaporation temperatures of CO TCO = 19–20 K and a constant abundance Xin = 5.0 × 10−5.

2.2 Chemical model

Deuterium is incorporated into the gas chemistry mainly through the following ion molecule reactions (Gerner et al. 2015; Turner 2001) $${\text{H}}_{\text{3}}^{\text{+}}+\text{HD\hspace{0.17em}}\leftrightharpoons {\text{H}}_{\text{2}}{\text{D}}^{+}{+\text{H}}_{\text{2}}+\text{230\hspace{0.17em}K},$$(1a) $${\text{CH}}_{\text{3}}^{\text{+}}+\text{HD}\leftrightharpoons {\text{CH}}_{\text{2}}{\text{D}}^{+}{+H}_{\text{2}}+\text{370\hspace{0.17em}K},$$(1b) $${\text{C}}_{\text{2}}{\text{H}}_{\text{2}}^{+}+\text{HD}\leftrightharpoons {\text{C}}_{\text{2}}{\text{HD}}^{+}{+\text{H}}_{\text{2}}+\text{550\hspace{0.17em}K}.$$(1c)

The right-to-left reactions of Eqs. (1a), (1b), and (1c) are endothermic and effectively enhance deuterium fractionation in low temperature environments. Equation (1a) corresponds to the CD channel and are active at temperatures ranging from 10 to 30 K (Millar et al. 1989; Albertsson et al. 2013). Equations (1b) and (1c), involving light hydrocarbons, correspond to the WD channel and are active at warmer temperatures ranging from 10 to 80 K.

DCO⁺ is formed in the gas phase involving both the CD and WD channels through the reactions (Watson 1976; Wootten 1987; Favre et al. 2015) $${\text{H}}_{\text{2}}{\text{D}}^{+}+\text{CO}\to {\text{DCO}}^{\text{+}}{+\text{H}}_{\text{2}},$$(2a) $${\text{CH}}_{\text{2}}{\text{D}}^{+}+\text{O}\to {\text{DCO}}^{+}{+\text{H}}_{\text{2}},$$(2b)

or with products of CH₂ D⁺ such as CH₄D⁺ $${\text{CH}}_{\text{4}}{\text{D}}^{+}+\text{CO}\to {\text{DCO}}^{+}{+\text{CH}}_{\text{4}}.$$(3)

We model the CD channel using a simple chemical network in 2D, and regard the WD channel as a constant abundance (X_WD) contribution that occurs at temperatures lower than an effective temperature (T_eff). The DCO⁺ chemical network involving the CD channel can be boiled down to only ten chemical reactions. This system can be solved analytically, as proposed by Murillo et al. (2015). We use their prescription and simplified chemical network of the CD channel, and apply it to HD 163296. The input parameters are the gas density (n(H₂)), the gas temperature (T_gas), the CO gas abundance (X(CO)), and the HD gas abundance (X(HD)). The ionization rate (ζ) is constant throughout the disk and equal to 1.3 × 10⁻¹⁷ s⁻¹. The ortho-to-para ratio of H₂ is considered to be in thermal equilibrium (LTE) and is approximated by the following expression: $$\frac{o}{p}=9\text{exp}\left(-\frac{\text{170\hspace{0.17em}K}}{T}\right).$$(4)

We use a lower limit of 10⁻³ at low temperatures as described in Murillo et al. (2015). They contrasted their results to a full chemical network including gas-grain balance (freeze-out, thermal desorption, and cosmic-ray-induced photodesorption) confirming the general trend found by the simplified network. The advantage of using this simple network over a full chemical calculation is that it allows us to more easily investigate the dependency of the DCO⁺ emission on individual model characteristics. Detailed chemical modeling, including the reactions shown in Eqs. (1b) and (1c), would be needed to further investigate the DCO⁺ radial distribution, but this is beyond the scope of this paper.

2.3 Implementation

We adopt the gas density and dust temperature from the physical model used by Mathews et al. (2013) as inputs for the simple chemical network describe above. This parametric model is an approximation of the model used by Qi et al. (2011), which fits both the SED and their millimeter observations. The density structure is defined by $${\Sigma }_d(R)=\{\begin{array}{ll}{\Sigma }_{\text{C}} {\left(\frac{R}{R_{\text{c}}}\right)}^{-1}\text{exp}\left[-\left(\frac{R}{R_{\text{c}}}\right)\right],\hfill & \text{if }R\ge R_{\text{in}},\hfill \\ 0,\hfill & \text{if }R<R_{\text{in}},\hfill \end{array}$$

where Σ_C is determined by the total disk mass M_disk (0.089 M_⊙), R_C (150 AU) is the characteristic radius, and R_in (0.6 AU) is the inner rim of the disk. The vertical structure is treated as a Gaussian distribution with an angular scale height defined by $$h(R)=h_{\text{C}}{\left(\frac{R}{R_{\text{C}}}\right)}^{\psi },$$

where ψ (0.066) is the flaring power of the disk and h_C is the angular scale height at thecharacteristic radius R_C that can take different values for the gas and the dust distribution (see Appendix A; Mathews et al. 2013). The dust temperatureprofile was computed by the 2D radiative transfer code RADMC (Dullemond & Dominik 2004) and is shown in Fig. 2 along with the gas density profile.

In addition to the gas density and dust temperature structure, CO and HD abundance profiles are requiredas inputs for the simple chemical network. We use a CO abundance profile described by three parameters: an effective dust temperature where CO starts to evaporate and becomes optimal for DCO⁺ formation (T_CO), and two constant abundances for gas-phase CO inside (X_low) and outside (X_high) the freeze-out zone. We set X_low = 0.01X_high for all our models. We assume that the dust temperature (T_dust) equals the gas temperature(T_gas), which is a reasonable assumption in the dense regions that we focus on. We set the HD abundance, with respect to the nuclear density n_H = 2n_{H_2}, constant through the entire disk and equal to the cosmic D/H ratio ~10⁻⁵ (Vidal-Madjar 1991).

Fig. 2 Visualization of the adopted physical model. Left panel: model density profile with white contours at 108 cm−3 and 109 cm−3. Right panel: gas temperature structured clipped at 80 K to enhance the color gradient in the region of interest and white contours at 20, 30, and 80 K.

2.4 Radiative transfer

We used LIME (v1.5), a 3D radiative transfer code in non-LTE (Brinch & Hogerheijde 2010) to compare our DCO⁺ observations to the DCO⁺ abundance model obtained using the prescription described above. LIME can produce line and continuum radiation from a physical model, an abundance distribution, and the rate coefficients for a given molecular transition. We use the rate coefficients from the Leiden Atomic and Molecular Database (Schöier et al. 2005)¹ for DCO⁺. These are the same collision rates as those listed for HCO⁺ (Flower 1999). The Einstein A coefficients taken are from CDMS and JPL. We use a grid of 100 000 points that are created applying a weighted random selection in R using a logarithmic scale. This weighted random selection favors denser gas regions and will always select a random point where the DCO⁺ abundance of the model is higher than 10⁻¹⁵. Establishing a convergence criterion that encompasses all of the grid points is difficult. We manually set the number of iterations to 20 and confirm convergence by comparing consecutive iterations.

The resulting model cube is continuum subtracted using the first channel as a continuum estimator. Then each spectral plane of the model is convolved with the synthesized beam of the DCO⁺ data cube. This is equivalent to simulating the visibilities from a sky model since the uv-space is well sampled (see the result in Appendix A).

3 Results

3.1 Standard model (CD)

We chose a standard CO abundance model with a constant abundance X_high = 5.0 × 10⁻⁵ at temperatures above an evaporation temperature of T_CO = 19 K. This evaporation temperature corresponds to ~150 AU at the midplane in our adopted temperature structure and correlates well with the second emission ring thought to be produced by the CD channel (Salinas et al. 2017). The adopted CO abundance is similar to previous models of CO isotopologues (Qi et al. 2015; Carney et al. 2017). We also include a radial cutoff at 300 AU where the DCO⁺ abundance drops to zero because the emissiondisappears at ~300 AU. We discuss this value further in Sect. 4.1.

The resulting DCO⁺ J = 3–2 emission radial profile for the standard model is shown in Fig. 1. The figure also shows our standard model with T_CO = 20 K keeping X_high at the same value to illustrate the model dependency on this parameter. We note that by changing the CO abundance, the model is capable of compensating the difference in evaporation temperature because these two parameters are degenerate. A higher CO abundance will produce less DCO⁺ and vice versa. The adopted CO abundance of 5 × 10⁻⁵ corresponds to a moderate carbon depletion consistent with other warmer disks such as HD 100546 (Kama et al. 2016).

We can obtain an estimate of the emission contribution from the WD channel by subtracting the standard model from the data. We perform a channel-by-channel subtraction from the data cube and the convolved model cube to calculate a residual cube (see example in Appendix A). Figure 3 shows the residual radial curve from the standard models with T_CO = 19 K and T_CO = 20 K and their correspondent abundance estimate. The residual radial profile is obtained applying a Keplerian mask (see Appendix in Salinas et al. 2017) to the residual cube, and each of the radial bins corresponds to a concentric ellipsoid projected to be equidistant to the central star. At R ≳ 180 AU, we can only provide an upper limit of about a few 10⁻¹². At R ≲ 180 AU, the abundance is a few 10⁻¹³ and starts declining at ~50 AU.

The radial abundance estimate is calculated assuming LTE and that the emission is optically thin. If we regard the emission as coming from an isothermal medium in LTE, we can calculate the corresponding column density at different radii assuming an excitation temperature by using the analytical formula of Remijan et al. (2003). This 1D analysis waspreviously performed by Salinas et al. (2017), and we use the same prescription and excitation temperature profile for DCO⁺. Finally, we divide the column density estimate by the surface density profile in Eq. 2.3 to get a vertically averaged abundance. The radial abundance profile only provides a lower limit to the actual DCO⁺ because it emits from a layer set by the activation temperature of the CD and WD channel and the CO evaporation temperature constraining its vertical extent. The simple chemical model of the CD channel can reproduce the DCO⁺ emission in the disk around HD 163296 at large radii (R > 180 AU), but requires additional DCO⁺, of a few 10⁻¹² in abundance to account for the inner emission produced by the WD channel and a radial cutoff at ~300 AU.

Fig. 3 Top panel: residual radial profiles from the models in Fig. 1 with an evaporation temperature of CO TCO = 19–20 K and a constant abundance Xin = 5.0 × 10−5. Bottom panel: vertically averaged abundance estimate of the radial curves of the top panel using the same excitation temperatureprofile and prescription proposed by Salinas et al. (2017) to convert emission to column densities. The 50 AU bar corresponds to the semiminor axis of the synthesized beam and serves as a measure of the spatial resolution.

3.2 Thermal inversion model (CD+TI)

The standard model uses a radial cutoff to suppress the emission of DCO⁺ in the outer disk, but a thermal inversion (TI) could prevent the CD channel from being active (Cleeves 2016). Temperatures higher than the activation barrier for the CD (and WD) channel would prevent the reactions shown in Eqs. (1a) and (1b), reducing DCO⁺ formation. A thermal inversion is produced by the dust evolution modeling of Facchini et al. (2017) as a direct consequence of grain growth and settling. Considering both processes results in a radial decrease in the scale height of the dust because the turbulence is less capable of stirring up the grains as the density drops. This leads to a radial temperature drop as less stellar light is intercepted. At larger radii, where almost all ofthe big grains have migrated inward, small dust is stirred to the disk surface out of the shadow cast in intermediate radii, intercepting more radiation and raising the temperature.

We modify the temperature structure of our standard model to mimic this effect using the following parametrization $$T\prime (R)=25\text{\hspace{0.17em}K}\left(1-\text{exp}{\left(-\frac{R-R_1}{R_2-R_1}\right)}^2\right)+T(R),$$(5)

where R₁ = 240 AU is the inflection point at which theTI occurs and R₂ = 300 AU a characteristic radius. We choose these values so that the temperature at ~290 AU reaches >30 K and effectively blocks the CD channel. We also use a CO abundance parameter of X_high = 2.0 × 10⁻⁷, keeping T_CO at 19 K to better match the shape of DCO⁺ at larger radii. This value is much lower than the 5 × 10⁻⁵ of the fiducial model.

The adopted simple chemical network reaches a maximum in DCO⁺ production for CO abundances of 1 × 10⁻⁵. Higher values result in lower DCO⁺ abundances because the gas-phase CO abundance surpasses the assumed HD abundance blocking its reaction with H${}_3^{+}$ (Eq. (1a)). Similarly, values lower than 1 × 10⁻⁵ result in lower DCO⁺ abundances because CO in the gas phase is required for Eq. (2a) to proceed. If the gradient of the CO abundance is steep, and if the resolution is lower than or equal to ours, we can think of the X_high parameter of our simple step abundance model as an effective CO abundance that accounts for the production gradient of DCO⁺. The thermal inversion parametrization region in our model has a steep temperature gradient that results in a steep second desorption front of CO in the midplane, at ~240–300 AU, comparable in size to our resolution. Our preferred value of X_high = 2.0 × 10⁻⁷ is an effective CO abundance that reproduces the DCO⁺ emission for such a steep gradient.

The first CO desorption front occurs much farther in. Our CD+TI model places this desorption front at about ~150 AU, corresponding to a temperature of ~19 K in our adopted model, where the CO starts to evaporate with modest abundances of ~2 × 10⁻⁷. Observations of C¹⁸O and N₂ H⁺ place the CO snowline at ~90 AU (Qi et al. 2015), corresponding to a temperature of ~22 K in our adopted model, where the bulk of CO is released into the gas phase with abundances of ~5 × 10⁻⁵. Coincidentally, CO abundances of ~5 × 10⁻⁵ and ~2 × 10⁻⁷ produce the same amount of DCO⁺ at temperatures where the CD channel operates. This explains why in Sect. 3.1 the CD model reproduces the DCO⁺ emission with a CO abundance of ~5 × 10⁻⁵ and a desorption front at ~150 AU. We find that to reconcile the CO snowline location at 90 AU, with an abundance of 5 × 10⁻⁵, and the observed DCO⁺ emission our CO abundance profile requires a zone between 90 and 150 AU with an effective abundance of 2 × 10⁻⁷ which we interpret as the onset of the thermal release of CO into the gas phase.

The resulting DCO⁺ abundance, gas density, and emission profile are shown in the middle panels of Fig. 4. DCO⁺ is confined between 30 and 19 K, corresponding to the temperature where the CD channel starts to be efficient and to the temperature where CO freeze-out passes 99%, respectively. The DCO⁺ abundance increases radially as a consequence of the radial increase of the ionization fraction which is $\propto \frac{1}{\sqrt{n}}$ if the ionization rate ζ is constant. The thermal inversion parametrization of this model, coupled with the simple chemical network for CD channel, can reproduce the radial shape of the DCO⁺ in the disk around HD 163296 at R > 180 AU and the lack of emission at R > 300 AU.

Fig. 4 Top panels: resulting radial profile of the DCO+ integrated intensity map of both models and data. The error bars on the data correspond to 3σ, where σ is calculated as the standard deviation of a single ellipsoid divided by the square root of the number of beams.Middle panels: models’ temperature profiles with contours at 19 and 30 K. Bottom panels: DCO+ abundance profiles over-plotted with temperature contours at 19 and 30 K. The hatched region in the DCO+ abundance from model CD marks the radial cutoff used in our standard model to reproduce the absence of emission at R > 300 AU.

3.3 The CD+WD+TI model

Finally, we add a WD channel component to our CD+TI model by modifying the obtained DCO⁺ abundance in the following parametrization: $$X_{{\text{DCO}}^{+}}(r,t)=\left\{\begin{array}{ll}{X}_{\text{CD+TI}}(r,t)+X_{\text{WD}},\hfill & \text{if }T(r,t)\le T_{\text{eff}}\hfill \\ X_{\text{CD+TI}}(r,t),\hfill & \text{if }T(r,t)>T_{\text{eff}}\hfill \end{array}\right\},$$(6)

where X_WD is a constant abundance parameter describing the WD channel contribution and T_eff the effective temperature for the WD channel at which the WD channel switches on. The T_eff is not the temperature where WD starts to operate (which is 80 K), but rather where it reaches its full effectiveness. Just like the CD channel, the WD channel does not switch on abruptly, but gradually increases toward lower temperatures. The right column of Fig. 4 shows the modified DCO⁺ abundance, temperature, and emission profile of the CD+WD+TI model using X_WD = 3.2 × 10⁻¹² and T_eff = 32 K and the parameters from the CD+TI model. These values are taken to match the inner DCO⁺ emission at R ≲ 150 AU. We note that the parameter T_eff = 32 K traces the location of an effective isothermal surface where WD starts to operate. This location is effectively constrained from the data and corresponds to ~40 AU in the midplane. The exact value of the parameter T_eff cannot be constrained without a gas temperature model and only represents the location of the isothermal surface. On the other hand, the X_WD parameter is much better constrained by our modeling because it does not depend strongly on the assumed gas temperature. The CD+WD+TI model reproduces simultaneously the inner and outer features of the DCO⁺ radial profile. The residual (and model) channel maps can be seen in Appendix A.

For comparison, the left column of Fig. 4 shows the DCO⁺ abundance, temperature, and emission profile of model CD. The emission profile of model CD does not include the radial cutoff of our standard model. This radial cutoff is shown as a hatched region in the DCO⁺ abundance profile of model CD. Without reducing the amount of DCO⁺ at larger radii, the CD channel alone overproduces the observed emission. The CD+WD+TI model effectively reproduces the DCO⁺ radial emission profile, at large and small radii, of the disk around HD 163296 including the lack of emission at R > 300 AU.

4 Discussion

4.1 DCO⁺ outer radius

Our standard model uses a radial cutoff in the DCO⁺ abundance to reproduce the absence of DCO⁺ emission at R ≳ 300 AU. In this section we discuss several possibilities to explain such a drop: a drop in total gas density, a local increase in the ortho-to-para ratio of H₂, photodesorption or thermal desorption of CO, and a higher electron fraction.

First, the emission drop could be the consequence of a drop in total gas density. Observations of CO isotoplogues at very high angular resolution (Isella et al. 2016) show that the ¹²CO, ¹³CO, and C¹⁸O emission abruptly diminish at different radii, ~630, ~510, and ~360 AU, respectively. If C¹⁸O (the least abundant CO isotopologue) traces the total gas density, the radial cutoff of DCO⁺ might be due to a lack of total gas density at a similar radius at which the C¹⁸O emission drops. However, isotope selective dissociation of C¹⁸O is a more likely explanation of its drop off at 360 AU (Visser et al. 2009; Miotello et al. 2014). If we adopt a plausible model of the CO column density profile, such as that obtained by Facchini et al. (2017; left panel of their Fig. 8), the column density of C¹⁸O at R ≳ 360 AU is not high enough to self-shield from UV radiation assuming the canonical ISM C¹⁸O-to-¹²CO ratio of 550. This would be consistent with the different radii at which the emission of ¹²CO, ¹³CO, and C¹⁸O drop. It is therefore unlikely that the drop off of DCO⁺ is due to an overall drop in gas surface density. Nevertheless, planet–disk interactions, as probed by the high resolution observations of HD 163296 in the millimeter continuum and CO isotopologues (Isella et al. 2016) can produce features in the dust and gas density profile with spatial scales that are similar to those seen in the DCO⁺ ringlike structures.

Second, the difference in the zero-point energies of ortho and para H₂. Since o-H₂ has more energy than p-H₂, the endothermic reaction that introduces deuterium via the latter has a lower energy barrier than the former, enhancing the deuterium fractionation more efficiently at cold temperatures for a thermal o∕p ratio of H₂. This means that a local increase in the o∕p ratio of H₂ could result in a drop in DCO⁺ abundance. However, it is unclear why H₂ would have a spin temperature much higher (80 K) than the region where DCO⁺ forms (<30 K).

Third, desorption of CO at large radii, either thermally or through photodesorption, can raise the CO abundance above the value where the formation of DCO⁺ is quenched. The distinctive feature at ~260 AU, where the DCO⁺ emission shows a bump, is a natural consequence of this. Both in the model with a thermal inversion and in the model with increased photodesorption through increased UV penetration, the layer of DCO⁺ folds back to the midplane. Increased excitation and column density lead to a maximum in emission; therefore, the presenceof this bump in our data strongly supports this scenario. Small-scale structures in CO abundance, gas temperature, or H₂ o∕p ratio could also produce the observed emission excess at ~250 AU, but not as naturally as a thermal inversion. An extreme case of this scenario can also explain the DCO⁺ emission rings seen in other disks such as IM Lup (Öberg et al. 2015). If the inflection point of the temperature profile occurs at large heights, at low gas densities the emission could hide below the noise level and come back again at larger radii where the DCO⁺ comes back to the midplane. This would give the illusion of two DCO⁺ rings at low sensitivities.

The shape of our modified temperature profile resembles qualitatively the models of Facchini et al. (2017), but their results are heavily model dependent, in particular on the turbulent parameter α. The thermal inversion effect is most pronounced at lower alpha parameters (10⁻³–10⁻⁴). The temperature inversion in these models is smoother than our proposed parametrization. The CO ice in their models is confined from ~200 to ~400 AU, with α =10⁻⁴, whereas our models confine CO ice from ~150 to ~260 AU. A different temperature structure with a slightly hotter disk could also shape the CO ice region differently. Our constraining CD+WD+TI model applied to the DCO⁺ data reveals the location of the thermal inversion and the temperature structure at large radii.

Instead of thermal desorption, photodesorption of CO by increased UV penetration can yield a similar cutoff to the DCO⁺ emission profile. This has been invoked for the DCO⁺ ring seen in IM Lup, and full chemical models of typical disks around T Tauri stars show that this is a plausible scenario (Öberg et al. 2015). Our modeling cannot distinguish this from thermal desorption. Additional observations are needed, either constraining the temperature (e.g., via multiple transitions of the rotational lines of optically thin species such as H₂ CO (Carney et al. 2017) or through UV tracers and/or modeling of the dust. Both thermal and photodesorption could be occurring, potentially creating even more complex structures.

Finally, a higher electron density at radii R > 300 AU cannot significantly destroy DCO⁺ resulting in the observed absence of emission. Our DCO⁺ models cannot constrain the electron density directly because it is degenerate with gas-phase CO abundance. However, if no thermal inversion is invoked and CO remains in the gas phase above 19 K (model CD), the required ionization rate to completely quench the DCO⁺ production at R > 300 AU is at least 8 orders of magnitude higher than the assumed canonical value of 1.36 × 10⁻¹⁷ s⁻¹. This implies an electron abundance of a few 10⁻⁴, much higherthan expected for the warm molecular layer.

4.2 Low-versus high-temperature deuteration pathways

Our preferred CD+WD+TI model uses a constant abundance X_WD = 3.2 × 10⁻¹² and an effective temperatureof T_eff = 32 K to parametrize the contribution and onset of the WD channel. From Fig. 4 the abundance of model CD is about 1.2 × 10⁻¹¹ at larger radii corresponding to a ratio between the DCO⁺ column densities produced by the WD and the CD channel of 0.2. This ratio agrees well with the detailed chemical modeling of Favre et al. (2015), which is on average about 0.25 outside the CO snowline. Since the amount of DCO⁺ produced by the CD channel depends on the CO abundance, an appropriate CO gradient can maximize the DCO⁺ produced via the CD channel and minimize the contribution of the WD channel to ~10%. This limit is only slightly lower than the 20% found by our CD+WD+TI model, which is therefore a robust number. The innermost drop-off of the DCO⁺ emission is probably caused by the dust becoming optically thick in the inner 50 AU (Isella et al. 2016), and not by the WD channel becoming fully operative. This means that our T_eff does not trace the onset of the WD channel, but only the radius where the dust becomes optically thick. On the other hand, our X_WD parameter is independent from this effect. Provided that the contribution of the WD channel can be constrained, DCO⁺ is a promising tracer of the temperature structure and CO snowline of a protoplanetary disk through their CD channel formation pathway.

5 Summary

In this work, we implemented a simple chemical network for cold deuteration and a parametrized treatment of warm deuteration, and carry out a 2D modeling of the DCO⁺ emission in the disk around HD 163296. The following points summarize the conclusions of this work:

• We found that a simple chemical model of the CD channel using a CO constant abundance of 2.0 × 10⁻⁷ above an effective dust temperature of 19 K where CO starts to evaporate, coupled with thermal inversion at around ~260 AU can reproduce the DCO⁺ emission in the outer regions of the disk surrounding HD 163296.

• In addition, modeling the contribution of the WD channel with constant abundance of 3.2 × 10⁻¹² and an effective temperature of 32 K, describing an isothermal surface corresponding to a midplane radius of ~40 AU, reproduces the DCO⁺ emission in the inner disk where the CD channel is not yet active. The ratio of the amount of DCO⁺ produced by the WD to that produced by the CD channels outside the CO snowline in this model is 0.2, consistent with previous full chemical models of DCO⁺.

• With the CD channel tracing the CO abundance at 2.0 × 10⁻⁷ for radii <150 AU, DCO⁺ is a tracer of the onset of CO evaporation. Full evaporation of CO occurs at radii <90 AU, where temperatures are too high for the CD channel to be efficient. This opens possibilities to probe the binding energy of CO ice and its evaporation process.

• We conclude that the formation and destruction mechanisms of DCO⁺ are very temperature sensitive, both through the efficiency of the CD channel and the CO abundance. With proper treatment of the DCO⁺ production through the WD channel, DCO⁺ can be used as a tracer of the location of the CO snowline and the temperature structure, and specifically its gradient, in the outer disk. Furthermore, since the temperature structure and CO abundance both depend on the dust size distribution and spatial distribution, DCO⁺ is a powerful probe of the dust evolution.

Acknowledgements

The authors acknowledge support by Allegro, the European ALMA Regional Center node in The Netherlands, and expert advice from Luke Maud in particular. We also thank Prof. Karin Öberg and Dr. Stefano Facchini for their very useful discussions that helped improve this paper. This work was partially supported by grants from the Netherlands Organization for Scientific Research (NWO) and the Netherlands Research School for Astronomy (NOVA) This paper makes use of the following ALMA data: ADS/JAO.ALMA# 2013.1.01268.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ.

Appendix A Residual channel maps

Fig. A.1 Toppanel: DCO+ J = 3–2 data channel maps. Middle panel: CD+WD+TI model. Bottom panel: residual channels. Contours are 3σ, where σ corresponds to the standard deviation of a line free channel. The channel maps presented here are binned in velocity to enhance S/N.

References

All Figures

Fig. 1 Observed DCO+ radial profile (black continuous line). These profiles are obtained by averaging the values of concentric ellipsoids of the integrated intensity maps. The error bars of the observed DCO+ radial profilecorresponds to 3σ, where σ is the standard deviation of the values contained in the ellipsoid divided by the square root of the number of beams. The 50 AU bar corresponds to the semiminor axis of the synthesized beam and serves as a measure of the spatial resolution. The dashed lines show radial profiles of the emission from models with two different evaporation temperatures of CO TCO = 19–20 K and a constant abundance Xin = 5.0 × 10−5.

In the text

	Fig. 2 Visualization of the adopted physical model. Left panel: model density profile with white contours at 108 cm−3 and 109 cm−3. Right panel: gas temperature structured clipped at 80 K to enhance the color gradient in the region of interest and white contours at 20, 30, and 80 K.
In the text

Fig. 3 Top panel: residual radial profiles from the models in Fig. 1 with an evaporation temperature of CO TCO = 19–20 K and a constant abundance Xin = 5.0 × 10−5. Bottom panel: vertically averaged abundance estimate of the radial curves of the top panel using the same excitation temperatureprofile and prescription proposed by Salinas et al. (2017) to convert emission to column densities. The 50 AU bar corresponds to the semiminor axis of the synthesized beam and serves as a measure of the spatial resolution.

In the text

Fig. 4 Top panels: resulting radial profile of the DCO+ integrated intensity map of both models and data. The error bars on the data correspond to 3σ, where σ is calculated as the standard deviation of a single ellipsoid divided by the square root of the number of beams.Middle panels: models’ temperature profiles with contours at 19 and 30 K. Bottom panels: DCO+ abundance profiles over-plotted with temperature contours at 19 and 30 K. The hatched region in the DCO+ abundance from model CD marks the radial cutoff used in our standard model to reproduce the absence of emission at R > 300 AU.

In the text

	Fig. A.1 Toppanel: DCO+ J = 3–2 data channel maps. Middle panel: CD+WD+TI model. Bottom panel: residual channels. Contours are 3σ, where σ corresponds to the standard deviation of a line free channel. The channel maps presented here are binned in velocity to enhance S/N.
In the text