# Molecules with ALMA at Planet-forming Scales (MAPS). X. Studying Deuteration at High Angular Resolution toward Protoplanetary Disks

Gianni Cataldi, Yoshihide Yamato, Yuri Aikawa, Jennifer B. Bergner, Kenji Furuya, Viviana V. Guzmán, Jane Huang, Ryan A. Loomis, Chunhua Qi, Sean M. Andrews, et al.

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Gianni Cataldi, Yoshihide Yamato, Yuri Aikawa, Jennifer B. Bergner, Kenji Furuya, et al. Molecules with ALMA at Planet-forming Scales (MAPS). X. Studying Deuteration at High Angular Resolution toward Protoplanetary Disks. The Astrophysical Journal Supplement Series, 2021, 257, 10.3847/1538-4365/ac143d. insu-03672370

HAL Id: insu-03672370

https://insu.hal.science/insu-03672370v1

Submitted on 19 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Molecules with ALMA at Planet-forming Scales (MAPS). X. Studying Deuteration at High Angular Resolution toward Protoplanetary Disks

Gianni Cataldi<sup>1,2</sup>, Yoshihide Yamato<sup>2</sup>, Yuri Aikawa<sup>2</sup>, Jennifer B. Bergner<sup>3,22</sup>, Kenji Furuya<sup>1</sup>, Viviana V. Guzmán<sup>4</sup>, Jane Huang<sup>5,6,22</sup>, Ryan A. Loomis<sup>7</sup>, Chunhua Qi<sup>6</sup>, Sean M. Andrews<sup>6</sup>, Edwin A. Bergin<sup>5</sup>, Alice S. Booth<sup>8,9</sup>, Arthur D. Bosman<sup>5</sup>, L. Ilsedore Cleeves<sup>10</sup>, Ian Czekala<sup>11,12,13,14,15,22</sup>, John D. Ilee<sup>9</sup>, Charles J. Law<sup>6</sup>, Romane Le Gal<sup>6,16,17,18</sup>, Yao Liu<sup>19</sup>, Feng Long<sup>6</sup>, François Ménard<sup>16</sup>, Hideko Nomura<sup>1</sup>, Karin I. Öberg<sup>6</sup>, Kamber R. Schwarz<sup>20,22</sup>, Richard Teague<sup>6</sup>, Takashi Tsukagoshi<sup>1</sup>, Catherine Walsh<sup>9</sup>, David J. Wilner<sup>6</sup>, and Ke Zhang<sup>5,21,23</sup>

<sup>1</sup> National Astronomical Observatory of Japan, Osawa 2-21-1, Mitaka, Tokyo 181-8588, Japan; gianni.cataldi@nao.ac.jp

<sup>2</sup> Department of Astronomy, Graduate School of Science, The University of Tokyo, Tokyo 113-0033, Japan

<sup>3</sup> University of Chicago, Department of the Geophysical Sciences, Chicago, IL 60637, USA

<sup>4</sup> Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile

<sup>5</sup> Department of Astronomy, University of Michigan, 323 West Hall, 1085 S. University Ave., Ann Arbor, MI 48109, USA

<sup>6</sup> Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA

<sup>7</sup> National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA 22903, USA

<sup>8</sup> Leiden Observatory, Leiden University, 2300 RA Leiden, The Netherlands

<sup>9</sup> School of Physics and Astronomy, University of Leeds, Leeds LS2 9JT, UK

<sup>10</sup> Department of Astronomy, University of Virginia, 530 McCormick Rd., Charlottesville, VA 22904, USA

<sup>11</sup> Department of Astronomy and Astrophysics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA

<sup>12</sup> Center for Exoplanets and Habitable Worlds, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA

<sup>13</sup> Center for Astrostatistics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA

<sup>14</sup> Institute for Computational & Data Sciences, The Pennsylvania State University, University Park, PA 16802, USA

<sup>15</sup> Department of Astronomy, 501 Campbell Hall, University of California, Berkeley, CA 94720-3411, USA

<sup>16</sup> Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

<sup>17</sup> IRAP, Université de Toulouse, CNRS, CNES, UT3, 31400 Toulouse, France

<sup>18</sup> IRAM, 300 rue de la piscine, F-38406 Saint-Martin d'Hères, France

<sup>19</sup> Purple Mountain Observatory & Key Laboratory for Radio Astronomy, Chinese Academy of Sciences, Nanjing 210023, People's Republic of China

<sup>20</sup> Lunar and Planetary Laboratory, University of Arizona, 1629 E. University Blvd., Tucson, AZ 85721, USA

<sup>21</sup> Department of Astronomy, University of Wisconsin-Madison, 475 N. Charter St., Madison, WI 53706, USA

*Received 2021 February 23; revised 2021 July 2; accepted 2021 July 7; published 2021 November 3*

## Abstract

Deuterium fractionation is dependent on various physical and chemical parameters. Thus, the formation location and thermal history of material in the solar system is often studied by measuring its D/H ratio. This requires knowledge about the deuteration processes operating during the planet formation era. We aim to study these processes by radially resolving the DCN/HCN (at 0.′3 resolution) and N₂D⁺/N₂H⁺ (~0.′3–0.′9) column density ratios toward the five protoplanetary disks observed by the Molecules with ALMA at Planet-forming scales (MAPS) Large Program. DCN is detected in all five sources, with one newly reported detection. N₂D⁺ is detected in four sources, two of which are newly reported detections. We derive column density profiles that allow us to study the spatial variation of the DCN/ HCN and N₂D⁺/N₂H⁺ ratios at high resolution. DCN/HCN varies considerably for different parts of the disks, ranging from 10⁻³ to 10⁻¹. In particular, the inner-disk regions generally show significantly lower HCN deuteration compared with the outer disk. In addition, our analysis confirms that two deuterium fractionation channels are active, which can alter the D/H ratio within the pool of organic molecules. N₂D⁺ is found in the cold outer regions beyond ~50 au, with N₂D⁺/N₂H⁺ ranging between 10⁻² and 1 across the disk sample. This is consistent with the theoretical expectation that N₂H⁺ deuteration proceeds via the low-temperature channel only. This paper is part of the MAPS special issue of the Astrophysical Journal Supplement.

*Unified Astronomy Thesaurus concepts: Astrochemistry (75); Isotopic abundances (867); Interferometry (808); Protoplanetary disks (1300); Planet formation (1241); Planetary system formation (1257); Millimeter astronomy (1061); Aperture synthesis (53)*

## 1. Introduction

It is clear that the physical and chemical properties of a protoplanetary disk heavily influence the properties and evolution of the bodies (planets, asteroids, and comets) emerging from it. In particular, the composition of nascent planets is set by the composition of the disk at the location of formation (e.g., Öberg & Bergin 2021). Varying physical parameters such as temperature,

gas density, or UV exposure result in a radially and vertically varying composition of the disk (e.g., Henning & Semenov 2013). Spatially resolving the chemical composition of protoplanetary disks is thus a cornerstone of planet formation studies. In this paper, we consider the radial distribution of two deuterated molecules: DCN and N₂D⁺. We are particularly interested in the radial variation of the deuteration fraction: DCN/HCN and N₂D⁺/N₂H⁺. This will allow us to connect deuterium chemistry in disks to measurements of deuterium fractionation in solar system bodies.

<footnote n="22"><sup>22</sup> NASA Hubble Fellowship Program Sagan Fellow.</footnote>

<footnote n="23"><sup>23</sup> NASA Hubble Fellow.</footnote>
In the interstellar medium (ISM), the elemental D/H ratio (by number) is of the order of $1.5 \times 10^{-5}$ (e.g., Hébrard et al. 2005; Linsky et al. 2006, and references therein), consistent with predictions of big bang nucleosynthesis (e.g., Cyburt et al. 2016). However, the deuteration fraction, that is, the ratio of a deuterated molecule to its nondeuterated isotopologue, can exceed the elemental ratio by several orders of magnitude. This phenomenon is known as deuterium fractionation and depends on chemical and physical parameters. For example, the degree of fractionation is generally enhanced at low temperatures. As a consequence, the deuteration fraction carries information about the physical properties of the environment in which the molecules were formed and is often used to infer the formation location and thermal history of material in the ISM or the solar system. For example, comparison of the deuteration fractions of the water in Earth's oceans and in solar system asteroids or comets is used to study the possibility of asteroidal or cometary water delivery to Earth during terrestrial planet formation (e.g., Alexander 2017; O'Brien et al. 2018; Lis et al. 2019).

As in molecular clouds, deuteration in protoplanetary disks starts from a number of exchange reactions involving the main reservoir of deuterium: HD (e.g., Millar et al. 1989; Turner 2001). At low temperature ($\lesssim 30$ K), the most relevant reaction is

$$HD + H_3^+ \rightleftharpoons H_2 + H_2D^+. \quad (1)$$

$H_2D^+$ can then propagate the D to other molecules such as $N_2D^+$ and DCN by proton transfer and subsequent recombination: $N_2 + H_2D^+ \rightarrow N_2D^+ + H_2$, HNC + $H_2D^+ \rightarrow DCNH^+ + H_2$, and DCNH$^+$ + $e^-$ $\rightarrow$ DCN + H. Reaction (1) is exothermic in the forward direction ($\Delta E \approx 230$ K). Therefore, at gas temperatures below $\sim 30$ K, the inverse reaction is suppressed and significant fractionation can occur. In addition, at these low temperatures, molecules that can destroy $H_2D^+$ are frozen out, further enhancing the fractionation (e.g., Roberts & Millar 2000). It is also important to consider the spin state of $H_2$. Ortho-$H_2$ has a higher internal energy than para-$H_2$, meaning that it can more easily drive the reverse reaction. Thus, a higher ortho-to-para ratio means less efficient deuteration (e.g., Willacy et al. 2015). Compared to molecular clouds, the ortho/para ratio of $H_2$ is more easily thermalized by ion-molecule reactions and grain-surface conversions in protoplanetary disks (Aikawa et al. 2015; Furuya et al. 2019). The strong thermal gradients in protoplanetary disks then likely lead to a gradient of the ortho-para ratio, with a higher abundance of ortho-$H_2$ in the warm inner disk.

Reaction (1) is often referred to as the low temperature deuteration channel. At temperatures $\gtrsim 30$ K, the most relevant fractionation reactions are instead expected to be (Millar et al. 1989)

$$HD + CH_3^+ \rightleftharpoons H_2 + CH_2D^+, \quad (2)$$

$$HD + C_2H_2^+ \rightleftharpoons H_2 + C_2HD^+. \quad (3)$$

These reactions are more exothermic than reaction (1), that is, $\Delta E \approx 500$ K (e.g., Roberts & Millar 2000; Roueff et al. 2013; Nyman & Yu 2019). Thus, they could dominate the fractionation for temperatures $\gtrsim 30$ K and are often called the high-temperature deuteration channel. DCN, for example, can be formed by the reaction of an N atom with CHD, which is formed by the dissociative recombination of $CH_2D^+$. The different deuteration pathways are therefore expected to operate in different parts of

the disk (Aikawa et al. 2018): in the midplane of the cold outer disk, fractionation should mainly be initiated by the low-temperature channel, while in the warmer inner region and the disk atmosphere, the high-temperature channel should dominate.

While HCN is mainly destroyed by photodissociation and various ion-molecule reactions (e.g., Aikawa et al. 1999; Willacy & Langer 2000), $N_2H^+$ is destroyed by CO and is thus expected to be abundant in cold regions beyond the CO snow line (e.g., Qi et al. 2013). Therefore, $N_2H^+$ deuteration should proceed via the low-temperature channel only (Millar et al. 1989). On the other hand, both the high- and low-temperature channels are expected to contribute to DCN formation (e.g., Huang et al. 2017; Salinas et al. 2017; Aikawa et al. 2018). This picture can be observationally tested by radially resolving the distribution of DCN and $N_2D^+$ (e.g., Huang et al. 2017; Salinas et al. 2017; Öberg et al. 2021). $DCO^+$ can be formed by both pathways and thus is also used to study the two deuteration pathways (e.g., Öberg et al. 2012, 2021; Salinas et al. 2017; Carney et al. 2018). In this work, we build on these previous efforts by using high-resolution DCN and $N_2D^+$ data obtained by the Molecules with ALMA at Planet-forming Scales (MAPS) Large Program$^{24}$ (Öberg et al. 2021). MAPS targeted five protoplanetary disks (Öberg et al. 2021): three disks around T Tauri stars (IM Lup, GM Aur, AS 209) and two disks around Herbig Ae stars (HD 163296 and MWC 480). These systems cover a range of stellar masses from $1.1 M_\odot$ for IM Lup and GM Aur (Teague et al. 2021) to $2.1 M_\odot$ for MWC 480 (Simon et al. 2019), a range of stellar luminosities from $1.2 L_\odot$ for GM Aur (Macías et al. 2018) to $21.9 L_\odot$ for MWC 480 (Montesinos et al. 2009), and a range of ages from $\sim 1$ Myr for IM Lup (Mawet et al. 2012) and AS 209 (Andrews et al. 2018) to $\sim 7$ Myr for MWC 480 (Simon et al. 2000; Montesinos et al. 2009). All disks have dust substructures in the form of rings and gaps (Andrews et al. 2018; Huang et al. 2018, 2020; Long et al. 2018). GM Aur is the only disk with a central dust cavity. The $^{12}$CO 2–1 disk size, which encloses 90% of the flux (Law et al. 2021a), is largest for IM Lup ($\sim 480$ au) and smallest for AS 209 ($\sim 200$ au). Low signal-to-noise ratio (S/N) $^{12}$CO 2–1 emission often extends to considerably larger radii (e.g., out to $\sim 800–900$ au for IM Lup, Law et al. 2021a). A more detailed description of the MAPS targets is given by Öberg et al. (2021). The MAPS data allow us to study the deuteration at an unprecedented spatial resolution and sensitivity. We aim to compare radial deuteration profiles to model predictions and to study the temperature dependence of the deuteration fraction.

