# The Chemical Inventory of the Planet-hosting Disk PDS 70

Stefano Facchini, Richard Teague, Jaehan Bae, Myriam Benisty, Miriam Keppler, Andrea Isella

► To cite this version:

Stefano Facchini, Richard Teague, Jaehan Bae, Myriam Benisty, Miriam Keppler, et al.. The Chemical Inventory of the Planet-hosting Disk PDS 70. The Astronomical Journal, 2021, 162, 10.3847/1538-3881/abf0a4. insu-03705282

HAL Id: insu-03705282

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

Submitted on 3 Aug 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
<figure><img src="image_2.png" /></figure>

# The Chemical Inventory of the Planet-hosting Disk PDS 70

Stefano Facchini<fnref n="1" />, Richard Teague<fnref n="2" />, Jaehan Bae<fnref n="3" />,<fnref n="8" />, Myriam Benisty<fnref n="4" />,<fnref n="5" />, Miriam Keppler<fnref n="6" />, and Andrea Isella<fnref n="7" />

## Abstract

As host to two accreting planets, PDS 70 provides a unique opportunity to probe the chemical complexity of atmosphere-forming material. We present ALMA Band 6 observations of the PDS 70 disk and report the first chemical inventory of the system. With a spatial resolution of 0.′′4–0.′′5 (~50 au), 12 species are detected, including CO isotopologs and formaldehyde, small hydrocarbons, HCN and HCO⁺ isotopologs, and S-bearing molecules. SO and CH₃OH are not detected. All lines show a large cavity at the center of the disk, indicative of the deep gap carved by the massive planets. The radial profiles of the line emission are compared to the (sub)millimeter continuum and infrared scattered light intensity profiles. Different molecular transitions peak at different radii, revealing the complex interplay between density, temperature, and chemistry in setting molecular abundances. Column densities and optical depth profiles are derived for all detected molecules, and upper limits obtained for the nondetections. Excitation temperature is obtained for H₂CO. Deuteration and nitrogen fractionation profiles from the hydrocyanide lines show radially increasing fractionation levels. Comparison of the disk chemical inventory to grids of chemical models from the literature strongly suggests a disk molecular layer hosting a carbon-to-oxygen ratio C/O > 1, thus providing for the first time compelling evidence of planets actively accreting high C/O ratio gas at present time.

*Unified Astronomy Thesaurus concepts: Astrochemistry (75); Protoplanetary disks (1300); Exoplanet formation (492); Interferometry (808); Chemical abundances (224); Circumstellar gas (238)*

## 1. Introduction

The frequency and characteristics of exoplanets detected in the last two decades indicate that planet formation is a robust and efficient mechanism that can lead to an outstanding diversity in the orbital architectures of planet systems and the physical properties of the planets within them (Winn & Fabrycky 2015). The properties and architecture of exoplanetary systems are likely inherited from the physical and chemical processes that led to their formation in their birth environment, the protoplanetary disks. Key properties such as mass, radius, orbital period, and multiplicity are likely determined by the formation mechanisms and subsequent evolution of planetary systems (Mordasini 2018). Possible signatures of forming planets, and of planet–disk interactions, are routinely detected as substructures in high-resolution continuum observations of protoplanetary disks (Andrews 2020) and as local deviations from Keplerian rotation in the gas velocity field (Disk Dynamics Collaboration et al. 2020). The tremendous advancement in observations of protoplanetary disks can now allow us to link the properties of the mature exoplanet population with the physical and chemical properties of their natal environments.

The chemical composition of planetary atmospheres, in particular, is thought to hold an imprint of the planet formation history, as the elemental composition directly traces the elemental abundances that have been accreted onto the planet surface through its formation phases (Öberg et al. 2011; Öberg & Bergin 2016). Atmospheric retrieval techniques have been

used to determine the molecular and atomic abundances of a few species of massive exoplanets. With the aid of sophisticated chemical models, transmission and thermal emission spectra of super-Jovian planets allow an inference of the elemental ratio present in their atmosphere, in particular the C/O ratio, which is critical in determining the chemical processes that occur within it (e.g., Madhusudhan 2019; Mollière et al. 2020). At the same time, molecular line observations of protoplanetary disks have been used to derive volatile C/O ratios. C and O abundances in the warm molecular layer have been determined in a handful of objects, showing that elements are reprocessed between ice and volatile phases through physical and chemical processes (Bergin et al. 2016; Cleeves et al. 2018).

The presence of massive planets embedded in their natal protoplanetary disks directly affects the physical and chemical properties of the material that they can accrete. On one hand, giant planets induce surface density and opacity variations across the disk, which change the thermal and ionization structure of the disk (e.g., Facchini et al. 2018; Isella & Turner 2018; Favre et al. 2019). On the other hand, protoplanetary disks evolve both physically and chemically (e.g., Eistrup et al. 2018), with a complex evolutionary path that is bound to determine the chemical abundances of gas and dust material that is accreted by embedded planets. In order to establish the chemical properties of the material being accreted by massive planets, high-angular-resolution observations probing the chemical abundances of the hosting disks in the proximity of planets are needed.

<footnote n="1">¹European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching, Germany; stefano.facchini@eso.org</footnote>

<footnote n="2">²Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

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

<footnote n="4">⁴Unidad Mixta Internacional Franco-Chilena de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile</footnote>

<footnote n="5">⁵Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

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

<footnote n="7">⁷Department of Physics and Astronomy, Rice University, 6100 Main Street, MS-108, Houston, TX 77005, USA</footnote>

<footnote>Received 2021 January 20; revised 2021 March 8; accepted 2021 March 18; published 2021 August 11</footnote>

<footnote n="8">⁸NASA Hubble Fellowship Program Sagan Fellow.</footnote>
<table><caption>**Table 1**
Log of ALMA Observations Used in This Paper</caption><thead><tr><th>Frequency Range</th><th>Date</th><th>Antennas<sup>a</sup></th><th>Baselines</th><th>Time on Source</th><th>PWV</th><th>Bandpass/flux Calibrator</th><th>Phase Calibrator</th></tr></thead><tbody><tr><td rowspan="2">217.194–233.955 GHz</td><td>06 Mar 2020</td><td>43</td><td>15–783m</td><td>42 minutes</td><td>2.2 mm</td><td>J1427–4206</td><td>J1407–4302</td></tr><tr><td>16 Mar 2020</td><td>46</td><td>15–969m</td><td>42 minutes</td><td>2.7 mm</td><td>J1427–4206</td><td>J1407–4302</td></tr><tr><td rowspan="2">262.081 GHz</td><td>02 Mar 2020</td><td>44</td><td>15–783 m</td><td>38 minutes</td><td>1.6 mm</td><td>J1427–4206</td><td>J1407–4302</td></tr><tr><td>06 Mar 2020</td><td>44</td><td>15–783 m</td><td>38 minutes</td><td>2.4 mm</td><td>J1427–4206</td><td>J1407–4302</td></tr></tbody></table>

**Notes.** The spectral setup is specified in Table 2.

<sup>a</sup> Number of antennas after flagging.

This paper presents the first chemical inventory of the only disk known to harbor massive protoplanets, PDS 70, and as such lays a foundation stone in establishing the chemicophysical connection between massive planets and their host environment. The ~5 Myr old young star PDS 70 hosts two giant planets, PDS 70b and PDS 70 c, recently directly imaged at multiple wavelengths (e.g., Keppler et al. 2018; Müller et al. 2018; Wagner et al. 2018; Haffert et al. 2019; Christiaens et al. 2019; Mesa et al. 2019). The planets are still accreting material from their host disk, as indicated by their strong emission in the Hα line (Wagner et al. 2018; Haffert et al. 2019). The two planets are in a dynamically stable configuration, near 2:1 mean motion resonance (Bae et al. 2019; Wang et al. 2021), with PDS 70b following an eccentric orbit and PDS 70 c in a nearly circular one (Wang et al. 2021). Other physical properties of the planets such as their mass are still rather unconstrained due to the likely presence of circumplanetary dust (e.g., Isella et al. 2019; Stolker et al. 2020) affecting their spectra. First attempts in characterizing their atmospheric properties are being performed, but their C/O ratio is still unconstrained (Wang et al. 2021). The planets are observed to sculpt the disk structure in both its gas and dust components. They are carving a large cavity in the dust traced by continuum millimeter emission, with a deep gap traced by CO emission at the orbital radius of the two planets (Hashimoto et al. 2012; Dong et al. 2012; Hashimoto et al. 2015; Keppler et al. 2019).

PDS 70 thus provides the best candidate for where to establish a link between massive protoplanets and their host environment, as the disk gas and dust structure, planetary accretion rates, masses, orbital radii, and dynamical stability have been determined observationally. In this paper, we present a spatially resolved chemical inventory of the disk around PDS 70 using observations obtained with the Atacama Large Millimeter Array (ALMA). We demonstrate how the embedded planets affect the radial dependence of the molecular properties of their hosting disk, while the disk chemistry is bound to determine the chemical abundances of the building planetary atmospheres. In Section 2, we describe the observations and data reduction. In Section 3, we present the first observational results of the chemical survey, while in Section 4 we derive physical properties for the detected molecular transitions. Finally, in Section 5 we discuss our findings and summarize our results in Section 6.

## 2. Observations and Data Reduction

This paper mostly presents results from ALMA Band 6 observations taken in 2020 March as part of the ALMA program #2019.1.01619.S (PI S. Facchini). These data represent the short spacing of a program aiming at imaging molecular lines at 0′′1 resolution in the system, for which the long-baseline data have not been taken yet. A summary of observation dates and properties is

reported in Table 1. Two spectral setups were used to target several molecular species (see Table 2), with two execution blocks of ~40 minutes per spectral setup. In particular, the lower-frequency setup in the 217.194–233.955 GHz spectral range consisted of seven spectral windows, while the higher frequency setup had nine spectral windows. In both cases, one of the spectral windows has a bandwidth of 1.875 GHz to increase the continuum signal-to-noise ratio (S/N) and to allow for good self-calibration of the data. These spectral windows are set to Frequency Division Mode (FDM) to allow for serendipitous line detection.

The data were pipeline calibrated using the CASA package, version 5.6. After minor manual flagging, we imaged the continuum of each execution block separately and noticed that likely due to poor weather conditions the pointing was offset between execution blocks by a significant fraction of the beam. While this can be corrected for with self-calibration, we preferred to recenter each execution block beforehand in order to construct the best model for self-calibration from the first step (Andrews et al. 2018). Because PDS 70 is a transition disk with a large cavity and that no clear emission from an inner disk is present in the single execution blocks, we were unable to align the individual execution blocks by fitting the peak continuum intensity at the disk center, as is often done. We thus fitted the outer ring using a Gaussian ring model in the visibility plane, with seven free parameters (intensity normalization, radius and width of the ring, inclination, position angle, and disk center). Each execution was fitted separately using the GALARIO code (Tazzari et al. 2018), sampling the posterior distributions with the emcee package (Foreman-Mackey et al. 2013), a Python implementation of the Metropolis-Hastings ensemble sampler, with 60 walkers and 2000 steps after burn-in. Flat priors were used, with the parameters initialized from the fit to the ALMA Band 7 continuum reported in Keppler et al. (2019). Baselines were selected to be <650 kλ to have the most uniform (u, v) coverage across different execution blocks. We then shifted the phase center of each execution block using the fixvis task and realigned all the execution blocks to the same coordinate system using the fixplanet task.

