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

# Chemistry in the GG Tau A Disk: Constraints from H₂D⁺, N₂H⁺, and DCO⁺ High Angular Resolution ALMA Observations

Parashmoni Kashyap<fnref n="1" />,<fnref n="2" /> ID, Liton Majumdar<fnref n="1" />,<fnref n="2" /> ID, Anne Dutrey<fnref n="3" /> ID, Stéphane Guilloteau<fnref n="3" /> ID, Karen Willacy⁴ ID, Edwige Chapillon<fnref n="3" />,<fnref n="5" /> ID, Richard Teague<fnref n="6" /> ID, Dmitry Semenov<fnref n="7" />,<fnref n="8" /> ID, Thomas Henning<fnref n="7" /> ID, Neal Turner⁴ ID, Raghvendra Sahai⁴ ID, Ágnes Kóspál<fnref n="7" />,<fnref n="9" />,<fnref n="10" /> ID, Audrey Coutens<fnref n="11" /> ID, V. Piétu<fnref n="5" /> ID, Pierre Gratier<fnref n="3" /> ID, Maxime Ruaud<fnref n="12" />,<fnref n="13" /> ID, N. T. Phuong<fnref n="14" /> ID, E. Di Folco<fnref n="3" /> ID, Chin-Fei Lee¹⁵ ID, and Y.-W. Tang¹⁵ ID

Received 2024 February 17; revised 2024 September 23; accepted 2024 September 26; published 2024 November 27

## Abstract

Resolved molecular line observations are essential for gaining insight into the physical and chemical structure of protoplanetary disks, particularly in cold, dense regions where planets form and acquire their chemical compositions. However, tracing these regions is challenging because most molecules freeze onto grain surfaces and are not observable in the gas phase. We investigated cold molecular chemistry in the triple stellar T Tauri disk GG Tau A, which harbours a massive gas and dust ring and an outer disk, using Atacama Large Millimeter/ submillimeter Array (ALMA) Band 7 observations. We present high angular resolution maps of N₂H⁺ and DCO⁺ emission, with upper limits reported for H₂D⁺, ¹³CS, and SO₂. The radial intensity profile of N₂H⁺ shows most emission near the ring's outer edge, while DCO⁺ exhibits a double peak, one near the ring's inner edge and the other in the outer disk. With complementary observations of lower-lying transitions, we constrained the molecular surface densities and rotation temperatures. We compared the derived quantities with model predictions across different cosmic-ray ionization (CRI) rates, carbon-to-oxygen (C/O) ratios, and stellar UV fluxes. Cold molecular chemistry, affecting the N₂H⁺, DCO⁺, and H₂D⁺ abundances, is most sensitive to the CRI rate, while the stellar UV fluxes and C/O ratios have minimal impact on these three ions. Our best model requires a low CRI rate of 10⁻¹⁸ s⁻<fnref n="1" />. However, it fails to match the low temperatures derived from N₂H⁺ and DCO⁺, 12–16 K, which are much lower than the CO freezing temperature.

*Unified Astronomy Thesaurus concepts: Protoplanetary disks (1300); Planet formation (1241); Astrochemistry (75); Isotopic abundances (867)*

## 1. Introduction

A protoplanetary disk is a crucial stage in the evolution of cosmic matter, where the interplay of inherited material from a parent cloud and the in situ physics and chemistry establishes the conditions for planet formation. Protoplanetary disks are characterized by complex chemical stratification owing to radial and vertical variations in density, temperature, ionization, and dissociative radiation (Y. Aikawa et al. 2002). The uppermost layers usually allow the formation of simple atoms, ions, and photo-stable radicals, while the intermediate warm

molecular layer produces most of the gas-phase molecules traced by Atacama Large Millimeter/submillimeter Array (ALMA) observations. There is minimal penetration of stellar radiation in the innermost regions near the midplane in the vertical direction, making them very cold and poorly ionized. In this dense region, most of the molecules freeze onto grain surfaces. Hence, midplane regions are very difficult to probe. N₂H⁺, DCO⁺, and H₂D⁺ are among the very few molecules which remain in the gas phase in the midplane region (C. Ceccarelli et al. 2004; K. I. Öberg et al. 2011).

N₂H⁺ is expected to be abundant in regions where CO is depleted. This anticorrelation is observed in both dense cores (e.g., P. Caselli et al. 1999; E. A. Bergin et al. 2002) as well as in disks (C. Walsh et al. 2012). N₂H⁺ is produced through the protonation of N₂ and rapidly destroyed by CO to form HCO⁺.

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Exoplanets and Planetary Formation Group, School of Earth and Planetary Sciences, National Institute of Science Education and Research, Jatni 752050, Odisha, India; liton@niser.ac.in, dr.liton.majumdar@gmail.com</footnote>

<footnote n="2">² Homi Bhabha National Institute, Training School Complex, Anushaktinagar, Mumbai 400094, India</footnote>

<footnote n="3">³ Laboratoire d'Astrophysique de Bordeaux, Université de Bordeaux, CNRS, B18N, Allée Geoffroy Saint-Hilaire, F-33615 Pessac, France
⁴ Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr., Pasadena, CA 91109, USA</footnote>

<footnote n="5">⁵ Institut de Radioastronomie Millimétrique (IRAM), 300 rue de la Piscine, F-38406 Saint-Martin d'Hères, France</footnote>

<footnote n="6">⁶ Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, USA</footnote>

<footnote n="7">⁷ Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="8">⁸ Department of Chemistry, Ludwig Maximilian University, Butenandtstr. 5-13, D-81377 Munich, Germany</footnote>

<footnote n="9">⁹ Konkoly Observatory, HUN-REN Research Centre for Astronomy and Earth Sciences, CSFK, MTA Centre of Excellence, Konkoly-Thege Miklós út 15-17, 1121 Budapest, Hungary</footnote>

<footnote n="10">¹⁰ Institute of Physics and Astronomy, ELTE Eötvös Loránd University, Pázmány Péter sétány 1/A, 1117 Budapest, Hungary</footnote>

<footnote n="11">¹¹ Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse, UT3-PS, CNRS, CNES, 9 av. du Colonel Roche, 31028 Toulouse Cedex 4, France</footnote>

<footnote n="12">¹² NASA Ames Research Center, Moffett Field, CA 94035, USA</footnote>

<footnote n="13">¹³ Carl Sagan Center, SETI Institute, Mountain View, CA 94035, USA</footnote>

<footnote n="14">¹⁴ Department of Astrophysics, Vietnam National Space Center, Vietnam Academy of Science and Technology, 18 Hoang Quoc Viet, Cau Giay, Hanoi, Vietnam
¹⁵ Academia Sinica Institute of Astronomy and Astrophysics, PO Box 23-141, Taipei 106, Taiwan</footnote>
<table><caption>**Table 1**

