Using HCO$^+$ isotopologues as tracers of gas depletion in protoplanetary disk gaps

The widespread rings and gaps seen in the dust continuum in protoplanetary disks are sometimes accompanied by similar substructures seen in molecular line emission. One example is the outer gap at 100 au in AS 209, which shows that the H$_{13}$CO$^+$ and C$_{18}$O emission intensities decrease along with the continuum in the gap, while the DCO$^+$ emission increases inside the gap. We aim to study the behavior of DCO$^+$/H$_{13}$CO$^+$ and DCO$^+$/HCO$^+$ ratios in protoplanetary disk gaps assuming the two scenarios: the gas depletion follows the dust depletion and only the dust is depleted. We first modeled the physical disk structure using the thermo-chemical model ANDES. This 1+1D steady-state disk model calculates the thermal balance of gas and dust and includes the FUV, X-rays, cosmic rays, and other ionization sources together with the reduced chemical network for molecular coolants. Afterward, this physical structure was adopted for calculations of molecular abundances with the extended gas-grain chemical network with deuterium fractionation. Ideal synthetic spectra and 0th-moment maps were produced with LIME. We are able to qualitatively reproduce the increase in the DCO$^+$ intensity and the decrease in the H$_{13}$CO$^+$ and C$_{18}$O intensities inside the disk gap, which is qualitatively similar to what is observed in the outer AS 209 gap. The corresponding disk model assumes that both the gas and dust are depleted in the gap. The model with the gas-rich gap, where only the dust is depleted, produces emission that is too bright in all HCO$^+$ isotopologues and C$_{18}$O. The DCO$^+$/H$_{13}$CO$^+$ line ratio can be used to probe gas depletion in dust continuum gaps outside of the CO snow line. The DCO$^+$/C$_{18}$O line ratio shows a similar, albeit weaker, effect; however, these species can be observed simultaneously with a single ALMA or NOEMA setup.


