# Molecules with ALMA at Planet-forming Scales (MAPS). XIII. HCO+ and Disk Ionization Structure

Yuri Aikawa, Gianni Cataldi, Yoshihide Yamato, Ke Zhang, Alice S. Booth, Kenji Furuya, Sean M. Andrews, Jaehan Bae, Edwin A. Bergin, Jennifer B. Bergner, et al.

► To cite this version:

Yuri Aikawa, Gianni Cataldi, Yoshihide Yamato, Ke Zhang, Alice S. Booth, et al., Molecules with ALMA at Planet-forming Scales (MAPS). XIII. HCO+ and Disk Ionization Structure. The Astrophysical Journal Supplement Series, 2021, 257, 10.3847/1538-4365/ac143c. insu-03672367

HAL Id: insu-03672367

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

Submitted on 29 Jul 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Molecules with ALMA at Planet-forming Scales (MAPS) XIII: HCO$^+$ and disk ionization structure

YURI AIKAWA,<fnref n="1" /> GIANNI CATALDI,<fnref n="1" />,<fnref n="2" /> YOSHIHIDE YAMATO,<fnref n="1" /> KE ZHANG,<fnref n="3" />,<fnref n="4" />,* ALICE S. BOOTH,<fnref n="1" />,<fnref n="5" />,<fnref n="6" />

KENJI FURUYA,<fnref n="1" />,<fnref n="2" /> SEAN M. ANDREWS,<fnref n="1" />,<fnref n="7" /> JAEHAN BAE,<fnref n="1" />,<fnref n="8" />,<fnref n="9" />,† EDWIN A. BERGIN,<fnref n="1" />,<fnref n="4" /> JENNIFER B. BERGNER,<fnref n="1" />,<fnref n="10" />,†

ARTHUR D. BOSMAN,<fnref n="1" />,<fnref n="4" /> L. ILSEDORE CLEEVES,<fnref n="11" /> IAN CZEKALA,<fnref n="1" />,<fnref n="2" />,<fnref n="13" />,<fnref n="14" />,<fnref n="15" />,<fnref n="16" />,† VIVIANA V. GUZMÁN,<fnref n="1" />,<fnref n="17" />

JANE HUANG,<fnref n="1" />,<fnref n="4" />,<fnref n="7" />,† JOHN D. ILEE,<fnref n="1" />,<fnref n="6" /> CHARLES J. LAW,<fnref n="1" />,<fnref n="7" /> ROMANE LE GAL,<fnref n="1" />,<fnref n="7" />,<fnref n="18" />,<fnref n="19" />,<fnref n="20" /> RYAN A. LOOMIS,<fnref n="1" />,<fnref n="21" />

FRANÇOIS MÉNARD,<fnref n="1" />,<fnref n="19" /> HIDEKO NOMURA,<fnref n="1" />,<fnref n="2" /> KARIN I. ÖBERG,<fnref n="1" />,<fnref n="7" /> CHUNHUA QI,<fnref n="1" />,<fnref n="7" /> KAMBER R. SCHWARZ,<fnref n="1" />,<fnref n="22" />,†

RICHARD TEAGUE,<fnref n="1" />,<fnref n="7" /> TAKASHI TSUKAGOSHI,<fnref n="1" />,<fnref n="2" /> CATHERINE WALSH,<fnref n="1" />,<fnref n="6" /> AND DAVID J. WILNER,<fnref n="1" />,<fnref n="7" />

(Received March 3, 2021; Revised June 16, 2021; Accepted June 17, 2021)

Submitted to ApJ

ABSTRACT

We observed HCO$^+$ *J* = 1–0 and H$^{13}$CO$^+$ *J* = 1–0 emission towards the five protoplanetary disks around IM Lup, GM Aur, AS 209, HD 163296, and MWC 480 as part of the MAPS project. HCO$^+$ is detected and mapped at 0.3′′ resolution in all five disks, while H$^{13}$CO$^+$ is detected (SNR> 6σ) towards GM Aur and HD 163296 and tentatively detected (SNR> 3σ) towards the other disks by a matched filter analysis. Inside a radius of *R* ~ 100 au, the HCO$^+$ column density is flat or shows a central dip. At outer radii (≥ 100 au), the HCO$^+$ column density decreases outwards, while the column density ratio of HCO$^+/CO$ is mostly in the range of ~ 10$^{-5}$ – 10$^{-4}$. We derived the HCO$^+$ abundance in the warm CO-rich layer, where HCO$^+$ is expected to be the dominant molecular ion. At *R* ≥ 100 au, the HCO$^+$ abundance is ~ 3 × 10$^{-11}$ – 3 × 10$^{-10}$, which is consistent with a template disk model with X-ray ionization. At the smaller radii, the abundance decreases inwards, which indicates that the

<footnote n="1">¹Department of Astronomy, Graduate School of Science, The University of Tokyo, 113-0033, Tokyo, Japan</footnote>

<footnote n="2">²National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="3">³Department of Astronomy, University of Wisconsin-Madison, 475 N Charter St, Madison, WI 53706</footnote>

<footnote n="4">⁴Department of Astronomy, University of Michigan, 323 West Hall, 1085 S. University Avenue, Ann Arbor, MI 48109, USA</footnote>

<footnote n="5">⁵Leiden Observatory, Leiden University, 2300 RA Leiden, the Netherlands</footnote>

<footnote n="6">⁶School of Physics and Astronomy, University of Leeds, Leeds, LS2 9JT, UK</footnote>

<footnote n="7">⁷Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA</footnote>

<footnote n="8">⁸Earth and Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Road NW, Washington, DC 20015, USA</footnote>

<footnote n="9">⁹Department of Astronomy, University of Florida, Gainesville, FL 32611, USA</footnote>

<footnote n="10">¹⁰University of Chicago Department of the Geophysical Sciences, Chicago, IL 60637, USA</footnote>

<footnote n="11">¹¹Department of Astronomy, University of Virginia, Charlottesville, VA 22904, USA</footnote>

<footnote n="12">¹²Department of Astronomy and Astrophysics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="13">¹³Center for Exoplanets and Habitable Worlds, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="14">¹⁴Center for Astrostatistics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="15">¹⁵Institute for Computational & Data Sciences, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="16">¹⁶Department of Astronomy, 501 Campbell Hall, University of California, Berkeley, CA 94720-3411, USA</footnote>

<footnote n="17">¹⁷Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 7820436 Macul, Santiago, Chile</footnote>

<footnote n="18">¹⁸IRAP, Université de Toulouse, CNRS, CNES, UT3, 31400 Toulouse, France</footnote>

<footnote n="19">¹⁹Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="20">²⁰IRAM, 300 rue de la piscine, F-38406 Saint-Martin d'Hères, France</footnote>

<footnote n="21">²¹National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA 22903, USA</footnote>

<footnote n="22">²²Lunar and Planetary Laboratory, University of Arizona, 1629 E. University Blvd, Tucson, AZ 85721, USA</footnote>
ionization degree is lower in denser gas, especially inside the CO snow line, where the CO-rich layer is in the midplane. Comparison of template disk models with the column densities of HCO$^+$, N$_2$H$^+$, and N$_2$D$^+$ indicates that the midplane ionization rate is $\gtrsim 10^{-18}$ s$^{-1}$ for the disks around IM Lup, AS 209, and HD 163296. We also find hints of an increased HCO$^+$ abundance around the location of dust continuum gaps in AS 209, HD 163296, and MWC 480. This paper is part of the MAPS special issue of the Astrophysical Journal Supplement.

**Keywords:** Astrochemistry — Exoplanet formation — Interferometry — Millimeter astronomy — Protoplanetary disks

# 1. INTRODUCTION

Protoplanetary disks are the birth site of planetary systems including our Solar System. The rate and degree of ionization are important parameters for both the physical and chemical evolution of the disk, and thus for the formation of planetary systems. The disk gas is much denser, and thus its ionization degree $x_i$, which is the relative abundance of electrons to H$_2$, is much lower than in parental molecular clouds (e.g. Umebayashi & Nakano 1988)$^1$. Yet the ionization degree is expected to be sufficient for the gas to be partially coupled with the magnetic fields present in the disk, which could induce magneto-hydrodynamic instabilities, disk winds, and thus the angular momentum transfer needed for mass accretion (e.g. Suzuki & Inutsuka 2009; Bai & Stone 2013). More specifically, the ionization degree, and thus the coupling between gas and magnetic fields (i.e. non-ideal MHD effects) should vary spatially within the disk, which affects the physical structure of the disk (e.g. Béthune et al. 2017). Ionization also triggers ion-molecule reactions in the disk. For example, carbon monoxide, which is the main tracer of disk gas, can be converted to other species by ion-molecule reactions within the typical evolutionary timescale of the disk (10$^6$ yr), if the ionization rate $\zeta$ is similar to or higher than the cosmic-ray ionization rate in the ISM (10$^{-17}$ s$^{-1}$; Furuya & Aikawa 2014; Bergin et al. 2014; Bosman et al. 2018; Schwarz et al. 2018).

The ionization degree is determined by the ionization rate, gas density, and number density and size distribution of dust grains. For protoplanetary disks, there are several possible ionization sources: X-rays from the central star, cosmic rays (CRs) and stellar energetic particles, and the decay of short-lived radio active nuclei (SLR) (e.g. Umebayashi & Nakano 1981; Glassgold et al. 1997; Umebayashi & Nakano 2009; Rab et al. 2017). X-rays are the dominant ionization source in the disk surface. However, X-rays are significantly attenuated before reaching the disk midplane. The attenuation length corresponds to a hydrogen column density of $N_H \sim 10^{22}$ cm$^{-2}$ at 1 keV, while the hydrogen column densities of our target disks are $\sim 10^{23} - 10^{25}$ cm$^{-2}$ at the radius of 100 au (Zhang et al. 2021). Although the attenuation length is larger at higher energies, the X-ray ionization rate is still expected to be $\lesssim 10^{-18}$ s$^{-1}$ in the midplane (e.g. §4) (see also Rab et al. 2018). The cosmic-ray ionization rate is $\sim 5 \times 10^{-17}$ s$^{-1}$ in molecular clouds (Dalgarno 2006), and its attenuation length is much larger than that of X-rays. While Umebayashi & Nakano (1981) evaluated the attenuation length to be 96 g cm$^{-2}$, more recent work by Padovani et al. (2018) reassessed the propagation of CR particles and obtained an even larger attenuation length. But CRs can be scattered by magnetized stellar winds and/or magnetic fields in the disk (e.g. Umebayashi & Nakano 1981; Cleeves et al. 2014; Padovani et al. 2018). Stellar energetic particles could also be disturbed by magnetic fields. Finally, the ionization rate from the decay of $^{26}$Al is estimated to be $10^{-18}$ s$^{-1}$ in the primordial Solar System based on meteorite analysis (e.g. Umebayashi & Nakano 2009). While the abundance of SLRs should vary among star-forming regions, recent chemo-hydrodynamic simulations of the Milky Way Galaxy predict that the abundance of $^{26}$Al in the primordial Solar System may be typical for other star-forming regions (Fujimoto et al. 2018). In summary, the dominant ionization source and ionization degree should vary spatially within the disk and could be different between disks.

The importance of understanding ionization processes motivates the observation of molecular ions in protoplanetary disks. Theoretical models show that the major molecular ions in disks are H$_3^+$, HCO$^+$, N$_2$H$^+$, and their deuterated isotopologues (e.g. Aikawa & Herbst 2001; Bergin et al. 2007; Willacy 2007; Aikawa et al. 2015, 2018). The most abundant molecular ion varies both radially and vertically in disks (Figure 1). The disk surface is the Photon Dominated Region (PDR); photoionization makes atomic ions such as C$^+$ and S$^+$ the main charge carriers. Deeper in the disk,