Spectral Setup of the ALMA Band 7 Observations of GG Tau A</caption><thead><tr><th rowspan="2">Spectral Window</th><th rowspan="2">Central Frequency (GHz)</th><th colspan="2">Bandwidth</th><th colspan="2">Resolution</th><th rowspan="2">No. of Channels</th></tr><tr><th>Frequency (MHz)</th><th>Velocity (km s<sup>-1</sup>)</th><th>Frequency (MHz)</th><th>Velocity (km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>25</td><td>358.000000</td><td>1875.00</td><td>1570.2</td><td>1.129</td><td>0.945</td><td>3840</td></tr><tr><td>27</td><td>359.770685</td><td>117.19</td><td>97.7</td><td>0.141</td><td>0.118</td><td>1920</td></tr><tr><td>29</td><td>360.169780</td><td>117.19</td><td>97.5</td><td>0.141</td><td>0.117</td><td>1920</td></tr><tr><td>31</td><td>369.908554</td><td>117.19</td><td>95.0</td><td>0.141</td><td>0.114</td><td>1920</td></tr><tr><td>33</td><td>372.421385</td><td>117.19</td><td>94.3</td><td>0.122</td><td>0.098</td><td>1920</td></tr><tr><td>35</td><td>372.672493</td><td>58.59</td><td>47.1</td><td>0.141</td><td>0.114</td><td>960</td></tr></tbody></table>

Hence, $N_2H^+$ can persist in the gas phase at temperatures a few kelvin below the CO freeze-out temperature (K. I. Öberg et al. 2005). $N_2H^+$ has been observed in multiple disks (e.g., A. Dutrey et al. 2007; C. Qi et al. 2013; N. T. Phuong et al. 2021; D. E. Anderson et al. 2022), but is often characterized by comparatively weak emission.

An enhanced abundance of $H_2D^+$ is predicted by chemical models at low temperatures ($T < 20$ K), with the major formation pathway being protonation of HD (H. Roberts & T. J. Millar 2000). $H_2D^+$ can only be observed in high-density regions with low CO and $N_2$ content where $H_2^+$ has not already been destroyed. $H_2D^+$ has already been detected in prestellar cores (e.g., P. Caselli et al. 2008), but there is no clear detection in any of the protoplanetary disks (e.g., C. Ceccarelli et al. 2004; W. F. Thi et al. 2004; C. Qi et al. 2008; E. Chapillon et al. 2011).

$DCO^+$ serves as an effective indicator of low-temperature CO freeze-out regions within the disk and is predominantly formed through the interaction of CO with $H_2D^+$ (A. Wootten 1987). On the one hand, the formation of $H_2D^+$ mainly occurs through HD protonation and is hindered by the destruction of $H_2^+$ in the presence of CO. On the other hand, the presence of CO is essential for $DCO^+$ formation. Consequently, $DCO^+$ emission is expected to peak in a narrow region around the CO snow line where both the parent molecules are in the gas phase in small amounts (Y. Aikawa et al. 2002; K. Willacy 2007; G. S. Mathews et al. 2013).

In summary, $N_2H^+$, $DCO^+$, and $H_2D^+$ are excellent tracers of low-temperature, high-density regions of the disk, specifically the midplane region in outer disks. Simultaneous observations of these molecules allow us to locate the CO snow line. Snow lines are thought to play an important role in planet formation because the ice coating on grains increases the solid mass, potentially promotes the coagulation of grains into larger particles, and/or the larger bodies formed from these grains require a higher collision velocity for destruction. These effects are particularly prominent just outside the snow line (S. Okuzumi et al. 2012). Additionally, such observations provide insights into ionization fractions within the relevant disk areas, offering constraints that can impact magnetorotational instability. In instances of excessively low ionization fractions, “dead zones” (C. F. Gammie 1996) form, which are conducive to efficient grain growth and subsequent planet formation (e.g., K. I. Öberg et al. 2011).

This work reports spatially resolved emission from $N_2H^+$ (4–3) and $DCO^+$ (5–4) and upper limits for $O-H_2D^+$ (1(1, 0)–1(1, 1)), $^{13}$CS (8–7), and $SO_2$ 19(4, 16)–19(3, 17) in the disk around a triple stellar T Tauri system, GG Tau A, utilizing the ALMA interferometer. This represents the most sensitive observation of these molecules in GG Tau A to date. The system (age 1–5 Myr) is located at a distance of 150 pc (using

the Gaia parallax measurements of the GG Tau Ba star of this hierarchical quintuple system, because those of GG Tau A are contaminated by its triple nature; Gaia Collaboration et al. 2016, 2018) from us in the Taurus–Auriga star-forming region. This circumtertiary disk is characterized by a dense ring of gas and dust between 193 and 285 au (from here on, we simply use the term “ring” to denote this region) where 90% of the continuum emission arises and an outer disk extended up to 800 au in $^{13}$CO observations (S. Guilloteau et al. 1999; A. Dutrey et al. 2014). The majority of the total disk mass (0.15 $M_⊙$) is concentrated in the ring region (0.13 $M_⊙$; S. Guilloteau et al. 1999; S. M. Andrews et al. 2014). A. Dutrey et al. (2014) reported the presence of a hot spot located at the outer edge of the dense ring where a partial CO gap is observed (Y.-W. Tang et al. 2016), suggesting the presence of an embedded planet. This hypothesis is reinforced by the detection of a CO spiral emanating from the hot spot (N. T. Phuong et al. 2020a).

This disk is very cold, with a midplane temperature of ~14 K (A. Dutrey et al. 2014) and an atmospheric temperature of ~30 K (S. Guilloteau et al. 1999; N. T. Phuong et al. 2020b), at a radial distance of 214 au from the center. The first ever detection of $H_2S$ in the outer disk of GG Tau A by N. T. Phuong et al. (2018) further confirms how dense and massive this system is. Given its substantial mass, high density (which should ease the detection of weak molecular lines), and low temperature, GG Tau A emerges as an ideal candidate for unraveling the intricacies of cold molecular chemistry within a protoplanetary disk.

We present our observational details in Section 2 of this paper. We explain the schemes we undertake to reduce the data in Section 2.1 and the observational results in Section 2.2, including integrated intensity maps, radial profiles, and Keplerian deprojected spectra. We also constrain the surface density of the detected molecules and the upper limits of surface density for the undetected molecules, as explained in Section 3. We perform a chemical analysis to understand our observations better using our astrochemical model, detailed in Section 4. We discuss our results in Section 5 and summarize the key understandings of our study in Section 6.

## 2. Observations

The GG Tau A system (International Celestial Reference System 04:32:30.3460, 17:31:40.642) was observed as part of the ALMA project #2021.1.00342.S (PI: Liton Majumdar) targeting to investigate the presence of Band 7 transitions, $O-H_2D^+$ (1(1, 0)–1(1, 1)), $N_2H^+$ (4–3), $DCO^+$ (5–4), $^{13}$CS (8–7), and $SO_2$ 19(4, 16)–19(3, 17). The observations were made over two executions using 47 antennas with baselines ranging from 14 to 456 m for a total on-source integration time of 2.7 hr. The spectral setup details can be found in Table 1.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Integrated intensity maps of the targeted lines N<sub>2</sub>H<sup>+</sup> (4–3), DCO<sup>+</sup> (5–4), and o–H<sub>2</sub>D<sup>+</sup> (1(1, 0)–1(1, 1)). The overplotted elliptical ring contours in white are at 193 au (1.′′28) and 285 au (1.′′9) to indicate the continuum ring region. Both N<sub>2</sub>H<sup>+</sup> and DCO<sup>+</sup> emissions are extended beyond the ring region. As far as H<sub>2</sub>D<sup>+</sup> is concerned, there is no significant emission from the ring region. The figure at the lower right corner is the continuum image at 358 GHz.</figcaption></figure>

indicate that the N<sub>2</sub>H<sup>+</sup> (4–3) emission exhibits a ring-like radial distribution, peaking near the outer edge of the continuum ring. In contrast, the radial distribution of the DCO<sup>+</sup> (5–4) peak brightness temperature (Figure B1) shows a double-ringed structure. The inner DCO<sup>+</sup> (5–4) emission arises from the continuum ring, while an outer, brighter DCO<sup>+</sup> (5–4) emission originates from the outer disk region. This suggests radial variations in the deuteration process across the disk. However, this double-ringed structure is less prominent in the velocity-integrated radial profile (right panel in Figure 2), possibly due to the lower angular resolution. Other targeted species, including H<sub>2</sub>D<sup>+</sup>, <sup>13</sup>CS, and SO<sub>2</sub>, are not detected (Figure A5).

**3. Derivation of the Disk's Physical Parameters**

**3.1. DISKFIT Modeling**

To derive disk properties, we perform a least squares fit to the observed visibilities to adjust an empirical parametric disk model to the data, using the DISKFIT tool from V. Piétu et al. (2007).

This approach allows an accurate derivation of the geometrical parameters. Furthermore, when several transitions of the same molecules are available, it allows a quantitative estimate of the molecular surface densities and excitation temperatures. It also allows for consistent estimates of the errors on the disk parameters (within the framework of the adopted disk model). This can be done either from the covariance matrix (when parameter coupling is limited) or through a more elaborate Markov Chain Monte Carlo method.

The disk model is similar to that used in N. T. Phuong et al. (2020b). It uses a flared disk geometry with a Gaussian vertical profile and a radial power law for the scale height. Molecules are assumed to follow the same vertical distribution for simplicity. The model assumes level populations are governed by the Boltzmann law, with the temperature being a power law as a function of radius. For a simple comparison with the observed data, we assume that the intrinsic line width remains constant with radius. In fact, fitting line width as a power law instead suggests a small but not significant decrease of intrinsic
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Azimuthally averaged velocity-integrated radial profiles of the detected lines generated from the integrated intensity maps shown in Figure 1. The shaded region in the above figure corresponds to $1\sigma$. The dashed vertical lines are to designate the location of the ring around GG Tau A, spanning from 193 to 285 au. The gray horizontal line at the top right corner represents the beam major axis.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Keplerian deprojected and stacked spectra of the detected lines $N_2H^+$ (4–3) and $DCO^+$ (5–4). The shaded region is the region over which we have done a Gaussian fitting (orange curve) and calculated the integrated flux densities. The stacked spectra are peaking at the systemic velocity, $v_{LSRK}$ of 6.4 km s$^{-1}$, confirming the detections. The location of the systemic velocity is shown by the vertical dashed lines. Integrated flux densities and matched filter responses are tabulated in Table 2.</figcaption></figure>

<table><caption>**Table 2**

Properties of the Lines Observed</caption><thead><tr><th>Species</th><th>Transition</th><th>Frequency (GHz)</th><th>E<sub>u</sub> (K)</th><th>log<sub>10</sub>(A<sub>ul</sub>) (s<sup>-1</sup>)</th><th>g<sub>u</sub></th><th>Filter Response (σ)</th><th>Integrated Intensity (mJy km s<sup>-1</sup>)</th><th>Line Width (km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>N<sub>2</sub>H<sup>+</sup></td><td>4–3</td><td>372.6724808</td><td>44.71</td><td>-2.50934</td><td>81</td><td>43</td><td>1570 ± 40</td><td>0.55 ± 0.01</td></tr><tr><td>DCO<sup>+</sup></td><td>5–4</td><td>360.1697783</td><td>51.86</td><td>-2.42480</td><td>11</td><td>20</td><td>290 ± 20</td><td>0.46 ± 0.01</td></tr><tr><td>H<sub>2</sub>D<sup>+</sup></td><td>1(1, 0)–1(1, 1)</td><td>372.4213558</td><td>104.20</td><td>-3.96567</td><td>9</td><td>...</td><td>&lt;58</td><td>...</td></tr><tr><td><sup>13</sup>CS</td><td>8–7</td><td>369.9085505</td><td>79.89</td><td>-2.97351</td><td>34</td><td>...</td><td>&lt;70</td><td>...</td></tr><tr><td>SO<sub>2</sub></td><td>19(4, 16)–19(3, 17)</td><td>359.7706846</td><td>214.26</td><td>-3.41473</td><td>39</td><td>...</td><td>&lt;35</td><td>...</td></tr></tbody></table>

**Note.** Here, $E_u$ is the upper energy level of the transition, $A_{ul}$ is the Einstein coefficient for spontaneous transition, and $g_u$ is the upper state degeneracy. The error bars on the integrated intensity for the detected molecules are $1\sigma$. We have calculated a $3\sigma$ upper limit on flux densities for the undetected molecules for an integration area between radii of 150 and 450 au and a line width of 0.55 km s$^{-1}$ (Appendix C). The line width comes from the FWHM of the Gaussian fitting in Figure 3 (Section 2.2).

line width with radius that does not affect the other disk parameters. The emerging line intensity is derived using simple ray tracing, considering the geometry and dynamics (position, orientation, inclination, velocity, and Keplerian rotation pattern) and level populations computed at every sampled

point along the line of sight. The hyperfine structure of the observed lines is also considered at this stage. This approach implicitly accounts for line optical depth. Visibilities are then computed on the same (u, v) points as the data from the 3D data cube. A least squares minimization is performed on the
<table><caption>**Table 3**

Geometric and Dynamic Parameters Determined through DISKFIT Modeling</caption><thead><tr><th>Quantity</th><th>N<sub>2</sub>H<sup>+</sup> (4–3)</th><th>DCO<sup>+</sup> (5–4)</th><th>Adopted</th></tr></thead><tbody><tr><td>X<sub>0</sub> (arcsec)</td><td>-0.007 ± 0.011</td><td>0.019 ± 0.016</td><td>0</td></tr><tr><td>Y<sub>0</sub> (arcsec)</td><td>-0.001 ± 0.007</td><td>0.004 ± 0.011</td><td>0</td></tr><tr><td>PA (deg)</td><td>8.4 ± 0.5</td><td>8.1 ± 0.5</td><td>8</td></tr><tr><td>i (deg)</td><td>36.5 ± 0.4</td><td>37.0 ± 0.4</td><td>37</td></tr><tr><td>V<sub>sys</sub> (km s<sup>-1</sup>)</td><td>6.45 ± 0.03</td><td>6.43 ± 0.03</td><td>6.43</td></tr><tr><td>V<sub>100</sub> (km s<sup>-1</sup>)</td><td>3.61 ± 0.03</td><td>3.54 ± 0.03</td><td>3.55</td></tr></tbody></table>

**Note.** Best-fit geometric parameters from the observed visibilities. Offsets ($X_0$, $Y_0$) are from the ring center. $V_{100}$ is the Keplerian velocity at 100 au. A PA of 8° in DISKFIT is equivalent to 278° in GoFish. In DISKFIT, the PA is measured with respect to the disk minor axis, while in GoFish, this angle is measured with respect to the disk major axis. The “Adopted” column indicates the fixed parameter values used for the multitransition model fit in Table 4.

different visibilities using a modified Levenberg-Marquardt algorithm. Guided by the radial distribution of velocity-integrated brightness (Figure 2) and apparent brightness temperature (Figure B1), we represented the radial surface density profiles for each molecule by the sum of two Gaussians, truncated to inner and outer radii of 150 and 550 au, respectively. The rotation temperature is assumed to be a simple power law, defined by its value at 250 au, $T_{250}$, and exponent $q$. The model thus has 15 possible parameters, six due to geometry and nine to represent the molecular distribution and excitation conditions (see Tables 3–4 for their designation).

Error bars are derived from the covariance matrix. The six geometric parameters, the line width $\delta V$, and temperature parameters exhibit very little coupling among them and with the other ones. However, the Gaussian distribution parameters are strongly coupled, particularly the pair ($N_i$, $W_i$) for distribution $i$, so their error bars should be treated with caution.

The parameter fitting process involves a two-stage approach. In the first stage, the geometric parameters are individually fitted for every observed transition. The geometric parameters derived from $N_2H^+$ and $DCO^+$ are given in Table 3, using the distance of 150 pc determined from Gaia. These transitions are relatively optically thin and, thus, more reliable tracers of the disk inclination than those previously used. Derived inclinations and orientation agree (but offer better precision) with previous determinations from CO isotopologues (A. Dutrey et al. 2014; N. T. Phuong et al. 2018).

The second stage entails averaging the corresponding values of the six geometric parameters across both spectral lines, fixing these averaged values, and proceeding to fit the remaining nine parameters.

To derive the surface densities and rotation temperatures, we used complementary data, which are made of the interferometric visibilities obtained with NOEMA for $N_2H^+$ (1–0) and $DCO^+$ (1–0) from N. T. Phuong et al. (2021), and $DCO^+$ (3–2) from N. T. Phuong et al. (2018). Geometric parameters from the “Adopted” column of Table 3 were used at this stage. Results are presented in Table 4. Figure 4 shows the constrained radial distributions of the surface densities. Agreement with the observations can be seen in Figure B1, where the best-fit radial profile (reconstructed by imaging the model visibilities and applying the Kepler deprojection) is compared with the data.

While this process is appropriate for $DCO^+$, the $N_2H^+$ distribution can almost equally well be represented by a single

<table><caption>**Table 4**

Best-fit Parameters for the Molecular Distributions and Excitation Conditions</caption><thead><tr><th>Quantity</th><th>N<sub>2</sub>H<sup>+</sup></th><th>DCO<sup>+</sup></th><th>Units</th></tr></thead><tbody><tr><td>δV</td><td>0.149 ± 0.005</td><td>0.150 ± 0.005</td><td>km s<sup>-1</sup></td></tr><tr><td>T<sub>250</sub></td><td>12.2 ± 0.7</td><td>16.0 ± 1.2</td><td>K</td></tr><tr><td>q</td><td>-0.32 ± 0.19</td><td>0.39 ± 0.15</td><td>...</td></tr><tr><td>N<sub>1</sub></td><td>22 ± 8</td><td>2.8 ± 0.6</td><td>10<sup>11</sup> cm<sup>-2</sup></td></tr><tr><td>R<sub>1</sub></td><td>252 ± 9</td><td>210 ± 50</td><td>au</td></tr><tr><td>W<sub>1</sub></td><td>95 ± 115</td><td>90 ± 360</td><td>au</td></tr><tr><td>N<sub>2</sub></td><td>12 ± 2</td><td>5.3 ± 1.0</td><td>10<sup>11</sup> cm<sup>-2</sup></td></tr><tr><td>R<sub>2</sub></td><td>336 ± 20</td><td>375 ± 10</td><td>au</td></tr><tr><td>W<sub>2</sub></td><td>145 ± 15</td><td>90 ± 15</td><td>au</td></tr></tbody></table>

**Note.** Parameters derived from the best fit of a double Gaussian radial distribution model. $N_i$ is the peak surface density of Gaussian distribution $i$, $R_i$ is its peak position, and $W_i$ is its width (FWHM). $\delta V$ is the 1/e line width. $T_{250}$ is the temperature at 250 au and $q$ is the exponent of the temperature power law (a positive exponent indicates a quantity decreasing with radius).

Gaussian or a truncated power law. This, however, does not affect the temperature derivation. The temperature that best represents $N_2H^+$ is almost constant, slightly increasing with radius, and quite low: 12 K. On the contrary, $DCO^+$ is better represented by higher temperatures, decreasing with radius from ~18 K at 200 au to ~12 K at 500 au.

As expected, the FWHM line width derived from the uv plane analysis with DISKFIT (0.25 km s⁻¹ using the 1/e values in Table 4) is smaller than those in the recentered spectra from GoFish (0.55 km s⁻¹, Figure 3) and KEPLER (0.37 km s⁻¹, Figure B1) because of the remaining Keplerian shear within the synthesized beams used in these image plane reconstructions.

For $H_2D^+$, we simply assumed the same geometrical parameters (“Adopted” column in Table 3), a constant temperature of 15 K, and a uniform surface density between 150 and 450 au. The best-fit surface density of o–$H_2D^+$ from DISKFIT is then $(0.7 \pm 0.7) \times 10^{10} \text{ cm}^{-2}$. This gives the $3\sigma$ upper limit on o–$H_2D^+$ to be $2.1 \times 10^{10} \text{ cm}^{-2}$.

## 3.2. Results

Figure 4 displays the radial distributions derived from the DISKFIT modeling. The results from a classical local thermodynamic equilibrium approach on the disk-averaged quantities (presented in Appendix C) are also indicated. The discrepancy between the radial profile and disk-averaged value for $DCO^+$ at large radii can be attributed to the assumption of constant $T_{rot}$ = 16 K while the DISKFIT model shows it drops to 12 K at 500 au, leading to a stronger Boltzmann distribution correction factor.

We did not detect the lines o–$H_2D^+$ (1(1, 0)–1(1, 1)), $^{13}$CS (8–7), and $SO_2$ 19(4, 16)–19(3, 17). We calculated upper limits on the disk-averaged surface densities for these molecules using Equation (C3), while assuming optically thin emission. The calculations considered $3\sigma$ values integrated over the line width of $N_2H^+$ (coming from the Gaussian fit in Figure 3). Here, the noise $\sigma = f \times \delta S \times \sqrt{\delta v \times \Delta v}$, where $f$ is the spectral correlation factor of 1.6, $\delta S$ is the rms noise in the flux, $\delta v$ is the channel width, and $\Delta v$ is the line width over which the noise is calculated. The excitation temperature was fixed at 15 K for this purpose. The integration area $\Omega$ is considered for a radial region ranging from 150 to 450 au. The constrained values are tabulated in Table 5.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Radial distribution of the surface density derived from the DISKFIT modeling. The solid purple lines represent the total surface density while the dashed lines show the contribution from the individual Gaussians fitted in DISKFIT. The gray horizontal lines indicate the derived disk-averaged values using the simple LTE approach.</figcaption></figure>

Our simple LTE approach reports an upper limit on the total (summation of both the ortho and para species) H₂D⁺ surface density of 7.48 × 10¹² cm⁻², while our DISKFIT analysis quotes an upper limit on the ortho–H₂D⁺ surface density of 2.1 × 10¹⁰ cm⁻². Converting this to a total H₂D⁺ surface density requires an assumption on the ortho-to-para ratio. Under the hypothesis of thermalization, this ratio can be very small since the ortho–H₂D⁺ ground state is 86 K above the para ground state, leading to a huge Boltzmann correction factor (from 70 to 300 for temperatures ranging from 20 down to 15 K, respectively). At 15 K, the 3σ upper limit on the total H₂D⁺ surface density, derived from the DISKFIT ortho–H₂D⁺ value, is 6.3 × 10¹² cm⁻².

Note that the low-level detection of ¹³CS (2–1) with NOEMA, a low energy transition reported by N. T. Phuong et al. (2021), implies a surface density of 3 × 10¹¹ cm⁻² for T<sub>rot</sub> = 15 K, below our upper limit.

## 4. Astrochemical Modeling

### 4.1. Disk Physical Model

Our physical model is constructed for a radial range of 190–290 au, which primarily includes the ring region around the central stellar system GG Tau A. Only 10%–20% of the total disk mass is contained in the outer disk region, which, therefore, exhibits a lower density than in the ring. Moreover, the outer disk, under the shadow of the dense ring, is very cold (A. Dutrey et al. 2014; Y.-W. Tang et al. 2016; R. Brauer et al. 2019) with poor constraints on the physical properties. We chose not to focus our present modeling effort on this outer region, which has a more complex structure (Y.-W. Tang et al. 2016; N. T. Phuong et al. 2020b).

In the ring, the radial distribution of the midplane (T<sub>mid</sub>) and the atmospheric temperature (T<sub>atm</sub>) are taken as power laws, constrained in previous studies by A. Dutrey et al. (2014) and S. Guilloteau et al. (1999), respectively. They are of the form

$$T_{\text{mid}}(r) = T_{\text{mid},R_{\text{ref}}} \left( \frac{r}{R_{\text{ref}}} \right)^{-q}, \quad (1)$$

<table><caption>**Table 5**
Disk-averaged Surface Densities of the Detected Molecules and Upper Limits on the Undetected Molecules Under the Local Thermodynamic Equilibrium Assumption</caption><thead><tr><th>Species</th><th>Transition</th><th>Surface Density, N<sub>T</sub> (cm<sup>-2</sup>)</th><th>Excitation Temperature, T<sub>ex</sub> (K)</th><th>Optical Depth, τ</th></tr></thead><tbody><tr><td>N<sub>2</sub>H<sup>+</sup></td><td>4–3</td><td>(1.19 ± 0.04) × 10<sup>12</sup></td><td>12 (fixed)</td><td>0.42</td></tr><tr><td>DCO<sup>+</sup></td><td>5–4</td><td>(1.27 ± 0.08) × 10<sup>11</sup></td><td>16 (fixed)</td><td>0.04</td></tr><tr><td>H<sub>2</sub>D<sup>+</sup></td><td>1(1, 0)–1 (1, 1)</td><td>&lt;7.48 × 10<sup>12</sup></td><td>15 (fixed)</td><td>...</td></tr><tr><td><sup>13</sup>CS</td><td>8–7</td><td>&lt;1.12 × 10<sup>12</sup></td><td>15 (fixed)</td><td>...</td></tr><tr><td>SO<sub>2</sub></td><td>19(4, 16)–19(3, 17)</td><td>&lt;2.46 × 10<sup>16</sup></td><td>15 (fixed)</td><td>...</td></tr></tbody></table>

**Notes.** (a) Surface density errors are coming from 16th and 84th percentiles of the posterior distribution. (b) The upper limits on the disk-averaged surface densities of the undetected molecules (H₂D⁺, ¹³CS, and SO₂) are constrained using the 3σ upper limit on the line flux as explained in Section C. (c) In this table, we have reported the upper limit on the total (sum of the ortho and para spin states) H₂D⁺ surface density. The 3σ upper limit for o–H₂D⁺ from the DISKFIT analysis is 2.1 × 10¹⁰ cm⁻².

$$T_{\text{atm}}(r) = T_{\text{atm},R_{\text{ref}}} \left( \frac{r}{R_{\text{ref}}} \right)^{-q} \quad (2)$$

The values of the parameters can be found in Table 6. The radial dependency of the scale height, H<sub>r</sub>, at the disk midplane is calculated assuming hydrostatic equilibrium using the midplane temperature (Equation (1)), the mass of the stellar system, M<sub>*</sub> = 1.37 M<sub>☉</sub> (rescaled from S. Guilloteau et al. 1999 considering a distance of 150 pc), the mean molecular weight of the gas per H nuclei μ = 2.37, and the atomic mass unit, m<sub>H</sub>, as follows

$$H_r = \sqrt{\frac{k_B T_{\text{mid}} r^3}{\mu m_{\text{H}} G M_*}} \quad (3)$$

Here k<sub>B</sub> and G are the Boltzmann and gravitational constants, respectively. We constructed the 1D vertical structure at each
<table><caption>**Table 6**
Parameter Prescription of the Best-fit Model Considered</caption><thead><tr><th>Parameter Description</th><th>Values</th><th>Units</th></tr></thead><tbody><tr><td colspan="3"><strong>Fixed Parameters</strong></td></tr><tr><td>Stellar mass: M<sub>*</sub></td><td>1.37</td><td>M<sub>☉</sub></td></tr><tr><td>Ring mass: M<sub>g</sub></td><td>0.13<sup>a</sup></td><td>M<sub>☉</sub></td></tr><tr><td>Ring inner radius: R<sub>in</sub></td><td>193</td><td>au</td></tr><tr><td>Ring outer radius: R<sub>out</sub></td><td>285</td><td>au</td></tr><tr><td>Reference radius: R<sub>ref</sub></td><td>214</td><td>au</td></tr><tr><td>Midplane temperature at R<sub>ref</sub>: T<sub>mid, R<sub>ref</sub></sub></td><td>14</td><td>K</td></tr><tr><td>Atmospheric temperature at R<sub>ref</sub>: T<sub>atm, R<sub>ref</sub></sub></td><td>30</td><td>K</td></tr><tr><td>Surface density at R<sub>ref</sub>: Σ<sub>g,ref</sub></td><td>9.7</td><td>g cm<sup>-2</sup></td></tr><tr><td>Temperature power-law index: q</td><td>1.1</td><td></td></tr><tr><td>Temperature vertical gradient index: σ</td><td>0.5</td><td></td></tr><tr><td>Surface density power-law index: γ</td><td>1.4</td><td></td></tr><tr><td>Grain size: r<sub>d</sub></td><td>0.1</td><td>μm</td></tr><tr><td colspan="3"><strong>Adjusted Parameters</strong></td></tr><tr><td>UV flux: f<sub>UV,R<sub>ref</sub></sub></td><td>375</td><td>Draine's unit (χ<sub>0</sub>)</td></tr><tr><td>Cosmic-ray ionization (CRI) rate, ζ<sub>CR</sub></td><td>10<sup>-18</sup></td><td>s<sup>-1</sup></td></tr><tr><td>Dust settling</td><td>Yes</td><td></td></tr><tr><td>C/O ratio</td><td>1.0</td><td></td></tr><tr><td>Age of the parent molecular cloud</td><td>10<sup>6</sup></td><td>yr</td></tr><tr><td>Age of the disk</td><td>10<sup>6</sup></td><td>yr</td></tr></tbody></table>

**Note.** All the fixed values, except $\sigma$ and $\Sigma_{\text{ref}}$, are taken from S. Guilloteau et al. (1999), A. Dutrey et al. (2014), and $^a$S. M. Andrews et al. (2014) and rescaled considering a distance of 150 pc. $\sigma$ is expected to be small as the inner edge of the ring casts a shadow on the rest of the structure (Y.-W. Tang et al. 2016). $\Sigma_{\text{ref}}$ is calculated using Equation (5) assuming 0.13 $M_{\odot}$ of total disk mass is contained in the ring region.

radial point up to a $4 H_r$ height. The vertical temperature structure is adopted from E. Dartois et al. (2003) and J. P. Williams & W. M. J. Best (2014), which is given by

$$T(z) = T_{\text{mid}} + (T_{\text{atm}} - T_{\text{mid}}) \left[ \sin\left(\frac{\pi z}{2z_q}\right) \right]^{2\sigma} \quad (4)$$

$T_{\text{mid}}$ and $T_{\text{atm}}$ are defined at the midplane ($z=0$) and the upper end ($z=z_q$) of the atmosphere assuming Equations (1) and (2), respectively. $z_q$ is set at $4 H_r$ in our case. The parameter $\sigma$ denotes the stiffness of the vertical temperature distribution. As pointed out by Y.-W. Tang et al. (2016), the inner edge of the GG Tau ring obstructs a substantial amount of starlight, minimizing the vertical temperature variation within. Therefore, we have chosen a value of 0.5 for $\sigma$. We have kept the dust temperature equal to the gas temperature, as dust and gas are relatively well coupled for such densities ($>10^8 \text{ cm}^{-3}$).

The gas surface density, $\Sigma_g(r)$ [units of $g \text{ cm}^{-2}$] is parameterized as a power-law profile suggested by D. Lynden-Bell & J. E. Pringle (1974), given by

$$\Sigma_g(r) = \Sigma_{g, R_{\text{ref}}} \left( \frac{r}{R_{\text{ref}}} \right)^{-\gamma} \quad (5)$$

The power-law index $\gamma$ is taken to be 1.4 following A. Dutrey et al. (2014). The gas surface density at the reference radius $R_{\text{ref}}$, $\Sigma_{g, R_{\text{ref}}}$ is calculated considering the mass contained in the circumtertiary ring, $M_g = 0.13 M_{\odot}$ (S. Guilloteau et al. 1999;

S. M. Andrews et al. 2014). The calculations are as follows

$$M_g = 2\pi \int_{R_{\text{in}}}^{R_{\text{out}}} \Sigma_g(r) r dr. \quad (6)$$

Substituting Equation (5) in Equation (6), for $R_{\text{ref}} = 214$ au, we get

$$\Sigma_{g, R_{\text{ref}}} = \frac{(2 - \gamma) M_g}{2\pi R_{\text{ref}}^{\gamma}} [R_{\text{out}}^{2-\gamma} - R_{\text{in}}^{2-\gamma}]^{-1} = 9.7 \text{ g cm}^{-2}. \quad (7)$$

We then calculate the $H_2$ number density [units of $\text{cm}^{-3}$] at the midplane with

$$n_{H_2, \text{midplane}} = \frac{\Sigma_g(r)}{\mu \text{m}_H H_r \sqrt{2\pi}}. \quad (8)$$

Using the midplane density as the lower limit, we integrate the equation of hydrostatic equilibrium to obtain the vertical density structure (L. Reboussin et al. 2015)

$$\frac{\partial \ln n_{H_2}(z)}{\partial z} = - \left[ \left( \frac{GM_*z}{r^3} \right) \left( \frac{\mu m_H}{k_B T} \right) + \frac{\partial \ln T}{\partial z} \right]. \quad (9)$$

The local visual extinction is determined from the hydrostatic density structure, assuming a conversion factor of $(A_v/N_H)_0 = 6.25 \times 10^{-22}$ (R. Wagenblast & T. W. Hartquist 1989) where $N_H = 2N_{H_2}$ is the vertical hydrogen surface density. We assume $N_H$, and consequently $A_v$, is zero above 4 $H_r$. To account for the impact of the grain size ($r_d$) and dust-to-gas mass ratio ($\epsilon$), we scale the conversion factor as follows (V. Wakelam et al. 2019)

$$\frac{A_v}{N_H} = \left(\frac{A_v}{N_H}\right)_0 \frac{\epsilon}{10^{-2}} \frac{10^{-5}}{r_d(\text{cm})}. \quad (10)$$

The stellar UV flux factor $f_{\text{UV}}$ follows the inverse-squared law of distance and acts on the top of the disk atmosphere, only half of which is scattered toward the midplane, the rest being absorbed or scattered upwards (V. Wakelam et al. 2016)

$$f_{\text{UV}}(r) = \frac{1}{2} \frac{f_{\text{UV},R_{\text{ref}}}}{\left(\frac{r}{R_{\text{ref}}}\right)^2 + \left(\frac{4H_r}{R_{\text{ref}}}\right)^2}. \quad (11)$$

$f_{\text{UV},R_{\text{ref}}}$ is defined at the surface of the disk in units of interstellar Draine field spectrum, $\chi_0$ (B. T. Draine 1978) and reported in Table 6. $H_r$ is calculated using Equation (3). The vertical distribution of UV flux is obtained by multiplying $f_{\text{UV}}$ with $e^{-A_v/1.086}$, where $A_v$ is the visual extinction distribution (F. Du & E. A. Bergin 2014).

In our models, we have not implemented grain growth; instead, we considered a uniform grain size of 0.1 $\mu$m (N. T. Phuong et al. 2018) settling below 1 $H_r$. To incorporate dust settling, we adjust $\epsilon$ to $10^{-3}$ above 1 $H_r$, while maintaining the standard value of $\epsilon$ of $10^{-2}$ below 1 $H_r$. Figure 5 illustrates the corresponding structure of the best physical model.

## 4.2. Chemical Model and Network

The gas-grain astrochemical model DNAUTILUS, introduced in L. Majumdar et al. (2017), is used to compute the 1+1D chemical composition in the ring of GG Tau A by incorporating all the updates from DNAUTILUS 2.0 outlined in K. Taniguchi et al. (2024). DNAUTILUS 2.0 has the capability to investigate deuterium fractionation in both two-phase (gas and grain surface) and three-phase (gas, grain surface, and grain bulk) modes. DNAUTILUS 2.0 includes 1606 gas species, 1472 grain-
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Physical structure of the ring region in our model derived using the parameters tabulated in Table 6. The left panel indicates the temperature distribution, the middle one is the gas density distribution in log scale, and the one on the right is the visual extinction in log scale. The description of the model is explained in Section 4.</figcaption></figure>

surface species, and 1472 grain-mantle species, connected by 83,715 gas-phase reactions, 10,967 reactions on grain surfaces, and 9431 reactions in the grain mantles.

Although including spin-state chemistry of hydrogenated species like H₂, H⁺, H₃⁺, and their isotopologues would have allowed us to directly investigate o–H₂D⁺ chemistry, we chose not to do so. The primary reason is that the newly proposed state-to-species and species-to-state reaction rates for the H₃⁺ + H₂ systems by O. Sipilä et al. (2017) vary depending on density and temperature regimes. Applying them to a protoplanetary disk would require rigorous benchmarking and testing across various star-forming environments, and such extensive validation is beyond the scope of this paper.

### 4.3. Grid of Models

To find the model that best represents our observationally constrained radial surface density distributions of N₂H⁺ and DCO⁺, and the upper limits on H₂D⁺, SO₂, and ¹³CS (including the constraint from the *J* = 2–1 transition by N. T. Phuong et al. 2021), we explored a grid of models utilizing the least squares method (refer to Appendix D). These models span a range of initial carbon-to-oxygen (C/O) ratios, CRI rates, and UV flux values (see Table 8 in Appendix D). The model with the least disagreement with observations (χ²<sub>red</sub> = 2.78) is an inheritance disk model with initial abundances characterized from a starless dense molecular cloud, illuminated by stellar UV radiation of 375 χ₀ at a reference radius of 214 au. First, we allowed the molecular cloud to evolve chemically for 10⁶ yr from initial atomic abundances (Table 7) with typical physical conditions such as a gas and dust temperature of 10 K, a total gas density of 2 × 10⁴ cm⁻³, a visual extinction of 15 mag, and a CRI rate of 1.3 × 10⁻¹⁷ s⁻¹. We then computed the temporal evolution of chemistry for the inherited material from the cloud in the ring region around GG Tau A for 10⁶ yr. The assumption that protoplanetary disks form from a cold, starless, dense molecular cloud with the mentioned physical conditions is based on detailed discussions from V. Wakelam et al. (2019), without taking into account the chemical evolution between the cold core phase and the protoplanetary disk itself (M. N. Drozdovskaya et al. 2016).

<table><caption>**Table 7**
Initial Abundances Used in Our Cloud Model</caption><thead><tr><th>Element</th><th>Abundance Relative to H</th><th>References</th></tr></thead><tbody><tr><td>H<sub>2</sub></td><td>0.5</td><td></td></tr><tr><td>He</td><td>9.00 × 10<sup>-2</sup></td><td>(1)</td></tr><tr><td>N</td><td>6.20 × 10<sup>-5</sup></td><td>(2)</td></tr><tr><td>O</td><td>(1.40–3.3) × 10<sup>-4a</sup></td><td>(3)</td></tr><tr><td>C<sup>+</sup></td><td>1.70 × 10<sup>-4</sup></td><td>(2)</td></tr><tr><td>S<sup>+</sup></td><td>8.00 × 10<sup>-8</sup></td><td>(4)</td></tr><tr><td>Si<sup>+</sup></td><td>8.00 × 10<sup>-9</sup></td><td>(4)</td></tr><tr><td>Fe<sup>+</sup></td><td>3.00 × 10<sup>-9</sup></td><td>(4)</td></tr><tr><td>Na<sup>+</sup></td><td>2.00 × 10<sup>-9</sup></td><td>(4)</td></tr><tr><td>Mg<sup>+</sup></td><td>7.00 × 10<sup>-9</sup></td><td>(4)</td></tr><tr><td>P<sup>+</sup></td><td>2.00 × 10<sup>-10</sup></td><td>(4)</td></tr><tr><td>Cl<sup>+</sup></td><td>1.00 × 10<sup>-9</sup></td><td>(4)</td></tr><tr><td>F</td><td>6.68 × 10<sup>-9</sup></td><td>(5)</td></tr><tr><td>HD</td><td>1.60 × 10<sup>-5</sup></td><td>(6)</td></tr></tbody></table>

**Note.** (a) O initial abundance is adjusted for a C/O ratio ranging from 0.5 to 1.2.

**References.** (1) V. Wakelam & E. Herbst (2008); (2) E. B. Jenkins (2009); (3) L. Reboussin et al. (2015); (4) low metal abundances from T. E. Graedel et al. (1982); (5) depleted value from D. A. Neufeld et al. (2005); and (6) L. Majumdar et al. (2017).

It is crucial to highlight that in order to reproduce the observed surface densities in our disk simulation, a CRI rate (ζ<sub>CR</sub>) of about 10⁻¹⁸ s⁻¹ and an initial C/O ratio of 1.0 are necessary. Table 6 presents the parameters describing the best-fit model. Figure 6 showcases the modeled number density distributions within the ring region.

## 5. Discussion

### 5.1. Comparison with Previous Observations

o–H₂D⁺ (1(1, 0)–1(1, 1)) has been observed in multiple dense starless cores and prestellar cores (see C. Ceccarelli et al. 2014 for a review). P. Caselli et al. (2008) report an o–H₂D⁺ surface density ~10¹²–10¹³ cm⁻² in regions characterized by a gas density ~10⁵–10⁶ cm⁻³ and kinetic temperatures around
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Number density distributions (log scale) from our best-fit model for CO, N<sub>2</sub>H<sup>+</sup>, H<sub>2</sub>D<sup>+</sup>, and DCO<sup>+</sup> from left to right and top to bottom. We have only modeled the ring region around GG Tau A, spanning a radial distance of 193–285 au from the central stellar system. The model disk atmosphere extends up to four scale heights.</figcaption></figure>

10 K, with the most prominent detections in the densest and centrally concentrated cores. However, no clear detection has been obtained in protoplanetary disks. E. Chapillon et al. (2011) derived 3σ upper limits on the o–H<sub>2</sub>D<sup>+</sup> surface density for the T Tauri disks DM Tau and TW Hya. Depending on the physical model, DM Tau’s upper limit lies between 4.5 × 10<sup>11</sup> and 1.9 × 10<sup>12</sup> cm<sup>-2</sup> and for TW Hya, it is between 9.0 × 10<sup>11</sup> and 1.4 × 10<sup>12</sup> cm<sup>-2</sup>. Our upper limit for o–H<sub>2</sub>D<sup>+</sup> in GG Tau is 2.1 × 10<sup>10</sup> cm<sup>-2</sup>. Note that E. Chapillon et al. (2011)’s upper limits are from single-dish observations (1σ noise level of

0.33 Jy km s<sup>-1</sup>). We have achieved a much higher sensitivity (1σ noise level of 0.019 Jy km s<sup>-1</sup>), making our observation the best available constraint on the upper limit of the o–H<sub>2</sub>D<sup>+</sup> surface density in any disk to date.

The N<sub>2</sub>H<sup>+</sup> (3–2) transition has been observed in five disks (three around T Tauri stars: IM Lup, GM Aur, and AS 209, and two around Herbig Ae stars: HD 163296 and MWC 480) as part of the Molecules with ALMA at Planet-forming Scales (MAPS) program (G. Cataldi et al. 2021). Despite the diverse continuum structures of the disks with multiple rings and gaps,
$N_2H^+$ emission showcases a ring morphology (Figure 4 in G. Cataldi et al. 2021). Our $N_2H^+$ emission is also characterized by a single-ring structure (Figure 2). C. Qi et al. (2013) have previously reported $N_2H^+$ detection toward TW Hya, and the single-ring emission structure is also evident there. The $N_2H^+$ ring morphology across different disks can be attributed to radial and vertical temperature variations (C. Qi et al. 2013). In the MAPS survey by G. Cataldi et al. (2021), two different $N_2H^+$ surface density values are derived by fixing the excitation temperature, $T_{ex}$, to (a) the observationally constrained CO freeze-out temperature, ~20 K and (b) the midplane temperatures, $T_{mid}$, of the respective disks. For the T Tauri disks (IM Lup, GM Aur, and AS 209), the $N_2H^+$ surface density peaks around $10^{13} \text{ cm}^{-2}$, while for the Herbig Ae disks (HD 163296 and MWC 480), it reaches up to $10^{12} \text{ cm}^{-2}$. Interestingly, GG Tau, despite being a T Tauri disk, exhibits a $N_2H^+$ surface density of approximately $10^{12} \text{ cm}^{-2}$. The surface density in TW Hya reported by C. Qi et al. (2013) ranges from $4 \times 10^{12}$ to $2 \times 10^{15} \text{ cm}^{-2}$ for different assumed physical structures. G. Cataldi et al. (2021) indicate that temperatures as low as 12 K in the $N_2H^+$-emitting region can also explain the observed flux of the 3-2 transition in the protoplanetary disks (e.g., GM Aur, HD 163296, and MWC 480) under the LTE assumption. Previously, in N. T. Phuong et al. (2021)'s study, $N_2H^+$ in GG Tau A was assumed to be arising from the same layer as CO, following a radial temperature variation, $T(r) = 20(r/250 \text{ au})^{-1}$ K. Our direct measure of $T(r) = 12$ K for the $N_2H^+$-emitting region (Table 4), based on a multiline excitation analysis, suggests that we are tracing a region much closer to the midplane than initially assumed.

A ring-like radial distribution for $DCO^+$ has been commonly observed in multiple disks, including the TW Hya, DM Tau, AS 209, V4046 Sgr, MWC 480, and HD 163296 (C. Qi et al. 2008; G. S. Mathews et al. 2013; R. Teague et al. 2015; J. Huang et al. 2017) and is facilitated by reaction between $H_2D^+$ and CO (e.g., K. Willacy 2007). Contrary to the common scenario, we observe a partially overlapping double-ring radial distribution of $DCO^+$ emission in GG Tau A (refer to Figure 2), which can likely be attributed to different chemical effects in the ring and the outer disk regions owing to the difference in densities between the two regions. This is discussed in detail in Sections 5.2 and 5.3. Notably, K. I. Öberg et al. (2015) observe double $DCO^+$ rings in the disk IM Lup. They found that the $DCO^+$ abundance first decreases radially due to CO depletion, then another ring forms due to nonthermal desorption of CO in low surface density outer regions as a result of higher UV penetration. N. T. Phuong et al. (2018, 2021) have previously assumed that $DCO^+$ emission originates from the same region as $N_2H^+$ and CO, characterized by the radial temperature, $T(r) = 20(r/250 \text{ au})^{-1}$ K in GG Tau A. In this work, we found that $DCO^+$ emission region is best represented by $T(r) = 16(r/250 \text{ au})^{-0.39}$ K (Table 4). Hence, we expect that $DCO^+$ emission region lies vertically between the CO-emitting region and $N_2H^+$-emitting region and closer to the midplane than previously thought.

## 5.2. Chemical Trends in the Ring Region

We present a subset of chemical models for a selected set of input parameters in Figures 7 and 8 from the full model grids presented in Table 8. Only the CRI rates and the C/O ratios were found to have any significant impact on the modeled

surface densities of the targeted molecules. Other parameters cannot be constrained by the current data.

CRI. Figure 7 shows the variation of the modeled abundances with CRI rate while keeping the C/O ratio constant at 0.7, and the remaining parameters are set to the best-fit values. The $N_2H^+$, $DCO^+$, and $H_2D^+$ abundances are quite sensitive to the CRI rate, while the UV flux values have no impact whatsoever. This is expected because these molecules probe a high-density, low-temperature, CO-depleted region near the midplane where UV penetration is minimal. Hence, cosmic rays are the primary ionization source for this cold molecular chemistry. However, the CRI rate has barely any impact on S-bearing species. It is worth noting that the standard interstellar CRI rate of $10^{-17} \text{ s}^{-1}$ predicts an $H_2D^+$ surface density above the upper limit and a subinterstellar CRI rate of $10^{-18} \text{ s}^{-1}$ explains the observationally constrained surface densities the best. This low CRI rate can be explained by a rise in the number of collisions experienced by cosmic-ray particles due to the greater effective surface density as they gyrate toward the disk midplane (M. Padovani et al. 2013). The attenuation is especially prominent in our case, as our chemistry focuses on high density near midplane regions. L. I. Cleeves et al. (2015) have previously inferred a very low CRI ($10^{-19} \text{ s}^{-1}$) in the protoplanetary disk TW Hya, which they suggest is due to an exclusion by a young stellar wind. Our current astrochemical model only includes a constant CRI rate across the disk. The inclusion of short-lived radionuclides (L. I. Cleeves et al. 2013) and a vertical distribution of CRI rates (M. Padovani et al. 2018) may provide a better picture.

C/O ratio. We have also modeled the abundances of our targeted molecules for different elemental C/O ratios from 0.5 to 1.2 by changing the initial O abundances, as shown in Figure 8. Changing the C/O ratio from 0.5 to 1.0 decreases the $N_2H^+$ surface density by a factor of approximately 2.5, while the $DCO^+$ surface density increases by at most a factor of 1–2. However, when the C/O ratio was further increased to 1.2, changes in the surface densities of $N_2H^+$ and $DCO^+$ became significant (by a factor of approximately 10). The $H_2D^+$ surface density remains unaffected for any change in the C/O ratio from 0.5 to 1.2. The chemistry of $^{13}$CS and $SO_2$, for which we obtained upper limits from our observations, is quite sensitive to changes in the C/O ratio. The general trend is that O-poor chemistry favors the formation of the C–S bond, and for $SO_2$, it is the other way around (see the discussions by D. Semenov et al. 2018; R. Le Gal et al. 2021, for further details). Our observations are best represented by a CRI rate of $10^{-18} \text{ s}^{-1}$ and an elemental C/O ratio of 1.0. At the commencement of the evolution of a molecular cloud, a significant amount of important elements, including carbon and oxygen, exist in a refractory state, which is expected to shape the chemical budget of the protoplanetary disk formed from that cloud. Based on A. C. A. Boogert et al. (2015)'s observation of ice abundances, K. I. Öberg & E. A. Bergin (2021) report that about a 40% O chemical budget was unaccounted for at the onset of star and planet formation, leading to a poor constraint on the initial C and O gas-phase abundances. Our elemental C/O ratio of 1.0 suggests an oxygen-depleted environment in the dense ring. Similar inferences of high C/O ratios (0.9–2.0) were also found in more typical protoplanetary disks (e.g., R. Le Gal et al. 2021; M. Ruaud et al. 2022), which may have consequences on the formation of planets in the inner regions.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Comparison between models with different CRI rates ($\zeta$) represented by the solid colored lines in the figure. Here, the “x” points represent the observationally contained surface density values from our DISKFIT analysis. The dashed horizontal lines represent observationally constrained upper limits on the surface density for the undetected molecules. The gray horizontal patch denotes the surface density constraint from a different transition of $^{13}$CS ($J=2-1$) by N. T. Phuong et al. (2021). Best model: CRI rate = $10^{-18}$ s$^{-1}$ ($\chi^2_{red}$ = 5.01). Here, C/O = 0.7 and the rest of the parameters are set to the best-fit parameters. The upper limit on the $^{13}$CS surface density is calculated from the CS surface density upper limit by assuming $N_{CS}/N_{^{13}CS}$ = 100 (N. T. Phuong et al. 2021).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Comparison between models with different initial C/O ratios represented by the solid colored lines in the figure. The rest of the denotations are the same as Figure 7. Best model: C/O = 1.0 ($\chi^2_{red}$ = 2.78). All the parameters other than C/O are set to the best-fit parameters. The upper limit on the $^{13}$CS surface density is calculated from the CS surface density upper limit by assuming $N_{CS}/N_{^{13}CS}$ = 100 (N. T. Phuong et al. 2021).</figcaption></figure>

**Spatial distributions.** Figure 6 represents the molecular number density distributions for our best model. We see a DCO$^+$ distribution in those regions where CO is available in the gas phase in small amounts. Observationally, DCO$^+$ follows a ring-like emission distribution in the ring region (refer to Figures 2 and B1). In our modeled abundances, DCO$^+$ is primarily peaking

in temperature regions representative of the warm molecular layer (2.5–3.5 $H_r$), somewhat showcasing the ring-like distribution. One of the dominant ions in the warm molecular layer is HCO$^+$. It essentially forms through the protonation of CO and, in turn, produces DCO$^+$ through isotope exchange. The major reaction for DCO$^+$ production contributing to our observation above 2.5 $H_r$ is
$HCO^+ + D \rightarrow DCO^+ + H$. We also see some $DCO^+$ arising below the warm molecular layer (~0.5–2.5 $H_r$) produced through a different channel, $CO + N_2D^+ \rightarrow N_2 + DCO^+$. Our DISKFIT analysis locates the peak of the $DCO^+$ ring at a radial distance of 210 ± 50 au ($R_1$ in Table 4), and this can be attributed to the fact that both $DCO^+$ formation channels are contributing to the $DCO^+$ abundances in this radial range. Near the disk midplane, $DCO^+$ starts depleting because its parent molecule, CO, efficiently freezes onto the grain surfaces owing to the high density and low temperature.

The $N_2H^+$ modeled distribution lies just below the $DCO^+$ layer, with its number density peaking in regions where $DCO^+$ starts depleting. Our observationally constrained representative temperatures of 16 and 12 K for the $DCO^+$- and $N_2H^+$-emitting regions, respectively (Table 4), suggest the same. $N_2H^+$ is abundant in the midplane at the inner edge of the ring as well as in the warm molecular layer (up to $3.0 H_r$) at the outer edge of the ring. All throughout the ring region, the major formation pathway of $N_2H^+$ is the protonation of $N_2$, i.e., $N_2 + H_3^+ \rightarrow H_2 + N_2H^+$. So, the chemistry is mainly controlled by the gas-phase abundances of CO, $N_2$, HD, and $H_3^+$. As long as CO remains in the gas phase, it destroys $N_2H^+$ to form $HCO^+$. When CO starts depleting, $N_2H^+$ start increasing, as $N_2$ is still in the gas phase owing to its lower binding energy with respect to CO, forming $N_2H^+$, and there is less destruction of $N_2H^+$ because of the reduced CO presence. Our model predicts a decreasing $N_2H^+$ surface density with increasing radius, contrary to the observed increase within the continuum ring. This discrepancy arises primarily from the model's overestimation of $N_2$ depletion from the gas phase. A similar trend was observed in the protoplanetary disk around DM Tau, where V. Wakelam et al. (2019) conducted sensitivity analyses on various model parameters. One limitation of our model is the absence of a grain growth mechanism and the assumption that grains are perfectly spherical. A. Kataoka et al. (2013) proposed a dust growth mechanism in protoplanetary disks, which would create inhomogeneities in the dust distribution and potentially allow more UV penetration toward the disk's midplane. UV penetration, in turn, facilitates photo-desorption of $N_2$. Chemically, nonspherical grains would have a larger cross-sectional area for collisions with gas-phase species and a larger surface area, increasing the likelihood of $N_2$ desorption. Additionally, our models do not include X-ray chemistry, which could also affect the abundance of $N_2H^+$ (L. I. Cleeves et al. 2015; R. Teague et al. 2015).

Our astrochemical model suggests that $H_2D^+$ is abundant in the midplane region, situated just below the $N_2H^+$ layer. The major formation pathway predicted by our model is $H_3^+ + HD \rightarrow H_2D^+ + H_2$. As $N_2$ begins to freeze onto grain surfaces, more $H_3^+$ become available for HD to capture, and we see a rise in $H_2D^+$ abundance and a fall in $N_2H^+$ abundance.

We did not detect $H_2D^+$. In our effort to replicate the observed surface density for the detected molecules, we realized that a lower CRI rate may provide the best explanation for our observations. $H_2D^+$ abundances are regulated by the availability of $H_3^+$, which in turn is dependent on the rate of ionization of $H_2$. A lower CRI rate results in reduced $H_2^+$ production ($H_2 + \text{cosmic ray} \rightarrow H_2^+ + e^-$), leading to a lower abundance of $H_3^+$ ($H_2 + H_2^+ \rightarrow H_3^+ + H$). Hence, a lower ionization rate at the midplane, owing to the very high density in the ring region, may explain the nondetection of $\text{o-H}_2D^+ (1(1, 0)-1(1, 1))$. Indeed, our model using a CRI rate ($\zeta_{CR}$) of $10^{-18} \text{ s}^{-1}$ predicts an $H_2D^+$ surface density of the order of $\sim 10^{12} \text{ cm}^{-2}$ close to the current upper limit reported in Table 5. This limit suggests we are close to detection, and increasing the

observation time should assure our chance of detection for $H_2D^+$. However, the upper limit on total $H_2D^+$ is based on the assumption of thermalized ortho-to-para ratio (LTE), with a high correction Boltzmann factor (see Section 3.2). While thermalization is facilitated by the very high midplane density, spin-dependent reaction rates may lead to a different prediction (e.g., E. Chapillon et al. 2011).

### 5.3. Possible Chemistry in the Outer Disk

$N_2H^+$ (4–3) and $DCO^+$ (5–4) emission is also spatially extended beyond the ring, reaching up to a radial distance of approximately 550 au (Figures 1, 2). The velocity-integrated radial profiles of $N_2H^+$ (4–3) and $DCO^+$ (5–4) (Figure 2) indicate that the former is decreasing as we go radially outward in the outer disk, while the latter shows another peak at around ~375 au.

We have not modeled the outer disk due to its complex structure (Y.-W. Tang et al. 2016; N. T. Phuong et al. 2020b). However, based on our current understanding, the outer disk is characterized by a lower density and colder temperatures as compared to the ring (refer to Section 4.1). Lower densities outside the dense ring allow higher cosmic-ray penetration, leading to nonthermal desorption of CO and $N_2$. Consequently, an increased presence of CO and $N_2$ in the gas phase is anticipated in the outer disk, suggesting that the primary pathway for $DCO^+$ production would be through $H_2D^+ + CO$ (N. T. Phuong et al. 2018). On the other hand, $N_2H^+$ formed through protonation of $N_2$ will be readily destroyed in the presence of CO, which will cause a rapid decrease in the $N_2H^+$ abundance as seen in Figure 2.

### 5.4. Temperature Trends in the Ring Region

The DISKFIT analysis of our observations of $N_2H^+$ (4–3) and $DCO^+$ (5–4) along with the NOEMA observations of $N_2H^+$ (1–0) and $DCO^+$ (1–0, 3–2) by N. T. Phuong et al. (2018, 2021) indicate that $N_2H^+$ and $DCO^+$ are best represented by average temperatures of 12 and 16 K (Table 4), respectively. These rotation temperatures are consistent with $N_2H^+$ being located closer to the cold disk midplane than $DCO^+$ and a fortiori, CO and its isotopologues previously studied by N. T. Phuong et al. (2020b).

However, in our chemical model, while the distributions of $N_2H^+$ and $DCO^+$ extend to regions with temperatures of 12 and 16 K, their density-weighted average temperatures are 21 K and 24 K, respectively. The discrepancy with the observed temperatures can be attributed to our limited understanding of the physical structure of the ring region and the simplifications applied in the model to derive the physical structure. $DCO^+$ abundances are sensitive to grain size, and the dust-to-gas mass ratio significantly influences gas-phase molecular abundances. S. Gavino et al. (2021, 2023) recently investigated the impact of the dependence of the grain size on their temperature. They found that large grains exhibit lower temperatures than smaller ones, with an effect on chemistry that cannot be properly represented by a single grain size. This can influence the structure of the CO distribution (with more CO near the midplane) and, consequently, the gas-phase abundances of $DCO^+$ and $N_2H^+$ in the densest, coldest regions.

## 6. Conclusions

In this paper, we have presented high-sensitivity ALMA Band 7 observations of the circumstellar disk around the triple star system GG Tau A, which highlight CO-regulated cold molecular chemistry in the shielded high-density regions near
the disk midplane. Our key findings and conclusions are summarized below.

1. With Band 7 single-pointing observation of GG Tau A, we achieved a spectral resolution of 0.1 km s⁻¹ and an angular resolution of 0.′′8–0.′′9. We spatially resolved the emission from N₂H⁺ (4–3) and DCO⁺ (5–4) while o–H₂D⁺ (1, 0)–1 (1, 1), ¹³CS (8–7), and SO₂ 19(4, 16)–19(3, 17) show no detection.

2. Comparing the predictions of a parametric disk model to our observed visibilities along with the lower excitation transition observations obtained with NOEMA by N. T. Phuong et al. (2018) and N. T. Phuong et al. (2020b), we constrain the average temperature of the N₂H⁺- and DCO⁺-emitting regions to be around 12 K and 16 K, respectively. These are at odds with the predictions from all the chemical models we explored and all the previous studies. Chemical models suggest DCO⁺ peaks at 20–25 K, where some CO remains in the gas phase, and N₂H⁺ peaks around 20 K, where CO is fully frozen onto grains. On the contrary, our observational constraints indicate that both molecules trace regions below the CO layer (N. T. Phuong et al. 2020a), much closer to the midplane than previously thought, with N₂H⁺ being the closest.

3. Our radial distribution of DCO⁺ (5–4) emission reveals an atypical partially overlapping double-ring structure. The inner ring peaks at 210 ± 50 au within the continuum ring region, while the more intense outer ring peaks at 375 ± 10 au in the outer disk. We hypothesize that DCO⁺ forms via two different pathways in the ring and the outer disk regions depending on the gas-phase CO abundance. In the ring region, DCO⁺ forms through isotope exchange after HCO⁺ arises from CO protonation. However, constraining DCO⁺ formation in the outer disk regions requires future work that takes into account the complex physical structure.

4. N₂H⁺ (4–3) emission shows a ring-like radial distribution facilitated through protonation of N₂ with its emission layer situated right below the DCO⁺ layer. The anticorrelation between the N₂H⁺ distribution and the CO distribution is clearly visible in our modeled simulations.

5. We ran numerous chemical models varying the CRI rates, C/O ratios, and stellar UV fluxes. These models suggest that the chemistry of H₂D⁺, N₂H⁺, and DCO⁺ in cold, dense regions is primarily influenced by the gas-phase CO abundance and the dominant midplane ion, H₃⁺. These regions are shielded from external UV radiation, making CRI the most crucial factor in determining the CO and H₃⁺ abundances and, consequently, the distributions of N₂H⁺, DCO⁺, and H₂D⁺. On the other hand, the chemistry of ¹³CS and SO₂ is very sensitive to the C/O ratio.

6. After running a grid of astrochemical models to simultaneously fit the surface density constraints from the observed transitions along with the upper limits provided by the nondetections, our best-fit model, with a C/O ratio of 1.0 and a subinterstellar CRI rate of 10⁻¹⁸ s⁻¹, successfully reproduced the radial trends in surface density for DCO⁺ but slightly overestimated both the N₂H⁺ content and the expected rotation temperatures. The CRI attenuation in the high-density GG Tau A ring can be attributed to an increased number of collisions experienced by cosmic-ray particles as they gyrate toward the midplane.

7. The nondetection of o–H₂D⁺ (1(1, 0)–1(1, 1)) suggests that ionization is much lower than the standard interstellar rate in the disk midplane due to its very high density. In the cold midplane region, H₂D⁺ is directly dependent on the abundance of H₃⁺, which, in turn, is regulated by the ionization rate. The nondetection of H₂D⁺ aligns with our best model, where the expected surface density is below our observationally determined 3σ upper limit only for a subinterstellar CRI rate.

The GG Tau A disk is very complex, but it is one of the nearest (hence larger) and more massive disks found around a low-mass star system. Such properties should facilitate the detection of o–H₂D⁺ (1(1, 0)–1(1, 1)). However, its modeled surface density being close to the observationally constrained 3σ upper limit suggests that the detection of H₂D⁺ is within ALMA’s capabilities, although at the cost of a significantly longer integration time.

## Acknowledgments

This paper makes use of the following ALMA data: ADS/JAO ALMA #2021.1.00342.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. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under a cooperative agreement by Associated Universities, Inc. L.M. acknowledges the financial support of DAE and DST-SERB research grants (SRG/2021/002116 and MTR/2021/000864) of the Government of India. A.D. and S.G. thank the French CNRS programs PNP, PNPS, and PCMI. This research was carried out in part at the Jet Propulsion Laboratory, which is operated for NASA by the California Institute of Technology. K.W. acknowledges the financial support from the NASA Emerging Worlds grant 18-EW-182-0083. Th.H. and D.S. acknowledge support from the European Research Council under the Horizon 2020 Framework Program via the ERC Advanced Grant Origins 83 24 28 (PI: Th. Henning). N.T. P. acknowledges the financial support of Vingroup Innovation Foundation (VINIF) under project code VINIF.2023.DA.057". A.C. received financial support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (ERC Starting Grant “Chemtrip,” grant agreement No 949278). R.S. thanks NASA for financial support through various Astrophysics Data Program awards. A.K. acknowledges support from the NKFIH excellence grant TKP2021-NKTA-64. We would like to thank the anonymous referees for their constructive comments, which helped improve the manuscript.

Facility: ALMA.

Software: CASA (J. P. McMullin et al. 2007), kepler-ian_mask (R. Teague 2020), GoFish (R. Teague 2019), VISIBLE (R. A. Loomis et al. 2018b), dynasty (J. S. Speagle 2020), IMAGER (https://imager.oasu.ubordeaux.fr), DNAUTILUS (L. Majumdar et al. 2017), numpy (C. R. Harris et al. 2020), astropy (Astropy Collaboration et al. 2013, 2018, 2022), matplotlib (J. D. Hunter 2007).

## Appendix A

### Channel Maps, Moment Maps, and Radial Profiles

The generated channel maps for the detected transitions, N₂H⁺ (4–3) and DCO⁺ (5–4), are shown in Figures A1 and A2,
respectively. Figure A3 shows integrated intensity maps between the velocity range 5.0 and 8.0 km s⁻¹ for all the targeted transitions. Sigma clipping and masking are not performed while generating these integrated intensity maps. Figure A4 shows the radial distribution of azimuthally averaged emission for the

targeted transitions and the continuum generated from the integrated intensity maps shown in Figure A3. We have displayed the matched filter responses for all the targeted transitions in Figure A5. The methods for generating all these figures are explained in Section 2.2.

<figure><img src="image_11.png" /><figcaption>**Figure A1.** Channel maps of N₂H⁺ (4-3) imaged with the tclean task in CASA as explained Section 2. The channel velocities (LSRK) are displayed in the top right corner of each channel in units of km s⁻¹. The solid contours are [3, 5, 7, 9,...] × σ levels, where σ = 14.6 mJy beam⁻¹. The ellipse at the lower left corner indicates the beam size of 0.′93 × 0.′77 with a PA of 68.76. The color bar is stretched nonlinearly to showcase the line emissions more prominently.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure A2.** Channel maps of DCO$^{+}$ (5–4) imaged with the tclean task in CASA as explained Section 2. The channel velocities (LSRK) are displayed in the top right corner of each channel in units of km s$^{-1}$. The solid contours are [3, 5, 7, 9,...] × σ levels, where σ = 5.6 mJy beam$^{-1}$. The ellipse at the lower left corner indicates the beam size of 0.′′91 × 0.′′73 with a PA of 99$^{◦}$20. The color bar is stretched nonlinearly to showcase the line emissions more prominently.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure A3.** Integrated intensity maps of the targeted transitions generated without sigma clipping and masking. The overplotted elliptical ring contours in the white are at 193 and 285 au to indicate the continuum ring region where 90% of the circumtertiary emission arises (S. Guilloteau et al. 1999). The panel at the bottom right corner is the continuum image at 358 GHz.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure A4.** The dashed dot and solid curves represent azimuthally averaged velocity-integrated radial profiles generated from integrated intensity maps with and without Keplerian masking, respectively. The shaded region in the above figure corresponds to 1σ rms. The gray horizontal line at the top right denotes the beam major axis. The dashed vertical lines are to designate the location of the ring around GG Tau A, spanning from 193 to 285 au. The subplot at the bottom right corner is the radial distribution of the azimuthally averaged continuum emission at 358 GHz.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure A5.** Matched filter responses to the emission for the targeted transitions, written at the upper left corner. The vertical gray dashed line is plotted at source velocity $v_{LSRK} = 6.4 \text{ km s}^{-1}$ for each subplot. Only $N_2H^+$ (4–3) and DCO$^+$ (5–4) show clear detections; the rest of them do not.</figcaption></figure>

## Appendix B
IMAGER and DISKFIT Analysis

Figure B1 shows the Keplerian deprojection produced by the IMAGER command *KEPLER*. Images used for this purpose were produced with IMAGER and have an angular resolution of 0.′′66 × 0.′′53 at a PA of 83°. Spectral resolution prior to Keplerian velocity correction was about 0.1 km s⁻¹. This

correction requires resampling of the spectra on a different grid, which results in correlated channels; the correlation factor is estimated from the noise autocorrelation spectrum. The best-fit radial profile derived from DISKFIT is always slightly above the measured one because the deconvolution with a limited signal-to-noise ratio cannot recover the full line brightness, while the uv plane fitting method is not biased in this respect.
<figure><img src="image_16.png" /><figcaption>**Figure B1.** Kepler plots generated for the observed transitions of N₂H⁺ (4–3) (top) and DCO⁺ (5–4) (bottom). For each molecule, the map is the brightness temperature (in units of K) as a function of velocity and radius (“teardrop” plot). The left panel shows the peak brightness radial profile, with error bars indicated, and the red profile is the one derived from the best-fit DISKFIT model. The top panel shows the integrated spectrum over the range indicated by the blue line in the radial profile, with the best-fit Gaussian spectrum superimposed in red. The peak (S, in units of Jy) and integrated (∫S, dV, in units of Jy km s⁻¹) flux densities, central velocity (V, in units of km s⁻¹), and FWHM line width (ΔV, in units of km s⁻¹) obtained from a Gaussian fit to the integrated spectrum are given in the upper left corner of each panel.</figcaption></figure>
Appendix C

Simple Local Thermodynamic Equilibrium Approach

We have undertaken a simpler approach to calculate disk-averaged surface densities where we assume LTE considering the disk densities are relatively high compared to critical densities for the studied transitions. This method requires prior knowledge of rotational temperature for a single transition.

For optically thin emission, the surface density of molecules in the upper state of a specific transition is expressed as

$$N_u^{\text{thin}} = \frac{4\pi}{hcA_{ul}} \frac{\left(\int S_{\nu} dv\right)}{\Omega}, \qquad (C1)$$

$$= \frac{8\pi k \nu^2}{hc^3 A_{ul}} \left( \int T_b dv \right), \qquad (C2)$$

where $\int S_{\nu} dv$ is the velocity-integrated flux density, $A_{ul}$ is the Einstein coefficient, $\Omega$ denotes the solid angle subtended by the emission area, and $h$ and $c$ stand for the Planck constant and the speed of light in a vacuum, respectively (P. F. Goldsmith & W. D. Langer 1999). In the alternate form, $\nu$ is the line frequency and $T_b$ is the line brightness temperature.

The total surface density can be derived from the upper state surface density, assuming a Boltzmann distribution

$$N_T^{\text{thin}} = N_u^{\text{thin}} \frac{Q(T_{\text{rot}})}{g_u} \exp\left(\frac{E_u}{T_{\text{rot}}}\right), \qquad (C3)$$

where $g_u$ is the degeneracy, $E_u$ is the energy of the upper energy level in units of K, and $Q(T_{\text{rot}})$ is the partition function at the rotation temperature $T_{\text{rot}}$, which can be computed by summation over the energy levels or by interpolating tabulated $Q$ values for discrete $T_{\text{rot}}$ values found in the Cologne Database for Molecular Spectroscopy (C. P. Endres et al. 2016).

If the transition is optically thick, the above value gives a strict lower limit to the observed surface density. Opacity corrections are possible in a purely homogeneous medium with only turbulent or thermal line broadening

(P. F. Goldsmith & W. D. Langer 1999)

$$N_T = N_T^{\text{thin}} \frac{\tau}{1 - e^{-\tau}}, \qquad (C4)$$

where we can self-consistently calculate the optical depth, $\tau$ assuming a Gaussian-like line profile $\sim 1/\Delta\nu$ at the line center (P. F. Goldsmith & W. D. Langer 1999), $\Delta\nu = 0.25 \text{ km s}^{-1}$ being the intrinsic FWHM line width constrained from the DISKFIT analysis

$$\tau = \frac{A_{ul} c^3}{8\pi \nu^3 \Delta\nu} N_u \left[ \exp\left(\frac{h\nu}{kT_{\text{rot}}}\right) - 1 \right], \qquad (C5)$$

Unfortunately, in our case, the abovementioned correction does not hold while calculating the radially resolved surface density, even if we apply this to radially dependent emission. Because of the Keplerian shear, there are significant regions in our images where the line-of-sight opacity varies both as a function of azimuth and radius, with both directions being insufficiently well resolved at our angular resolution (about 100 au). Accordingly, the azimuthally averaged brightness does not originate from a unique opacity below about 400 au.

However, the effect is expected to be small. The peak brightness of $N_2H^+$ (see Figure B1) is 2 K, while the brightness temperature of an optically thick line would be $J_{\nu}(T_{\text{rot}}) = 5.3$ K for $T_{\text{rot}} = 12$ K at the $J = 4-3$ line frequency. This indicates a peak line opacity of about 0.47, leading to a $\sim 22\%$ opacity correction at most. For DCO$^+$, whose emission is about 3 times fainter, the correction is even smaller, as for the (undetected) emission from $H_2D^+$ and other molecules. We thus expect the constraint given by the optically thin approximation to be a reasonable estimate for the mean surface density.

We explore the likely ranges of disk-averaged surface density using nested sampling while assuming lognormal priors for both the molecules in the range $10^5 < N_T < 10^{18} \text{ cm}^{-2}$ and fixing the rotational temperature to the values constrained in the DISKFIT analysis. The rotation temperatures for $N_2H^+$ and DCO$^+$ were set at 12 K and 16 K, respectively. A value of 15 K was used for the undetected molecules. The Python

<figure><img src="image_17.png" /><figcaption>Figure C1. Posterior distributions for disk-averaged surface density calculations using the simple LTE approach. The left panel is for the molecule $N_2H^+$, and the right one is for DCO$^+$. The gray dashed lines represent 16th, 50th, and 84th percentiles of the distribution. The median value with $1\sigma$ error bars (68% confidence interval) are reported at the top right.</figcaption></figure>
package utilized to sample the posterior distribution is *dynesty* (J. S. Speagle 2020). We consider the median of the posterior as the best-fit value, while the 16th and 84th percentile values account for the 1σ uncertainties. Posterior distributions for these averaged surface density derivations are shown in Figure C1.

## Appendix D

### Grids of Astrochemical Models and Finding the Best-fit Model

We ran a large number of astrochemical models with various combinations of parameters. The grids of models are summarized in Table 8. Note that we have not implemented any grain growth in our models due to poor constraints on larger grain sizes. We explored the parameter space using the least squares method while simultaneously considering the observational constraints on the surface densities for N₂H⁺ and DCO⁺ from our observations and ¹³CS as reported by N. T. Phuong et al. (2021). We ensured that the modeled values for undetected molecules remained below their respective upper limits. To quantify the goodness of fit, we calculated reduced chi-squared values. This metric evaluates

the agreement between the observed and modeled surface densities across multiple molecules and radii. The reduced chi-squared calculation is defined as follows

$$ \chi_{\text{red}}^2 = \frac{1}{\nu} \sum_j \left( \sum_i \left( \frac{N_{\text{obs}}(i, j) - N_{\text{model}}(i, j)}{\sigma_{\text{obs}}(i, j) + \sigma_{\text{model}}(i, j)} \right)^2 \right). \quad (\text{D1}) $$

In the equation, $\nu$, the number of degrees of freedom is 25. $N_{\text{obs}}(i, j)$ and $N_{\text{model}}(i, j)$ denote the observed and modeled surface densities of the $j$th molecule at the $i$th radial point, respectively. The terms $\sigma_{\text{obs}}(i, j)$ and $\sigma_{\text{model}}(i, j)$ represent the associated error estimates for the observed and modeled values, respectively. Observational errors were derived from the error bars in the peak column density values in Table 4, while the modeled errors were assumed to be 20% to account for uncertainties in the chemical models due to rate coefficients and other chemical parameters. The best-fit model is the one with the lowest $\chi_{\text{red}}^2$ value. Although $\chi_{\text{red}}^2$ depends on the assumed errors, the best-fit solution was found to be independent of this assumption.

<table><caption>**Table 8**
Parameter Space Explored with DNAUTILUS</caption><thead><tr><th>Initial Abundances</th><th>Settling</th><th>Grain Growth</th><th>f<sub>UV, 214 au</sub> (χ<sub>0</sub>)<sup>a</sup></th><th>C/O<sup>b</sup></th><th>ζ<sub>CR</sub> <sup>c</sup><br>(s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Atomic<sup>d</sup></td><td>no</td><td>no</td><td>375, 1500, 3000</td><td>0.5, 0.7, 0.9, 1.0, 1.2</td><td>10<sup>-16</sup>, 10<sup>-17</sup>, 10<sup>-18</sup>, 10<sup>-19</sup>, 10<sup>-20</sup></td></tr><tr><td>10<sup>6</sup> yr Molecular Cloud<sup>e</sup></td><td>no</td><td>no</td><td>375, 1500, 3000</td><td>0.5, 0.7, 0.9, 1.0, 1.2</td><td>10<sup>-16</sup>, 10<sup>-17</sup>, 10<sup>-18</sup>, 10<sup>-19</sup>, 10<sup>-20</sup></td></tr><tr><td>Atomic<sup>d</sup></td><td>yes</td><td>no</td><td>375, 1500, 3000</td><td>0.5, 0.7, 0.9, 1.0, 1.2</td><td>10<sup>-16</sup>, 10<sup>-17</sup>, 10<sup>-18</sup>, 10<sup>-19</sup>, 10<sup>-20</sup></td></tr><tr><td>10<sup>6</sup> yr Molecular Cloud<sup>e</sup></td><td>yes</td><td>no</td><td>375, 1500, 3000</td><td>0.5, 0.7, 0.9, 1.0, 1.2</td><td>10<sup>-16</sup>, 10<sup>-17</sup>, 10<sup>-18</sup>, 10<sup>-19</sup>, 10<sup>-20</sup></td></tr></tbody></table>

**Notes.**

<footnote><sup>a</sup> *f*<sub>UV</sub> is the stellar UV flux at the reference radius, 214 au in the units of the Draine interstellar UV field, χ<sub>0</sub>.</footnote>

<footnote><sup>b</sup> C/O ratios have been considered around the standard value of 0.7 to account for low depletion and high depletion of oxygen (L. Reboussin et al. 2015).</footnote>

<footnote><sup>c</sup> Different CRI rates ζ<sub>CR</sub> have been explored to account for how they impact our chemistry.</footnote>

<footnote><sup>d</sup> In the reset scenario, we have taken our initial abundances in atomic forms along with H₂ and HD and run our astrochemical model for 10<sup>6</sup> yr.</footnote>

<footnote><sup>e</sup> We have also considered the inheritance scenario as another possibility, where we have run a cloud model with atomic initial abundances, and the final abundances after 10<sup>6</sup> yr are used as the initial abundance in the disk model. The disk model was then run for 10<sup>6</sup> yr.</footnote>
ORCID iDs

Parashmoni Kashyap https://orcid.org/0009-0008-2350-5876

Liton Majumdar https://orcid.org/0000-0001-7031-8039

Anne Dutrey https://orcid.org/0000-0003-2341-5922

Stéphane Guilloteau https://orcid.org/0000-0003-3773-1870

Karen Willacy https://orcid.org/0000-0001-6124-5974

Edwige Chapillon https://orcid.org/0009-0003-6480-8084

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

Dmitry Semenov https://orcid.org/0000-0002-3913-7114

Thomas Henning https://orcid.org/0000-0002-1493-300X

Neal Turner https://orcid.org/0000-0001-8292-1943

Raghvendra Sahai https://orcid.org/0000-0002-6858-5063

Ágnes Kóspál https://orcid.org/0000-0001-7157-6275

Audrey Coutens https://orcid.org/0000-0003-1805-3920

V. Piétu https://orcid.org/0009-0006-3497-397X

Pierre Gratier https://orcid.org/0000-0002-6636-4304

Maxime Ruaud https://orcid.org/0000-0003-0522-5789

N. T. Phuong https://orcid.org/0000-0002-4372-5509

E. Di Folco https://orcid.org/0009-0009-9618-4927

Chin-Fei Lee https://orcid.org/0000-0002-3024-5864

Y.-W. Tang https://orcid.org/0000-0002-0675-276X

References

Aikawa, Y., van Zadelhoff, G. J., van Dishoeck, E. F., & Herbst, E. 2002, A&A, 386, 622

Anderson, D. E., Cleeves, L. I., Blake, G. A., et al. 2022, ApJ, 927, 229

Andrews, S. M., Chandler, C. J., Isella, A., et al. 2014, ApJ, 787, 148

Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

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

Bergin, E. A., Alves, J., Huard, T., & Lada, C. J. 2002, ApJL, 570, L101

Boogert, A. C. A., Gerakines, P. A., & Whittet, D. C. B. 2015, ARA&A, 53, 541

Brauer, R., Pantin, E., Di Folco, E., et al. 2019, A&A, 628, A88

Caselli, P., Vastel, C., Ceccarelli, C., et al. 2008, A&A, 492, 703

Caselli, P., Walmsley, C. M., Tafalla, M., Dore, L., & Myers, P. C. 1999, ApJL, 523, L165

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

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

Ceccarelli, C., Dominik, C., Lefloch, B., Caselli, P., & Caux, E. 2004, ApJL, 607, L51

Chapillon, E., Parise, B., Guilloteau, S., & Du, F. 2011, A&A, 533, A143

Cleeves, L. I., Adams, F. C., & Bergin, E. A. 2013, ApJ, 772, 5

Cleeves, L. I., Bergin, E. A., Qi, C., Adams, F. C., & Öberg, K. I. 2015, ApJ, 799, 204

Dartois, E., Dutrey, A., & Guilloteau, S. 2003, A&A, 399, 773

Draine, B. T. 1978, ApJS, 36, 595

Drozdovskaya, M. N., Walsh, C., van Dishoeck, E. F., et al. 2016, MNRAS, 462, 977

Du, F., & Bergin, E. A. 2014, ApJ, 792, 2

Dutrey, A., di Folco, E., Guilloteau, S., et al. 2014, Natur, 514, 600

Dutrey, A., Henning, T., Guilloteau, S., et al. 2007, A&A, 464, 615

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

Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1

Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1

Gammie, C. F. 1996, ApJ, 457, 355

Gavino, S., Dutrey, A., Wakelam, V., et al. 2021, A&A, 654, A65

Gavino, S., Kobus, J., Dutrey, A., et al. 2023, A&A, 680, A59

Goldsmith, P. F., & Langer, W. D. 1999, ApJ, 517, 209

Graedel, T. E., Langer, W. D., & Frerking, M. A. 1982, ApJS, 48, 321

Guilloteau, S., Dutrey, A., & Simon, M. 1999, A&A, 348, 570

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

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

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

Jenkins, E. B. 2009, ApJ, 700, 1299

Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013, A&A, 554, A4

Le Gal, R., Öberg, K. I., Teague, R., et al. 2021, ApJS, 257, 12

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

Loomis, R. A., Öberg, K. I., Andrews, S. M., et al., 2018b VISIBLE: VISIBILITY Based Line Extraction, Astrophysics1133 Source Code Library, ascl:1802.006

Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603

Majumdar, L., Gratier, P., Ruaud, M., et al. 2017, MNRAS, 466, 4470

Mathews, G. S., Klaassen, P. D., Juhász, A., et al. 2013, A&A, 557, A132

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

Neufeld, D. A., Wolfire, M. G., & Schilke, P. 2005, ApJ, 628, 260

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

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

Öberg, K. I., Qi, C., Wilner, D. J., & Andrews, S. M. 2011, ApJ, 743, 152

Öberg, K. I., van Broekhuizen, F., Fraser, H. J., et al. 2005, ApJL, 621, L33

Okuzumi, S., Tanaka, H., Kobayashi, H., & Wada, K. 2012, ApJ, 752, 106

Padovani, M., Hennebelle, P., & Galli, D. 2013, A&A, 560, A114

Padovani, M., Ivlev, A. V., Galli, D., & Caselli, P. 2018, A&A, 614, A111

Phuong, N. T., Chapillon, E., Majumdar, L., et al. 2018, A&A, 616, L5

Phuong, N. T., Dutrey, A., Chapillon, E., et al. 2021, A&A, 653, L5

Phuong, N. T., Dutrey, A., Di Folco, E., et al. 2020a, A&A, 635, L9

Phuong, N. T., Dutrey, A., Diep, P. N., et al. 2020b, A&A, 635, A12

Piétu, V., Dutrey, A., & Guilloteau, S. 2007, A&A, 467, 163

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

Qi, C., Wilner, D. J., Aikawa, Y., Blake, G. A., & Hogerheijde, M. R. 2008, ApJ, 681, 1396

Reboussin, L., Wakelam, V., Guilloteau, S., Hersant, F., & Dutrey, A. 2015, A&A, 579, A82

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

Ruaud, M., Gorti, U., & Hollenbach, D. J. 2022, ApJ, 925, 49

Semenov, D., Favre, C., Fedele, D., et al. 2018, A&A, 617, A28

Sipilä, O., Harju, J., & Caselli, P. 2017, A&A, 607, A26

Speagle, J. S. 2020, MNRAS, 493, 3132

Tang, Y.-W., Dutrey, A., Guilloteau, S., et al. 2016, ApJ, 820, 19

Taniguchi, K., Rayalacheruvu, P., Yonetsu, T., et al. 2024, ApJ, 963, 12

Teague, R. 2019, JOSS, 4, 1632

Teague, R. 2020, richteague/keplerian_mask: Initial Release v1.0, Zenodo, doi:10.5281/ZENODO.4321137

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

Teague, R., Semenov, D., Guilloteau, S., et al. 2015, A&A, 574, A137

Thi, W. F., van Zadelhoff, G. J., & van Dishoeck, E. F. 2004, A&A, 425, 955

Wagenblast, R., & Hartquist, T. W. 1989, MNRAS, 237, 1019

Wakelam, V., Chapillon, E., Dutrey, A., et al. 2019, MNRAS, 484, 1563

Wakelam, V., & Herbst, E. 2008, ApJ, 680, 371

Wakelam, V., Ruaud, M., Hersant, F., et al. 2016, A&A, 594, A35

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

Willacy, K. 2007, ApJ, 660, 441

Williams, J. P., & Best, W. M. J. 2014, ApJ, 788, 59

Wootten, A. 1987, in IAU Symp. 120, Astrochemistry, ed. M. Vardya & S. Tarafdar (Dordrecht: Reidel), 311