Using HCO + isotopologues as tracers of gas depletion in protoplanetary disk gaps
Grigorii V. Smirnov-Pinchukov, Dmitry A. Semenov, Vitaly V. Akimkin, Thomas Henning
AAstronomy & Astrophysics manuscript no. 38572corr c (cid:13)
ESO 2020September 22, 2020
Using HCO + isotopologues as tracers of gas depletion inprotoplanetary disk gaps Grigorii V. Smirnov-Pinchukov , Dmitry A. Semenov , , Vitaly V. Akimkin , and Thomas Henning Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germanye-mail: [email protected] Department of Chemistry, Ludwig Maximilian University, Butenandtstr. 5-13, D-81377 Munich, Germany Institute of Astronomy, Russian Academy of Sciences, Pyatnitskaya str. 48, Moscow, 119017, RussiaReceived September 22, 2020; accepted
ABSTRACT
Context.
The widespread rings and gaps seen in the dust continuum in protoplanetary disks are sometimes accompanied by similarsubstructures seen in molecular line emission. One example is the outer gap at ∼
100 au in AS 209, which shows that the H CO + andC O emission intensities decrease along with the continuum in the gap, while the DCO + emission increases inside the gap. Aims.
We aim to study the behavior of DCO + / H CO + and DCO + / HCO + ratios in protoplanetary disk gaps assuming the two scenar-ios: (A) the gas depletion follows the dust depletion and (B) only the dust is depleted. Methods.
We first modeled the physical disk structure using the thermo-chemical model ANDES. This 1 +
1D steady-state disk modelcalculates the thermal balance of gas and dust and includes the far ultraviolet (FUV), X-rays, cosmic rays, and other ionization sourcestogether with the reduced chemical network for molecular coolants. Afterward, this physical structure was adopted for calculationsof molecular abundances with the extended gas-grain chemical network with deuterium fractionation. Ideal synthetic spectra and0th-moment maps were produced with the LIne Modeling Engine (LIME).
Results.
We are able to qualitatively reproduce the increase in the DCO + intensity and the decrease in the H CO + and C O intensitiesinside the disk gap, which is qualitatively similar to what is observed in the outer AS 209 gap. The corresponding disk model (A)assumes that both the gas and dust are depleted in the gap. The model (B) with the gas-rich gap, where only the dust is depleted,produces emission that is too bright in all HCO + isotopologues and C O. Conclusions.
The DCO + / H CO + line ratio can be used to probe gas depletion in dust continuum gaps outside of the CO snow line.The DCO + / C O line ratio shows a similar, albeit weaker, e ff ect; however, these species can be observed simultaneously with a singleALMA or NOEMA setup. Key words. astrochemistry – methods: numerical – protoplanetary disks – stars: pre-main sequence – ISM: molecules – submillime-ter: planetary systems
1. Introduction
Protoplanetary disks are believed to be the birthplaces of plan-etary systems. While more and more theoretical studies of theplanet formation in disks appear in the literature, they need ob-servational constraints regarding physical conditions and chem-ical composition in various disks and disk locations. Now, inthe Atacama Large (sub)Millimeter Array (ALMA), NOrth-ern Extended Millimeter Array (NOEMA), Very Large Tele-scope / Spectro-Polarimetric High-contrast Exoplanet REsearch(VLT / SPHERE), and Gemini Planet Imager (GPI) era, the spa-tially resolved observations of dust and gas in protoplanetarydisks in nearby star-forming regions have become a routine.High-resolution studies in dust continuum and scattered lightwith an angular resolution up to 0 . (cid:48)(cid:48) − . (cid:48)(cid:48) reveal com-plex substructures with multiple gaps, inner holes, rings, andspirals (DSHARP: Andrews et al. 2018; Huang et al. 2018b;Long et al. 2019; Huang et al. 2020; van Boekel et al. 2017).The high-resolution observations in the CO isotopologues andother molecules also show the presence of substructures in thedisk gas (see, e.g., Isella et al. 2016; Fedele et al. 2017; Teagueet al. 2017; Huang et al. 2018a). Often disk substructures thatare detected in the scattered light at near-infrared wavelengths do not coincide with the substructures visible in the submillime-ter dust continuum and gas emission lines (Keppler et al. 2019).One of the peculiar disks showing very wide gaps is AS 209.Huang et al. (2017) have studied continuum and line emission inthis disk at the ∼ . (cid:48)(cid:48) resolution with ALMA. They have founda notable di ff erence between H CO + and DCO + emission radialprofiles, with the DCO + emission peak extending farther out upto radii of (cid:46) −
100 au. The follow-up higher resolution ob-servations have shown that the AS 209 disk has a dust gap ataround 55–120 au (Fedele et al. 2018; Favre et al. 2019). At thisgap location, the optically thin C O J = CO J = + J = + isotopologues could be used as a probe forgas depletion in disk gaps that have been detected before in dustcontinuum or CO isotopologues.There are several major factors that can cause molecular lineemission in disks to show gap-like or ring-like substructures.Firstly, a circular gap or a ring in the line emission could bedriven by a local, rather azimuthally symmetric deviation in the Article number, page 1 of 12 a r X i v : . [ a s t r o - ph . E P ] S e p & A proofs: manuscript no. 38572corr physical structure of the disk, either in density or temperature(e.g., due to planet-disk interactions, snowlines, a change in dustproperties, etc. Lin & Papaloizou 1993; Guzmán et al. 2018;Huang et al. 2020). In the case of a surface density drop, the totalmolecular column density could become smaller or higher (e.g.,for ions), which would lead to weaker or stronger line emission,as long as the line remains optically thin. Similarly, a local in-crease or decrease in the gas’s kinetic temperature or backgroundradiation a ff ect molecular line excitation and hence emission in-tensity, assuming that the line is thermalized.Secondly, molecular emission could show rings or gapscaused primarily by various chemical e ff ects. The most promi-nent ones are ultraviolet (UV)-driven production or (selective)photodissociation in the disk atmosphere, C / O variations, andthe freeze-out of molecules onto the dust grain surfaces in thedisk midplane (see, e.g., Dutrey et al. 2017; Teague et al. 2017;Cazzoletti et al. 2018; Miotello et al. 2018, 2019; van Terwisgaet al. 2019; Garufi et al. 2020). For example, c-C H , CCH, CN,H CO, and CS emission often shows a ring-like appearance indisks. Whereas for hydrocarbons and CN emission, the C / O vari-ations and high-energy UV or X-ray irradiation are the most im-portant factors; for CS and H CO emission, the surface chem-istry processes and freeze-out play the major role (Bergin et al.2016; Cazzoletti et al. 2018; Miotello et al. 2019; van Terwisgaet al. 2019). The ring-like appearance of N D + and DCO + emis-sion in disks is mainly due to low-temperature deuterium frac-tionation when gaseous CO freezes out in the disk midplane (Qiet al. 2013; Ceccarelli et al. 2014; Huang & Öberg 2015; Salinaset al. 2018; Garufi et al. 2020). In addition, temporal variationsin physical conditions can also a ff ect chemical abundances, forexample, via X-ray flares or episodic accretion outbursts (e.g.,Cleeves et al. 2017; Molyarova et al. 2018; Lee et al. 2019).Finally, a third relevant e ff ect for molecular line excita-tion and emission intensity is the departure from local ther-mal equilibrium (LTE). It depends on the molecular propertiesand the location of the emitting molecular layer in the disk(Pavlyuchenkov et al. 2007). The most notable examples are theCN radical emitting from upper disk layer, where excitation isdominated by radiative processes, and CH OH, which have acomplex energy level structure with torsional and coupled en-ergy levels (e.g., Parfenov et al. 2016; Cazzoletti et al. 2018;Teague & Loomis 2020). Another e ff ect that may a ff ect the lineexcitation in disks locally is the localized deviation in the phys-ical structure of the disk, which does not a ff ect the global struc-ture of the disk otherwise (e.g., local heating in a circumplane-tary disk Cleeves et al. 2015). This shows the complex interplaybetween the physical structure, chemistry, and excitation con-ditions in defining the strength of line emission and its spatialdistribution in a disk.The paper is structured as follows. In the Section 2, wepresent the tools we use and the assumptions we make. In Sec-tion 3 we present the results for each model (3.1-3.3). In Sec-tion 4 we discuss the model uncertainties (4.1) and the compar-ison with previous observational gas gap studies (4.2). We pro-vide the final conclusion in Section 5.
2. Methods and models
In this section we list all of the steps involved in our modeling.Firstly, we describe ANDES, the thermo-chemical code we usedto simulate 1 +
1D disk gas and dust thermal structure (Subsec-tion 2.1). Secondly, we explain how we combined the simulatedvertical columns into three 2D disk models, with one referencemodel and two models with and without gas depletion in the dust gap. We recalculated the radial ionization and photodissoci-ation self-shielding factors for those 2D physical structures (2.2).Thirdly, we briefly overview the chemical model ALCHEMICwith the deuterated network, which we used to calculate thechemical abundances (2.3). Finally, we describe how we com-puted the synthetic molecular line emission maps with LIME(2.4).
The standard practice of chemical modeling of protoplanetarydisks includes some density structure setup, dust radiative trans-fer simulations for the thermal structure of the dust, gas thermalbalance, and chemical kinetics calculations using the resultingstructure (Dust And LInes (DALI) (Bruderer et al. 2009, 2014),PRotoplanetary DIsk MOdel (ProDiMo) (Woitke et al. 2009),ANDES (Akimkin et al. 2013), also (Gorti & Hollenbach 2004;Cleeves et al. 2013; Du & Bergin 2014; Walsh et al. 2015; Sali-nas et al. 2018)). The ANDES model allows one to simulatethe chemical evolution of a protoplanetary disk with a detailedtreatment of gas and dust thermal balance. Its original version(Akimkin et al. 2013) is based on a 1 +
1D approach, where thevertical disk structure at each radius is simulated independently.In the upper, low-density regions of the protoplanetary disk, thegas temperature is not coupled to the dust temperature. The gasdensity becomes too low for thermal accommodation with thedust, and the gas usually becomes hotter than the dust as it coolsitself rather ine ffi ciently (Röllig et al. 2007). The gas density,temperature, and chemical composition depend on each other,so iterative modeling of the vertical disk structure is performed.After the initial chemical assumptions are made, the thermal bal-ance for this predefined chemical structure is calculated. Thechemical kinetics calculations are iterated with thermal balanceuntil convergence in further iterations. In the utilized versionof the ANDES code, the dust temperature step is recalculatedagain for the hydrostatic solution that fits the vertical temper-ature structure, assuming well-mixed small dust particles and,therefore, a constant dust-to-gas mass ratio.The recent version of the ANDES code is based on afully 2D approach and it simultaneously iterates the vertical diskstructure and the chemical evolution, which allows one to accu-rately simulate some time-dependent e ff ects such as luminosityoutbursts in the FU Ori-type systems (Molyarova et al. 2017,2018). However, this recent version does not feature the detailedgas thermal balance. In this study, we use the original versionof ANDES with the reduced chemical network for molecularcoolants and modifications to the modeling of the X-ray and cos-mic ray ionization rates. The reduced chemical network consistsof 63 species and 480 reactions to accurately model the abun-dances of C, C + , H O, CO, CO , OH, H , and other simplespecies. It has been obtained from the full network (described inSubsection 2.3) by using the iterative reaction-based reductionmethod followed by manual tuning (Wiebe et al. 2003; Semenovet al. 2004).Using the ANDES model, we first computed a set of diskthermo-chemical models to obtain 1 +
1D temperature and den-sity structures. We assumed the stellar and disk parameters, aslisted in Table 1. We selected a T Tauri star with a mass of M = . M (cid:12) , and we derived its e ff ective temperature T e ff = R = . R (cid:12) from the evolutionary track model as-suming an age of 2 Myr (Yorke & Bodenheimer 2008, privatecommunication). We did not aim to qualitatively model the K5e-type AS 209 star, rather we modeled a generic T Tauri protoplan-etary disk. For the UV excess, we assumed the blackbody emis- Article number, page 2 of 12rigorii V. Smirnov-Pinchukov et al.: Disk gaps traced by HCO + isotopologues Table 1.
ANDES input parameters.
Parameter Value M star . M (cid:12) T star R star . R (cid:12) UV excess B ν ( T e ff = . × − M (cid:12) yr − Dust opacity Draine & Lee (1984)Dust grain size 10 − cmDust-to-gas mass ratio 0.01, 0.001Grazing angle 0.05Gas column density at 1 au 10, 100 g cm − Density power-law slope − T e ff = T e ff = ff ective temperature into account. Forthe dust radiative transfer and surface chemistry, we used non-evolving single-size silicate grains with a size of 10 − cm. We setthe dust-to-gas mass ratio to 0.01 in our reference disk model.The disk surface density was parameterized using the follow-ing power-law distribution, without taking tapering at outer radiiinto account: Σ ( r ) = Σ (cid:18) r (cid:19) − . (1)This approach leads to a denser outer disk and a ff ects the in-tegral disk properties. On the other hand, it does not a ff ect thee ff ect of the local substructure much. We also separately com-puted radii between 80 and 250 au with a 10 times lower surfacedensity, and with a 10 times lower dust-to-gas ratio for the gapmodels, which is described in the next section. The summary ofthe adopted stellar and disk parameters is presented in Table 1. After that, we combined the vertical structures computed by AN-DES for each disk radii into the three protoplanetary disk mod-els (see Table 2). The first model, the so-called reference model,hereafter referred to as (R), is the disk model without any dustor gas gaps. For the other two disk models, we used the sepa-rately computed vertical structures for the radii between 80 and250 au, and those from the reference model (R) outside of thisgap. For the second model, the so-called gas-poor gap model,hereafter referred to as (A), we lowered gas and dust surfacedensities by a factor of 10 between 80 and 250 au, such that thedust-to-gas mass ratio remains the same inside the gap as in thereference model (R). For the third model, the so-called gas-richgap model, hereafter referred to as (B), we only lowered the dustsurface density by a factor of 10 between 80 and 250 au, suchthat the gas surface density remains the same inside the gap as inthe reference model. This setup leads to a much lower dust-to-gas mass ratio of 0.001 inside the gap in model (B). Our 1 + ff erent radii,which are treated independently; furthermore, the physical con-ditions and abundances at the gap’s edges may not be completelyfeasible. We present the resulting density and thermal disk struc-tures in Fig. 1. The width of the gap used in our modeling is larger than a typical gap width of ∼ −
20 au observed in thedust continuum, but it is similar to the wide gap found in theAS 209 disk. The e ff ects of narrower gaps of a similar depth lo-cated anywhere between 80 and 250 au on the disk physics andchemistry would be similar because of the 1 +
1D nature of ourmodeling tools.
Other than cosmic rays, stellar X-ray radiation is one of the pri-mary ionization mechanisms in protoplanetary disks. X-rays areproduced by the accreting gas and its interaction with magneticfields (see Rab et al. (2018) and references therein). In additionto X-rays, it has recently been recognized that the high-energystellar energetic particles (SEP) could also play an important rolein the disk ionization (Rab et al. 2017). We added the SEP ion-ization using the active T Tauri model of Rab et al. (2017). Whilefar ultraviolet (FUV) radiation dominates ionization in the diskatmosphere (Röllig et al. 2007), both X-rays and SEP are theprimary ionization sources in the intermediate layers of the disk(Semenov et al. 2004). In the midplane, the ionization is dom-inated by the cosmic rays or the decay of short-lived radionu-clides (SLRs).The X-ray ionization rate at a specific disk location primarilydepends on the X-ray flux and the hardness of the X-ray energyspectrum. To compute the X-ray flux at a given point in the disk,one should perform X-ray radiative transfer modeling. The ad-dition of the Compton scattering is necessary for more realisticcalculations of the X-ray ionization rates in the disk. However,the Monte-Carlo X-ray radiative transfer is computationally ex-pensive.Instead, a more conventional approach is the computation ofthe X-ray ray tracing on simplified 1D or 2D geometries and thefurther usage of parametrization to approximate the disk struc-ture. In recent versions of the ProDiMo code, the X-ray radiativetransfer was performed for every disk structure (Rab et al. 2018).In the DALI code, the Compton scattering is neglected, and theone-dimension ray tracing from a point-like stellar source is as-sumed (Bruderer et al. 2009).The parametrization of the disk’s X-ray irradiation givenin Bai & Goodman (2009) is based on the radiative transfermodel of two bremsstrahlung-emitting coronal rings located atthe height of 10 R (cid:12) from the rotation axis and a similar dis-tance above and below the disk midplane. This model also takesthe Compton scattering and photoionization absorption into ac-count. Bai & Goodman (2009) followed the approach of Igea &Glassgold (1999), who studied the X-ray ionization of the inner(on the order of a few au) regions of protoplanetary disks. Igea& Glassgold (1999) found that the ionization structure can beparameterized as a function of only vertical column density fora given radius by scaling down an unattenuated X-ray flux at thedisk atmosphere. We adopted the simplified approach of Bai &Goodman (2009) in our ANDES 1 +
1D iterative modeling of thedisk structure. After the iterative calculations of the disk werefinished, we utilized the more rigorous method of Bruderer et al.(2009) and recalculated the disk ionization structure, assuming asingle-point X-ray source associated with the central star. The X-ray shielding of the gap by the inner disk regions was computedby the 2D ionizing radiation ray tracing without scattering.Finally, Galactic cosmic rays (CRs) dominate the ioniza-tion rate in the disk midplanes when the local gas surface den-sity does not exceed about 100 g cm − or when the stellar windand magnetic mirroring are not that strong (Dolginov & Stepin-ski 1994; Gammie 1996; Cleeves et al. 2013; Rab et al. 2017; Article number, page 3 of 12 & A proofs: manuscript no. 38572corr
Model (R) Model (A) Model (B)
100 200 300 400 500Radius [au]10 − − C o l u m n d e n s i t y [ g c m − ] Gas100 × Dust 100 200 300 400 500Radius [au]10 − − C o l u m n d e n s i t y [ g c m − ] Gas100 × Dust 100 200 300 400 500Radius [au]10 − − C o l u m n d e n s i t y [ g c m − ] Gas100 × Dust10 Radius [au]0.50.00.5 H e i g h t / R a d i u s × Dust densityGas density 10 − − − − − − V o l u m e d e n s i t y [ g c m − ] Radius [au]0.50.00.5 H e i g h t / R a d i u s × Dust densityGas density 10 − − − − − − V o l u m e d e n s i t y [ g c m − ] Radius [au]0.50.00.5 H e i g h t / R a d i u s × Dust densityGas density 10 − − − − − − V o l u m e d e n s i t y [ g c m − ] Radius [au]0.50.00.5 H e i g h t / R a d i u s Dust temperatureGas temperature T = K T = K T e m p e r a t u r e [ K ] Radius [au]0.50.00.5 H e i g h t / R a d i u s Dust temperatureGas temperature T = K T = K T e m p e r a t u r e [ K ] Radius [au]0.50.00.5 H e i g h t / R a d i u s Dust temperatureGas temperature T = K T = K T e m p e r a t u r e [ K ] Fig. 1.
Physical structure of the three protoplanetary disk models. The left panels are for the reference model (R), the middle panels are for thegas-poor gap model (A), and the right panels are for the gas-rich gap model (B). In the upper row, the column densities of gas and dust arepresented. The dust column density is multiplied by 100 for comparison purposes with the gas and the typical dust-to-gas mass ratio of 0.01. Inthe second row, the upper half of each panel is the gas density, and the bottom half is the dust density, multiplied by 100 for comparison. In thebottom row, the upper half of each panel is the gas temperature, and the bottom half is the dust temperature. The location of the gap is marked bya hatched rectangle. The 20 and 40 K gas temperature isotherms are shown as blue and red contour lines, respectively.
Table 2.
Properties of the disk models.
Model Gas gap Gas column density Dust gap Dust-to-gas mass ratio[au] (inside the gap) [au] (inside the gap)(R) Reference disk no 100 g cm − (cid:16) r (cid:17) − no 0.01(A) Gas- and dust-poor gap 80 −
250 0 . × (Reference) 80 −
250 0.01(B) Gas-rich, dust-poor gap no (Reference) 80 −
250 0.001Padovani et al. 2018). Otherwise, in the densest disk regionsshielded from CRs, the decay of SLRs plays a major role. Weused the unattenuated CR ionization rate of 1 . × − s − andthe SLR ionization rate of 6 . × − s − (Semenov et al. 2004).In the studied regions of our adopted disk models, neither the gassurface density nor the stellar wind is strong enough to e ffi cientlyblock the CRs, so the impact of SLRs on the disk midplane ion-ization is negligible.Last but not least, for each iterative modeling step, the radialCO and H column densities were calculated. They were used to compute the UV shielding factors, using an interpolation tablefrom Lee et al. (1996) (see also Semenov & Wiebe 2011; Visseret al. 2009). Next, we used the computed physical structures of the diskand performed time-dependent chemical modeling as post-processing with the ALCHEMIC code (Semenov et al. 2010).The gas-grain network is based on the osu.2007 rate file with
Article number, page 4 of 12rigorii V. Smirnov-Pinchukov et al.: Disk gaps traced by HCO + isotopologues updates as of June 2013 from the kida database (see Wake-lam et al. (2012)). The network is supplied with a set of ap-proximately 1 000 reactions with high-temperature barriers fromHarada et al. (2010) and Harada et al. (2012). This network wasextended by using a statistical branching approach with a setof deuterium fractionation reactions and it includes up to triply-deuterated species (see Albertsson et al. 2013, 2014). Primal iso-tope exchange reactions for H + , CH + , and C H + from Roberts& Millar (2000), Gerlich et al. (2002), Roberts et al. (2004), andRoue ff et al. (2005) were included. Ortho and para forms of H ,H + , and H + isotopologues were considered. The correspondingnuclear spin-state exchange processes were added from experi-mental and theoretical studies (Gerlich 1990; Flower et al. 2004;Walmsley et al. 2004; Flower et al. 2006; Pagani et al. 2009;Hugo et al. 2009; Honvault et al. 2011; Sipilä et al. 2013).The assumed grains are spherical nonporous silicate parti-cles with the material density of 3 g cm − and a radius of 0 . µ m.Each grain provides around 1 . × surface sites (Biham et al.2001). The gas-grain interactions include the freeze-out of neu-tral species and electrons to dust grains with a 100% stickingprobability. In addition, the dissociative recombination and ra-diative neutralization of molecular ions on charged grains andgrain recharging are taken into account. As desorption mecha-nisms for ices, we considered thermal, CR-, UV-driven, and re-active desorption processes. An FUV photodesorption yield of1 × − was adopted. In addition, FUV-driven photodissocia-tion reactions inside ice mantles were implemented as in Gar-rod & Herbst (2006) and Semenov & Wiebe (2011). Surface re-combination is assumed to proceed via the classical Langmuir-Hinshelwood mechanism (e.g., Hasegawa et al. 1992). The ratiobetween di ff usion and desorption energies of surface reactants isassumed to be 0.4. Hydrogen tunneling through rectangular reac-tion barriers with a thickness of 1 Å was computed using Eq. (6)from (Hasegawa et al. 1992). For each act of surface recombina-tion, the e ffi ciency of reactive desorption of the reaction productsdirectly into the gas phase was assumed to be 1% (Garrod et al.2007; Vasyunin & Herbst 2013).Our chemical network consists of 1268 species made of 13elements and 38812 reactions. As initial abundances, we usedthe so-called low metals set of mainly atomic abundances (ex-cept H and HD) from (Lee et al. 1998). Table 3 summarizesthe initial relative abundances (wrt the total amount of hydrogennuclei). The ALCHEMIC code and the chemical network wereused to model the chemical evolution of the disk over an evolu-tionary time span of 1 Myr for all three physical structures thatwere considered.With our large gas-grain chemistry network with deuteriumfractionation, it was not practical to add another full set of Cand O reactions. That is why we opted for simple scaling rela-tions when deriving abundances of C and O-bearing isotopo-logues of CO and HCO + . We assumed that the H CO + abun-dances are equal to 1% of the HCO + abundances, while C Oabundances are equal to 0.2% of the CO abundances (Wil-son & Rood 1994). This simplistic approach does not take se-lective photodissociation of CO isotopologues and isotopic C-exchange reactions into account, which may lead to uncertaintiesin the adopted H CO + and C O abundances by a factor of 3–5(e.g., Visser et al. 2009; Miotello et al. 2014, 2018). Since wemainly focus on analyzing the behavior of abundance and lineratios, rather than absolute values, such a discrepancy should nota ff ect the major underlying trends. http://kida.obs.u-bordeaux1.fr Table 3.
Initial Chemical Abundances.
Species Abundance Species Abundance n (X) / n (H) n (X) / n (H)ortho-H . × − S 9 . × − para-H . × − Si 9 . × − He 9 . × − Fe 2 . × − O 1 . × − Na 2 . × − C 7 . × − Mg 1 . × − N 2 . × − Cl 1 . × − HD 1 . × − P 2 . × − Table 4.
LIME input parameters.
Parameter Internal name ValueMinimal spatial scale minScale 1 auNumber of grid points pIntensity 10 Number of sink points sinkPoints 3000Velocity resolution velres 20 m s − Distance to the source distance 100 pc
To test whether the radial variations in the chemical structurescould become visible in the emission images, we have used thecomputed physical and chemical structures of the disk and per-formed LTE line radiative transfer with LIME v1.9.4 (LIne Mod-eling Engine) by Brinch & Hogerheijde (2010). We used the cor-responding spectroscopic and collisional rate data from the Lei-den Atomic and Molecular Database (LAMDA: Schöier et al.2005). Other key input parameters for LIME are listed in Table4. We tested the LTE line radiative transfer calculations, bothwith and without the dust opacities taken into account, andwe did not notice any significant di ff erences in the continuum-subtracted line emission for the low-J HCO + and CO isotopo-logue lines in the studied low-density regions of the disk. Asthe dust continuum is usually subtracted from the science-readyobservational data and since the LIME computations, which in-clude dust opacities, take much longer, we only show the resultsobtained for the pure gas emission case. For each synthetic im-age, we ran 50 instances of LIME and then averaged them toreduce the noise. As global emission distributions look similarfor the J = =
3. Results
The reference model does not have any density gaps. Its gas sur-face density is set by the the power law (Eq. 1) with Σ =
100 g cm − . This type of monotonous surface density profile stillleads to rather nonmonotonous abundance structures (see Fig. 2).This nonhomogeneous abundance distribution is typical of diskmodels with gas-grain chemistry, which are sensitive to the localvariations in temperature, density, ionization, and high-energyradiation intensities (e.g., Semenov & Wiebe 2011).In our reference disk model, the self-shielded CO extends farinto the disk atmosphere, until the height-to-radius ratio of ∼ . Article number, page 5 of 12 & A proofs: manuscript no. 38572corr
Model (R) Model (A) Model (B) Radius [au]0.50.00.5 H e i g h t / R a d i u s CO iceCO T = K T = K − − N u m b e r d e n s i t y [ c m − ] Radius [au]0.50.00.5 H e i g h t / R a d i u s CO iceCO T = K T = K − − N u m b e r d e n s i t y [ c m − ] Radius [au]0.50.00.5 H e i g h t / R a d i u s CO iceCO T = K T = K − − N u m b e r d e n s i t y [ c m − ] Radius [au]0.50.00.5 H e i g h t / R a d i u s D in H + H + − − − − − N u m b e r d e n s i t y [ c m − ] Radius [au]0.50.00.5 H e i g h t / R a d i u s D in H + H + − − − − − N u m b e r d e n s i t y [ c m − ] Radius [au]0.50.00.5 H e i g h t / R a d i u s D in H + H + − − − − − N u m b e r d e n s i t y [ c m − ] Radius [au]0.50.00.5 H e i g h t / R a d i u s DCO + HCO + − − − − − N u m b e r d e n s i t y [ c m − ] Radius [au]0.50.00.5 H e i g h t / R a d i u s DCO + HCO + − − − − − N u m b e r d e n s i t y [ c m − ] Radius [au]0.50.00.5 H e i g h t / R a d i u s DCO + HCO + − − − − − N u m b e r d e n s i t y [ c m − ] Fig. 2.
Chemical structures of the three protoplanetary disk models. The left panels are for the reference model (R), the middle panels are for thegas-poor gap model (A), and the right panels are for the gas-rich gap model (B). Di ff erent species are shown in the top and bottom halves of eachpanel. D in H + is the total number of D atoms in H + isotopologues. The location of the gap is shown as a hatched background. The gas isothermsof 20 and 40 K are shown in blue and red in the upper panels, respectively, highlighting the CO snowline location and shape. and above (Figs. 2 & 3). The CO snowline is located at about30 au in the midplane. Inside the radius of 20 au, the gas-phaseCO exists all the way through the disk down to the midplane.The CO molecules are e ffi ciently converted to CO in a regionbetween ∼ −
30 au, causing the gas-phase CO abundances todecrease at height-to-radius ratio below ∼ .
2. Outside of 30 au,the gas-phase CO resides in the molecular layer above the heightand radius scales of ∼ . − .
3. The CO ice is mainly concen-trated in the midplane at r (cid:38)
30 au, and the height-to-radius ratiobelow ∼ . + isotopologues sensitively dependson the local ionization rate and the local gas density, since itsprimary production mechanism involves the ionization of H followed by the ion-molecule reactions of H + with H . Conse-quently, in the dense midplane, where ionization rates are low,the abundances of H + isotopologues decrease, but they neverdisappear entirely as they do for the gas-phase CO. A key dif-ference between the main H + isotopologue and its minor D-isotopologues is that the molecular layer of the deuterated H + ions is shifted closer to the cold disk midplane compared to thatof H + (Fig. 3, left panels). This is because unlocking D from HD via deuterium exchange reactions most rapidly proceeds at tem-peratures below ∼ −
100 K and is particularly e ffi cient whenCO is frozen out and when H mainly exists in the para form(e.g., Flower et al. 2006; Albertsson et al. 2013; Sipilä et al.2013; Teague et al. 2015; Huang et al. 2017). This leads to asituation when abundances of, for example, H D + can becomecomparable to the H + abundances in some disk locations, whiletheir vertically-integrated column density would still di ff er by afactor of ∼ −
100 (Fig. 2 and 4, bottom left panels). Panelswith the label "D in H + " show the total number of D atoms inH + isotopologues (H D + , HD + , and D + ).The global distribution of the HCO + isotopologues is deter-mined by the distribution of their parental species, the gaseousCO and the H + isotopologues. In general, the HCO + and DCO + abundances do follow the H + and CO gas-phase distributionsand they are absent in the CO depletion region. However, thereare still major di ff erences between these ions. First, the DCO + molecular layer is, by necessity, co-spatial with the upper part ofthe H + deuterated isotopologue layer, where CO is still not com-pletely frozen out (Fig. 3, top left panel). Second, because ofthat, the DCO + layer is located beneath the HCO + layer and thus Article number, page 6 of 12rigorii V. Smirnov-Pinchukov et al.: Disk gaps traced by HCO + isotopologues Model (R) Model (A) Model (B)
100 200 300 400 500Radius [au]0.00.20.40.60.8 Z / R COD in H + DCO +
100 200 300 400 500Radius [au]0.00.20.40.60.8 Z / R COD in H + DCO +
100 200 300 400 500Radius [au]0.00.20.40.60.8 Z / R COD in H + DCO +
100 200 300 400 500Radius [au]0.00.20.40.60.8 Z / R COH + HCO +
100 200 300 400 500Radius [au]0.00.20.40.60.8 Z / R COH + HCO +
100 200 300 400 500Radius [au]0.00.20.40.60.8 Z / R COH + HCO +
100 200 300 400 500Radius [au]0.00.20.40.60.8 Z / R HCO + DCO +
100 200 300 400 500Radius [au]0.00.20.40.60.8 Z / R HCO + DCO +
100 200 300 400 500Radius [au]0.00.20.40.60.8 Z / R HCO + DCO + Fig. 3.
Vertical mass distribution of H + , deuterium in H + isotopologues, CO, HCO + , and DCO + as a function of the radius. The color-filledstripes show the vertical location of the 25, 50, and 75 mass percentiles (bottom border, median line, and top border, respectively). Half of the totalmolecular gas mass is located within the corresponding stripe. The dotted contour lines show the gas isodensities from 1 × − to 1 × − g cm − with a logarithmic step of 10 . , which is the same as in the top panel of Fig. 1. The left panel is the reference model (R), the middle panel isthe gas-poor gap model (A), and the right panel is the gas-rich gap model (B). The location of the gap is shown as a hatched background. Whenmolecules are co-spatial, their mass distributions overlap (e.g., as in the case of deuterated isotopologues of H + in the model (A)). probes slightly di ff erent physical conditions in the disk (Fig. 3,bottom left panel).As a result of localized variations in the chemical structure,the vertical column densities of the HCO + isotopologues andother key species in the gap-free reference model do show thepresence of weak chemical gaps at outer radii or r >
100 au(Fig. 4). Furthermore, these chemical gaps a ff ect the line excita-tion and appear as emission gaps on the ideal synthetic spectrafor the reference disk model with the monotonous global physi-cal structure (see Fig. 5). The appearance of these chemical gapsis di ff erent for various HCO + isotopologues and C O, and itdepends on their line optical depths and / or the location of theemitting layer. Thus, when taking molecular emission gaps atface value, one could misinterpret the data as being indicative ofreal physical gaps in the disk gas, while it may not be the case. The lower amount of dust and gas leads to higher temperatures,larger pressure scale heights, and hence lower densities in thegas-poor gap model compared to the reference model (Fig. 1,middle panels). The reduced opaqueness of the disk matter inthe gas-poor gap results in a vertical shift of the molecular lay-ers down toward the midplane (Figs. 2 and 3, middle panels).While the vertical column density of H and H + and the total(ice and gas) CO concentration decrease by a factor of ∼
10 in-side the gap, the vertical column density of gaseous CO almostremains intact. The balance between the CO photodissociationand freeze-out is similar in the shifted CO layer inside the gap,just as in the gap-free reference model (Fig. 4, top right panel).Naturally, the total amount of CO ice, which dominates CO con-centration outside the CO snowline, is lower in the gas-poor gapmodel (A) compared to the reference model (R).This vertical shift of the CO molecular layer in the gas-poor gap brings it closer to the molecular layer of the H + D-isotopologues, but moves it away from the main H + molecular Article number, page 7 of 12 & A proofs: manuscript no. 38572corr
100 200 300 400 500Radius [au]10 C o l u m n d e n s i t y [ c m − ] para-H ortho-H RAB 100 200 300 400 500Radius [au]10 C o l u m n d e n s i t y [ c m − ] CO iceCO100 200 300 400 500Radius [au]10 C o l u m n d e n s i t y [ c m − ] HCO + DCO +
100 200 300 400 500Radius [au]10 C o l u m n d e n s i t y [ c m − ] H + D in H + Fig. 4.
Radial profiles of the vertical column densities of the selected molecules. D in H + is the total number of D atoms in H + isotopologues. Thesolid lines correspond to the reference model (R), the dash-dotted lines to the gas-poor model (A), and the dashed lines to the gas-rich model (B).The location of the gap is shown by the gray rectangle. layer (Fig. 3, middle panels). It leads to an increase in the rateof the DCO + -forming reaction, which is proportional to a prod-uct of deuterated isotopologues of H + and CO volume densities.However, it also leads to an opposite e ff ect for the HCO + abun-dances. Consequently, the column density of HCO + decreases bya factor of ∼
3, whereas the column density of DCO + increasesby a similar factor (Fig. 4, bottom left panel).These chemical e ff ects remain visible in the emission maps.Therefore, the gas-poor gap model leads to weaker H CO + emission inside the gap, while the chemically-una ff ected, op-tically thin C O and optically-thick HCO + emission remainrelatively una ff ected. In contrast, the DCO + emission increasesinside the gas-poor gap (Fig. 5). This e ff ect is very similar tothe observed behavior of DCO + , H CO + , and C O in AS 209(Huang et al. 2017; Favre et al. 2019).
Similarly to the gas-poor model (A), the disk model where thegap is only devoid from dust, model (B) shows higher tempera-tures and slightly lower volume densities compared to the refer-ence model (R) (Fig. 1, right panels). The smaller dust columndensity in the gas-rich gap leads to lower extinction of the ion-izing and dissociating radiation and makes molecular freeze-outless e ffi cient. This also shifts the molecular layers closer to the midplane, as in the previous gap model (Figs. 2 and 3, right pan-els).The dust-depleted gap with the same amount of gas as inthe reference model leads to about a 10 times higher columndensity of gaseous CO and a ∼ − + D-isotopologues faster, while a higher concentration of CO in themolecular layer also boosts the production of HCO + . Thus, themodel with the dust-poor gap leads to a similar increase in thecolumn densities of HCO + , DCO + , and CO. This strong increasein the molecular concentrations is imprinted into the correspond-ing increase in the HCO + , DCO + , and C O emission in thedust-poor gap (Fig. 5). Thus, using the single ALMA Band 6or NOEMA Band 3 spectral setups with DCO + and the CO iso-topologue or the two setups with the HCO + , H CO + , and DCO + ions, one could distinguish between gas-poor and dust-poor diskgaps and estimate the overall depletion of gas and / or dust there.
4. Discussion
The version of ANDES used in this study simulates each disk ra-dius independently, thus limiting us in the modeling approachesfor ionizing radiation transport that is consistent with the thermalstructure. The X-rays and stellar energetic particles from the star
Article number, page 8 of 12rigorii V. Smirnov-Pinchukov et al.: Disk gaps traced by HCO + isotopologues
500 0 5005002500250500 O ff c e t , [ a u ] HCO + , (R) 500 0 5005002500250500 DCO + , (R) 500 0 5005002500250500 H CO + , (R) 500 0 5005002500250500 C O, (R)500 0 5005002500250500 O ff c e t , [ a u ] HCO + , (A) 500 0 5005002500250500 DCO + , (A) 500 0 5005002500250500 H CO + , (A) 500 0 5005002500250500 C O, (A)500 0 500Offcet, [au]5002500250500 O ff c e t , [ a u ] HCO + , (B)10 Intensity, K km s
500 0 500Offcet, [au]5002500250500 DCO + , (B)10 Intensity, K km s
500 0 500Offcet, [au]5002500250500 H CO + , (B)10 Intensity, K km s
500 0 500Offcet, [au]5002500250500 C O, (B)10 Intensity, K km s
500 0 5005002500250500 O ff c e t , [ a u ] HCO + (A)/(R) 500 0 5005002500250500 DCO + (A)/(R) 500 0 5005002500250500 H CO + (A)/(R) 500 0 5005002500250500 C O (A)/(R)500 0 500Offcet, [au]5002500250500 O ff c e t , [ a u ] HCO + (B)/(R)10 Intensity ratio 500 0 500Offcet, [au]5002500250500 DCO + (B)/(R)10 Intensity ratio 500 0 500Offcet, [au]5002500250500 H CO + (B)/(R)10 Intensity ratio 500 0 500Offcet, [au]5002500250500 C O (B)/(R)10 Intensity ratio
Fig. 5.
Integrated intensity 0th-moment maps. From left to right: HCO + , DCO + , H CO + , and C O, all J = / (R), and the ratio (B) / (R). The gap boundaries aremarked by the two dotted lines. The optically thick HCO + emission is not significantly a ff ected by the presence of the gas-poor gap. In contrast,the C O and particularly H CO + intensities decrease inside the gas-poor gap, while the DCO + intensity increases inside this gap. cause significant heating in the gas, a ff ecting its vertical structureand ionization fraction. The 2D nature of a gap must have an im-pact on the underlying chemical structure and thus it should alsoa ff ect the gas thermal balance. Consistent 2D modeling with ac- curate modeling of the gas thermal balance will be implementedin the future 2D version of ANDES , first presented in Mol-yarova et al. (2017, 2018), which is being actively developed. Article number, page 9 of 12 & A proofs: manuscript no. 38572corr
From the chemical point of view, the calculated abundancesof simple species in our model su ff er from intrinsic uncertaintiesdue to the reaction rate uncertainties, which are on the order of afactor of 3 for HCO + isotopologues and less than a factor of 2 forthe gas-phase CO (Vasyunin et al. 2008; Albertsson et al. 2013;Iqbal et al. 2018). In addition, our chemical network does notinclude reactions involving isotopologues except for deuterium,and hence we had to rescale the abundances of HCO + and CO toget the H CO + and C O abundances. This simplification mayintroduce additional chemical uncertainty by a factor of 2-3 forboth these minor isotopologues. The main reason for that is thecomputational demands required for our model with deuteriumfractionation. A further increase in the number of species from ∼ > − C- and O-isotopologues wouldincrease the computation time of the chemical evolution by atleast a factor of 10, and we decided to leave the more detailedanalysis for future studies.
Carney et al. (2018) have studied the gap at 40–60 au be-tween the two dusty rings in the protoplanetary disk aroundHD 169142. In contrast to our case, DCO + J = + radial profiles, the behavior ofthe DCO + emission alone is neither consistent with the DCO + observations of AS 209 nor with our modeling. The first reasonfor this inconsistency could be a stronger gas density depletion,estimated as a factor of 40 in their model, in comparison to ourfactor of 10 depletion. In their case, the DCO + / H ratio insidethe gap is lower by a factor of ∼
20 in comparison with the outerdust ring. This is somewhat similar to the increase we have inthe gas-poor model (A). The second major di ff erence is that thisis a Herbig Ae disk (the stellar mass is 1.65 M (cid:12) ). Compared toour cold T Tauri disk with the CO snowline located at ∼
30 au,the CO snowline in the HD 169142 disk is likely located beyond50 au. Thus, the gap at 40–60 au in HD 169142 could be insidethe CO snowline and not outside, as in our model.Teague et al. (2017) have studied the dip in the CS emissionat the radius of r ∼ −
90 au in TW Hya, which coincideswith one of the (sub)mm dust gaps in this disk. They computeda set of models with di ff erent maximal depths of the Gaussiangap ( Σ x / Σ A = . B ) , . C ) , D )). They considered severalcases of the gap, both with and without gas inside. The modelthat qualitatively reproduces the observed CS radial profile, es-pecially its second derivative, was model C with an intermediategap depth. The behavior of the CS emission in TW Hya is similarto that of HCO + in our model. In Fig. 6 we show the CS columndensity derived from our model in comparison with HCO + andDCO + . CS is more sensitive to the lowering of the dust-to-gasratio because it unlocks more carbon from CO by photodisso-ciation and due to less e ffi cient freeze out. In the gas-poor gapmodel (A), there is no significant di ff erence in the CS abundancewith respect to the reference model (R). On the other hand, in thecase of the gas-rich gap model (B), the CS abundance increasesstronger than the abundance of HCO + with respect to the refer-ence model (R). It is also consistent with the results of Teagueet al. (2017) for their gas-rich model Cd . Thus, CS, together with
100 200 300 400 500Radius [au]10 C o l u m n d e n s i t y [ c m − ] HCO + CS RAB100 200 300 400 500Radius [au]10 C o l u m n d e n s i t y [ c m − ] DCO + CS Fig. 6.
Column densities of CS (purple) in comparison with the HCO + (blue) and DCO + (green) column densities. The solid line denotes thereference model (R), the dash-dotted line is the gas-poor model (A),and the dashed line is the gas-rich model (B). The location of the gap isshown by the gray box. DCO + , can also serve as a probe for gas depletion in disk gaps ifH CO + data are not available.Finally, Favre et al. (2019) performed the thermo-chemicalmodeling to fit the observed CO isotopologues emission inAS 209. For the outer gap at ∼
100 au, they modeled the de-pletion of both gas and dust. Unlike this work, they assumed thatthe dust surface density is depleted by a factor of 100, while thegas surface density is depleted by only a factor of 2 . −
10. Ourresults are consistent with their modeling results, as the strongergas depletion reproduces the CO observations in the gap better.While they have shown the ALMA data for the DCO + J = + and C O emission line profiles and the location of the assumedgaps are shown in Fig. 7. The profiles were normalized to havethe same slope in the outer disk. The increase in the DCO + fluxand the decrease in C O in comparison to their outer disk slopeat the location of the gaps are clearly visible, as in our gas-poormodel (A) (Fig. 5, second and fourth rows).However, our modeling does not fit the AS 209 observationsquantitatively. The stellar mass of 0 . M (cid:12) and thus the relativelylow luminosity that we adopted in our models are typical forlow-mass T Tauri stars. However, the type K4-5 AS 209 star ismore luminous and massive and its disk is likely warmer than ourdisk models (Teague et al. 2018). Compared to our disk models, Article number, page 10 of 12rigorii V. Smirnov-Pinchukov et al.: Disk gaps traced by HCO + isotopologues the CO snowline in the AS 209 disk is located at ∼ −
40 au(or about 10 au farther out from the central star). Still, it does notchange the main conclusion of our study because the outer gap inAS 209 is located outside of the CO snowline, as in our models,leading to similar chemical e ff ects for the HCO + isotopologues. Fig. 7.
Azimuthally averaged profiles of DCO + (green) and C O (gray)emission line fluxes from Favre et al. (2019). The profiles were normal-ized to have the same slope outside of the gap (120 - 250 au); the slopeis shown in blue. There is an increase in DCO + emission and a decreasein C O emission intensity at the gap location in comparison with it. Inthe bottom left corner of the image, the beam size of 0.23” is shown.The locations of the gaps are shown by the hatched boxes.
5. Conclusion
We present detailed physical-chemical and line radiative trans-fer simulations of protoplanetary disks, both with and without awide gap. For that, we used the ANDES thermo-chemical modelto calculate 1 +
1D thermal and density disk structures and in-cluded various ionization sources (FUV, X-rays, SEP, CRs, andSLRs). Next, we used the gas-grain chemical code ALCHEMICwith a deuterated network to calculate molecular abundances inthe disk. Finally, we utilized LIME for the line radiative transferto produce the ideal synthetic images in molecular lines of theHCO + isotopologues and C O. Three disk models were consid-ered: the reference gap-free model, the model with the gas-poorgap where both gas and dust surface densities are depleted by afactor of 10, and the gas-rich gap model where solely the dustsurface density is depleted by a factor of 10.Our simulations reveal that in using the HCO + and CO iso-topologues, one could observationally distinguish between vari-ous types of the disk gaps (gas-poor versus gas-rich). The mostconvenient proxy of the gas or dust depletion in the disk gapsis either the DCO + / H CO + and DCO + / HCO + line ratios (us-ing two ALMA Band 6 or NOEMA Band 3 setups) or theDCO + / C O line ratio (using single ALMA Band 6 or NOEMABand 3 spectral setups). Together with high-resolution dust con-tinuum data, these ratios allow one to estimate the degree of thegas or dust depletion in the disk gaps. Namely, if most of theemission lines appear brighter inside the (sub)mm dust gap, es-pecially H CO + , it means that the gas in the gap is not stronglydepleted. In contrast, if C O or H CO + emission decreaseswhile the DCO + emission increases inside the gap, it means thatthe disk gap is substantially depleted in both gas and dust. Fi-nally, the results of our study using the gas-poor disk model are in accordance with the ALMA observations of AS 209, whereC O and DCO + emission inside the outer ∼
100 au gap showsthe opposite behavior. Thus, using intensity ratios of DCO + andC O or the HCO + isotopologues, one might get invaluable in-formation about the disk gap physics and the gap clearing mech-anisms. Acknowledgements.
The authors thank Cecile Favre for sharing the azimuthallyaveraged emission line profiles for AS 209 ALMA observations. The authors ac-knowledge the Python open-source community. This work was done using theopen-source packages NumPy, SciPy, AstroPy, and Matplotlib. D.S. acknowl-edges support by the Deutsche Forschungsgemeinschaft through SPP 1833:“Building a Habitable Earth” (SE 1962 / References
Akimkin, V., Zhukovska, S., Wiebe, D., et al. 2013, ApJ, 766, 8Albertsson, T., Semenov, D., & Henning, T. 2014, ApJ, 784, 39Albertsson, T., Semenov, D. A., Vasyunin, A. I., Henning, T., & Herbst, E. 2013,ApJS, 207, 27Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41Bai, X.-N. & Goodman, J. 2009, ApJ, 701, 737Bergin, E. A., Du, F., Cleeves, L. I., et al. 2016, ApJ, 831, 101Biham, O., Furman, I., Pirronello, V., & Vidali, G. 2001, ApJ, 553, 595Brinch, C. & Hogerheijde, M. R. 2010, A&A, 523, A25Bruderer, S., Doty, S. D., & Benz, A. O. 2009, ApJS, 183, 179Bruderer, S., van der Marel, N., van Dishoeck, E. F., & van Kempen, T. A. 2014,A&A, 562, A26Carney, M. T., Fedele, D., Hogerheijde, M. R., et al. 2018, A&A, 614, A106Cazzoletti, P., van Dishoeck, E. F., Visser, R., Facchini, S., & Bruderer, S. 2018,A&A, 609, A93Ceccarelli, C., Caselli, P., Bockelée-Morvan, D., et al. 2014, in Protostars andPlanets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning,859Cleeves, L. I., Adams, F. C., & Bergin, E. A. 2013, ApJ, 772, 5Cleeves, L. I., Bergin, E. A., & Harries, T. J. 2015, ApJ, 807, 2Cleeves, L. I., Bergin, E. A., Öberg, K. I., et al. 2017, ApJ, 843, L3Dolginov, A. Z. & Stepinski, T. F. 1994, ApJ, 427, 377Draine, B. T. & Lee, H. M. 1984, ApJ, 285, 89Du, F. & Bergin, E. A. 2014, ApJ, 792, 2Dutrey, A., Guilloteau, S., Piétu, V., et al. 2017, A&A, 607, A130Favre, C., Fedele, D., Maud, L., et al. 2019, ApJ, 871, 107Fedele, D., Carney, M., Hogerheijde, M. R., et al. 2017, A&A, 600, A72Fedele, D., Tazzari, M., Booth, R., et al. 2018, A&A, 610, A24Flower, D. R., Pineau des Forêts, G., & Walmsley, C. M. 2004, A&A, 427, 887Flower, D. R., Pineau Des Forêts, G., & Walmsley, C. M. 2006, A&A, 449, 621Gammie, C. F. 1996, ApJ, 457, 355Garrod, R. T. & Herbst, E. 2006, A&A, 457, 927Garrod, R. T., Wakelam, V., & Herbst, E. 2007, A&A, 467, 1103Garufi, A., Podio, L., Codella, C., et al. 2020, A&A, 636, A65Gerlich, D. 1990, J. Chem. Phys., 92, 2377Gerlich, D., Herbst, E., & Roue ff , E. 2002, Planet. Space Sci., 50, 1275Gorti, U. & Hollenbach, D. 2004, ApJ, 613, 424Guzmán, V. V., Huang, J., Andrews, S. M., et al. 2018, ApJ, 869, L48Harada, N., Herbst, E., & Wakelam, V. 2010, Astrophys. J.„ 721, 1570Harada, N., Herbst, E., & Wakelam, V. 2012, ApJ, 756, 104Hasegawa, T. I., Herbst, E., & Leung, C. M. 1992, ApJS, 82, 167Honvault, P., Jorfi, M., González-Lezana, T., Faure, A., & Pagani, L. 2011, Phys-ical Review Letters, 107, 023201Huang, J., Andrews, S. M., Cleeves, L. I., et al. 2018a, ApJ, 852, 122Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018b, ApJ, 869, L42Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2020, ApJ, 891, 48Huang, J. & Öberg, K. I. 2015, ApJ, 809, L26Huang, J., Öberg, K. I., Qi, C., et al. 2017, AJ, 835Hugo, E., Asvany, O., & Schlemmer, S. 2009, J. Chem. Phys., 130, 164302Igea, J. & Glassgold, A. E. 1999, ApJ, 518, 848Iqbal, W., Wakelam, V., & Gratier, P. 2018, A&A, 620, A109Isella, A., Guidi, G., Testi, L., et al. 2016, Phys. Rev. Lett., 117, 251101Keppler, M., Teague, R., Bae, J., et al. 2019, A&A, 625, A118 Article number, page 11 of 12 & A proofs: manuscript no. 38572corr
Lee, H. H., Herbst, E., Pineau des Forets, G., Roue ff , E., & Le Bourlot, J. 1996,A&A, 311, 690Lee, H.-H., Roue ff , E., Pineau des Forets, G., et al. 1998, A&A, 334, 1047Lee, J.-E., Lee, S., Baek, G., et al. 2019, Nature Astronomy, 3, 314Lin, D. N. C. & Papaloizou, J. C. B. 1993, in Protostars and Planets III, ed. E. H.Levy & J. I. Lunine, 749Long, F., Herczeg, G. J., Harsono, D., et al. 2019, ApJ, 882, 49Miotello, A., Bruderer, S., & van Dishoeck, E. F. 2014, A&A, 572, A96Miotello, A., Facchini, S., van Dishoeck, E. F., & Bruderer, S. 2018, A&A, 619,A113Miotello, A., Facchini, S., van Dishoeck, E. F., et al. 2019, A&A, 631, A69Molyarova, T., Akimkin, V., Semenov, D., et al. 2018, ApJ, 866, 46Molyarova, T., Akimkin, V., Semenov, D., et al. 2017, ApJ, 849, 130Padovani, M., Ivlev, A. V., Galli, D., & Caselli, P. 2018, A&A, 614, A111Pagani, L., Vastel, C., Hugo, E., et al. 2009, A&A, 494, 623Parfenov, S. Y., Semenov, D. A., Sobolev, A. M., & Gray, M. D. 2016, MNRAS,460, 2648Pavlyuchenkov, Y., Semenov, D., Henning, T., et al. 2007, ApJ, 669, 1262Qi, C., Öberg, K. I., & Wilner, D. J. 2013, ApJ, 765, 34Rab, C., Güdel, M., Padovani, M., et al. 2017, A&A, 603, 96Rab, C., Güdel, M., Woitke, P., et al. 2018, A&A, 609, A91Roberts, H., Herbst, E., & Millar, T. J. 2004, A&A, 424, 905Roberts, H. & Millar, T. J. 2000, A&A, 361, 388Röllig, M., Abel, N. P., Bell, T., et al. 2007, A&A, 467, 187Roue ff , E., Lis, D. C., van der Tak, F. F. S., Gerin, M., & Goldsmith, P. F. 2005,A&A, 438, 585Salinas, V. N., Hogerheijde, M. R., Murillo, N. M., et al. 2018, A&A, 616, A45Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005,A&A, 432, 369Semenov, D., Hersant, F., Wakelam, V., et al. 2010, A&A, 522, A42Semenov, D. & Wiebe, D. 2011, ApJ, 196, 25Semenov, D., Wiebe, D., & Henning, T. 2004, A&A, 417, 93Sipilä, O., Caselli, P., & Harju, J. 2013, A&A, 554, A92Teague, R., Bae, J., Birnstiel, T., & Bergin, E. A. 2018, ApJ, 868, 113Teague, R. & Loomis, R. 2020, arXiv e-prints, arXiv:2007.11906Teague, R., Semenov, D., Gorti, U., et al. 2017, ApJ, 835, 228Teague, R., Semenov, D., Guilloteau, S., et al. 2015, A&A, 574, A137van Boekel, R., Henning, T., Menu, J., et al. 2017, ApJ, 837, 132van Terwisga, S. E., van Dishoeck, E. F., Cazzoletti, P., et al. 2019, A&A, 623,A150Vasyunin, A. I. & Herbst, E. 2013, ApJ, 769, 34Vasyunin, A. I., Semenov, D., Henning, T., et al. 2008, ApJ, 672, 629Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, Astron. Astrophys, 503, 323Wakelam, V., Herbst, E., Loison, J.-C., et al. 2012, Astrophys. J., Suppl. Ser, 199,21Walmsley, C. M., Flower, D. R., & Pineau des Forêts, G. 2004, A&A, 418, 1035Walsh, C., Nomura, H., & van Dishoeck, E. 2015, A&A, 582, A88Wiebe, D., Semenov, D., & Henning, T. 2003, A&A, 399, 197Wilson, T. L. & Rood, R. 1994, ARA&A, 32, 191Woitke, P., Kamp, I., & Thi, W. F. 2009, A&A, 501, 383Yorke, H. W. & Bodenheimer, P. 2008, in Astronomical Society of the PacificConference Series, Vol. 387, Massive Star Formation: Observations ConfrontTheory, ed. H. Beuther, H. Linz, & T. Henning, 189, E., Lis, D. C., van der Tak, F. F. S., Gerin, M., & Goldsmith, P. F. 2005,A&A, 438, 585Salinas, V. N., Hogerheijde, M. R., Murillo, N. M., et al. 2018, A&A, 616, A45Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005,A&A, 432, 369Semenov, D., Hersant, F., Wakelam, V., et al. 2010, A&A, 522, A42Semenov, D. & Wiebe, D. 2011, ApJ, 196, 25Semenov, D., Wiebe, D., & Henning, T. 2004, A&A, 417, 93Sipilä, O., Caselli, P., & Harju, J. 2013, A&A, 554, A92Teague, R., Bae, J., Birnstiel, T., & Bergin, E. A. 2018, ApJ, 868, 113Teague, R. & Loomis, R. 2020, arXiv e-prints, arXiv:2007.11906Teague, R., Semenov, D., Gorti, U., et al. 2017, ApJ, 835, 228Teague, R., Semenov, D., Guilloteau, S., et al. 2015, A&A, 574, A137van Boekel, R., Henning, T., Menu, J., et al. 2017, ApJ, 837, 132van Terwisga, S. E., van Dishoeck, E. F., Cazzoletti, P., et al. 2019, A&A, 623,A150Vasyunin, A. I. & Herbst, E. 2013, ApJ, 769, 34Vasyunin, A. I., Semenov, D., Henning, T., et al. 2008, ApJ, 672, 629Visser, R., van Dishoeck, E. F., & Black, J. H. 2009, Astron. Astrophys, 503, 323Wakelam, V., Herbst, E., Loison, J.-C., et al. 2012, Astrophys. J., Suppl. Ser, 199,21Walmsley, C. M., Flower, D. R., & Pineau des Forêts, G. 2004, A&A, 418, 1035Walsh, C., Nomura, H., & van Dishoeck, E. 2015, A&A, 582, A88Wiebe, D., Semenov, D., & Henning, T. 2003, A&A, 399, 197Wilson, T. L. & Rood, R. 1994, ARA&A, 32, 191Woitke, P., Kamp, I., & Thi, W. F. 2009, A&A, 501, 383Yorke, H. W. & Bodenheimer, P. 2008, in Astronomical Society of the PacificConference Series, Vol. 387, Massive Star Formation: Observations ConfrontTheory, ed. H. Beuther, H. Linz, & T. Henning, 189