<footnote>* NASA Hubble Fellow</footnote>

<footnote>† NASA Hubble Fellowship Program Sagan Fellow</footnote>

<footnote n="1">¹ Considering neutrality, the electron abundance should be the same as the total abundance of cations. In dense midplane regions, grain particles can be the dominant charge carrier.</footnote>
$H_3^+$ becomes the most abundant ion, while $HCO^+$ becomes dominant when the abundance ratio $n(CO)/n(e)$ is higher than $\sim 10^3$, where $n(i)$ denotes the number density of species $i$ (Aikawa et al. 2015). In the regions where $T < 20$ K (i.e. below the CO snow surface and outside the CO snow line), CO freezes out onto grains, which enhances the $N_2H^+$ abundance. As temperature declines towards deeper layers, the dominant ion changes from $N_2H^+$ to $H_3^+$ and its deuterated isotopologues. Unfortunately, $H_3^+$ cannot be observed at millimeter wavelengths, and its deuterated isotopologues have not been detected in disks so far. $N_2D^+$ is considered to be an alternative probe to constrain the ionization degree in the cold midplane (Cleeves et al. 2014).

While the rotational transitions of $HCO^+$, $N_2H^+$, and their isotopologues have been observed towards several disks, a quantitative evaluation of their column densities and the disk ionization degree has been reported only in a limited number of references. Öberg et al. (2011a) used the IRAM 30 m telescope to observe $H^{13}CO^+ J = 3-2$ in the disk of DM Tau. They combined their observation with previous Submillimeter Array (SMA) data of $N_2H^+ J = 3-2$, $HCO^+ J = 3-2$, and $DCO^+ J = 3-2$ to estimate the ionization degree in three temperature regions of the disk: in the upper warm molecular layer ($T > 20$ K), the ionization degree $x_i$ is estimated to be $4 \times 10^{-10}$ based on the $HCO^+$ data. In the cooler molecular layer ($T = 16-20$ K), where $N_2H^+$ and $DCO^+$ would be abundant, $x_i$ is derived to be $3 \times 10^{-11}$, while in the cold, dense midplane ($T < 16$ K), the non-detection of $H_2D^+$ constrains $x_i$ to be $< 3 \times 10^{-10}$. Teague et al. (2015) observed $HCO^+ J = 1-0$ and $J = 3-2$, and $DCO^+ J = 3-2$ in DM Tau at $\sim 1.5''$ resolution using the Plateau de Bure Interferometer (PdBI). The column densities of $HCO^+$ and $DCO^+$ at a radius of 100 au are derived to be $9.8 \times 10^{12}$ and $1.2 \times 10^{12}$ cm$^{-2}$, respectively. They derive an ionization degree of $\sim 10^{-7}$ from the abundance ratio of $DCO^+/HCO^+$, assuming steady state balance between $H_3^+$ + HD $\rightarrow$ $H_2D^+$ + $H_2$, $H_2D^+$ + CO $\rightarrow$ $DCO^+$ + $H_2$, and the destruction of $DCO^+$ (Caselli et al. 2002). Cleeves et al. (2015) calculated chemical models of the disk around TW Hya for various cosmic-ray ionization rates $\zeta_{CR}$ and X-ray spectra. Instead of evaluating molecular column densities from the observational data, they calculated the disk-integrated flux of molecular lines from the models to compare with their own $HCO^+$ and $H^{13}CO^+$ $J = 3-2$ observations of TW Hya, and $HCO^+$, $H^{13}CO^+$, and $N_2H^+$ data from the literature (e.g. Qi et al. 2013a,b). They concluded that the observations are best reproduced by the model with low CR ionization rate $\zeta_{CR} \lesssim 10^{-19}$ s$^{-1}$ and modest X-ray spectra, in which the ionization degree is $\sim 10^{-11} - 10^{-10}$ near the midplane outside of a radius of $\sim 100$ au. Since the ionization degree varies spatially and between disks, observations with higher spatial resolution and towards more targets are highly desirable.

In light of the ring-gap structures observed in the millimeter dust continuum of many disks in recent years (e.g. Andrews et al. 2018; Huang et al. 2018; Cieza et al. 2021), it is also interesting to investigate whether and how the molecular ion abundances and ionization degree correlate with the dust substructures. Theoretical studies predict that $HCO^+$ is sensitive to, and thus could be a good probe of, gas density and/or small dust abundance. In the region where $HCO^+$ is the most abundant positive charge carrier, its abundance (i.e. ionization degree) should be proportional to $n_H^{-1/2}$ (e.g. Duley & Williams 1984). In the gap region, the ionization degree could also be enhanced due to deeper penetration of X-rays. Aikawa & Nomura (2006), on the other hand, showed that the $HCO^+$ column density declines as the abundance of sub-micron sized dust grains decreases with grain growth. This is because the PDR, in which atomic ions dominate over molecular ions, extends deeper into the disk when the abundance of sub-micron dust grains decreases. Similar results are obtained in models of disks around Herbig Ae stars (Jonkheid et al. 2007) (see also Wakelam et al. 2019). While these models do not include any ring-gap structures, they indicate that the $HCO^+$ abundance could vary over the substructures. The $HCO^+$ abundance would be higher in the gap than in the ring, if the gas density is lower but the UV attenuation is sufficient in the gap. On the other hand, the $HCO^+$ abundance would be lower in the gap if the UV penetration is significant (see also Smirnov-Pinchukov et al. 2020).

In the present work, we analyze the $J = 1-0$ transition of $HCO^+$ and $H^{13}CO^+$ towards five disks around IM Lup, GM Aur, AS 209, HD 163296, and MWC 480 observed as a part of the Molecules with ALMA at Planet-forming Scales (MAPS) Large Program (Oberg et al. 2021, project code 2018.1.01055.L). For a given column density of $HCO^+$ at a lukewarm temperature (a few tens of K), which corresponds to the CO rich molecular layer in the disk, the $J = 1-0$ line tends to be optically thinner than higher transitions, and thus useful to derive the column density. Our spatial resolution, 0.3'', is high enough to investigate the radial distribution of $HCO^+$ and to marginally resolve the continuum gap regions. We also combine our results with the analysis of CO isotopologues, $N_2H^+$, and $N_2D^+$ in the companion MAPS papers by Zhang et al. (2021) and Cataldi et al. (2021) to estimate the radial variation in the ionization degree and the ionization rate in the disks. Various data products of the MAPS project, including the present work, such as reduced observational data, zeroth-moment maps, as well as derived $HCO^+$ column densities can be downloaded at www.alma-maps.info.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematic distribution of major ions in protoplanetary disks.</figcaption></figure>

<table><caption>**Table 1.** Molecular data of the targeted lines taken from the CDMS database.</caption><thead><tr><th>transition</th><th>frequency<br>[GHz]</th><th>A<sub>ij</sub><sup>a</sup><br>[s<sup>-1</sup>]</th><th>E<sub>u</sub><sup>b</sup><br>[K]</th></tr></thead><tbody><tr><td>HCO<sup>+</sup> J = 1–0</td><td>89.1885247</td><td>4.2512 × 10<sup>-5</sup></td><td>4.28</td></tr><tr><td>H<sup>13</sup>CO<sup>+</sup> J = 1–0</td><td>86.7542884</td><td>3.8534 × 10<sup>-5</sup></td><td>4.16</td></tr></tbody></table>

<sup>a</sup>Einstein A coefficient.

<sup>b</sup>Upper state energy.

The outline of this paper is as follows. We briefly describe the observations in §2. §3 describes the observational results: the zeroth moment maps, radial emission profiles, azimuthally averaged spectra, and molecular column densities. We derive the HCO⁺ abundance, which corresponds to a lower limit of the ionization degree in the warm molecular layer, and compare our results with a template disk chemistry model in §4. Our conclusions are presented in §5.

## 2. OBSERVATIONS

The observational setup chosen by MAPS covered four spectral setups: two in Band 3 at ~3 mm and two in Band 6 at ~1 mm. The *J* = 1–0 transitions of HCO⁺ and H<sup>13</sup>CO⁺ are in one of the Band 3 setups. The spectral resolution is 0.237 km s<sup>-1</sup> for HCO⁺ and 0.488 km s<sup>-1</sup> for H<sup>13</sup>CO⁺. Molecular data for the targeted lines is taken from the CDMS database<fnref n="2" /> (Müller et al. 2001; Endres et al. 2016) and summarized in table 1.

Two array configurations were used: a short baseline configuration that was sensitive to extended structure, and a long baseline configuration to achieve high angular resolution. The details of the data calibration are described in Oberg et al. (2021). The MAPS collaboration produced images from the calibrated visibilities using the CLEAN algorithm implemented in the CASA tclean task. Details of the imaging procedure are described in Czekala et al. (2021). We use the fiducial images provided by MAPS, which have a 0.′3 circular beam. This beam size corresponds to a physical size of 30.3 au at the distance of HD 163296, which is the closest target, and 48.6 au at the distance of MWC 480, which is the most distant. Basic properties of the data cubes such as noise levels are summarized in Table 11 and 12 in Oberg et al. (2021).

## 3. OBSERVATIONAL RESULTS

### 3.1. Zeroth moment maps, disk-integrated fluxes, and radial emission profiles

The MAPS collaboration produced zeroth moment maps by applying a Keplerian mask to the data cubes and integrating over the velocity axis (Law et al. 2021a). The mask parameters are conservatively chosen to incorporate all <sup>13</sup>CO *J* = 2–1 emission, which is the most widespread emission of any species except for CO *J* = 2–1. The mask parameters for the five disk are summarized in Table 1 in Czekala et al. (2021). The stellar masses adopted for the

<footnote n="2">² https://cdms.astro.uni-koeln.de/</footnote>
generation of the masks are listed in Table 1 in Oberg et al. (2021). These maps were used for all scientific analysis, in particular to derive radial emission profiles (Law et al. 2021a). A Keplerian mask, however, introduces spatial variance in noise distribution, since the number of channels summed to create the zeroth moment maps vary spatially. The discontinuous noise distribution occasionally causes arc-like artefacts, which could be mistaken as real substructure, in the central regions. To mitigate such artefacts, MAPS also produced “hybrid” zeroth moment maps by combining a Keplerian mask and a smoothed $\sigma$-clip mask (thus their name “hybrid”, see Law et al. 2021a, Appendix A). The latter removes the pixels with a S/N below a given threshold. We emphasize that by using a clipping mask with a threshold larger than $0\sigma$, some emission is inevitably lost. Therefore, the hybrid zeroth moment maps are for presentational purposes only, and are not used for any quantitative analysis.

In Figure 2, we show hybrid zeroth moment maps for $HCO^+ J = 1-0$ and $H^{13}CO^+ J = 1-0$. We choose the best-looking $\sigma$-clip value by visual inspection: $1\sigma$ for $HCO^+ J = 1-0$ towards IM Lup, and $0\sigma$ for the other maps. Hybrid zeroth moment maps with a $0''5$ tapered beam are shown in Figure 13 in Appendix A; the faint emission is more easily seen in the lower spatial resolution images.