Introduction
Protoplanetary disks are believed to be the birthplaces of planetary systems. While more and more theoretical studies of the planet formation in disks appear in the literature, they need observational constraints regarding physical conditions and chemical composition in various disks and disk locations. Now, in the Atacama Large (sub)Millimeter Array (ALMA), NOrthern Extended Millimeter Array (NOEMA), Very Large Telescope / Spectro-Polarimetric High-contrast Exoplanet REsearch (VLT/SPHERE), and Gemini Planet Imager (GPI) era, the spatially resolved observations of dust and gas in protoplanetary disks in nearby star-forming regions have become a routine. High-resolution studies in dust continuum and scattered light with an angular resolution up to 0.025 − 0.04 reveal complex substructures with multiple gaps, inner holes, rings, and spirals (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 and other molecules also show the presence of substructures in the disk gas (see, e.g., Isella et al. 2016;Fedele et al. 2017;Huang et al. 2018a). Often disk substructures that are detected in the scattered light at near-infrared wavelengths do not coincide with the substructures visible in the submillimeter 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 in this disk at the ∼ 0.4 resolution with ALMA. They have found a notable difference between H 13 CO + and DCO + emission radial profiles, with the DCO + emission peak extending farther out up to radii of 80 − 100 au. The follow-up higher resolution observations have shown that the AS 209 disk has a dust gap at around 55-120 au Favre et al. 2019). At this gap location, the optically thin C 18 O J=2-1 emission decreases slightly in intensity, while optically thick 12 CO J=2-1 does not show any intensity decreases, and the DCO + J=3-2 line reaches its peak intensity. Inspired by these puzzling observational findings, we decided to study whether a combination of optically thin and thick HCO + isotopologues could be used as a probe for gas depletion in disk gaps that have been detected before in dust continuum or CO isotopologues.
There are several major factors that can cause molecular line emission in disks to show gap-like or ring-like substructures. Firstly, a circular gap or a ring in the line emission could be driven by a local, rather azimuthally symmetric deviation in the Article number, page 1 of 12 arXiv:2009.09962v1 [astro-ph.EP] 21 Sep 2020 physical structure of the disk, either in density or temperature (e.g., due to planet-disk interactions, snowlines, a change in dust properties, etc. Lin & Papaloizou 1993;Guzmán et al. 2018;Huang et al. 2020). In the case of a surface density drop, the total molecular 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 increase or decrease in the gas's kinetic temperature or background radiation affect molecular line excitation and hence emission intensity, assuming that the line is thermalized.
Secondly, molecular emission could show rings or gaps caused primarily by various chemical effects. The most prominent ones are ultraviolet (UV)-driven production or (selective) photodissociation in the disk atmosphere, C/O variations, and the freeze-out of molecules onto the dust grain surfaces in the disk midplane (see, e.g., Dutrey et al. 2017;Cazzoletti et al. 2018;Miotello et al. 2018Miotello et al. , 2019van Terwisga et al. 2019;Garufi et al. 2020). For example, c-C 3 H 2 , CCH, CN, H 2 CO, and CS emission often shows a ring-like appearance in disks. Whereas for hydrocarbons and CN emission, the C/O variations and high-energy UV or X-ray irradiation are the most important factors; for CS and H 2 CO emission, the surface chemistry processes and freeze-out play the major role (Bergin et al. 2016;Cazzoletti et al. 2018;Miotello et al. 2019;van Terwisga et al. 2019). The ring-like appearance of N 2 D + and DCO + emission in disks is mainly due to low-temperature deuterium fractionation when gaseous CO freezes out in the disk midplane (Qi et al. 2013;Ceccarelli et al. 2014;Huang & Öberg 2015;Salinas et al. 2018;Garufi et al. 2020). In addition, temporal variations in physical conditions can also affect chemical abundances, for example, 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 effect for molecular line excitation and emission intensity is the departure from local thermal equilibrium (LTE). It depends on the molecular properties and the location of the emitting molecular layer in the disk (Pavlyuchenkov et al. 2007). The most notable examples are the CN radical emitting from upper disk layer, where excitation is dominated by radiative processes, and CH 3 OH, which have a complex energy level structure with torsional and coupled energy levels (e.g., Parfenov et al. 2016;Cazzoletti et al. 2018;Teague & Loomis 2020). Another effect that may affect the line excitation in disks locally is the localized deviation in the physical structure of the disk, which does not affect the global structure of the disk otherwise (e.g., local heating in a circumplanetary disk Cleeves et al. 2015). This shows the complex interplay between the physical structure, chemistry, and excitation conditions in defining the strength of line emission and its spatial distribution in a disk.
The paper is structured as follows. In the Section 2, we present the tools we use and the assumptions we make. In Section 3 we present the results for each model (3.1-3.3). In Section 4 we discuss the model uncertainties (4.1) and the comparison with previous observational gas gap studies (4.2). We provide the final conclusion in Section 5.

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 used to simulate 1+1D disk gas and dust thermal structure (Subsection 2.1). Secondly, we explain how we combined the simulated vertical columns into three 2D disk models, with one reference model and two models with and without gas depletion in the dust gap. We recalculated the radial ionization and photodissociation self-shielding factors for those 2D physical structures (2.2). Thirdly, we briefly overview the chemical model ALCHEMIC with the deuterated network, which we used to calculate the chemical abundances (2.3). Finally, we describe how we computed the synthetic molecular line emission maps with LIME (2.4).