Besides studying deuteration, we aim to investigate the relation of the CO snow line with $N_2D^+$. As mentioned above, the abundance of both $N_2H^+$ and $N_2D^+$ is expected to anticorrelate with the CO abundance, because $N_2H^+$ destruction is enhanced in the presence of CO (e.g., Bergin et al. 2001; Qi et al. 2013). This was used by Qi et al. (2013, 2015, 2019) to estimate the radius of the CO snow line by matching it to the inner edge of the observed $N_2H^+$ emission. But theoretical models show that, besides the midplane, $N_2H^+$ can be moderately abundant in the warm molecular layer as well (e.g., van't Hoff et al. 2017). This can complicate the inference of the CO snow line. On the other hand, $N_2D^+$ is expected to mainly trace the midplane (Aikawa et al. 2018). Therefore, we will test whether $N_2D^+$ can be used as an alternative to $N_2H^+$ to trace the CO snow line.

<footnote>$^{24}$ https://www.alma-maps.info</footnote>
<table><caption>**Table 1**

Overview of Properties of HCN, DCN, and H$^{13}$CN Data Cubes Provided by MAPS and Used in This Study</caption><thead><tr><th>Line</th><th>Rms<sup>a</sup><br>(mJy beam<sup>-1</sup>)</th><th>JvM</th><th>Spectral Resolution<br>(km s<sup>-1</sup>)</th><th>Channel Width<br>(km s<sup>-1</sup>)</th><th>Beam Size</th><th>Usage<sup>c</sup></th></tr></thead><tbody><tr><td colspan="7" style="text-align: center;">IM Lup</td></tr><tr><td>HCN 1-0</td><td>0.7</td><td>0.72</td><td>0.24</td><td>0.5</td><td>0"30 × 0"30</td><td rowspan="5">MOM0, PROF<br>FLUX, AASPEC, COLDENS</td></tr><tr><td>HCN 3-2</td><td>0.3</td><td>0.25</td><td>0.16</td><td>0.2</td><td>0"15 × 0"15</td></tr><tr><td>HCN 3-2</td><td>0.6</td><td>0.49</td><td>0.16</td><td>0.2</td><td>0"30 × 0"30</td></tr><tr><td>DCN 3-2</td><td>1.1</td><td>0.79</td><td>0.20</td><td>0.2</td><td>0"30 × 0"30</td></tr><tr><td>H<sup>13</sup>CN 1-0</td><td>0.7</td><td>0.73</td><td>0.49</td><td>0.5</td><td>0"30 × 0"30</td></tr><tr><td colspan="7" style="text-align: center;">GM Aur</td></tr><tr><td>HCN 1-0</td><td>1.8</td><td>1.00</td><td>0.24</td><td>0.5</td><td>0"31 × 0"31</td><td rowspan="5">MOM0, PROF<br>FLUX, AASPEC, COLDENS</td></tr><tr><td>HCN 3-2</td><td>0.8</td><td>0.57</td><td>0.16</td><td>0.2</td><td>0"15 × 0"15</td></tr><tr><td>HCN 3-2</td><td>1.0</td><td>0.74</td><td>0.16</td><td>0.2</td><td>0"30 × 0"30</td></tr><tr><td>DCN 3-2</td><td>0.7</td><td>0.63</td><td>0.20</td><td>0.2</td><td>0"30 × 0"30</td></tr><tr><td>H<sup>13</sup>CN 1-0</td><td>1.5</td><td>1.00</td><td>0.49</td><td>0.5</td><td>0"33 × 0"33</td></tr><tr><td colspan="7" style="text-align: center;">AS 209</td></tr><tr><td>HCN 1-0</td><td>0.8</td><td>0.83</td><td>0.24</td><td>0.5</td><td>0"30 × 0"29</td><td rowspan="5">MOM0, PROF<br>FLUX, AASPEC, COLDENS</td></tr><tr><td>HCN 3-2</td><td>0.5</td><td>0.28</td><td>0.16</td><td>0.2</td><td>0"15 × 0"15</td></tr><tr><td>HCN 3-2</td><td>0.8</td><td>0.51</td><td>0.16</td><td>0.2</td><td>0"30 × 0"30</td></tr><tr><td>DCN 3-2</td><td>0.7</td><td>0.58</td><td>0.20</td><td>0.2</td><td>0"30 × 0"30</td></tr><tr><td>H<sup>13</sup>CN 1-0</td><td>0.8</td><td>0.90</td><td>0.49</td><td>0.5</td><td>0"30 × 0"30</td></tr><tr><td colspan="7" style="text-align: center;">HD 163296</td></tr><tr><td>HCN 1-0</td><td>0.8</td><td>0.93</td><td>0.24</td><td>0.5</td><td>0"30 × 0"30</td><td rowspan="5">MOM0, PROF<br>FLUX, AASPEC, COLDENS</td></tr><tr><td>HCN 3-2</td><td>0.4</td><td>0.30</td><td>0.16</td><td>0.2</td><td>0"15 × 0"15</td></tr><tr><td>HCN 3-2</td><td>0.7</td><td>0.53</td><td>0.16</td><td>0.2</td><td>0"30 × 0"30</td></tr><tr><td>DCN 3-2</td><td>0.9</td><td>0.72</td><td>0.20</td><td>0.2</td><td>0"30 × 0"30</td></tr><tr><td>H<sup>13</sup>CN 1-0</td><td>0.7</td><td>0.95</td><td>0.49</td><td>0.5</td><td>0"30 × 0"30</td></tr><tr><td colspan="7" style="text-align: center;">MWC 480</td></tr><tr><td>HCN 1-0</td><td>1.7</td><td>0.99</td><td>0.24</td><td>0.5</td><td>0"31 × 0"31</td><td rowspan="5">MOM0, PROF<br>FLUX, AASPEC, COLDENS</td></tr><tr><td>HCN 3-2</td><td>0.8</td><td>0.54</td><td>0.16</td><td>0.2</td><td>0"15 × 0"15</td></tr><tr><td>HCN 3-2</td><td>1.0</td><td>0.74</td><td>0.16</td><td>0.2</td><td>0"30 × 0"30</td></tr><tr><td>DCN 3-2</td><td>0.8</td><td>0.73</td><td>0.20</td><td>0.2</td><td>0"30 × 0"30</td></tr><tr><td>H<sup>13</sup>CN 1-0</td><td>1.5</td><td>1.00</td><td>0.49</td><td>0.5</td><td>0"31 × 0"31</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Measured in an emission-free region of the non-primary-beam-corrected data cube.

<sup>b</sup> Ratio of the CLEAN beam area to the dirty beam area used in the JvM correction. See Appendix A and Czekala et al. (2021) for details.

<sup>c</sup> For HCN 3-2, the following flags indicate the usage: FLUX: disk-integrated flux; MOM0: zeroth moment; PROF: radial emission profile; AASPEC: azimuthally averaged spectra; COLDENS: column density calculation. For the moment 0 maps, the radial emission profiles, and the disk-integrated fluxes, we used the non-primary-beam-corrected images, while primary-beam-corrected images are used for all other analyses.

In summary, we observed DCN and N$_{2}$D$^{+}$ with high sensitivity and angular resolution to study the deuteration chemistry of five protoplanetary disks. In Section 2, we give an overview of the data. In Section 3, we discuss the emission morphology inferred from zeroth moment maps and radial emission profiles. We then proceed to derive radial column density profiles and deuteration profiles in Section 4. We discuss the implications of our results in Section 5 and summarize our conclusions in Section 6.

## 2. Observations

### 2.1. MAPS Data: HCN, H$^{13}$CN, DCN, N$_{2}$D$^{+}$

The MAPS Large Program (project code 2018.1.01055.L) used four spectral setups: two in ALMA Band 3 at a wavelength of ~3 mm and two in ALMA Band 6 at ~1.3 mm (Öberg et al. 2021). Two array configurations were used: a short baseline configuration and a long baseline configuration, resulting in baselines between 15 and 3638 m. The details of the data calibration are described in Öberg et al. (2021).

#### 2.1.1. CLEANing and JvM Correction

As described in detail in Czekala et al. (2021), the MAPS collaboration produced images from the calibrated visibilities using the CLEAN deconvolution algorithm (Högbom 1974) implemented in the CASA task *tclean*, where the “multiscale” version of the algorithm was used (*deconvolver* = ‘multiscale’). Furthermore, MAPS employed a flux scale correction to the CLEANed images that was first discussed by Jorsater & van Moorsel (1995, hereafter JvM correction). A brief summary of the correction is given in Appendix A, while the full details can be found in Czekala et al. (2021).

#### 2.1.2. Adopted MAPS Imaging Products

In this paper, we use the following lines observed by MAPS: the *J* = 3–2 transitions of HCN, DCN, and N$_{2}$D$^{+}$ in Band 6 and the *J* = 1–0 transitions of HCN and H$^{13}$CN in Band 3. For each emission line, images with several different beam sizes were produced by MAPS. For our Band 3 data, we use the MAPS fiducial images with a circular 0"/3 beam. For the Band 6 data, we use images tapered to the same circular 0"/3 beam
<table><caption>**Table 2**

Overview of Properties and Usage of N₂H⁺ and N₂D⁺ Data Cubes</caption><thead><tr><th>Line</th><th>Rms<sup>a</sup><br>(mJy beam<sup>-1</sup>)</th><th>JvM ϵ<sup>b</sup></th><th>δv<sup>c</sup><br>(km s<sup>-1</sup>)</th><th>Channel Width<br>(km s<sup>-1</sup>)</th><th>Beam Size</th><th>Beam PA<sup>d</sup><br>(deg)</th><th>Data Source<sup>e</sup></th><th>Usage<sup>f</sup></th></tr></thead><tbody><tr><td colspan="9" style="text-align:center; font-weight:bold;">IM Lup</td></tr><tr><td>N₂H⁺ 3–2</td><td>3.6</td><td>0.71</td><td>0.15</td><td>0.15</td><td>0″42 × 0″35</td><td>76</td><td>Qi19</td><td>FLUX, MOM0, PROF</td></tr><tr><td>N₂D⁺ 3–2</td><td>1.6</td><td>0.78</td><td>0.09</td><td>0.2</td><td>0″3 × 0″3</td><td>N/A</td><td>Öberg21</td><td>AASPEC, COLDENS</td></tr><tr><td>N₂D⁺ 3–2</td><td>1.8</td><td>0.78</td><td>0.09</td><td>0.2</td><td>0″42 × 0″35</td><td>76</td><td>Öberg21</td><td>AASPEC, COLDENS</td></tr><tr><td colspan="9" style="text-align:center; font-weight:bold;">GM Aur</td></tr><tr><td>N₂H⁺ 3–2</td><td>2.4</td><td>0.75</td><td>0.15</td><td>0.15</td><td>0″31 × 0″20</td><td>-2</td><td>Qi19</td><td>FLUX, MOM0, PROF</td></tr><tr><td>N₂H⁺ 3–2</td><td>2.4</td><td>0.75</td><td>0.15</td><td>0.15</td><td>0″32 × 0″32</td><td>N/A</td><td>Qi19</td><td>AASPEC, COLDENS</td></tr><tr><td>N₂D⁺ 3–2</td><td>1.2</td><td>0.73</td><td>0.09</td><td>0.2</td><td>0″3 × 0″3</td><td>N/A</td><td>Öberg21</td><td>FLUX, MOM0, PROF</td></tr><tr><td>N₂D⁺ 3–2</td><td>1.3</td><td>0.73</td><td>0.09</td><td>0.2</td><td>0″32 × 0″32</td><td>N/A</td><td>Öberg21</td><td>AASPEC, COLDENS</td></tr><tr><td colspan="9" style="text-align:center; font-weight:bold;">AS 209</td></tr><tr><td>N₂H⁺ 3–2</td><td>2.1</td><td>0.67</td><td>0.15</td><td>0.15</td><td>0″34 × 0″24</td><td>-74</td><td>Qi19</td><td>FLUX, MOM0, PROF</td></tr><tr><td>N₂H⁺ 3–2</td><td>2.1</td><td>0.67</td><td>0.15</td><td>0.15</td><td>0″35 × 0″35</td><td>N/A</td><td>Qi19</td><td>AASPEC, COLDENS</td></tr><tr><td>N₂D⁺ 3–2</td><td>1.1</td><td>0.58</td><td>0.09</td><td>0.2</td><td>0″3 × 0″3</td><td>N/A</td><td>Öberg21</td><td>FLUX, MOM0, PROF</td></tr><tr><td>N₂D⁺ 3–2</td><td>1.2</td><td>0.58</td><td>0.09</td><td>0.2</td><td>0″35 × 0″35</td><td>N/A</td><td>Öberg21</td><td>AASPEC, COLDENS</td></tr><tr><td colspan="9" style="text-align:center; font-weight:bold;">HD 163296</td></tr><tr><td>N₂H⁺ 3–2</td><td>2.8</td><td>0.84</td><td>0.26</td><td>0.2</td><td>0″45 × 0″35</td><td>-88</td><td>Qi15</td><td>FLUX, MOM0, PROF</td></tr><tr><td>N₂D⁺ 3–2</td><td>1.2</td><td>0.69</td><td>0.09</td><td>0.2</td><td>0″3 × 0″3</td><td>N/A</td><td>Öberg21</td><td>AASPEC, COLDENS</td></tr><tr><td>N₂D⁺ 3–2</td><td>1.4</td><td>0.69</td><td>0.09</td><td>0.2</td><td>0″45 × 0″35</td><td>-88</td><td>Öberg21</td><td>AASPEC, COLDENS</td></tr><tr><td colspan="9" style="text-align:center; font-weight:bold;">MWC 480</td></tr><tr><td>N₂H⁺ 3–2</td><td>1.8</td><td>0.95</td><td>1.21</td><td>1.05</td><td>0″93 × 0″48</td><td>-29</td><td>Loomis20</td><td>FLUX, MOM0, PROF</td></tr><tr><td>N₂D⁺ 3–2</td><td>1.3</td><td>0.74</td><td>0.09</td><td>0.2</td><td>0″3 × 0″3</td><td>N/A</td><td>Öberg21</td><td>AASPEC, COLDENS</td></tr><tr><td>N₂D⁺ 3–2</td><td>2.2</td><td>0.74</td><td>0.09</td><td>0.2</td><td>0″93 × 0″48</td><td>-29</td><td>Öberg21</td><td>AASPEC, COLDENS</td></tr></tbody></table>

**Notes.**

a Measured in an emission-free region of the non-primary-beam-corrected data cube.

b Ratio of the CLEAN beam area to the dirty beam area used in the JvM correction. See Appendix A and Czekala et al. (2021) for details.

c Spectral resolution.

d Position angle of the beam major axis, measured north through east. For circular beams, position angle is undefined.

e References: Qi19: Qi et al. (2019), Öberg21: Öberg et al. (2021), Qi15: Qi et al. (2015), Loomis20: Loomis et al. (2020).

f For lines with several image cubes, the following flags indicate their usage: FLUX: disk-integrated flux; MOM0: zeroth moment; PROF: radial emission profile; AASPEC: azimuthally averaged spectra; COLDENS: column density calculation. For the moment 0 maps, the radial emission profiles, and the disk-integrated fluxes, we used the non-primary-beam-corrected images, while primary-beam-corrected images are used for all other analyses.

(Czekala et al. 2021) for the weak DCN and N₂D⁺ lines to improve the sensitivity. For the strong HCN 3–2 line in Band 6, we use both the images with 0″15 and 0″3 circular beams, as listed in Table 1. In particular, when fitting for the column densities, we use the image cubes with 0″3 beams in order to have uniform spatial resolution for all the emission lines included in the fit. Since MAPS provided individual HCN and H¹³CN image cubes with spectral axes centered onto individual hyperfine components (Öberg et al. 2021), we used the CASA tasks regrid and imageconcat to combine these cubes into a single cube. The basic properties of the MAPS data cubes we used are listed in Tables 1 (for HCN, DCN, and H¹³CN) and 2 (for N₂D⁺).

## 2.2. N₂H⁺ Archival Data

No transitions of N₂H⁺ were observed by MAPS. Fortunately, the J = 3–2 transition of this molecule has been targeted by previous observations in ALMA Band 7 (λ ~ 1 mm) for all five of our targets. For IM Lup, GM Aur, and AS 209, we use data presented in Qi et al. (2019, project code 2015.1.00678.S, ~0″3–0″4 resolution). For HD 163296, we use the data by Qi et al. (2015, project code 2012.1.00681.S, ~0″5 resolution). Finally, for MWC 480, we use data presented in Loomis et al. (2020, project code 2015.1.00657.S, ~1″ resolution). The

calibration of the data is described in the original papers. Using CASA 5.6, we imaged the visibilities in the same way as for the MAPS data: first, Keplerian masks (i.e., masks that select the regions of the image cube where emission is expected based on the Keplerian rotation of the disk) were constructed following the procedure described in Czekala et al. (2021), taking into account the hyperfine components listed in Table 7. Using the Keplerian masks, the data were imaged with the multiscale CLEAN deconvolver implemented in the CASA task tclean with scales of [0, 5, 15, 25] pixels and a Briggs parameter of 0.5. Finally, to yield the correct flux, the JvM correction was applied to all images following Czekala et al. (2021). For most sources, this process results in N₂H⁺ image cubes with fairly circular beams. The exception is MWC 480, where the beam is 0″93 × 0″48. However, the beam sizes of the resulting N₂H⁺ images do not match the beam sizes of the MAPS N₂D⁺ images. Thus, we produced additional N₂H⁺ and N₂D⁺ data cubes with matching beams using the CASA task imsmooth. These additional cubes are used for the derivation of azimuthally averaged spectra and column densities. For IM Lup, HD 163296, and MWC 480, we only had to smooth the N₂D⁺ data. For AS 209 and GM Aur, both the N₂H⁺ and N₂D⁺ data required smoothing to achieve matching beam sizes. Table 2 lists the properties and the usage of the various N₂H⁺ and N₂D⁺ image cubes.
<table><caption>**Table 3**
Disk-integrated Fluxes (Derived from Nonprimary-beam-corrected Images)</caption><thead><tr><th rowspan="2"></th><th colspan="3">IM Lup</th><th colspan="3">GM Aur</th><th colspan="3">AS 209</th><th colspan="3">HD 163296</th><th colspan="3">MWC 480</th></tr><tr><th>r<sub>min</sub><br>(au)</th><th>r<sub>max</sub><br>(au)</th><th>Flux<br>(mJy km s<sup>-1</sup>)</th><th>r<sub>min</sub><br>(au)</th><th>r<sub>max</sub><br>(au)</th><th>Flux<br>(mJy km s<sup>-1</sup>)</th><th>r<sub>min</sub><br>(au)</th><th>r<sub>max</sub><br>(au)</th><th>Flux<br>(mJy km s<sup>-1</sup>)</th><th>r<sub>min</sub><br>(au)</th><th>r<sub>max</sub><br>(au)</th><th>Flux<br>(mJy km s<sup>-1</sup>)</th><th>r<sub>min</sub><br>(au)</th><th>r<sub>max</sub><br>(au)</th><th>Flux<br>(mJy km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>HCN 1–0</td><td>0</td><td>600</td><td>240 ± 27</td><td>0</td><td>400</td><td>125 ± 16</td><td>0</td><td>250</td><td>198 ± 21</td><td>0</td><td>500</td><td>658 ± 68</td><td>0</td><td>300</td><td>155 ± 18</td></tr><tr><td>HCN 3–2</td><td>0</td><td>600</td><td>2409 ± 242</td><td>0</td><td>400</td><td>1792 ± 180</td><td>0</td><td>300</td><td>2965 ± 297</td><td>0</td><td>500</td><td>7346 ± 736</td><td>0</td><td>300</td><td>2469 ± 248</td></tr><tr><td>DCN 3–2</td><td>0</td><td>500</td><td>88 ± 13</td><td>0</td><td>400</td><td>37 ± 6</td><td>0</td><td>200</td><td>235 ± 24</td><td>0</td><td>400</td><td>115 ± 15</td><td>0</td><td>400</td><td>74 ± 10</td></tr><tr><td>H<sup>13</sup>CN 1–0</td><td>0</td><td>600</td><td>&lt;16</td><td>0</td><td>400</td><td>&lt;36</td><td>0</td><td>250</td><td>&lt;13</td><td>0</td><td>500</td><td>&lt;15</td><td>0</td><td>300</td><td>&lt;20</td></tr><tr><td>N<sub>2</sub>H<sup>+</sup> 3–2</td><td>0</td><td>500</td><td>1483 ± 150</td><td>0</td><td>400</td><td>1025 ± 103</td><td>0</td><td>300</td><td>698 ± 71</td><td>0</td><td>350</td><td>450 ± 51</td><td>0</td><td>350</td><td>279 ± 31</td></tr><tr><td>N<sub>2</sub>D<sup>+</sup> 3–2</td><td>0</td><td>400</td><td>96 ± 14</td><td>0</td><td>400</td><td>24 ± 7</td><td>0</td><td>250</td><td>74 ± 9</td><td>50</td><td>250</td><td>135 ± 15</td><td>60</td><td>250</td><td>28 ± 6</td></tr></tbody></table>

**Notes.** Upper limits are at 3σ significance.

a Minimum radius of the Keplerian mask.

b Maximum radius of the Keplerian mask.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Gallery of hybrid zeroth moment maps for HCN 1–0, HCN 3–2, DCN 3–2, and H$^{13}$CN 1–0, generated by combining a Keplerian mask and a smoothed σ-clip mask (Law et al. 2021a). The use of a σ-clip mask means that some emission is inevitably lost, except for a threshold of 0 σ. Thus, the flux scale can be unreliable, and these maps should be used for presentational purposes only. The clip values employed are indicated in the upper right of each map. The color scales employ either linear or arcsinh stretches, with the lower end saturating at 0 mJy beam$^{-1}$ km s$^{-1}$. The beam is shown by the white ellipse. Note that the noise is not uniform over these maps due to the use of masks. Text in the lower right of the panels marks lines not or only tentatively detected in total flux or a matched uv-plane filter.</figcaption></figure>

in HCN 3–2, while the outer DCN ring corresponds to a shoulder in the HCN 3–2 profile. For GM Aur, both HCN and DCN show a bright inner ring and a faint outer ring. For AS 209, both HCN and DCN are in a ring, while for HD 163296, both species show a double ring. MWC 480 is a notable exception in that HCN and DCN show different morphologies: while HCN is centrally peaked, DCN is in a ring.

There are also a few associations of the HCN and DCN line emission structure with dust substructures (Law et al. 2021a). For example, the inner HCN and DCN rings in IM Lup are associated with a continuum ring-gap structure at ~125 au. The most prominent feature is the association of the outer DCN rings in IM Lup and GM Aur with the edge of the dust continuum disk, as discussed in Section 5.3. For a detailed analysis of the relationships between the radial structures seen

in the MAPS data for HCN, DCN, as well as other molecules and the dust, see Law et al. (2021a).

### 3.3. Radial Emission Profiles of $N_2H^+$ and $N_2D^+$

Figure 4 shows the radial emission profiles of $N_2H^+$ 3–2 and $N_2D^+$ 3–2. These profiles are also produced by azimuthally averaging a zeroth moment map generated with a Keplerian mask only (i.e., without a σ-clipping mask). The $N_2H^+$ 3–2 radial emission profiles of IM Lup, GM Aur, AS 209, and MWC 480 are consistent with the ones presented in Qi et al. (2019) and Loomis et al. (2020).

$N_2H^+$ and $N_2D^+$ show similar ring emission structures, although their detailed structure differs from source to source. Qi et al. (2019) proposed that the morphology of the $N_2H^+$ emission reflects the vertical temperature structure of the disk. A narrow ring with extended tenuous emission (GM Aur) is
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Gallery of hybrid zeroth moment maps for N₂H⁺ 3–2 and N₂D⁺ 3–2, generated by combining a Keplerian mask and a smoothed σ-clip mask (Law et al. 2021a). The use of a σ-clip mask means that some emission is inevitably lost, except for a threshold of 0σ. Thus, the flux scale can be unreliable, and these maps should be used for presentational purposes only. The clip values employed are indicated in the upper right of each map. The color scales saturate at 0 mJy beam⁻¹ km s⁻¹ at the lower end and employ a linear stretch, except for N₂H⁺ toward GM Aur (arcsinh stretch). The ellipses in the lower left indicate the beam size. N₂D⁺ 3–2 toward GM Aur is only tentatively detected in total flux as well as with a matched filter analysis in the uv plane.</figcaption></figure>

expected for a disk that is vertically isothermal up to substantial heights ($z/r \approx 0.2$) above the midplane. Conversely, disks with a vertical temperature gradient above the midplane are expected to present broad rings (IM Lup, AS 209). For HD 163296 and MWC 480, observations with higher S/N and, in the case of MWC 480, higher angular resolution would be useful to firmly distinguish between these two cases.

For the N₂D⁺ emission, the four sources with a clear detection show ring structures. Most interestingly, IM Lup clearly shows a double ring structure peaking at ~100 and ~330 au, similar to the rings detected in DCN (Section 3.2) and DCO⁺ (Öberg et al. 2015). N₂H⁺ in IM Lup also shows a subtle shoulder at the location of the outer ring of N₂D⁺. The relation between the structures seen for DCN, N₂D⁺, and DCO⁺ will be discussed more in detail in Section 5.3.

We note the presence of negative N₂D⁺ emission toward the center of the HD 163296 and MWC 480 disks. This is probably due to the difficulty of achieving a precise continuum subtraction for this line, as discussed in Appendix B.

### 3.4. Azimuthally Averaged Spectra

In order to compute radial column density profiles, we will model the azimuthally averaged spectra of equally spaced, deprojected annuli (radial bins). However, at each spatial pixel of a data cube, the spectrum is shifted with respect to the systemic velocity due to the Keplerian rotation of the gas. Therefore, in the azimuthally averaged spectrum, the emission is spread over a broad range of velocities. This results in a suboptimal S/N. Furthermore, the hyperfine structure of HCN and N₂H⁺ is not resolved in such a broad spectrum and cannot be used to constrain the column density. Therefore, we apply the following procedure prior to azimuthally averaging (Teague et al. 2016; Yen et al. 2016; Matrà et al. 2017): we shift each spectrum by the projected Keplerian velocity,

which is given by

$$v_{\text{Kep}}(r, \theta) = \sqrt{\frac{GM_{\star}}{r}} \sin i \cos \theta, \quad (4)$$

where r and θ are the (deprojected) radial and azimuthal coordinates of the disk, M$_{\star}$ is the (dynamical) stellar mass, and i is the inclination (see Öberg et al. 2021, for the adopted values). Here we are assuming that emission originates from the midplane. The center of the disk is assigned to the proper motion-corrected stellar position. This procedure centers each spectrum at the systemic velocity. When averaging azimuthally, the S/N is increased and the hyperfine structure remains spectrally resolved.

The calculation of the error bars of the averaged spectra is described in Appendix D.1. The size of the radial bins was chosen to be equal to half of the mean of the beam major and minor axes.

The west side of the AS 209 disk is known to be affected by foreground cloud contamination (Öberg et al. 2011; Huang et al. 2016; Guzmán et al. 2018), which should be most relevant for Band 3 data. Thus, for the 1–0 transitions of HCN and H¹³CN, we follow Teague et al. (2018a) and extract all spectra from a ±55° wedge that encompasses the uncontaminated eastern side of the disk. The other disks are not affected by cloud contamination and thus are averaged over the full azimuth.

Figures 5 and 6 show examples of extracted spectra. The full gallery of spectra can be found in Appendix D.2. In Figure 5, the hyperfine structure of HCN 1–0 and 3–2 is readily visible, except in the innermost region, where the finite spatial resolution causes considerable broadening of the spectra. Hyperfine structure is also seen for N₂H⁺ in Figure 6.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Radial emission profiles of HCN 1–0, HCN 3–2, DCN 3–2, and H$^{13}$CN 1–0. The horizontal dashed line marks the zero flux level. The shaded area shows the ±1σ error. If the profile has been scaled by a constant factor, the scaling is indicated in the upper right of the panel. The beam major axis is shown as a horizontal line in the upper right. H$^{13}$CN 1–0 is marked as undetected or tentatively detected based on the total flux measurement and the matched filter analysis.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Radial emission profiles of N$_2$H$^+$ 3–2 and N$_2$D$^+$ 3–2. The horizontal dashed line indicates the zero flux level. The shaded area shows the ±1σ error. If the profile has been scaled by a constant factor, the scaling is indicated in the upper right of the panel. The beam major axis is shown as horizontal black line in the upper right of each panel. N$_2$D$^+$ 3–2 toward GM Aur is marked as tentatively detected based on the measurement of the disk-integrated flux and the matched filter analysis.</figcaption></figure>

# 4. Analysis

## 4.1. Radial Column Density Profiles of HCN and DCN

To derive the radial column density profiles of HCN and DCN, we use the 3–2 transitions of HCN and DCN covered

in Band 6 and the 1–0 transitions of HCN and H$^{13}$CN covered in Band 3. Our modeling includes the hyperfine satellite lines, which help to derive reliable column densities even if the main components are optically thick. Table 6 provides an overview of the hyperfine lines used in the
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Examples of HCN, DCN, and H$^{13}$CN model spectra fit to azimuthally averaged spectra of AS 209 for a few radial bins. Spectra are centered on the systemic velocity. The colored lines show the data, with the shaded regions corresponding to the $1\sigma$ uncertainty. The black curves show 50 randomly selected models drawn from the Markov Chain Monte Carlo (MCMC), with the selection probability proportional to the posterior probability of the model. The small black vertical lines mark the hyperfine components. Spectra that have been scaled show the corresponding scaling factor on their right. Spectra are vertically offset for clarity.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Examples of N$_2$H$^+$ and N$_2$D$^+$ model spectra fit to azimuthally averaged spectra of IM Lup for a few radial bins. Spectra are centered on the systemic velocity. The colored lines show the data, with the shaded regions corresponding to the $1\sigma$ uncertainty. The black curves show 50 randomly selected models drawn from the MCMC with the selection probability proportional to the posterior probability of the model. The small black vertical lines mark the hyperfine components. The N$_2$D$^+$ spectra have been scaled and show the corresponding scaling factor on their right. Spectra are vertically offset for clarity.</figcaption></figure>
<table><caption>**Table 4**
Free Parameters for the Fitting of Azimuthally Averaged Spectra to Derive HCN and DCN Column Densities</caption><thead><tr><th>Parameter</th><th>Prior Low<sup>a</sup></th><th>Prior High<sup>b</sup></th><th>Unit</th><th>LTE/ Non-LTE<sup>c</sup></th></tr></thead><tbody><tr><td>T<sub>ex</sub></td><td>10</td><td>100</td><td>[K]</td><td>LTE</td></tr><tr><td>log<sub>10</sub>N<sub>HCN</sub></td><td>3</td><td>18 (LTE), 16 (non-LTE)</td><td>log<sub>10</sub>([cm<sup>-2</sup>])</td><td></td></tr><tr><td>log<sub>10</sub>N<sub>DCN</sub></td><td>3</td><td>18 (LTE), 16 (non-LTE)</td><td>log<sub>10</sub>([cm<sup>-2</sup>])</td><td></td></tr><tr><td>FWHM<sub>B3</sub></td><td>0.5<sup>d</sup></td><td>variable<sup>e</sup></td><td>[km s<sup>-1</sup>]</td><td></td></tr><tr><td>FWHM<sub>B6</sub></td><td>0.2<sup>d</sup></td><td>variable<sup>e</sup></td><td>[km s<sup>-1</sup>]</td><td></td></tr><tr><td>Δν<sub>HCN 1-0</sub></td><td>-0.3</td><td>0.3</td><td>[km s<sup>-1</sup>]</td><td></td></tr><tr><td>Δν<sub>H13CN 1-0</sub></td><td>-0.3</td><td>0.3</td><td>[km s<sup>-1</sup>]</td><td></td></tr><tr><td>Δν<sub>HCN 3-2</sub></td><td>-0.2</td><td>0.2</td><td>[km s<sup>-1</sup>]</td><td></td></tr><tr><td>Δν<sub>DCN 3-2</sub></td><td>-0.2</td><td>0.2</td><td>[km s<sup>-1</sup>]</td><td></td></tr><tr><td>T<sub>kin</sub></td><td>10</td><td>90</td><td>[K]</td><td>non-LTE</td></tr><tr><td>log<sub>10</sub>n<sub>H2</sub></td><td>3</td><td>10</td><td>log<sub>10</sub>([cm<sup>-3</sup>])</td><td>non-LTE</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Lower bound of flat prior.

<sup>b</sup> Upper bound of flat prior.

<sup>c</sup> Marks parameters used exclusively in the LTE or the non-LTE run.

<sup>d</sup> Equal to the channel width.

<sup>e</sup> Initial value for innermost radial bin is 20 km s<sup>-1</sup>. Dynamically adjusted as larger and larger radii are fitted (see Figure 29 and text).

analysis. Combining the 3–2 and 1–0 transitions allows us to constrain the gas temperature.

For each radial bin, we fit the azimuthally averaged spectra calculated in Section 3.4 of all four lines simultaneously, using the image cubes with a circular 0.′′3 beam. We start by assuming local thermodynamic equilibrium (LTE), that is, the excitation temperature $T_{ex}$ is the same for all transitions and equals the kinetic gas temperature $T_{kin}$. The free parameters are the excitation temperature, the HCN column density, the DCN column density, and additional parameters describing the line width and velocity offsets of the spectra. This results in a total of nine free parameters (see Table 4). We assume that HCN, DCN, and H<sup>13</sup>CN have the same temperature, that is, that they are spatial. Although DCN might be more concentrated toward the midplane compared with HCN, their theoretically expected spatial distributions (Aikawa et al. 2018) are similar enough to justify this first-order approximation. The general similarity of the radial emission profiles of HCN 3–2 and DCN 3–2 (Figure 3) also supports this assumption (Huang et al. 2017).

We fix the H<sup>13</sup>CN/HCN ratio to the ISM value of 13C/12C = 1/68 (Milam et al. 2005). To explore the parameter space, we use the MCMC method implemented in the emcee package (Foreman-Mackey et al. 2013). The details of the model and the fitting procedure are described in Appendix E. A few example fits are shown in Figure 5 for AS 209, and the full gallery of fits is shown in Appendix D.2.

As can be seen in Figure 5, in the innermost region (roughly within one beam FWHM from the disk center, i.e., two radial bins, corresponding to 30–47 au depending on the disk), the lines are strongly broadened by the velocity gradient within the beam and the hyperfine structure is not resolved (see also Figures 18, 19, 20, 21, and 22 in Appendix D.2). We caution that the strong broadening in the inner two radial bins might introduce additional uncertainties for the inferred column densities that are not reflected by our error bars.

Figures 7 and 8 show the derived optical depths, temperatures, and column densities for all five sources. The HCN 3–2 transition is optically thick ($\tau_0 > 1$) in the inner ~100 au for all sources except IM Lup.

We next performed a fit without assuming LTE, that is, the excitation temperatures are not necessarily equal to the kinetic gas temperature. In this case, we fit for the kinetic gas temperature $T_{kin}$ and the H<sub>2</sub> number density. Details are again given in Appendix E. The results are shown in Figure 9. The H<sub>2</sub> number density is constrained to $\gtrsim 10^6$ cm<sup>-3</sup> for all disks for $r \lesssim 200–400$ au. For comparison, the critical density is $2 \times 10^5$ cm<sup>-3</sup> and $4 \times 10^6$ cm<sup>-3</sup> for HCN 1–0 and 3–2, respectively (Dutrey et al. 1997). The fitted H<sub>2</sub> number density is mostly consistent with the midplane H<sub>2</sub> number density from the MAPS reference disk models (Zhang et al. 2021), shown with black dotted lines. The column densities are not significantly different from the LTE case, as can be seen in Figure 8, where the non-LTE results are overplotted for comparison. This suggests that the emission is indeed in LTE. The exception is the region around 500 au of the disks around IM Lup and HD 163296, where non-LTE conditions might prevail (see Section 5.1.1).

## 4.2. Radial Column Density Profiles for N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup>

We derive column densities of N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup> with the same procedure as in Section 4.1, assuming LTE. However, since we observed only one transition for each molecule, we did not fit the excitation temperature. Theoretical models predict that the N<sub>2</sub>H<sup>+</sup> abundance peaks around the threshold temperature of CO freeze-out (e.g., Aikawa et al. 2015, 2018), which is constrained by observations to ~20 K (Qi et al. 2011; Schwarz et al. 2016; Pinte et al. 2018). Thus, for simplicity we fixed the excitation temperature of both molecules to 20 K for our fiducial fits. This results in a total of six free parameters that are fitted for each radial bin: the column densities, the velocity offsets, and the parameter describing line broadening for each of the two lines. These parameters and their priors are listed in Table 5. The hyperfine structure is taken into account in the same way as for HCN and DCN. Tables 7 and 8 list the 28 and 25 hyperfine components considered for the fitting of N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup>, respectively. Figure 6 shows example fits for IM Lup for a few selected radial bins. The complete gallery of fits is shown in Figures 23–27. Similar to the HCN and DCN column densities, we caution that the column densities of the two innermost radial bins might be less reliable than those in the outer regions because of strong line broadening (see, e.g., Figures 23 and 24 for extreme examples of broadening of N<sub>2</sub>H<sup>+</sup> 3–2).

Figure 10 shows the derived optical depth and column density profiles for all five sources. The N<sub>2</sub>H<sup>+</sup> 3–2 transition is optically thick around the emission peak for IM Lup, GM Aur, and AS 209. The column densities of N<sub>2</sub>H<sup>+</sup> are close to $10^{13}$ cm<sup>-2</sup> for these sources, while they only reach $\lesssim 10^{12}$ cm<sup>-2</sup> for the two warmest disks in the sample around the Herbig stars HD 163296 and MWC 480. For all sources, the N<sub>2</sub>D<sup>+</sup> 3–2 transition is optically thin throughout the whole disk. The N<sub>2</sub>D<sup>+</sup> peak column densities exceed $10^{11}$ cm<sup>-2</sup> except for GM Aur.

In order to check the dependence of the column density on the assumed excitation temperature, we also run fits assuming that the excitation temperature equals the midplane gas temperature $T_{mid}$ shown in Figure 30 of Appendix E.1. These
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Optical depths from the MCMC calculation assuming LTE (i.e., $T_{ex} = T_{kin}$, see Appendix E.1) for the HCN 1–0, HCN 3–2, DCN 3–2, and H$^{13}$CN 1–0 lines. The black solid line shows the median. The blue shaded regions encompass the 16th to 84th, 2.3th to 97.7th, and 0.15th to 99.85th percentile regions. The horizontal dotted line marks an optical depth of 1. The beam major axis is shown as a horizontal black line in the top right of each panel.</figcaption></figure>

temperature profiles are extracted from the MAPS reference models (Zhang et al. 2021). As can be seen in Figure 10, there are some areas where the two assumptions about the excitation temperature result in different column densities. The strongest difference is seen for IM Lup, where the N$_2$H$^+$ column density increases by more than an order of magnitude if assuming $T_{ex} = T_{mid}$. This is because for IM Lup, the model midplane temperature is as low as ~7 K. However, even with such a high column density, the fit actually underpredicts the flux of the main N$_2$H$^+$ component by a factor of ~2, while at the same time overpredicting the flux of the hyperfine components (see Figure 23). This strongly suggests that the N$_2$H$^+$ excitation temperature is actually ≳10 K. Beyond ~300 au where the model midplane temperature rises to ~15 K, the column densities derived for the two excitation temperature choices agree again. Compared to N$_2$H$^+$, the N$_2$D$^+$ optical depth and column density of IM Lup show a less strong increase if we assume $T_{ex} = T_{mid}$.

Similarly, when assuming $T_{ex} = T_{mid}$, the N$_2$H$^+$ column density toward AS 209 is higher by a factor ~4 around ~75 au (where $T_{mid} \approx 11$ K; see Figure 30), but again the N$_2$H$^+$ data are not fit well by this model (see Figure 25), in contrast to the fit where $T_{ex} = 20$ K. Finally, Figure 10 shows that, when assuming $T_{ex} = T_{mid}$, the N$_2$H$^+$ column density toward GM Aur is increased by up to a factor of ~3 in the region between ~50 au and ~250 au. In this region, $T_{mid} \approx 12$ K. However, in contrast to IM Lup and AS 209, the model assuming $T_{ex} = T_{mid}$ fits the N$_2$H$^+$ data as good as the model

assuming $T_{ex} = 20$ K (see Figure 24). Thus, if one believes that the temperature in the emitting region is indeed as low as 12 K, the N$_2$H$^+$ column density toward GM Aur would be a factor of ~3 larger. From a theory point of view, this low temperature is not unreasonable; the models by Aikawa et al. (2015; see their Figure 10) show that, while the N$_2$H$^+$ abundance should peak around the CO freeze-out temperature (~20 K), it can remain substantial even below the N$_2$ freeze-out temperature of ~17 K (see also Aikawa et al. 2021).

## 4.3. Radial Profiles of the Deuteration Fraction

### 4.3.1. DCN/HCN

Using the column density profiles of HCN and DCN calculated in Section 4.1, we now derive the radial profile of the deuteration fraction of HCN. Figure 11 (upper panel) shows the radial profiles of the DCN/HCN ratio. The results from the LTE and non-LTE fits agree well. The inferred DCN/HCN ratios range from ~10$^{-3}$ to ~10$^{-1}$. For IM Lup, the outer DCN ring at ~350 au is more strongly deuterated than the inner ring at ~100 au. The data also suggest that the deuteration further decreases inward of the inner ring, although not at high significance. For the disk of GM Aur, the outer ring at ~300 au is more strongly deuterated than the inner ~100 au of the disk. AS 209 shows a weak decreasing trend of the deuteration toward the inner parts of the disk. Finally, for the two disks around the Herbig stars HD 163296 and MWC 480, the
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Temperatures and HCN and DCN column densities derived from MCMC fits of azimuthally averaged spectra assuming LTE (i.e., $T_{ex} = T_{kin}$, see Appendix E.1). The solid black line shows the median. The blue shaded regions encompass the 16th to 84th, 2.3th to 97.7th, and 0.15th to 99.85th percentile regions. For comparison, the median values (thick dotted lines) and 16th and 84th percentiles (thin dotted lines) of the kinetic temperature and column densities derived from non-LTE fits are shown. The beam major axis is shown as a horizontal black line in the upper right of each panel.</figcaption></figure>

deuteration in the innermost region is reduced by almost two orders of magnitude compared to the deuteration at ~150 au. Thus, while there is some diversity in the deuteration profiles toward our five targets, generally the outer disk regions are more strongly deuterated in HCN compared with the inner disk. In the inner 100 au of the disk around GM Aur, the deuteration seems to be increasing toward the star instead, although higher S/N and angular resolution data would be needed for confirmation. These trends are further discussed in Section 5.1.4.

### 4.3.2. $N_2D^+/N_2H^+$

The lower panels of Figure 11 show radial profiles of the $N_2D^+/N_2H^+$ column density ratio. We find a $N_2D^+/N_2H^+$ ratio typically between a few times $10^{-2}$ and 1. The profiles do not show a strong dependence on the assumed excitation temperature. The radial distribution of the $N_2D^+/N_2H^+$ ratio varies among sources. For IM Lup, the ratio has a minimum located at ~220 au; this structure reflects the gap in the column density profile of $N_2D^+$. For GM Aur, $N_2D^+$ remains undetected while strong emission of $N_2H^+$ is seen, resulting in upper limits on the ratio. For AS 209, the ratio monotonically increases with radius from ~50 to 170 au. For the disk around HD 163296, the profile is rather flat. For MWC 480, the $N_2D^+/N_2H^+$ ratios is a few times $10^{-1}$, but the low angular resolution (beam major axis of 0.94) of the $N_2H^+$ data precludes an accurate determination of the $N_2D^+/N_2H^+$ profile.

# 5. Discussion

## 5.1. HCN and DCN

### 5.1.1. Comparing LTE and Non-LTE

It is interesting to compare the HCN gas temperatures derived from the LTE and non-LTE fits (Figure 8 top row). In general, the temperature is not well constrained, although there are some disk regions where information on the temperature can be extracted. For example, $T \gtrsim 30$ K around 150 au in the disk around IM Lup, and $T \approx 25$ K in the inner ~100 au of the disk around AS 209. Despite the large uncertainties, we find significant differences between LTE and non-LTE temperature estimates in some disk regions: from 300 to 550 au for IM Lup, at ~110 au for AS 290, inward of 100 au and between 300 and 500 au for HD 163296, and inward of 70 au for MWC 480. We find that in these regions, non-LTE generally provides a better fit to the data. In particular, the LTE fits tend to underpredict the HCN 1–0 emission. As an example, this can clearly be seen in Figure 18 around 500 au for IM Lup and in Figure 21 around 450 au for HD 163296. For these two outer disk regions, our fits thus suggest that the gas density is low enough for non-LTE conditions to prevail. This is further supported by Figure 9 that shows that the $H_2$ density might be below the critical density of HCN 3–2. The HCN column density derived from the non-LTE fit is higher in those two regions compared with the LTE case.

However, non-LTE seems more unlikely an explanation for the disk regions $\lesssim 100$ au mentioned above, where the gas density is expected to be higher. Instead, for those regions, the difference between the LTE and non-LTE temperatures might indicate that our assumption of a single temperature being able
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Kinetic gas temperature, HCN and DCN column densities, and H₂ number density derived from non-LTE MCMC fits of azimuthally averaged spectra. The solid black line shows the median. The blue shaded regions encompass the 16th to 84th, 2.3th to 97.7th, and 0.15th to 99.85th percentile regions. The beam major axis is shown as a horizontal black line in the upper right of each panel. The black dashed curves in the last row show the H₂ number density in the midplane of the MAPS reference disk models by Zhang et al. (2021), while the horizontal dotted and dashed lines show the critical density of HCN 1–0 and 3–2, respectively.</figcaption></figure>

<table><caption>**Table 5**
Free Parameters for the Fitting of Azimuthally Averaged Spectra to Derive N₂H⁺ and N₂D⁺ Column Densities</caption><thead><tr><th>Parameter</th><th>Prior Low<sup>a</sup></th><th>Prior High<sup>b</sup></th><th>Unit</th></tr></thead><tbody><tr><td>log<sub>10</sub>N<sub>N2H+</sub></td><td>3</td><td>18</td><td>log<sub>10</sub>([cm<sup>-2</sup>])</td></tr><tr><td>log<sub>10</sub>N<sub>N2D+</sub></td><td>3</td><td>18</td><td>log<sub>10</sub>([cm<sup>-2</sup>])</td></tr><tr><td>FWHM<sub>N2H+</sub> 3–2</td><td>variable<sup>c</sup></td><td>variable<sup>d</sup></td><td>[km s<sup>-1</sup>]</td></tr><tr><td>FWHM<sub>N2D+</sub> 3–2</td><td>0.2<sup>c</sup></td><td>variable<sup>d</sup></td><td>[km s<sup>-1</sup>]</td></tr><tr><td>Δν<sub>N2H+</sub> 3–2</td><td>-0.2</td><td>0.2</td><td>[km s<sup>-1</sup>]</td></tr><tr><td>Δν<sub>N2D+</sub> 3–2</td><td>-0.2</td><td>0.2</td><td>[km s<sup>-1</sup>]</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Lower bound of flat prior.

<sup>b</sup> Upper bound of flat prior.

<sup>c</sup> Equal to the channel width (see Table 2).

<sup>d</sup> Initial value for innermost radial bin is 20 km s<sup>-1</sup>. Dynamically adjusted as larger and larger radii are fitted (see Figure 29 and text).

to describe the emission of all lines is invalid. For example, since the HCN emission is optically thick in those regions (Figure 7), it could be that HCN 3–2 is emitting from a different vertical layer with a different temperature than HCN 1–0 and DCN 3–2. Similarly, it could be that the hyperfine components trace a different temperature compared with the main component due to optical depth. Another possibility is

effects due to dust optical depth that might affect the 3–2 and 1–0 transitions differently. Source-specific radiative transfer modeling would be necessary to investigate these ideas. Fortunately, the column densities do not differ strongly between the LTE and non-LTE fits in those inner ~100 au regions.

### 5.1.2. Comparison with Previous Observations

Bergner et al. (2019, 2020a) derived HCN column density profiles toward AS 209, HD 163296, and MWC 480 by using observations of the 3–2 transition of HCN and H<sup>13</sup>CN. Their data were of significantly lower resolution: ~0.′′5 for AS 209 and HD 163296 and ~2.′′7 for MWC 480. Our HCN column density profiles show a much greater dynamical range compared with Bergner et al. (2019, 2020a), but they show order of magnitude agreement in terms of mean column densities. Bergner et al. (2019, 2020a) derive excitation temperatures mostly below 20 K, that is, lower than the temperatures we derive. To test whether these differences arise because of the higher spatial resolution of our data (i.e., less beam dilution), we convolved our data to match the resolution of the Bergner et al. (2019) data and repeated the fitting. We find that even in that case, significant differences in column density and excitation temperature are present. This suggests that the lower resolution of the Bergner et al. (2019, 2020a) data is not the primary reason for the observed differences but rather that their analysis did not
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Optical depths and column density profiles of $N_2H^+$ and $N_2D^+$. The results from the fiducial fit ($T_{ex}$ = 20 K) are shown with the solid black line (median) and blue shaded regions encompassing the 16th to 84th, 2.3th to 97.7th, and 0.15th to 99.85th percentile regions. For comparison, the median (thick dotted lines) and 16th and 84th percentile (thin dotted lines) from the alternative fit ($T_{ex}$ = $T_{mid}$) are shown. In the first two rows, the horizontal dotted line marks an optical depth of 1. The beam size is shown as a horizontal black line in the upper right of each panel.</figcaption></figure>

include the HCN 1–0 transition and that they were unable to use the hyperfine structure of HCN 3–2.

Guzmán et al. (2021) and Bergner et al. (2021) also fitted spectra extracted from the MAPS data to derive HCN column density profiles. Bergner et al. (2021) only used the HCN 3–2 data, while Guzmán et al. (2021) used HCN 3–2 and 1–0. Neither of them used $H^{13}CN$ 1–0. In general, the profiles show good agreement with ours. Some differences are observed in the inner 50–100 au, where Guzmán et al. (2021) and Bergner et al. (2021) derive up to an order of magnitude higher HCN column densities for GM Aur, AS 209, HD 163296, and MWC 480. Adopting those values would further strengthen the trend of low DCN/HCN toward the disk center. We refer to Appendix C of Guzmán et al. (2021) for a comparison of the derived HCN column densities.

The disk-integrated DCN/HCN ratios derived by Huang et al. (2017) and Salinas et al. (2017) (orange shaded areas in Figure 11) are in good agreement with our results. Concerning the radial variation of DCN/HCN, our results are in line with trends of enhanced deuteration in the outer disks seen in previous work. Huang et al. (2017) discuss the possibility of enhanced deuteration in the outer disks of V4046 Sgr, LkCa 15, and HD 163296 by comparing the emission morphologies of $H^{13}CN$ and DCN. More recently, Facchini et al. (2021) presented a radial profile of the DCN/HCN column density ratio

for the PDS 70 disk at 0.42 (47 au) resolution. They also identify an increasing HCN deuteration with increasing radius. Similarly, Öberg et al. (2021) found that the DCN emission in the TW Hya disk drops toward the center, while the HCN emission is centrally peaked (Hily-Blant et al. 2019). While a determination of the column density DCN/HCN ratio is pending, TW Hya may well show the same trend of an increasing DCN/HCN ratio with radius.

### 5.1.3. Comparison of HCN and DCN Column Densities to Model Predictions

In this section, we compare the derived HCN and DCN column densities to some models presented in the literature (Willacy 2007; Walsh et al. 2012; Aikawa et al. 2018; Cleeves et al. 2018). Among these models, the Cleeves et al. (2018) model is the only one tailored to a specific MAPS source: IM Lup. The others are generic models of disks around T Tauri stars. This should especially be kept in mind when comparing these models to the disks around the Herbig Ae stars HD 163296 and MWC 480.

We first consider the HCN column density. The models by Willacy (2007), Walsh et al. (2012), and Aikawa et al. (2018) predict a flat HCN profile (except for the innermost few astronomical units). In contrast, IM Lup shows an HCN profile slowly rising inward (except for the innermost 100 au), while
the disks around AS 209, GM Aur, HD 163296, and MWC 480 all show HCN column densities that are steeply rising toward the star inward of 100–200 au (Figure 8). Cleeves et al. (2018) modeled the HCN column density for a range of elemental C/O ratios and cosmic-ray ionization rates. Their models show increasing HCN column densities toward the star. Therefore, the Cleeves et al. (2018) models actually better reproduce the morphology of the HCN column density profiles for all disks (depending on the assumed C/O), although they were tailored to IM Lup. The Cleeves et al. (2018) models also reproduce reasonably well the order of magnitude values of the HCN column densities we derive. One possible reason might be that Cleeves et al. (2018) assumed an elevated interstellar radiation field strength ($G_0 = 4$). Also, for the UV radiation from the central star, they scaled the TW Hya UV spectrum to be consistent with the flux density of IM Lup. These points might matter, since UV photons can dissociate HCN, and the UV radiation field is thus an important parameter influencing the HCN column density profile (e.g., Bergin et al. 2003; Chapillon et al. 2012; Bergner et al. 2021). The Cleeves et al. (2018) study also demonstrated that the HCN column density depends strongly on the C/O ratio. Thus, future work might be able to constrain the C/O ratio using observed HCN column density profiles. This could then be compared to the C/O ratio inferred from, for example, the study of C$_2$H (e.g., Bergin et al. 2016; Alarcón et al. 2021; Bosman et al. 2021a) or sulfur-bearing molecules (Le Gal et al. 2021).

In the innermost ~5 au, the models by Willacy & Woods (2009), Walsh et al. (2012), and Aikawa et al. (2018) all predict an up to six orders of magnitude increase of the HCN column density due to thermal desorption. The angular resolution of our data is insufficient to test this prediction directly. A detailed analysis of the line kinematics may be able to probe this innermost region (e.g., Bosman et al. 2021b), but this is beyond the scope of the present paper.

Finally, we compare our DCN column density profiles to the models by Willacy (2007) and Aikawa et al. (2018). They typically predict DCN column densities of ~$10^{12}$ cm$^{-2}$, that is, the order of magnitude agrees with the column densities we derive (see Figure 8).

In summary, the HCN and DCN column density profiles derived here provide useful benchmarks for disk chemistry models. In particular, we expect that parameters such as the UV field or the C/O ratio can be constrained in future modeling work.

### 5.1.4. The Radial Variation of the DCN/HCN Ratio

A major result of our work is the variation of the DCN/HCN ratio within the disks: generally, DCN/HCN decreases when moving from the outer disk regions ($r \gtrsim 100$ au) toward the disk center (Figure 11). For IM Lup, this decrease is seen when comparing the outer ring to the inner ring, as well as when comparing the inner ring to the disk center. For GM Aur, the decrease is seen by comparing the outer ring at ~300 au to the disk's central ~100 au, although DCN/HCN might be increasing from 100 au toward the disk center. For AS 209, the decrease is quite weak from ~150 au toward the disk center. For all these T Tauri disks, the decrease of DCN/HCN is not more than an order of magnitude. On the other hand, for the disks around the Herbig stars HD 163296 and MWC 480,

the ratio decreases by almost two orders of magnitude, from $10^{-1}$ to almost $10^{-3}$, when moving inward from ~150 au.

Generally speaking, we interpret the decrease of the DCN/HCN ratio toward the disk center as a confirmation that HCN fractionation is due to in situ exothermic exchange reactions such as those shown in Equations (1)–(3). Indeed, the temperature is expected to increase toward the disk center, which in turn increases the efficiency of the reverse reactions, thus lowering the deuteration. To get further insight, we also compare our DCN/HCN profiles to models from the literature. Interestingly, the disk chemistry model (assuming a T Tauri host star) by Favre et al. (2015) predicts an increasing DCN/HCN ratio from 60 au inward. Unfortunately, the DCN/HCN inward of 60 au is not constrained well for our disks, and the model does not go to larger radii. A comparison is therefore difficult, but we can at least say that the qualitative trend of the DCN/HCN model profile is different from the data, with the notable exception of the inner ~100 au of the disk around GM Aur. On the other hand, the disk model (assuming a T Tauri host star) by Aikawa et al. (2018) predicts a decreasing DCN/HCN ratio when moving inward from 300 to 10 au (see their Figure 5(d)). One of the updates of the Aikawa et al. (2018) model compared with Favre et al. (2015) is the inclusion of the ortho/para ratio of H$_2$. An elevated ortho/para ratio of H$_2$ in the inner disk compared with the outer regions can reduce the efficiency of the deuteration.

For a more detailed interpretation of the DCN/HCN radial profiles, source-specific modeling will be needed. In any case, if the gas we observe is the main source of HCN for forming icy bodies in the disk, we expect the formation location to be imprinted in the DCN/HCN ratio of the body, that is, the higher the DCN/HCN ratio, the further out in the disk the body formed. However, there remains a caveat that the gas we observe may not be the primary HCN source of forming comets (see Section 5.1.7).

### 5.1.5. DCN/HCN as a Function of Temperature

In order to study the contributions of the low- and high-temperature pathways to HCN deuteration, we now look at the deuteration fraction versus the temperature estimated from the multiline analysis. In Figure 12, we show the temperature as a function of radius (top) and the DCN/HCN ratio as a function of temperature (bottom). We only include points with a relative error on the temperature smaller than 20%. This is a somewhat arbitrary choice, but we found that a stricter cutoff results in only a few points available for analysis, while a more generous cutoff results in a large number of data points with poorly constrained temperatures.

If, for a given radius, the LTE and non-LTE fits give similar results (relative difference smaller than 15% for both the fitted temperature and DCN/HCN ratio), we only consider the LTE fit. Data points within half a beam FWHM from the disk center are excluded, because the additional line broadening makes the inference of column densities and temperatures more uncertain. We further define that the upper or lower bound of the fitted DCN/HCN is unconstrained whenever the 99th percentile is larger than 10 or the 1st percentile is smaller than $10^{-4}$, respectively. We exclude data points where both bounds are unconstrained by this definition and plot upper and lower limits whenever only the lower or upper bound is unconstrained, respectively.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Radial profiles of the HCN/DCN (top) and N<sub>2</sub>D<sup>+</sup>/N<sub>2</sub>H<sup>+</sup> column density ratios. The solid black lines show the median of the LTE fits for DCN/HCN and of the fits with T<sub>ex</sub> = 20 K for N<sub>2</sub>D<sup>+</sup>/N<sub>2</sub>H<sup>+</sup>. The blue shaded regions encompass the 16th to 84th, 2.3th to 97.7th, and 0.15th to 99.85th percentile regions. Medians (thick dotted line) and the 16th and 84th percentiles (thin dotted lines) of the alternative fits (DCN/HCN: non-LTE; N<sub>2</sub>D<sup>+</sup>/N<sub>2</sub>H<sup>+</sup>: T<sub>ex</sub> = T<sub>mid</sub>) are shown for comparison. The horizontal orange shaded areas show the disk-integrated HCN/DCN ratios inferred by Huang et al. (2017, IM Lup, AS 209, HD 163296, MWC 480) and Salinas et al. (2017, HD 163296) and disk-integrated N<sub>2</sub>D<sup>+</sup>/N<sub>2</sub>H<sup>+</sup> ratios by Huang & Öberg (2015, AS 209) and Salinas et al. (2017, HD 163296). The beam major axis is shown with the horizontal black line in the upper left.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Top: the gas temperature as a function of radius. The solid lines show power-law fits to the <sup>12</sup>CO 2–1 brightness temperature derived by Law et al. (2021b). Bottom: the DCN/HCN column density ratio as a function of gas temperature. Upper and lower triangles mark upper and lower limits corresponding to the 99th and 1st percentiles, respectively. The black and red lines are upper limits on the HCN deuteration from the low- and high-temperature pathways, respectively. They show the H<sub>2</sub>D<sup>+</sup>/H<sub>3</sub><sup>+</sup> and CH<sub>2</sub>D<sup>+</sup>/CH<sub>3</sub><sup>+</sup> ratios calculated by balancing Equations (1) and (2), respectively. For CH<sub>2</sub>D<sup>+</sup>/CH<sub>3</sub><sup>+</sup>, we show ratios calculated using the exothermicity of Nyman & Yu (2019, solid) and Roueff et al. (2013, dotted). In both panels, only data points with a relative error on the temperature smaller than 20% are included. See the main text for more details on the data point selection. Data points from disk regions where the LTE and non-LTE temperatures disagree are plotted with open symbols to indicate the possibility of systematic uncertainty. Data points exceeding the CO temperature are plotted fainter.</figcaption></figure>

In Section 5.1.1, we noticed that the gas temperatures derived from the LTE and non-LTE fits differ in certain regions of the disks. For the regions beyond 300 au in IM Lup and HD 163296, these differences can be explained by non-LTE conditions, that is, the excitation temperature does not equal the gas temperature. As a consequence, we only consider the gas temperature derived from the non-LTE fits for r > 300 au for IM Lup and HD 163296, but it turns out that no non-LTE data

points from those regions pass the temperature error bar cutoff. For the inner regions of AS 209, HD 163296, and MWC 480 where the LTE and non-LTE temperatures disagree as well, non-LTE seems unlikely, making the inferred temperatures uncertain. Thus, we plot data points from those regions (90 to 140 au for AS 209, 0 to 130 au for HD 163296, and 0 to 70 au for MWC 480) with open symbols to indicate the possibility of systematic uncertainty. We also plot in Figure 12 power-law
<figure><img src="image_16.png" /><figcaption>**Figure 13.** The column density profiles of N<sub>2</sub>H<sup>+</sup> (blue), N<sub>2</sub>D<sup>+</sup> (orange), and CO (gray, scaled by 10<sup>-8</sup> for HD 163296 and MWC 480 and 10<sup>-6</sup> for the other targets). The CO column densities are from Zhang et al. (2021). The color-shaded regions extend from the 16th to 84th percentile of the MCMC fits for N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup> and over the 1σ uncertainty for CO. The vertical green bars indicate the CO snow line locations inferred from the MAPS reference models (Zhang et al. 2021). The vertical gray-shaded regions indicate the CO snow line locations estimated from N<sub>2</sub>H<sup>+</sup> observations by Qi et al. (2015, 2019).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Comparison of the radial emission profiles of DCN 3–2 (blue), N<sub>2</sub>D<sup>+</sup> 3–2 (orange), and DCO<sup>+</sup> (green). For the latter, we show the 4–3 transition in GM Aur and the 3–2 transition in the other disks. The vertical gray lines indicate the outer edge of the millimeter dust continuum emission determined by Law et al. (2021a). The size of the beam major axis for each emission line is indicated by a colored horizontal line in the upper right corner of each panel.</figcaption></figure>

fits to the <sup>12</sup>CO 2–1 brightness temperature (Law et al. 2021b) for guidance, since the HCN temperature is unlikely to exceed the CO temperature.

The bottom panel of Figure 12 also shows the H<sub>2</sub>D<sup>+</sup>/H<sub>3</sub><sup>+</sup> and CH<sub>2</sub>D<sup>+</sup>/CH<sub>3</sub><sup>+</sup> ratios calculated by balancing Equations (1) (low-temperature pathway) and (2) (high-temperature pathway), respectively. We assumed a thermalized ortho-to-para ratio of H<sub>2</sub>, which is theoretically expected for protoplanetary disks (Aikawa et al. 2018; Furuya et al. 2019). For the CH<sub>2</sub>D<sup>+</sup>/CH<sub>3</sub><sup>+</sup> ratio, we show two curves, corresponding to the exothermicities by Roueff et al. (2013) and Nyman & Yu (2019). The H<sub>2</sub>D<sup>+</sup>/H<sub>3</sub><sup>+</sup> and CH<sub>2</sub>D<sup>+</sup>/CH<sub>3</sub><sup>+</sup> ratios are upper limits on the HCN deuteration by the low- and high- temperature pathways, respectively, since, for example, the destruction of H<sub>2</sub>D<sup>+</sup> by CO or the destruction of CH<sub>2</sub>D<sup>+</sup> by electrons is neglected.<fnref n="25" /> We see that for a considerable number of points, the low-temperature pathway alone cannot explain the observed deuteration fraction. Thus, a contribution from the high-temperature pathway is needed.

Naively, we might expect that the deuteration should decrease with increasing temperature. There are suggestions of such a trend with temperature in the bottom panel of Figure 12 for the data points from AS 209 and MWC 480. However, the presence of two deuteration pathways as well as possible systematic uncertainties in the derived temperatures complicates the interpretation. More data points with a well-constrained temperature and DCN/HCN

ratio would be necessary to study the temperature dependence of the deuteration in more detail.

### 5.1.6. Is the Drop of DCN/HCN in the Inner Disks an Artefact Caused by High Dust Optical Depth?

Here we consider the possibility that the observed decrease of the DCN/HCN ratio toward the disk centers is merely an artifact caused by high dust optical depth. Optically thick dust can block line photons from escaping (e.g., Cleeves et al. 2016). Furthermore, continuum subtraction can lead to underestimated line flux (Weaver et al. 2018). Naively, we might expect that taking the ratio of DCN/HCN at least partially cancels out such effects, but since DCN and HCN are not necessarily perfectly cospacial, there is a possibility that DCN is more strongly affected.

Both of the effects described above would occur in regions of high dust optical depth. Thus, we consider the dust optical depth at ALMA Band 6 wavelengths (~1.3 mm), which should be higher than in Band 3 (~3 mm). Previous studies suggested that the dust is mostly optically thin in ALMA Band 6, with the exception of the inner 10 au and 50 au for IM Lup and MWC 480, respectively (Huang et al. 2018, 2020; Liu et al. 2019). However, these studies did not consider the effects of scattering. Zhu et al. (2019) found that dust scattering can considerably reduce the emission from an optically thick region, which leads to underestimation of the dust optical depth. Therefore, we consider radial profiles of the vertical<fnref n="26" />

<footnote n="25">²⁵ The effect of the destruction of H<sub>2</sub>D<sup>+</sup> by CO is small (a factor of 2 at most), because at temperatures above 20 K and with a thermalized H<sub>2</sub> ortho/para ratio, the main destruction partner of H<sub>2</sub>D<sup>+</sup> is H<sub>2</sub> as long as CO/H<sub>2</sub> is lower than 10<sup>-4</sup>. Quantifying the CH<sub>2</sub>D<sup>+</sup> destruction by electrons is not straightforward, because we do not know the electron abundances in disks.</footnote>

<footnote n="26">²⁶ All MAPS disk have fairly low inclination, ranging from 35° for AS 209 to 53° for GM Aur (Öberg et al. 2021).</footnote>
dust optical depth derived from models including scattering by Sierra et al. (2021). These models show that the dust emission is optically thin for $r \gtrsim 60$ au. The exception is the disk around IM Lup that could be optically thick even for $r > 150$ au. Generally, we find that the DCN/HCN ratio starts to decrease at radii well beyond 60 au, even when considering the finite beam size. In addition, there is no clear correlation between the DCN/HCN profiles and the dust optical depth profiles by Sierra et al. (2021). This suggests that dust optical depth is not the dominant cause for the decrease of DCN/HCN toward the disk centers.

### 5.1.7. Comparison with the DCN/HCN Ratio in Comets

In this section, we attempt to connect our results to the formation history of cometary bodies. Unfortunately, the number of comets with a measured DCN/HCN ratio is limited (e.g., Bockelée-Morvan et al. 2015). For comet Hale-Bopp, DCN/HCN = (2.3 ± 0.4) × 10⁻³ (Meier et al. 1998; Crovisier et al. 2004). Upper limits of DCN/HCN < 10⁻² are reported for C/1996 B2 (Hyakutake, Bockelée-Morvan et al. 1998) and 103P/Hartley 2 (Gicquel et al. 2014). The DCN/HCN ratios measured in our sample of protoplanetary disks are typically ≳10⁻², with the exception of the inner ~50 au of the HD 163296 and MWC 480 disks (see Figure 11). The higher DCN/HCN ratios we derive compared with Hale-Bopp suggest that comets do not directly form from the bulk of the gas reservoir we observe. The simplest explanation might be that comets form from a gas reservoir located closer to the star than we probe here (Ceccarelli et al. 2014). Indeed, in the framework of the Nice model, both the Oort cloud (the reservoir of long-period comets) and the scattered disk (the reservoir of short-period, that is, Jupiter-family comets) originate from dynamical scattering of objects originally located in the Uranus–Neptune zone (at ~30 au) by the migrating giant planets (Brasser & Morbidelli 2013). Although the resolution of our data is insufficient to put strong constraints on the DCN/HCN in these inner regions of the disks, the observed trends of decreasing DCN/HCN toward the disk centers of IM Lup, AS 209, HD 163296, and MWC 480 are consistent with this picture. On the other hand, in the inner ~100 au of the disk around GM Aur, DCN/HCN is increasing toward the star, making such a scenario less likely for this target.

An alternative interpretation would be that comets do actually form in the outer disk regions but not from the gas we are observing here. Instead, they might form from ISM ices that were incorporated into the disk. The origin of cometary material (direct inheritance of ISM ice versus condensation in the protoplanetary disk or both) is still actively debated in the literature (e.g., Bockelée-Morvan et al. 2015; Willacy et al. 2015; Rubin et al. 2020, and references therein). HCN and DCN are chemically stable and have relatively high desorption energies (Noble et al. 2013). Thus, a significant amount of HCN and DCN ice formed in molecular clouds could survive in the disk and be incorporated into comets.

We expect that radially resolved deuteration profiles as presented in this work will be useful to constrain models of comet formation in the solar nebula such as those presented by Mousis et al. (2000).

### 5.1.8. Comparison with the DCN/HCN Ratio in Earlier Evolutionary Stages

The DCN/HCN column density ratios inferred for class I young stellar objects by Le Gal et al. (2020) and Bergner et al. (2020b) are within the range of values found for our sources. The DCN/HCN ratios inferred for even earlier evolutionary stages of star formation are also comparable to the results from our work: infrared dark clouds and hot molecular cores (Gerner et al. 2015) as well as high-mass star-forming clumps (Feng et al. 2019) show DCN/HCN ratios roughly of the order of 10⁻². Considering the relatively high desorption energy of HCN (3370 K, Noble et al. 2013), these single-dish observations should mostly probe HCN and DCN formed in the gas phase rather than sublimated from ice (the sublimation zone is small compared with the beam). In summary, this suggests that the gas-phase deuteration chemistry is similar between a wide range of ISM environments and protoplanetary disks.

## 5.2. N₂H⁺ and N₂D⁺

### 5.2.1. Comparison with Previous Observation of N₂D⁺ and Theoretical Predictions

N₂D⁺ has previously been detected toward only two protoplanetary disks (Huang & Öberg 2015; Salinas et al. 2017). Huang & Öberg (2015) observed N₂D⁺ 3–2 toward the disk around AS 209 with ALMA and compared it with a Submillimeter Array observation of N₂H⁺ by Öberg et al. (2011). Because the resolution (~1") and the sensitivity were not high enough to resolve the disk structure, they derived disk-averaged column densities for two fixed excitation temperatures: 10 and 25 K. They find column densities of N₂D⁺ and N₂H⁺ of (1.4–2.1) × 10¹¹ cm⁻² and (3.1–6.3) × 10¹¹ cm⁻², respectively. The N₂D⁺/N₂H⁺ column density ratio is then estimated as 0.3–0.5. While we derive higher column densities, our N₂H⁺ deuteration ratio agrees with their disk-averaged value (Figure 11).

Salinas et al. (2017) observed the N₂D⁺ 3–2 transition toward the HD 163296 disk with ALMA and derived a disk-averaged N₂D⁺ column density of (2.5 ± 0.3) × 10¹¹ cm⁻² for the assumed excitation temperature of 10 K (for 25 K, they derive (1.6 ± 0.2) × 10¹¹ cm⁻²). They also estimated the N₂H⁺ column density from the N₂H⁺ 3–2 observation by Qi et al. (2015) and derived an N₂D⁺/N₂H⁺ column density ratio of 0.34 ± 0.15 and 0.45 ± 0.21 for the assumed excitation temperatures of 10 and 25 K, respectively. We derive a slightly higher N₂D⁺ column density and N₂D⁺/N₂H⁺ ratio (Figures 10 and 11). One reason could be that their values are averaged over the whole disk, including the inner disk, which shows negative emission, while our higher-resolution observations allow us to ignore that inner region. Indeed, their integrated N₂D⁺ 3–2 flux (from which they estimate the column density) is significantly lower compared with our value (61.6 ± 7.5 mJy km s⁻¹, compared with our 135 ± 15 mJy km s⁻¹, where the latter is derived by excluding the inner 50 au, see Table 3).

The radial emission profile of N₂D⁺ 3–2 in Salinas et al. (2017) also shows a morphology similar to our data (see Figure 4): a ring extending from ~100 to 300 au. The line flux is negative in the inner-disk region in both Salinas et al. (2017) and our profile. This negative emission is probably due to the difficulty of achieving an accurate continuum subtraction, possibly caused by an atmospheric absorption feature. We confirmed, however, that our N₂D⁺ profile
is not affected by this issue outside of 50 au. See Appendix B for a
more detailed discussion.

Based on the chemical disk modeling by Aikawa et al. (2018), $N_2D^+$ is expected to be present mainly in the cold outer regions (≥50 au) of the disk. We find that for all of our sources where $N_2D^+$ is detected, it is indeed distributed in the outer regions, peaking beyond ~50 au. This is consistent with the theoretical expectation that the deuteration of $N_2H^+$ proceeds via the low-temperature channel. In addition to the radial profile, the $N_2D^+/N_2H^+$ ratio derived from the observations is in reasonable order of magnitude agreement with theoretical models. In Aikawa et al. (2018), the ratio reaches almost unity in the outer disk, while it declines to ≤0.1 inside a radius of ~100 au.

In summary, we find that both previous observations and predictions by chemical disk models are broadly consistent with our analysis of the $N_2H^+$ and $N_2D^+$ data. The overall picture of efficient $N_2H^+$ deuteration in the cold midplane of the outer disk is confirmed.

### 5.2.2. $N_2D^+$ as CO Snow Line Tracer

In this section, we discuss the prospects of using $N_2D^+$ as a CO snow line tracer. Determining the CO snow line from observations of CO emission lines is generally challenging, because even outside of the CO snow line, there is a warm molecular layer with CO gas. This can be seen in Figure 13, which shows the CO column density profiles for the MAPS disks determined from C$^{18}$O observations (Zhang et al. 2021). The vertical green bars indicate the model estimates for the CO snow line, that is, the radii where the abundances of CO gas and CO ice become the same in the MAPS reference models (Zhang et al. 2021). In the two disks around the Herbig stars, HD 163296 and MWC 480, the CO column density profiles rapidly increase inward across the CO snow lines, tracing CO evaporation (Zhang et al. 2021). However, no clear trend is seen for the other disks, that is, the snow line cannot be determined from the CO column density profile.

Instead, the snow line is usually determined from indirect tracers. In particular, $N_2H^+$ has been used to trace the CO snow line (Qi et al. 2013, 2015, 2019). Since $N_2H^+$ is destroyed by CO ($N_2H^+ + CO \rightarrow HCO^+ + N_2$), $N_2H^+$ is expected to be abundant outside of the CO snow line. Therefore, the inner edge of the $N_2H^+$ emission can be used as an estimate of the CO snow line. However, models by van't Hoff et al. (2017) and Aikawa et al. (2018) predict that $N_2H^+$ can, in addition to the midplane region, also form in a surface layer of the disk. This surface layer can make the inference of the CO snow line location considerably more uncertain (van't Hoff et al. 2017). Thus, $N_2D^+$ might be a better snow line tracer, since no surface layer is expected to form (Aikawa et al. 2018). Here we test how well $N_2D^+$ traces the CO snow line.

In Figure 13, we plot the column density profiles of $N_2H^+$ and $N_2D^+$ (remember that the column densities might carry additional systematic uncertainties for $r \lesssim 50$ au due to severe broadening; see Section 4.2) and the CO snow line estimate derived from $N_2H^+$ observations (gray-shaded area, Qi et al. 2015, 2019) and from the MAPS reference models (green bar, Zhang et al. 2021). Except for GM Aur, where the S/N of $N_2D^+$ is low, the estimates for the CO snow line approximately correspond to the inner edges of the $N_2D^+$ column density profiles. Unfortunately, an accurate measurement of the inner edge of $N_2D^+$ can be challenging since it is difficult to achieve

an accurate continuum subtraction for the $N_2D^+$ 3–2 line, which can result in unphysical negative emission in the inner-disk region (see Appendix B). Furthermore, $N_2D^+$ emission is generally weaker than $N_2H^+$. Still, our results suggest that $N_2D^+$ can in principle be used as an alternative or complement to $N_2H^+$ for constraining the CO snow line.

### 5.3. Outer Emission Rings of DCN, $N_2D^+$, and DCO$^+$

In this section, we discuss ring structures seen in the radial emission profiles of the deuterated molecules DCN, $N_2D^+$, and DCO$^+$. In Figure 14, we show DCN 3–2 and $N_2D^+$ 3–2 from this work (0.3 resolution) together with DCO$^+$ 4–3 for GM Aur (PI C. Qi, project code 2015.1.00678.S, 0.35 resolution) and DCO$^+$ 3–2 for the other disks (Huang et al. 2017; Favre et al. 2019, 0.26–0.7 resolution). The DCO$^+$ emission profiles were generated from the archival data in the same way as for the MAPS data, that is, by azimuthal averaging of a zeroth moment map generated with a Keplerian mask (Law et al. 2021a).

The disks around the T Tauri targets IM Lup, GM Aur, and AS 209 all show outer rings (i.e., a ring beyond the most central emission component) in at least two of the three deuterated molecules. IM Lup shows a double-ring structure in all three deuterated molecules. For GM Aur, both DCN and DCO$^+$ have an outer ring at ~300 au. For AS 209, $N_2D^+$ and DCO$^+$ show an outer ring at ~150 au. In addition, by employing a deconvolution technique, Flaherty et al. (2017) revealed an outer DCO$^+$ ring at 215 au$^{27}$ in HD 163296, in good agreement with the shoulder of the $N_2D^+$ emission profile (see also Salinas et al. 2018). Interestingly, all of these outer rings coincide with the edge of the millimeter dust continuum, which is indicated by the vertical gray line in Figure 14. There is also a shoulder in the $N_2D^+$ profile of MWC 480 just inside the millimeter dust edge.

The MAPS disks are not the only sources where rings of deuterated molecules are observed. A DCN ring at the millimeter dust edge is also seen in the disk around LkCa 15 (Huang et al. 2017). Furthermore, nondeuterated molecules can show similar associations to the dust edge: Law et al. (2021a) find such behavior among some of the MAPS sources for HCN and $H_2CO$. Carney et al. (2017), Pegues et al. (2020), and Guzmán et al. (2021) also report examples of enhanced $H_2CO$ emission at the dust edge.

Öberg et al. (2015) suggested that the outer DCO$^+$ ring in the IM Lup disk results from nonthermal desorption of CO ice by UV photons at the edge of the dust disk where the UV penetration depth increases. The newly supplied CO gas then reacts with $H_2D^+$, which is abundant due to the low temperature in these outer regions, to form DCO$^+$. In addition, CO might return to the gas phase at these large radii as a result of a thermal inversion: generic disk models by Cleeves (2016) and Facchini et al. (2017) predict that the removal of large grains due to radial drift leads to an increase in dust temperature at the dust disk edge.

Similar to DCO$^+$ forming from photodesorbed CO, an $N_2D^+$ ring might be formed from photodesorbed $N_2$ (Öberg et al. 2009) that reacts with $H_2D^+$. To get a ring in both $N_2D^+$ and DCO$^+$, the CO gas phase abundance should be neither too high (otherwise $N_2D^+$ gets destroyed) nor too low (otherwise DCO$^+$ will not form efficiently, e.g., Pagani et al. 2011). As for DCN, in cold gas it forms from $DCNH^+ + e^- \rightarrow DCN + H$. $DCNH^+$

<footnote>$^{27}$ We scaled the radial distance of the ring to account for the different distance of HD 163296 assumed in our study (101 pc) compared with Flaherty et al. (2017, 122 pc).</footnote>
is produced by the reaction of HNC with either DCO$^+$ or D$_3^+$ (Willacy 2007). Like the CO required for DCO$^+$ formation, the HNC might originate from UV photodesorption from icy grains or alternatively be formed from other species (e.g., HCN) that are photodesorbed. No direct laboratory measurements of HNC or HCN photodesorption are available, but naively we might expect that the photodesorption yield is roughly similar to other molecules such as CO, H$_2$O, or N$_2$ for which those measurements exist (Cuppen et al. 2017). A thermal inversion could also prompt the desorption of N$_2$ and other precursor molecules of the deuterated species discussed here.

Not all deuterated molecules show rings at the dust disk edge in all five disks. For example, there is no DCN ring at the millimeter dust edges of AS 209 and MWC 480. This suggests that the mechanism(s) for forming deuterated molecules at the dust disk edge depend on disk properties.

## 6. Summary

In this work, we study the distribution of the deuterated molecules DCN and N$_2$D$^+$ in five protoplanetary disks observed as part of the ALMA MAPS Large Program. We combine the observations of DCN 3–2 and N$_2$D$^+$ 3–2 with additional observations of HCN 1–0 and 3–2 and H$^{13}$CN 1–0 from MAPS and N$_2$H$^+$ 3–2 from the archive. We model azimuthally averaged spectra (including hyperfine structure) to derive radial column density profiles, from which we determine the profiles of the DCN/HCN and N$_2$D$^+/N_2H^+$ column density ratios. Our results are summarized as follows:

1. DCN 3–2 is detected in all five sources, with the detection toward GM Aur reported for the first time. N$_2$D$^+$ 3–2 is detected in four sources (IM Lup, AS 209, HD 163296, and MWC 480), two of which (IM Lup and MWC 480) are newly reported detections. This brings the total number of protoplanetary disks with N$_2$D$^+$ detected to four. In addition, we tentatively detect N$_2$D$^+$ 3–2 in GM Aur.

2. We find a variety of DCN emission profile morphologies. This hints at a complex dependence of HCN deuteration on individual disk properties such as the temperature or ionization structure.

3. N$_2$D$^+$ is expected to be abundant in the cold midplane beyond the CO snow line. This is reflected in the observed radial emission profiles that all show ring-like structures. For IM Lup, we discover a double ring that corresponds well to double rings seen for DCN and DCO$^+$. This is the first time the radial distribution of the N$_2$D$^+$ column density is derived in protoplanetary disks.

4. The DCN/HCN ratio shows a wide range of values across our sample, extending from a few times $10^{-3}$ in the inner regions of HD 163296 and MWC 480 to ~0.2 in the outer ring of IM Lup. We observe a general trend of decreasing DCN/HCN when moving from the outer disk regions toward the disk center. The strongest decrease of almost two orders of magnitude is seen for the disks around the Herbig stars HD 163296 and MWC 480. We interpret this trend as a confirmation that deuteration proceeds via in situ exchange reactions, because the backward endothermic reactions should become more efficient at the higher temperatures expected in the inner disk.

5. Even for temperatures $\ge 30$ K, we find DCN/HCN column density ratios above $10^{-2}$, which cannot be explained by the

low-temperature deuteration channel alone. We conclude that both the low- and high-temperature deuteration channels are active.

6. The derived DCN/HCN column density ratios are mostly larger than observed for comet Hale-Bopp (Meier et al. 1998). This suggests that Hale-Bopp either originated from smaller radii than probed here (~30 au) or inherited its HCN from the ice of the prestellar cloud.

7. We find N$_2$D$^+/N_2H^+$ column density ratios typically between a few times $10^{-2}$ and 1, consistent with model predictions of efficient N$_2$H$^+$ deuteration in the cold midplane of the outer disk. The disks around IM Lup and HC 163296 show a relatively flat N$_2$D$^+/N_2H^+$ profile, while the ratio strongly increases outward between ~50 and ~170 au for AS 209.

8. The inner edge of the observed N$_2$D$^+$ rings is consistent with the snow line locations inferred from N$_2$H$^+$ observations (Qi et al. 2015, 2019) and predicted by the MAPS reference models (Zhang et al. 2021). This suggests that N$_2$D$^+$ might be a useful tracer of the CO snow line. However, a careful continuum subtraction is essential to derive the inner edge of N$_2$D$^+$ 3–2 (see Appendix B).

9. All disks show rings or shoulders in the emission profiles of at least one of the deuterated molecules DCN, N$_2$D$^+$, and DCO$^+$ at the edge of the millimeter dust disks. These rings are probably linked to an increased UV flux due to the decreasing dust optical depth, with UV photodesorption releasing CO, N$_2$, and other precursor molecules. In addition, a thermal inversion at the dust disk edge might prompt thermal desorption of ices that results in the observed rings.

We would like to thank the anonymous referee for an exceptionally careful review of our manuscript that helped us to clarify many aspects of our paper. We are grateful to Satoshi Yamamoto and Yoko Oya for helpful discussions on spectroscopic properties of the observed lines. We thank Cécile Favre for sharing the DCO$^+$ data of AS 209. We also thank Marie-Lise Dubernet and Yaye Awa Ba for help with SPECTCOL and Holger S. P. Müller for help with the CDMS. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.01055.L, ADS/JAO.ALMA#2015.1.00678.S, ADS/JAO.ALMA#2012.1.00681.S, ADS/JAO.ALMA#2015.1.00657.S, ADS/JAO.ALMA#2013.1.00226.S, ADS/JAO.ALMA#2015.1.00486.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This research has made use of NASA's Astrophysics Data System and the SIMBAD database, operated at CDS, Strasbourg, France.

G.C. is supported by the NAOJ ALMA Scientific Research grant code 2019-13B. Y.Y. is supported by IGPEES, WINGS Program, the University of Tokyo. Y.A. acknowledges support by NAOJ ALMA Scientific Research grant code 2019-13B and Grant-in-Aid for Scientific Research (S) 18H05222, and Grant-in-Aid for Transformative Research Areas (A) 20H05844 and 20H05847. J.B.B., J.H., I.C., K.R.S., and K.Z. acknowledge the support of NASA through Hubble Fellowship grants HST-HF2-51429.001-A, HST-HF2-51460.001-A, HST-HF2-51405.001-A, HST-HF2-51419.001, and HST-HF2-51401.001, awarded by the
Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. J.H. and S.M.A. acknowledge funding support from the National Aeronautics and Space Administration under grant No. 17-XRP17 2-0012 issued through the Exoplanets Research Program. V.V.G. acknowledges support from FONDECYT Iniciación 11180904 and ANID project Basal AFB-170002. E.A.B. and A.D.B. acknowledge support from NSF AAG grant No. 1907653. L.I.C. gratefully acknowledges support from the David and Lucille Packard Foundation and Johnson & Johnson's WiSTEM2D Program. A.S.B. acknowledges the studentship funded by the Science and Technology Facilities Council of the United Kingdom (STFC). J. D.I. acknowledges support from STFC under ST/T000287/1. C. J.L. acknowledges funding from the National Science Foundation Graduate Research Fellowship under grant No. DGE1745303. R. L.G. acknowledges support from a CNES fellowship grant. Y.L. acknowledges the financial support by the Natural Science Foundation of China (grant No. 11973090). F.M. acknowledges support from ANR of France under contract ANR-16-CE31-0013 (Planet-Forming-Disks) and ANR-15-IDEX-02 (through CDP "Origins of Life"). F.L. and R.T. acknowledge support from the Smithsonian Institution as Submillimeter Array (SMA) Fellows. H.N. acknowledges support by NAOJ ALMA Scientific Research grant code 2018-10B and Grant-in-Aid for Scientific Research 18H05441. K.I.Ö. acknowledges support from the Simons Foundation (SCOL #321183) and an NSF AAG Grant (#1907653). T.T. is supported by JSPS KAKENHI grant Nos. JP17K14244 and JP20K04017. C.W. acknowledges financial support from the University of Leeds, STFC, and UKRI (grant Nos. ST/R000549/1, ST/T000287/1, MR/T040726/1). K.Z. acknowledges the support of the Office of the Vice Chancellor for Research and Graduate Education at the University of Wisconsin—Madison with funding from the Wisconsin Alumni Research Foundation.

*Facility: ALMA.*

*Software: astropy (Astropy Collaboration et al. 2013; Price-Whelan et al. 2018), bettermoments (Teague & Foreman-Mackey 2018), CASA (McMullin et al. 2007), CBcolors.py (https://gist.github.com/thriveth/8560036), corner.py (Foreman-Mackey 2016), emcee (Foreman-Mackey et al. 2013), gofish (Teague 2019), matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011), python-radex (https://github.com/gica3618/pythonradex), SciPy (Virtanen et al. 2020), SPECTCOL (VAMDC Consortium, http://www.vamdc.org), VISIBLE (Loomis et al. 2018).

# Appendix A
## Brief Summary of the JvM Correction

We first briefly summarize the CLEAN algorithm in order to set the stage for the JvM correction. CLEAN first initializes a “residual image” equal to the dirty image (i.e., the Fourier transform of the calibrated visibilities) and a “CLEAN model” equal to a blank image. Then a CLEAN component is placed into the CLEAN model at the position of the current peak in the residual map. At the same time, the newly placed CLEAN component is deconvolved from the residual map by subtracting its convolution with the dirty beam. This procedure is repeated until a stopping criterion is reached, which in the case of MAPS is a threshold of the peak S/N in the residual map. Then the CLEAN model is convolved with the CLEAN beam, which is a Gaussian fit to the dirty beam. Lastly, the residual

image is added to the convolved CLEAN model to obtain the final image cube. This last step is necessary because the residual image still contains astrophysical flux below the noise threshold. A more detailed discussion as well as visualization of the CLEAN algorithm is presented by Czekala et al. (2021).

A problem can arise in the final step of the workflow described above: the CLEAN model has units of flux/(CLEAN beam), while the residual map has units of flux/(dirty beam). Adding together these two maps produces a final image with inconsistent flux units. This is not a concern as long as the difference between the dirty and CLEAN beams is small. However, the combination of visibilities from short and long baseline antenna configurations (such as for MAPS) can lead to a significantly non-Gaussian dirty beam and therefore an inconsistent flux scale. Thus, MAPS employed a correction first discussed by Jorsater & van Moorsel (1995, hereafter JvM correction) and also used by, for example, Walter & Brinks (1999), Walter et al. (2008), Koda et al. (2019), and Pinte et al. (2020). MAPS implemented the correction by changing the units of the residual map to flux/(CLEAN beam) prior of combining it with the CLEAN model (Czekala et al. 2021). This is achieved by scaling the residual map by the ratio of the beam areas,

$$ \epsilon = \frac{\Omega_{\text{CLEAN}}}{\Omega_{\text{dirty}}} \qquad (A1) $$

The more $\epsilon$ deviates from 1, the more severe the mismatch between the CLEAN and dirty beam is, and the more important the JvM correction becomes. We list the $\epsilon$ factors for the image cubes used in this study in Tables 1 and 2.

The noise level of the final image cube is set by the residual map. Without the JvM correction, the noise level in the final cube is wrong by a factor $\epsilon$. Furthermore, for weak lines (such as those from DCN or $N_2D^+$), a significant fraction of the total flux is below the noise threshold and thus resides in the residual map. For those cases, the JvM correction is crucial to determine the correct line flux. As an illustration, Czekala et al. (2021) consider a channel of DCN 3–2 toward AS 209. About half of the total flux resides in the residual map. With an $\epsilon = 0.36$, this implies a ~50% difference between the total fluxes inferred from the corrected and uncorrected images, respectively.

# Appendix B
## Negative Emission of $N_2D^+$ 3–2

The radial profiles of the $N_2D^+$ 3–2 emission are strongly negative inside of ~50 au for HD 163296 and MWC 480. The negative emission is much deeper than the statistical error bars (see Figure 4), implying that this feature is highly unlikely to be just random noise. Furthermore, the $N_2D^+$ 3–2 profile for HD 163296 derived by Salinas et al. (2017) from an independent ALMA data set shows similarly negative emission, further suggesting that this is a systematic effect. The fact that both sources show strong continuum emission within ~50 au (Huang et al. 2018; Liu et al. 2019; Sierra et al. 2021) suggests that the effect is related to the continuum subtraction (Salinas et al. 2017). The difficulties with the continuum subtraction of $N_2D^+$ 3–2 are probably due to the fact that the line is located at the edge of a strong atmospheric feature, which may affect the bandpass calibration.

We inspected $N_2D^+$ spectra extracted from the central disk region from a data cube where the continuum has not been
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Radial emission profiles of $N_2D^+$ 3–2 for different orders of the polynomials used to subtract the continuum in the *uv* plane. In this example, the polynomials were only fitted over wavelengths without line emission. Strong variation is seen in the inner region of GM Aur, HD 163296, and MWC 480, depending on the adopted fit order. Similar variations are seen when the polynomial fit is performed using the whole spectral window.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Spectra of $N_2D^+$ 3–2 without continuum subtraction toward MWC 480 for different EBs, obtained by integrating over the central 0.′′3 region. The observations consist of a total of eight EBs, three with short baselines (i.e., a compact array configuration) and five with long baselines (i.e., an extended array configuration, Öberg et al. 2021). Spectra are vertically offset for clarity.</figcaption></figure>

subtracted. We found that the continuum indeed shows a wave-like structure, making accurate continuum subtraction difficult. The calibration procedure adapted by the MAPS team uses the CASA task *uvcontsub* to subtract the continuum in the *uv* plane by fitting first-order polynomials to those wavelengths without line emission (Öberg et al. 2021). For the $N_2D^+$ 3–2 data, this results in an overprediction of the continuum level at the wavelength of the line. In order to investigate the effect of the continuum subtraction on the radial emission profiles, we subtracted the continuum using *uvcontsub* with polynomials of order 1 to 4. We tested two cases: (1) fitting the polynomial to wavelengths outside of the line and (2) fitting the polynomial over all wavelengths of the spectral window, including the channels were line emission is expected. The latter approach is motivated by the fact that the $N_2D^+$ 3–2 emission is weak and might not influence the fit significantly. This might allow a more precise continuum subtraction at the wavelengths of the

line emission. For both cases, we find that in the inner ~50 au of MWC 480 and HD 163296 (and also GM Aur), the radial profiles strongly vary depending on the adopted fit order and the wavelength range chosen for the fit (Figure 15). On the other hand, in the outer region, the profiles are independent of the adopted fit order. This supports the hypothesis that the negative emission is due to oversubtraction of the continuum. Importantly, it also tells us that the radial profile is not affected by this problem beyond ~50 au.

For MWC 480, we also investigated whether the problem comes from a single execution block (EB), for example, due to poor weather. We imaged the continuum-subtracted visibilities for each EB separately and derived the radial emission profiles. We did not find a correlation between the strength of the negative emission and the precipitable water vapor for individual EBs, suggesting that variable weather conditions are not the primary cause of the issue. However, we find that
the short baseline (SB) executions show the strongest negative emission. We then imaged the visibilities without the continuum subtracted for each EB separately. Figure 16 shows spectra of MWC 480 extracted from the innermost circular region (~beam size) of these image cubes for each EB. The spectra from the SB executions indeed show the most obvious wavy structure of the continuum. We experimented with continuum subtractions tailored to individual EBs (e.g., excluding the high-velocity region of SB1 and SB2 for the continuum fit; see Figure 16). However, the inner region of the emission profile still showed strong variations depending on the adopted fit order and fitted wavelength range.

In summary, our tests suggest that nonoptimal continuum subtraction makes the N<sub>2</sub>D<sup>+</sup> 3–2 radial profiles of HD 163296 and MWC 480 unreliable within the inner ~50 au. The disk beyond ~50 au is not affected by this problem, and therefore our main results still hold. However, the negative emission means that the inner edge of the N<sub>2</sub>D<sup>+</sup> 3–2 emission is difficult to determine, which could reduce its usefulness as a CO snow line tracer.

visibilities are calculated from a Keplerian disk model or from a model based on another, strong emission line that is believed to have a spatial distribution similar to the weak line. The model visibilities are then cross-correlated with the observed visibilities, producing a response for each channel of the data. We applied this procedure for a series of Keplerian disk models with varying radial extents using the VISIBLE software (Loomis et al. 2018). In Figure 17 we present the matched filter responses for DCN 3–2, H<sup>13</sup>CN 1–0, and N<sub>2</sub>D<sup>+</sup> 3–2. For each disk and emission line, we selected the Keplerian mask that maximizes the filter response at the systemic velocity. DCN is detected above 10σ in all sources. N<sub>2</sub>D<sup>+</sup> is detected above 5σ in all sources except GM Aur. This is consistent with our analysis in the image plane. There is a marginal N<sub>2</sub>D<sup>+</sup> detection in GM Aur at 4σ. However, several peaks of similar magnitude in the spectrum cast some doubt on its significance, so we consider it only tentatively detected. H<sup>13</sup>CN is tentatively detected at 3σ for GM Aur and not detected in any of the other disks.

**Appendix C**
Matched Filter Analysis in the uv Plane

Loomis et al. (2018) proposed the application of a matched
filter in the uv plane to detect weak lines. With this method,

We recall that the first step for calculating azimuthally
averaged spectra is to shift the spectrum at each spatial position

**Appendix D**
Azimuthally Averaged Spectra

D.1. *Calculation of Errors*

<figure><img src="image_20.png" /><figcaption>**Figure 17.** Matched filter response for DCN 3–2, H<sup>13</sup>CN 1–0, and N<sub>2</sub>D<sup>+</sup> 3–2. For each disk and emission line, the filter with the Keplerian mask that maximizes the S/N was chosen, as indicated in the upper right of each panel. The x-axis is centered on the systemic velocity. The y-axis is in units of the noise in the filter response.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 18. Azimuthally averaged spectra of HCN 1–0, HCN 3–2, DCN 3–2, and H$^{13}$CN 1–0 toward IM Lup. Spectra are vertically offset for clarity and centered on the systemic velocity. If a scaling is applied to the spectra of a transition, it is indicated at the top of the column. The shaded region marks the 1σ error. The horizontal dotted line marks the zero flux level. For each spectrum, 20 randomly selected models from the MCMC chain (with the selection probability proportional to the model’s posterior probability) are shown by the blue lines (LTE) and the red lines (non-LTE), respectively. For most spectra, the models overlap closely so that only the blue lines are visible.</figcaption></figure>

of the data cube according to the projected Keplerian velocity (see Section 3.4). For each velocity $v_i$, the averaged flux is then computed as an azimuthal average over a deprojected ring within channel $i$ of the (shifted) cube. The error $\epsilon_i$ on the averaged flux is given by

$$ \epsilon_i = \frac{\sigma_i}{\sqrt{N_{\text{indep}}^i}}, \qquad (\text{D1}) $$

where $\sigma_i$ is the standard deviation of the pixels that were averaged and $N_{\text{indep}}^i$ is the number of independent samples
<figure><img src="image_22.png" /><figcaption>Figure 18. (Continued.)</figcaption></figure>

represented by those pixels. It is calculated as

$$N_{\text{indep}}^{i} = \sum_{k} \frac{1}{n_{\text{corr}}^{i,k}}, \qquad (\text{D2})$$

where the index k runs over all pixels included in the average and $n_{corr}^{i,k}$ is the number of pixels correlated with pixel k (where only pixels included in the average are considered). For example, consider the case where all pixels within a beam are

correlated. Then $n_{corr}^{i,k} = N_{beam}$ where $N_{beam}$ is the number of pixels per beam. In that case, $N_{indep}^{i} = N/N_{beam}$ with N being the total number of pixels included in the average. In other words, $N_{indep}^{i}$ equals the number of beams included in the average. However, in general $n_{corr}^{i,k} < N_{beam}$ because (1) the deprojected ring over which the average runs can be narrower than the beam, so that for a given pixel, some neighboring pixels are not included in the average, even though they are
<figure><img src="image_23.png" /><figcaption>Figure 19. Same as Figure 18, but for GM Aur.</figcaption></figure>

within a beam distance, and (2) the spectral shifting procedure applied to the data cube prior to azimuthally averaging means that pixels within a beam can be uncorrelated if their relative spectral shift exceeds the spectral resolution (see also Section 2.2 in Yen et al. 2016). Therefore, for a given pixel $k$, we calculate the number of correlated pixels as the number of pixels that satisfy all of the following conditions: they are in the same channel of the shifted data cube, are within the

deprojected ring, are within a beam centered onto pixel $k$, and have a spectral shift relative to pixel $k$ smaller than the spectral resolution.

As an alternative, we also computed the standard deviation of the azimuthally averaged spectrum measured at velocities without line emission, which we shall denote $\delta$. We find that generally, $\delta \approx \epsilon_i$. Note that $\epsilon_i$ depends on velocity, while $\delta$ does not. Most of the time, $\epsilon_i$ is slightly larger than $\delta$. The exception is the innermost disk region, where $\epsilon_i$ at times is significantly
<figure><img src="image_24.png" /><figcaption>**Figure 20.** Same as Figure 18, but for AS 209. Note the hyperfine structure of HCN 1–0 and 3–2.</figcaption></figure>

smaller than $\delta$. This points to an important limitation of estimating the error bar with $\epsilon_i$: when most averaged pixels are correlated, $N_{\text{indep}}^i \approx 1$, and thus $\epsilon_i \approx \sigma_i$. But $\sigma_i$ is small if most pixels included in the average are correlated, leading to an underestimation of the error bar. This problem only occurs in the innermost disk where all pixels within the deprojected ring under consideration are part of the same beam. To be conservative, we adopt $\max(\epsilon_i, \delta)$ as our final error bar.

## D.2. Gallery of Azimuthally Averaged Spectra

Figures 18–22 show the azimuthally averaged spectra of HCN 1–0, HCN 3–2, DCN 3–2, and H$^{13}$CN 1–0 as described in Section 3.4, together with LTE and non-LTE model spectra. Figures 23–27 show the corresponding spectra for N$_2$H$^+$ and N$_2$D$^+$, together with model spectra where $T_{\text{ex}} = 20$ K and $T_{\text{ex}} = T_{\text{mid}}$. In Figure 28, we show the azimuthally averaged spectra over the whole disk for N$_2$D$^+$ 3–2. The line is clearly
<figure><img src="image_25.png" /><figcaption>**Figure 21.** Same as Figure 18, but for HD 163296.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Figure 21. (Continued.)</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Figure 21. (Continued.)</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Figure 22. Same as Figure 18, but for MWC 480.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>**Figure 23.** Azimuthally averaged spectra for IM Lup of the 3–2 transitions of N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup>. Spectra are vertically offset for clarity and centered on the systemic velocity. If a scaling is applied to the spectra of a transition, it is indicated at the top of the column. The shaded region marks the 1σ error. The horizontal dotted line marks the zero flux level. For each spectrum, 20 randomly selected models from the MCMC chain (with the selection probability proportional to the model’s posterior probability) are shown by the blue lines (T<sub>ex</sub> = 20 K) and red lines (T<sub>ex</sub> = T<sub>mid</sub>), respectively. For most spectra, the models overlap closely so that only the blue lines are visible.</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>Figure 24. Same as Figure 23, but for GM Aur.</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>**Figure 25.** Same as Figure 23, but for AS 209.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>**Figure 26.** Same as Figure 23, but for HD 163296.</figcaption></figure>
<figure><img src="image_33.png" /><figcaption>**Figure 27.** Same as Figure 23, but for MWC 480.</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>**Figure 28.** Disk-integrated spectra of N<sub>2</sub>D<sup>+</sup> 3–2 centered on the systemic velocity. The radial integration range is indicated in the title of each panel. The shaded region denotes the 1σ uncertainty.</figcaption></figure>

detected toward IM Lup, AS 209, HD 163296, and MWC 480 and tentatively detected toward GM Aur. This is consistent with the disk-integrated flux (Section 3.1) and the results of the matched filter analysis (Appendix C).

## Appendix E
### Derivation of HCN and DCN Column Densities

In this appendix, we describe how the HCN and DCN column densities were derived by simultaneously fitting the shifted and azimuthally averaged spectra of HCN 1–0 and 3–2, H<sup>13</sup>CN 1–0, and DCN 3–2. An analogous approach was used to derive the column densities of N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup>.

#### E.1. LTE

Here we consider the column density calculation under the LTE assumption where the excitation temperature T<sub>ex</sub> is the same for all transitions of a molecule and equals the gas kinetic temperature T<sub>kin</sub>. We further assume that the different molecular species are sufficiently spatial to share the same temperature. Our model has the following free parameters:

1. The excitation temperature T<sub>ex</sub>.

2. The logarithms of the column densities: log<sub>10</sub>N<sub>HCN</sub> and log<sub>10</sub>N<sub>DCN</sub>. The column density of H<sup>13</sup>CN is fixed to the ISM ratio of <sup>13</sup>C/<sup>12</sup>C = 1/68 of the HCN column density.

3. The FWHM of the Gaussian convolution kernel used to model instrumental broadening of the lines (see below): FWHM<sub>B3</sub> (for Band 3 data, i.e., the 1–0 transitions) and FWHM<sub>B6</sub> (for Band 6 data, i.e., the 3–2 transitions).

4. For each line, a velocity offset with respect to the systemic velocity: Δv<sub>HCN1–0</sub>, Δv<sub>HCN3–2</sub>, Δv<sub>DCN3–2</sub>, and Δv<sub>H13CN1–0</sub>.

This sums up to nine free parameters. For given values of the parameters, model spectra (including hyperfine structure) are calculated by following an approach described in Appendix A of Teague & Loomis (2020). These model spectra can then be compared to the data. The intrinsic emission spectrum is given by

$$I_{\nu}(v) = (B_{\nu}(T_{\text{ex}}) - B_{\nu}(T_{\text{CMB}})) \cdot (1 - e^{-\tau_{\nu}}), \quad (\text{E1})$$

with B<sub>ν</sub> the Planck function and T<sub>CMB</sub> = 2.73 K the cosmic microwave background (CMB) temperature (i.e., the background radiation is assumed to be the CMB). The frequency-dependent optical depth τ<sub>ν</sub> is given by

$$\tau_{\nu} = \tau_0 \sum_i r_i \exp \left( - \frac{(v - v_{\text{sys}} - \delta v_i - \Delta v_X)^2}{2\sigma_v^2} \right), \quad (\text{E2})$$

with v<sub>sys</sub> the systemic velocity (taken from Öberg et al. 2021) and Δv<sub>X</sub> the free parameter describing the velocity offset of the observed emission line X. The sum goes over all hyperfine components, with r<sub>i</sub> the strength and δv<sub>i</sub> the velocity offset of each component. The strength r<sub>i</sub> is given by

$$r_i = \frac{g_{u,i} A_{ul,i}}{\sum_j g_{u,j} A_{ul,j}}, \quad (\text{E3})$$

with g<sub>u</sub> the weight of the upper level and A<sub>ul</sub> the Einstein coefficient for spontaneous emission of the hyperfine line. Tables 6, 7, and 8 list the values of molecular parameters used
<table><caption>**Table 6**
Spectral Line Data for HCN 1–0 and 3–2, DCN 3–2 and H$^{13}$CN 1–0 Used in This Study</caption><thead><tr><th>Molecule</th><th>J</th><th>F</th><th>ν<sub>0</sub><br>(GHz)</th><th>g<sub>u</sub></th><th>A<sub>ul</sub><br>(s<sup>-1</sup>)</th><th>r<sub>i</sub><sup>a</sup></th><th>δν<sup>b</sup><br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td colspan="8" style="text-align: center; font-weight: bold;">Hyperfine-resolved transitions</td></tr><tr><td>HCN</td><td>1–0</td><td>0–1</td><td>88.6339357</td><td>1</td><td>2.41 × 10<sup>-5</sup></td><td>0.1111</td><td>-7.063</td></tr><tr><td>HCN</td><td>1–0</td><td>2–1</td><td>88.6318475</td><td>5</td><td>2.41 × 10<sup>-5</sup></td><td>0.5555</td><td>0.000</td></tr><tr><td>HCN</td><td>1–0</td><td>1–1</td><td>88.6304156</td><td>3</td><td>2.41 × 10<sup>-5</sup></td><td>0.3333</td><td>4.843</td></tr><tr><td>HCN</td><td>3–2</td><td>2–2</td><td>265.8885221</td><td>5</td><td>1.30 × 10<sup>-4</sup></td><td>0.0370</td><td>-2.354</td></tr><tr><td>HCN</td><td>3–2</td><td>2–3</td><td>265.8869793</td><td>5</td><td>3.71 × 10<sup>-6</sup></td><td>0.0011</td><td>-0.615</td></tr><tr><td>HCN</td><td>3–2</td><td>4–3</td><td>265.8864999</td><td>9</td><td>8.36 × 10<sup>-4</sup></td><td>0.4286</td><td>-0.074</td></tr><tr><td>HCN</td><td>3–2</td><td>3–2</td><td>265.8864339</td><td>7</td><td>7.43 × 10<sup>-4</sup></td><td>0.2963</td><td>0.000</td></tr><tr><td>HCN</td><td>3–2</td><td>2–1</td><td>265.8861886</td><td>5</td><td>7.02 × 10<sup>-4</sup></td><td>0.2000</td><td>0.277</td></tr><tr><td>HCN</td><td>3–2</td><td>3–3</td><td>265.8848912</td><td>7</td><td>9.29 × 10<sup>-5</sup></td><td>0.0370</td><td>1.739</td></tr><tr><td>DCN</td><td>3–2</td><td>2–2</td><td>217.240622</td><td>5</td><td>7.12 × 10<sup>-5</sup></td><td>0.0370</td><td>-2.876</td></tr><tr><td>DCN</td><td>3–2</td><td>2–3</td><td>217.239079</td><td>5</td><td>2.03 × 10<sup>-6</sup></td><td>0.0011</td><td>-0.747</td></tr><tr><td>DCN</td><td>3–2</td><td>4–3</td><td>217.238612</td><td>9</td><td>4.57 × 10<sup>-4</sup></td><td>0.4286</td><td>-0.102</td></tr><tr><td>DCN</td><td>3–2</td><td>3–2</td><td>217.238555</td><td>7</td><td>4.07 × 10<sup>-4</sup></td><td>0.2963</td><td>-0.024</td></tr><tr><td>DCN</td><td>3–2</td><td>2–1</td><td>217.238300</td><td>5</td><td>3.84 × 10<sup>-4</sup></td><td>0.2000</td><td>0.328</td></tr><tr><td>DCN</td><td>3–2</td><td>3–3</td><td>217.236999</td><td>7</td><td>5.08 × 10<sup>-5</sup></td><td>0.0370</td><td>2.124</td></tr><tr><td>H<sup>13</sup>CN</td><td>1–0</td><td>0–1</td><td>86.3422543</td><td>1</td><td>2.23 × 10<sup>-5</sup></td><td>0.11</td><td>-7.25</td></tr><tr><td>H<sup>13</sup>CN</td><td>1–0</td><td>2–1</td><td>86.3401666</td><td>5</td><td>2.23 × 10<sup>-5</sup></td><td>0.56</td><td>0.00</td></tr><tr><td>H<sup>13</sup>CN</td><td>1–0</td><td>1–1</td><td>86.3387352</td><td>3</td><td>2.23 × 10<sup>-5</sup></td><td>0.33</td><td>4.97</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Hyperfine-unresolved transitions</td></tr><tr><td>HCN</td><td>1–0</td><td></td><td>88.6318475</td><td>9</td><td>2.41 × 10<sup>-5</sup></td><td></td><td></td></tr><tr><td>HCN</td><td>3–2</td><td></td><td>265.8864339</td><td>21</td><td>8.36 × 10<sup>-4</sup></td><td></td><td></td></tr><tr><td>DCN</td><td>3–2</td><td></td><td>217.2385378</td><td>21</td><td>4.57 × 10<sup>-4</sup></td><td></td><td></td></tr><tr><td>H<sup>13</sup>CN</td><td>1–0</td><td></td><td>86.3401666</td><td>9</td><td>2.23 × 10<sup>-5</sup></td><td></td><td></td></tr></tbody></table>

**Notes.**

<sup>a</sup> Relative strength of hyperfine component.

<sup>b</sup> Velocity offset with respect to the hyperfine-unresolved transition.

in this work. $\tau_0$ is the optical depth of the hyperfine-unresolved transition at the line center, given by

$$ \tau_0 = \frac{hc}{4\pi} (x_l B_{lu} - x_u B_{ul}) N \frac{1}{\sqrt{2\pi} \sigma_v}, \quad (E4) $$

with N the column density. Here the Einstein B coefficients and the fraction of molecules in the upper or lower level ($x_u$ and $x_l$) refer to the hyperfine-unresolved transition. The latter are given by

$$ x_u = g_u \frac{e^{-E_u/(kT_{ex})}}{Z(T_{ex})}, \quad (E5) $$

$$ x_l = g_l \frac{e^{-E_l/(kT_{ex})}}{Z(T_{ex})}, \quad (E6) $$

where the partition function Z is taken from the CDMS database.<fnref n="28" />

In contrast to Teague & Loomis (2020), the width of the intrinsic spectrum is fixed to the thermal width, that is,

$$ \sigma_v = \sqrt{\frac{kT_{kin}}{m}}, \quad (E7) $$

where the kinetic temperature $T_{kin} = T_{ex}$ in LTE and m is the mass of the molecule. We ignore turbulent broadening given

recent results suggesting low turbulence (e.g., Teague et al. 2018b). Finally, the spectrum given by Equation (E1) is convolved with a Gaussian kernel with FWHM given by the free parameter FWHM$_{Bx}$. This models any broadening of the line by effects such as beam smearing or imprecise Keplerian shifting. This broadening is highest in the inner region of the disk and gradually becomes smaller at larger radii (see Figures 18–27).

Table 4 shows an overview of the free parameters. We explore the parameter space using the MCMC method implemented in the `emcee` package (Foreman-Mackey et al. 2013) for each radial bin separately, starting with the innermost bin. We run 200 walkers, each with a total of 5000 steps. The walkers of radial bin $n+1$ are initialized using the result from radial bin $n$. The first 2500 steps are discarded prior to analyzing the results. We choose flat priors over the ranges detailed in Table 4.

To take into account the correlation between neighboring wavelengths in the azimuthally averaged spectra, we rescale the error bars in the spectra by multiplying by the square root of the correlation length (Booth et al. 2017). We define the latter as $max(1, D_{corr})$, where $D_{corr}$ is the FWHM (in units of spectral pixels) of the autocorrelation of the spectrum. Before calculating the autocorrelation, we remove the region of the spectrum with line emission and then subtract the mean value such that the line-free spectrum has a mean value of zero. The correlation length is typically $\lesssim 3$ spectral pixels, meaning that

<footnote n="28">²⁸ https://cdms.astro.uni-koeln.de/cdms/portal/</footnote>
<table><caption>**Table 7**
Spectral Line Data of N₂H⁺ 3–2 Used in This Study</caption><thead><tr><th>J</th><th>F<sub>1</sub></th><th>F</th><th>ν<sub>0</sub><br>(GHz)</th><th>g<sub>u</sub></th><th>A<sub>ul</sub><br>(s<sup>-1</sup>)</th><th>r<sub>i</sub><sup>a</sup></th><th>δν<sup>b</sup><br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td colspan="8">Hyperfine-resolved transitions</td></tr><tr><td>3–2</td><td>2–2</td><td>1–2</td><td>279.5147240</td><td>3</td><td>4.88 × 10<sup>-5</sup></td><td>0.0018</td><td>-3.191</td></tr><tr><td>3–2</td><td>2–2</td><td>3–2</td><td>279.5145858</td><td>7</td><td>2.06 × 10<sup>-5</sup></td><td>0.0018</td><td>-3.043</td></tr><tr><td>3–2</td><td>2–2</td><td>1–1</td><td>279.5143988</td><td>3</td><td>9.32 × 10<sup>-5</sup></td><td>0.0035</td><td>-2.842</td></tr><tr><td>3–2</td><td>2–2</td><td>3–3</td><td>279.5143328</td><td>7</td><td>1.60 × 10<sup>-4</sup></td><td>0.0141</td><td>-2.771</td></tr><tr><td>3–2</td><td>2–2</td><td>2–2</td><td>279.5142143</td><td>5</td><td>1.82 × 10<sup>-4</sup></td><td>0.0115</td><td>-2.644</td></tr><tr><td>3–2</td><td>2–2</td><td>2–3</td><td>279.5139612</td><td>5</td><td>3.50 × 10<sup>-5</sup></td><td>0.0022</td><td>-2.373</td></tr><tr><td>3–2</td><td>2–2</td><td>2–1</td><td>279.5138890</td><td>5</td><td>3.42 × 10<sup>-5</sup></td><td>0.0022</td><td>-2.295</td></tr><tr><td>3–2</td><td>4–3</td><td>3–3</td><td>279.5123144</td><td>7</td><td>1.16 × 10<sup>-4</sup></td><td>0.0103</td><td>-0.606</td></tr><tr><td>3–2</td><td>3–2</td><td>2–2</td><td>279.5121160</td><td>5</td><td>1.18 × 10<sup>-4</sup></td><td>0.0074</td><td>-0.394</td></tr><tr><td>3–2</td><td>4–3</td><td>5–4</td><td>279.5118523</td><td>11</td><td>1.26 × 10<sup>-3</sup></td><td>0.1752</td><td>-0.111</td></tr><tr><td>3–2</td><td>2–1</td><td>1–1</td><td>279.5118447</td><td>3</td><td>4.95 × 10<sup>-4</sup></td><td>0.0188</td><td>-0.103</td></tr><tr><td>3–2</td><td>4–3</td><td>4–3</td><td>279.5118028</td><td>9</td><td>1.20 × 10<sup>-3</sup></td><td>0.1360</td><td>-0.058</td></tr><tr><td>3–2</td><td>4–3</td><td>3–2</td><td>279.5118022</td><td>7</td><td>1.14 × 10<sup>-3</sup></td><td>0.1008</td><td>-0.057</td></tr><tr><td>3–2</td><td>3–2</td><td>4–3</td><td>279.5117986</td><td>9</td><td>1.10 × 10<sup>-3</sup></td><td>0.1253</td><td>-0.053</td></tr><tr><td>3–2</td><td>3–2</td><td>2–1</td><td>279.5117908</td><td>5</td><td>9.68 × 10<sup>-4</sup></td><td>0.0612</td><td>-0.045</td></tr><tr><td>3–2</td><td>3–2</td><td>3–2</td><td>279.5116770</td><td>7</td><td>1.00 × 10<sup>-3</sup></td><td>0.0888</td><td>0.077</td></tr><tr><td>3–2</td><td>2–1</td><td>1–2</td><td>279.5116372</td><td>3</td><td>2.87 × 10<sup>-5</sup></td><td>0.0011</td><td>0.120</td></tr><tr><td>3–2</td><td>2–1</td><td>3–2</td><td>279.5114990</td><td>7</td><td>1.07 × 10<sup>-3</sup></td><td>0.0948</td><td>0.268</td></tr><tr><td>3–2</td><td>3–2</td><td>3–3</td><td>279.5114240</td><td>7</td><td>1.59 × 10<sup>-4</sup></td><td>0.0141</td><td>0.349</td></tr><tr><td>3–2</td><td>4–3</td><td>4–4</td><td>279.5114066</td><td>9</td><td>6.35 × 10<sup>-5</sup></td><td>0.0072</td><td>0.367</td></tr><tr><td>3–2</td><td>2–1</td><td>1–0</td><td>279.5113769</td><td>3</td><td>5.88 × 10<sup>-4</sup></td><td>0.0223</td><td>0.399</td></tr><tr><td>3–2</td><td>2–1</td><td>2–1</td><td>279.5113349</td><td>5</td><td>7.60 × 10<sup>-4</sup></td><td>0.0481</td><td>0.444</td></tr><tr><td>3–2</td><td>2–1</td><td>2–2</td><td>279.5111274</td><td>5</td><td>2.44 × 10<sup>-4</sup></td><td>0.0155</td><td>0.667</td></tr><tr><td>3–2</td><td>3–3</td><td>4–3</td><td>279.5102649</td><td>9</td><td>1.21 × 10<sup>-5</sup></td><td>0.0014</td><td>1.592</td></tr><tr><td>3–2</td><td>3–3</td><td>3–3</td><td>279.5098903</td><td>7</td><td>8.05 × 10<sup>-5</sup></td><td>0.0071</td><td>1.994</td></tr><tr><td>3–2</td><td>3–3</td><td>4–4</td><td>279.5098686</td><td>9</td><td>1.46 × 10<sup>-4</sup></td><td>0.0167</td><td>2.017</td></tr><tr><td>3–2</td><td>3–3</td><td>2–2</td><td>279.5098171</td><td>5</td><td>1.55 × 10<sup>-4</sup></td><td>0.0098</td><td>2.072</td></tr><tr><td>3–2</td><td>3–3</td><td>3–4</td><td>279.5094940</td><td>7</td><td>1.13 × 10<sup>-5</sup></td><td>0.0010</td><td>2.419</td></tr><tr><td colspan="8">Hyperfine-unresolved transition</td></tr><tr><td>3–2</td><td></td><td></td><td>279.5117491</td><td>63</td><td>1.26 × 10<sup>-3</sup></td><td></td><td></td></tr></tbody></table>

**Notes.**

a Relative strength of hyperfine component.

b Velocity offset with respect to the hyperfine-unresolved transition.

the rescaling of the error bars is marginal. However, we find that the correlation length can be much larger (up to ~20 spectral pixels) for the azimuthally averaged spectra extracted from the innermost beam. This is probably because the individual spectra that were azimuthally averaged are all spatially correlated. This spatial correlation is transformed to spectral correlation by the Keplerian shifting procedure.

The correlation length is difficult to determine for spectra extracted from inner-disk regions because the line emission becomes very broad. Thus, no line-free region is available. To circumvent this problem, for the innermost regions (within 1.5 beams), we adopt the correlation lengths determined from the weakest lines: H¹³CN 1–0 for HCN 1–0, HCN 3–2, and DCN 3–2, and N₂D⁺ 3–2 for N₂H⁺ 3–2.

As we move outward in radius, the upper boundaries of the priors on FWHM<sub>B3</sub> and FWHM<sub>B6</sub> are adjusted according to the following scheme. First, we calculate the S/N of FWHM<sub>Bx</sub>, defined as the median of the MCMC (“signal”) divided by the mean of the 16th and 84th percentiles (“error”). If the S/N is higher than a set threshold (we choose between 3 and 8), then the upper boundary of the prior for the next radial bin is set to the median plus five times the error. This procedure reflects the observed sharp decrease of the line width with increasing

radius and allows us to set tighter priors in the outer regions that show little signal. In the case where only the B6 or B3 FWHM exceeds the S/N required for adjustment, we still adjust both FWHM by using the higher S/N FWHM as a basis to estimate a conservative upper boundary for the lower S/N FWHM, taking into account the difference in channel size if necessary. Figure 29 shows an example of the fitted FWHM together with the radially adjusted prior boundaries.

To derive the column densities of N₂H⁺ and N₂D⁺, we used an analogous approach and simultaneously fitted the 3–2 transitions. However, the excitation temperature was fixed to 20 K. In addition, to test the dependence on the assumed temperature, we also run fits with the excitation temperature fixed to the midplane temperature extracted from the MAPS reference models by Zhang et al. (2021; see our Figure 30). These source-specific thermochemical models fit the spectral energy distribution, (sub)millimeter continuum images, and vertical locations of the CO-emitting layers reported by Law et al. (2021b). We note that the increase of temperature in the outer disk (particularly clearly seen for IM Lup at ~300 au) coincides with the truncation radius of the pebble disk, which makes the midplane warmer.
<figure><img src="image_35.png" /><figcaption>**Figure 29.** Example of the fitted FWHM of the Gaussian kernel used to convolve the model spectrum for MWC 480. The black line marks the median, while the shaded regions encompass the 16th to 84th, 2.3th to 97.7th, and 0.15th to 99.85th percentile regions. The red dashed lines show the upper and lower bound of the flat prior. Sequentially fitting larger and larger radii, the upper bound is adjusted by using the fit results from the previous radius (see text for details).</figcaption></figure>

<table><caption>**Table 8**

Spectral Line Data of N<sub>2</sub>D<sup>+</sup> 3–2 Used in This Study</caption><thead><tr><th>J</th><th>F<sub>1</sub></th><th>F</th><th>ν<sub>0</sub><br>(GHz)</th><th>g<sub>u</sub></th><th>A<sub>ul</sub><br>(s<sup>-1</sup>)</th><th>r<sub>i</sub><sup>a</sup></th><th>δν<sup>b</sup><br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td colspan="8" style="text-align: center; font-weight: bold;">Hyperfine-resolved transitions</td></tr><tr><td>3–2</td><td>2–2</td><td>1–2</td><td>231.3248131</td><td>3</td><td>2.76 × 10<sup>-5</sup></td><td>0.0019</td><td>-3.868</td></tr><tr><td>3–2</td><td>2–2</td><td>3–2</td><td>231.3246669</td><td>7</td><td>1.17 × 10<sup>-5</sup></td><td>0.0018</td><td>-3.679</td></tr><tr><td>3–2</td><td>2–2</td><td>1–1</td><td>231.3244843</td><td>3</td><td>5.25 × 10<sup>-5</sup></td><td>0.0035</td><td>-3.442</td></tr><tr><td>3–2</td><td>2–2</td><td>3–3</td><td>231.3244147</td><td>7</td><td>9.04 × 10<sup>-5</sup></td><td>0.0142</td><td>-3.352</td></tr><tr><td>3–2</td><td>2–2</td><td>2–2</td><td>231.3242977</td><td>5</td><td>1.03 × 10<sup>-4</sup></td><td>0.0115</td><td>-3.200</td></tr><tr><td>3–2</td><td>2–2</td><td>2–3</td><td>231.3240455</td><td>5</td><td>1.99 × 10<sup>-5</sup></td><td>0.0022</td><td>-2.873</td></tr><tr><td>3–2</td><td>2–2</td><td>2–1</td><td>231.3239689</td><td>5</td><td>1.94 × 10<sup>-5</sup></td><td>0.0022</td><td>-2.774</td></tr><tr><td>3–2</td><td>4–3</td><td>3–3</td><td>231.3223976</td><td>7</td><td>6.57 × 10<sup>-5</sup></td><td>0.0103</td><td>-0.738</td></tr><tr><td>3–2</td><td>3–2</td><td>2–2</td><td>231.3222005</td><td>5</td><td>6.64 × 10<sup>-5</sup></td><td>0.0074</td><td>-0.482</td></tr><tr><td>3–2</td><td>2–1</td><td>1–1</td><td>231.3219314</td><td>3</td><td>2.81 × 10<sup>-4</sup></td><td>0.0188</td><td>-0.134</td></tr><tr><td>3–2</td><td>4–3</td><td>5–4</td><td>231.3219281</td><td>11</td><td>7.14 × 10<sup>-4</sup></td><td>0.1758</td><td>-0.129</td></tr><tr><td>3–2</td><td>4–3</td><td>3–2</td><td>231.3218798</td><td>7</td><td>6.46 × 10<sup>-4</sup></td><td>0.1012</td><td>-0.067</td></tr><tr><td>3–2</td><td>4–3</td><td>4–3</td><td>231.3218794</td><td>9</td><td>6.78 × 10<sup>-4</sup></td><td>0.1365</td><td>-0.066</td></tr><tr><td>3–2</td><td>3–2</td><td>4–3</td><td>231.3218770</td><td>9</td><td>6.24 × 10<sup>-4</sup></td><td>0.1257</td><td>-0.063</td></tr><tr><td>3–2</td><td>3–2</td><td>2–1</td><td>231.3218717</td><td>5</td><td>5.49 × 10<sup>-4</sup></td><td>0.0614</td><td>-0.056</td></tr><tr><td>3–2</td><td>3–2</td><td>3–2</td><td>231.3217560</td><td>7</td><td>5.69 × 10<sup>-4</sup></td><td>0.0891</td><td>0.094</td></tr><tr><td>3–2</td><td>2–1</td><td>3–2</td><td>231.3215801</td><td>7</td><td>6.07 × 10<sup>-4</sup></td><td>0.0951</td><td>0.322</td></tr><tr><td>3–2</td><td>3–2</td><td>3–3</td><td>231.3215039</td><td>7</td><td>9.01 × 10<sup>-5</sup></td><td>0.0141</td><td>0.420</td></tr><tr><td>3–2</td><td>4–3</td><td>4–4</td><td>231.3214836</td><td>9</td><td>3.59 × 10<sup>-5</sup></td><td>0.0072</td><td>0.447</td></tr><tr><td>3–2</td><td>2–1</td><td>1–0</td><td>231.3214595</td><td>3</td><td>3.33 × 10<sup>-4</sup></td><td>0.0224</td><td>0.478</td></tr><tr><td>3–2</td><td>2–1</td><td>2–1</td><td>231.3214160</td><td>5</td><td>4.31 × 10<sup>-4</sup></td><td>0.0482</td><td>0.534</td></tr><tr><td>3–2</td><td>2–1</td><td>2–2</td><td>231.3212109</td><td>5</td><td>1.39 × 10<sup>-4</sup></td><td>0.0155</td><td>0.800</td></tr><tr><td>3–2</td><td>3–3</td><td>3–3</td><td>231.3199797</td><td>7</td><td>4.55 × 10<sup>-5</sup></td><td>0.0071</td><td>2.396</td></tr><tr><td>3–2</td><td>3–3</td><td>4–4</td><td>231.3199570</td><td>9</td><td>8.30 × 10<sup>-5</sup></td><td>0.0167</td><td>2.425</td></tr><tr><td>3–2</td><td>3–3</td><td>2–2</td><td>231.3199064</td><td>5</td><td>8.82 × 10<sup>-5</sup></td><td>0.0099</td><td>2.491</td></tr><tr><td colspan="8" style="text-align: center; font-weight: bold;">Hyperfine-unresolved transition</td></tr><tr><td>3–2</td><td></td><td></td><td>231.3218283</td><td>63</td><td>4.38 × 10<sup>-4</sup></td><td></td><td></td></tr></tbody></table>

**Notes.**

<footnote><sup>a</sup> Relative strength of hyperfine component.</footnote>

<footnote><sup>b</sup> Velocity offset with respect to the hyperfine-unresolved transition.</footnote>
<figure><img src="image_36.png" /><figcaption>**Figure 30.** Profiles of the midplane gas temperature $T_{\text{mid}}$. These profiles were extracted from the MAPS reference models (Zhang et al. 2021) and are used for fitting the $N_2H^+$ and $N_2D^+$ column densities under the assumption that $T_{\text{ex}} = T_{\text{mid}}$ (Section 4.2).</figcaption></figure>

Since the $N_2H^+$ and $N_2D^+$ data come from different observations, two separate parameters for the FWHM of the Gaussian kernels were used, one for each line. The free parameters and their priors are listed in Table 5.

The method described here is essentially the same as employed by Bergner et al. (2021) to calculate HCN and CN column densities and by Guzmán et al. (2021) to calculate HCN and $C_2H$ column densities. One difference is that Bergner et al. (2021) and Guzmán et al. (2021) do not convolve the intrinsic spectrum but simply redistribute the total flux of the intrinsic spectrum in a model spectrum with broader components. They also do not include $H^{13}\text{CN}$ in their fits of the HCN column density. Furthermore, Guzmán et al. (2021) employ a two-step fitting method, where they first fit HCN 1–0 and 3–2 at 0.′3 resolution to constrain the excitation temperature and then fit HCN 3–2 at higher 0.′15 resolution, using the constraints on $T_{\text{ex}}$ from the first step. The HCN column densities derived by Bergner et al. (2021) and Guzmán et al. (2021) are generally consistent with our results, with differences only seen inward of 50–100 au (Guzmán et al. 2021, Appendix C).

## E.2. Non-LTE

We repeated the calculation of the HCN and DCN column densities described in Appendix E.1 without the assumption of LTE, that is, $T_{\text{ex}}$ is not necessarily equal to $T_{\text{kin}}$ and can now take different values for each of the emission lines. We now fit the kinetic gas temperature $T_{\text{kin}}$ and the $H_2$ number density $n_{H_2}$ with priors detailed in Table 4. Instead of using Equations (E5) and (E6), the fractional-level populations are now determined by solving the statistical equilibrium of the excitation and deexcitation processes. To this end, we use the radiative transfer code `pythonradex`.<fnref n="29" /> For the purpose of the radiative transfer, we ignore the hyperfine structure: we assume that all hyperfine lines of a given emission line share the same excitation temperature (Teague & Loomis 2020) and we do not consider cross-excitation. The excitation temperature for a given transition is then computed from the fractional-level populations by using

$$ \frac{x_u}{x_l} = \frac{g_u}{g_l} e^{-\Delta E/(kT_{\text{ex}})}, \qquad (\text{E8}) $$

where $\Delta E = E_u - E_l$. The rest of the calculation is identical to the LTE case. In practice, we precompute a grid of fractional-level populations and excitation temperatures in the $T_{\text{kin}}-n_{H_2}$-space and interpolate on that grid during the MCMC run

for computational efficiency. We adopt the scaled HCN-He collision rates by Dumouchel et al. (2010) available in the LAMDA database<fnref n="30" /> for all three species (HCN, DCN, $H^{13}\text{CN}$). We find that the calculation of the non-LTE level populations fails for column densities above ~$10^{16} \text{ cm}^{-2}$. Thus, we reduce the upper boundary of the column density priors compared to the LTE case (Table 4). Because we find column densities below $10^{16} \text{ cm}^{-2}$ with both the LTE and non-LTE runs, this should not have a significant effect on our results.

## Appendix F: Molecular Data

We used spectral line data from the CDMS database (Müller et al. 2001, 2005; Endres et al. 2016), downloaded with the SPECTCOL software. CDMS gives the following sources for the data:

1. HCN: Ebenstein & Muenter (1984), Maki et al. (2000), Ahrens et al. (2002), Thorwirth et al. (2003), Lapinov (2006)

2. DCN: DeLeon & Muenter (1984), Möllmann et al. (2002), Brünken et al. (2004)

3. $H^{13}\text{CN}$: Winnewisser & Vogt (1978), Preusser & Maki (1993), Maki et al. (2000), Maiwald et al. (2000), Fuchs et al. (2004), Cazzoli & Puzzarini (2005)

4. $N_2H^+$: Verhoeve et al. (1990), Caselli et al. (1995), Amano et al. (2005), Cazzoli et al. (2012)

5. $N_2D^+$: Dore et al. (2004), Amano et al. (2005), Pagani et al. (2009).

Table 6 lists the line data for HCN 1–0 and 3–2, DCN 3–2, and $H^{13}\text{CN}$ 1–0. Tables 7 and 8 show corresponding data for $N_2H^+$ 3–2 and $N_2D^+$ 3–2, respectively.

## Appendix G: Data Available for Download

The MAPS data products (calibrated visibilities, image cubes, moment maps, radial emission profiles, and emission surfaces), as well as the scripts used to generate the data products are available for download from the ALMA archive<fnref n="31" /> and from the MAPS website.<fnref n="32" /> In addition, the following data, specific to this paper, are also available for download:

<footnote n="29">²⁹ https://github.com/gica3618/pythonradex</footnote>

<footnote n="30">³⁰ https://home.strw.leidenuniv.nl/~moldata/, accessed 2021 March 1.</footnote>

<footnote n="31">³¹ https://almascience.nrao.edu/alma-data/lp/maps</footnote>

<footnote n="32">³² http://www.alma-maps.info</footnote>
1. Image cubes of HCN 1-0 and 3-2 as well as H$^{13}$CN 1-0 covering all hyperfine transitions (Section 2.1.2).

2. Visibilities, (smoothed) image cubes, zeroth-moment maps, and radial emission profiles of the archival N$_2$H$^+$ 3-2 data (Sections 2.2, 3.1, and 3.3).

3. Azimuthally averaged spectra of the HCN, DCN, H$^{13}$CN, N$_2$H$^+$, and N$_2$D$^+$ emission lines analysed in this paper (Section 3.4).

4. Radial profiles of column density, optical depth, temperature, and deuteration fraction (Section 4).

5. Image cubes, zeroth-moment maps, and radial emission profiles of the archival DCO$^+$ data (Section 5.3).

6. Visibilities, image cubes, radial emission profiles, and associated scripts used for studying the negative emission of N$_2$D$^+$ 3-2 (Appendix B).

7. Disk-integrated spectra of N$_2$H$^+$ 3-2 (Appendix D.2).

8. Molecular data and partition functions (Appendix F).

## ORCID iDs

Gianni Cataldi https://orcid.org/0000-0002-2700-9676

Yoshihide Yamato https://orcid.org/0000-0003-4099-6941

Yuri Aikawa https://orcid.org/0000-0003-3283-6884

Jennifer B. Bergner https://orcid.org/0000-0002-8716-0482

Kenji Furuya https://orcid.org/0000-0002-2026-8157

Viviana V. Guzmán https://orcid.org/0000-0003-4784-3040

Jane Huang https://orcid.org/0000-0001-6947-6072

Ryan A. Loomis https://orcid.org/0000-0002-8932-1219

Chunhua Qi https://orcid.org/0000-0001-8642-1786

Sean M. Andrews https://orcid.org/0000-0003-2253-2270

Edwin A. Bergin https://orcid.org/0000-0003-4179-6394

Alice S. Booth https://orcid.org/0000-0003-2014-2121

Arthur D. Bosman https://orcid.org/0000-0003-4001-3589

L. Ismedore Cleeves https://orcid.org/0000-0003-2076-8001

Ian Czekala https://orcid.org/0000-0002-1483-8811

John D. Ilee https://orcid.org/0000-0003-1008-1142

Charles J. Law https://orcid.org/0000-0003-1413-1776

Romane Le Gal https://orcid.org/0000-0003-1837-3772

Yao Liu https://orcid.org/0000-0002-7616-666X

Feng Long https://orcid.org/0000-0002-7607-719X

François Ménard https://orcid.org/0000-0002-1637-7393

Hideko Nomura https://orcid.org/0000-0002-7058-7682

Karin I. Öberg https://orcid.org/0000-0001-8798-1347

Kamber R. Schwarz https://orcid.org/0000-0002-6429-9457

Richard Teague https://orcid.org/0000-0003-1534-5186

Takashi Tsukagoshi https://orcid.org/0000-0002-6034-2892

Catherine Walsh https://orcid.org/0000-0001-6078-786X

David J. Wilner https://orcid.org/0000-0003-1526-7587

Ke Zhang https://orcid.org/0000-0002-0661-7517

## References

Ahrens, V., Lewen, F., Takano, S., et al. 2002, *ZNatA*, 57, 669

Aikawa, Y., Cataldi, G., Yamato, Y., et al. 2021, *ApJS*, 257, 13

Aikawa, Y., Furuya, K., Hincelin, U., & Herbst, E. 2018, *ApJ*, 855, 119

Aikawa, Y., Furuya, K., Nomura, H., & Qi, C. 2015, *ApJ*, 807, 120

Aikawa, Y., Umebayashi, T., Nakano, T., & Miyama, S. M. 1999, *ApJ*, 519, 705

Alarcón, F., Bosman, A. D., Bergin, E. A., et al. 2021, *ApJS*, 257, 8

Alexander, C. M. O. 2017, *RSPTA*, 375, 20150384

Amano, T., Hirao, T., & Takano, J. 2005, *JMoSp*, 234, 170

Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, *ApJL*, 869, L41

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, *A&A*, 558, A33

Bergin, E., Calvet, N., D'Alessio, P., & Herczeg, G. J. 2003, *ApJL*, 591, L159

Bergin, E. A., Ciardi, D. R., Lada, C. J., Alves, J., & Lada, E. A. 2001, *ApJ*, 557, 209

Bergin, E. A., Du, F., Cleeves, L. I., et al. 2016, *ApJ*, 831, 101

Bergner, J. B., Öberg, K. I., Bergin, E. A., et al. 2019, *ApJ*, 876, 25

Bergner, J. B., Öberg, K. I., Bergin, E. A., et al. 2020a, *ApJ*, 896, 176

Bergner, J. B., Öberg, K. I., Bergin, E. A., et al. 2020b, *ApJ*, 898, 97

Bergner, J. B., Öberg, K. I., Guzmán, V. V., et al. 2021, *ApJS*, 257, 11

Bockelée-Morvan, D., Calmonte, U., Charnley, S., et al. 2015, *SSRv*, 197, 47

Bockelée-Morvan, D., Gautier, D., Lis, D. C., et al. 1998, *Icar*, 133, 147

Booth, M., Dent, W. R. F., Jordán, A., et al. 2017, *MNRAS*, 469, 3200

Bosman, A. D., Alarcón, F., Bergin, E. A., et al. 2021a, *ApJS*, 257, 7

Bosman, A. D., Bergin, E. A., Loomis, R. A., et al. 2021b, *ApJS*, 257, 15

Brasser, R., & Morbidelli, A. 2013, *Icar*, 225, 40

Brüken, S., Fuchs, U., Lewen, F., et al. 2004, *JMoSp*, 225, 152

Carney, M. T., Fedele, D., Hogerheijde, M. R., et al. 2018, *A&A*, 614, A106

Carney, M. T., Hogerheijde, M. R., Loomis, R. A., et al. 2017, *A&A*, 605, A21

Caselli, P., Myers, P. C., & Thaddeus, P. 1995, *ApJL*, 455, L77

Cazzoli, G., Cludi, L., Buffa, G., & Puzzarini, C. 2012, *ApJS*, 203, 11

Cazzoli, G., & Puzzarini, C. 2005, *JMoSp*, 233, 280

Ceccarelli, C., Caselli, P., Bockelée-Morvan, D., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 859

Chapillon, E., Guilloteau, S., Dutrey, A., Piétu, V., & Guélin, M. 2012, *A&A*, 537, A60

Cleeves, L. I. 2016, *ApJL*, 816, L21

Cleeves, L. I., Öberg, K. I., Wilner, D. J., et al. 2016, *ApJ*, 832, 110

Cleeves, L. I., Öberg, K. I., Wilner, D. J., et al. 2018, *ApJ*, 865, 155

Crovisier, J., Bockelée-Morvan, D., Colom, P., et al. 2004, *A&A*, 418, 1141

Cuppen, H. M., Walsh, C., Lamberts, T., et al. 2017, *SSRv*, 212, 1

Cyburt, R. H., Fields, B. D., Olive, K. A., & Yeh, T.-H. 2016, *RvMP*, 88, 015004

Czekala, I., Loomis, R. A., Teague, R., et al. 2021, *ApJS*, 257, 2

DeLeon, R. L., & Muenter, J. S. 1984, *JChPh*, 80, 6092

Dore, L., Caselli, P., Beninati, S., et al. 2004, *A&A*, 413, 1177

Dumouchel, F., Faure, A., & Lique, F. 2010, *MNRAS*, 406, 2488

Dutrey, A., Guilloteau, S., & Guelin, M. 1997, *A&A*, 317, L55

Ebenstein, W. L., & Muenter, J. S. 1984, *JChPh*, 80, 3989

Endres, C. P., Schlemmer, S., Schilke, P., Stutzki, J., & Müller, H. S. P. 2016, *JMoSp*, 327, 95

Facchini, S., Birnstiel, T., Bruderer, S., & van Dishoeck, E. F. 2017, *A&A*, 605, A16

Facchini, S., Teague, R., Bae, J., et al. 2021, *AJ*, 162, 99

Favre, C., Bergin, E. A., Cleeves, L. I., et al. 2015, *ApJL*, 802, L23

Favre, C., Fedele, D., Maud, L., et al. 2019, *ApJ*, 871, 107

Feng, S., Caselli, P., Wang, K., et al. 2019, *ApJ*, 883, 202

Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, *ApJ*, 843, 150

Foreman-Mackey, D. 2016, *JOSS*, 1, 24

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, *PASP*, 125, 306

Fuchs, U., Bruenken, S., Fuchs, G. W., et al. 2004, *ZNatA*, 59, 861

Furuya, K., Aikawa, Y., Hama, T., & Watanabe, N. 2019, *ApJ*, 882, 172

Gerner, T., Shirley, Y. L., Beuther, H., et al. 2015, *A&A*, 579, A80

Gicquel, A., Milam, S. N., Villanueva, G. L., et al. 2014, *ApJ*, 794, 1

Guzmán, V. V., Bergner, J. B., Law, C. J., et al. 2021, *ApJS*, 257, 6

Guzmán, V. V., Huang, J., Andrews, S. M., et al. 2018, *ApJL*, 869, L48

Hébrard, G., Tripp, T. M., Chayer, P., et al. 2005, *ApJ*, 635, 1136

Henning, T., & Semenov, D. 2013, *Chrv*, 113, 9016

Hily-Blant, P., Magalhaes de Souza, V., Kastner, J., & Forveille, T. 2019, *A&A*, 632, L12

Högbom, J. A. 1974, *A&AS*, 15, 417

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, *ApJL*, 869, L42

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2020, *ApJ*, 891, 48

Huang, J., & Öberg, K. I. 2015, *ApJL*, 809, L26

Huang, J., Öberg, K. I., & Andrews, S. M. 2016, *ApJL*, 823, L18

Huang, J., Öberg, K. I., Qi, C., et al. 2017, *ApJ*, 835, 231

Hunter, J. D. 2007, *CSE*, 9, 90

Jorsater, S., & van Moorsel, G. A. 1995, *AJ*, 110, 2037

Koda, J., Teuben, P., Sawada, T., Plunkett, A., & Fomalont, E. 2019, *PASP*, 131, 054505

Lapinov, A. V. 2006, *Proc. SPIE*, 6580, 658001

Law, C., Loomis, R. A., Teague, R., et al. 2021a, *ApJS*, 257, 3

Law, C., Teague, R., Loomis, R. A., et al. 2021b, *ApJS*, 257, 4

Le Gal, R., Öberg, K. I., Huang, J., et al. 2020, *ApJ*, 898, 131

Le Gal, R., Öberg, K. I., Teague, R., et al. 2021, *ApJS*, 257, 12
Linsky, J. L., Draine, B. T., Moos, H. W., et al. 2006, ApJ, 647, 1106

Lis, D. C., Bockelée-Morvan, D., Güsten, R., et al. 2019, A&A, 625, L5

Liu, Y., Dipierro, G., Ragusa, E., et al. 2019, A&A, 622, A75

Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17

Loomis, R. A., Öberg, K. I., Andrews, S. M., et al. 2018, AJ, 155, 182

Loomis, R. A., Öberg, K. I., Andrews, S. M., et al. 2020, ApJ, 893, 101

Macías, E., Espaillat, C. C., Ribas, Á, et al. 2018, ApJ, 865, 37

Maiwald, F., Lewen, F., Ahrens, V., et al. 2000, JMoSp, 202, 166

Maki, A. G., Mellau, G. C., Klee, S., Winnewisser, M., & Quapp, W. 2000, JMoSp, 202, 67

Matrà, L., MacGregor, M. A., Kalas, P., et al. 2017, ApJ, 842, 9

Mawet, D., Absil, O., Montagnier, G., et al. 2012, A&A, 544, A131

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in ASP Conf. Ser. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell (San Francisco, CA: ASP), 127

Meier, R., Owen, T. C., Jewitt, D. C., et al. 1998, Sci, 279, 1707

Milam, S. N., Savage, C., Brewster, M. A., Ziurys, L. M., & Wyckoff, S. 2005, ApJ, 634, 1126

Millar, T. J., Bennett, A., & Herbst, E. 1989, ApJ, 340, 906

Möllmann, E., Maki, A. G., Winnewisser, M., Winnewisser, B. P., & Quapp, W. 2002, JMoSp, 212, 22

Montesinos, B., Eiroa, C., Mora, A., & Merín, B. 2009, A&A, 495, 901

Mousis, O., Gautier, D., Bockelée-Morvan, D., et al. 2000, Icar, 148, 513

Müller, H. S. P., Schlöder, F., Stutzki, J., & Winnewisser, G. 2005, JMoSt, 742, 215

Müller, H. S. P., Thorwirth, S., Roth, D. A., & Winnewisser, G. 2001, A&A, 370, L49

Noble, J. A., Theule, P., Borget, F., et al. 2013, MNRAS, 428, 3262

Nyman, G., & Yu, H.-G. 2019, AIPA, 9, 095017

Öberg, K. I., & Bergin, E. A. 2021, PhR, 893, 1

Öberg, K. I., Cleeves, L. I., Bergner, J. B., et al. 2021, AJ, 161, 38

Öberg, K. I., Furuya, K., Loomis, R., et al. 2015, ApJ, 810, 112

Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, ApJS, 257, 1

Öberg, K. I., Qi, C., Fogel, J. K. J., et al. 2011, ApJ, 734, 98

Öberg, K. I., Qi, C., Wilner, D. J., & Hogerheijde, M. R. 2012, ApJ, 749, 162

Öberg, K. I., van Dishoeck, E. F., & Linnartz, H. 2009, A&A, 496, 281

O'Brien, D. P., Izidoro, A., Jacobson, S. A., Raymond, S. N., & Rubie, D. C. 2018, SSRv, 214, 47

Pagani, L., Daniel, F., & Dubernet, M. L. 2009, A&A, 494, 719

Pagani, L., Roueff, E., & Lesaffre, P. 2011, ApJL, 739, L35

Pegues, J., Öberg, K. I., Bergner, J. B., et al. 2020, ApJ, 890, 142

Pinte, C., Ménard, F., Duchêne, G., et al. 2018, A&A, 609, A47

Pinte, C., Price, D. J., Ménard, F., et al. 2020, ApJL, 890, L9

Preusser, J., & Maki, A. G. 1993, JMoSp, 162, 484

Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123

Qi, C., D'Alessio, P., Öberg, K. I., et al. 2011, ApJ, 740, 84

Qi, C., Öberg, K. I., Andrews, S. M., et al. 2015, ApJ, 813, 128

Qi, C., Öberg, K. I., Espaillat, C. C., et al. 2019, ApJ, 882, 160

Qi, C., Öberg, K. I., Wilner, D. J., et al. 2013, Sci, 341, 630

Remijan, A., Biggs, A., Cortes, P., et al. 2021, ALMA Technical Handbook ALMA Doc. 8.5, https://almascience.eso.org/documents-and-tools/cycle8/alma-technical-handbook

Roberts, H., & Millar, T. J. 2000, A&A, 361, 388

Roueff, E., Gerin, M., Lis, D. C., et al. 2013, JPCA, 117, 9959

Rubin, M., Engrand, C., Snodgrass, C., et al. 2020, SSRv, 216, 102

Salinas, V. N., Hogerheijde, M. R., Mathews, G. S., et al. 2017, A&A, 606, A125

Salinas, V. N., Hogerheijde, M. R., Murillo, N. M., et al. 2018, A&A, 616, A45

Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2016, ApJ, 823, 91

Sierra, A., Pérez, L. M., Zhang, K., et al. 2021, ApJS, 257, 14

Simon, M., Dutrey, A., & Guilloteau, S. 2000, ApJ, 545, 1034

Simon, M., Guilloteau, S., Beck, T. L., et al. 2019, ApJ, 884, 42

Teague, R. 2019, JOSS, 4, 1632

Teague, R., Bae, J., Aikawa, Y., et al. 2021, ApJS, 257, 18

Teague, R., Bae, J., Birnstiel, T., & Bergin, E. A. 2018a, ApJ, 868, 113

Teague, R., & Foreman-Mackey, D. 2018, RNAAS, 2, 173

Teague, R., & Loomis, R. 2020, ApJ, 899, 157

Teague, R., Guilloteau, S., Semenov, D., et al. 2016, A&A, 592, A49

Teague, R., Henning, T., Guilloteau, S., et al. 2018b, ApJ, 864, 133

Thorwirth, S., Müller, H. S. P., Lewen, F., et al. 2003, ApJL, 585, L163

Turner, B. E. 2001, ApJS, 136, 579

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, CSE, 13, 22

van 't Hoff, M. L. R., Walsh, C., Kama, M., Facchini, S., & van Dishoeck, E. F. 2017, A&A, 599, A101

Verhoeve, P., Zwart, E., Versluis, M., et al. 1990, RScI, 61, 1612

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Walsh, C., Nomura, H., Millar, T. J., & Aikawa, Y. 2012, ApJ, 747, 114

Walter, F., & Brinks, E. 1999, AJ, 118, 273

Walter, F., Brinks, E., de Blok, W. J. G., et al. 2008, AJ, 136, 2563

Weaver, E., Isella, A., & Boehler, Y. 2018, ApJ, 853, 113

Willacy, K. 2007, ApJ, 660, 441

Willacy, K., & Langer, W. D. 2000, ApJ, 544, 903

Willacy, K., & Woods, P. M. 2009, ApJ, 703, 479

Willacy, K., Alexander, C., Ali-Dib, M., et al. 2015, SSRv, 197, 151

Winnewisser, M., & Vogt, J. 1978, ZNatA, 33, 1323

Yen, H.-W., Koch, P. M., Liu, H. B., et al. 2016, ApJ, 832, 204

Zhang, K., Booth, A. S., Law, C. J., et al. 2021, ApJS, 257, 5

Zhu, Z., Zhang, S., Jiang, Y.-F., et al. 2019, ApJL, 877, L18