The $HCO^+ J = 1-0$ line is clearly detected in all five disks. The emission is relatively bright inside a radius of ~200 au, while its diffuse emission extends to $\gtrsim 400$ au, except for the disk of AS 209, in which the emission extends only up to $R \sim 300$ au. In AS 209, the emission from the west side of the disk is absorbed by a foreground cloud, which happens to have a similar line-of-sight velocity (Oberg et al. 2011b). We thus performed all analysis for AS 209 only on the un-obscured east side of the disk within a $\pm 55$ degree wedge centered on the semi-major axis (Teague et al. 2018b; Law et al. 2021a). The emission of $H^{13}CO^+ J = 1-0$ is very weak and is not clearly seen in the zeroth moment maps.

Disk-integrated fluxes of $HCO^+ J = 1-0$ and $H^{13}CO^+ J = 1-0$ are presented in Table 2. These were calculated by integrating the flux within a Keplerian mask. The radial extent of the Keplerian mask is different from that used for the zeroth moment maps. For $HCO^+$, the radial extent of the mask is visually determined from the extent of the emission seen in the radial emission profile (see Figure 3 and Table 2). For $H^{13}CO^+$, we measured fluxes by using the same radial extents as for $HCO^+$. In an attempt to increase the SNR, we computed additional fluxes by adopting the same radial extent as the Keplerian models that maximize the SNR in the matched filter analysis in the uv plane (Appendix B, Figure 14), as shown in the right column of Table 2. Errors were calculated by repeating the flux measurement procedure at off-source positions and taking the standard deviation of the off-source fluxes. A 10% flux calibration error was added in quadrature. If the resulting SNR was smaller than 3, Table 2 reports the $3\sigma$ upper limit. To calculate these disk-integrated fluxes and errors, we used the image cubes that are not corrected for the primary beam. This is because for primary beam corrected images, the noise increases towards the edges of the image, making our approach to estimate the error from off-source positions invalid. However, we verified that the difference to fluxes extracted from primary beam corrected images is negligible. For deriving radially resolved column density profiles, we use the primary beam corrected images.

$H^{13}CO^+$ is seen with SNR> 3 for IM Lup, GM Aur, and HD 163296, if we use the masks with the extent informed by the matched filter analysis. The matched filter analysis itself (Appendix B) gives a response larger than $6\sigma$ for GM Aur and HD 163296, and larger than $3\sigma$ for the other disks. Thus, we consider $H^{13}CO^+ J = 1-0$ detected in GM Aur and HD 163296, and tentatively detected for the other three disks. We note that the values in Table 2 are different from those in Table 12 in Oberg et al. (2021), who adopted the same Keplerian masks to generate zeroth moment maps and to estimate disk-integrated fluxes.

Figure 3 shows the radial emission profiles of the $J = 1-0$ transitions of $HCO^+$ and $H^{13}CO^+$. These profiles are produced by azimuthal averaging the zeroth moment maps that were produced with a Keplerian mask only (i.e. without a $\sigma$-clip mask). The inclination and position angles of each disk are listed in Table 1 in Oberg et al. (2021). Since $HCO^+$ is expected to be formed by $CO + H_3^+$ and to be a dominant molecular ion in CO gas rich layers, we also plotted the radial emission profile of $C^{18}O J = 1-0$ derived by Law et al. (2021a). The vertical lines indicate the radius of rings (solid), gaps (dashed), and the dust disk edge (dotted) observed in the millimeter dust continuum (Huang et al. 2018; Long et al. 2018; Liu et al. 2019; Sierra et al. 2021; Law et al. 2021a). Both $HCO^+$ and $C^{18}O$ emission extend out to several 100 au. In the IM Lup disk, the radial emission profile of $HCO^+ J = 1-0$ is similar to that of $C^{18}O$, showing a ring-like distribution, although the peak position of the $HCO^+$ emission is ~36 au outside of the $C^{18}O$ peak. In GM Aur, on the other hand, both the $C^{18}O$ and $HCO^+$ emission are centrally peaked. The radial emission profiles of $C^{18}O$ and $HCO^+$ are quite different in AS 209: around a radius of ~77 au, the $HCO^+$ emission has a local maximum, while the radial gradient of $C^{18}O$ changes from negative to plateau. Interestingly, the peak
<table><caption>**Table 2.** Disk-integrated fluxes. Upper limits are at $3\sigma$ significance.</caption><thead><tr><th rowspan="3"></th><th colspan="3">HCO<sup>+</sup> 1-0</th><th colspan="3">H<sup>13</sup>CO<sup>+</sup> 1-0<sup>a</sup></th><th colspan="3">H<sup>13</sup>CO<sup>+</sup> 1-0<sup>b</sup></th></tr><tr><th>r<sub>min</sub> <sup>c</sup></th><th>r<sub>max</sub> <sup>d</sup></th><th>flux</th><th>r<sub>min</sub></th><th>r<sub>max</sub></th><th>flux</th><th>r<sub>min</sub></th><th>r<sub>max</sub></th><th>flux</th></tr><tr><th>[au]</th><th>[au]</th><th>[mJy km s<sup>-1</sup>]</th><th>[au]</th><th>[au]</th><th>[mJy km s<sup>-1</sup>]</th><th>[au]</th><th>[au]</th><th>[mJy km s<sup>-1</sup>]</th></tr></thead><tbody><tr><td>IM Lup</td><td>0</td><td>700</td><td>519 ± 53</td><td>0</td><td>700</td><td>&lt; 46</td><td>0</td><td>400</td><td>23 ± 6</td></tr><tr><td>GM Aur</td><td>0</td><td>500</td><td>412 ± 42</td><td>0</td><td>500</td><td>&lt; 43</td><td>0</td><td>400</td><td>27 ± 7</td></tr><tr><td>AS 209</td><td>0</td><td>300</td><td>178 ± 19</td><td>0</td><td>300</td><td>&lt; 33</td><td>0</td><td>200</td><td>&lt; 26</td></tr><tr><td>HD 163296</td><td>0</td><td>550</td><td>1002 ± 101</td><td>0</td><td>550</td><td>&lt; 54</td><td>50</td><td>400</td><td>27 ± 7</td></tr><tr><td>MWC 480</td><td>0</td><td>550</td><td>323 ± 34</td><td>0</td><td>550</td><td>&lt; 34</td><td>50</td><td>100</td><td>&lt; 18</td></tr></tbody></table>

<sup>a</sup>For a Keplerian mask with the same radial extent as used for HCO<sup>+</sup> 1-0.

<sup>b</sup>For a Keplerian mask with the same radial extent as the matched filter that maximises the SNR of the filter response (see Appendix B and Fig. 14).

<sup>c</sup>Minimum radius of the Keplerian mask.

<sup>d</sup>Maximum radius of the Keplerian mask.

position of the HCO<sup>+</sup> emission coincides with that of a dust continuum ring. In the disks around the Herbig Ae stars, HD 163296 and MWC 480, the HCO<sup>+</sup> emission shows a central dip, while C<sup>18</sup>O is flat for HD 163296 and centrally peaked for MWC 480.

The intensity ratio of the C<sup>18</sup>O and HCO<sup>+</sup> lines varies significantly with radius and among objects. For example, the ratio reaches ~5 towards the center of the MWC 480 disk, while HCO<sup>+</sup> is brighter than C<sup>18</sup>O around R ~ 77 au in AS 209. Although this indicates that the abundance (column density) ratio of HCO<sup>+</sup> to C<sup>18</sup>O also varies, we need to consider the line optical depth. Under LTE conditions with an excitation temperature of 30 K, which is a typical temperature in the warm molecular layers, the optical depth of C<sup>18</sup>O J = 1-0 reaches unity for a C<sup>16</sup>O column density of ~6 × 10<sup>18</sup> cm<sup>-2</sup>, assuming a <sup>16</sup>O/<sup>18</sup>O isotope ratio of 557 (Wilson 1999). Zhang et al. (2021) derived the radial column density distribution of CO by analyzing several transitions of CO and its isotopologues. The CO column density exceeds 6 × 10<sup>18</sup> cm<sup>-2</sup> inside a radius of ~100 au in the GM Aur, HD 163296, and MWC 480 disks. At least for these inner radii, C<sup>18</sup>O J = 1-0 is optically thick, and thus the intensity ratio of C<sup>18</sup>O to HCO<sup>+</sup> does not reflect their column density ratio. We also need to evaluate the optical depth of HCO<sup>+</sup> to derive its column density; the optical depth of HCO<sup>+</sup> J = 1-0 reaches unity for a HCO<sup>+</sup> column density of ~1 × 10<sup>13</sup> cm<sup>-2</sup> under LTE conditions at a temperature of 30 K.

## 3.2. Azimuthally averaged spectra and radial profiles of the HCO<sup>+</sup> column density

In order to evaluate the radial profile of the optical depth and column density of HCO<sup>+</sup>, we first derive the azimuthally averaged spectrum for radial bins with a width of half the beam size. Due to the Keplerian rotation of the disk, at each spatial location in the data cube, the spectrum is shifted with respect to the systemic velocity. We thus shift each spectrum by the Keplerian velocity projected to the line of sight before averaging the spectra azimuthally. This results in all individual spectra being centered at the systemic velocity, which increases the S/N ratio of the azimuthally averaged spectrum (Yen et al. 2016; Teague et al. 2016; Matrà et al. 2017). The calculation of the error bars of the averaged spectrum is described in detail in Appendix D of Cataldi et al. (2021). Briefly, we calculate both the standard deviation of the spectrum in regions without line emission, and an analytical error bar based on the number of independent samples included in the average. We then adopt the maximum of the two as our final error bar. Figure 4 shows examples of extracted spectra of HD 163296. The full gallery of spectra is found in Appendix C.

We then calculate the column density of HCO<sup>+</sup>, N(HCO<sup>+</sup>), at each radial bin by simultaneously fitting the azimuthally averaged spectra of HCO<sup>+</sup> and H<sup>13</sup>CO<sup>+</sup>. We employ the same fitting procedure as Cataldi et al. (2021). Briefly, for a given HCO<sup>+</sup> column density, we compute HCO<sup>+</sup> J = 1-0 and H<sup>13</sup>CO<sup>+</sup> J = 1-0 model spectra that can be compared to the data. There are five free parameters: the logarithm of the HCO<sup>+</sup> column density, the excitation temperature, an offset Δv of the line center with respect to the systemic velocity for each line, and the full width at half maximum (FWHM) of the Gaussian kernel with which the model spectra are convolved to mimic observational broadening. The H<sup>13</sup>CO<sup>+</sup> column density is fixed to 1/68 of the HCO<sup>+</sup> column density (Milam et al. 2005). We use an MCMC method implemented in the emcee package (Foreman-Mackey et al. 2013) to explore the parameter
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Gallery of hybrid zeroth moment maps for HCO⁺ *J* = 1–0 (top row) and H¹³CO⁺ *J* = 1–0 (bottom row) for the MAPS sample, ordered from left to right by increasing stellar mass (see Oberg et al. 2021, table 1). These maps were generated by combining a Keplerian mask with a smoothed σ-clip mask (Appendix A of Law et al. 2021a). The color scales employ either linear or arcsinh stretches, with the lower end saturating at 0 mJy beam⁻¹ km s⁻¹. The synthesized beam is shown by the white ellipse in the lower left of each panel. Due to the use of a mask, the noise level is not constant across a map. Lines that are only tentatively detected in total flux or matched filter are marked in the lower right of the panel.</figcaption></figure>

space and assume flat priors. We employ 200 walkers taking 5000 steps each, and discard the first 2500 steps from the analysis. The free parameters and the prior boundaries are listed in Table 3. The full details of the fitting procedure can be found in Cataldi et al. (2021). Figure 4 shows examples of model spectra fitted to the HD 163296 data. The full gallery of models is shown in Appendix C.

