Multiple paths of deuterium fractionation in protoplanetary disks
aa r X i v : . [ a s t r o - ph . S R ] M a r MULTIPLE PATHS OF DEUTERIUM FRACTIONATION INPROTOPLANETARY DISKS
Yuri Aikawa
Department of Astronomy, The University of Tokyo [email protected]
Kenji Furuya
Center for Computational Sciences, University of Tsukuba, Japan
Ugo Hincelin
Department of Chemistry, The University of Virginia, USA andEric Herbst
Department of Chemistry, The University of Virginia, USA
ABSTRACT
We investigate deuterium chemistry coupled with the nuclear spin-state chem-istry of H and H +3 in protoplanetary disks. Multiple paths of deuterium frac-tionation are found; exchange reactions with D atoms, such as HCO + + D, areeffective in addition to those with HD. In a disk model with grain sizes appropri-ate for dark clouds, the freeze-out of molecules is severe in the outer midplane,while the disk surface is shielded from UV radiation. Gaseous molecules, includ-ing DCO + , thus become abundant at the disk surface, which tends to make theircolumn density distribution relatively flat. If the dust grains have grown to mil-limeter size, the freeze-out rate of neutral species is reduced, and the abundancesof gaseous molecules, including DCO + and N D + , are enhanced in the cold mid-plane. Turbulent diffusion transports D atoms and radicals at the disk surfaceto the midplane, and stable ice species in the midplane to the disk surface. Theeffects of turbulence on chemistry are thus multifold; while DCO + and N D + abundances increase or decrease depending on the regions, HCN and DCN in thegas and ice are much reduced at the innermost radii, compared with the model 2 –without turbulence. When cosmic rays penetrate the disk, the ortho-to-para ra-tio (OPR) of H is found to be thermal in the disk, except in the cold ( . +3 and H D + , as well as the mainreactions of H D + , DCO + , and N D + to analytically derive their abundances inthe cold midplane. Subject headings: astrochemistry — star-formation — protoplanetary disks
1. INTRODUCTION
In the primordial material in the Solar System, such as comets and meteorites, molecularD/H ratios are often higher than the elemental D/H ratio, which is 2 × − (Geiss & Gloeckler1998). For example, the D/H ratio of water is found to be (1 . − × − in comets (e.g.Mumma & Charnley 2011; Altwegg et al. 2015). Since the zero-point energy of D-bearingspecies is lower than that of the normal isotope by up to a few hundred K, the significantdeuterium enrichment is considered to originate in chemistry at low temperatures.There are two possible sites of such low-temperature chemistry: molecular clouds priorto star formation and outer regions of protoplanetary disks. In molecular clouds, the tem-perature is ∼
10 K and high D/H ratios are found for various molecules. For example,even a doubly deuterated molecule, D CO, has been detected in prestellar cores, and the n (D CO)/ n (H CO) ratio has been derived to be 0 . − . n ( i ) denotes the number density of species i . In such dense cold cores, where HD is the pri-mary reservoir of deuterium, deuterated H +3 is produced via exothermic exchange reactions;e.g., H +3 + HD → H D + + H . (1)The high D/H ratio of H +3 propagates to other molecules via ion-molecule reactions in thegas phase. The atomic D/H ratio is also enhanced through the electron recombination ofH +3 and H D + . The high atomic D/H ratio propagates to icy molecules through hydrogena-tion/deuteration of atoms and molecules on grain surfaces. The ions H D + and D H + areactually detected at the center of the prestellar core L1544 (Caselli et al. 2003; Vastel et al.2006). In such a dense cold region, the D/H ratio of H +3 is further enhanced by the freeze-outof CO, which is otherwise the major reactant of H +3 and its isotopologs.Protoplanetary disks are also partially ionized by X-rays and/or cosmic rays, and thetemperature is .
20 K in the outer midplane region, so that deuterium enrichment canproceed in disks, as well. In fact, DCN, DCO + , and N D + have been detected in sev-eral disks (e.g. van Dishoeck et al. 2003; Qi et al. 2008; Huang & ¨Oberg 2015). While a 3 –stable neutral species such as DCN can originate in interstellar ice, then delivered to anddesorbed in the disk, molecular ions are apparently formed in situ, considering their shortdestruction timescales. van Dishoeck et al. (2003) used the JCMT telescope to derive adisk-averaged n (DCO + )/ n (HCO + ) ratio of 0.035 in TW Hya, providing evidence of ongo-ing deuterium enrichment in protoplanetary disks. Although the deuterium enrichment inthe primordial material in the Solar System can originate in interstellar chemistry, at leastpartially, it has been debated if and to what extent the molecular D/H ratios are modifiedin disks (Aikawa & Herbst 1999; Willacy 2007; Cleeves et al. 2014a). Spatially resolvedobservations of deuterated ions can specify the region of active deuteration in protoplanetarydisks.The emission of DCO + is spatially resolved in the disks around TW Hya (Qi et al.2008) and HD 163296 (Mathews et al. 2013). In both disks, the DCO + emission shows aring structure, which is considered to reflect the production of DCO + by the reactionH D + + CO → DCO + + H . (2)The exchange reaction (1) to form H D + is exothermic only by .
250 K, and thus itsbackward reaction becomes active at T &
30 K. Mathews et al. (2013) found that DCO + emission shows a ring structure of radius r = 110 −
160 au, and argued that the outerradius of the DCO + ring corresponds to the CO snow line. The ion DCO + is expected to beabundant in the region with temperature T ∼ −
21 K; at higher temperatures ( T & (see below) would become abundant enough to destroy the precursormolecules, H +3 and H D + , while at T .
19 K, CO abundance would be too low to makeDCO + abundant (Mathews et al. 2013). The DCN emission in TW Hya disk, on the otherhand, is centrally peaked. This molecule is expected to form primarily via the reaction of Natoms with CHD, which is formed by the dissociative recombination of CH D + . Since theexothermicity of the exchange reactionCH +3 + HD → CH D + + H (3)is higher than that of reaction (1), DCN can be abundant even at T &
30 K (Millar, Bennet, & Herbst1989; ¨Oberg et al. 2012).Recent observations of deuterated species in disks, however, challenge the above sce-nario. Qi et al. (2015) observed DCO + in HD 163296 with a higher spatial resolution, andfound that the inner edge of the DCO + ring is at 40 AU, which is much closer to the centralstar than derived by Mathews et al. (2013). Huang et al. (2017) observed six protoplan-etary disks, spatially resolving DCO + , H CO + , H CN, and DCN in most of them. Whilethe DCO + emission tends to be spatially more extended than the DCN emission, the relativedistributions of DCO + and DCN vary among the disks. While the DCN emission is more 4 –compact than the DCO + emission in HD 163296, the radial intensity profiles of DCO + andDCN are similar in AS 209, for example. Such variations indicate that multiple paths todeuteration occur in the disks.On the theoretical side, deuterium chemistry in protoplanetary disks has long beeninvestigated (e.g. Aikawa & Herbst 1999; Willacy 2007; Cleeves et al. 2014a). Recentprogress in the study of deuterium fractionation includes the evaluation of the state-to-staterate coefficients of the H +3 + H system (Hugo et al. 2009) and the re-evaluation of theexothermicity of reaction (3) (Roueff et al. 2013). Since the ground-state energy of ortho-H (o-H ) is higher than that of para-H (p-H ), the backward reactions of exchange reactionssuch as (1) are less endothermic with o-H than with p-H . The deuterium enrichmentwould thus be less efficient if o-H is abundant. Although several groups have presented thedeuterium chemistry in disks coupled with ortho-para chemistry, they concentrated mostly onmolecular species such as HDO rather than DCO + and DCN, or made an assumption that theortho-para chemistry is locally thermalized (Albertsson et al. 2014; Cleeves et al. 2014a,2016). One of the exceptions is the paper by Teague et al. (2015); these authors observedHCO + and DCO + lines in DM Tau, and calculated molecular evolution in a disk model tobe compared with their observations. While they solved for the spin states of H and H +3 ,they did not seem to adopt the updated exothermicity for reaction (3). Roueff et al. (2013)calculated the exothermicity of reaction (3) for p-H to be 654 K, which is higher than theprevious estimate of 370 K (Smith et al. 1982). Favre et al. (2015) adopted this updatedenergy difference for reaction (3) in their chemical model, and showed that it enhances theabundance of DCO + in the warm surface layers of the disk. Their model, however, does notinclude the ortho-para chemistry, and implicitly assumes that all H is in the para state.The assumption would not be appropriate in the warm layers.In the present work, we investigate the deuteration paths of HCO + , N H + , and HCN inprotoplanetary disks using an updated chemical reaction network with deuterium fractiona-tion and ortho-para chemistry. Instead of constructing best-fit models of observed disks, weinvestigate the chemistry in template disk models. The effects of elemental C/O ratio, grainsize, vertical mixing, and exclusion of cosmic rays are also studied. The rest of the paperis organized as follows. Our chemical reaction network and protoplanetary disk models aredescribed in §
2. Section 3 presents the results of our numerical calculations; i.e., abundanceand column density distributions of HCO + , N H + , HCN and their deuterated isotopologs indisk models. The model results are qualitatively compared with observations in §
4, and ourconclusions are summarized in §
5. In appendices, we also present an analysis of the OPRsof H , H +3 , and H D + , and analytical formulae of abundances of H D + , DCO + and N D + inthe cold midplane. 5 –
2. MODELS
We adopt basically the same disk structure and chemical reaction network as in Aikawa et al.(2015). Although basic descriptions of the model and network are given in this section,more detailed descriptions can be found in Aikawa et al. (2015), Furuya et al. (2013), andAikawa & Nomura (2006).
In the present work, we investigate deuterium chemistry in template disk models, ratherthan constructing a model for a specific object. The radial distribution of the gas columndensity in our models is thus rather arbitrary; it is determined so that the mass accretion rate(10 − M ⊙ /yr) is constant at all radii for a constant viscosity parameter α (see § M ∗ =0.5 M ⊙ , surface temperature T ∗ = 4000 K, and radius R ∗ = 0 . R ⊙ . We refer to the UVand X-ray spectrum of TW Hya, and adopt the luminosities of L UV = 10 erg s − and L X =10 erg s − (Herczeg et al. 2002; Kastner et al. 2002)(see also Nomura & Millar 2005;Nomura et al. 2007). The vertical density distribution is set by hydrostatic equilibrium.The temperatures of gas and dust are obtained by solving two-dimensional radiative transferand the balance between cooling and heating. The temperature and density structures aresolved self-consistently. The gas and dust temperatures are the same in the midplane, whilegas is warmer than dust at the disk surface, which is basically a dense photon-dominatedregion.In our fiducial model, the grains have grown to the maximum size a max of 1 mm size,while the minimum size is a min = 0 . µ m, with the dust-to-gass mass ratio of 0.01, all overthe disk. A grain size distribution is assumed to follow the power law dn ( a ) /da ∝ a − . ,where n ( a ) is the number density of grains with size a . For comparison, we also adopt a diskmodel with dust appropriate to dark clouds, which we label dark cloud dust. The maximumgrain size of the latter model is 10 µ m (Weingartner & Draine 2001). The density andtemperature distributions of the two disk models are shown in Figure 1. The temperatureis generally higher in the model with dark cloud dust, because it absorbs stellar radiationmore efficiently. While we assume that the disk is static (i.e. no diffusion or accretion), wealso calculate models with vertical mixing (Furuya et al. 2013), in order to investigate itseffect on deuterium chemistry, 6 – The gas-grain chemical reaction network is based on Garrod & Herbst (2006), but hasbeen updated to include reactions that are effective at high temperatures ( T &
100 K)(Harada et al. 2010, 2012). Our network includes up to triply deuterated species, and nu-clear spin states of H , H +3 , and their isotopologs (Hugo et al. 2009; Hincelin et al. 2014;Coutens et al. 2014; Furuya et al. 2015). We refer to Roueff et al. (2013) for the reac-tion rate coefficients and exothermicities/endothermcities of the exchange reaction (3) andanalogous reactions with multi-deuteration and spin states.Ultra violet radiation from the central star and interstellar radiation field causes pho-todissociation and photoionization in the gas phase; the rate coefficients are calculated byconvolving the dissociation/ionization cross sections of molecules and the UV spectrum ateach position in the disk. Self- and mutual shielding of H , HD, CO, N , and C atoms aretaken into account (Draine & Bertoldi 1996; Visser et al. 2009; Wolcott-Green & Haiman2011; Li et al. 2013; Kamp & Bertoldi 2000). The shielding of D is not considered in thepresent model, but we confirmed that it does not affect our results, because D is muchless abundant than HD at the disk surface, where photodissociation is effective. Photo-dissociation of icy molecules is also considered in our model (e.g. Furuya et al. 2013). Theionization sources in the disk are X-rays, cosmic rays, and the decay of radioactive nuclei.We assume a cosmic-ray ionization rate of 5 × − s − with an attenuation length of 96g (Umebayashi & Nakano 1981). Considering that the penetration of cosmic rays can beprohibited by stellar winds, we also run a model without cosmic rays (Aikawa et al. 1999;Cleeves et al. 20014b); the disk midplane is then mainly ionized by the decay of radioactivenuclei with a rate of 10 − s − . Layers above the midplane are mainly ionized by stellarX-rays; the ionization rate could reach 10 − s − at the disk surface, for example (see Figure1 in Aikawa et al. 2015).We adopt a two-phase model, which consists of a gas phase and an undifferentiatedgrain ice mantle. The sticking probability is assumed to be unity when a neutral atom ormolecule collides with a grain surface, except for H and D atoms, for which we adopt thetemperature-dependent sticking probability of Tielens & Hollenbach (1985). Adsorptionenergies of molecules on gain surfaces are generally taken from Garrod & Herbst (2006).The adsorption energies of deuterated species are generally set to be the same as those ofnormal species. One exception is the D atom, the adsorption energy of which is set to be 21K higher than that of the H atom, which is assumed to be 600 K (Caselli et al. 2002). Theadsorption energies of CO, N and HCN, the abundance of which we present in the followingsection, are set to be 1150 K, 1000 K, and 3370 K, respectively. The value for HCN isadopted from the Temperature Programmed Desorption (TPD) experiment by Noble et al. 7 –(2013). The dependence of molecular abundances on the adsorption energies of CO and N are presented in Aikawa et al. (2015).In addition to thermal desorption, we take into account three non-thermal desorptionprocesses: photodesorption, stochastic heating by cosmic rays, and reactive desorption (e.g.¨Oberg et al. 2009a; Hasegawa & Herbst 1993; Garrod et al. 2007). The efficiency of chem-ical desorption is set to be 10 − − − depending of the reactions, referring to Garrod et al.(2007). Recent laboratory experiments show that the efficiency depends on the molecularcomposition of the grain surface (Minissale et al. 2016). In order to include such an effect,we need to adopt the three-phase model, discriminating the ice surface from the bulk icemantle, which we postpone to future work. We assume that the surface reactions occur viathe Langmuir-Hinshelwood mechanism; i.e. the adsorbed species diffuse on grain surfacesvia thermal hopping, and react with each other when they meet, before they desorb. Themodified rate of Caselli et al. (1998) is adopted to ensure that the rate of hydrogenation(deuteration) is not higher than the adsorption rate of H atoms (D atoms) when the numberof such atoms on a grain is small ( . E diff ) is set toone-half of the adsorption energy ( E ads ). Lower values from 0.3 to 0.4 have also been used inastrochemical models (e.g. Ruaud et al. 2016; Penteado et al. 2017); the surface reactionsbecome more rapid with the lower ratio of E diff /E ads . The efficiency of H formation andthose of its isotopologs on granular surfaces become lower in warmer regions (e.g. T & a few10 K) in the disk, since the physisorbed H (and D) atoms are desorbed more promptly. Whenthe temperature is high enough, however, a bare grain surface appears, on which H atomscan be chemisorbed. At such high temperatures, we assume that the formation rate of H is 0.2 times the sticking rate of H atom onto grain surfaces, referring to Cazaux & Tielens(2004, 2010) (see eq. 19 of Furuya et al. 2013).The elemental abundance of deuterium is set to be 1 . × − relative to hydrogen(Linsky 2003). The initial molecular abundances are determined by considering the molecu-lar evolution from a cloud formation stage to a collapse phase to form a protostar. Furuya et al.(2015) first calculated molecular evolution behind the shock front of colliding HI gas. Furuya et al.(2016) then calculated the subsequent molecular evolution in the hydrostatic prestellar coreand in its collapse phase to form a protostar, as in Aikawa et al. (2012). We adopt thefiducial model of Furuya et al. (2016) for our initial abundances; the shock model is rununtil the column density of the post shock gas reaches a visual extinction of 2 mag, andthe duration of the hydrostatic prestellar phase is set to be 10 yr. Then the core collapsesto form a protostar in 2 . × yr, and the molecular evolution in the infalling envelope iscalculated until the protostar is 9 . × years old. While Furuya et al. (2016) calculateda spatial distribution of molecular abundances in a protostellar core, our initial abundancesfor the present work are adopted from the infalling fluid parcel which reaches a radius of 60 8 –au in the core. The initial abundances of major species are presented in Table 1; since thetemperature ( ∼
144 K) and density ( ∼ × cm − ) are high in the infalling fluid parcelat 60 au, ices and ions are not abundant. While the OPR of H is initially 3 . × − , itreaches thermal equilibrium or steady state in a relatively short timescale (see Appendix A).Our results thus do not significantly depend on the initial OPR of H .
3. RESULTS
In the following, we consider the spatial distributions of HCO + , DCO + , N H + , N D + ,HCN, and DCN, and describe the major formation pathways of the deuterated species.We present the results at 3 × yr, which is earlier than a typical age of T Tauri stars( ∼ + and N H + (and their isotopologs) are strongly coupled with COand N , which are gradually converted to less volatile species such as CO and NH , anddepleted from the gas phase. The timescale of this conversion in the gas phase is, however, ∼ × (cid:0) ζ He . × − s − (cid:1) − yr, where ζ He is the ionization rate of He atom. It is comparable to orlonger than the accretion timescale of the disk at r .
100 au (Furuya & Aikawa 2014)(seealso Furuya et al. 2013). In other words, gaseous CO and N can be supplied from theouter radius by viscous accretion and the radial drift of dust grains with ice mantles. Theabundances of CO and N are thus underestimated, if we choose 1 Myr in our static models.The temporal variation of the abundance of each molecular species is also described in thefollowing. + and DCO + Let us first consider HCO + and DCO + . The former is mainly formed by the reactionCO + H +3 . As a reference, distributions of gaseous CO abundance and its column densityare shown in Figure 2 (a) (b), while those of HCO + and DCO + are shown in Figure 2, panels(c)(d) and (f). The horizontal axis depicts the radius ( r ) of the disk. In panels (a) (c) (d)and (e), the vertical axis is the distance from the midplane z , normalized by the radius;i.e. z/r , while the vertical axis is the column density for panels (b) and (f). In the panelsshowing the abundance distributions (a, c and d), the dotted lines depict the positions wherethe X-ray ionization rate is equal to the cosmic-ray ionization rate (5 × − s − ) (the upperdotted line) and to the ionization rate by decay of radio active nuclei (1 × − s − ) (the 9 –lower dotted line). The long dashed line indicates the CO snow surface; it indicates thepositions where CO abundances in the gas phase and in the ice mantle should be equal, ifthe abundances are determined simply by adsorption onto grains and thermal desorption.In the layer below the long dashed line, CO is expected to be mainly in the ice mantle, whileit is expected to be abundant in the gas phase, elsewhere. The snow line is defined as theradius at which the snow surface crosses the disk midplane.In our model, CO is gradually converted to less volatile molecules such as CO ice andCH OH ice and thus depleted even in the regions inside the CO snow line ( ∼
20 au) andin the layer above the CO snow surface (Furuya & Aikawa 2014). Although we presentthe results at 3 × yr, which is shorter than the timescale of chemical conversion in thegas phase, chemical conversion on grain surfaces is more efficient than in the gas phase. Inthe midplane at 40 au . r .
100 au, CO (both in the gas phase and in ice mantle), andthus HCO + , are deficient ( n ( i ) /n (H) . − ), and such regions with low CO and HCO + abundances extend inwards over time. For example, the gaseous CO abundance is as low as5 × − and the HCO + abundance is 1 × − in the midplane at r = 10 au at t = 1 Myr.It should be noted that ALMA observations recently revealed such depletion of gaseous COinside the CO snow line in TW Hya (Nomura et al. 2016; Schwarz et al. 2016) (see alsoZhang et al. 2017). In the midplane at inner radii ( r .
10 au), the HCO + abundance islimited by the low ionization degree, which occurs at high density. At r .
70 au, molecularions such as H +3 and HCO + are deficient at z/r ∼ .
15, because photoionization makesatomic ions (e.g. S + ) the dominant positive charge carrier in such upper layers. The highelectron abundance then reacts quickly with molecular ions.Figure 2 (e) depicts the major formation pathways of DCO + at each position in thedisk. Comparing panels (d) and (e), we can see that DCO + is abundant in the layer belowthe CO snow surface (i.e. at the lower z ) due to reaction (2), despite the relatively low COabundance in the gas phase. In the midplane, deuterated-H +3 increases towards the outerradius, where the temperature is lower and the chemical conversion of HD is less efficient(see Appendix B). The chemical conversion of CO is also less efficient at outer radii ( & + abundance thus increases outwards at r &
100 AU. In the upper layers( z/r & . + is mainly formed viaHCO + + D → DCO + + H + 796K . (4)Since the reaction is exothermic by 796 K (Adams & Smith 1985), the backward reactionis less efficient than the forward reaction even at warm temperatures.These results are basically the same as described in ¨Oberg et al. (2015). One update isthat, following Roueff et al. (2013) we set the exothermicities of reaction (3) and its multi-deuterated counterparts to values higher than those previously determined (Smith et al. 10 –1982) and used in ¨Oberg et al. (2015). For example, the exothermicity of CH +3 + p-H andCH +3 + o-H are set to 660 K and 489 K, respectively. Although we do not distinguish thenuclear spin state of CH +3 , this omission does not affect our results, since the difference inexothermicities among the reactions with ortho- and para-CH +3 is only 6 K (Roueff et al.2013). Although DCO + is also formed via the reaction between CO and CH D + , wherethe ion is produced via reaction (3) and subsequent radiative association with H , thispath is the major formation path of DCO + only in limited regions, i.e. the green regionsin panel (e). It is in contrast to the model of Favre et al. (2015), in which reaction (3)significantly contributes to the DCO + formation in warm surface layers and thus enhancesthe DCO + column density. Although the direct comparison between our model results andthose of Favre et al. (2015) is not straightforward due to the differences in disk physicalstructure, knowledge of the spin state of H would be a key to differentiate between the tworesults. While Favre et al. (2015) assumed all H in the para form, and used a constantexothermicity of 654 K for reaction (3), we found that the OPR of H is almost thermal ateach position in the disk (see Appendix A); the effective exothermicity of reaction (3) thusdecreases as o-H increases with temperature. Figure 3 shows the ratio of the backwardto forward reaction rate coefficients referring to Roueff et al. (2013). For the solid line, weassumed that the OPR of H is thermal. The dashed line, on the other hand, shows the ratiowhen all H is in para state, so that the exothermicity of the forward reaction is constantat 654 K. The dotted line depicts the ratio of the former (i.e. the ratio shown with thesolid line) to the latter (the dashed line); it reaches a maximum at the temperature of ∼ + claimed by Favre et al.(2015). It should be noted, however, that the contributions of reaction (3) become moresignificant when the elemental C/O ratio is higher than unity ( § + (dashed lines)and DCO + (solid lines). The blue, green, and red lines depict the values at t = 1 × yr, 3 × yr (the fiducial value), and 9 . × yr, respectively. Both HCO + and DCO + increase inwards from r ∼
40 au, which is slightly outside the CO snow line ( ∼
20 AU).In the midplane inside r ∼
40 au, thermal desorption of CO becomes non-negligible, whichmakes HCO + the dominant charge carrier, while H +3 and its isotoplogues dominate in theouter radius (see eq. (23) in Aikawa et al. 2015). Inside a radius of ∼ + abundance is low in the midplane; the warm temperature lowers the atomic D/H ratio, andthe D atom abundance is also lowered by reaction with HS. Outside 50 au, the columndensity of DCO + increases outwards due to its enhancement in the outer midplane, whilethe HCO + column density is relatively flat. Both column densities gradually decrease withtime, as CO is converted to less volatile species. The decline of the DCO + column densityis more significant than that of HCO + , since DCO + has its abundance peak in the midplane 11 –at outer radii, where CO conversion is efficient. H + and N D + Figure 4 shows distributions of molecular abundances, column densities, and majorformation pathways of deuterium isotopologs as in Figure 2, but for N , N H + , and N D + .N H + is formed by N + H +3 , and is destroyed by recombination with electrons and protontransfer to CO. Thus its abundance depends on N , CO, H +3 and electrons (Aikawa et al.2015); the abundance has a peak around the CO snow surface, where the abundance ratioof CO to electron is ∼ , and in the layer above the CO snow surface, where gaseous COis converted to less volatile species to be frozen onto grains.In the midplane, N is gradually converted to less volatile species such as NH ice,and thus the region of low N H + abundance basically expands with time, although N H + becomes abundant at t = 9 . × yr around the radius of several au, where gaseous COis reduced by the conversion effect. The formation pathways of N D + are similar to thoseof DCO + ; N D + in the midplane is formed by the reaction of N with deuterated H +3 , whileN D + is formed by N H + + D → N D + + H + 550K (5)above the CO snow surface. In the upper layers, the ratio of N D + /N H + is lower than thatof DCO + /HCO + , because the exothermicity of reaction (5) (550 K) is lower than that ofreaction (4). The column densities of both N H + and N D + have a peak at 20 au . r .
50 au. The inner boundary corresponds to the snow line of CO, the main reactant withN H + , while the outer boundary is slightly outside the N snow line, for a similar reasonused for HCO + . Outside r ∼
150 au, the N D + column density increases outwards, since itis abundant in the outer midplane. Figure 5 shows the distributions of HCN and DCN, the formation pathways of DCN,and the column densities of HCN and DCN. Since the desorption energy of HCN is high(3370 K), the abundance of gaseous HCN is . − outside a radius of a few au, wherethe temperature is .
100 K. HCN in the ice mantle, on the other hand, is as abundant as10 − − − in the midplane. Towards the outer ( r &
30 au) midplane, HCN forms by thedissociative recombination of H CN + , which forms by H +3 + CN → HCN + + H followedby HCN + + H → H CN + + H. In the upper layers and inner radii, it forms by H + H CN 12 –and N + HCO. The precursor molecules H CN and HCO are produced by CH + N andCH + O, respectively. HCN also is produced by reactive desorption after the grain surfacereaction of H + CN.The major formation paths of DCN are mostly the deuterated version of the abovereactions. One exception is HCN + D → DCN + H , (6)which is effective at the disk surface, depicted by the green region in Figure 5 (c). In ourmodel, the activation barrier of both the forward and backward reaction of (6) is set to be 500K, following Schilke et al. (1992). A more recent quantum chemical calculation (Kayanumain private communication) evaluates this barrier to be 3271 K, but also predicts that theeffect of the barrier would be significantly lowered by tunneling. Our result does not change,even if we assume a higher barrier, e.g. 1800 K, for reaction (6). Another exception isC H D + + N, which dominates in the midplane at r ∼ < T <
50 K), hydrocarbons are also deuterated viaC H +2 + HD → C HD + + H , (7)which is exothermic by 550 K. The major exothermic exchange reaction that initiates theenhancement of the DCN/HCN ratio is thus reaction (1) in the outer midplane, while reac-tions (3) and (7) dominate in the upper layers and in inner radii. It should be noted thatHCND + and DCNH + are not distinguished in our reaction network, which means that thenetwork of DCN and DNC is partially mixed in regions where the recombination of HDCN + is their major formation path.Due to the contribution from the outer midplane, the column density of DCN increasesoutwards at r &
20 au. Inside a radius of a few au, both the HCN and DCN column densitiesare high, since they are relatively abundant in the z/r ∼ .
06 and z/r ∼ C N serve as reservoirs of carbonand nitrogen, while the reactants of C-bearing species, such as O atoms, are deficient in themidplane with high density.
Observations in recent years suggest that the surface and molecular layers are deficientin elemental carbon and oxygen, especially in relatively cold disks. Hogerheijde et al. (2011)detected ground-state rotational emission lines of H O in the disk around TW Hya usingthe Herschel Space Observatory to find that the emissions are significantly weaker thanexpected from disk models. Du et al. (2017) obtained and analyzed the water emission of 13 –13 protoplanetary disks. They compared the observational data with disk models to showthat the abundance of gas-phase oxygen needs to be reduced by a factor of at least ∼ O taken by
Spitzer to find that water vapor is significantly depletedin the disk surface beyond the radius of ∼ − − − ) even in the warm molecular layer (above the CO snow surface) in the diskof TW Hya by comparing C O ( J = 2 −
1) and HD ( J = 1 −
0) emission line intensities.Kama et al. (2016) combined various emission lines including CI and OI towards TW Hyaand HD100546; carbon and oxygen were found to be strongly depleted from the gas phase inthe disk of TW Hya, while the depletion is moderate in the disk of HD100546. Kama et al.(2016) also presented an analytical model of the vertical cold finger effect to show that acombination of turbulent mixing and settling of large dust grains deplete C- and O-bearingvolatiles from the surface and molecular layers of disks by locking them in the ices in themidplane (see also Krijt et al. 2016; Xu et al. 2017). Since H O is less volatile than CO,such a depletion mechanism would be more efficient for oxygen than carbon, which couldenhance the C/O ratio in the surface and molecular layers. Kama et al. (2016) found theC/O ratio to be higher than unity in the disk of TW Hya. Bergin et al. (2016) observed C Hemission, which is bright in a ring region, in TW Hya and DM Tau. In order to reproducethe bright C H emission in disk models, a high C/O ratio ( >
1) is required together with astrong UV field.In order to investigate the depletion of elemental carbon and oxygen on deuteriumchemistry, we performed a calculation of our fiducial disk model as in § O to zero and reducethe CO abundance by an order of magnitude, while the abundances of other species are thesame as in our fiducial model. Even with this reduced abundance, CO is still the majorcarbon carrier, although CH is the most abundant among C-bearing species (see Table 1).Major oxygen carrier in the initial condition is CO, H CO, and CO . The elemental C/Oratio is 1.43. Figure 6 shows the distributions of CO, HCO + and DCO + abundances (panelsa, c, and d), their column densities (panels b and f) and major formation pathways of DCO + (panel e). While the spatial distributions of the molecular abundances and column densitiesof HCO + and DCO + are basically similar to the fiducial model (Figure 2), there are severalnotable differences. Firstly, CO depletion via chemical conversion to CO ice is less effectiveabove the CO snow surface in the outer radius ( r &
20 au) than in the fiducial model, which 14 –is natural considering the reduced oxygen abundance. HCO + and DCO + thus becomes moreabundant in the layer above the CO snow surface, in which DCO + is mainly formed byreaction (3). The high C/O ratio enhances the abundance of hydrocarbons and thus theimportance of reaction (3). In spite of the reduced CO abundance, the column densities ofHCO + and DCO + are similar to those in the fiducial model, except that the depression ofthe DCO + column density at r ∼
50 au is more modest and the HCO + column density isreduced at the central region ( r . a few au), in which HCO + exists mostly in the surfacelayer.Figure 7 (a-f) shows the distributions of N , N H + , and N D + , their column densities,and the major formation pathways of N D + . The distributions of N and N D + are similarto those in our fiducial model (Figure 4), while N H + is redued in the layer above the COsnow surface, in which the gaseous CO abundance exceeds that in the fiducial model.Figure 7 (g-j) shows the distributions of HCN and DCN, their column densities, and themajor formation pathways of DCN. Their abundances and column densities are significantlyhigher than in the fiducial model, since the high C/O ratio enhances the abundance ofhydrocarbons, which react with N atoms to form HCN.We also calculated a model in which the initial abundance of H O is totally depletedbut CO is not (the C/O ratio then becomes 1.10); the results are quite similar to the modeldescribed above, except that CO is more abundant inside the CO snow line, which lowersthe N H + abundance . Figure 8 shows distributions of gaseous CO, HCO + , and DCO + (panels a, c and d), themain formation paths of DCO + (panel e), and their column densities (panels b and f) in themodel with dark cloud dust. The initial abundances are the same as in our fiducial model.Since the total surface area of grains is larger than in the fiducial model, the freeze-out ofCO and subsequent conversion to other molecules is more efficient in the dark cloud dustmodel than in the mm-sized grain model, especially in the midplane regions. In the layerabove the midplane, on the other hand, UV radiation is more efficiently attenuated, so thatthe molecular layer extends to larger z than in Figure 2.At early time (e.g. t = 1 × yr), HCO + is abundant in the layer above the COsnow surface, and in the midplane inside the CO snow line. Since the contribution from theupper layers is significant, the radial distribution of the HCO + column density is rather flat;although it slightly increases inwards around the CO snow line ( r ∼
100 au), the increment 15 –is not significant compared with that in Figure 2 at r ∼
40 au. HCO + in the midplanedecreases with time, as CO is converted to less volatile species. Its column density, however,varies by less than a factor of three, due to its constantly high abundance in the upper( z/r & .
3) layers.The DCO + ion is mainly formed by reaction (4), since CO is depleted in the midplane,where reaction (2) should be efficient. At the early time, D atoms are abundant in the upper( z/r & .
2) warm layers, where reformation of HD is inefficient. DCO + at that stage is thusabundant in the upper layers, and its column density distribution is relatively flat at r & + decreases as D atoms areincorporated into hydrocarbons, and CO inside the snow line is converted to less volatilesspecies. Its column density thus decreases significantly with time. Even at 9 . × yr,DCO + in the midplane does not significantly contribute to its column density except arounda radius of a few tens of au, where the column density has a sharp peak. This is in contrastto the model with mm dust grains, in which DCO + in the midplane mostly determines itscolumn density distribution.In order to check the significance of the updated exothermicity of reaction (3), we alsoran a model in which its exothermicity is set to 370 K. Compared with the model with thisold value, the DCO + abundance is enhanced at 0 . . z/r . . + there.Figure 9 shows the distributions of N , N H + , HCN and their deuterated isotopologs,their column densities, and the major formation pathways of the deuterated isotopologs inthe model with dark cloud dust. N H + is abundant in the upper layers, where H +3 is thedominant ion, and in the midplane region between the snow lines of CO ( ∼
120 au) andN ( ∼
230 au). Its abundance in the midplane, however, decreases as N is converted toNH ice in a few × yr (Furuya & Aikawa 2014). The distribution of N D + is similar tothat of N H + , but its abundance in the upper ( z/r & .
3) layer is lower than the maximumabundance in the midplane. The column density of N D + , therefore, decreases significantlyas N is converted to NH ice in the midplane. The major deuteration path is via H D + inthe midplane, and via D atoms (reaction 5) in the upper layers.Outside a radius of r ∼
10 au, gaseous HCN and DCN are distributed mostly in theupper ( z/r & .
2) layers, where they form via the recombination of H CN + (HDCN + ) andreaction H CN (HDCN) + H. The radial distributions of their column densities graduallyincrease outwards at r &
10 au. At the innermost radii ( r . &
100 K) to desorb HCN and DCN from ice mantles. 16 –We also calculated molecular abundances in the disk model with the maximum grain(pebble) size of 10 cm, and found that the distributions of gaseous molecular abundancesare qualitatively similar to those in our fiducial model (i.e. the model with mm-sized dust).Some notable differences are as follows. Firstly, the midplane temperature is slightly higherin the model with cm dust. Secondly, CO is converted to CO ice more efficiently in themodel with cm dust inside the CO snow line (around a radius of a few tens of au), becausethe OH radical, which reacts with CO to produce CO , is produced from H O ice due tothe higher UV flux. Thirdly, gaseous HCN and DCN are more abundant in the model withcm dust outside a radius of a few au; a deeper penetration of UV and slightly warmertemperature enhance their precursors, N atoms and hydrocarbons, in the gas phase. Insidetheir snow line ( . Protoplanetary disks are considered to be turbulent, most probably due to magneto-rotational instability (Balbus & Hawley 1991). The direct measurement of non-thermalvelocity dispersion v has been one of the major challenges in radio observations of disks. Thedispersion is basically subsonic with a Mach number M = v/c s ∼ . − . c s is the sound speed. While a very low velocity dispersion M < .
03 isderived in HD163296 (Flaherty et al. 2015), Teague et al. (2016) took into account theuncertainties in flux calibration to derive an upper limit of
M ∼ .
16. Since the diskhas vertical temperature gradients, and since the mixing timescale in the vertical directionis shorter than that in the radial direction (e.g. Aikawa et al. 1996), the vertical mixingcould alter the molecular D/H ratios (Furuya et al. 2013; Albertsson et al. 2014). In thissubsection, we investigate the effect of vertical turbulent mixing on chemistry including theD/H ratios.The diffusion coefficient is of the same order as the kinematic viscosity coefficient, αc s H ,where H is the scale height of the disk. Although there could be a slight difference betweenthe two values (e.g. Johansen & Klahr 2001), we assume that the values of the two coef-ficients are the same in the present work. The non dimensional parameter α is equal to( v/c s )( l/H ) ≈ ( v/c s ) , where l is the size of the turbulent eddy. The value of α is thusestimated to be . − in protoplanetary disks. Figure 10 and Figure 11 show distributionsof neutral and ionic species, respectively, in models with a diffusion coefficient α = 10 − (toprow) and 10 − (middle row). The bottom panels show the molecular column densities for themodel without diffusion (solid lines) and with a diffusion coefficient of α = 10 − (dashed) 17 –and 10 − (dotted). CO and N are also shown in Figure 10 in addition to HCN and DCN,since DCO + and N D + are chemically coupled to these precursor species.The effect of turbulence is twofold. First, it transports ices to the disk surface, wherethe ices are thermally desorbed and photodissociated. Secondly, it transports H atomsand radicals from the disk surface to the midplane, where they contribute to grain-surfacereactions. In the static model, CO is efficiently converted to CO via the grain-surfacereaction of CO + OH in the layer above the CO snow surface ( § ice is transported to the disk surfaceto be photodissociated. Thus the column density of gaseous CO is slightly higher in themodels with diffusion at r &
60 au. In the midplane at 10 au . r .
50 au, on the otherhand, CO is more efficiently converted to CO ( α = 10 − ) and CH OH ( α = 10 − ) in themodels with diffusion, due to the enhanced abundances of OH and H atoms. Inside a radiusof ∼
40 au, HCO + is mostly on the disk surface, and is increased by the diffusion of CO inthe turbulent disk.The major nitrogen carriers are N and NH in our models. In the models with diffusion,the N column density is reduced at r &
50 AU, while it is enhanced at smaller radii,compared with the model without diffusion. In the cold outer midplane region, NH is moreefficiently formed via H atoms coming from the disk surface, while in the inner warm regions,midplane temperatures are too warm to (re)form NH via grain-surface hydrogenation, andNH is transported to the disk surface to be photodissociated (Furuya & Aikawa 2014). Thehigh abundance of N in the midplane enhances the abundance of N H + at 10 au . r . D ices (see Appendix B) less efficientvia the transport of the deuterated ices to the disk surface, where they can be dissociated.Thus in the model with mixing, HD and D atoms are more abundant in the midplane at r &
40 au. Secondly, D atoms are transported towards the midplane and enhance deuterationvia reactions such as (4). 18 –
So far we have assumed that cosmic rays provide a minimum ionization rate of 5 × − s − in the midplane, where X-ray penetration is significantly attenuated. But the cosmicrays may be excluded by stellar winds with magnetic fields (e.g. Umebayashi & Nakano1981; Cleeves et al. 2013a). In order to check the effect of the exclusion of cosmic rays,we investigated the model without cosmic-ray ionization. Here, the model disk is ionized bystellar X-rays, which dominate in and above the molecular layer, and the decay of radioactivenuclei, which sets a minimum ionization rate of 1 × − s − in the midplane (below the lowerdotted line in Figure 12), given the Al abundance derived from the analysis of meteoritesfor the formation stage of the Solar System (Umebayashi & Nakano 2009). Although theionization rate could be smaller depending on the abundances of radioactive nuclei and thesurface density of the disk (Umebayashi & Nakano 2009; Cleeves et al. 2013b), the effect ofa low ionization rate is already apparent in the model presented here, and it is straightforwardto extrapolate the results to an even lower ionization rate, at least qualitatively.Figure 12 shows the distributions and column densities of HCO + , N H + , HCN and theirdeuterated isotopologs at t = 3 × yr. The solid lines depict the column density in themodel without cosmic-ray ionization, while the dashed lines depict our fiducial model. It isnatural that the abundances of ionic molecules are suppressed by the low ionization rate.It should be noted that the dependence of molecular ion column densities on cosmic-rayionization rate would vary among disk models. In the model with dark cloud dust, themidplane abundances of these molecular ions are lower due to the more efficient freeze-out ofC-bearing and N-bearing molecules, and thus the X-ray dominated layer contributes more tothe column density than in the model with mm-sized grains. Thus the decline of molecularion column densities due to the exclusion of cosmic rays would be less significant in the diskswith dark cloud dust. We found that the column density of HCO + in the model withoutcosmic-ray ionization, for example, is similar to that in Figure 8 in the dark cloud dustmodel. The column density of DCO + , on the other hand, is reduced by about an order ofmagnitude at 10 au . r .
100 au compared with Figure 8.Comparing the right column in Figure 12 with that in Figure 5, we can see that bothHCN and DCN are reduced in the outer ( r &
10 au) midplane in the model without cosmic-ray ionization. Cosmic rays play a key role in forming N atoms and hydrocarbons such asCH from N and CO, respectively. N atoms are formed by the recombination of N H + witha small branching ratio, which is produced by protonation of N . CO reacts with He + to formC + , which goes through successive reactions with H and electrons to form hydrocarbons.Inside a radius of ∼ is suppressed by the low ionization rate. In thedisk with dark cloud dust, both HCN and DCN are mostly abundant in the warm X-raydominated layers, and thus their column densities do not depend much on the midplaneionization rate, except at r .
4. Discussion
The motivation of the present work is to investigate the major deuteration paths andtheir efficiency in disk models with an updated gas-grain chemical network, rather thanconstructing a best fit model for a specific object. It is, however, worth comparing ourmodel results with recent observations of deuterated species and their normal isotopologs.Table 2 summarizes the morphologies of integrated intensity maps of H CO + , DCO + , N H + ,N D + , H CN, and DCN in two full disks around T Tauri stars (AS 209 and IM Lup), twotransition disks around T Tauri stars (V4046 and LkCa 15) and two disks around HerbigAe stars (MWC 480 and HD 163296). It is mainly based on Huang et al. (2017), and issupplemented by ¨Oberg et al. (20011), Huang & ¨Oberg (2015), Salinas et al. (2017), andFlaherty et al. (2017). Comparison should be qualitative, rather than quantitative, sincedisk physical structures are expected to vary among objects, including the radial distributionof midplane temperature, which sets the location of the CO snow line. While we comparethe distribution of molecular lines with the estimated position of the CO snow line andtemperature distribution for some disks, the derivation of temperature distributions in disksis not straightforward due to the opacity effect and temperature gradient in the verticaldirection. CO sublimation temperature, and thus the position of CO snow line could alsodepend on the ice composition, i.e. whether its surface is water-rich or not, while desorptionenergy of a molecule is set to be constant in our model. It should also be noted that wecompare the radial profiles (distributions) of molecular column density in our models withthe observed distributions of line emissions. Huang et al. (2017) observed the J = 3 − + , DCN, H CO + , and H CN, while Huang & ¨Oberg (2015) observedthe J = 3 − D + . Using the RADEX code (van der Tak et al. 2007), we estimatethat the opacity of the J = 3 − + , DCO + , N H + and HCN reachesunity, when the molecular column density is a few 10 cm − , with a gas temperature of 30K. These emission lines from our model disks are thus expected to be optically thin exceptfor limited regions where the molecular column densities have a peak; the radial profile ofthe emission (e.g. ring emission) basically reflects the column density distributions. 20 – In IM Lup, DCO + emission shows a double-ring structure with radii of 110 au and 310au, while H CO + has a single ring structure located at ∼
130 au. ¨Oberg et al. (2015)attributed the H CO + ring and the inner ring of DCO + to the CO snow line. The centraldip (i.e. a small region of low brightness) in the molecular emission lines is suggested tobe caused by the subtraction of the optically thick dust continuum. Thus we cannot tellif H CO + and DCO + decrease inwards at r .
95 au. The location of the DCO + outerring coincides with the edge of the disk observed in the millimeter continuum. ¨Oberg et al.(2015) argued that in this outermost radius the dust opacity is reduced, which enhances thephotodesorption of CO and thus the DCO + abundance. Our model with mm-sized grains isconsistent with this observation. The HCO + and DCO + column density increases inwardsaround (slightly outside) the CO snow line. Beyond the CO snow line, on the other hand,the DCO + column density increases outwards, as in the outer ring in IM Lup, althoughdesorption via cosmic-rays is more efficient than photodesorption in the cold midplane inour model.Another T Tauri disk, AS 209, shows more compact emission of H CO + and DCO + than does IM Lup. While H CO + shows a ring of r ∼
50 au, the radial size of the DCO + emission is similar to or slightly extended than 90 au, at which continuum emission shows abreak. DCO + also has a central dip, which is shallower than that of H CO + (Huang et al.2017). The midplane dust temperature distribution is estimated from continuum observa-tions (Andrews et al. 2009); it is about 20 K at r ∼
60 au. Hence the ring-like emission ofH CO + could coincide with the CO snow line. Our models indicate that DCO + inside theCO snow line could be formed via reactions (4) and CH D + + CO. Perez et al. (2012) foundthat the dust opacity index β increases outwards, which suggests that grains have grown (atleast) to mm-sizes inside 80 au, while grains are still small at the outer radius. The absenceof an outer DCO + ring could be due to efficient freeze out (and chemical conversion) of COmolecules on small dust grains. C O emission, however, shows a local peak at r ∼ CO + is more clearthan that of DCO + . Although the HCO + abundance in the midplane could decline at smallradii due to a low ionization degree, or sublimation of H CO, H S and H O, which havehigher proton affinities than CO, HCO + is the dominant charge carrier at the disk surface,which makes the distribution of HCO + column density relatively flat. In our model with COand H O depletion, the HCO + column density shows a ring-like structure with a depressionat r . a few au, which indicates that the disk surface might be deficient in carbon at thedip radius. A disk structure model specific for AS 209 and a chemical model with isotopeselective CO photodissociation might also be necessary to account for the central dip of 21 –H CO + .In addition to HCO + , N D + is detected in AS 209. Its emission is offset from the diskcenter, and has a peak around the outer edge of the CO emission, which is consistent withour models.LkCa15 is a transient disk with a central hole of 40 −
50 au in the dust continuum.Distributions of H CO + and DCO + emission in LkCa15 are qualitatively similar to thosein IM Lup; the H CO + emission is rather compact, peaking at ∼
40 au, while the DCO + emission is more extended than the dust disk of radius ∼
200 AU, although the ring-likestructure is much less clear in LkCa 15. Pi´etu et al. (2006) derived a disk temperature of22 ± r = 100 au. Thus DCO + emission comes from both inside and outside the COsnow line. The ring of H CO + emission, on the other hand, could be linked to the hole seenin the dust continuum.Another transient disk, V4046, has a central hole of radius ∼
29 au in the dust con-tinuum. While H CO + emission is diffuse extending to ∼
200 au, DCO + emission has aring-like structure with its peak at ∼
70 au (Huang et al. 2017). Rosenfeld et al. (2012) es-timated the temperature distribution to be T ( r ) = 115( r/ − . K; i.e. the temperatureis 27 K at r = 100 au. The DCO + ring thus could be caused by reactions (4) and CH D + +CO. It is consistent with our model that the H CO + distribution is rather flat and extendsinwards compared with that of DCO + . Rosenfeld et al. (2013) determined that most of thedust mass is confined to a ring with a peak at r = 37 au and FWHM of 16 au. The diskmay not extend beyond the CO snowline, where DCO + could have another (outer) emissionpeak.In MWC 480, both H CO + and DCO + emissions have a peak at ∼
40 au (Huang et al.2017). Despite the relatively high luminosity of the central star (11.5 L ⊙ ), the disk temper-ature is rather low. Pi´etu et al. (2006) constrained the disk temperature to be ∼
20 K at aradius of 20 −
30 au. Akiyama et al. (2013) observed CO ( J = 3 − CO ( J = 1 − CO ( J = 1 − ∼ −
15 K for the layer tracedby CO ( J = 1 −
0) at the radius of 100 au. HCO + and DCO + emission thus seems tooriginate in the region around and slightly outside the CO snow line.In HD 163296, the H CO + ( J = 3 −
2) line features an emission ring peaking at r ∼ r ∼
200 au. The DCO + emission shows a rather broad ring with a peak at r ∼
70 au (Huang et al. 2017). Thetemperature distribution in the disk of HD 163296 was derived by Isella et al. (2016) fromCO observations; the midplane temperature is 23 K at r ∼
80 au. Qi et al. (2015), on theother hand, derived the location of the CO snow line to be r ∼
90 au from observations 22 –of N H + and CO isotopologs. HCO + and DCO + emission thus originates both inside andoutside of the CO snow line. Recently, Salinas et al. (2017) reported spatially resolvedemissions of DCO + , N D + , and DCN. Their analysis shows that DCO + emission can bedivided to three ring-like components at 70 au, 150 au, and 260 au (see also Flaherty et al.2017). The innermost and sencond ring coincide with the emission peaks of DCN and N D + ,respectively, while the outermost ring could arise from the non-thermal desorption of CO,as in the outer ring of IM Lup (Salinas et al. 2017). In our fiducial model, the radialdistribution of the DCO + column density can be divided to three components at r ∼ severalau, 10-30 au, and an outermost radius ( &
100 au). The innermost peak coincides with thelocal peak of DCN, while DCO + originates in non-thermal desorption of CO in the midplaneat the outermost radius, which could be consistent with Salinas et al. (2017). The secondcomponent at r ∼ −
30 au, however, is located inside the local peak of the N D + columndensity ( & −
100 au).So far, the DCO + emission detected in disks always has a central hole or dip (Huang et al.2017; Qi et al. 2015). In our models, the DCO + distribution does have a dip. Right outsidethe dip, DCO + is mainly formed by reaction (4), while D atoms are destroyed mainly bythe reaction with HS inside the dip. Although the sublimation of S-bearing species such asH S seems to be a key for the decline of D atom in our model, a derivation of the physicalcondition for the decline of D atoms and thus of DCO + may not be straightforward, sinceD atoms are chemically active. It should also be noted that a central dip could be causedby subtraction of optically thick dust continuum, rather than by a chemical effect, in somedisks (Huang et al. 2017; Salinas et al. 2017). In our models, both in models with mm-sized dust and dark cloud dust, the radialdistributions of the HCN and DCN column densities are centrally peaked, although theDCN/HCN ratio is low (10 − − − ) at the central region ( r .
10 au) due to high temper-ature. Their column densities drop sharply outside this high temperature region, and thenslightly increase outwards, with the gradient steeper for DCN. In the models with verticaldiffusion, the central peak disappears, since hydrocarbons, from which HCN and DCN areformed, are transported to the disk surface to be destroyed.In the disk of IM Lup, H CN is not detected and DCN emission is weak and diffuse(Huang et al. 2017), while both H CN and DCN are clearly detected and have a centraldip (Huang et al. 2017) in the disk of AS 209. In these disks, HCN and DCN are notcentrally peaked, possibly because the temperature is not high enough to desorb HCN at 23 –the radius traced by current observations (beam size of ∼
60 au), or because turbulentmixing is at work. V4046 and HD163296, on the other hand, show centrally peaked H CNemission and a central dip in DCN emission. If these features simply reflect their columndensity distributions, it is difficult to account for in our fiducial model, in which both HCNand DCN are most abundant in the innermost radii. Our models with turbulent diffusion,however, might produce the observed feature; while DCN is abundant only in the outer coldregions, HCN is abundant at the disk surface even inside 20 au, which could be bright dueto high temperatures.
5. SUMMARY
We investigated deuterium chemistry coupled with the nuclear spin-state chemistry ofH and H +3 and their isotopologs in protoplanetary disks. Our principal findings are asfollows:1. We have found multiple paths for deuterium enrichment. The exchange reactions withD atoms, such as HCO + + D, are found to be effective, in addition to H +3 + HD, CH +3 + HD, and C H +2 + HD, which had been considered in previous studies.2. As discussed in Appendix A, the OPR of H is found to be almost thermal, as long asthe cosmic rays ionize the disk with a rate ∼ − s − . In the cold ( ∼
10 K) midplane,however, the OPR reaches its minimum value, which is higher than the thermal value.The minimum value is determined by the balance between the rates of H formation,which sets the OPR to be 3, and spin conversion via ion-molecule reactions involvingprotons and H +3 mainly. The OPR could reach the thermal value at such cold denseregions, if we take into account the spin conversion on grain surfaces, which has recentlybeen found in laboratory experiments (Ueta et al. 2016). We have also analyzed theOPR of H +3 and H D + , and the abundances of H D + , DCO + , and N D + in the coldmidplane, in part by the use of derived analytical formulae.3. In our models, the contribution of the exchange reaction CH +3 + HD is found to beless significant than that described in Favre et al. (2015). The increasing OPR of H helps the backward reaction at T &
20 K.4. In the disk model with mm-sized grains, reduced freeze-out rates enhance the gaseousmolecular abundances in the cold midplane. DCO + , N D + , and DCN in the outermidplane thus contribute significantly to their column densities. The radial distributionof the DCO + column density has a double-ring structure, similar to the DCO + emission 24 –observed in IM Lup. While the outer ring is caused by the enhanced deuteration of H +3 and less efficient chemical conversion of HD and CO in the outer radii in our model,the inner ring is linked to the CO snow line and depletion of D atoms due to reactionswith HS for example. N D + , on the other hand, is more abundant outside the CO snowline. Gaseous DCN decreases inward, except in the central hot region ( T &
100 K)where it is thermally desorbed.5. If the elemental C/O ratio is higher than unity due to sedimentation of H O ice,hydrocarbons become abundant. The exchange reaction CH +3 + HD, which eventuallyleads to the ion CH D + , thus contributes to form DCO + via the reaction CH D + +CO in the warm molecular layer. The column densities of gaseous HCN and DCNare also enhanced by an order of magnitude compared with the fiducial model. Thespatial distributions of molecular abundances and radial profiles of molecular columndensities are, however, qualitatively similar to those in our fiducial model.6. In the disk model with dark cloud dust, the freeze out of molecules is more severe inthe outer midplane, while the disk surface is better shielded from UV radiation than inthe model with mm-sized grains. The disk surface area thus harbors abundant gaseousmolecules and contributes to the column densities (and emissions) of HCO + , DCO + ,N H + , HCN and DCN. One exception is N D + , which is not abundant in the disksurface.7. Turbulence helps to prevent chemical conversion of molecules to less volatile speciesby transporting ices from the midplane to the disk surface, where ices are desorbedand photodissociated, but it also enhances the formation of saturated or less volatilemolecules in the midplane by transporting H atoms and radicals from the disk surface.For example, NH ice formation is hampered and thus N and N H + abundances tendto be enhanced inside the N snow line. Turbulence also transports D atoms in thedisk surface to the lower layers, which helps the formation of DCO + and N D + . Atthe innermost radii, the abundances of HCN and DCN are significantly reduced by theturbulence; their icy components and hydrocarbons are transported to the disk surfaceand destroyed, while O atoms are transported to the midplane to react with C-bearingspecies.8. If the penetration of cosmic rays is hampered by the stellar wind, the midplane ioniza-tion rate decreases by more than one order of magnitude. Column densities of molecularions decline, although the decrement varies with radius, species, and the dust grainsizes in the disk model. HCN and DCN also decrease at r > several au, since thecosmic-ray ionization is needed to form their precursors, N atoms and hydrocarbons,from N and CO, respectively. 25 –This work is supported by JSPS KAKENHI Grant Numbers 23540266, 16H00931, and17K14245. We would like to thank Hideko Nomura for providing the disk physical models,and to Jane Huang, Megumi Kayanuma, and Liton Majumdar for helpful discussions. Weare grateful to the anonymous referee for helpful comments, which improved the manuscript.E. H. wishes to acknowledge the support of the National Science Foundation through grantAST-1514844. A. H OPR
In our chemical reaction network, the OPR is assumed to be 3 for newly formed H molecules on grain surfaces. Neglecting rapid thermalization on the grain, they are desorbedinto the gas phase, where the spin state changes via proton exchange reactions with ionssuch as H +3 . Figure 13 (a) shows the 2D distribution of gas temperature, and Figure 13(b) shows the 2D distribution of the H OPR at t = 3 × yr in our fiducial model. Thedistribution is the same at t = 1 × yr, 3 × yr and 9 . × , because the H OPRreaches steady state on a relatively short timescale ( < yr). In Figure 13 (a) and (b) thehorizontal axis ( r ) is shown in linear scale, since we discuss only the low temperature regionsin this Appendix. The OPR apparently follows the gas temperature distribution. In orderto study the dependence of OPR on gas temperature more quantitatively, we plot the OPRas a function of gas temperature in Figure 13 (c). At low temperatures ( . several tens ofK), where only the lowest two rotational levels need be considered, the OPR is mostly equalto the value of thermal equilibrium, n (o-H ) /n (p-H ) = 9 exp( − . /T ) , (A1)which is depicted by the solid line, while at higher temperatures the OPR asymptoticallyreaches the statistical equilibrium value of 3. At T ∼
10 K, however, the OPR obtained inthe numerical calculation reaches minimum value, which is higher than the thermal value.The minimum value can be explained as follows.The steady state value of the H OPR is given byOPR st = β + b β − b ) β (A2) β = τ op /τ po β = τ op /τ H2 , where b (= 0 .
75) is the branching ratio to form o-H for H formation on grains, τ H2 is thetimescale of H formation, and τ op and τ po represent the timescale of spin conversion via 26 –proton exchange from ortho to para and from para to ortho, respectively (e.g. Furuya et al.2015). When the abundance of ions such as H +3 is high enough, β dominates in (A2), andthe OPR reaches the thermal value. If the ion abundance is low, on the other hand, the OPRis determined by β ; i.e. the balance between the H formation, which sets the OPR to 3,and the spin conversion via ion-molecule reactions. In the disk model, the low temperature( T ∼
10 K) region corresponds to the outer midplane; due to the low ionization degree, theOPR is set by β : β = πa q kTπm n (grain) n (H) k op n (XH + ) n (H ) , (A3)where k is the Boltzmann constant, m is the mass of a hydrogen atom, n (XH + ) is the numberdensity of ions which convert ortho- to para-H through proton exchange reactions with arate coefficient k op . Considering the balance between the formation and destruction of H ,the numerator, i.e. the formation rate of H , can be replaced by the H destruction rate viacosmic-ray/X-ray ionization; β = ζ n (H ) k op n (XH + ) n (H ) = ζk op n (XH + ) , (A4)where ζ is the H ionization rate. Then the minimum value is almost equivalent to the OPRgiven by Walmsley et al. (2004). It should also be noted that the minimum value of OPR inour disk model is much lower than ∼ − , which is often assumed in the model of molecularclouds (e.g. Flower et al. 2006; Faure et al. 2013), because n (XH + ) in eq. (A4) dependson the ionization rate and gas density.Recent laboratory experiments found that the OPR of H is thermalized via the spinconversion on amorphous water ice (Ueta et al. 2016). The conversion rate is higher thanthe desorption rate at .
12 K. Thus OPR could reach the thermal value even at 10 K, if weinclude the conversion on grain surfaces.Figure 13 (c) shows that there are some deviations of the H OPR from the thermal valueat warm temperatures ( T &
25 K), as well. Such deviations are caused by the relatively lowionization degree in the warm dense midplane. The OPR cannot simply be calculated from β (eq. A4) in these regions, since β and β are comparable, unlike the low temperaturelimit ( T ∼
10 K) discussed above.Finally it should be noted that, in the model without cosmic-ray ionization, the regionwith non-thermal H OPR extends both in the radial and vertical directions (Figure 13 (d)).The ionization rate is too low around the midplane for β to dominate in equation (A2). 27 – B. Analytical formula of the OPR values and abundances of molecular ions inthe outer midplane In §
3, we showed that the major formation paths of deuterated molecular ions varyspatially. It is thus difficult to derive analytical formulae of abundances of these ions whichare applicable to the whole disk. Derivation of the analytical formulas is possible, however,if we restrict ourselves to the cold outer midplane. Since molecules such as H D + and N D + have their peak abundances in the outer midplane, and since they are often considered tobe tracers of ionization degree, the analytical formulae could be useful.First, we derive the OPR of H +3 . The analytical formula for the H +3 abundance is givenin eq (18) in Aikawa et al. (2015); it is determined by the formation via cosmic-ray/X-rayionization of H and the destruction by reactions with CO, N , electrons, and negativelycharged grains. Its spin state changes via reactions with H (Oka 2004):o-H +3 + p-H → p-H +3 + o-H (B1)o-H +3 + o-H → p-H +3 + o-H (B2)p-H +3 + o-H → o-H +3 + p-H (B3)p-H +3 + o-H → o-H +3 + o-H . (B4)The OPR of H +3 is thus given at steady state by n (o-H +3 ) n (p-H +3 ) = ( k B3 + k B4 ) n (o-H ) k B1 n (p-H ) + k B2 n (o-H ) . (B5)Figure 14 (a) shows the H +3 OPR in our fiducial disk model (green crosses) and the valuegiven by equation (B5) (red line). We fix the H OPR of (A1) and its minimum value of1 × − , for simplicity. We see reasonable agreement between (B5) and the OPR of H +3 in the numerical calculation, while there are some deviations at T &
25 K, which originatefrom the non-thermal OPR of H (Figure 13 c). The deviation stands out in Figure 14 (a),in which the vertical axis is shown in linear scale, while it is shown in logarithmic scale inFigure 13 (c).The major formation and destruction paths of H D + areo-H +3 + HD → o-H D + + o-H (B6)p-H +3 + HD → o-H D + + p-H (B7)p-H +3 + HD → p-H D + + o-H (B8)p-H +3 + HD → p-H D + + p-H (B9)o-H D + + o-H → o-H +3 + HD (B10) 28 –p-H D + + o-H → o-H +3 + HD (B11)p-H D + + o-H → p-H +3 + HD (B12)o-H D + + CO → HCO + + HD (B13)p-H D + + CO → HCO + + HD (B14)o-H D + + CO → DCO + + o-H (B15)p-H D + + CO → DCO + + p-H . (B16)The balance between these reactions is described as n (HD) { kB n (o-H +3 ) + ( kB kB kB n (p-H +3 ) } = (B17) n (o-H ) { kB n (o-H D + ) + ( kB
11 + kB n (p-H D + ) } +( kB
13 + kB n (o-H D + ) n (CO) + ( kB
14 + kB n (p-H D + ) n (CO) . The spin state of H D + is converted viap-H D + + o-H ⇋ o-H D + + p-H (B18)o-H D + + o-H ⇋ p-H D + + p-H (B19)o-H D + + o-H ⇋ p-H D + + o-H (B20)in a shorter timescale than that of H D + /H +3 conversion. The o/p ratio of H D + can thusbe described by n (o-H D + ) n (p-H D + ) = ( k + B
18 + k − B × n (o-H ) /n (p-H ) + k − B k + B
19 + k + B × n (o-H ) /n (p-H ) + k − B
18 (B21)(Gerlich et al. 2002; Br¨unken et al. 2014). The H D + abundance is then given by n (H D + ) = n (HD) { kB n (o-H +3 ) + ( kB kB kB n (p − H +3 ) } kB n (o − H ) xx +1 + ( kB
11 + kB n (o-H ) x +1 + ( kB
13 + kB n (CO) , (B22)where x is the o/p ratio of H D + (B21). Note that the rate coefficient of H D + + CO is thesame for o-H D + and p-H D + . If the H D + abundance derived above is comparable to thatof main isotopolog, H +3 , it should be replaced by n (H D + ) = n (H +3 ) n ′ (H D + ) n (H +3 ) + n ′ (H D + ) , (B23)where n ′ (H D + ) denotes the value obtained in (B22).Reasonable agreement between the analytical value and the numerical results for the o/pratio of H D + is found in Figure 14 (b). In our numerical models, multiply-deuterated H +3
29 –becomes more abundant than H D + in the outermost ( r &
250 au) midplane. The analysisof multi-deuterated species is, however, beyond the scope of this paper. Since we neglectthe multiply deuterated isotopologs of H +3 , H D + represents the total deuterated-H +3 in ouranalysis. In the following, we therefore compare the analytical value of the H D + abundancewith the sum of the abundances of deuterated H +3 ( n (D +3 ) /n H + n (D H + ) /n H + n (H D + ) /n H )obtained in our fiducial model. Figure 14 (c) shows the sum of deuterated H +3 abundancein our fiducial model (green), which is compared with the analytical value (red). At warmregions ( T &
20 K), the abundances of H +3 and its deuterated isotopologs are not simplefunctions of temperature; they depend on abundances of other species such as CO, N andelectrons (Aikawa et al. 2015). We use the latter abundances from our fiducial model inthe calculation of the analytical formula of H D + abundance. Equation (B22) also needsthe abundance of HD as an input. While its canonical value is 1 . × − , HD abundance isdecreased down to ∼ − via conversion to D-bearing ices such as NH D and HDO in thecold midplane in our numerical model (see also Sipil¨a et al. 2013; Teague et al. 2015). Wethus use the HD abundance from our fiducial model, as well.Figures 15 (a) and (b) show the 2D distribution of deuterated H +3 in our fiducial model,and the H D + abundance derived from eq. (B22), respectively. We can see a reasonableagreement between the numerical model and the analytical formula. For a comparison,Figure 15 (c) shows the distribution of H D + abundance derived from eq. (B22) assumingthat CO, N and HD are not converted to other molecules; HD abundance is set to beconstant ( n (HD)/ n H = 1 . × − ), and CO and N abundances are determined simply bythe balance between adsorption and desorption (Aikawa et al. 2015). H D + abundance ishigher than the numerical model, mainly due to a higher HD abundance.Once we obtain the analytical formula of the H D + abundance, we can calculate theDCO + and N D + abundances, by simply replacing H +3 with H D + and modifying the corre-sponding reaction rate coefficients in equations (22) and (21) in Aikawa et al. (2015). Theresultant 2D distributions are shown in Figure 15 (d-i). Panels (e) and (h) show distributionsof DCO + and N D + , respectively, derived from the analytical formulas using the abundancesof HD, CO, N , and electrons from our fiducial model, while the abundances of HD, CO andN are given analytically for panels (f) and (i). The agreements between the numerical model(panels d and g) and analytical values (panels e and h) confirm that DCO + and N D + aremainly formed via deuterated H +3 in the outer midplane. REFERENCES
Adams, N. G. & Smith, D. 1985, ApJ, 294, L63 30 –Table 1: Initial abundances of major speciesspecies abundanceN 1.8 (-8) a O 8.8 (-9)o-H CO 3.5 (-6)CO OH 8.8 (-7)N D 2.0 (-7)H O 1.1 (-4)HDO 8.5 (-8) a a ( b ) means a × b .
31 –Aikawa, Y., Miyama, S. M., Nakano, T., & Umebayashi, T. 1996, ApJ, 467, 684Aikawa, Y., Umebayashi, T., Nakano, T. & Miyama, S. M. 1997, ApJ, 486, L51Aikawa, Y. & Herbst, E. 1999, ApJ, 526, 314Aikawa, Y., Umebayashi, T., Nakano, T., & Miyama, S.M. 1999, ApJ, 519, 705Aikawa, Y. & Nomura, H. 2006, ApJ, 642, 1152Aikawa, Y., Wakelam, V., Hersant, F., Garrod, R.T. & Herbst, E. 2012, ApJ, 760, 40Aikawa, Y., Furuya, K., Nomura, H. & Qi, C. 2015, ApJ, 807, 120Akiyama, E., Momose, M., Kitamura, Y. et al. 2013, PASJ, 65, id. 123Albertsson, T., Semenov, D., & Henning, Th. 2014, ApJ, 748, 39Altwegg, K., Balsiger, H., Bar-Nun, A. et al. 2015, Science, 347, id. 1261952Andrews, S. M., Wilner, D. J., Hughes, A. M., Qi, C., & Dullemond, C.P. 2009, ApJ, 700,1502Bacmann, A., Lefloch, B., Ceccarelli, C. et al. 2003, ApJ, 585, L55Balbus, S.A. & Hawley, J. F. 1991, ApJ, 376, 214Bergin, E. A., Du, F., Cleeves, L. I. et al. 2016, ApJ, 831, 101Br¨unken, S., Sipil¨a, O., Chambers, E.T. et al. 2014, Nature, 516, 219Caselli, P., Stancheva, T., Shalabiea, O., Shematovich, V. I., & Herbst, E. 2002, P&SS, 50,1257Caselli, P., Hasegawa, T.I., & Herbst, E. 1998, ApJ, 495, 309Caselli, P., van der Tak, F.F.S., Ceccarelli, C., & Bacmann, A. 2003, A&A, 403, L37Cazzaux, S. & Tielens, A.G.G.M. 2004, ApJ, 604, 222Cazzaux, S. & Tielens, A.G.G.M. 2010, ApJ, 715, 698Cleeves, L.I., Adams, F.C. , & Bergin, E. A., 2013, ApJ, 772, 5Cleeves, L. I., Adams, F.C., Bergin, E. A., & Visser, R. 2013, ApJ, 777, 28Cleeves, L. I., Bergin, E.A., Alexander, C.M.O’D. et al. 2014, Science, 345, 1590 32 –Cleeves, L. I., Bergin, E. A., & Adams, F. C. 2014, ApJ, 794, 123Cleeves, L.I., Bergin, E.A., Alexander, C.M. O’D. et al. 2016, ApJ, 819, 13Coutens, A., Vastel, C., Hincelin, U. et al. 2014, MNRAS, 445, 1299Draine, B.T. & Bertoldi, F. 1996, ApJ, 468, 269Du, F., Bergin, E. A., Hogerheijde, M. et al. 2017, ApJ, 842, 98Favre, C., Cleeves, L. I., Bergin, E. A., Qi, C., & Blake, G. A. 2013, ApJ, 776, L38Favre, C., Bergin, E. A., Cleeves, L. I. et al. 2015, ApJ, 802, L23Flaherty, K. M., Hughes, A. M., Rosenfeld, K. A. et al. 2015, ApJ, 813, 99Flaherty, K. M., Hughes, A. M., Rose, S. C. et al. 2017, ApJ, 843, 150Flower, D. R., Pineau des Forˆets, & Walmsley, C. M. 2006, A&A, 449, 621Faure, A., Hily-Blant, P., Le Gal, R., Rist, C., & Pineau des Forˆets 2013, ApJ, 770, L2Furuya, K., Aikawa, Y., Nomura, H., Hersant, F., & Wakelam, V. 2013, ApJ, 779, 11Furuya, K. & Aikawa, Y. 2014, ApJ, 790, 97Furuya, K., Aikawa, Y., Hincelin, U. et al. 2015, A&A, 584, A124Furuya, K., van Dishoeck, E. F., & Aikawa, Y. 2016, A&A, 586, A127Garrod, R.T. & Herbst, E. 2006, A&A, 457, 927Garrod, R.T., Wakelam, V., & Herbst, E. 2007, A&A, 467, 1103Geiss, J. & Gloeckler 1998, SSRv, 84, 239Gerlich, D., Herbst, E., & Roueff, E. 2002, P&SS, 50, 1275Harada, N., Herbst, E., & Wakelam, V. 2010, ApJ, 721, 1570Harada, N., Herbst, E., & Wakelam, V. 2012, ApJ, 756, 104Hasegawa, T.I. & Herbst, E. 1993, MNRAS, 261, 83Herczeg, G. J., Linsky, J. L., Valenti, J. A., Johns-Krull, C. M., & Wood, B. E. 2002, ApJ,572, 310 33 –Hincelin, U., Herbst, E., Chang, Q. et al. 2014, 69th International Symposium on MolecularSpectroscopy, MF09Hogerheijde, M. R., Bergin, E. A., Brinch, C. et al. 2011, Science, 334, 338Huang, J., & ¨Oberg, K.I. 2015, ApJ, 809, L26Huang, J., ¨Oberg, K.I., & Andrews, S. 2016, ApJ, 823, L18Huang, J., ¨Oberg, K.I., Qi, C. et al. 2017, ApJ, 835, 231Hugo, E., Asvany, O., Schlemmer, S. 2009, J. Ch. Ph., 130, 164302Isella, A., Guidi, G., Testi, L. et al. 2016, PRL, 117, 251101Johansen, A. & Klahr, H. 2005, ApJ, 634, 1353Kama, M., Bruderer, S., van Dishoeck, E. F. et al. 2016, A&A, 592, 83Kamp, I. & Bertoldi, F. 2000, A&A, 353, 276Kastner, J. H., Huenemoerder, D. P., Schulz, N. S., Canizares, C. R., & Weintraub, D. A.2002, ApJ, 567, 434Krijt, S., Ciesla, F. J., & Bergin, E. A., 2016, ApJ, 833, 285Li, X., Heays, A. N., Visser, R. et al. 2013, A&A, 555, A14Linsky, J.L. 2003, SSRv, 106, 49Mathews, G.S., Klaassen, P.D., Juh´asz, A. et al. 2013, A&A, 557, A132Meijerink, R., Pontoppidan, K. M., Blake, G.A., Poelman, D.R., & Dullemond, C. P. 2009,ApJ, 704, 1471Millar, T.J., Bennet, A, & Herbst, E. 1989, ApJ, 340, 906Minissale, M., Dulieu, F., Cazaus, S. & Hocuk, S. 2016, A&A, 585, A24Mumma, M.J. & Charnley, S.B. 2011, ARAA, 49, 471Noble, J.A., Theule, P., Borget, F. et al. 2013, MNRAS, 428, 3262Nomura, H. & Millar, T. 2005, A&A, 438, 923Nomura, H., Aikawa, Y., Tsujimoto, M., Nakagawa, Y., & Millar, T., 2007, ApJ, 661, 334 34 –Nomura, H., Tsukagoshi, T., Kawabe, R. et al. 2016, ApJ, 819, L7¨Oberg, K.I., van Dishoeck, E. F., & Linnartz, H. 2009, A&A, 496, 281¨Oberg, K.I., Qi, C., Fogel. K.J. et al. 2011, ApJ, 734, 98¨Oberg, K.I, Qi, C., Wilner, D. J., & Hogerheijde, M. R. 2012, ApJ, 749, 162¨Oberg, K.I., Furuya, K., Loomis, R. et al. 2015, ApJ, 810, 112Oka, T., 2004, JMoSp, 228, 635Penteado, E. M., Walsh, C. & Cuppen, H. M. 2017, ApJ, 844, 71P´erez. L. et al. 2012, ApJ, 760, L17Pi´etu, V., Dutrey, A., Guilloteau, S., Chapillon, E., & Pety, J. A&A 460, L43-L47Qi, C., Wilner, D.J., Aikawa, Y., Blake, G.A., & Hogerheijde, M. R. 2008, ApJ, 681, 1396Qi, C., ¨Oberg, K.I., Andrews, S. M. et al. 2015, ApJ, 813, 128Roberts, H. & Millar, T.J. 2000, A&A, 361, 388Rosenfeld, K.A., Andrews, S. M., Wilner, D. J., & Stempels, H. C. 2012, ApJ, 759, 119Rosenfeld, K. A., Andrews, S. M., Wilner, D. J., Kastner, J. H., & McClure, M. K. 2013,ApJ, 775, 136Roueff, E., Gerin, M., Lis, D.C. et al. 2013, J. Phys. Chem. A, 117, 9959Ruaud, M., Wakelam, V. & Hersant, F. 2016, MNRAS, 459, 3756Salinas, V. N., Hogerheijde, M. R., Matthews, G. S. et al. 2017, A&A, 606, A125Schilke, P., Walmsley, C.M., Pineau des Forets, G. et al. 1992, A&A, 256, 595Schwarz, K. R., Bergin, E. A., Cleeves, L. I. et al. 2016, ApJ, 823, 91Sipil¨a, O., Caselli, P., & Harju, J. 2013, A&A, 554, A92Smith, D., Adams, N.G., & Alge, E. 1982, ApJ, 263, 123Teague, R., Semenov, D., Guilloteau, S. et al. 2015, A&A, 574, A137Teague, R., Guilloteau, S., Semenov, D. et al. 2016, A&A, 592, A49 35 –Tielens, A.G.G.M. & Hollenbach, D., 1985, ApJ, 291, 722Umebayashi, T., & Nakano, T. 1981, PASJ, 33, 617Umebayashi, T. & Nakano, T. 2009, ApJ, 690, 69Ueta, H., Watanabe, N., Hama, T. & Kouchi, A. 2016, PRL, 116, 253201van der Tak, F.F.S., Black, J.H., Sch´’oier, F.L., Jansen, D.J., van Dishoeck, E.F. 2007, A&A468, 627van Dishoeck, E.F., Thi, W. -F., & Zadelhoff, G.-J. 2003, A&A, 2003, 400, L1Vastel, C., Caselli, P., Ceccarelli, C. et al. ApJ, 645, 1198Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, A&A, 503, 323Wolcott-Green, J. & Haiman, Z. 2011, MNRAS, 412, 2603Walmsley, C. M., Flower, D. R., & Forˆets 2004, A&A, 418, 1035Weingartner, J. C. & Draine, B. T., 2001, ApJ, 548, 296Willacy, K. 2007, ApJ, 660, 441Xu, R., Bai, X.-N., & ¨Oberg, K. I. 2017, ApJ, 835, 162Zhang, K., Bergin, E. A., Blake, G. A., Cleeves, L. I., Schwarz, K. R. 2017, NatAs, 1, 0130
This preprint was prepared with the AAS L A TEX macros v5.2.
36 –Table 2: Observations morphologies of molecular lines intensities aObject Class H CO + DCO + N H + N D + H CN DCNIM Lup full disk (T Tauri) ring-like multiple rings offset b - ring-like c extended? d AS 209 full disk (T Tauri) ring-like ring-like offset b offset e ring-like ring-likeLkCa 15 transition disk (T Tauri) ring-like diffuse - - compact? d multiple rings? d V4046 Sgr transition disk (T Tauri) diffuse ring-like offset b - centrally peaked ring-likeMWC 480 Herbig Ae ring-like ring-like - - centrally peaked compact? d HD 163296 Herbig Ae multiple rings multiple rings f ,g - ring-like g centrally peaked ring-like a Based on Huang et al. (2017) unless otherwise stated. b ¨Oberg et al. (20011) c H CN is not detected. H CN shows a ring-like feature. d Detected but not categorized in Huang et al. (2017) due to relatively low S/N. e Huang & ¨Oberg (2015) f Flaherty et al. (2017) g Salinas et al. (2017)
37 –Table 3: Rate coefficients of major reactions of deuterated molecular ions a reaction rate coefficients [cm s − ]B1 1 . × − exp( − /T ) b B2 4 . × − exp(0 . /T )B3 3 . × − exp(0 . /T )B4 8 . × − exp( − . /T )B6 1 . × − exp( − . /T )B7 6 . × − exp(1 . /T )B8 4 . × − exp(0 . /T )B9 3 . × − exp(0 . /T )B10 1 . × − exp( − . /T )B11 9 . × − exp( − . /T )B12 1 . × − exp( − . /T )B13 1 . × − B15 5 . × − B18 + . × − exp( − . /T )B18 − . × − exp( − . /T )B19 + . × − exp(0 . /T )B19 − . × − exp( − /T )B20 + . × − exp(0 . /T )B20 − . × − exp( − . /T ) a Oka (2004), Hugo et al. (2009), Coutens et al. (2014) and references therein. b The unit of temperature ( T ) is the Kelvin.
38 –Fig. 1.— Distributions of gas density (top), dust temperature (middle), and gas temperature(bottom) in the disk model with mm-sized grains (left) and dark cloud dust (right). Theupper dotted line depicts the positions where the X-ray ionization rate is equal to the cosmic-ray ionization rate (5 × − s − ), while the X-ray ionization rate is 1 × − s − at thelower dotted line. 39 –Fig. 2.— Panels (a), (c), and (d) depict distributions of CO, HCO + , and DCO + in the gasphase in the disk model with mm-sized grains at 3 × yr. The long dashed lines depict theCO snow surface, while the dotted lines are the same as in Figure 1. Panel (e) is color codedreferring to the major formation pathways of DCO + at each position in the disk. Panels(b) and (f) show the radial distribution of molecular column densities at 1 × yr (blue),3 × yr (green), and 9 . × yr (red). 40 –Fig. 3.— The ratio of backward to forward reaction rate coefficients of CH +3 + HD → CH D + + H (reaction 3) as a function of gas temperature given in Roueff et al. (2013). The solidline depicts the value with the thermal OPR of H , while the dashed line depicts the valuewith a constant exothermicity of 654 K assuming all H in the para form. The ratio of theformer to the latter is shown by the dotted line. 41 –Fig. 4.— Distributions of molecular abundances (a, c, and d) and molecular column densities(b and f) as in Figure 2, but for N , N H + , and N D + in the gas phase. The long dashedlines in panels (a), (c) and (d) depict the N snow surface. Panel (e) is color coded referringto the major formation pathways of N D + at each position in the disk.Fig. 5.— Distribution of molecular abundances (a-b) and molecular column densities (d)as in Figure 2, but for HCN, and DCN in the gas phase. The long dashed lines in panels(a) and (b) depict the HCN snow surface. Panel (c) is color coded referring to the majorformation pathways of DCN at each position in the disk. 42 –Fig. 6.— Same as Figure 2, but for the disk model with depletions of CO and H O. Inpanel (e), DCO + is formed mainly via H D + + CO in the pink region, while the reactionsof multi-deuterated H +3 with CO dominate in the white region. 43 –Fig. 7.— Distribution of molecular abundances (a, c, d, g and h), and molecular columndensities (b, f, and j) as in Figure 2, but for N , N H + , N D + , HCN and DCN in the gasphase in the disk model with depletions of CO and H O. The long dashed lines depict thesnow surface of the mother molecule. Panels (e) and (i) are color coded referring to themajor formation pathways of N D + and DCN, respectively, at each position in the disk. Inpanel (e), the reactions of multi-deuterated H +3 with N dominate in the white region, whilethe reaction of H D + + N dominates in the purple region. 44 –Fig. 8.— Same as Figure 2, but for the disk model with dark cloud dust.Fig. 9.— Same as Figure 7, but for the disk model with dark cloud dust. 45 –Fig. 10.— Distributions of gaseous CO, N , HCN and DCN in the mm-sized grain model withvertical mixing at t = 3 × yr. The diffusion coefficient is α = 10 − in the top panels, and α = 10 − in the middle panels. The bottom panels show the radial distribution of molecularcolumn densities of models without diffusion (solid) and with a diffusion coefficient of 10 − (dashed), and 10 − (dotted) at t = 3 × yr. 46 –Fig. 11.— Same as Figure 10, but for ionic molecules. 47 – (cid:16116)(cid:16111)(cid:16123)(cid:16087)(cid:16112)(cid:16111)(cid:16123)(cid:16087) (cid:16122)(cid:16094)(cid:16116)(cid:16087)(cid:16122)(cid:16094)(cid:16112)(cid:16087) (cid:16116)(cid:16111)(cid:16122)(cid:16112)(cid:16111)(cid:16122) Fig. 12.— Distributions of abundances and column densities of HCO + , N H + , HCN andtheir deuterated isotopologs in the mm-sized grain model without cosmic-ray ionization at t = 3 × yr. The solid lines depict the molecular column density without cosmic-rayionization, while the dashed lines depict the column density in our fiducial model. 48 –Fig. 13.— (a) Distribution of gas temperature in our fiducial disk model. (b) Distributionof H OPR in our fiducial disk model. (c) H OPR in our fiducial model (red crosses) andthermal equilibrium value (blue) as functions of temperature. (d) Same as (b), but for themodel without cosmic-ray ionization. The upper and lower dotted lines depict the positionswhere the X-ray ionization rate is equal to 5 × − s − and 1 × − s − , respectively. 49 –Fig. 14.— (a) OPR of H +3 in our fiducial model (green crosses) and thermal equilibriumvalue (red) as a function of temperature. (b) The o/p ratio of H D + in our fiducial model(green) and the analytical value (red). (c) H D + abundance in our fiducial model (green)and the value obtained by the analytical formula (red). 50 –Fig. 15.— Panels (a), (d), and (g) show distributions of deuterated H +3 , DCO + and N D + in our fiducial model at t = 3 × yr. Panels (b), (e), and (h) show distributions of H D + ,DCO + and N D + calculated by the analytical formulas using the abundances of HD, CO, N ,and electrons from the fiducial model. Panels (c), (f), and (i) show distributions of H D + ,DCO + and N D + calculated by the analytical formulas with the abundances of HD, CO andN2