Thermo-chemical modeling
The standard practice of chemical modeling of protoplanetary disks includes some density structure setup, dust radiative transfer simulations for the thermal structure of the dust, gas thermal balance, and chemical kinetics calculations using the resulting structure (Dust And LInes (DALI) (Bruderer et al. 2009, 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;Salinas et al. 2018)). The ANDES model allows one to simulate the chemical evolution of a protoplanetary disk with a detailed treatment of gas and dust thermal balance. Its original version (Akimkin et al. 2013) is based on a 1+1D approach, where the vertical disk structure at each radius is simulated independently. In the upper, low-density regions of the protoplanetary disk, the gas temperature is not coupled to the dust temperature. The gas density becomes too low for thermal accommodation with the dust, and the gas usually becomes hotter than the dust as it cools itself rather inefficiently (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 balance for this predefined chemical structure is calculated. The chemical kinetics calculations are iterated with thermal balance until convergence in further iterations. In the utilized version of the ANDES code, the dust temperature step is recalculated again for the hydrostatic solution that fits the vertical temperature structure, assuming well-mixed small dust particles and, therefore, a constant dust-to-gas mass ratio.
The recent version of the ANDES 2D code is based on a fully 2D approach and it simultaneously iterates the vertical disk structure and the chemical evolution, which allows one to accurately simulate some time-dependent effects such as luminosity outbursts in the FU Ori-type systems (Molyarova et al. 2017(Molyarova et al. , 2018. However, this recent version does not feature the detailed gas thermal balance. In this study, we use the original version of ANDES with the reduced chemical network for molecular coolants and modifications to the modeling of the X-ray and cosmic ray ionization rates. The reduced chemical network consists of 63 species and 480 reactions to accurately model the abundances of C, C + , H 2 O, CO, CO 2 , OH, H 2 , and other simple species. It has been obtained from the full network (described in Subsection 2.3) by using the iterative reaction-based reduction method followed by manual tuning (Wiebe et al. 2003;Semenov et al. 2004).
Using the ANDES model, we first computed a set of disk thermo-chemical models to obtain 1+1D temperature and density structures. We assumed the stellar and disk parameters, as listed in Table 1. We selected a T Tauri star with a mass of M = 0.5 M , and we derived its effective temperature T eff = 3820 K and radius R = 1.47 R from the evolutionary track model assuming an age of 2 Myr (Yorke & Bodenheimer 2008, private communication). We did not aim to qualitatively model the K5etype AS 209 star, rather we modeled a generic T Tauri protoplanetary disk. For the UV excess, we assumed the blackbody emis- Dust opacity Draine & Lee (1984) Dust grain size 10 −5 cm Dust-to-gas mass ratio 0.01, 0.001 Grazing angle 0.05 Gas column density at 1 au 10, 100 g cm −2 Density power-law slope −1 sion with T eff = 20000 K, normalized to the total flux of blackbody radiation with T eff = 4000 K. We did not take the evolution of the stellar radius and effective temperature into account. For the dust radiative transfer and surface chemistry, we used nonevolving single-size silicate grains with a size of 10 −5 cm. We set the dust-to-gas mass ratio to 0.01 in our reference disk model. The disk surface density was parameterized using the following power-law distribution, without taking tapering at outer radii into account: (1) This approach leads to a denser outer disk and affects the integral disk properties. On the other hand, it does not affect the effect of the local substructure much. We also separately computed radii between 80 and 250 au with a 10 times lower surface density, and with a 10 times lower dust-to-gas ratio for the gap models, which is described in the next section. The summary of the adopted stellar and disk parameters is presented in Table 1.

Models with gaps
After that, we combined the vertical structures computed by AN-DES for each disk radii into the three protoplanetary disk models (see Table 2). The first model, the so-called reference model, hereafter referred to as (R), is the disk model without any dust or gas gaps. For the other two disk models, we used the separately computed vertical structures for the radii between 80 and 250 au, and those from the reference model (R) outside of this gap. For the second model, the so-called gas-poor gap model, hereafter referred to as (A), we lowered gas and dust surface densities by a factor of 10 between 80 and 250 au, such that the dust-to-gas mass ratio remains the same inside the gap as in the reference model (R). For the third model, the so-called gas-rich gap model, hereafter referred to as (B), we only lowered the dust surface density by a factor of 10 between 80 and 250 au, such that the gas surface density remains the same inside the gap as in the reference model. This setup leads to a much lower dust-togas mass ratio of 0.001 inside the gap in model (B). Our 1+1D approach does not include the interaction between different radii, which are treated independently; furthermore, the physical conditions and abundances at the gap's edges may not be completely feasible. We present the resulting density and thermal disk structures in Fig. 1. The width of the gap used in our modeling is larger than a typical gap width of ∼ 1 − 20 au observed in the dust continuum, but it is similar to the wide gap found in the AS 209 disk. The effects of narrower gaps of a similar depth located anywhere between 80 and 250 au on the disk physics and chemistry would be similar because of the 1+1D nature of our modeling tools.

X-rays, cosmic rays, and stellar wind
Other than cosmic rays, stellar X-ray radiation is one of the primary ionization mechanisms in protoplanetary disks. X-rays are produced by the accreting gas and its interaction with magnetic fields (see Rab et al. (2018) and references therein). In addition to X-rays, it has recently been recognized that the high-energy stellar energetic particles (SEP) could also play an important role in the disk ionization (Rab et al. 2017). We added the SEP ionization using the active T Tauri model of Rab et al. (2017). While far ultraviolet (FUV) radiation dominates ionization in the disk atmosphere (Röllig et al. 2007), both X-rays and SEP are the primary ionization sources in the intermediate layers of the disk (Semenov et al. 2004). In the midplane, the ionization is dominated by the cosmic rays or the decay of short-lived radionuclides (SLRs).
The X-ray ionization rate at a specific disk location primarily depends on the X-ray flux and the hardness of the X-ray energy spectrum. To compute the X-ray flux at a given point in the disk, one should perform X-ray radiative transfer modeling. The addition of the Compton scattering is necessary for more realistic calculations of the X-ray ionization rates in the disk. However, the Monte-Carlo X-ray radiative transfer is computationally expensive.
Instead, a more conventional approach is the computation of the X-ray ray tracing on simplified 1D or 2D geometries and the further usage of parametrization to approximate the disk structure. In recent versions of the ProDiMo code, the X-ray radiative transfer was performed for every disk structure (Rab et al. 2018). In the DALI code, the Compton scattering is neglected, and the one-dimension ray tracing from a point-like stellar source is assumed (Bruderer et al. 2009).
The parametrization of the disk's X-ray irradiation given in Bai & Goodman (2009) is based on the radiative transfer model of two bremsstrahlung-emitting coronal rings located at the height of 10 R from the rotation axis and a similar distance above and below the disk midplane. This model also takes the Compton scattering and photoionization absorption into account. 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 be parameterized as a function of only vertical column density for a given radius by scaling down an unattenuated X-ray flux at the disk atmosphere. We adopted the simplified approach of Bai & Goodman (2009) in our ANDES 1+1D iterative modeling of the disk structure. After the iterative calculations of the disk were finished, we utilized the more rigorous method of Bruderer et al. (2009) and recalculated the disk ionization structure, assuming a single-point X-ray source associated with the central star. The Xray shielding of the gap by the inner disk regions was computed by the 2D ionizing radiation ray tracing without scattering.
Finally, Galactic cosmic rays (CRs) dominate the ionization rate in the disk midplanes when the local gas surface density does not exceed about 100 g cm −2 or when the stellar wind and magnetic mirroring are not that strong (Dolginov & Stepinski 1994;Gammie 1996;Cleeves et al. 2013;Rab et al. 2017;  In the upper row, the column densities of gas and dust are presented. 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. In the 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 the bottom 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 by a hatched rectangle. The 20 and 40 K gas temperature isotherms are shown as blue and red contour lines, respectively.  Padovani et al. 2018). Otherwise, in the densest disk regions shielded from CRs, the decay of SLRs plays a major role. We used the unattenuated CR ionization rate of 1.3 × 10 −17 s −1 and the SLR ionization rate of 6.5 × 10 −19 s −1 (Semenov et al. 2004).
In the studied regions of our adopted disk models, neither the gas surface density nor the stellar wind is strong enough to efficiently block the CRs, so the impact of SLRs on the disk midplane ionization is negligible. Last but not least, for each iterative modeling step, the radial CO and H 2 column densities were calculated. They were used to compute the UV shielding factors, using an interpolation table from Lee et al. (1996) (see also Semenov & Wiebe 2011;Visser et al. 2009).

Chemical model
Next, we used the computed physical structures of the disk and performed time-dependent chemical modeling as postprocessing with the ALCHEMIC code (Semenov et al. 2010). The gas-grain network is based on the osu.2007 rate file with updates as of June 2013 from the kida database 1 (see Wakelam et al. (2012)). The network is supplied with a set of approximately 1 000 reactions with high-temperature barriers from Harada et al. (2010) and Harada et al. (2012). This network was extended by using a statistical branching approach with a set of deuterium fractionation reactions and it includes up to triplydeuterated species (see Albertsson et al. 2013Albertsson et al. , 2014. Primal isotope exchange reactions for H + 3 , CH + 3 , and C 2 H + 2 from Roberts & Millar (2000), Gerlich et al. (2002), Roberts et al. (2004), and Roueff et al. (2005) were included. Ortho and para forms of H 2 , H + 2 , and H + 3 isotopologues were considered. The corresponding nuclear spin-state exchange processes were added from experimental 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 particles with the material density of 3 g cm −3 and a radius of 0.1 µm. Each grain provides around 1.88 × 10 6 surface sites (Biham et al. 2001). The gas-grain interactions include the freeze-out of neutral species and electrons to dust grains with a 100% sticking probability. In addition, the dissociative recombination and radiative neutralization of molecular ions on charged grains and grain recharging are taken into account. As desorption mechanisms for ices, we considered thermal, CR-, UV-driven, and reactive desorption processes. An FUV photodesorption yield of 1 × 10 −5 was adopted. In addition, FUV-driven photodissociation reactions inside ice mantles were implemented as in Garrod & Herbst (2006) and Semenov & Wiebe (2011). Surface recombination is assumed to proceed via the classical Langmuir-Hinshelwood mechanism (e.g., Hasegawa et al. 1992). The ratio between diffusion and desorption energies of surface reactants is assumed to be 0.4. Hydrogen tunneling through rectangular reaction barriers with a thickness of 1 Å was computed using Eq. (6) from (Hasegawa et al. 1992). For each act of surface recombination, the efficiency of reactive desorption of the reaction products directly into the gas phase was assumed to be 1% (Garrod et al. 2007;. Our chemical network consists of 1268 species made of 13 elements and 38812 reactions. As initial abundances, we used the so-called low metals set of mainly atomic abundances (except H 2 and HD) from (Lee et al. 1998). Table 3 summarizes the initial relative abundances (wrt the total amount of hydrogen nuclei). The ALCHEMIC code and the chemical network were used to model the chemical evolution of the disk over an evolutionary time span of 1 Myr for all three physical structures that were considered.
With our large gas-grain chemistry network with deuterium fractionation, it was not practical to add another full set of 13 C and 18 O reactions. That is why we opted for simple scaling relations when deriving abundances of 13 C and 18 O-bearing isotopologues of CO and HCO + . We assumed that the H 13 CO + abundances are equal to 1% of the HCO + abundances, while C 18 O abundances are equal to 0.2% of the 12 CO abundances (Wilson & Rood 1994). This simplistic approach does not take selective photodissociation of CO isotopologues and isotopic 13 Cexchange reactions into account, which may lead to uncertainties in the adopted H 13 CO + and C 18 O abundances by a factor of 3-5 (e.g., Visser et al. 2009;Miotello et al. 2014Miotello et al. , 2018. Since we mainly focus on analyzing the behavior of abundance and line ratios, rather than absolute values, such a discrepancy should not affect the major underlying trends. 1 http://kida.obs.u-bordeaux1.fr

Line radiative transfer
To test whether the radial variations in the chemical structures could become visible in the emission images, we have used the computed physical and chemical structures of the disk and performed LTE line radiative transfer with LIME v1.9.4 (LIne Modeling Engine) by Brinch & Hogerheijde (2010). We used the corresponding spectroscopic and collisional rate data from the Leiden Atomic and Molecular Database (LAMDA: Schöier et al. 2005). Other key input parameters for LIME are listed in Table  4. We tested the LTE line radiative transfer calculations, both with and without the dust opacities taken into account, and we did not notice any significant differences in the continuumsubtracted line emission for the low-J HCO + and CO isotopologue lines in the studied low-density regions of the disk. As the dust continuum is usually subtracted from the science-ready observational data and since the LIME computations, which include dust opacities, take much longer, we only show the results obtained for the pure gas emission case. For each synthetic image, we ran 50 instances of LIME and then averaged them to reduce the noise. As global emission distributions look similar for the J=5-4, 4-3, and 3-2 transitions, we focus on the J=3-2 results.

Reference model (R)
The reference model does not have any density gaps. Its gas surface density is set by the the power law (Eq. 1) with Σ 1au = 100 g cm −1 . This type of monotonous surface density profile still leads to rather nonmonotonous abundance structures (see Fig. 2). This nonhomogeneous abundance distribution is typical of disk models with gas-grain chemistry, which are sensitive to the local variations in temperature, density, ionization, and high-energy radiation intensities (e.g., Semenov & Wiebe 2011).
In our reference disk model, the self-shielded CO extends far into the disk atmosphere, until the height-to-radius ratio of ∼ 0.5 Article number, page 5 of 12

D in H +
3 is the total number of D atoms in H + 3 isotopologues. The location of the gap is shown as a hatched background. The gas isotherms of 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 about 30 au in the midplane. Inside the radius of 20 au, the gas-phase CO exists all the way through the disk down to the midplane. The CO molecules are efficiently converted to CO 2 in a region between ∼ 20 − 30 au, causing the gas-phase CO abundances to decrease at height-to-radius ratio below ∼ 0.2. Outside of 30 au, the gas-phase CO resides in the molecular layer above the height and radius scales of ∼ 0.1 − 0.3. The CO ice is mainly concentrated in the midplane at r 30 au, and the height-to-radius ratio below ∼ 0.3. The distribution of the H + 3 isotopologues sensitively depends on the local ionization rate and the local gas density, since its primary production mechanism involves the ionization of H 2 followed by the ion-molecule reactions of H + 2 with H 2 . Consequently, in the dense midplane, where ionization rates are low, the abundances of H + 3 isotopologues decrease, but they never disappear entirely as they do for the gas-phase CO. A key difference between the main H + 3 isotopologue and its minor Disotopologues is that the molecular layer of the deuterated H + 3 ions is shifted closer to the cold disk midplane compared to that of H + 3 (Fig. 3, left panels). This is because unlocking D from HD via deuterium exchange reactions most rapidly proceeds at temperatures below ∼ 20 − 100 K and is particularly efficient when CO is frozen out and when H 2 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 a situation when abundances of, for example, H 2 D + can become comparable to the H + 3 abundances in some disk locations, while their vertically-integrated column density would still differ by a factor of ∼ 10 − 100 ( Fig. 2 and 4, bottom left panels). Panels with the label "D in H + 3 " show the total number of D atoms in H + 3 isotopologues (H 2 D + , HD + 2 , and D + 3 ).
The global distribution of the HCO + isotopologues is determined by the distribution of their parental species, the gaseous CO and the H + 3 isotopologues. In general, the HCO + and DCO + abundances do follow the H + 3 and CO gas-phase distributions and they are absent in the CO depletion region. However, there are still major differences between these ions. First, the DCO + molecular layer is, by necessity, co-spatial with the upper part of the H +  Fig. 3. Vertical mass distribution of H + 3 , deuterium in H + 3 isotopologues, CO, HCO + , and DCO + as a function of the radius. The color-filled stripes show the vertical location of the 25, 50, and 75 mass percentiles (bottom border, median line, and top border, respectively). Half of the total molecular gas mass is located within the corresponding stripe. The dotted contour lines show the gas isodensities from 1 × 10 −18 to 1 × 10 −14 g cm −2 with a logarithmic step of 10 0.5 , which is the same as in the top panel of Fig. 1. The left panel is the reference model (R), the middle panel is the 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. When molecules are co-spatial, their mass distributions overlap (e.g., as in the case of deuterated isotopologues of H + 3 in the model (A)).
probes slightly different 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 and other key species in the gap-free reference model do show the presence of weak chemical gaps at outer radii or r > 100 au (Fig. 4). Furthermore, these chemical gaps affect the line excitation and appear as emission gaps on the ideal synthetic spectra for the reference disk model with the monotonous global physical structure (see Fig. 5). The appearance of these chemical gaps is different for various HCO + isotopologues and C 18 O, and it depends on their line optical depths and/or the location of the emitting layer. Thus, when taking molecular emission gaps at face value, one could misinterpret the data as being indicative of real physical gaps in the disk gas, while it may not be the case.

Gas-poor model (A)
The lower amount of dust and gas leads to higher temperatures, larger pressure scale heights, and hence lower densities in the gas-poor gap model compared to the reference model (Fig. 1,  middle panels). The reduced opaqueness of the disk matter in the gas-poor gap results in a vertical shift of the molecular layers down toward the midplane (Figs. 2 and 3, middle panels). While the vertical column density of H 2 and H + 3 and the total (ice and gas) CO concentration decrease by a factor of ∼ 10 inside the gap, the vertical column density of gaseous CO almost remains intact. The balance between the CO photodissociation and 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 concentration outside the CO snowline, is lower in the gas-poor gap model (A) compared to the reference model (R).
This vertical shift of the CO molecular layer in the gaspoor gap brings it closer to the molecular layer of the H + 3 Disotopologues, but moves it away from the main H + 3 molecular  Fig. 4. Radial profiles of the vertical column densities of the selected molecules. D in H + 3 is the total number of D atoms in H + 3 isotopologues. The solid 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 rate of the DCO + -forming reaction, which is proportional to a product of deuterated isotopologues of H + 3 and CO volume densities. However, it also leads to an opposite effect for the HCO + abundances. Consequently, the column density of HCO + decreases by a factor of ∼ 3, whereas the column density of DCO + increases by a similar factor (Fig. 4, bottom left panel).
These chemical effects remain visible in the emission maps. Therefore, the gas-poor gap model leads to weaker H 13 CO + emission inside the gap, while the chemically-unaffected, optically thin C 18 O and optically-thick HCO + emission remain relatively unaffected. In contrast, the DCO + emission increases inside the gas-poor gap (Fig. 5). This effect is very similar to the observed behavior of DCO + , H 13 CO + , and C 18 O in AS 209 (Huang et al. 2017;Favre et al. 2019).

Gas-rich, dust-poor model (B)
Similarly to the gas-poor model (A), the disk model where the gap is only devoid from dust, model (B) shows higher temperatures and slightly lower volume densities compared to the reference model (R) (Fig. 1, right panels). The smaller dust column density in the gas-rich gap leads to lower extinction of the ionizing and dissociating radiation and makes molecular freeze-out less efficient. This also shifts the molecular layers closer to the midplane, as in the previous gap model (Figs. 2 and 3, right panels).
The dust-depleted gap with the same amount of gas as in the reference model leads to about a 10 times higher column density of gaseous CO and a ∼ 2 − 3 times lower column density of the CO ice (Fig. 4). The higher concentration of gaseous CO in the cold midplane region makes reactions with the H + The gap boundaries are marked by the two dotted lines. The optically thick HCO + emission is not significantly affected by the presence of the gas-poor gap. In contrast, the C 18 O and particularly H 13 CO + intensities decrease inside the gas-poor gap, while the DCO + intensity increases inside this gap. cause significant heating in the gas, affecting its vertical structure and ionization fraction. The 2D nature of a gap must have an impact on the underlying chemical structure and thus it should also affect the gas thermal balance. Consistent 2D modeling with ac-curate modeling of the gas thermal balance will be implemented in the future 2D version of ANDES 2D , first presented in Molyarova et al. (2017Molyarova et al. ( , 2018, which is being actively developed. From the chemical point of view, the calculated abundances of simple species in our model suffer from intrinsic uncertainties due to the reaction rate uncertainties, which are on the order of a factor of 3 for HCO + isotopologues and less than a factor of 2 for the gas-phase CO (Vasyunin et al. 2008;Albertsson et al. 2013;Iqbal et al. 2018). In addition, our chemical network does not include reactions involving isotopologues except for deuterium, and hence we had to rescale the abundances of HCO + and CO to get the H 13 CO + and C 18 O abundances. This simplification may introduce additional chemical uncertainty by a factor of 2-3 for both these minor isotopologues. The main reason for that is the computational demands required for our model with deuterium fractionation. A further increase in the number of species from ∼ 1250 in the current version to > 2000−5000 species in the network that also includes key 13 C-and 18 O-isotopologues would increase the computation time of the chemical evolution by at least a factor of 10, and we decided to leave the more detailed analysis for future studies. Carney et al. (2018) have studied the gap at 40-60 au between the two dusty rings in the protoplanetary disk around HD 169142. In contrast to our case, DCO + J=3-2 emission intensity decreases in this gap. The thermo-chemical structure of the HD 169142 disk was simulated with the DALI code . To reproduce the observed behavior, Bruderer et al. (2014) assumed the low gas densities in both the gap and the inner dust ring, similar to our gas-poor model (A). While the authors do not show the HCO + radial profiles, the behavior of the DCO + emission alone is neither consistent with the DCO + observations of AS 209 nor with our modeling. The first reason for this inconsistency could be a stronger gas density depletion, estimated as a factor of 40 in their model, in comparison to our factor of 10 depletion. In their case, the DCO + / H 2 ratio inside the gap is lower by a factor of ∼ 20 in comparison with the outer dust ring. This is somewhat similar to the increase we have in the gas-poor model (A). The second major difference is that this is a Herbig Ae disk (the stellar mass is 1.65 M ). Compared to our cold T Tauri disk with the CO snowline located at ∼ 30 au, the CO snowline in the HD 169142 disk is likely located beyond 50 au. Thus, the gap at 40-60 au in HD 169142 could be inside the CO snowline and not outside, as in our model.  have studied the dip in the CS emission at the radius of r ∼ 80 − 90 au in TW Hya, which coincides with one of the (sub)mm dust gaps in this disk. They computed a set of models with different maximal depths of the Gaussian gap (Σ x /Σ A = 0.3(B), 0.55(C), 1(D)). They considered several cases of the gap, both with and without gas inside. The model that qualitatively reproduces the observed CS radial profile, especially its second derivative, was model C with an intermediate gap depth. The behavior of the CS emission in TW Hya is similar to that of HCO + in our model. In Fig. 6 we show the CS column density derived from our model in comparison with HCO + and DCO + . CS is more sensitive to the lowering of the dust-to-gas ratio because it unlocks more carbon from CO by photodissociation and due to less efficient freeze out. In the gas-poor gap model (A), there is no significant difference in the CS abundance with respect to the reference model (R). On the other hand, in the case of the gas-rich gap model (B), the CS abundance increases stronger than the abundance of HCO + with respect to the reference model (R). It is also consistent with the results of    DCO + , can also serve as a probe for gas depletion in disk gaps if H 13 CO + data are not available. Finally, Favre et al. (2019) performed the thermo-chemical modeling to fit the observed CO isotopologues emission in AS 209. For the outer gap at ∼ 100 au, they modeled the depletion of both gas and dust. Unlike this work, they assumed that the dust surface density is depleted by a factor of 100, while the gas surface density is depleted by only a factor of 2.5 − 10. Our results are consistent with their modeling results, as the stronger gas depletion reproduces the CO observations in the gap better. While they have shown the ALMA data for the DCO + J=3-2 emission, the modeling involving deuterated chemistry has been left for the forthcoming paper. The azimuthally averaged DCO + and C 18 O emission line profiles and the location of the assumed gaps are shown in Fig. 7. The profiles were normalized to have the same slope in the outer disk. The increase in the DCO + flux and the decrease in C 18 O in comparison to their outer disk slope at the location of the gaps are clearly visible, as in our gas-poor model (A) (Fig. 5, second and fourth rows).

Comparison with previous studies
However, our modeling does not fit the AS 209 observations quantitatively. The stellar mass of 0.5 M and thus the relatively low luminosity that we adopted in our models are typical for low-mass T Tauri stars. However, the type K4-5 AS 209 star is more luminous and massive and its disk is likely warmer than our disk models (Teague et al. 2018). Compared to our disk models, the CO snowline in the AS 209 disk is located at ∼ 30 − 40 au (or about 10 au farther out from the central star). Still, it does not change the main conclusion of our study because the outer gap in AS 209 is located outside of the CO snowline, as in our models, leading to similar chemical effects for the HCO + isotopologues.

Conclusion
We present detailed physical-chemical and line radiative transfer simulations of protoplanetary disks, both with and without a wide gap. For that, we used the ANDES thermo-chemical model to calculate 1+1D thermal and density disk structures and included various ionization sources (FUV, X-rays, SEP, CRs, and SLRs). Next, we used the gas-grain chemical code ALCHEMIC with a deuterated network to calculate molecular abundances in the disk. Finally, we utilized LIME for the line radiative transfer to produce the ideal synthetic images in molecular lines of the HCO + isotopologues and C 18 O. Three disk models were considered: the reference gap-free model, the model with the gas-poor gap where both gas and dust surface densities are depleted by a factor of 10, and the gas-rich gap model where solely the dust surface density is depleted by a factor of 10.
Our simulations reveal that in using the HCO + and CO isotopologues, one could observationally distinguish between various types of the disk gaps (gas-poor versus gas-rich). The most convenient proxy of the gas or dust depletion in the disk gaps is either the DCO + /H 13 CO + and DCO + /HCO + line ratios (using two ALMA Band 6 or NOEMA Band 3 setups) or the DCO + /C 18 O line ratio (using single ALMA Band 6 or NOEMA Band 3 spectral setups). Together with high-resolution dust continuum data, these ratios allow one to estimate the degree of the gas or dust depletion in the disk gaps. Namely, if most of the emission lines appear brighter inside the (sub)mm dust gap, especially H 13 CO + , it means that the gas in the gap is not strongly depleted. In contrast, if C 18 O or H 13 CO + emission decreases while the DCO + emission increases inside the gap, it means that the disk gap is substantially depleted in both gas and dust. Finally, the results of our study using the gas-poor disk model are in accordance with the ALMA observations of AS 209, where C 18 O and DCO + emission inside the outer ∼ 100 au gap shows the opposite behavior. Thus, using intensity ratios of DCO + and C 18 O or the HCO + isotopologues, one might get invaluable information about the disk gap physics and the gap clearing mechanisms.