Figure 5 and 6 show the optical depths and excitation temperature derived by fitting the azimuthally averaged spectra. We find that HCO⁺ *J* = 1–0 is optically thick at the column density peaks. However, we are still able to constrain the HCO⁺ column density because the H¹³CO⁺ *J* = 1–0 transition is in the optically thin regime. The top panels in Figure 7 present the derived HCO⁺ column densities.

The excitation temperature is poorly constrained for large radial regions of the disks, as shown in Fig. 6. However, at least a lower limit on the excitation temperature can be placed in some disk locations. The innermost ~50 au of GM Aur stand out with an apparently well constrained *T*<sub>ex</sub> ≈ 90 K. This region also requires a large FWHM of 17 km s⁻¹. These rather extreme parameters might indicate that the fit is not reliable. In fact, considerable broadening is seen for all disks for the two innermost radial bins (i.e. within one beam FWHM from the disk center) due to beam smearing of the large velocity gradient in the inner disk (see Figures 16a). Therefore, the column density estimates for the two innermost radial bins should be considered with caution.

In Appendix D, we show optical depth and column density profiles for additional fits where the excitation temperature has been fixed to 30 K, which is a typical temperature for the warm molecular layer where HCO⁺ is expected to be present. The corresponding models are shown in Figure 16a. The column densities generally agree well with the fits where *T*<sub>ex</sub> is a free parameter, as can be seen in Fig. 17. Furthermore, for both fits, the model spectra fit the data well (Figure 16a). One exception is the region inwards of ~150 au towards GM Aur where the models with *T*<sub>ex</sub> as a free parameter generally provide a better fit, especially for the two innermost radial bins discussed in the previous
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Radial emission profiles of HCO$^+$ *J* = 1–0, H$^{13}$CO$^+$ *J* = 1–0, and C$^{18}$O *J* = 1–0 for IM Lup, GM Aur, AS 209, HD 163296, and MWC 480. The error bar shows the ±1σ error. The size of the radial bins is 0.075'', which is a quarter of the beam size. The vertical lines indicate the positions of rings (solid), gaps (dashed), and the edge of millimeter dust continuum (dotted) referring to Huang et al. (2018); Long et al. (2018); Liu et al. (2019); Law et al. (2021a) (see also Sierra et al. 2021).</figcaption></figure>

<table><caption>**Table 3.** Free parameters for fitting of azimuthally averaged spectra to derive the HCO$^+$ column density.</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 N<sub>HCO+</sub></td><td>1</td><td>16</td><td>log([cm<sup>-2</sup>])</td></tr><tr><td>T<sub>ex</sub></td><td>10</td><td>100</td><td>[K]</td></tr><tr><td>FWHM<sub>kernel</sub></td><td>0.5<sup>c</sup></td><td>variable<sup>d</sup></td><td>[km s<sup>-1</sup>]</td></tr><tr><td>Δv<sub>HCO+</sub><sup>e</sup></td><td>-0.3</td><td>0.3</td><td>[km s<sup>-1</sup>]</td></tr><tr><td>Δv<sub>H<sup>13</sup>CO+</sub></td><td>-0.3</td><td>0.3</td><td>[km s<sup>-1</sup>]</td></tr></tbody></table>

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

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

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

<sup>d</sup>Initial value for innermost radial bin is 20 km s<sup>-1</sup>. Dynamically adjusted when sequentially fitting larger and larger radii (see Cataldi et al. 2021, for details).

<sup>e</sup>Offset of the line center with respect to the systemic velocity.

paragraph (Fig. 16a). The other exception is the region inwards of 76 au towards HD 163296, where the fits with T<sub>ex</sub> = 30 K underpredict the HCO<sup>+</sup> emission and overpredict the H<sup>13</sup>CO<sup>+</sup> emission (Fig. 16a). Here, the models with T<sub>ex</sub> as a free parameter provide a better fit and predict a column density smaller by a factor of a few compared to the fits where T<sub>ex</sub> = 30 K.

In Figure 7, we also plot the column densities of CO which are derived from C$^{18}$O *J* = 2–1 emission in Zhang et al. (2021). For each disk, they constructed a thermo-chemical model which reproduces the SED. The model provides the 2D (*R*, *Z*) distributions of temperature and CO abundance. Then they introduced a CO depletion factor, which is varied from 0.001 to 50 with a step-size of 1.1 in a logarithmic scale, at each radial bin, and calculated the CO flux for grid of models to find the best fit depletion factor. The uncertainty in the CO column density is thus 10 % in most of the disk radii. Exceptions are the innermost and outermost radii; at these radii, the C$^{18}$O line fluxes have relatively
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Examples of HCO$^+$ and H$^{13}$CO$^+$ *J* = 1–0 model spectra fit to azimuthally averaged spectra of HD 163296 for a few radial bins. The orange and blue solid lines show the data, with the shaded regions corresponding to the 1σ uncertainty. The black curves show 50 randomly selected models drawn from the Monte Carlo chain, with the selection probability being proportional to the posterior probability of the model. Spectra are vertically offset for clarity.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Optical depth profiles of HCO$^+$ (*J* = 1–0) and H$^{13}$CO$^+$ (*J* = 1–0) derived from fitting azimuthally averaged spectra. The black lines show the median of the posterior probability, while the shaded regions extend between the 16th to 84th and the 2.3th to 97.7th percentiles. The beam size is shown as a horizontal line in the upper right.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** The excitation temperature derived from fitting azimuthally averaged spectra. The black lines show the median of the posterior probability, while the shaded regions extend between the 16th to 84th and the 2.3th to 97.7th percentiles. The beam size is shown as a horizontal line in the upper right.</figcaption></figure>

large uncertainties, and the error of the CO column density is 1 sigma. They also derived the CO column densities using the *J* = 2 − 1 transition of ¹³CO, and *J* = 1 − 0 transition of ¹³CO, C¹⁸O, and C¹⁷O, which are in reasonable agreement with the value based on C¹⁸O (*J* = 2 − 1).