Particular care was also given to the flux calibration of the observations. For all execution blocks, the Band 6 flux was corrected by deriving a spectral index from a Band 3 and Band 7 observation of the flux calibrator J1427–4207. While the Band 3 observation of the flux calibrator was always within 3 days from our science observations, the closest observation in Band 7 was from 2020 February 27, 18 days before the latest observation of PDS 70. We thus decided to correct for flux offsets between execution blocks with the same spectral setup by using the flux from the observation that was closest in time to the Band 7 calibrator estimate. To do so, we followed the procedure used in the DSHARP collaboration (Andrews et al. 2018), correcting for
<table><caption>**Table 2**

Spectral Windows of the Two Spectral Settings of the ALMA Band 6 Observations</caption><thead><tr><th>Central Rest Frequency (MHz)</th><th>Frequency Resolution (kHz)</th><th>Velocity Resolution (km s<sup>-1</sup>)</th><th>#channels</th><th>Main Lines Included in spw</th></tr></thead><tbody><tr><td>217238.530</td><td>244.141</td><td>0.389</td><td>480</td><td>DCN J = 3-2</td></tr><tr><td>218222.192</td><td>244.141</td><td>0.388</td><td>480</td><td>H<sub>2</sub>CO 3<sub>03</sub>-2<sub>02</sub></td></tr><tr><td>218732.732</td><td>244.141</td><td>0.387</td><td>960</td><td>H<sub>2</sub>CO 3<sub>21</sub>-2<sub>20</sub><br>c-C<sub>3</sub>H<sub>2</sub> 7<sub>16</sub>-7<sub>07</sub></td></tr><tr><td>219560.358</td><td>122.070</td><td>0.193</td><td>960</td><td>C<sup>18</sup>O J = 2-1</td></tr><tr><td>220398.684</td><td>122.070</td><td>0.192</td><td>960</td><td><sup>13</sup>CO J = 2-1</td></tr><tr><td>230538.000</td><td>61.035</td><td>0.092</td><td>1920</td><td><sup>12</sup>CO J = 2-1</td></tr><tr><td>233002.432</td><td>976.562</td><td>1.453</td><td>1920</td><td>CH<sub>3</sub>OH 10<sub>2</sub>-9<sub>3</sub>(A<sup>+</sup>)</td></tr><tr><td>244500.000</td><td>976.562</td><td>1.384</td><td>1920</td><td>CS J = 5-4<br>CH<sub>3</sub>OH 5<sub>1</sub>-4<sub>1</sub>(A<sup>-</sup>)<br>c-C<sub>3</sub>H<sub>2</sub> 3<sub>21</sub>-2<sub>12</sub></td></tr><tr><td>246404.588</td><td>244.141</td><td>0.343</td><td>960</td><td>SO 2<sub>3</sub>-3<sub>2</sub></td></tr><tr><td>246663.470</td><td>244.141</td><td>0.343</td><td>960</td><td><sup>34</sup>SO 5<sub>6</sub>-4<sub>5</sub></td></tr><tr><td>258156.996</td><td>244.141</td><td>0.328</td><td>480</td><td>HC<sup>15</sup>N J = 3-2</td></tr><tr><td>258255.825</td><td>244.141</td><td>0.328</td><td>480</td><td>SO 6<sub>6</sub>-5<sub>5</sub></td></tr><tr><td>259011.787</td><td>244.141</td><td>0.327</td><td>960</td><td>H<sup>13</sup>CN J = 3-2</td></tr><tr><td>260255.342</td><td>122.070</td><td>0.163</td><td>480</td><td>H<sup>13</sup>CO<sup>+</sup> J = 3-2</td></tr><tr><td>261843.721</td><td>122.070</td><td>0.162</td><td>480</td><td>SO 7<sub>6</sub>-6<sub>5</sub></td></tr><tr><td>262004.260</td><td>122.070</td><td>0.161</td><td>960</td><td>C<sub>2</sub>H N = 3-2, J = <sup>7</sup>&frasl;<sub>2</sub> - <sup>5</sup>&frasl;<sub>2</sub>, F = 4-3<br>C<sub>2</sub>H N = 3-2, J = <sup>7</sup>&frasl;<sub>2</sub> - <sup>5</sup>&frasl;<sub>2</sub>, F = 3-2<br>C<sub>2</sub>H N = 3-2, J = <sup>5</sup>&frasl;<sub>2</sub> - <sup>3</sup>&frasl;<sub>2</sub>, F = 3-2</td></tr></tbody></table>

Note. Frequency and velocity resolutions are obtained after Hanning smoothing being performed directly at the correlator.

flux offsets using the *gaincal* task. We emphasize that for the low-frequency observations the closest flux calibrator measurement in Band 7 was 8 days apart from the science observation of the first execution block.

The data sets with same spectral setup were combined using the *concat* task, and then self-calibration was performed on the two frequency settings separately. All spectral windows were included in the self-calibration, after flagging channels within ±20 km s⁻¹ (~14.5–17.5 MHz depending on the central rest frequency) from all the lines listed in Table 2. The resulting continuum was channel averaged with channels <500 MHz to avoid bandwidth smearing. At every step of the self-calibration, the model was imaged cleaning down to 5σ of the initial data, with a Briggs robust parameter of 0.5 and an elliptical mask with a position angle (PA) of 160°, semimajor axis of 1′′7, and a ratio of major to minor axis given by an inclination of 51° (Keppler et al. 2019). The rms is estimated over an elliptical annulus with the same PA and axes ratio between 2′′55 and 4′′25. We performed three rounds of phase self-calibration in both frequency setups. The first one accounted for polarization offsets with scan length as solution interval, and then two rounds combining polarizations and spectral windows with 120 and 60 s solution intervals for the low-frequency setting, and 120 and 30 s solution intervals for the high-frequency one. A final step of amplitude self-calibration at scan length was performed in both cases. The overall improvement was 119% and 54%, for a peak S/N of 434 and 503, for the low- and high-frequency setups respectively. Image reconstruction was performed with the *tclean* task by using the *multiscale* deconvolver, cleaning down to the 1σ level. More details about the continuum imaging are reported in Section 3.1. The bandwidth-weighted mean continuum frequencies are 225.56 GHz (1.33 mm) and 252.80 GHz (1.19 mm).

The self-calibration tables were then applied to the unflagged, unaveraged visibilities. Continuum subtraction was performed

using first-order fitting with the *uvcontsub* task. Imaging of the spectral lines was performed with *tclean* with the *multiscale* deconvolver, allowing for point sources in the emission (particularly important for the spatially unresolved high-velocity channels). In the cleaning process, we apply so-called Keplerian masks, i.e., exploit the rotation of the disk well fitted by a Keplerian model in the <sup>12</sup>CO J = 3–2 line to have a prior on where the molecular lines originate. We use the following geometrical parameters: $M_{\text{star}} = 0.875 M_{\odot}$, i = 51°,7, and PA = 160°,4 (Keppler et al. 2019). For bright optically thick lines (<sup>12</sup>CO and HCO<sup>+</sup>), we accounted for an elevated disk emission surface in creating the maps, using the best-fit parameters by Kepler et al. (2019): $z(r)'' = 0.33 \times (r/1')^{0.76}$, with z and r being the vertical and radial cylindrical coordinates. The channel masks were created with the *keplerian_mask.py* tool<fnref n="9" />.

As a first attempt, all lines listed in Table 2 were imaged between −6 and +17 km s⁻¹ with both natural and Briggs weighting with a robust parameter equaling zero. Given the very different line intensities spanning two orders of magnitudes, multiple imaging runs have been performed on each line to obtain the best compromise between angular resolution, intensity, and channel spacing. The choice of the imaging parameters for each line cube depends on the scientific question that is being addressed, e.g., natural weighting and uv-tapering are used to maximize line detection, while Briggs weighting with a robust parameter of 0 is used to maximize the spatial resolution to extract intensity profiles of bright lines. Table 3 summarizes the fiducial imaging parameters chosen for the lines included in this study.

Band 7 observations of HCO<sup>+</sup> J = 4–3 and <sup>12</sup>CO J = 3–2 are included in the analysis. In particular, data from programs #2015.1.00888.S and #2017.A.00006.S include the <sup>12</sup>CO

<footnote n="9">⁹ https://github.com/richteague/keplerian_mask</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1. Intensity maps of the continuum emission at 1.19 mm (left) and 1.33 mm (right), imaged with superuniform weighting.</figcaption></figure>

transition disks being on the high side within the protoplanetary disk sample (Pinilla et al. 2014; Ansdell et al. 2018; Tazzari et al. 2020). By using the upper limit on the inner disk flux density at 1.33 mm, we obtain a lower limit on the 0.85–1.33 mm spectral index of the inner disk of $1.61 \pm 0.44$, indicative that the central emission is indeed affected by the outer disk, which artificially increases the flux density at longer wavelengths, thus reducing the spectral index. Dust thermal emission at warm temperatures in the proximity of the star would imply spectral indices $\ge 2$ (e.g., Testi et al. 2014).

## 3.2. Line Fluxes

The spectral setup was built to target simple molecules because complex organic molecules (COMs) are difficult to detect in protoplanetary disks due to the low temperatures and large partition functions of the numerous transitions. As described in Section 1, the molecules are divided into five groups: CO isotopologs and hydrogenated versions of CO ($^{12}$CO, $^{13}$CO, $C^{18}$O, $H_2$CO, $CH_3$OH), hydrocarbons ($C_2$H and $c-C_3$H$_2$), cyanides ($H^{13}$CN, $HC^{15}$N, DCN), ions ($HCO^+$, $H^{13}CO^+$), and S-bearing molecules (CS, SO). A list of the main molecular lines covered by the spectral setup is reported in Table 2.

After imaging all lines with weighting and channel widths as reported in Table 3, we employ a stacking method to boost the S/N of the lines. To do so, we use the same technique used in Teague et al. (2016), similar to the one developed by Yen et al. (2016) and Matrà et al. (2017). Exploiting the kinematical fit of the Keplerian disk obtained by Keppler et al. (2019), at every spatial location we deproject the emission and then shift the local spectrum to the systemic velocity using the GoFish package (Teague 2019). We then integrate the spectra within $3''$ from the disk center ($4''$ for the $^{12}$CO $J=2-1$ line). Figure 2 shows the integrated spectrum of the $^{12}$CO $J=2-1$ and of the $C_2$H $N=3-2$, $J=\frac{7}{2}-\frac{5}{2}$ fine-structure line as examples, wherein the latter two hyperfine components $F=4-3$ and $F=3-2$ are clearly distinguishable. To obtain line fluxes, we spectrally integrate the spatially integrated stacked spectrum between 3 and 8 km s$^{-1}$, with the systemic velocity being $5505 \pm 2$ m s$^{-1}$ (Keppler et al. 2019). The two $^{12}$CO and the $HCO^+$ lines are integrated between 1.5 and 9.5 km s$^{-1}$, as high-velocity

emission is observed in the channel maps (Section 3.3). Similarly, the $^{13}$CO, $H^{13}$CN, and $C_2$H lines are integrated between 2 and 9 km s$^{-1}$. The $C_2$H $J=\frac{7}{2}-\frac{5}{2}$ hyperfine structure lines are partially blended (see Figure 2), and fluxes have been estimated by fitting a double Gaussian to the profile. Only the $F=3-2$ hyperfine component of the $C_2$H $J=\frac{5}{2}-\frac{3}{2}$ is covered by the spectral setup. Because only half of the line is detected at the edge of the spectral window, we consider only the part of the disk within 80° of the redshifted major axis and then assume azimuthal symmetry to extract the total flux. The flux uncertainty for all lines is estimated as the standard deviation of the velocity bins not included in the velocity integration, times the square root of the number of bins over which we integrate. The line fluxes and upper limits are reported in Table 3 and shown in Figure 3.

We detect 14 lines in our Band 6 (~217–262 GHz) program, and with the two additional $^{12}$CO $J=3-2$ and $HCO^+$ lines from the Band 7 (~345–357 GHz) data, we detect a total of 16 lines, from 12 species (including isotopologs). We detect multiple transitions of $^{12}$CO, $H_2$CO, and $C_2$H. The detected lines span a large range of fluxes, with more than two orders of magnitude between the bright $^{12}$CO 3–2 line (12.3 ± 0.1 Jy km s$^{-1}$) and the faint $c-C_3H_2\ 3_{21}-2_{12}$ line (60 ± 14 mJy km s$^{-1}$).

The $c-C_3H_2\ 3_{21}-2_{12}$ and $H_2$CO $3_{21}-2_{20}$ lines are only marginally detected, with an S/N of 4.3 and 4.4, respectively. In order to increase the significance of the detections, we image the lines with an additional uv-tapering of 0.7. While this reduces the sensitivity to point sources, the downweight of long baselines can increase the sensitivity to extended structure. For $c-C_3H_2$, the resulting beam is $1.06 \times 0.92$ (PA = 104.6) with an rms of 0.9 mJy beam$^{-1}$ over one 1.38 km s$^{-1}$ channel, while for $H_2$CO, the synthesized beam is $1.03 \times 0.93$ (PA = 92.0°), with an rms of 1.2 mJy beam$^{-1}$ over one 1 km s$^{-1}$ channel. The $c-C_3H_2\ 3_{21}-2_{12}$ line shows a 5.4σ detection in the channel map at 6.45 km s$^{-1}$, with another adjacent channel showing a 4σ detection. The $H_2$CO $3_{21}-2_{20}$ presents an 8.3σ detection at 8 km s$^{-1}$, with four adjacent channels showing a signal at $\ge 4$σ. Both channel maps are included in Appendix B (Figure 12). As a final test to assess the robustness of the detection of the two lines, we divide the projected emission in concentric annuli, separated by one-quarter of the major axis of the beam. Along
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Spectra (flux density) of $^{12}$CO J = 2–1 (left) and C$_2$H N = 3–2, J = $\frac{7}{2}-\frac{5}{2}$, where the shifted and stacked spectra are integrated over an elliptical region with the disk geometrical parameters and outer radius of 4'' and 3'', respectively. The intensity maps are produced with the parameters specified in Table 3. The gray lines show the original disk spectra extracted from the same area, with the typical double horn profile apparent in the $^{12}$CO line. Two hyperfine components are clearly visible in the C$_2$H emission. The systemic velocity is 5.5 km s$^{-1}$ (Keppler et al. 2019).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Line fluxes and upper limits, portrayed as triangles. The corresponding values are reported in Table 3. Error bars represent the 1$\sigma$ intervals, while upper limits are shown at 3$\sigma$, where $\sigma$ is computed as the rms of the spectrum times the square root of the number of velocity channels used in the integration. Flux calibration uncertainties are not included. In most cases, the uncertainty is smaller than the marker size. Colors indicate the five different groups of molecules included in this study: CO isotopologs and hydrogenated versions of CO ($^{12}$CO, $^{13}$CO, $^{18}$O, H$_2$CO, CH$_3$OH), hydrocarbons (C$_2$H and C-C$_3$H$_2$), cyanides (H$^{13}$CN, HC$^{15}$N, DCN), ions (HCO$^+$, H$^{13}$CO$^+$), and S-bearing molecules (CS, SO).</figcaption></figure>

these annuli, we compute the azimuthally averaged spectrum after shifting the velocities to the systemic one by exploiting the kinematical information. This is performed using the GoFish package (Teague 2019). The azimuthally averaged spectra are shown in Figure 4, together with the azimuthally averaged spectrum of $^{12}$CO J = 2–1 for comparison. Both lines clearly show an emission peak at the systemic velocity, indicating again that both lines are robustly detected. It is worth noting that thanks to this method we manage to obtain a detection at the 50 mK level in brightness temperature units.

We do not detect SO or methanol. The three lines of sulfur monoxide are SO 6$_6$–5$_5$, 7$_6$–6$_5$, and 2$_3$–3$_2$, with upper energies of 56.5, 47.5, and 21.0 K (see Table 5), with the last being fainter by two orders of magnitude. None of the lines are detected in the channel maps or in the azimuthally averaged spectra of the natural images with and without uv-tapering. One line of $^{34}$SO is also included in the spectral setup at 246.6636 GHz, but is not detected, as expected from the nondetection of the main isotopolog. The two brightest methanol lines included in the spectral setup are both for A-type states (analog to ortho-ammonia
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Azimuthally averaged spectra of three molecular lines, showing the detection of H₂CO 3₂₁–2₂₀ and c-C₃H₂ 3₂₁–2₁₂, and the ¹²CO J = 2–1 line as comparison. The two weak lines are imaged with natural weighting, and an additional uv-tapering to increase sensitivity to extended emission. The 2D spectra are shown in brightness temperature, with the conversion performed under the Rayleigh–Jeans approximation. The spectra are computed after deprojecting and shifting the velocities using the kinematical and geometrical model by Keppler et al. (2019) and then averaging the spectra along the disk annuli. The H₂CO and c-C₃H₂ lines have been supersampled by a factor of 4 in velocity. Dashed white contours show the 3 K level of the ¹²CO emission.</figcaption></figure>

in the spin state of the methyl group). The transitions are 10<sub>2</sub>–9<sub>3</sub>(A<sup>+</sup>) at 232.4185 GHz and 5<sub>1</sub>–4<sub>1</sub>(A<sup>−</sup>) at 243.9157 GHz. Neither of the two lines are detected with the methods described above. The last notable nondetection is for the c-C<sub>3</sub>H<sub>2</sub> 7<sub>16</sub>–7<sub>07</sub> at 218.7327 GHz. Finally, a matched filter analysis is performed on the whole spectral setup in the uv-plane, with particular attention given to the wide basebands. To do so, we exploit the VISIBLE package by Loomis et al. (2018a), which has been successfully applied to protoplanetary disk observations (e.g., Loomis et al. 2018b, 2020). This method requires a model image as matched filter, and we thus run the analysis onto the continuum subtracted visibilities with the CS and H<sup>f</sup>CO<sup>+</sup> channel maps as filters. No new additional line is detected with this method.

### 3.3. Integrated Maps

Integrated intensity maps of the 14 detected transitions are derived with the same imaging parameters used in Section 3.2. The c-C<sub>3</sub>H<sub>2</sub> 3<sub>21</sub>–2<sub>12</sub> line is imaged with natural weighting and an additional uv-tapering, as specified in Section 3.2, to increase brightness temperature sensitivity. Integrated intensity maps of the 14 detected lines are then created integrating over the whole cubes with a 2.5σ clipping. In order to increase the contrast of the images for faint lines, we adopt a Keplerian masking technique, by masking any signal that does not originate from the Keplerian mask used in the cleaning (see Section 2). This technique has already been exploited by different groups to boost the S/N in their data (e.g., Salinas et al. 2017; Ansdell et al. 2018; Bergner et al. 2018). The cases where the image quality improves significantly are CS and H₂CO 3₂₁–2₂₀. The emission morphology does not vary significantly between the different techniques for the other lines.

The integrated intensity maps are shown in Figure 5, where the molecular emission is compared to the high-resolution

855 µm continuum emission from Isella et al. (2019; cleaned with Briggs weighting, robust = 0.3). All lines show a bright ring of emission outside the orbit of the two massive planets, which are highlighted as white dots in the figure. While the general spatial morphology is similar among different lines, they clearly exhibit different structures, together with very different intensities. For example, hydrogen-cyanide isotopologs and H<sup>13</sup>CO<sup>+</sup> show a very defined narrow ring similarly to the continuum emission, while H₂CO and C₂H exhibit more extended emission. We analyze the radial profiles of the emission lines in Section 3.4. Most lines seem to show very azimuthally symmetric emission, with the only possible exception of H<sup>13</sup>CO<sup>+</sup>.

The high-resolution images of the Band 7 data of ¹²CO J = 3–2 and HCO<sup>+</sup>J = 4–3 show an inner disk, inside the deep gap carved by the planets (Long et al. 2018; Keppler et al. 2019). Similarly to the Band 6 continuum emission, it is difficult to distinguish any inner disk in the other molecular tracers, due to the lower angular resolution. However, we can explore the contribution of molecular gas within the orbits of the planets in the high-velocity channels, associated with emission from small radii. From Wang et al. (2021), planet c has an orbit with semimajor axis of 33.2 au. Using the geometrical parameters by Keppler et al. (2019), in particular the disk inclination of 51°7 and a stellar mass of 0.87 M<sub>⊕</sub>, any Keplerian gas with projected velocity higher than 3.75 km s<sup>-1</sup> from the systemic velocity must emit from within the orbit of PDS 70 c. With a systemic velocity of 5.5 km s<sup>-1</sup>, this includes emission at <1.75 km s<sup>-1</sup> and >9.25 km s<sup>-1</sup>. The ¹²CO J = 2–1 line shows emission at high velocities (at ~0 and ~10 km s<sup>-1</sup>; see Figure 2), indicating a contribution from small radii, as is the case for the ¹²CO J = 3–2 and the HCO<sup>+</sup>J = 4–3 lines. Emission lines of C<sub>2</sub>H, <sup>13</sup>CO, and H<sup>13</sup>CN start to show detectable emission at ~1.85 km s<sup>-1</sup>, thus indicating that no
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Integrated intensity maps of the 14 detected lines and of the 885 µm continuum from Isella et al. (2019) in the bottom right in comparison. The imaging parameters of all lines are listed in Table 4. The CS and H₂CO 3₂₁–2₂₀ channel maps have been collapsed using a Keplerian mask, while all others have been produced with a 2.5σ clip.</figcaption></figure>
emission is likely seen within the orbit of the outermost detected planet at the sensitivity of our observations.

### 3.4. Radial Profiles

In order to extract the radial profiles of the integrated intensities, we azimuthally average the deprojected integrated intensity maps shown in Figure 5. Error bars are estimated from the standard deviations of the emission within the annuli, divided by the square root of the number of independent beams across the azimuthal range. The width of the annuli is a quarter of the beam major axis. The resulting integrated intensity profiles are shown in Figure 6.

A deep gap colocated with the orbital radius of PDS 70b is observed in HCO$^+$ and $^{12}$CO $J=3-2$ (Long et al. 2018; Keppler et al. 2019). For all other lines, no emission is seen from an inner disk in the integrated intensity profiles, in line with the discussion of Section 3.3. The brightest lines clearly exhibit emission at radii $>2''$, in particular H$_2$CO, $^{12}$CO, $^{13}$CO, and C$_2$H. The $^{12}$CO $J=2-1$ line, which has a higher sensitivity than the $J=3-2$ transition, shows emission out to radii $>3''$. The H$_2$CO line shows a shoulder of emission outside the continuum ring between $1''$ and $1''5$, which is also hinted in the DCN line. This may be due to enhanced gas temperature outside the dust ring where photons penetration depths are higher due to the lower dust column densities (e.g., Öberg et al. 2015; Cleeves 2016; Facchini et al. 2017; Öberg et al. 2021). This would induce a stronger thermal and photodesorption of CO, leading to the H$_2$CO shoulder.

## 4. Derivation of Disk Physical Parameters

In the following section, disk fluxes and integrated intensity radial profiles are exploited to derive physical and chemical parameters of the PDS 70 disk, mainly column densities and deuteration and fractionation profiles, without requiring any chemical modeling.

### 4.1. Column Densities and Optical Depths

In the optically thin limit, the column density of an emitting molecule, $N$, can be related to the integrated flux density, $S_νΔν$, and excitation temperature, $T_{ex}$, via

$$N_{\text{thin}} = \frac{4\pi}{A_{\text{ul}}} \frac{S_{\nu} \Delta\nu}{\Omega} \frac{Q(T_{\text{ex}})}{g_{\text{u}}} \exp\left(\frac{E_{\text{u}}}{T_{\text{ex}}}\right), \quad (1)$$

where $Ω$ is the area over which the flux was integrated, $A_{ul}$ is the Einstein coefficient of the transition, $E_u$ is the upper-state energy, $h$ is the Planck constant, $c$ the speed of light in vacuum, $g_u$ the upper-state degeneracy, and $Q$ the partition function. When the optical depth, $τ$, is no longer negligible, this can be corrected by

$$N = \frac{N_{\text{thin}} \tau}{1 - \exp(-\tau)} = N_{\text{thin}} C_{\tau}, \quad (2)$$

following, for example, Goldsmith & Langer (1999). In turn, the optical depth at the line center (assuming a top-hat line profile) can be calculated via

$$\tau = \frac{A_{\text{ul}} c^3}{8 \pi \nu^3 \Delta \nu} N g_{\text{u}} \frac{\exp(-E_{\text{u}}/T_{\text{ex}})}{Q(T_{\text{ex}})} \left[ \exp\left(\frac{h\nu}{kT_{\text{ex}}}\right) - 1 \right], \quad (3)$$

where $ν$ is the frequency of the transition. In the following, we use the molecular data provided in Table 5, taken from the CDMS catalog (Müller et al. 2001), and interpolate the partition function, $Q$, with a cubic spline.

### 4.1.1. H$_2$CO

As we detected two transitions of H$_2$CO, $3_{03}-2_{02}$ and $3_{21}-2_{20}$, we are able to constrain both $N$ and $T_{ex}$ simultaneously (e.g., Salinas et al. 2017; Bergner et al. 2019; Öberg et al. 2021, for similar analyses). For the integrated fluxes, we radially bin the integrated intensity profiles shown in Figure 6 into $0''1$ wide annuli to measure $S_νΔν/Ω$.

Following, e.g., Loomis et al. (2018b), we use emcee (Foreman-Mackey et al. 2013) to explore the posterior distributions for {$N(H_2CO), T_{ex}$} while self-consistently solving for the optical depth of the two transitions. For each radial bin, we use 128 walkers, each taking 1000 steps to burn in, and then an additional 500 steps to sample the posterior distributions. As the posteriors were found to be Gaussian, we take the median value for each parameter as the most likely value, and the 16th to 84th percentile as a quantification of the uncertainty on this. The resulting radial profiles of $N(H_2CO)$ and $T_{ex}$ are found in Figure 7.

The excitation temperature ranges between ~10 and 40 K, with the temperature decreasing in the outer disk regions. Higher temperatures are found for the inner regions, as expected for the cavity wall being strongly irradiated; however, the $S/N$ of the data is low, as reflected by the large inferred uncertainties. Column densities across the $0''5-2''$ interval range between $10^{13}-10^{14} \text{cm}^{-2}$, in line with other protoplanetary disks (e.g., Pegues et al. 2020; Terwisscha van Scheltinga et al. 2021).

### 4.1.2. C$_2$H

Three hyperfine transitions of C$_2$H were detected; however, the two $3-2$ $J = \frac{7}{2} - \frac{5}{2}$ transitions were sufficiently blended that typical hyperfine-fitting methods, e.g., Hily-Blant et al. (2013b, 2017) and Teague & Loomis (2020), could not be used. Similarly, there was only a single $3-2$ $J = \frac{5}{2} - \frac{3}{2}$ fine-structure group component detected, precluding typical hyperfine analyses. Instead we follow a similar approach to that used for H$_2$CO, but for each individual component separately, which implicitly assumes that local thermodynamic equilibrium (LTE) is valid.

As the $J = \frac{7}{2} - \frac{5}{2}$ transitions were blended, we used GoFish to calculate azimuthally averaged, velocity-corrected spectra for annuli spaced $0''1$ apart (as in the right panel of Figure 3). To each of these spectra, we fit two Gaussian components using the scipy.optimize.curve_fit routine, fixing the line centers to the rest frequency of the two transitions, and assuming that the lines are optically thin such that the combined spectrum is a simple sum of Gaussian components. For the single $J = \frac{5}{2} - \frac{3}{2}$ transition, only part of the disk emission was covered (see discussion in Section 3.2). Thus, the azimuthally averaged, velocity-corrected spectra only considered emission within $80^\circ$ of the redshifted major axis. To this spectrum, a single Gaussian component was fit in order to measure the line flux.

As only a single transition is considered, $T_{ex}$ was held constant at values of {20, 30, 50} K, guided by previous analyses of C$_2$H
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Azimuthally averaged radial profiles of the 14 detected lines, and of the 885 µm continuum by Isella et al. (2019) in the bottom right in comparison. The imaging parameters of the respective channel maps are listed in Table 4, with integrated map images shown in Figure 5. The ribbon shows the 1σ rms, excluding the uncertainty from absolute flux calibration. The gray line at the top right of all panels shows the beam major axis. The brightness temperature conversion was done under the Rayleigh-Jeans approximation.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Left: H₂CO excitation temperature derived from rotational diagram analysis of the 3<sub>03</sub>−2<sub>02</sub> and 3<sub>21</sub>−2<sub>20</sub> transitions. The gray line at the top right shows the beam major axis. Right: H₂CO column density. Integrated intensity maps with natural weighting were used in the analysis to allow for enough S/N in the weaker 3<sub>21</sub>−2<sub>20</sub> transition.</figcaption></figure>

in other sources (e.g., Bergner et al. 2019) in order to explore
likely ranges of N(C₂H). Again, emcee was used to sample the
posterior distribution of N(C₂H), with τ being calculated self-
consistently. The 128 walkers took 1000 steps to burn in and an
additional 500 to sample the posteriors. The resulting radial
profiles of N(C₂H) are shown for the three hyperfine components
in Figure 8, where the solid lines show the result for T<sub>ex</sub> = 30 K,
and the shaded regions show the ranges when considering lower
and higher T<sub>ex</sub> values. The concordance between these three
transitions suggests that the C₂H column density is well
constrained and that LTE is a reasonable assumption.

4.1.3. Other Molecules

We follow a similar procedure to that used for C₂H for all
other molecules that only have a single transition detected. The
integrated intensities are calculated from the azimuthally
averaged profiles of the total intensity maps shown in
Figure 5 and binned into annuli of 0.1′′1. As with C₂H, three
T<sub>ex</sub> values were considered: 20, 30, and 50 K. We add a prior
on the optical depth to be <1 for all molecules, because all of
the considered transitions show a brightness temperature that is
well below the assumed excitation temperature range. The
inferred column density profiles and optical depths are shown
in Figure 9.

4.1.4. *c-C<sub>3</sub>H<sub>2</sub> and Upper Limits on SO and CH<sub>3</sub>OH*

The detected c-C₃H₂ line exhibits a low S/N, and we thus
compute a disk-averaged column density by using Equation (1)
and assuming optically thin emission. A radius of 3'' is used to
compute Ω. A column density of 5.7<sup>+4.0</sup><sub>-1.6</sub> × 10<sup>12</sup> cm<sup>-2</sup> is
obtained, including the uncertainty from the unknown excita-
tion temperature (ranged again between 20 and 50 K). The
relatively high column, associated with the weak molecular
emission, is easily explained by the large partition function.

From the upper limits on SO and CH3OH, we can derive
upper limits on the disk-averaged column density following the
same approach. The SO 76–65 line is the most stringent of the
two considered in this work, with a 3σ upper limit on the SO
disk-averaged column density of 1.5 × 10¹¹ cm⁻², where the
uncertainty of the excitation temperature has been considered
as above. As for CH3OH, the 51–41(A⁻) is the most
constraining transition, with a relative upper limit on the
column density of 6.5 × 10¹² cm⁻². The 102–93(A⁺) transition
leads to much weaker constraint, due to the high upper-level
energy and lower Einstein coefficient.

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Derived column densities and optical depths for the three hyperfine C₂H lines included in the spectral setup. The assumed excitation temperature Tₐₓ is varied between 20 and 50 K, with the solid lines showing the 30 K case. The ribbons include the scatter in the obtained column densities from the posterior distributions between the 16th and 84th percentiles after accounting for the 10% flux calibration uncertainty and the systematic uncertainty driven by the unknown Tₐₓ. The hatched area at r < 0.1''/2 indicates the region where beam smearing dominates the column density reconstruction.</figcaption></figure>

4.2. Deuteration and N Fractionation

The hydrogen-cyanide isotopologs can be used to infer the
level of deuteration and nitrogen fractionation in PDS 70. The
H<sup>13</sup>CN/DCN and H<sup>13</sup>CN/HC<sup>15</sup>N ratios have been already
used in the literature to extract disk-average deuteration and
fractionation values in a sample of disks comprising about 10
sources (e.g., Huang et al. 2017; Guzmá et al. 2017; Bergner
et al. 2019, 2020). We follow the same procedure used in
Bergner et al. (2020) to extract the column density ratio
between the three isotopologs. In order to derive HCN column
densities, a <sup>12</sup>C/<sup>13</sup>C ratio of 68 ± 15 is assumed (Milam et al.
2005). While the carbon fractionation is still not well
established in the majority of disks, Hily-Blant et al. (2019)
derived a uniform carbon fractionation across the TW Hya disk
using HCN and H<sup>13</sup>CN lines, with a value compatible with the
one we assume. The column density ratio between two
isotopologs is derived using Equation (1) and taking the ratio
between two different isotopologs. As noticed by Guzmá et al.
(2017), the similar upper-state energy and Q/g<sub>u</sub> ratios of the
three isotopologs lead to a column density ratio that is very
similar to a flux ratio. We assume that the emission from the
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Column densities and optical depths of all detected lines, excluding C<sub>2</sub>H and H<sub>2</sub>CO, which are reported in Figures 7–8, and c-C<sub>3</sub>H<sub>2</sub>, which presents a very weak emission line. The method used to derive both column densities and optical depths is described in Section 4.1. The ribbon indicates statistical uncertainties including a 10% flux calibration uncertainty and systematic uncertainties driven by the unknown excitation temperatures, which have been varied between 20 and 50 K. The solid lines represent the profiles for T<sub>ex</sub> = 30 K.</figcaption></figure>

three isotopologs is at the same excitation temperature, which we fix at 30 K. As in Section 4.1, the temperature is varied between 20 and 50 K in the computation of the uncertainties on the ratio. The systematics on the excitation temperature leads to relatively small uncertainties in the ratio, between 7% for the H/D ratio and 0.1% for the <sup>14</sup>N/<sup>15</sup>N ratio (where we assumed the molecules to be at the same temperature in the uncertainty calculation). Other sources of uncertainty are included in the error estimate (flux rms, <sup>12</sup>C/<sup>13</sup>C, and absolute flux calibration), and in all cases dominate the uncertainties over the excitation temperature systematics. The absolute flux

calibration uncertainty is not included for H<sup>13</sup>CN/HC<sup>15</sup>N, because the two transitions are very close in frequency in the spectral setup (we ignore any uncertainty in the bandpass calibration). The derived disk-averaged H/D and <sup>14</sup>N/<sup>15</sup>N in HCN isotopologs are 60 ± 17 and 111 ± 27, respectively, where a radius of 3" has been used to extract the disk-averaged values.

The two values show a significant D and <sup>15</sup>N enrichment in HCN isotopologs relative to protosolar elemental ratios. In Jupiter atmosphere, D/H ~ 2.5 × 10<sup>-5</sup> from the HD/H<sub>2</sub> ratio (Linsky et al. 2006), while solar winds indicate a ratio <sup>14</sup>N/<sup>15</sup>N
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Left: H/D ratio calculated as N(HCN)/N(DCN), where N(HCN) has been computed from the H<sup>13</sup>CN column density. Right: <sup>14</sup>N/<sup>15</sup>N ratio calculated as N(HCN)/N(HC<sup>15</sup>N). The gray lines at the top right show the beam major axis, after smoothing the images used in the ratios to the same beam. The ribbons indicate 1σ uncertainties, which include the standard deviation and rms divided by the number of beams along the azimuthal ring (left), the uncertainty on the <sup>12</sup>C/<sup>13</sup>C ratio, on the excitation temperature, and for the N(HCN)/N(DCN) ratio on the absolute flux due to calibration uncertainties.</figcaption></figure>

~ 441 ± 5 (Marty et al. 2011). Fractionation levels of HCN isotopologs in both nitrogen and hydrogen show consistently lower values in protoplanetary disks. In particular, Huang et al. (2017) showed a H/D ranging between ~15 and 200 in a sample of six disks, based on DCN and H<sup>13</sup>CN lines. Guzmá et al. (2017) obtained a <sup>14</sup>N/<sup>15</sup>N between ~80 and 160 in the same sample by analyzing H<sup>13</sup>CN and HC<sup>15</sup>H transitions. PDS 70 falls well within these distributions. The nitrogen fractionation of HCN in PDS 70 is in line with values observed in solar system bodies, e.g., comets showing HCN/HC<sup>15</sup>N ≈ 100–250 (Mumma & Charnley 2011). Such low ratios are comparable with HCN isotopolog observations of prestellar cores (e.g., Hily-Blant et al. 2013a), even though new observations suggest that there may be a large variation between different cores (Magalhães et al. 2018), in particular in the high-mass regime (Colzi et al. 2018). Thus, the high nitrogen fractionation inferred in PDS 70 may be inherited from its progenitor core.

We derive the spatial dependence of H/D and $^{14}$N/$^{15}$N using column densities of the three HCN isotopologs. Purely inherited material would show uniform fractionation and deuteration across the radial extent of the disk (e.g., Guzmá et al. 2017; Hily-Blant et al. 2019), while radial gradients would be indicative of continued chemical processing. The radial profiles of the column densities are computed as described in Section 4.1 from integrated intensity maps created with matching synthesized beams. In particular, the H$^{13}$CN and DCN lines are compared with maps obtained with a robust = 0 parameter as in Section 3.3 and then smoothed to a common circular beam of 0.42 with the CASA *imsmooth* task. The H$^{13}$CN and HC$^{15}$N lines are instead compared from images with natural weighting, which lead to the same beam. An excitation temperature of 30 K is assumed, and uncertainties on T<sub>ex</sub> are considered in the error propagation. Other sources of uncertainties are considered, including the uncertainty on $^{12}$C/$^{13}$C, as for the disk-average values. The column density ratios are shown in Figure 10. The figures include the derived radial profiles of N(HCN)/N(DCN) and N(HCN)/N(HC$^{15}$N), with N(HCN) being 68 × N(H$^{13}$CN). The analysis is limited to r ≲ 1.8 due to the low S/N at larger radii.

Both D/H and $^{14}$N/$^{15}$N profiles show some radial dependence, with enhanced deuteration and nitrogen fractionation in

outer regions of the disk. To our knowledge, this is the first time that a radial profile of the deuteration of HCN has been measured in a protoplanetary disk. Huang et al. (2017) suggest that LkCa 15 may present a similar behavior, with enhanced deuteration in the outer regions of the disk, outside the millimeter ring of this transition disk. Because the formation of DCN is mediated by CH<sub>2</sub>D<sup>+</sup> (e.g., Millar et al. 1989; Turner 2001), a high level of deuterium fractionation may indicate favorable deuterium exchange between HD and C<sub>2</sub>H<sub>2</sub><sup>+</sup> at temperatures <80 K outside the pebble disk at lower temperatures (see Section 3.4). Finally, the HCN deuteration profile indicates that the dust ring at 0.7–0.8 may affect the deuteration process, with the N(H$^{13}$CN)/N(DCN) profile showing a peak colocated with the dust millimeter thermal emission. This coincides with the temperature increase in the H<sub>2</sub>CO excitation temperature profile (see Figure 7), thus directly tracing a less efficient deuterium fractionation at higher temperatures in the disk, as expected from chemical models. The correlation of the N(H$^{13}$CN)/N(DCN) peak with the submillimeter continuum ring highlights the direct impact that dust substructure can have on chemical processes in protoplanetary disks, including isotope exchange reactions.

The radial profile of the nitrogen fractionation is particularly surprising. Guzmán et al. (2017) and Hily-Blant et al. (2019) observed radially decreasing fractionation levels in HCN isotopologs in the V4046 Sgr and TW Hya disks, respectively. Hily-Blant et al. (2019) observed tentative evidence of increasing fractionation levels in the very outer disk, but the low S/N of the observations precluded a robust conclusion. A radially decreasing level of fractionation can be explained by enhanced fractionation in the inner regions of the disk due to isotope selective photodissociation of the self-shielding molecule N<sub>2</sub> (Heays et al. 2014; Visser et al. 2018). Higher UV fluxes in the inner regions can selectively photodissociate N<sup>15</sup>N, with the <sup>15</sup>N atom free to enter gas-phase chemical reactions and lead to formation of HC<sup>15</sup>N. However, the radially decreasing fractionation in PDS 70 suggests that in the outer disk the main fractionation mechanism is low-temperature isotope exchange reactions (e.g., Terzieva & Herbst 2000; Roueff et al. 2015; Wiström & Charnley 2018), which dominates over the isotope selective photodissociation pathway. It is not clear why the behavior in PDS 70 is so different
from the TW Hya and V4046 Sgr disks. The most likely
possibility is that the inner cavity and the bright dust wall at the
cavity edge significantly affects the temperature and radiation
structure across the whole disk, thus leading to favorable low-
temperature reactions outside the cavity wall. The radial
polarized intensity profile in the J band (see Figure 2 in
Keppler et al. 2018 and Section 5.2 in this paper) shows a
pronounced peak at ~0.48 from the central star, suggesting
that the outer regions may be shadowed by the warm wall of
material at that cavity edge, as predicted by theoretical models
(Isella & Turner 2018).

As specified earlier, the radial profiles assume a uniform $N(HCN)/N(H^{13}CN)$ ratio. However, Visser et al. (2018) noted that from thermochemical models, high $N(HCN)/N(H^{13}CN)$ are found at large radii ($r > 150$ au in their model). We thus caution that our conclusions about deuteration and nitrogen fractionation may be affected by the carbon fractionation profile, were this not uniform (e.g., Colzi et al. 2020). Observations of the main HCN isotopolog are needed to break this degeneracy (e.g., Hily-Blant et al. 2019).

Our observations indicate that organic material in the disk of PDS 70 possesses high nitrogen and hydrogen fractionation levels. Thus, organic material accreted by the two planets will possess high fractionation levels, which may alter the isotope ratios of the planetary atmospheres.

# 5. Discussion

## 5.1. C/O Ratio

The elemental carbon-to-oxygen (C/O) ratio in the atmosphere of giant planets has been proposed as a probe of the accretion history of massive protoplanets during their coevolution with their host protoplanetary disk (e.g., Pudritz et al. 2018; Lammer & Blanc 2018; Madhusudhan 2019; Öberg & Bergin 2021; van Dishoeck & Bergin 2020, for some recent reviews). While molecular phase transitions at specific temperatures across the disk can affect the elemental ratio of gas-phase material (e.g., Öberg et al. 2011), additional chemical evolution within the disk can alter the C/O ratio and thus the chemical outcome of the planet formation process (e.g., Eistrup et al. 2016, 2018; Cridland et al. 2019, 2020). Probing the C/O ratio in disks is thus key in accessing whether the combination of both chemical and physical evolution of gas and dust in disks has affected the elemental abundances available in the disk gas-phase material (e.g., Kama et al. 2016; Zhang et al. 2019; Krijt et al. 2020).

To determine the C/O ratio in the vertical layers of disks that can be probed by ALMA observations, emission lines of small hydrocarbons are extremely valuable. In particular, the abundance of C₂H increases very rapidly as soon as C/O > 1, where free C atoms are not locked up in CO, thereby quenching complex organic chemistry. For this reason, the large C₂H fluxes of rotational transitions have been used to infer C/O ratios above the solar 0.54 value (in some cases above unity) in a few protoplanetary disks (e.g., Bergin et al. 2016; Cleeves et al. 2018; Miotello et al. 2019).

Comparing the C₂H luminosities of disks from the surveys by Bergner et al. (2019) and Miotello et al. (2019), PDS 70 is on the high flux range of the Lupus survey (accounting for source distances; see Figure 4 in Miotello et al. 2019), with a total flux of the $J = \frac{7}{2} - \frac{5}{2}$ transition of ~1327 ± 18 mJy km s⁻¹ (see Table 3). In the thermochemical models by Miotello et al. (2019), such a

high flux cannot be obtained with a solar C/O ratio, and elemental
ratios close to ~1.5 are required. Similarly, Fedele & Favre (2020)
proposed to normalize line fluxes to the ¹³CO J = 2–1 line to ease
comparisons between sources. By exploring part of the wide
parameter space that disks can exhibit in their physical and stellar
properties, they showed that C/O > 1 is required when the ratio
of the C₂H N = 3–2 J = $\frac{7}{2} - \frac{5}{2}$ F = 4–3 line to the ¹³CO J = 2–1
line is $\gtrsim$0.2. For PDS 70, this ratio ranges between 0.25–0.35
throughout the whole disk out to r ~ 2'', indicating again a C/O
ratio above unity for this particular disk. However, we note that all
these models have been fine-tuned for disks that do not exhibit a
large cavity, as is the case for PDS 70. Accurate C/O ratio
retrievals need a forward model of the PDS 70 disk, which will be
the focus of future works.

An alternative way to qualitatively determine the C/O ratio
in PDS 70 is by looking at the N(CS)/N(SO) ratio, as suggested
by Semenov et al. (2018) for the DM Tau disk (see also Dutrey
et al. 2011; Guilloteau et al. 2016 for a larger disk sample).
With a grid of thermochemical models, the authors show that
the only models able to reproduce their observations of sulfur-
bearing molecules (in particular CS and SO) are with a C/O
ratio $\gtrsim$1. Their reference model produces a disk-averaged
column density ratio N(CS)/N(SO)~0.03, while observations
indicate a ratio with a lower limit >1. Other chemical and
physical parameters can affect the column density ratio, but not
as significantly as the C/O ratio.

We apply the same qualitative argument to the PDS 70 disk. By using the $3\sigma$ upper limit on the SO column density and the disk-averaged CS column density of $1.8_{0.2}^{+0.4} \times 10^{12} \text{ cm}^{-2}$, we obtain a lower limit on the N(CS)/N(SO) ratio of ~109. This value is 4000 times higher than what Semenov et al. (2018) predicted for DM Tau with a solar C/O ratio. A C/O ratio of ~1 leads to an N(CS)/N(SO) ratio of ~10 Semenov et al. 2018. Our lower limit therefore strongly suggests that the upper molecular layers of PDS 70 have a C/O ratio above unity, with a significant amount of oxygen being sequestered in likely refractory form. Thermochemical models of PDS 70 are needed to extract a more quantitative conclusion.

In summary, a comparison of the PDS 70 inventory with
thermochemical models of protoplanetary disks from the
literature suggests a disk-average C/O ratio >1 using different
proxies, making this qualitative conclusion robust against
different methods.

## 5.2. Comparing Millimeter Continuum, IR Scattered Light, and Molecular Emission Profiles

In Section 3.4, we have shown that all molecular transitions
well detected in the ALMA data exhibit a ring-like morphology, with a high level of azimuthal symmetry. Molecular emission is observed to emit in ring-like morphologies in many
disks and with different tracers (e.g., Öberg et al. 2015; Teague
et al. 2016; Bergin et al. 2016), with the emission morphology
being primarily caused by chemical and excitation effects (e.g.,
Öberg et al. 2015; Cazzoletti et al. 2018; Long et al. 2021),
rather than the underlying dust density structure. However, the
gas structure can lead to ring-like morphologies in molecular
transitions more sensitive to gas densities, such as CO
isotopologs, as is clearly the case in extreme conditions as
for transition disks (e.g., van der Marel et al. 2016). The
presence of gaps and rings in dust is also expected to lead to
molecular abundance variations along the radial coordinate
from thermochemical models (e.g., Facchini et al. 2018;
<table><caption>**Table 4**

Images Used to Extract the Radial Profiles of the 14 Detected Lines</caption><thead><tr><th>Species</th><th>Transition</th><th>uv-tapering</th><th>Beam (PA)</th><th>Radius of Peak (")</th><th>Error<sup>a</sup> (")</th><th>Peak Integrated Intensity<sup>b</sup> (K km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td><sup>12</sup>CO</td><td>3–2</td><td>n</td><td>0.11 × 0.10 (90°3)</td><td>0.40</td><td>0.03</td><td>57.1</td></tr><tr><td><sup>12</sup>CO</td><td>2–1</td><td>n</td><td>0.37 × 0.33 (67°1)</td><td>0.42</td><td>0.09</td><td>46.0</td></tr><tr><td><sup>13</sup>CO</td><td>2–1</td><td>n</td><td>0.39 × 0.35 (69°6)</td><td>0.54</td><td>0.10</td><td>17.0</td></tr><tr><td>C<sup>18</sup>O</td><td>2–1</td><td>n</td><td>0.39 × 0.35 (68°3)</td><td>0.64</td><td>0.10</td><td>5.6</td></tr><tr><td>H<sub>2</sub>CO</td><td>3<sub>03</sub>–2<sub>02</sub></td><td>n</td><td>0.40 × 0.35 (69°0)</td><td>0.55</td><td>0.10</td><td>5.4</td></tr><tr><td>H<sub>2</sub>CO</td><td>3<sub>21</sub>–2<sub>20</sub></td><td>n</td><td>0.57 × 0.52 (100°8)</td><td>0.21</td><td>0.14</td><td>0.4</td></tr><tr><td>C<sub>2</sub>H<sup>c</sup></td><td>3–2 J = <sup>7</sup>&frasl;<sub>2</sub> <sup>5</sup>&frasl;<sub>2</sub></td><td>n</td><td>0.39 × 0.32 (107°7)</td><td>0.63</td><td>0.10</td><td>6.8</td></tr><tr><td>c-C<sub>3</sub>H<sub>2</sub></td><td>3<sub>21</sub>–2<sub>12</sub></td><td>y</td><td>1.06 × 0.92 (104°6)</td><td>0.93</td><td>0.27</td><td>0.1</td></tr><tr><td>H<sup>13</sup>CN</td><td>3–2</td><td>n</td><td>0.39 × 0.32 (107°1)</td><td>0.44</td><td>0.10</td><td>4.4</td></tr><tr><td>HC<sup>15</sup>N</td><td>3–2</td><td>n</td><td>0.56 × 0.46 (100°6)</td><td>0.49</td><td>0.14</td><td>1.1</td></tr><tr><td>DCN</td><td>3–2</td><td>n</td><td>0.40 × 0.35 (70°8)</td><td>0.45</td><td>0.10</td><td>2.7</td></tr><tr><td>HCO<sup>+</sup></td><td>4–3</td><td>n</td><td>0.21 × 0.18 (73°1)</td><td>0.45</td><td>0.05</td><td>23.4</td></tr><tr><td>H<sup>13</sup>CO<sup>+</sup></td><td>3–2</td><td>n</td><td>0.39 × 0.32 (108°8)</td><td>0.54</td><td>0.10</td><td>3.3</td></tr><tr><td>CS</td><td>5–4</td><td>n</td><td>0.42 × 0.34 (108°8)</td><td>1.09</td><td>0.10</td><td>2.7</td></tr></tbody></table>

**Notes.** The radius of the peak indicates the radius at which the radial profile of the integrated intensities has a maximum at $r > 0.2''$.

<sup>a</sup> Error is estimated as one-quarter of the major axis of the clean beam.

<sup>b</sup> Integrated intensity at the peak radius for $r > 0.2''$ in K km s<sup>-1</sup>, converted from intensity units assuming Rayleigh–Jeans regime.

<sup>c</sup> The two hyperfine components have been considered simultaneously.

<table><caption>**Table 5**

Molecular Line Data and Peak Column Densities and Optical Depths Derived in Section 4.1</caption><thead><tr><th>Species</th><th>Transition</th><th>E<sub>u</sub> (K)</th><th>log A<sub>ul</sub> (s<sup>-1</sup>)</th><th>g<sub>u</sub></th><th>Q(18.75 K)</th><th>Q(37.50 K)</th><th>Q(75 K)</th><th>N<sub>max</sub><sup>a</sup> (cm<sup>-2</sup>)</th><th>τ<sub>max</sub><sup>a</sup></th></tr></thead><tbody><tr><td>C<sup>18</sup>O</td><td>2–1</td><td>15.8</td><td>-6.2</td><td>5</td><td>7.46</td><td>14.57</td><td>28.81</td><td>6.2<sup>+0.2</sup><sub>-2.8</sub> × 10<sup>15</sup></td><td>9.8<sup>+0.02</sup><sub>-0.41</sub> × 10<sup>-1</sup></td></tr><tr><td>H<sub>2</sub>CO</td><td>3<sub>03</sub>–2<sub>02</sub></td><td>21.0</td><td>-3.55</td><td>7</td><td>44.68</td><td>128.65</td><td>361.72</td><td>8.4<sup>+1.8</sup><sub>-1.4</sub> × 10<sup>13</sup></td><td>8.8<sup>+11.3</sup><sub>-5.8</sub> × 10<sup>-1</sup></td></tr><tr><td>H<sub>2</sub>CO</td><td>3<sub>21</sub>–2<sub>20</sub></td><td>68.1</td><td>-3.80</td><td>7</td><td>44.68</td><td>128.65</td><td>361.72</td><td>8.4<sup>+1.8</sup><sub>-1.4</sub> × 10<sup>13</sup></td><td>1.7<sup>+0.3</sup><sub>-0.3</sub> × 10<sup>-1</sup></td></tr><tr><td>CH<sub>3</sub>OH</td><td>5<sub>1</sub> – 4<sub>1</sub>(A<sup>-</sup>)</td><td>49.7</td><td>-3.85</td><td>11</td><td>68.75</td><td>230.24</td><td>731.07</td><td>&lt;6.5 × 10<sup>12</sup></td><td>...</sup>c</td></tr><tr><td>CH<sub>3</sub>OH</td><td>10<sub>2</sub> – 9<sub>3</sub>(A<sup>+</sup>)</td><td>165.3</td><td>-4.22</td><td>21</td><td>68.75</td><td>230.24</td><td>731.07</td><td>&lt;1.2 × 10<sup>15</sup></td><td>...</sup>c</td></tr><tr><td>C<sub>2</sub>H</td><td>J = 7/2–5/2 F = 4–3</td><td>25.1</td><td>-4.12</td><td>9</td><td>37.14</td><td>72.91</td><td>144.71</td><td>2.0<sup>+0.5</sup><sub>-0.1</sub> × 10<sup>14</sup></td><td>3.5<sup>+0.2</sup><sub>-1.1</sub> × 10<sup>-2</sup></td></tr><tr><td>C<sub>2</sub>H</td><td>J = 7/2–5/2 F = 3–2</td><td>25.1</td><td>-4.13</td><td>7</td><td>37.14</td><td>72.91</td><td>144.71</td><td>3.4<sup>+2.0</sup><sub>-1.5</sub> × 10<sup>14</sup></td><td>7.4<sup>+1.7</sup><sub>-1.8</sub> × 10<sup>-2</sup></td></tr><tr><td>C<sub>2</sub>H</td><td>J = 5/2–3/2 F = 3–2</td><td>25.2</td><td>-4.15</td><td>7</td><td>37.14</td><td>72.91</td><td>144.71</td><td>1.8<sup>+0.3</sup><sub>-0.1</sub> × 10<sup>14</sup></td><td>3.8<sup>+2.2</sup><sub>-0.1</sub> × 10<sup>-2</sup></td></tr><tr><td>c-C<sub>3</sub>H<sub>2</sub></td><td>3<sub>21</sub>–2<sub>12</sub></td><td>18.2</td><td>-4.23</td><td>21</td><td>201.84</td><td>566.86</td><td>1597.89</td><td>5.7<sup>+4.0</sup><sub>-1.6</sub> × 10<sup>12</sup></td><td>...</sup>b</td></tr><tr><td>H<sup>13</sup>CN</td><td>3–2</td><td>24.9</td><td>-3.11</td><td>21</td><td>28.17</td><td>55.31</td><td>109.62</td><td>5.4<sup>+0.4</sup><sub>-1.3</sub> × 10<sup>12</sup></td><td>6.9<sup>+3.0</sup><sub>-3.0</sub> × 10<sup>-1</sup></td></tr><tr><td>HC<sup>15</sup>N</td><td>3–2</td><td>24.8</td><td>-3.12</td><td>7</td><td>9.42</td><td>18.50</td><td>36.66</td><td>1.1<sup>+0.2</sup><sub>-0.1</sub> × 10<sup>12</sup></td><td>1.4<sup>+1.1</sup><sub>-0.5</sub> × 10<sup>-1</sup></td></tr><tr><td>DCN</td><td>3–2</td><td>20.9</td><td>-3.34</td><td>21</td><td>33.39</td><td>65.75</td><td>130.51</td><td>3.7<sup>+0.5</sup><sub>-0.3</sub> × 10<sup>12</sup></td><td>3.6<sup>+3.7</sup><sub>-1.4</sub> × 10<sup>-1</sup></td></tr><tr><td>H<sup>13</sup>CO<sup>+</sup></td><td>3–2</td><td>25.0</td><td>-2.88</td><td>7</td><td>9.35</td><td>18.35</td><td>36.37</td><td>2.2<sup>+0.3</sup><sub>-0.1</sub> × 10<sup>12</sup></td><td>4.9<sup>+4.7</sup><sub>-1.9</sub> × 10<sup>-1</sup></td></tr><tr><td>CS</td><td>5–4</td><td>35.3</td><td>-3.53</td><td>11</td><td>16.28</td><td>32.24</td><td>64.15</td><td>1.0<sup>+2.6</sup><sub>-0.6</sub> × 10<sup>13</sup></td><td>5.0<sup>+4.8</sup><sub>-1.8</sub> × 10<sup>-1</sup></td></tr><tr><td>SO</td><td>7<sub>6</sub>–6<sub>5</sub></td><td>47.5</td><td>-2.11</td><td>15</td><td>38.88</td><td>90.35</td><td>197.51</td><td>&lt;1.5 × 10<sup>10</sup></td><td>...</sup>c</td></tr><tr><td>SO</td><td>6<sub>6</sub>–5<sub>5</sub></td><td>56.4</td><td>-2.21</td><td>13</td><td>38.88</td><td>90.35</td><td>197.51</td><td>&lt;7.1 × 10<sup>10</sup></td><td>...</sup>c</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Peak column densities and optical depths are derived for an excitation temperature of 30 K, with the exception of H<sub>2</sub>CO where T<sub>ex</sub> has been computed self-consistently. The error bars include the statistical uncertainties from the posterior distributions, a 10% flux calibration uncertainty, and for all lines but H<sub>2</sub>CO the systematic uncertainty in the excitation temperature.

<sup>b</sup> No radial profile was derived for c-C<sub>3</sub>H<sub>2</sub>. The column density represents the disk-averaged value, computed over a disk radius of 3" using the flux from Table 3.

<sup>c</sup> The 3σ upper limit on the line flux has been used to derive an upper limit on the disk-averaged column density, using the same approach as for c-C<sub>3</sub>H<sub>2</sub>.

Alarcón et al. 2020; Rab et al. 2020). To assess how different is the molecular emission morphology across different lines in PDS 70, and whether these are related to the underlying dust density distribution as traced by millimeter continuum thermal emission, the radial peak of the integrated intensity profile is analyzed. From the radial profiles of Section 3.4, the radius of the peak intensity is derived without any smoothing or fitting. An error of one-quarter of the beam major axis is assigned to the peak location. The measured values are reported in Table 4. For the two <sup>12</sup>CO lines and for HCO<sup>+</sup> only the off-center peak at $r > 0.2''$ is analyzed.

The summary plot comparing the peak locations to the 855 µm continuum emission presented in Isella et al. (2019) and the J-band polarized intensity by Keppler et al. (2018) is shown in Figure 11. The lines peak at different locations in the disk, with the large majority peaking within the radial maximum of the submillimeter continuum emission. The only two exceptions are CS and c-C<sub>3</sub>H<sub>2</sub>, peaking at $r \gtrsim 0.9''$ outside the continuum ring, similarly to what is observed in other systems (e.g., Teague et al. 2016; Cleeves et al. 2021). This general behavior indicates that the underlying gas density structure significantly affects most of the molecular lines'
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Top panel: normalized 855 µm continuum radial intensity profile (data by Isella et al. 2019) and normalized J-band polarized intensity profile (data by Keppler et al. 2018), with the two gray markers highlighting the orbital radius of PDS 70b and c (Wang et al. 2021). Bottom panel: the ALMA continuum intensity profile is shown in color code in physical units, where the synthesized beam is 67 × 50 mas. White dots indicate the radius of the integrated intensity profile peaks for all the detected species, listed on the y-axis. The C<sub>2</sub>H represents the *J* = 7/2-5/2 line including both detected hyperfine components; H<sub>2</sub>CO is shown using the brighter *3*<sub>03</sub>-*2*<sub>02</sub> line. The 1σ error bars represent one-quarter of the major axis of the beam of the respective images used to extract the radial profiles, which are listed in Table 4.</figcaption></figure>

emission. At the same time, the variety of peak locations within
the submillimeter cavity indicates how the flux of molecular
lines is determined by a combination of factors (density, kinetic
temperature, excitation temperature, optical depth, molecular
abundance, etc.). While reproducing all these emission lines
would require significant thermochemical modeling, some
trends can be identified heuristically. The most optically thick
lines, such as ¹²CO, peak at the innermost radii, indicating the
high temperatures at the edge of the cavity wall. Rarer CO
isotopologs, namely ¹³CO and C¹⁸O, progressively peak at
larger radii, with the lines better tracing bulk gas densities as
they become optically thinner. Interestingly, the C¹⁸O peak is
colocated with the peak in the dust continuum, likely indicating
a density maximum leading to efficient dust trapping. The HCN
isotopolog lines, which are found to have low optical depths
(see Section 4.2), peak at smaller radii than the optically thicker
C¹⁸O. This likely shows an excitation and abundance effect,
with HCN being mostly abundant at the gas cavity wall, where
the outer edge of the gap carved by the planet is directly
illuminated by UV photons from the central star triggering the
production of HCN via vibrationally excited H₂ (e.g.,
Cazzoletti et al. 2018; Visser et al. 2018). This is also indicated
by their emission peaking at the same radial range of the peak
in the IR scattered light image, which traces the cavity edge in
small grains. Similarly, HCO⁺ and H¹³CO⁺ peak at a similar
radial range, likely showing a high abundance of their
precursor H₃⁺, where X-ray photons from the central star
ionizing H₂ impinge onto the cavity wall (e.g., Cleeves et al.
2014, 2017). This indicates that the group of molecules for

which the emission is significantly determined by the
illumination pattern correlates better with IR scattered light
images, highlighting the illumination pattern onto the disk, than
with millimeter continuum images, which better highlight the
bulk dust density structure.

In summary, the radial profiles of the different molecules
show that the gas and dust densities, shaped by the massive
planets carving a large gap visible in ¹²CO, affect the molecular
emission. The large cavity of PDS 70 is an extreme case of a
gapped structure in the surface density profile, and it can
be used to constrain thermochemical models focusing on the
effect of gaps and rings in affecting chemical abundances and
molecular emission in protoplanetary disks (e.g., Facchini et al.
2018; Alarcón et al. 2020; Rab et al. 2020). In PDS 70, we have
shown that different molecules peak at a variety of radial
scales, with a separate group of molecules being sensitive to
different physicochemical factors. These observations will be at
the foundation of upcoming targeted thermochemical models of
PDS 70 to determine the physicochemical interplay between
dust and gas structure and molecular emission in transition
disks.

**6. Conclusions**

In this paper, we have presented the chemical inventory of
the planet-hosting disk PDS 70. The conclusions can be
summarized with the following:

1. Achieving a typical sensitivity of ~1 mJy beam⁻¹ in a 1 km s⁻¹ channel with an angular resolution of ~0.′′4–0.′′5,
we detect 16 transitions from 12 molecular species: ¹²CO, ¹³CO, C¹⁸O, H₂CO, C₂H, c-C₃H₂, H¹³CN, HC¹⁵N, DCN, HCO⁺, H¹³CO⁺, and CS. No emission from CH₃OH or SO is detected.

2. Integrated intensity maps and radial profiles of the integrated intensities show that the molecular emission originates primarily from a bright ring outside the orbit of the two detected planets. Only the two ¹²CO and the HCO⁺ lines are associated with emission from gas within the cavity carved by planets b and c.

3. Molecular lines peak at different radial locations, highlighting the complex interplay of temperature, density, and stellar irradiation in determining molecular abundances. The molecular structure indicates a chemically active layer at the edge of the cavity wall of the transition disk.

4. Column densities and optical depths are derived for all molecules with optically thin transitions. The excitation temperature is computed self-consistently for H₂CO, where two transitions are detected. For all other molecules a range of excitation temperatures is considered to estimate the uncertainty on the column density.

5. Deuteration and nitrogen fractionation profiles are derived from the hydrogen-cyanide isotopologs. Both show a high fractionation level, with fractionation increasing at larger radii, possibly suggesting a cold deuteration and fractionation pathway outside the millimeter-emitting ring.

6. The high C₂H/¹³CO flux ratio and a high lower limit on the N(CS)/N(SO) column density ratio indicate a disk-average C/O ratio above unity. This suggests that the actively accreting protoplanets are likely to be enriched with high C/O ratio gas material.

7. The emission of some molecular transitions, such as C¹⁸O J = 2–1, correlates well with the (sub)millimeter continuum emission profile, thus indicating that they are good tracers for the bulk disk density profile. Other groups of molecules, such as HCN and HCO⁺ isotopologs, peak in the same radial range of the IR scattered light intensity profile, highlighting the effect of direct illumination on the chemistry at the cavity edge.

PDS 70 provides a unique opportunity to witness how interactions between newly formed planets and their parental disk can sculpt the gas and dust distributions and hence influencing the chemical complexity of the atmosphere-building material that is actively accreted onto the planets. This work presents the first chemical inventory of the PDS 70 system, revealing a rich chemistry where radial variations in the line emission suggest a dynamical chemical structure that reacts to the planet-driven perturbations in the underlying physical structure. Upcoming high-resolution data, probing spatial

scales of ~10 au, will facilitate a comprehensive study of the chemical complexity at key locations in the disk, such as the highly irradiated cavity edge and along fast radial flows delivering the atmosphere-building material to the giant planets. Together with present and future projects trying to determine the chemical properties and C/O ratio in the atmospheres of PDS 70b and c (e.g., Wang et al. 2021), such work will provide the most direct link between the formation environment of an exoplanet and the elemental ratios found in its atmosphere, a key science case of current and upcoming facilities.

We thank the referee for their comments on the manuscript. We are thankful to Ewine F. van Dishoeck for early suggestions that helped to define this paper. We appreciate useful conversations with Ted Bergin, Dmitry Semenov, Romane Le Gal, and Teresa Paneque. We thank Mark Gurwell for performing the recalibration of the archival SMA data. This paper makes use of the following ALMA data:

ADS/JAO.ALMA#2015.1.00888.S,
ADS/JAO.ALMA#2017.A.00006.S,

ADS/JAO.ALMA#2019.1.01619.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), NSC 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. S.F. acknowledges an ESO Fellowship. This project has received funding from the European Unions Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 823823 (Dustbusters RISE project). J.B. acknowledges support by NASA through the NASA Hubble Fellowship grant #HST-HF2-51427.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555.

Facility: ALMA.

Software: GoFish (Teague 2019), VISIBLE (Loomis et al. 2018a), GALARIO (Tazzari et al. 2018), CASA (McMullin et al. 2007), Matplotlib (Hunter 2007), numpy (van der Walt et al. 2011), emcee (Foreman-Mackey et al. 2013).

## Appendix A
### Channel Maps of Weak Lines

The channel maps of the c-C₃H₂ 3₂₁–2₁₂ and H₂CO 3₂₁–2₂₀ transitions are reported in Figure 12. The images were produced with natural weighting and an additional uv-tapering of 0"/7. The properties of the images are listed in Section 3.2.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Channel maps of the c-C₃H₂ 3₂₁⁻₂₁₂ (top) and H₂CO 3₂₁⁻₂₂₀ (bottom) transitions, imaged with natural weighting and an additional uv-tapering of 0.′′7. The channel velocity is reported in the top-right corner of the panels in units of km s⁻¹. Dashed contours indicate the [3, 4, 5, 6, 7, 8]σ levels. The two white dots indicate the location of planets b and c.</figcaption></figure>

## Appendix B
### Radial Profiles with GoFish

In Section 3.4, radial profiles of the detected emission lines are computed by azimuthally averaging the integrated intensity maps shown in Figure 5. A different method to produce the intensity profiles has also been tested. The methodology is the same as described in Section 3.2 by the GoFish package, by deprojecting, shifting, and azimuthally averaging the spectra along radial annuli, as shown in Figure 4. The ¹²CO J = 3–2 and HCO⁺ lines are deprojected assuming the vertical component of the emission, with the same parameters derived by Keppler et al. (2019), and used to create the maps in the cleaning process (see parameters in Section 2). The 2D spectra are then collapsed along the velocity axis in order to obtain the

radial profile of the integrated intensity. The integral in velocity is performed on the following intervals: −1 to 12 km s⁻¹ for the ¹²CO lines, 0–11 km s⁻¹ for the HCO⁺ line, 0–12 km s⁻¹ for the C₂H J = 7/2, 5/2 line, and 1–10 km s⁻¹ for all the other lines. The profiles obtained with this method are shown in Figure 13, to be compared against the ones from the azimuthal average of the integrated intensity maps shown in Figure 6. For radii larger than the beam major axis, the agreement between the two methods is remarkable, whereas the inner regions of the stacking method lead to an underestimate of the integrated intensity due to beam and velocity smearing reducing the average intensity. Because most of the analysis in the paper focuses on the inner two to three beams, radial profiles from integrated intensity maps were preferred.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Azimuthally averaged radial profiles of the 14 detected lines and of the 885 µm continuum by Isella et al. (2019) in the bottom right in comparison. The profiles have been obtained with the GoFish package (see Section 3.4). The imaging parameters of the respective channel maps are listed in Table 4. The ribbon shows the 1σ rms, excluding the uncertainty from absolute flux calibration. The gray line at the top right of all panels shows the beam major axis. The brightness temperature conversion was done under the Rayleigh–Jeans approximation.</figcaption></figure>
ORCID iDs

Stefano Facchini https://orcid.org/0000-0003-4689-2684

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

Jaehan Bae https://orcid.org/0000-0001-7258-770X

Myriam Benisty https://orcid.org/0000-0002-7695-7605

Andrea Isella https://orcid.org/0000-0001-8061-2207

References

Alarcón, F., Teague, R., Zhang, K., Bergin, E. A., & Barraza-Alfaro, M. 2020, ApJ, 905, 68

Andrews, S. M. 2020, ARA&A, 58, 483

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

Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21

Bae, J., Zhu, Z., Baruteau, C., et al. 2019, ApJL, 884, L41

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

Bergner, J. B., Guzmán, V. G., Öberg, K. I., Loomis, R. A., & Pegues, J. 2018, ApJ, 857, 69

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

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

Cazzoletti, P., van Dishoeck, E. F., Visser, R., Facchini, S., & Bruderer, S. 2018, A&A, 609, A93

Christiaens, V., Casassus, S., Absil, O., et al. 2019, MNRAS, 486, 5819

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

Cleeves, L. I., Bergin, E. A., & Adams, F. C. 2014, ApJ, 794, 123

Cleeves, L. I., Bergin, E. A., Öberg, K. I., et al. 2017, ApJL, 843, L3

Cleeves, L. I., Loomis, R. A., Teague, R., et al. 2021, ApJ, 911, 29

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

Colzi, L., Fontani, F., Caselli, P., et al. 2018, A&A, 609, A129

Colzi, L., Sipilä, O., Roueff, E., Caselli, P., & Fontani, F. 2020, A&A, 640, A51

Cridland, A. J., Bosman, A. D., & van Dishoeck, E. F. 2020, A&A, 635, A68

Cridland, A. J., Eistrup, C., & van Dishoeck, E. F. 2019, A&A, 627, A127

Disk Dynamics Collaboration, Armitage, P. J., Bae, J., et al. 2020, arXiv:2009.04345

Dong, R., Hashimoto, J., Rafikov, R., et al. 2012, ApJ, 760, 111

Dutrey, A., Wakelam, V., Boehler, Y., et al. 2011, A&A, 535, A104

Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2016, A&A, 595, A83

Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2018, A&A, 613, A14

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

Facchini, S., Pinilla, P., van Dishoeck, E. F., & de Juan Ovelar, M. 2018, A&A, 612, A104

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

Fedele, D., & Favre, C. 2020, A&A, 638, A110

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

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

Guilloteau, S., Reboussin, L., Dutrey, A., et al. 2016, A&A, 592, A124

Guzmán, V. V., Öberg, K. I., Huang, J., Loomis, R., & Qi, C. 2017, ApJ, 836, 30

Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, NatAs, 3, 749

Hashimoto, J., Dong, R., Kudo, T., et al. 2012, ApJL, 758, L19

Hashimoto, J., Tsukagoshi, T., Brown, J. M., et al. 2015, ApJ, 799, 43

Heays, A. N., Visser, R., Gredel, R., et al. 2014, A&A, 562, A61

Hily-Blant, P., Bonal, L., Faure, A., & Quirico, E. 2013a, Icar, 223, 582

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

Hily-Blant, P., Magalhaes, V., Kastner, J., et al. 2017, A&A, 603, L6

Hily-Blant, P., Pineau des Fortès, G., Faure, A., Le Gal, R., & Padovani, M. 2013b, A&A, 557, A65

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

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

Isella, A., Benisty, M., Teague, R., et al. 2019, ApJL, 879, L25

Isella, A., & Turner, N. J. 2018, ApJ, 860, 27

Kama, M., Bruderer, S., van Dishoeck, E. F., et al. 2016, A&A, 592, A83

Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44

Keppler, M., Teague, R., Bae, J., et al. 2019, A&A, 625, A118

Krijt, S., Bosman, A. D., Zhang, K., et al. 2020, ApJ, 899, 134

Lammer, H., & Blanc, M. 2018, SSRv, 214, 60

Linsky, J. L., Draine, B. T., Moos, H. W., et al. 2006, ApJ, 647, 1106

Long, F., Bosman, A. D., Cazzoletti, P., et al. 2021, A&A, 647, A118

Long, Z. C., Akiyama, E., Sitko, M., et al. 2018, ApJ, 858, 112

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

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

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

Madhusudhan, N. 2019, ARA&A, 57, 617

Magalhães, V. S., Hily-Blant, P., Faure, A., Hernand ez-Vera, M., & Lique, F. 2018, A&A, 615, A52

Marty, B., Chaussidon, M., Wiens, R. C., Jurewicz, A. J. G., & Burnett, D. S. 2011, Sci, 332, 1533

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

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

Mesa, D., Keppler, M., Cantalloube, F., et al. 2019, A&A, 632, A25

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

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

Miotello, A., Facchini, S., van Dishoeck, E. F., et al. 2019, A&A, 631, A69

Mollière, P., Stolker, T., Lacour, S., et al. 2020, A&A, 640, A131

Mordasini, C. 2018, in Planetary Population Synthesis, ed. H. J. Deeg & J. A. Belmonte (Berlin: Springer), 143

Müller, A., Keppler, M., Henning, T., et al. 2018, A&A, 617, L2

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

Mumma, M. J., & Charnley, S. B. 2011, ARA&A, 49, 471

Öberg, K. I., & Bergin, E. A. 2016, ApJL, 831, L19

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

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

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

Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJL, 743, L16

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

Pinilla, P., Benisty, M., Birnstiel, T., et al. 2014, A&A, 564, A51

Pudritz, R. E., Cridland, A. J., & Alessi, M. 2018, in Connecting Planetary Composition with Formation, ed. H. J. Deeg & J. A. Belmonte (Berlin: Springer), 144

Rab, C., Kamp, I., Dominik, C., et al. 2020, A&A, 642, A165

Roueff, E., Loison, J. C., & Hickson, K. M. 2015, A&A, 576, A99

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

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

Stolker, T., Marleau, G. D., Cugno, G., et al. 2020, A&A, 644, A13

Tazzari, M., Beaujean, F., & Testi, L. 2018, MNRAS, 476, 4527

Tazzari, M., Testi, L., Natta, A., et al. 2020, arXiv:2010.02248

Teague, R. 2019, JOSS, 4, 1220

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

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

Terwisscha van Scheltinga, J., Hogerheijde, M. R., Cleeves, L. I., et al. 2021, ApJ, 906, 111

Terzieva, R., & Herbst, E. 2000, MNRAS, 317, 563

Testi, L., Birnstiel, T., Ricci, L., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 339

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

van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2016, A&A, 585, A58

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

van Dishoeck, E. F., & Bergin, E. A. 2020, arXiv:2012.01472

Visser, R., Bruderer, S., Cazzoletti, P., et al. 2018, A&A, 615, A75

Wagner, K., Follete, K. B., Close, L. M., et al. 2018, ApJL, 863, L8

Wang, J. J., Vigan, A., Lacour, S., et al. 2021, AJ, 161, 148

Winn, J. N., & Fabrycky, D. C. 2015, ARA&A, 53, 409

Wirström, E. S., & Charnley, S. B. 2018, MNRAS, 474, 3720

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

Zhang, K., Bergin, E. A., Schwarz, K., Krijt, S., & Ciesla, F. 2019, ApJ, 883, 98