Overall, the radial profile of the HCO⁺ column density is similar to that of CO. The panels in the middle row of Figure 7 show the column density ratio *N*(HCO⁺)/*N*(CO). Since the CO column density was calculated on a finer radial grid than HCO⁺, we interpolated the CO column density at the radii of HCO⁺. The error bars in the middle panels reflect only the error of HCO⁺ column density, which dominates over that of CO column density. The ratio is mostly within 10⁻⁵ − 10⁻⁴ at ∼ 100 − 500 au, which suggests a close chemical link between these species, as expected from chemical models (e.g. Aikawa et al. 2015; Teague et al. 2015).

In the disk of IM Lup, the column densities of both HCO⁺ and CO increase inwards from *R* ~ 500 au to 100 au, where the HCO⁺/CO column density ratio is remarkably constant. The HCO⁺ column density then declines at ~ 100 au, while the CO column density becomes flat and rises inwards at ~ 50 au. In the GM Aur disk, the HCO⁺/CO column density ratio exceeds 10⁻⁴ outside the dust continuum edge. Both the HCO⁺ and CO column densities increase inwards up to ~ 100 au, while they have small local humps. At *R* ~ 60 au, the CO column density has a local maximum, while that of HCO⁺ has a local minimum. This coincidence needs to be taken with caution, however, since the beam size of HCO⁺ data is two times larger than that of CO, and since the column density estimates at the innermost radii suffers line broadening due to the velocity gradient within a beam. In the disk of AS 209, the radial profile of HCO⁺ significantly differs from that of CO at *R* ~ 100 au. While the CO column density shows a broad depression at *R* ~ 45 – 120 au, the HCO⁺ column density increases inwards up to ~ 75 au, and stays constant at inner radii within the error bars. In the disk around HD 163296, the column density ratio of HCO⁺/CO decreases inwards up to *R* ~ 240 au, which corresponds to the edge of the millimeter dust continuum. From 240 au to 65 au, the ratio stays constant or has a shallow rise, while both the CO and HCO⁺ column densities increase inwards. Inside the radius of ~ 65 au, which coincides with the CO snow line in the thermo-chemical model, the HCO⁺ column density is roughly constant, while the CO column density increases towards the center by more than an order of magnitude. The sharp rise of CO column density and relatively flat distribution of HCO⁺ column density inside the CO snow line is also seen in MWC 480, while the overall column density ratio of HCO⁺/CO is slightly lower than that in HD 163296. We note that the X-ray spectrum of MWC 480 is significantly softer than that of HD 163296 (Dionatos et al. 2019). Since X-rays are the major ionization source in the molecular layer (see §4.3), the relatively low HCO⁺/CO column density ratio in MWC 480 could be due to the lack of high energy X-rays (≥ a few keV).

# 4. DISCUSSION

## 4.1. HCO⁺ abundance in the warm molecular layers

In the molecular layers of protoplanetary disks, the major molecular ions are H<sub>3</sub><sup>+</sup>, HCO<sup>+</sup> and N<sub>2</sub>H<sup>+</sup>, among which H<sub>3</sub><sup>+</sup> cannot be observed at millimeter wavelengths. Its deuterated counterpart o-H<sub>2</sub>D<sup>+</sup> has not been detected so far (e.g. Chapillon et al. 2011). Figure 8 shows the radial profiles of the HCO<sup>+</sup>, N<sub>2</sub>H<sup>+</sup>, and N<sub>2</sub>D<sup>+</sup> column densities obtained in this work and by Cataldi et al. (2021). We can see that HCO<sup>+</sup> has the largest column densities among the observable molecular ions in our five disks.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** (Top) Column densities of HCO⁺ (red) and CO column density scaled by 10⁻⁵ (black lines with gray error bars). The error bars correspond to 10% or 1 σ for CO (see §3.2), while they show the 16th and 84th percentile for HCO⁺. At radii where the median value of the molecular column density is lower than the value at 16th percentile by a factor of > 10, we plot the upper limits as inverted triangles (84 percentile for closed triangles and 98 percentile for open triangles). The vertical orange bars mark the CO snowline as inferred from the model by Zhang et al. (2021) with error of ±5 au. The blue solid and dashed lines depict the radius of rings and gaps, respectively, in dust continuum. The black dotted lines depict the outer edge of the millimeter dust continuum. (Middle) The column density ratio of HCO⁺ to CO. The yellow bars mark the HCO⁺/CO column density ratio of 10⁻⁵ – 10⁻⁴. (Bottom) Abundance of HCO⁺ relative to hydrogen nuclei in the CO gas rich layer. The yellow bars mark the abundance of 3 × 10⁻¹¹ – 3 × 10⁻¹⁰. In the middle and bottom panels, we consider only the error in HCO⁺ column density, since it dominates over the error of CO column density estimates. At radii where the median value of the HCO⁺ column density is lower than the value at 84th percentile by a factor of > 10, we plot the value which corresponds to the 1-σ upper limit (i.e. 84 percentile of HCO⁺ column density) with an inverted triangle.</figcaption></figure>

HCO⁺ is expected to be the major molecular ion in the warm molecular layer where CO gas is abundant. Thus, here we assume that the majority of HCO⁺ coexists with CO in the warm molecular layer, and derive the HCO⁺ abundance there, based on the HCO⁺ to CO column density ratio. Furthermore, we assume that the CO abundance relative to hydrogen nuclei is equal to its canonical abundance in molecular clouds, 5.0 × 10⁻⁵, multiplied by the CO depletion factor derived in Zhang et al. (2021). The bottom panels in Figure 7 show the estimated HCO⁺ abundance. Since the uncertainty (error) of the HCO⁺ column density dominates over the that of the CO column density at most disk radii, and since the methods of error estimation are different between HCO⁺ and CO, we evaluate the error of the HCO⁺ abundance considering only that of the HCO⁺ column density. Overall, the HCO⁺ abundance is in a range of 3 × 10⁻¹¹ – 3 × 10⁻¹⁰ outside ~100 au, and tends to decline towards the disk center at the inner radii.
<figure><img src="image_9.png" /><figcaption>Figure 7. (cont)</figcaption></figure>

In those inner regions, the CO-rich molecular layer extends towards the midplane, although detailed analysis by Zhang et al. (2021) shows that the gaseous CO fractional abundance does not necessarily recover its canonical value (i.e. 10⁻⁴ relative to H₂) even inside of the CO snow line. Closer to the midplane, the gas density is higher and the ionization degree, i.e. the HCO⁺ abundance, should be lower. Furthermore, HCO⁺ is suppressed by grain-surface recombination (i.e. collision with a negatively charged grain), which could be more efficient than the recombination in the gas phase at high densities. The HCO⁺ abundance averaged over the CO gas rich layer should thus be lower around and inside the CO snow line than at outer radii (see §4.3).

Although HCO⁺ is expected to be one of the major molecular ions in the warm molecular layer, it is not always the most abundant ion. The relative abundances of major molecular ions, H₃⁺, HCO⁺, and N₂H⁺, are determined by the balance between proton transfer and recombination of relevant species, i.e. H₂, CO, and N₂ (Aikawa et al. 2015). Furthermore, if the abundance of small grains is low enough to allow significant penetration of UV radiation, atomic ions such as C⁺ and S⁺ can dominate even in the CO rich layers (e.g. Aikawa et al. 2015). Therefore, the HCO⁺ abundance derived here should be considered as a lower limit of the ionization degree in the molecular layer.

## 4.2. HCO⁺ enhancement in gas gaps
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Column densities of HCO⁺ (red), N₂H⁺ (blue), and N₂D⁺ (light blue). Black lines with gray error bars depict the CO column density scaled by 10⁻⁵. The error bars correspond to 10% or 1 σ for CO (see §3.2), while they show the 16th and 84th percentile for other species. At radii where the median value of the molecular column density is lower than the value of the 16th percentile by a factor of > 10, we plot the upper limits as inverted triangles (84th percentile for closed triangles and 98th percentile for open triangles). The vertical orange bars mark the CO snowline as inferred from the model by Zhang et al. (2021).</figcaption></figure>

The HCO⁺ abundance exceeds 10⁻⁹ at 50 au ≲ R ≲ 150 au in AS 209, where the CO column density shows a depression. Favre et al. (2019) observed DCO⁺ *J* = 3 − 2 in AS 209 with high resolution (0.′′26 × 0.′′21) and found a ring-like emission with its peak at ~ 85 au. They suggest that the over-density of DCO⁺ is caused by the more efficient ionization at the radius of gas and dust depletion, which is carved by a planet of ≲ 0.3 M<sub>Jup</sub> at ~ 100 au. The high HCO⁺ abundance at 50 au ≲ R ≲ 150 au is consistent with their scenario.

Alarcón et al. (2021), on the other hand, calculated the thermal structure in two disk models for AS 209; model A assumes a smooth gas distribution and a CO abundance drop around *R* ~ 80 au, while model B assumes a constant CO abundance combined with a gas density drop around 80 au. They found that the radial pressure profile derived by Teague et al. (2018b) is shallower than in model B, and concluded that the CO column density drop at 80 au is caused mainly by a CO abundance drop rather than a gas density drop. Self-consistent modeling of the thermal structure and the HCO⁺ and DCO⁺ chemistry would be useful to distinguish between these scenarios.

We also note correlated local enhancements of HCO⁺ and DCO⁺ in HD 163296. Flaherty et al. (2017) observed DCO⁺ (*J* = 3 − 2) with a beam of 0.5″ × 0.59″ and found that DCO⁺ emission is confined to three concentric rings at 54, 124, and 214 au. The ring at 124 au is the brightest and coincides with the shoulder seen in HCO⁺ emission in Figure 3 and the local enhancement of the HCO⁺/CO column density ratio (Figure 7). Teague et al. (2018a) analyzed the disk gas kinematics around HD 163296 and found deviations from Keplerian rotation, which suggests the presence of a 1 M<sub>Jup</sub> planet at 83 au and a 1.3 M<sub>Jup</sub> planet at 137 au. Hydrodynamic models of the HD 163296 disk with these planets show gas gaps (i.e. gas density decrease) around the planet orbits. The local enhancement of DCO⁺ and HCO⁺ around 125 au could be related to such gas gaps.

In MWC 480, the CO column density shows a local minimum around the gap seen in dust continuum at ~ 76 au (Figure 7). We note that the HCO⁺ column density has a local maximum at this radius; we derived the column density distribution with radial grids of a quarter of the beam size to confirm the coincidence. It again suggests a correlation between gas gap and HCO⁺ enhancement. Observations of HCO⁺ with higher angular resolution are desirable for further studies.
## 4.3. Constraining the ionization structure with a template disk chemistry model

MAPS observations, as well as previous ALMA observations of disks in the last few years, show that detailed physical structure, e.g. size, flaring, and temperature structure vary significantly among disks (e.g. Law et al. 2021b). While source-specific modeling is desirable, it is out of the scope of the present work. However, we note that the derived abundance and its radial distribution of HCO$^+$ have common features; the HCO$^+$ abundance is ~3 × 10$^{-11}$–3 × 10$^{-10}$ at $\gtrsim$ 100 au, while it tends to decrease inwards at smaller radii. These features should reflect the basic chemistry of HCO$^+$, which is not very sensitive to the details of the disk structure. Therefore, we compare our results with some template disk models.

We adopt the disk chemistry model of Aikawa et al. (2018); the disk mass is $1.7 \times 10^{-2} M_{\odot}$ and the mass of the central star is $0.5 M_{\odot}$. Stellar UV and X-ray luminosities are $10^{31}$ erg s$^{-1}$ and $10^{30}$ erg s$^{-1}$, respectively. While the dust sedimentation is not considered, the maximum grain size is set to be 1 mm. Since the total gas-to-dust mass ratio is set to be the same as the interstellar value (~100), the abundance of small dust grains is reduced compared with the interstellar dust. The 2D ($R,Z$) distributions of the gas number density (i.e. number density of hydrogen nuclei $n_{\text{H}}$) and temperature are shown in Figure 19 in Appendix E. While Aikawa et al. (2018) varied several parameters in the calculation of the chemistry model, we here show three models: their fiducial model, a high C/O model, and a low $\zeta$ model, which is newly calculated for the present work. The fiducial model is a static disk model with a cosmic-ray ionization rate of $5 \times 10^{-17}$ s$^{-1}$. Two other ionization sources, X-rays and the decay of SLR are included in the model with the X-ray spectrum of TW Hya and the ionization rate by SLR of $\zeta_{\text{SLR}} = 1 \times 10^{-18}$ s$^{-1}$. The initial molecular abundance for the disk chemistry is set by calculating the chemical evolution from the molecular cloud formation stage to the collapse of a star-forming core. In the fiducial model, the elemental abundance of C and O are $7.82 \times 10^{-5}$ and $1.8 \times 10^{-4}$ relative to hydrogen nuclei (i.e. C/O=0.43). Disk observations in recent years, however, suggest low C/H ratios and high C/O ratios in the gas phase (e.g. Bergin et al. 2014, 2016; Zhang et al. 2021; Bosman et al. 2021). This is likely due to the removal of ice by the sedimentation and radial migration of ice-coated pebbles (Kama et al. 2016; Krijt et al. 2018). Thus, in the high C/O model, $H_2O$ is completely removed and CO is reduced by one order of magnitude in the initial molecular abundances, which results in a C/O ratio of 1.43. In addition to these two models, we run a low $\zeta$ model, in which the cosmic-ray ionization rate is set to zero. Ionization rates by X-ray and SLR are set to be the same as in the fiducial model, and the initial molecular abundance is the same as the high C/O model. Since we aim to constrain the ionization rate rather than identifying the main ionization source, the low $\zeta$ model can simply be considered as a model with a midplane ionization rate of $\zeta_{\text{mid}} = 10^{-18}$ s$^{-1}$, while $\zeta_{\text{mid}}$ is $5 \times 10^{-17}$ s$^{-1}$ in the fiducial model and high C/O model.

While molecular evolution in the disk is calculated up to $t = 1 \times 10^6$ yr, we compare our observational results with the molecular abundances at $t = 1 \times 10^5$ yr, which is comparable to the vertical mixing timescale at the radius of several tens of au assuming an $\alpha$ parameter of $10^{-3}$ for turbulent mixing (e.g. Aikawa et al. 1996). In gas-grain chemical networks, volatile species tend to be converted to less volatile ices. At $t = 10^6$ yr, a significant amount of CO and $N_2$ is converted to $CO_2$, $CH_3OH$, and $NH_3$, which could be an artifact of the static disk assumption (i.e. without mixing or accretion). The desorption rate and thermal diffusion rate (and thus the reaction rate) on the grain surface are very sensitive to temperature. In the static model, in which gas and dust stay at the same position, and thus are kept at a constant temperature, only a limited number of species can effectively diffuse and react, which results in an accumulation of specific products. At $t = 1 \times 10^5$ yr, such accumulations are less significant (see Aikawa et al. 2015, for more details).

Figure 9 shows the abundances of CO, HCO$^+$, and electrons relative to hydrogen nuclei, in the low $\zeta$ model. The abundances in the other models are shown in Figure 19 in Appendix E. We can see that CO and HCO$^+$ are mostly co-spatial and that the observed features are reproduced; the HCO$^+/\text{CO}$ abundance ratio is ~$10^{-5}$–$10^{-4}$, and the HCO$^+$ abundance is ~$10^{-10}$ in the lower part of the CO-rich layer, where the density is higher. In Figure 9, the dashed line depicts the CO snow surface and the dotted lines depict the position where the X-ray ionization rate is $5 \times 10^{-17}$ s$^{-1}$ (upper), $1 \times 10^{-18}$ s$^{-1}$ (middle), and $1 \times 10^{-19}$ s$^{-1}$ (lower). We note that the dashed line and dotted lines are roughly parallel at radii $> 50$ au, probably because both parameters are controlled by the radiation transfer from the central star (i.e. stellar irradiation and X-rays). Inside the CO snow line (where the dashed line crosses the midplane), CO is indeed abundant near the midplane, although its abundance is reduced by an order of magnitude compared to the initial value (~$10^{-5}$) due to the conversion to less volatile species such as $CO_2$ ice. In this region, the HCO$^+/\text{CO}$ abundance ratio is less than $10^{-5}$ in the low $\zeta$ model, which is also consistent with the observations. Comparing the HCO$^+$ and the electron abundances, we can see that the former is slightly lower than the latter even
<figure><img src="image_11.png" /><figcaption>**Figure 9.** The abundances of CO (left), HCO⁺ (middle), and electrons (right) relative to hydrogen nuclei in the low ζ model. The dashed line depicts the CO snow surface and the dotted lines show the positions where the X-ray ionization rate is 5 × 10⁻¹⁷ s⁻¹ (the upper dotted line), 1 × 10⁻¹⁸ s⁻¹ (the middle dotted line), and 1 × 10⁻¹⁹ s⁻¹ (the lower dotted line),</figcaption></figure>

in the layers where the HCO⁺ abundance peaks. This means that the ionization degree is actually higher than the HCO⁺ abundance; atomic ions are relatively abundant even in the CO-rich layer in the model.

Figure 10 shows the radial profiles of the column densities of CO (multiplied by 10⁻⁵), HCO⁺, N₂H⁺, and N₂D⁺ in the three models. As references, 2D distributions of absolute abundances (i.e. number densities) of major molecular ions in the low ζ model are shown in Figure 11. Similar plots for the other two disk models are shown in Figure 20 in Appendix E. The column densities in the models are almost constant outside ~100 au, while the observed values decline outwards. This discrepancy might arise because the model gas density is not exponentially tapered, but simply follows a power law up to a radius of 284 au. We thus only compare the model column densities inside ~200 au with the observed values. In all three models, the radial profile of the column density ratio N(HCO⁺)/N(CO) is similar to those in Figure 7; N(HCO⁺)/N(CO) is nearly constant (~10⁻⁵), but declines towards the disk center inside a certain radius, which corresponds to the CO snow line in the disk models.

Figure 10 shows that at radii outside the CO snow line (~20 au in the template model) the HCO⁺ column density is not sensitive to ζmid, since the contribution of the X-ray ionized region is significant in the CO-rich layer there (see Figure 11 and Figure 20). Inside the CO snow line, the CO-rich layer extends to the midplane. The HCO⁺ column density and the gas column density ratio N(HCO⁺)/N(CO) around and inside the CO snow line is thus more sensitive to ζmid.<fnref n="3" /> The peak value of the HCO⁺ column density is lower in the low ζ model than in the models with ζmid = 5 × 10⁻¹⁷ s⁻¹ by a factor of 4.6 (fiducial model) and 4.2 (high C/O model). The column density ratio of HCO⁺/CO at the HCO⁺ column density peak is also lower in the low ζ model accordingly. Around the CO snow lines predicted from the thermo-chemical model of Zhang et al. (2021), the HCO⁺ column density derived from the observation is ≥ 10¹⁴ cm⁻² with the column density ratio of HCO⁺/CO ~ 10⁻⁵ in our targeted disks, as in our fiducial model or high C/O model, except for those around IM Lup and MWC 480. It may suggest ζmid ~ 10⁻¹⁷ s⁻¹ around the CO snow line. We note, however, that the HCO⁺ column density would also depend on the total gas column density, which is not varied in the set of template models shown here. While the mass of the template disk model is within the range of estimated values for GM Aur and IM Lup, AS 209 is less massive and the disks around the Herbig Ae stars are more massive (Zhang et al. 2021) (see §4.4). Source-specific models are needed for further comparison.

### 4.4. Ionization rate traced by N₂H⁺ and N₂D⁺

The template disk models (Figures 10, 11, and 20) indicate that, in contrast to HCO⁺, the column densities of N₂H⁺ and N₂D⁺ are more sensitive to the ionization rate below the CO snow surface, where these molecular ions are relatively abundant. In the cold midplane, H₃⁺ and its deuterated isotopomers are the most abundant ions, but H₃⁺ and D₃⁺ are not observable at radio wavelengths, and H₂D⁺ and HD₂⁺ are not detected so far, possibly due to unfavorable o/p ratios (e.g. Chapillon et al. 2011). Figure 20 shows that among the observable molecular ions, N₂D⁺ is the best proxy of deuterated H₃⁺ and thus the best probe of ionization rate in the midplane. Although its abundance is lower than the sum of H₃⁺ and its deuterated isotopomers, the gas above the CO snow surface does not significantly

<footnote n="3">³ We note that the absolute abundance of HCO⁺ in the midplane is high even at radii slightly outside the CO snow line. For HCO⁺ to be the abundant ion, CO does not need to fully sublimate to the gas phase (Aikawa et al. 2015).</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Column densities of CO (multiplied by $10^{-5}$), HCO$^{+}$, N$_{2}$H$^{+}$, and N$_{2}$D$^{+}$ in the fiducial model (left), high C/O model (middle) and low $\zeta$ model (right).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** The absolute abundances (number densities) of HCO$^{+}$, N$_{2}$H$^{+}$, N$_{2}$D$^{+}$, and H$_{3}^{+}$ (and its deuterated isotopomers) relative to hydrogen nuclei in the low $\zeta$ model. The dashed line depicts the CO snow surface and the dotted lines depict the position where the X-ray ionization rate is $5 \times 10^{-17} \text{ s}^{-1}$ (the upper dotted line), $1 \times 10^{-18} \text{ s}^{-1}$ (the middle dotted line), and $1 \times 10^{-19} \text{ s}^{-1}$ (the lower dotted line),</figcaption></figure>

contribute to its column density in our template disk models. In the template disk models, the N$_{2}$H$^{+}$ and N$_{2}$D$^{+}$ column densities have a peak around R ~ 50 au, where these molecular ions are abundant in the midplane (see also Qi et al. 2019). This indicates that $\zeta_{\text{mid}}$ could be best investigated by their peak column densities, especially by the N$_{2}$D$^{+}$ column density.

The column densities of N$_{2}$H$^{+}$ and N$_{2}$D$^{+}$ in our target disks are derived by Cataldi et al. (2021) and are plotted in Figure 8. As expected from the template disk models, these column densities have a peak outside the CO snow line, except for GM Aur, for which the N$_{2}$H$^{+}$ column density is the highest at the innermost radial grid. The value at the innermost grid could, however, be overestimated, since these column densities are derived by assuming a fixed excitation temperature of 20 K. When the excitation temperature is set to be the same as the midplane temperature of the thermochemical model of Zhang et al. (2021), the N$_{2}$H$^{+}$ column density has a peak (~$10^{13} \text{ cm}^{-2}$) around R ~ 100 au (Cataldi et al. 2021).

Unlike HCO$^{+}$, the distributions of N$_{2}$H$^{+}$ and N$_{2}$D$^{+}$ are expected to anti-correlate with that of CO. We thus consider their column densities rather than the column density ratio to CO. In Figure 12, we plot the peak N$_{2}$H$^{+}$ and N$_{2}$D$^{+}$ column densities derived by Cataldi et al. (2021) versus disk mass to see if there is any trend. For each disk, a range of possible masses is shown, extending from the minimum masses to the best-fit mass given in Tables 3 and 2 of Zhang et al. (2021), respectively. While the estimated disk masses of IM Lup, GM Aur, HD 163296, and MWC 480 overlap with each other$^{4}$, their N$_{2}$H$^{+}$ and N$_{2}$D$^{+}$ peak column densities vary, which indicates that the midplane ionization rate varies among the disks. The gray squares and circles depict the peak column densities around R ~ 50 au in the fiducial

<footnote n="4">⁴ Fedele et al. (2018) reproduced their 1.3 mm dust continuum image with a disk model with a dust mass of $3.5 \times 10^{-4} \text{ M}_\odot$. If the gas/dust mass ratio is 100, the disk mass of AS 209 is also similar to the disk mass of IM Lup and GM Aur.</footnote>
<figure><img src="image_14.png" /></figure>

**Figure 12.** The N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup> column densities as a function of disk mass. For each disk, a mass range adopted from Zhang et al. (2021) is shown as horizontal error bar. The column densities correspond to the peak values, with the error bars indicating the 16th and 84th percentiles. For N<sub>2</sub>D<sup>+</sup> in GM Aur, we instead plot the 99.85th percentile at ~100 au as an upper limit. Predictions from the template models are shown by the gray points.

model and low $\zeta$ model, respectively. The gray diamonds show the values in the model with X-ray ionization only; i.e. the midplane ionization rate is $\lesssim 10^{-19}$ s<sup>-1</sup>. The peak column densities of the high C/O model are not plotted but are similar to those of the fiducial model.

The N<sub>2</sub>H<sup>+</sup> column densities of IM Lup, GM Aur, and AS 209 indicate $\zeta_{\text{mid}} \sim 10^{-18}$ s<sup>-1</sup>. For MWC 480 and HD 163296, the N<sub>2</sub>H<sup>+</sup> column density is lower than that of the X-ray only model, which indicates that the X-ray ionization is less effective in these disks than in our template models. We note that a fraction of N<sub>2</sub>H<sup>+</sup> exists in the layer with relatively high X-ray ionization rate (>$\sim 5 \times 10^{-17}$ s<sup>-1</sup>), which sets the floor value of the N<sub>2</sub>H<sup>+</sup> column density in our template models. For MWC 480, the low N<sub>2</sub>H<sup>+</sup> column density is consistent with its soft X-ray spectrum and relatively low HCO<sup>+</sup> column density. Alternative explanation for the variety of N<sub>2</sub>H<sup>+</sup> column density would be the temperature structure in the disk. The column density of N<sub>2</sub>H<sup>+</sup> should be lower if the midplane layer of the low temperature (>$\sim 20$ K) is thinner (Qi et al. 2019). This could be the case for MWC 480 disk, in which both N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup> column densities are low. HD 163296, on the other hand, has relatively high column density of N<sub>2</sub>D<sup>+</sup>, which indicate that there are plenty of cold gas in the midplane (see below).

The N<sub>2</sub>D<sup>+</sup> peak column densities of IM Lup, AS 209, and HD 163296 are similar to or slightly larger than the value from the low $\zeta$ model, which suggests $\zeta_{\text{mid}} \gtrsim 10^{-18}$ s<sup>-1</sup>. In the HD 163296 disk, the peak N<sub>2</sub>D<sup>+</sup> column density is relatively high, while the low N<sub>2</sub>H<sup>+</sup> column density indicates lower X-ray ionization rate than the template models. The midplane ionization rate, which is better traced by N<sub>2</sub>D<sup>+</sup>, would thus be set by SLR or high energy particles. In GM Aur, on the other hand, the upper limit of N<sub>2</sub>D<sup>+</sup> column density indicates a low ionization rate ($< 10^{-18}$ s<sup>-1</sup>) in the midplane.

Finally we note that the N<sub>2</sub>H<sup>+</sup> and N<sub>2</sub>D<sup>+</sup> column densities depend not only on $\zeta_{\text{mid}}$, but also on other parameters such as the vertical temperature distribution and the disk mass (Cleeves et al. 2014; Qi et al. 2019). Source specific models are needed to estimate the midplane ionization rate more quantitatively for each object.

## 5. CONCLUSIONS

We observed and analyzed HCO<sup>+</sup> *J* = 1–0 and H<sup>13</sup>CO<sup>+</sup> *J* = 1–0 lines towards the protoplanetary disks around IM Lup, GM Aur, AS 209, HD 163296, and MWC 480. H<sup>13</sup>CO<sup>+</sup> *J* = 1–0 was detected in all five disks, while H<sup>13</sup>CO<sup>+</sup> *J* = 1–0 was detected (SNR> 6σ) towards GM Aur and HD 163296 and tentatively detected (SNR> 3σ) towards the other disks by a matched filter analysis in the uv plane. The disk-integrated flux of H<sup>13</sup>CO<sup>+</sup> *J* = 1–0 is also above 3σ in IM Lup, GM Aur, and HD 163296.
We derived the HCO⁺ column density by fitting the azimuthally averaged spectra of HCO⁺ and H¹³CO⁺ simultaneously for radial bins of half the beam FWHM. In all five disks, the HCO⁺ column density increases inwards, but becomes flat or drops towards the center inside a radius of *R* ~ 100 au. The column density ratio of *N*(HCO⁺)/*N*(CO) is about 10⁻⁵ – 10⁻⁴ at *R* ≳ 100 au, except for the MWC 480 disk, in which the column density ratio is < 10⁻⁵ at *R* ≲ 200 au.

We derived the HCO⁺ abundance in the warm CO-rich layer, where HCO⁺ is expected to be the dominant molecular ion, via the column density ratio *N*(HCO⁺)/*N*(CO) using the radial profiles of CO column density and CO depletion factor from Zhang et al. (2021). Beyond ~100 au, the derived HCO⁺ abundance ranges from 3 × 10⁻¹¹ – 3 × 10⁻¹⁰ in the IM Lup, GM Aur, HD 163296 disks. The HCO⁺ abundance is lower in MWC 480, possibly due to the lack of high energy (> a few keV) X-rays. The HCO⁺ abundance tends to decline at towards the disk center for *R* ≲ 100 au. This can be explained by the lower ionization degree in denser gas, especially inside the CO snow line, where the CO-rich layer is in the midplane.

We find a hint of a correlation between the HCO⁺ abundance and the gap carved by a putative planet in AS 209: the HCO⁺ abundance exceeds 10⁻⁹ at *R* ~ 50 – 150 au, where the CO column density is depressed. This region seems to correspond to the gas gap, where ionization would be more efficient. In HD 163296, the HCO⁺ abundance shows a shallow bump at 100 – 150 au, which coincides with the radius of DCO⁺ enhancement found by Flaherty et al. (2017). This feature could also be related to the enhanced ionization around a gas gap. In MWC 480, the HCO⁺ column density has a local maximum at ~ 76 au, which coincides with the local depression of CO and the gap seen in the dust continuum.

Finally, we compared the column densities of HCO⁺, N₂H⁺, and N₂D⁺ with those of template disk models: a fiducial disk model with a midplane ionization rate ζₘid of 5 × 10⁻¹⁷ s⁻¹, a model with CO and H₂O depletion, and a model with ζₘid = 1 × 10⁻¹⁸ s⁻¹. The almost constant HCO⁺ abundance at *R* ≳ 100 au is explained by X-ray ionization in the CO-rich layer. The decline of the HCO⁺ abundance at the inner radii, on the other hand, can be explained by CO sublimation inside the CO snow line; the CO-rich layer then extends to the midplane, where the ionization degree is low due to high density. While the estimated disk mass range of IM Lup, GM Aur, HD 163296, and MWC 480 overlap with each other, their peak column densities of N₂H⁺ and N₂D⁺ vary, which may indicate that the midplane ionization rate varies among disks. The peak N₂D⁺ column density suggests a midplane ionization rate of ≳ 10⁻¹⁸ s⁻¹ for IM Lup, AS 209, and HD 163296, while the upper limit of N₂D⁺ column density indicates it is < 10⁻¹⁸ s⁻¹ for GM Aur. The peak column density of N₂H⁺ is lower in MWC 480 than in other disks, which is consistent with its soft X-ray spectrum and relatively low HCO⁺ column density. Alternatively, the low column densities of N₂H⁺ and N₂D⁺ in the MWC 480 disk could be due to its warmer temperature; i.e. cold miplane layer is thinner than in other disks. Source-specific models are needed for further evaluation of the midplane ionization rate.

ACKNOWLEDGMENTS

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. 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.

We would like to thank the anonymous referee for careful reading of our manuscript and for constructive comments. Y.A. acknowledges support by NAOJ ALMA Scientific Research Grant code 2019-13B, Grant-in-Aid for Scientific Research (S) 18H05222, and Grant-in-Aid for Transformative Research Areas (A) 20H05844 and 20H05847. 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. 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. S.M.A. and J.H. acknowledge funding support from the National Aeronautics and Space Administration under Grant No. 17-XRP17 2-0012 issued through the Exoplanets Research Program. E.A.B. and A.D.B. acknowledges support from NSF AAG Grant #1907653. L.I.C. gratefully acknowledges support from the David and Lucille Packard
Foundation and Johnson & Johnson's WiSTEM2D Program. V.V.G. acknowledges support from FONDECYT Iniciación 11180904 and ANID project Basal AFB-170002. J.D.I. acknowledges support from the Science and Technology Facilities Council of the United Kingdom (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. 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"). 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). R.T. acknowledges support from the Smithsonian Institution as a Submillimeter Array (SMA) Fellow. T.T. is supported by JSPS KAK-ENHI Grant Numbers JP17K14244 and JP20K04017. C.W. acknowledges financial support from the University of Leeds, STFC and UKRI (grant numbers ST/R000549/1, ST/T000287/1, MR/T040726/1). K.Z., K.R.S., J.H., J.B., J.B.B., and I.C. acknowledge the support of NASA through Hubble Fellowship grants HST-HF2-51401.001, HST-HF2-51419.001, HST-HF2-51460.001-A, HST-HF2-51427.001-A, HST-HF2-51429.001-A, and HST-HF2-51405.001-A 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.

*Facilities: ALMA*

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

# APPENDIX

## A. HYBRID ZEROTH MOMENT MAPS WITH A 0.′5 BEAM

While we used the data cubes tapered to a circular 0.′3 beam for the analysis, faint emission is better recognized in images with lower spatial resolution. Since Band 3 lines tend to be fainter than Band 6 lines, the MAPS collaboration also produced the data cube tapered to a circular 0.′5 beam for Band 3 lines (Oberg et al. 2021; Czekala et al. 2021). Figure 13 shows the hybrid zeroth moment maps of the $J = 1-0$ lines of HCO$^+$ and H$^{13}$CO$^+$ tapered to a circular 0.′5 beam. All maps were produced by combining a Keplerian mask and a smoothed $0\sigma$-clip mask (Law et al. 2021a).

## B. MATCHED FILTER ANALYSIS AND DISK-INTEGRATED SPECTRA FOR H$^{13}$CO$^+$ $J = 1-0$

The $J = 1-0$ line of H$^{13}$CO$^+$ is not clearly seen in our zeroth moment maps (Figure 2), radial emission profiles (Figure 3), or azimuthally averaged spectra (Figure 16a). To test whether the line is detected, we applied a matched filter analysis in the uv plane (Loomis et al. 2018) (see also Czekala et al. 2021). We calculate the visibilities of a Keplerian disk model, which are then cross-correlated with the observed visibilities to produce a response for each channel of the data. The response spectrum is then divided by the rms of the response spectrum at signal-free regions. Since we do not know the radial distribution of the emission, we computed the filter response for a series of Keplerian disk models with varying radial extents to find the Keplerian model that maximizes the SNR at the systemic velocity. Figure 14 shows the matched filter responses with the maximum SNR, while the radial extent of the corresponding Keplerian disk models is shown in the figure and also listed in Table 2. We note that the radial extent of the keplerian disk model roughly represents the radial extent of the emission. H$^{13}$CO$^+$ $J = 1-0$ is detected with a SNR above $6\sigma$ for GM Aur and HD 163296. The line is tentatively detected with a SNR above $3\sigma$ for IM Lup, AS 209 and MWC 480.

In Figure 15 we show the shifted and azimuthally averaged, disk-integrated H$^{13}$CO$^+$ $J = 1-0$ spectrum of each disk. The radial extent of the region over which we averaged is the same as for the matched filter. The disk-integrated spectra also point towards a detection in GM Aur and HD 163296, although at lower significance than the matched filter. In summary, we detect H$^{13}$CO$^+$ $J = 1-0$ towards GM Aur and IM Lup, and tentatively detect it towards the other three disks.
<figure><img src="image_15.png" /><figcaption>Figure 13. Same as figure 2, but for images tapered to a circular 0.5" beam.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 14. Matched filter response for H$^{13}$CO$^{+}$ J = 1 − 0. For each disk, we chose the Keplerian mask with the radius that maximized the SNR. The radial extent of that model is indicated in the upper right of each panel. The y-axis is in units of the noise in the filter response.</figcaption></figure>
<figure><img src="image_17.png" /></figure>

**Figure 15.** Disk-integrated spectra (over the radial range indicated in the title of each panel) of H<sup>13</sup>CO<sup>+</sup> *J* = 1 − 0. The spectrum at each spatial location of the data cube was shifted to the systemic velocity prior to integration. The shaded region corresponds to the 1σ error.

### C. AZIMUTHALLY AVERAGED SPECTRA

Figure 16 show the azimuthally averaged spectra of HCO<sup>+</sup> and H<sup>13</sup>CO<sup>+</sup> *J* = 1 − 0 together with model spectra from the MCMC fits.

### D. DEPENDENCE OF HCO<sup>+</sup> COLUMN DENSITY ON ASSUMED EXCITATION TEMPERATURE

In our fiducial fits, the excitation temperature is a free parameter. In order to explore the dependence of our results on the excitation temperature, we performed additional fits where the excitation temperature is fixed to 30 K. Figure 17 displays a comparison of the column density profiles derived from the different fits. The comparison is discussed in Section 3.2. Figure 18 shows the optical depth profiles derived assuming T<sub>ex</sub> = 30 K.

### E. TEMPLATE DISK CHEMISTRY MODELS

To top panels in Figure 19 show the 2D (*R*, *Z*) distributions of the number of hydrogen nuclei and gas temperature in the template disk model. The panels in the second row show the abundances of gaseous CO, HCO<sup>+</sup>, and electron relative to hydrogen nuclei in the fiducial model. The panels in the third and fourth row are for the high C/O model and the low *ζ* model.

Figure 20 shows the absolute abundances (i.e. number density) of HCO<sup>+</sup> (top row), N<sub>2</sub>H<sup>+</sup> (the second row), N<sub>2</sub>D<sup>+</sup> (the third row), H<sub>3</sub><sup>+</sup> and its deuterated isotopomers (the fourth row) in the fiducial model (left), the high C/O model (middle), and the low *ζ* model (right).

### REFERENCES

Aikawa, Y., Furuya, K., Hincelin, U., & Herbst, E. 2018, ApJ, 855, 119, doi: 10.3847/1538-4357/aaad6c

Aikawa, Y., Furuya, K., Nomura, H., & Qi, C. 2015, ApJ, 807, 120, doi: 10.1088/0004-637X/807/2/120
<figure><img src="image_18.png" /></figure>

**Figure 16a.** The black curves show azimuthally averaged spectra for IM Lup. Spectra are vertically offset for clarity. The shaded region marks the $1\sigma$ error. The horizontal dotted line marks the zero flux level. For each spectrum, we show 20 randomly selected model spectra from the MCMC with the blue curves ($T_{ex}$ as free parameter) and the red curves ($T_{ex}$ fixed to 30 K). The selection probability of each model was set proportional to its posterior probability. Most of the time, the model spectra overlap closely, so that only the blue curves are visible.
<figure><img src="image_19.png" /><figcaption>Figure 16b. Continuation of figure 16a.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 16c. Same as figure 16a, but for GM Aur.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure 16d. Same as figure 16a, but for AS 209.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure 16e. Same as figure 16a, but for HD 163296.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure 16f. Continuation of figure 16e.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Figure 16g. Same as figure 16a, but for MWC 480.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Figure 17.** The column density of HCO$^+$ assuming an excitation temperature of 30 K (top row) and for our fiducial fits where $T_{\text{ex}}$ is a free parameter (bottom row). The black solid lines mark the median, while the shading encompasses the regions from the 16th to 84th and the 2.3th to 97.7th percentile, respectively. To ease comparison, in the top panel, the thick and thin black dotted lines mark the median and the 16th and 84th percentile of the fits where $T_{\text{ex}}$ is a free parameter.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Figure 18.** Optical depth profiles of HCO$^+$ ($J = 1 - 0$) and H$^{13}$CO$^+$ ($J = 1 - 0$) derived from fitting azimuthally averaged spectra assuming $T_{\text{ex}} = 30$ K. The black lines mark the median, while the shading encompasses the regions from the 16th to 84th and the 2.3th to 97.7th percentile, respectively.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>**Figure 19.** Number density of hydrogen nuclei ($n_H$) and gas temperature of the template disk model are shown in the top row. Furthermore, we show the abundances of CO (left column), HCO$^+$ (middle column), and electrons (right column) relative to hydrogen nuclei in the fiducial model (second row), CO depleted model (third row), and the model without CR ionization (low $\zeta$ model, bottom row).</figcaption></figure>

Aikawa, Y., & Herbst, E. 2001, A&A, 371, 1107,
doi: 10.1051/0004-6361:20010416

Aikawa, Y., Miyama, S. M., Nakano, T., & Umebayashi, T.
1996, ApJ, 467, 684, doi: 10.1086/177644

Aikawa, Y., & Nomura, H. 2006, ApJ, 642, 1152,
doi: 10.1086/501114

Alarcón, F., Bosman, A., Bergin, E., et al. 2021, arXiv
e-prints, arXiv:2109.06263.
https://arxiv.org/abs/2109.06263

Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJL,
869, L41, doi: 10.3847/2041-8213/aaf741

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J.,
et al. 2013, A&A, 558, A33,
doi: 10.1051/0004-6361/201322068

Bai, X.-N., & Stone, J. M. 2013, ApJ, 769, 76,
doi: 10.1088/0004-637X/769/1/76

Bergin, E. A., Aikawa, Y., Blake, G. A., & van Dishoeck,
E. F. 2007, in Protostars and Planets V, ed. B. Reipurth,
D. Jewitt, & K. Keil, 751.
https://arxiv.org/abs/astro-ph/0603358

Bergin, E. A., Cleeves, L. I., Crockett, N., & Blake, G. A.
2014, Faraday Discussions, 168, 61,
doi: 10.1039/C4FD00003J
<figure><img src="image_28.png" /><figcaption>**Figure 20.** Absolute abundance (i.e. number density in units of cm⁻³) of major molecular ions in the fiducial disk model (left column), high C/O model (middle column), and low ζ model (right column).</figcaption></figure>

Bergin, E. A., Du, F., Cleeves, L. I., et al. 2016, ApJ, 831, 101, doi: 10.3847/0004-637X/831/1/101

Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75, doi: 10.1051/0004-6361/201630056

Bosman, A. D., Walsh, C., & van Dishoeck, E. F. 2018, A&A, 618, A182, doi: 10.1051/0004-6361/201833497

Bosman, A. D., Alarcón, F., Bergin, E. A., et al. 2021, arXiv e-prints, arXiv:2109.06221. https://arxiv.org/abs/2109.06221

Caselli, P., Walmsley, C. M., Zucconi, A., et al. 2002, ApJ, 565, 344, doi: 10.1086/324302

Cataldi, G., Yamato, Y., Aikawa, Y., et al. 2021, arXiv e-prints, arXiv:2109.06462. https://arxiv.org/abs/2109.06462

Chapillon, E., Parise, B., Guilloteau, S., & Du, F. 2011, A&A, 533, A143, doi: 10.1051/0004-6361/201117354

Cieza, L. A., González-Ruilova, C., Hales, A. S., et al. 2021, MNRAS, 501, 2934, doi: 10.1093/mnras/staa3787

Cleeves, L. I., Bergin, E. A., & Adams, F. C. 2014, ApJ, 794, 123, doi: 10.1088/0004-637X/794/2/123

Cleeves, L. I., Bergin, E. A., Qi, C., Adams, F. C., & Öberg, K. I. 2015, ApJ, 799, 204, doi: 10.1088/0004-637X/799/2/204

Czekala, I., Loomis, R. A., Teague, R., et al. 2021, arXiv e-prints, arXiv:2109.06188. https://arxiv.org/abs/2109.06188

Dalgarno, A. 2006, Proceedings of the National Academy of Science, 103, 12269, doi: 10.1073/pnas.0602117103
Dionatos, O., Woitke, P., Güdel, M., et al. 2019, A&A, 625, A66, doi: 10.1051/0004-6361/201832860

Duley, W. W., & Williams, Williams, A. 1984, Interstellar Chemistry (Academic Press)

Endres, C. P., Schlemmer, S., Schilke, P., Stutzki, J., & Müller, H. S. P. 2016, Journal of Molecular Spectroscopy, 327, 95, doi: 10.1016/j.jms.2016.03.005

Favre, C., Fedele, D., Maud, L., et al. 2019, ApJ, 871, 107, doi: 10.3847/1538-4357/aaf80c

Fedele, D., Tazzari, M., Booth, R., et al. 2018, A&A, 610, A24, doi: 10.1051/0004-6361/201731978

Flaherty, K. M., Hughes, A. M., Rose, S. C., et al. 2017, ApJ, 843, 150, doi: 10.3847/1538-4357/aa79f9

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067

Fujimoto, Y., Krumholz, M. R., & Tachibana, S. 2018, MNRAS, 480, 4025, doi: 10.1093/mnras/sty2132

Furuya, K., & Aikawa, Y. 2014, ApJ, 790, 97, doi: 10.1088/0004-637X/790/2/97

Glassgold, A. E., Najita, J., & Igea, J. 1997, ApJ, 480, 344, doi: 10.1086/303952

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJL, 869, L42, doi: 10.3847/2041-8213/aaf740

Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55

Jonkheid, B., Dullemond, C. P., Hogerheijde, M. R., & van Dishoeck, E. F. 2007, A&A, 463, 203, doi: 10.1051/0004-6361:20065668

Kama, M., Bruderer, S., van Dishoeck, E. F., et al. 2016, A&A, 592, A83, doi: 10.1051/0004-6361/201526991

Krijt, S., Schwarz, K. R., Bergin, E. A., & Ciesla, F. J. 2018, ApJ, 864, 78, doi: 10.3847/1538-4357/aad69b

Law, C. J., Loomis, R. A., Teague, R., et al. 2021a, arXiv e-prints, arXiv:2109.06210. https://arxiv.org/abs/2109.06210

Law, C. J., Teague, R., Loomis, R. A., et al. 2021b, arXiv e-prints, arXiv:2109.06217. https://arxiv.org/abs/2109.06217

Liu, Y., Dipierro, G., Ragusa, E., et al. 2019, A&A, 622, A75, doi: 10.1051/0004-6361/201834157

Long, F., Pinilla, P., Herczeg, G. J., et al. 2018, ApJ, 869, 17, doi: 10.3847/1538-4357/aae8e1

Loomis, R. A., Öberg, K. I., Andrews, S. M., et al. 2018, AJ, 155, 182, doi: 10.3847/1538-3881/aab604

Matrà, L., MacGregor, M. A., Kalas, P., et al. 2017, ApJ, 842, 9, doi: 10.3847/1538-4357/aa71b4

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 376, Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw, F. Hill, & D. J. Bell, 127

Milam, S. N., Savage, C., Brewster, M. A., Ziurys, L. M., & Wyckoff, S. 2005, ApJ, 634, 1126, doi: 10.1086/497123

Müller, H. S. P., Thorwirth, S., Roth, D. A., & Winnewisser, G. 2001, A&A, 370, L49, doi: 10.1051/0004-6361:20010367

Öberg, K. I., Qi, C., Wilner, D. J., & Andrews, S. M. 2011a, ApJ, 743, 152, doi: 10.1088/0004-637X/743/2/152

Öberg, K. I., Qi, C., Fogel, J. K. J., et al. 2011b, ApJ, 734, 98, doi: 10.1088/0004-637X/734/2/98

Oberg, K. I., Guzman, V. V., Walsh, C., et al. 2021, arXiv e-prints, arXiv:2109.06268. https://arxiv.org/abs/2109.06268

Padovani, M., Ivlev, A. V., Galli, D., & Caselli, P. 2018, A&A, 614, A111, doi: 10.1051/0004-6361/201732202

Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f

Qi, C., Öberg, K. I., & Wilner, D. J. 2013a, ApJ, 765, 34, doi: 10.1088/0004-637X/765/1/34

Qi, C., Öberg, K. I., Wilner, D. J., et al. 2013b, Science, 341, 630, doi: 10.1126/science.1239560

Qi, C., Öberg, K. I., Espaillat, C. C., et al. 2019, ApJ, 882, 160, doi: 10.3847/1538-4357/ab35d3

Rab, C., Güdel, M., Padovani, M., et al. 2017, A&A, 603, A96, doi: 10.1051/0004-6361/201630241

Rab, C., Güdel, M., Woitke, P., et al. 2018, A&A, 609, A91, doi: 10.1051/0004-6361/201731443

Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2018, ApJ, 856, 85, doi: 10.3847/1538-4357/aaae08

Sierra, A., Pérez, L. M., Zhang, K., et al. 2021, arXiv e-prints, arXiv:2109.06433. https://arxiv.org/abs/2109.06433

Smirnov-Pinchukov, G. V., Semenov, D. A., Akimkin, V. V., & Henning, T. 2020, A&A, 644, A4, doi: 10.1051/0004-6361/202038572

Suzuki, T. K., & Inutsuka, S.-i. 2009, ApJL, 691, L49, doi: 10.1088/0004-637X/691/1/L49

Teague, R. 2019, The Journal of Open Source Software, 4, 1632, doi: 10.21105/joss.01632

Teague, R., Bae, J., Bergin, E. A., Birnstiel, T., & Foreman-Mackey, D. 2018a, ApJL, 860, L12, doi: 10.3847/2041-8213/aac6d7

Teague, R., Bae, J., Birnstiel, T., & Bergin, E. A. 2018b, ApJ, 868, 113, doi: 10.3847/1538-4357/aae836
Teague, R., & Foreman-Mackey, D. 2018, Research Notes of
the American Astronomical Society, 2, 173,
doi: 10.3847/2515-5172/aae265

Teague, R., Semenov, D., Guilloteau, S., et al. 2015, A&A,
574, A137, doi: 10.1051/0004-6361/201425268

Teague, R., Guilloteau, S., Semenov, D., et al. 2016, A&A,
592, A49, doi: 10.1051/0004-6361/201628550

Umebayashi, T., & Nakano, T. 1981, PASJ, 33, 617

—. 1988, Progress of Theoretical Physics Supplement, 96,
151, doi: 10.1143/PTPS.96.151

—. 2009, ApJ, 690, 69, doi: 10.1088/0004-637X/690/1/69

van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011,
Computing in Science and Engineering, 13, 22,
doi: 10.1109/MCSE.2011.37

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020,
Nature Methods,
doi: https://doi.org/10.1038/s41592-019-0686-2

Wakelam, V., Chapillon, E., Dutrey, A., et al. 2019,
MNRAS, 484, 1563, doi: 10.1093/mnras/stz068

Willacy, K. 2007, ApJ, 660, 441, doi: 10.1086/512796

Wilson, T. L. 1999, Reports on Progress in Physics, 62,
143, doi: 10.1088/0034-4885/62/2/002

Yen, H.-W., Koch, P. M., Liu, H. B., et al. 2016, ApJ, 832,
204, doi: 10.3847/0004-637X/832/2/204

Zhang, K., Booth, A. S., Law, C. J., et al. 2021, arXiv
e-prints, arXiv:2109.06233.
https://arxiv.org/abs/2109.06233