ALMA high-resolution observations unveil planet formation shaping molecular emission in the PDS 70 disk

L. Rampinelli; S. Facchini; M. Leemker; J. Bae; M. Benisty; R. Teague; C. J. Law; K. I. Öberg; B. Portilla-Revelo; A. J. Cridland

doi:10.1051/0004-6361/202449698

Home

All issues

Volume 689 (September 2024)

A&A, 689 (2024) A65

Full HTML

Open Access

Issue		A&A Volume 689, September 2024


Article Number		A65
Number of page(s)		24
Section		Planets and planetary systems
DOI		https://doi.org/10.1051/0004-6361/202449698
Published online		03 September 2024

A&A, 689, A65 (2024)

ALMA high-resolution observations unveil planet formation shaping molecular emission in the PDS 70 disk^★

L. Rampinelli¹, S. Facchini¹, M. Leemker¹, J. Bae², M. Benisty³^,4, R. Teague⁵, C. J. Law⁶^,**, K. I. Öberg⁷, B. Portilla-Revelo⁸^,9 and A. J. Cridland¹⁰

¹ Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria 16, 20133 Milano, Italy
e-mail: luna.rampinelli@unimi.it
² Department of Astronomy, University of Florida, Gainesville, FL 32611, USA
³ Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange, 06304 Nice, France
⁴ Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
⁵ Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
⁶ Department of Astronomy, University of Virginia, Charlottesville, VA 22904, USA
⁷ Center for Astrophysics | Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA
⁸ Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands
⁹ Department of Astronomy & Astrophysics, The Pennsylvania State University, 525 Davey Laboratory, University Park, PA 16802, USA
¹⁰ Universitäts-Sternwarte München, Ludwig-Maximilians-Universität, Scheinerstr. 1, 81679 München, Germany

Received: 22 February 2024
Accepted: 2 July 2024

Abstract

With two directly detected protoplanets, the PDS 70 system is a unique source in which to study the complex interplay between forming planets and their natal environment. The large dust cavity carved by the two giant planets can affect the disk chemistry, and therefore the molecular emission morphology. On the other hand, chemical properties of the gas component of the disk are expected to leave an imprint on the planetary atmospheres. In this work, we reconstruct the emission morphology of a rich inventory of molecular tracers in the PDS 70 disk, and we look for possible chemical signatures of the two actively accreting protoplanets, PDS 70b and c. We leverage Atacama Large Millimeter/submillimeter Array (ALMA) band 6 high-angular-resolution and deep-sensitivity line emission observations, together with image and uv-plane techniques, to boost the detection of faint lines. We robustly detect ring-shaped emission from ¹²CO, ¹³CO, C¹⁸O, H¹³CN, HC¹⁵N, DCN, H₂CO, CS, C₂H, and H¹³CO⁺ lines in unprecedented detail. Most of the molecular tracers show a peak of the emission inside the millimeter dust peak. We interpret this as the direct impact of the effective irradiation of the cavity wall, as a result of the planet formation process. Moreover, we have found evidence of an O-poor gas reservoir in the outer disk, which is supported by the observations of bright C-rich molecules, the non-detection of SO, and a lower limit on the CS/SO ratio of ~1. Eventually, we provide the first detection of the c-C₃H₂ transitions at 218.73 GHz, and the marginal detection of an azimuthal asymmetry in the higher-energy H₂CO (3_2,1−2_2,0) line, which could be due to accretion heating near PDS 70b.

Key words: astrochemistry / protoplanetary disks / stars: individual: PDS 70

^★

The reduced datacubes are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/689/A65

^**

NASA Hubble Fellowship Program Sagan Fellow.

© The Authors 2024

Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

In the past decade, advances in exoplanet detection and characterization have allowed for demographic insights. They have revealed a broad diversity in exoplanet properties, and thus imply that the outcomes of planet formation are similarly diverse (Suzuki et al. 2016; Thompson et al. 2018; Fulton et al. 2021; Vigan et al. 2021; Currie et al. 2023; Lissauer et al. 2023, and references therein). Mature planetary systems store important markings of their formation history: for example, the composition of their atmosphere may be crucial in determining the chemical inheritance from their natal environment; that is, the protoplanetary disk (Öberg & Bergin 2021; Nomura et al. 2023; Öberg et al. 2023, and references therein). In recent years, detailed information about the chemical properties and composition of gas giant exoplanets has been inferred from transmission and thermal emission spectra thanks to powerful ground-based facilities, or space missions (Greene et al. 2016; Hinkley et al. 2023; Currie et al. 2023; Guillot et al. 2023). Analyzing the obtained spectra with radiative transfer thermo-chemical models has revealed unprecedented details about planetary atmospheres, constraining fundamental properties such as the C/O elemental ratio (Mollière et al. 2015) or carbon fractionation (Zhang et al. 2021).

On the one hand, planet formation is strongly influenced by the chemical and physical processing of the disk material during the disk evolution, which is expected to affect the location, timescale, mass budget, or chemical composition of forming planets. On the other hand, a forming planet is also expected to leave an imprint on the disk, such as substructures in the dust and gas distributions, chemical signatures, and kinematics perturbations, which have been observed and used to predict the planet’s properties (Bae et al. 2023; Benisty et al. 2023; Pinte et al. 2023, and references therein).

In this context, high spatial and spectral resolution observations of molecular lines with ALMA can be leveraged to detect both kinematic and chemical footprints of planet formation. The kinematic variations induced onto the Keplerian rotation pattern are observed as localized spatial and spectral features (Perez et al. 2015b; Pinte et al. 2018, 2019, 2023; Casassus & Pérez 2019; Teague et al. 2019; Izquierdo et al. 2022, 2023; Stadler et al. 2023). Chemical signatures induced by planetary accretion shocks can be targeted through specific molecular lines (Cleeves et al. 2015; Law et al. 2023), or chemical footprints of planet formation can be observed through substructures into the gas distribution (Öberg et al. 2021b; Facchini et al. 2021; Nomura et al. 2021; Law et al. 2021a; Bae et al. 2022; Jiang et al. 2023).

Planetary signatures in molecular line emission observations complement constraints on planet formation coming from studies on the dust distribution, which have been largely explored in the past few years (Andrews et al. 2018; Long et al. 2018a; Andrews 2020; Cieza et al. 2021; Sierra et al. 2021; Bae et al. 2023), resulting in the recent detection of a circumplanetary disk in the PDS 70 system through ALMA continuum observations (Benisty et al. 2021). Moreover, molecular line emission observations can provide a way to access open questions related to the astrochemical characterization of planet-forming disks, such as elemental composition (Bosman et al. 2021; Cataldi et al. 2021), evolution of organic chemistry (Guzmán et al. 2021; Ilee et al. 2021), and in particular sulfur chemistry (Le Gal et al. 2021), or snowline location.

In this context, the PDS 70 system is the best candidate to further explore the impact of planet formation on the chemical structure of the disk, and on the other hand how the chemical properties of the disk can result in the final atmospheric composition of the accreting planets (Cridland et al. 2023), as it is the first source for which a multiwavelength direct detection of two forming planets has been presented, in near-infrared (Keppler et al. 2019; Christiaens et al. 2019; Mesa et al. 2019), Hα line (Wagner et al. 2018; Haffert et al. 2019), ultraviolet (Zhou et al. 2021), and millimeter and submillimeter observations (Isella et al. 2019; Benisty et al. 2021). PDS 70 is a ~5 Myr old T Tauri star with a 0.87 M_⊙ mass (Paxton et al. 2011; Hashimoto et al. 2015; Choi et al. 2016; Keppler et al. 2018, 2019; Long et al. 2018b) at a distance of ~113 pc (Gaia Collaboration 2018) in the Upper Centaurus-Lupus subgroup. Facchini et al. (2021) presented the first chemical inventory of the PDS 70 disk, showing a complexity of morphologies and emitting properties, and evidence suggesting a gas-phase C/O > 1. Law et al. (2024) extracted the high-resolution vertical gas structure of the PDS 70 disk, showing strong evidence of a prominent cavity wall at the ¹²CO-emitting height, which could directly impact the molecular emission morphology. Moreover, Perotti et al. (2023) found water in the JWST spectrum of the inner disk of PDS 70, which could affect potential terrestrial planets being formed.

In this work, we analyze ALMA band 6 line emission observations of the planet-hosting disk around PDS 70. We present new high-angular-resolution observations of multiple molecular lines, providing an overview of the molecular complexity in the PDS 70 disk in unprecedented detail. In particular, we suggest that the effect of the strongly irradiated cavity wall is playing an important role in shaping the molecular emission. In Sect. 2, we describe the observations together with calibration and imaging procedures. In Sec. 3, we describe the image-plane and uv-plane techniques we applied to boost faint lines. In Sect. 4, we outline the analyses we performed in order to build a chemical overview of the observed molecules. In Sec. 5, we interpret the observed emission morphologies of the analyzed molecular tracers in light of the ongoing process of planet formation. We summarize the results and discussion in Sect. 6.

2 Observations and data reduction

In this paper, we present ALMA band 6 line emission observations from the ALMA program #2019.1.01619.S (PI: S. Facchini): we complement the results for the short baselines, previously shown by Facchini et al. (2021), with the long-baseline data. We performed the following data reduction through the CASA software (McMullin et al. 2007) v5.8.

The observations included two spectral setups. The calibration of the lower-frequency setup is reported by Law et al. (2024). We calibrated the high-frequency observations with the same procedure. These observations consist of two execution blocks (EBs) for the short spacings, and five EBs in an extended configuration. In summary, we first realigned the observations of each EB by shifting the phase center to the center of an ellipse fitted to the outer continuum ring of PDS 70. To perform the shift, we used the CASA tasks fixvis and fixplanets. After flagging the channels in the proximity of the lines detected by Facchini et al. (2021) in a ±15 km s⁻¹ range from the systemic velocity, we averaged the data into 250 MHz channels. We then self-calibrated the continuum data, the short baselines first, and then the concatenated short and long baselines (as in Andrews et al. 2018). For the short-baseline data, we used the same solution intervals as in Facchini et al. (2021), whereas for data combining short and long baselines we computed and applied a phase-only solution down to 120 s. We ran a final round of amplitude and phase self-calibration on EB-long intervals. The resulting continuum image reconstructed with the tclean task, using the multiscale deconvolver, with Briggs weighting and a robust parameter of 0.5 shows a root mean square (RMS) of 9.27 μJy. The peak signal-to-noise ratio (S/N) is 234, over a synthesized beam of $0_{\cdot}^{''} 14 × 0_{\cdot}^{''} 10$ $\[0^{\prime\prime}_\cdot 14\text { \times }0^{\prime\prime}_\cdot 10\]$ (with a PA of −85.3°). The total flux density of the disk at the weighted frequency of 252.8 GHz is 68.42 ± 0.10mJy, without accounting for 10% absolute flux calibration uncertainties, which dominate but only result in a systematic shift. The flux uncertainty was evaluated from standard deviations of fluxes from random ellipses with a minor axis of $1_{\cdot}^{''} 5 × \cos (i)$ $\[1^{\prime\prime}_\cdot 5\text { \times }\cos ~(i)\]$ (i is the disk inclination) and PA as the disk, taken outside the emitting area (the same procedure was applied for line fluxes, see Sec. 4.3 for more details). The gain solutions were then applied to the full spectral data. Finally, we subtracted the continuum emission using the uvcontsub task with a first-order polynomial.

We imaged the different lines using the tclean task, based on the CLEAN algorithm (Högbom 1974), using natural weighting. We applied the multiscale deconvolver, with scales (0, 5, 10, 20, 30)pix, where pix indicates the pixel dimension ( $0_{\cdot}^{''} 2 × 0_{\cdot}^{''} 2$ $\[0^{\prime\prime}_\cdot 2\text { \times }0^{\prime\prime}_\cdot 2\]$ ). We tested the effect of different imaging procedures on the final image, such as different flux thresholds to stop the cleaning (2σ and 3.5σ) and two different masking techniques: we built a first mask by cleaning down to 7σ, and a second one by leveraging the Keplerian rotation, using the keplerian_mask.py tool by Teague (2020), and we obtained consistent images. We fixed the geometrical parameters, M_* = 0.875 M_⊙, i = 51.7°, PA= 160.4°, and υ_sys = 5.5 km s⁻¹ (in LSRK frame), as they were obtained by Keppler et al. (2019), and a distance of 113 pc (Gaia Collaboration 2018). We subsequently convolved the mask with a Gaussian kernel with a full width at half maximum (FWHM) of $0_{\cdot}^{''} 4$ $\[0^{\prime\prime}_\cdot 4\]$ to reproduce the broadening effect of the finite resolution in the data. The systematic difference in the resulting channel maps due to the different masks are within the statistical uncertainties. We performed all the following analyses on cubes imaged using the 7σ mask.

In order to maximize the detection of the faint lines H₂CO (3_2,1−2_2,0), c-C₃H₂, SO, and CH₃OH, we imaged their data cubes with an additional uv-tapering of $0_{\cdot}^{''} 2$ $\[0^{\prime\prime}_\cdot 2\]$ to increase sensitivity to extended emission. We performed all the following analyses on the continuum subtracted cubes.

The spectral setup includes several molecular transitions, including CO isotopologs (¹²CO, C¹⁸O, ¹³CO), formaldehyde (H₂CO), cyanides (H¹³CN, HC¹⁵N, DCN), hydrocarbons (C₂H, c-C₃H₂), S-bearing molecules (CS, SO), ions (H¹³CO⁺), and methanol (CH₃OH). In particular, we observed two molecular transitions of c-C₃H₂, one at 244.2221 GHz (3_2,4−2_1,2) and another at 218.7327 GHz: the latter corresponds to the two ortho (7_1,6−7_0,7) and para (7_2,6−7_1,7) lines, which are blended. In the paper, we refer to the 3_2,1−2_1,2 transition as “c-C₃H₂ (244 GHz)” and to the blended one as “c-C₃H₂ (218 GHz).” An overview of the analyzed transitions, corresponding beam size, RMS, and line flux is presented in Table 1 (see also Sec. 4.3).

We accounted for the so-called JvM effect (Jorsater & van Moorsel 1995), following Czekala et al. (2021), by rescaling the residual map by the ϵ factor, which is the ratio between areas of the CLEAN and DIRTY beams. ϵ factors are listed in Table 1 for each line. Non-JvM corrected images show higher residuals, as they are not rescaled, resulting in an overestimated line flux, by a factor that we estimated to be ~10% for bright lines, reaching ~40% for faint lines. In this work, the analysis was performed on the JvM corrected cubes, except for undetected lines, not to overestimate the sensitivity of the observations (Casassus & Cárcamo 2022). For the same reason, we used non-JvM corrected images in the rotational diagram analysis of the H₂CO asymmetry (see Sec. 3.1), as it is performed on spatial scales equal to the beam size.

Table 1

Listed imaged lines, corresponding rest frequencies, imaging parameters (channel width, beam, RMS), and line fluxes.

3 Boosting faint lines detection

In this work, we applied boosting techniques both in the image and in the visibility plane, to the weak lines included in the dataset: H₂CO (3_2,1−2_2,0), c-C₃H₂ (244 GHz), c-C₃H₂ (218 GHz), SO (6₇−5₆), and SO (6₆−5₅). Techniques in the image plane leverage the disk rotation pattern, to stack the spatially integrated line spectrum, by shifting back the spectrum along each spatial pixel to the systemic velocity, and boost the S/N (Yen et al. 2016, and Fig. A.1). This method has been implemented in the GoFish package (Teague 2019), and allows one to reconstruct a shifted spectrum by correcting for the Keplerian rotation (see also Appendix A).

On the other hand, these techniques are sensitive to imaging artifacts, spatial covariance on the beam scale and/or spectral information loss due to spectral masks (Loomis et al. 2018). Directly applying boosting techniques in the visibility plane could prevent limitations linked to the image reconstruction. In this context, multiple lines of the same molecule can be stacked in the uv-plane (Walsh et al. 2016), or a matched filter analysis can be applied by spanning the sampled visibilities of the faint line with a model filter (Loomis et al. 2018).

Fig. 1

Integrated intensity maps of the detected lines, and of the 855 μm continuum (Isella et al. 2019), with the white contours showing the bright ring in the submillimeter continuum emission and the white dots marking the position of the two forming planets (Wang et al. 2021). Brightness temperatures were obtained under the assumption of Rayleigh-Jeans approximation. The white ellipse at the bottom left of each panel is the beam, while the white line at the bottom right of the last panel indicates the 100 au scale.

3.1 Detection of H₂CO (3_2,1–2_2,0)

Two molecular lines of H₂CO are included in the observed spectral window: the bright ring-shaped 3_0,3−2_0,2 line, which is detected at a high S/N (see Fig. 1) and the faint and asymmetric 3_2j1 −2₂,o transition. For the case of the faint H₂CO 3_2,1−2_2,0 line, the Keplerian mask and natural weighting during the CLEAN-ing procedure were enough to retrieve the Keplerian pattern in the channel maps. The emission distribution is represented in the integrated intensity map in Fig. 1, where the signal is spatially resolved and asymmetric.

To confirm the robustness of the detection, we used the GoFish package (Teague 2019) to extract the Keplerian shifted spatially integrated spectrum within a circle of 2″ of radius, which is represented in the top right plot in Fig. A.1. We also obtained the Keplerian shifted spectrum as a function of radius, in order to highlight the spatial location of the emitted signal, as it is represented in the bottom right panel in Fig. A.1. In this case, image-plane techniques were enough to confirm the detection, with an S/N of 9.5 (evaluated from the integrated flux, see Sec. 4.3), which confirms the marginal detection of 4.4σ presented by Facchini et al. (2021).

3.2 Detection of c-C₃H₂

We applied boosting techniques in the image plane to the two lines of the c-C₃H₂ molecule: we imaged the two lines with natural weighting and an additional uv-tapering of $0_{\cdot}^{''} 2$ $\[0^{\prime\prime}_\cdot 2\]$ , along with a Keplerian mask, as we described in Sec. 2. The emission is at least marginally detected in the channel maps of both lines, as it can be seen from the top and middle panels in Fig. B.1, where we highlight the 3σ contours in white. From the integrated intensity maps (see first column in Fig. 2), we also extracted the radial profile of the integrated intensity. The results are presented in the middle column of Fig. 2.

To assess the significance of the signal, we extracted the Keplerian shifted spectra as a function of radius (Teague 2019), which are shown in the right column of Fig. 2, with the dashed white lines defining the 3K level from the ¹²CO line. The two shifted spectra are centered on the systemic velocity, supporting the detection. As is also evident from the integrated intensity maps (left column Fig. 2), the two lines show a different emission morphology, with the 218 GHz ones being more compact. This could be due to the fact that the two transitions have different upper energies, and are sensitive to different environmental conditions, resulting in different emitting regions. This can be critical when extracting the rotational diagram analysis, as we discuss further in Sec. 5.5.

In summary, c-C₃H₂ 218 GHz lines are marginally detected with an S/N of 3.3 (undetected in the short-baseline data, Facchini et al. 2021), while the c-C₃H₂ 244 GHz transition is detected with an S/N of 8 (marginally detected by Facchini et al. 2021 with an S/N of 4.3). Signal-to-noise ratios were derived from disk-integrated fluxes, extracted after applying boosting techniques (see Sec. 4.3).

To confirm the detection, we applied the matched filter analysis directly in the uv-plane, using the visible package (Loomis et al. 2018). We built a model filter from a Keplerian mask using the Keplerian_mask.py tool developed by Teague (2020), assuming the geometrical and line parameters that we applied in the imaging, but varying the radial extent, in order to find the mask that maximizes the filter response, when spanning the kernel over the sampled visibilities. The best results for the filter responses to the c-C₃H₂ lines use a mask with an inner radius of 0″ and an outer radius of 1″ (218 GHz) and 2″ (244 GHz). The corresponding filter responses are shown in the left and middle panels of Fig. 3, as a function of velocity. We applied the analysis after running the ms tables through the cvel2 task, which corrects for Doppler shifts throughout the time of the observation, and transforms the reference frame from TOPO to LSRK, knowing the rest frequencies of the targeted lines. We therefore expect the peak of the filter response to be at the systemic velocity, as it is marked by the vertical red lines for the targeted faint c-C₃H₂ lines, in Fig. 3. The matched filter response peaks at the systemic velocity (5.5 km s⁻¹), and thus confirms the detection of the two lines, with a peak filter response of 5.7σ (for the 244 GHz line) and of 4.2σ (for the 218 GHz lines).

Fig. 2

Detection of the c-C₃H₂ transitions at 218 GHz (top row) and at 244 GHz (bottom row). Panels in the left column show the integrated intensity maps obtained from the imaged data cubes with natural weighting and an additional uv-tapering. The white contours show the band 7 submillimeter continuum bright ring, and the two white dots indicate the location of the two gas giants (Wang et al. 2021), while the white ellipses on the bottom left are the beams. Panels in the central column represent the radial profiles of the azimuthally averaged integrated intensity, obtained as described in Sec. 4.2, for c-C₃H₂ transitions (blue line), and for the 855 μm continuum emission (dashed black line, Isella et al. 2019). The hatched gray region indicates the beam major axis. Panels in the right column show the Keplerian shifted line spectra, in brightness temperature, as a function of radius, for each $0_{\cdot}^{''} 1$ $\[0^{\prime\prime}_\cdot 1\]$ wide annulus in which the disk has been divided. The dashed white lines are the 3K level of the ¹²CO.

3.3 uv-stacking of SO

The imaging techniques we applied for the other weak molecular transitions – that is, using natural weighting, and an additional uv-tapering, or Keplerian shifting – did not lead to a detection when applied to the two SO 6₇−5₆ and SO 6₆−5₅ lines. We therefore applied a stacking technique in the uv-plane, by combining the two individually undetected lines, as it has been previously presented by Walsh et al. (2016) for the case of the methanol detection in the TW Hya disk. As the two lines have different rest frequencies, we re-gridded the data using cvel2, such that the central frequency of each line refers to the systemic velocity, and they can subsequently be correctly stacked. We then combined the two measurement sets using the CASA task concat: Fig. C.1 shows a test of the uv-stacking performed on two molecular lines of H₂CO that are also individually detected, where the integrated intensity profile of the stacked cube (pink) lies between the higher and the lower profiles, representing the integrated intensity profiles of the bright (orange) and the faint (blue) lines, respectively.

We then imaged the stacked SO cube, using the same parameters applied for the re-gridding in cvel2 (channel width, spectral extent, starting channel). As we expected from the stacking procedure, the stacked data cube shows an RMS $\sim \sqrt{2}$ $\[\sim \sqrt{2}\]$ times lower than the original ones. We did not detect any signal in the SO channel maps of the uv-stacked cube above the 3σ level. The right panel in Fig. A.2 shows the Keplerian shifted spectra as a function of radius, for the original undetected lines (left and middle panels), and for the uv-stacked cube (right panel), without any clear boosted emission linked to the expected Keplerian pattern of the disk. The bottom panel in Fig. B.1 shows the channel maps obtained from the cleaned uv-stacked data cube, where the white contours show the 3σ level. The only marginal detection originates from the outer disk at 4.10 km s⁻¹ (see Sec. 5.6 for a detailed discussion).

We applied the matched filter analysis to the stacked visibilities, testing various model filters, built from Keplerian masks, or from the emission of other molecules. The strongest filter response is obtained using H₂CO 3_0,3−2_0,2 as a filter, as it shows extended emission, and is thus more sensitive to the 3σ feature observed in the outer disk. However, the result from the matched filter analysis is not significant, since it is below the 3σ threshold, as is shown in Fig. 3.

Fig. 3

Filter response as a function of velocity, as we obtained from the matched filter analysis performed on the c-C₃H₂ (218 GHz, left panel), and the C₃H₂ (244 GHz, middle panel) transitions, with a Keplerian mask as filter: the vertical dashed red line indicates the systemic velocity and the position of the peaks corresponding to the two c-C₃H₂ lines. Right panel: impulse response from the matched filter analysis performed on the SO stacked visibilities, using the clean image of the H₂CO (3_0,3−2_0,2) line as a filter. The horizontal red line shows the 3σ level.

4 Overview of the observed molecules

In this section, we present an overview of the analyzes we applied to reconstruct the emission morphology and line fluxes from the high resolution observations (see Facchini et al. 2021 for reference on the methodologies).

4.1 Disk-integrated spectra

After reconstructing the emission distribution of the various molecular transitions, we extracted line spectra, by applying the stacking method in the image plane, presented in Sec. 3. We used the GoFish package (Teague 2019), following Teague et al. (2016), under the assumption of an axisymmetric Keplerian disk and from the geometrical parameters and the star mass presented in Sec. 2. We assumed midplane emission for all molecules except for the most elevated ¹²CO and ¹³CO emission. For the latters, we included a tapered power law prescription of the emitting layer height z(r): $z (r) = {\begin{cases} z_{0} \times {(\frac{r - r_{cavity}}{1^{''}})}^{ϕ} \times \exp (- {[\frac{r}{r_{taper}}]}^{ψ}) & r > r_{cavity} \\ 0 & r \leq r_{cavity}, \end{cases}$ $\[z(r)= \begin{cases}z_0 \times\left(\frac{r-r_{\text {cavity }}}{1^{\prime \prime}}\right)^\phi \times \exp \left(-\left[\frac{r}{r_{\text {taper }}}\right]^\psi\right) & r>r_{\text {cavity }} \\ 0 & r \leq r_{\text {cavity }},\end{cases}\]$ (1)

with $z_{0} = 0_{\cdot}^{''} 41, r_{cavity} = 0_{\cdot}^{''} 37, ϕ = 0.51, r_{taper} = 1_{\cdot}^{''} 27, ψ = 5.74$ $\[z_0=0^{\prime \prime}_\cdot41, r_{\text {cavity }}=0^{\prime \prime}_\cdot37, \phi=0.51, r_{\text {taper }}=1^{\prime \prime}_\cdot27, \psi=5.74\]$ for $^{12} C O, and z_{0} = 0_{\cdot}^{''} 41, r_{cavity} = 0_{\cdot}^{''} 28, ϕ = 1.29, r_{taper} = 0_{\cdot}^{''} 74, ψ = 1.61$ $\[{ }^{12} \mathrm{CO} \text {, and } z_0=0^{\prime \prime}_\cdot41, ~r_{\text {cavity }}=0^{\prime \prime}_\cdot28, ~\phi=1.29, ~r_{\text {taper }}=~0^{\prime \prime}_\cdot74 \text {, }~\psi=1.61\]$ for ¹³CO (Law et al. 2024). We then extracted the integrated spectrum by spatially integrating the shifted spectra over a de-projected circle of radius 3″ for the lines with the most extended emission (¹²CO and ¹³CO), and 2″ for the other ones. Radial extents are chosen to be large enough not to exclude extended emission, but not to include excessive noise in the calculation. An example of the extracted disk-integrated spectra is shown in the top row of Fig. A.1.

4.2 Integrated intensity maps and radial profiles

We extracted the integrated intensity maps of the detected lines by spectrally integrating the JvM-corrected data cubes with the bettermoments tool (Teague & Foreman-Mackey 2019). We collapsed the cubes over the velocity ranges listed in Table 1, chosen after visual inspection of the data to include all channel maps showing signal, without applying any additional mask. Figure 1 shows the integrated intensity maps of the detected lines. The position of the two planets is marked with white dots (Wang et al. 2021), while the white contours show the bright ring in the ALMA band 7 submillimeter continuum emission (Isella et al. 2019). The azimuthal dark feature in the CS map is a channelization effect due to low spectral resolution. All molecules have a structured emission morphology, with a ringed-shape emission, not always colocated with the bright ring in the submillimeter continuum.

We extracted the radial profiles from the integrated intensity maps shown in Fig. 1, using the GoFish package (Teague 2019), by dividing each map into annuli with a radial extent of 1/4 of the beam major axis, and taking the azimuthal average of the integrated intensity for each annulus (see Facchini et al. 2021 for more details). The result is shown for each line in Fig. 4: the integrated intensity is expressed in integrated brightness temperature, assuming the Rayleigh-Jeans approximation. Uncertainties are obtained by taking the standard deviation along each annulus, divided by the square root of the number of independent beams in the annulus. For radii smaller than the beam major axis we took only the standard deviation along the annuli. The uncertainties may be underestimated, especially in the inner regions, inside the beam FWHM, as it does not account for the 2D covariance introduced by the beam.

The azimuthal average is performed after de-projecting the emission, by using the geometrical parameters presented in Sec. 2, assuming the surfaces from Eq. (1) for ¹²CO and ¹³CO, and midplane origin for the other molecules. All molecules in Fig. 4 have a peaked integrated intensity profile, with the ¹²CO, H¹³CN, DCN, and H₂CO (3_2,1−2_2,0) lines showing the innermost peak, inside the millimeter dust peak. Each molecule presents an inner cavity, with H¹³CN, H¹³CO⁺, and C₂H showing the steepest radial integrated intensity gradient, while H₂CO, C₂H, and ¹²CO also show an outer shoulder of emission. Most of the molecules show a peak of the emission inside the submillimeter continuum peak, which we discuss further in Sec. 5.1, in light of the efficiently irradiated cavity wall.

Fig. 4

Azimuthally averaged radial profiles of the integrated intensity for the analyzed molecular transitions. The ribbons show the standard deviation across each annulus, divided by the square root of the number of independent beams. The FWHM of the synthesized beam is highlighted by the hatched region, while the beam size of the continuum image is indicated on the top-right corner of the last panel. The dashed black line shows the azimuthally averaged radial profile of the 855 μm continuum emission (Isella et al. 2019).

4.3 Line fluxes

The line flux of each transition is obtained from the integrated intensity maps in Fig. 1, by spatially integrating the maps over de-projected circles of the same extent used for the integrated spectra, in Sec. 4.1. This procedure has been applied to all transitions except for the faintest two lines of c-C₃H₂. In the latter case, we extracted the line flux from the shifted spectra as a function of radius (see right panel in Fig. 2). We then applied a mask based on the ¹²CO 3 K brightness temperature level, which is marked by the dashed white line in the right panel in Fig. 2. We chose to define such a mask to select a region large enough to include all the line flux, but with a different radial extent for each velocity channel, to avoid including excessive noise. The 3 K level from ¹²CO is a good indicator of the extent of the emission, expected from a Keplerian disk (see also bottom-left panel in Fig. A.1). This procedure is similar to applying a Keplerian mask when extracting the integrated intensity map, but it prevents one from producing artifacts when creating the image. We tested this procedure on the two weakest lines, in order to robustly extract their line flux: their shifted spectra as a function of radius and the mask are shown in the right panel in Fig. 2.

Table 1 lists line fluxes and associated uncertainties for each line. The uncertainty is evaluated as the standard deviation of the flux measured in 26 de-projected circles, outside the emitting area of the integrated intensity maps, with the same radius used to extract the flux. We took the maximum number of circles that was possible to place on the image without overlap. This procedure has been performed on images which are not primary beam corrected, to ensure uniform noise. For the case of c-C₃H₂, we extracted the uncertainty associated with the line flux from the shifted spectra as a function of radius, by taking the standard deviation of 26 flux estimates, which we evaluated inside the 3K mask, randomly applied outside the spectral range used to evaluate the line flux.

For the undetected SO (6₇−5₆), SO (6₆−5₅), and CH₃OH (5_1,4,₀−4_1,3,0) lines, we provide an upper limit on the flux, evaluated from the integrated intensity maps, as 3 times the standard deviation of the flux of 26 de-projected circles, taken outside the possible emitting region. For the case of the C₂H (3_7/2−2_5/2) transition, as the two hyperfine components are spectrally resolved, we extracted the integrated flux of the individual components by fitting a double-Gaussian to the shifted-stacked spectrum (blue line in the top-middle panel in Fig. A.1), under the assumption of two optically thin components. The uncertainty on the flux estimate has been extracted as for the other lines from the integrated intensity map.

Table 2

Line parameters for the various species for which we estimated the column density.

4.4 Disk-averaged column densities

From the flux estimates listed in Table 1, we extracted the disk averaged molecular column density for all molecules (listed in Table 2), except for ¹²CO and ¹³CO, which are expected to be optically thick. We assumed local thermodynamic equilibrium (LTE), and optically thin emission, following Facchini et al. (2021). LTE is a reasonable assumption in disks, as density becomes high enough at the midplane: for the case of PDS 70 the midplane is reaching a density of ~10⁸ cm⁻³ (Portilla-Revelo et al. 2023). For the molecular transitions and the excitation temperatures (see below) we are considering, critical densities are typically between 5 × 10⁵ cm⁻³ and 5 × 10⁶ cm⁻³ (Shirley 2015): the LTE approximation is therefore generally valid for an emitting layer z/r ≲ 0.25 (Portilla-Revelo et al. 2023). The optically thin assumption is expected to be a good approximation for most of the considered molecules: Law et al. (2024) showed (see Fig. 12) that C¹⁸O is becoming marginally optically thick at its peak at the location of the bright millimeter dust ring. If we consider the peak of C¹⁸O in Fig. 4 as benchmark, we see that all the molecules are below the threshold, except for ¹²CO and ¹³CO, for which we did not derive the column density. We took C¹⁸O as reference, as it is expected to emit from the lowest layer of all molecular transitions discussed in this work, with the other molecules having a kinetic temperature T_kin ≳ T_kin(C¹⁸O), and therefore a lower optical depth than the C¹⁸O (2–1) transition, if their brightness temperature is lower than the C¹⁸O one. Both C₂H and H₂CO are becoming marginally optically thick, but only at their peak: we therefore included the two molecules to extract the column density, which needs to be taken as lower limit where the molecule is becoming marginally optically thick. We used the line parameters E_u (upper state energy), g_u (upper state degeneracy), Q(T) (partition function), and A_ul (Einstein coefficient) listed in Table 2, which we retrieved from the CDMS database (Pickett et al. 1998; Müller et al. 2001). We evaluated the partition functions by interpolating with a cubic spline, over tabulated values.

We extracted the average molecular column density for each line, inside the de-projected circle where we evaluated the line fluxes: the results are listed in Table 2, for an excitation temperature of 30 K. The errors include the propagation of the statistical uncertainty over the line flux and the systematic uncertainty from the unknown excitation temperature, which we varied between 20 and 50 K. This range is consistent with the excitation temperature range extracted for H₂CO (see following discussion). We therefore evaluated the column density N_i at different excitation temperatures within the assumed range, and the related statistical uncertainty propagated from the flux uncertainty dN_i. To evaluate the final uncertainty over the disk averaged column density, we took the largest value of N_i + dN_i and the lowest value of N_i − dN_i. We did not include the 10% absolute flux calibration uncertainty, as it would result only in a systematic shift. For the undetected lines, we provide an upper limit on the average molecular column density from the upper limit on the line flux (see Table 1).

If more than one line is detected for a certain molecule, both the excitation temperature T_ex and the total molecular column density N_t can be extracted by leveraging a rotational diagram analysis (Goldsmith & Langer 1999). We applied this methodology to the two lines of the H₂CO and the c-C₃H₂ molecules respectively. For the latter, we took into account the ortho- and para- spin isomers of c-C₃H₂, presenting a line at the same rest frequency (218.73 GHz) and with the same line parameters, except for the upper state degeneracy (g_u = 45 for ortho and g_u = 15 for para, assuming an ortho-to-para ratio of 3:1, see Table 2). We extracted the disk averaged rotational diagram inside a de-projected circular area, with a radius of 2″ for the more extended H₂CO and of $1_{\cdot}^{''} 5$ $\[1^{\prime\prime}_\cdot 5\]$ for the more compact c-C₃H₂, obtaining T_ex = (25 ± 2) K and N_t = (2.0 ± 0.2) × 10¹³ cm⁻² for H₂CO, T_ex = (15 ± 4) K and N_t = (6.8 ± 1.8) × 10¹² cm⁻² for c-C₃H₂, as it is shown in the two top panels of Fig. 5. The radius of $1_{\cdot}^{''} 5$ $\[1^{\prime\prime}_\cdot 5\]$ for the case of c-C₃H₂ has been chosen to balance the more extended line at 244 GHz and the more compact one at 218 GHz, as it is visible from Fig. 2. For c-C₃H₂ we also included the 10% absolute flux calibration uncertainty (sum in quadrature with flux statistical uncertainty), as the two lines are in different spectral setups, while we did not include it for H₂CO, as the two lines are in the same spectral setup, with absolute flux calibration uncertainty possibly resulting in the same systematic shift. We explored the posterior distribution through the MCMC method implemented in the emcee package (Foreman-Mackey et al. 2013), with 128 walkers, 5000 burn-in steps, and 500 steps. Error bars refer to the 16th and 84th percentiles of the posterior distribution. We considered a uniform prior with T_ex ∈ (0, 50) K, and N_t ∈ (10¹¹, 10¹⁴) cm⁻². Our result for the H₂CO molecule is consistent with what Facchini et al. (2021) presented, while they did not perform the rotational diagram analysis on c-C₃H₂, as the faint lines at 218 GHz were undetected.

Fig. 5

Rotational diagram analysis of H₂CO and c-C₃H₂. First row: disk-averaged rotational diagram analysis applied to the two lines of the H₂CO and c-C₃H₂ molecules, respectively. The average excitation temperatures and column densities are indicated on the top of the two panels. The dashed lines correspond to the best fit results for N_t and T_ex, while the solid lines are 100 random results from the posterior distribution of the emcee fit. Bottom row: radial profiles of the total column density (left panel) and of the excitation temperature (right panel) from H₂CO rotational diagram analysis, applied to $0_{\cdot}^{''} 1$ $\[0^{\prime\prime}_\cdot 1\]$ wide annuli. The solid lines show the median of the posterior distribution (which is representative of the best fit as the posterior distribution is well described by a Gaussian) from the emcee fit, for N_t and T_ex, while the ribbons represent the 16th and 84th percentiles. The horizontal dark blue line shows the results from the disk-average rotational diagram (from top left panel). Gray regions mark the radial range of the expected location of the CO snowline (Law et al. 2024).

Fig. 6

Radial profiles of the column density for the observed molecules. The solid blue lines refer to the value extracted for an excitation temperature of T_ex = 30 K, while the ribbons show the uncertainty, which takes into account the propagated uncertainty on the flux density and the systematic uncertainty on the excitation temperature, which has been varied between 20 and 50 K (Facchini et al. 2021). The 10% absolute flux calibration uncertainty is not included. The horizontal dark blue line shows the disk-averaged value, listed in Table 2. The dashed black line shows the azimuthally averaged radial profile of the 855μm continuum emission (Isella et al. 2019). Hatched regions indicate the beam major axis of the cubes.

4.5 Radial profiles of column density

From the radial profiles of the integrated intensity (Fig. 4), we extracted the molecular column density as a function of radius for all the molecules listed in Table 2 (see Fig. 6), by computing the column density along $0_{\cdot}^{''} 1$ $\[0^{\prime\prime}_\cdot 1\]$ wide annuli, using the line parameters listed in Table 2. We varied the excitation temperatures from 20 to 50 K, as for the disk-averaged column density. The solid blue lines in Fig. 6 show the radial profiles of the total column density obtained for the different molecules, at the assumed T_ex = 30 K.

For the H₂CO molecule we also applied a radial rotational diagram, to retrieve both T_ex(r) and N_t(r), where r is the cylindrical radius. We divided the disk into $0_{\cdot}^{''} 1$ $\[0^{\prime\prime}_\cdot 1\]$ wide annuli and we sampled the posterior distribution for each annulus using an MCMC method with the emcee package (Foreman-Mackey et al. 2013), with 128 walkers, 5000 burn-in steps and 500 steps. We considered a uniform prior with T_ex ∈ (0, 100) K, and N_t ∈ (10¹¹, 10¹⁶) cm⁻². Radially resolved H₂CO excitation temperature and column density are shown in the bottom panels in Fig. 5. Ribbons indicate the 16th and 84th percentiles of the posterior distribution, and the propagated uncertainty on the integrated intensity has been considered when generating the likelihood. The column density peaks at smaller radii with respect to the submillimeter continuum bright ring, which is consistent with the radial profile of the integrated intensity for the two H₂CO lines, in Fig. 4. The radial profile of the excitation temperature is almost flat, with T_ex between 20 and 30 K already from r = $0_{\cdot}^{''} 4$ $\[0^{\prime\prime}_\cdot 4\]$ inside the dust cavity (discussed in Sec. 5.4). The higher excitation temperature at smaller radii is likely due to the dominant contribution from the inner emission observed in the faint and higher energetic H₂CO (3_2,1−2_2,0) line, showing a hot spot at the cavity wall near PDS 70b (see Sec. 5.3 for deeper discussion).

Table 3

Radial position of the peak of the integrated intensity profiles, and radii enclosing 68 and 90% of the total flux.

5 Discussion

5.1 Molecular emission peaking at the cavity wall

From the integrated intensity maps in Fig. 1 and their azimuthally averaged radial profiles in Fig. 4, we found that the emission morphologies in high resolution data are in good agreement with Facchini et al. (2021).

To highlight the main emitting features, we list the radial location of the emission peak for each molecule in Table 3, and we compare them in Fig. 7 (bottom panel), together with the radial profile of the ALMA band 7 continuum intensity (top panel, blue line, Isella et al. 2019) and of the SPHERE J-band polarized intensity (top panel, pink line, Keppler et al. 2018). We extracted the radial position of the peaks by fitting a multiple-Gaussian profile (triple for ¹²CO and H₂CO, double for the other molecules) to each radial profile in Fig. 4, using scipy.optimize.curve_fit (see Fig. D.1). Uncertainties on the peaks position have been obtained from the scipy fit, taking into account also the uncertainties on the integrated intensity. For the faint c-C₃H₂ transitions we could not perform a Gaussian fitting, and we therefore provide the peak position along with 1/4 of the beam major axis as an associated uncertainty. In Table 3 we also list the radius corresponding to the 68% and 90% of the radius enclosing the whole line flux, as computed from a cumulative flux radial profile (Facchini et al. 2019; Sanchis et al. 2021).

As apparent from the bottom panel in Fig. 7, molecular emission does not necessarily correlate with the dust distribution. This is also consistent with what was shown by the MAPS Large Program (Öberg et al. 2021b) through high resolution line emission observations of five protoplanetary disks: all disks show a variety of chemical substructures which can critically differ from the dust ones (Law et al. 2021a,b).

For the specific case of the PDS 70 disk, high resolution line emission observations show that most of the molecules present a peak of the emission at smaller radii with respect to the submillimeter continuum bright ring (bottom panel of Fig. 7). Different tracers are affected by specific chemical and physical conditions, and thus peak at different radii, not always tracing the underlying gas density structure.

The emission of CO isotopologs is mainly affected by their different optical depth, with the optically thicker ¹²CO tracing the gas temperature, and thus peaks at an inner radius (~ $0_{\cdot}^{''} 43$ $\[0^{\prime\prime}_\cdot 43\]$ ) with respect to the submillimeter continuum bright ring (~ $0_{\cdot}^{''} 67$ $\[0^{\prime\prime}_\cdot 67\]$ ), while the optically thinner C¹⁸O better traces the bulk of the gas mass and peaks at ~ $0_{\cdot}^{''} 63$ $\[0^{\prime\prime}_\cdot 63\]$ (Portilla-Revelo et al. 2023). This trend in the peak of CO isotopologs is also consistent with planet-disk interaction models (Pinilla et al. 2012; Facchini et al. 2018) and with what it has been found for other disks (e.g., Perez et al. 2015a; van der Marel et al. 2016; Leemker et al. 2022).

Except for the CO isotopologs, another effect which can play a role in shaping the emission of the other molecular tracers is the irradiation field. As the two forming planets carved a large dust cavity, radiation from the central star can penetrate freely and efficiently irradiate the cavity wall. PDS 70 is the first transition disk with clear evidence of a prominent cavity wall, with a steep rise in emitting height and a vertical dip in the ¹²CO (3–2) line (Law et al. 2024). The emission of molecules whose production is enhanced in a strong irradiation field is expected to peak at the warm, strongly irradiated cavity wall rather than following the gas density distribution (Portilla-Revelo et al. 2023). Among these tracers, HCN isotopologs, C₂H, and H¹³CO⁺ peak inside the bright submillimeter ring, showing a correlation with the peak in the scattered light profile (see Fig. 7) which traces the rise of the density of small grains from the cavity wall (Keppler et al. 2018; Portilla-Revelo et al. 2023). In particular, C₂H is expected to be enhanced in a strong UV-field, and also HCN, as it is produced by UV-pumped H₂. Similarly, ions such as H¹³CO⁺ are formed in a strong X-ray field environment outside the water snowline. While the submillimeter continuum profile traces the distribution of large grains at the midplane, the NIR scattered light profile traces the disk illumination pattern. The correlation between the peak of the emission for the molecules listed above and the NIR peak is therefore consistent with their production being favored by photons from the star which can freely reach the cavity wall.

HCN isotopologs also show different emission morphologies, with H¹³CN being the most compact one, while HC¹⁵N and DCN present a more extended emission, possibly linked to active in-situ fractionation pathways. In particular, isotope selective photodissociation and the effective self-shielding of N2 could strongly affect H¹³CN and HC¹⁵N inside the UV-irradiated cavity (Heays et al. 2014; Visser et al. 2018; Hily-Blant et al. 2019). In particular, this process could explain the observed shift in the H¹³CN and HC¹⁵N peaks: while H¹³CN peaks at an inner radius, HC¹⁵N peaks at an outer radius (similar to C¹⁸O), which could reflect the fact that N₂ easily self-shields at a higher gas density, with respect to N¹⁵N. On the other hand, low temperature isotope exchange reactions may play a crucial role in producing the outer emission shoulders in the HC¹⁵N and DCN emission (Millar et al. 1989; Roueff et al. 2015; Huang et al. 2017; Öberg et al. 2021a; Cataldi et al. 2021; Muñoz-Romero et al. 2023). Fractionation processes presented above result in different HCN isotopologs peaking at different radii: ~ $0_{\cdot}^{''} 48$ $\[0^{\prime\prime}_\cdot 48\]$ for H¹³CN, ~ $0_{\cdot}^{''} 60$ $\[0^{\prime\prime}_\cdot 60\]$ for HC¹⁵N, and ~ $0_{\cdot}^{''} 47$ $\[0^{\prime\prime}_\cdot 47\]$ for DCN. A deeper characterization of the fractionation levels of the HCN molecule in the PDS 70 disk, their radial profiles and implications for the planet-forming environment will be investigated in a future work, as well as further modeling to constrain the role of the strongly irradiated cavity wall in setting the observed emission morphology.

In this picture, the bright ring in the submillimeter continuum emission appears to induce a chemical separation, between an inner warm and strongly irradiated cavity wall and an outer low-temperature and low-density environment. While most of the molecules show an inner bright ring, a few tracers also present outer peaks or shoulders of emission (see Fig. 4), which are due to different chemical reactions specifically unlocked by environmental conditions set beyond the continuum bright ring, as we previously outlined for the case of the DCN molecule. Similarly, also H₂CO shows an outer peak at ~1.27″, which we discuss further in Sec. 5.4.

In this context, we note that all C-rich molecules present bright emission beyond the submillimeter dust ring: CS peaks at ~1″, c-C₃H₂ (244 GHz) peaks at ~ $0_{\cdot}^{''} 83$ $\[0^{\prime\prime}_\cdot 83\]$ , and C₂H shows a bright outer shoulder beyond ~ $0_{\cdot}^{''} 87$ $\[0^{\prime\prime}_\cdot 87\]$ . This could be indicative of a carbon-rich chemistry boosted in the outer disk, and a corresponding radially increasing gas-phase C/O ratio, which we discuss further in Sect. 5.7. c-C₃H₂ (218 GHz) is the only exception, with centrally peaked compact emission: this is consistent with the higher upper state energy of the transition (E_u = 61 K), and thus better traces high temperature at smaller radii.

In Fig. D.3 we compare the position of the peaks of the emission extracted from the integrated intensity (in white) and from the peak intensity (in red) radial profiles, obtained by multiple-Gaussian fitting (see Fig. D.2). Radial profiles of the peak intensity were obtained by azimuthally averaging the peak intensity maps, which show the peak of the spectrum for each spatial pixel (see Fig. E.1). Peak intensity profiles better trace optical depth, as they do not depend on the line width like integrated intensity profiles, which include contribution at the same location both from the line center and line wings, where optical depth can change (Weaver et al. 2018; Rosotti et al. 2021). As line width increases at smaller radii, we expect the peaks from integrated intensity profiles to be at an inner radius with respect to peaks from peak intensity profiles: this is visible in Fig. D.3. However, the peak positions are consistent within the error bars for most of the molecules, and thus robustly confirm the observed trend of the peaks of line emission, which are inside the millimeter dust peak for most of the molecular tracers. This result is also consistent with the fact that the peaks in Fig. E.1 are spatially and spectrally resolved: beam dilution, which is expected to affect the peak position extracted from peak intensity profiles (Leemker et al. 2022), is thus negligible in this case.

We summarize the peak positions extracted from peak intensity profiles in Table D.1: the two only tracers with inconsistent peaks are C₂H and ¹³CO, with the peak position from peak intensity profiles being further out with respect to integrated intensity profiles, as expected. The most critical case is C₂H (R_peak = $0_{\cdot}^{''} 51$ $\[0^{\prime\prime}_\cdot 51\]$ from integrated intensity, R_peak = $0_{\cdot}^{''} 83$ $\[0^{\prime\prime}_\cdot 83\]$ from peak intensity): this can be due to the fact that the molecule presents two similar peaks at two different radii (see Fig. E.1), which weigh differently in the peak intensity and in the integrated intensity profiles, as explained above. Moreover, from the radial profile of the peak intensity in Fig. D.2, we notice that the H₂CO molecule is the only one showing a gap beyond the submillimeter bright ring, followed by an outer ring of emission (see also Sec. 5.6).

Fig. 7

Radial location of the peak of the emission of the various molecular tracers, compared to radial profiles of SPHERE scattered light and ALMA band 7 continuum intensity. Top panel: radial profiles of the ALMA band 7 continuum intensity (blue line, Isella et al. 2019) and SPHERE J-band polarized intensity (pink line, Keppler et al. 2018). The gray shadow shows the radius of the coronagraph in the J-band observations (Keppler et al. 2018). Bottom panel: radial position of the peak of the emission and its related uncertainty, for each detected line. The background color-bar refers to the radial profile of the ALMA band 7 continuum intensity. The two vertical dashed white lines indicate the planets position (Wang et al. 2021).

5.2 Comparison with other transition disks

Substructures in the line emission of molecular tracers are the product of a variety of physical and chemical pathways, and are therefore not always colocated with dust substructures (Öberg et al. 2021b, and references therein). We showed that the PDS 70 system fully agrees with this scenario (see Sec. 5.1), but in this case the efficient irradiation of the cavity wall could play an important role in enhancing molecular emission of specific tracers inside the dust cavity, which could be a specificity of transition disks.

In this context, we compare our result with two other transition disks around HD 169142 (Garg et al. 2022; Booth et al. 2023b) and HD 142527 (Temmink et al. 2023), for which a chemical overview has been provided through ALMA line emission observations. As different disks around different stars present different temperature structure, irradiation field, size and mass, in order to provide a consistent comparison and look for a possible correlation between molecular emission and the cavity wall, we extracted the ratio between the radius corresponding to the peak of the emission of various molecular tracers and the one of the innermost bright ring in the continuum emission (we call it “peak ratio” for simplicity).

CO isotopologs and CS show the most consistent trend among the three disks. CS always peaks outside the innermost continuum bright ring, with a peak ratio of ~1.5, 1.5 and 2.2 for the PDS 70, HD 142527, and HD 169142 disk respectively. In this context, Temmink et al. (2023) extracted a CS/SO ratio >1, as we obtained for the PDS 70 disk in this work (see Sec. 5.7), while Booth et al. (2023b) presented a radially increasing CS/SO ratio. These first attempts highlight the importance of better constraining the C/O ratio and its radial variations across transition disks directly from observations. The volatile C/O ratio has been shown to be strongly linked to planetary accretion (Jiang et al. 2023): on one hand, the C/O ratio of the gas accreted by the planet will be inherited by the planetary atmosphere, but on the other hand, a forming planet that is carving a gap in the disk is also expected to locally change the C/O ratio in the gas phase due to, for example, the sublimation of C-rich molecules.

The consistent trend shown by peaks of CO isotopologs among the three disks robustly confirm the dominant effect of their different optical depths, as we previously discussed in Sec. 5.1. In particular, the peak ratios increase with decreasing optical depths: for ¹²CO it is ~0.7, 0, and 0.6, for ¹³CO is ~0.8, 0.7 and 0.8, and for C¹⁸O is ~0.9, 1.0, and 1.0, for PDS 70, HD 142527, and HD 169142, respectively. In this context, a critical difference between the three sources is the presence of a gas gap/cavity: while both the PDS 70 and the HD 169142 disk present a gas gap/cavity clearly visible in the ¹²CO emission, ¹²CO extends towards the disk center for the HD 142527 disk.

On the other hand, except for CO isotopologs and CS, there are not any other consistent trends among the peak ratios for other molecular tracers. For example, for the disk around HD 169142 (Booth et al. 2023b), both H₂CO and c-C₃H₂ peak further out (peak ratios of ~5 and 4 respectively), with respect to the PDS 70 disk (peak ratios of ~0.8 and 1.2 respectively). Similarly, if we compare peak ratios of observed lines of HCN isotopologs, we obtain ~0.74 for H¹³CN in the PDS 70 disk and ~1.5 for HCN in the HD 142527 disk. Even if these are different isotopologs and could therefore be sensitive to different chemical pathways, we would have anyway expected a peak of the emission inside the inner cavity, if the irradiation of the cavity wall had played a major role. However, Cazzoletti et al. (2018) implemented a set of disk models through the thermo-chemical code DALI (Bruderer et al. 2012; Bruderer 2013), to show that rings in the CN emission can be traced back to its main formation pathway through UV-pumped H₂, which is mostly abundant in a ring-shaped region close to the surface layer of the disk. This process does not require the presence of a cavity, but results in a CN ring also for smooth disks. In particular, CN rings are at larger distances from the central star in disks around Herbig stars than T Tauri stars, due to their different UV fields. HCN is also produced by UV-pumped H₂, and we therefore expect it to be strongly affected by different UV fields, which could partially explain the observed difference between the radial position of the HCN isotopologs rings, as HD 169142 and HD 142527 are Herbig stars, while PDS 70 is a T Tauri. Moreover, PDS 70 shows the largest gas cavity in the ¹²CO emission (~43 au with respect to ~15 au for the HD 169142 disk). PDS 70 is the only disk for which a vertical dip at the cavity wall is clearly visible in ¹²CO line emission, revealing a prominent, warm, and strongly irradiated cavity wall (Law et al. 2024).

Nevertheless, the observed diversity in the molecular emission morphology throughout the three transition disks suggests that a further and deeper characterization of line emission observations and further model calculations in transition disks is needed in order to shed light on the effect of the inner cavity on the disk chemistry, and to assess if the observed chemical properties of the planet-hosting disk around PDS 70 make it a peculiar planet-forming environment.

Fig. 8

Azimuthal asymmetry of the H₂CO 3_2,1−2_2,0 line emission. Left panel: rotational diagram analysis applied to the two regions, highlighted in the integrated intensity maps of the two lines of the H₂CO molecule: the pink region corresponds to PA = 0°, while the blue one is at PA = 180° along the major axis. The dashed pink and blue lines correspond to the best fit N_t and T_ex (reported in the top right corner) obtained for the pink and blue region, respectively. Right panel: spatially integrated spectra of the H₂CO 3_2,1–2_2,0 line in the pink region and the blue region (lines in corresponding colors). The gray error bar in the top right of the panel is the uncertainty on the spatially integrated flux derived from the non-JvM corrected cube, which is the same for each velocity bin.

5.3 Azimuthally asymmetric emission of H₂CO

The two forming planets in the cavity of the PDS 70 disk have been shown to be actively accreting material from their environment, through ultraviolet (for PDS 70b) and H$alpha; line emission observations, which have been used to derive an estimate of the mass accretion rate over PDS 70b and c (Wagner et al. 2018; Aoyama & Ikoma 2019; Haffert et al. 2019; Zhou et al. 2021). Similarly, we could expect that the process of accretion onto the forming planets results in enhanced emission of particular molecules observed with ALMA (Cleeves et al. 2015). In this context, we looked for accretion signatures in the PDS 70 disk through high excitation lines, as they could be more sensitive to local heating mechanisms. We found an azimuthal asymmetry in the H₂CO (3_2,1−2_2,0) faint line emission, with a localized hot spot at the cavity wall at the same azimuth as PDS 70b, which could be linked to planetary accretion heating. This is particularly evident in the peak intensity map in Fig. E.1, but also in the integrated intensity map in Fig. 1. In order to test the significance of the azimuthal asymmetry, we extracted the spatially integrated spectrum within an area of the beam size localized on the emission feature, at an angular position of 0° along the major axis, and from the opposite angular position at 180°, as we show in the integrated intensity maps in the left panel of Fig. 8. The integrated spectra are shown in the right panel of Fig. 8, where the solid blue line refers to PA = 0°, while the pink one to PA = 180°. These do show an asymmetry between the blue and red-shifted parts of the disk. Nevertheless, the line showing the azimuthal asymmetry is faint, and the statistical significance of the spectral difference is low (2.5σ). We extracted the significance by taking the ratio of the difference between the two peaks and $\sqrt{2} \times$ $\[\sqrt{2} \times\]$ the spectrum uncertainty, which has been evaluated as the standard deviation of the fluxes from 50 random beam-sized regions outside the disk emitting area, in a random velocity channel. The uncertainty was extracted from the non-JvM corrected cube, since the JvM correction underestimates the uncertainty on the flux density of point sources by its correction factor, as we highlighted in Sect. 2.

To further explore the physical significance of the asymmetry, we applied the rotational diagram analysis in the two regions (pink and blue, as we previously described) to the two H₂CO (3₀_,3−2_0,2) and (3_2,1−2_2,0) lines, following the procedure we explained in Sect. 4.4. The result is presented in the left panel of Fig. 8, where the pink line refers to the results for the pink region, while the blue one refers to the blue region. As expected, we obtained a higher excitation temperature (~42 K) at the peak position, with respect to the opposite one (~22 K), but the two temperatures are consistent within the error bars, making the difference not statistically significant. We notice that we did not detect significant asymmetries in the H₂CO (3₀_,3−2_0,2) emission, or in any other accretion tracers. However, this result highlights the importance of further exploring planetary signatures in high energy transitions, as the only feature we identified is from the higher energy H₂CO line.

5.4 On the origin of the observed H₂CO emission

As we previously highlighted in Sects. 4.2 and 4.5, the emission of the bright H₂CO (3_0,3−2_0,2) line presents a peak at smaller radii with respect to the submillimeter continuum ring (see Fig. 4 and the bottom-left panel of Fig. 5), and an outer shoulder of emission extending beyond the outer edge of the millimeter continuum disk. The presence of an outer shoulder of emission for the H₂CO molecule has been observed in other sources (Loomis et al. 2015; Carney et al. 2017; Guzmán et al. 2021; Terwisscha van Scheltinga et al. 2021; Booth et al. 2023b; Temmink et al. 2023), but its interpretation is not straightforward, as it has been linked both to a possible gas-phase or solidstate origin. In particular, the dominant gas phase formation pathway is (Fockenberg & Preses 2002; Atkinson et al. 2006): ${C H}_{3} + O \to H_{2} C O + H,$ $\[\mathrm{CH}_3+\mathrm{O} \rightarrow \mathrm{H}_2 \mathrm{CO}+\mathrm{H},\]$ (2)

while the main grain surface formation route is through hydrogenation of CO ice (Hiraoka et al. 1994; 2002; Watanabe & Kouchi 2002; 2004; Hidaka et al. 2004; Fuchs et al. 2009): $C O \overset{H}{\to} H C O \overset{H}{\to} H_{2} C O .$ $\[\mathrm{CO} \xrightarrow{\mathrm{H}} \mathrm{HCO} \xrightarrow{\mathrm{H}} \mathrm{H}_2 \mathrm{CO}.\]$ (3)

If H₂CO originates from ice sublimation we would expect some correlation between H₂CO emission and dust grain locations, either through a correlation with the millimeter continuum bright ring (desorption from millimeter dust), or with the H₂CO emitting layer being consistent with or below the near infrared surface, where the μm dust is expected to reside. We therefore retrieved the H₂CO emitting layer, from the 2D temperature structure of the PDS 70 disk, which was extracted by Law et al. (2024), combining CO isotopologs and HCO⁺. The rotational diagram assumes that the excitation temperature T_ex(r) is equal to the kinetic temperature T(r, z) (LTE). We therefore extracted the emitting layer z(r) by inverting the T(r, z) relation (Law et al. 2024), as it has been previously performed by Ilee et al. (2021) for the case of complex organic molecules in the MAPS sample. We obtained that the H₂CO emission originates from z/r ~ 0.2 (see Fig. 9, solid black line), under the assumption of a geometrically thin emitting layer. We highlight that the LTE approximation is valid at the inferred layer, as the critical density of ~6 × 10⁵ cm⁻³ for H₂CO at the derived excitation temperatures is reached at z/r ≳ 0.3. The inferred layer for H₂CO is overlaid on the 2D temperature structure (Law et al. 2024) and compared to the emission surfaces of ¹²CO and ¹³CO (solid green and blue lines).

The trend of the H₂CO excitation temperature is almost flat, and between 20 and 30 K, which suggests a low-temperature, radical-radical gas-phase formation pathway. This scenario is also supported by the absence of an evident correlation between the H₂CO emission and the submillimeter continuum emission, and the detection of molecular emission from various molecules such as C₂H, c-C₃H₂, H¹³ CN whose gas-phase production is boosted in UV-exposed regions. Moreover, thermal-desorption from icy grains is an unlikely process to release H₂CO, as both the H₂CO excitation temperature and the ¹²CO kinetic temperature are lower than the typical H₂CO sublimation temperature of ~70 K (Noble et al. 2012; Fedoseev et al. 2015; van’t Hoff et al. 2020), for $r ≳ 0_{\cdot}^{''} 4$ $\[r \gtrsim 0^{\prime\prime}_\cdot 4\]$ .

Nevertheless, the grain-surface formation pathway cannot be completely ruled out, as H₂CO could be released also through nonthermal desorption, such as photo-desorption, in UV-exposed regions. This scenario is consistent with the low excitation temperature we inferred (~20–30 K), corresponding to CO freeze-out, and with the detection of bright UV-boosted molecules. Moreover, this formation route is supported by the fact that the radial profile of the H₂CO column density (see right panel, bottom row in Fig. 5) peaks at the position of the CO snowline (Öberg et al. 2017; Guzmán et al. 2021), which lays between ~ $0_{\cdot}^{''} 5$ $\[0^{\prime\prime}_\cdot 5\]$ and $0_{\cdot}^{''} 75$ $\[0^{\prime\prime}_\cdot 75\]$ (Law et al. 2024). In this context, desorption from μm grains cannot be ruled out, as the inferred H₂CO emitting layer is slightly above the NIR surface (dashed red line in Fig. 9, Keppler et al. 2018), but uncertainties do not allow one to distinguish the two scenarios. However, it is important to note that the efficiency of nonthermal desorption from icy grains is not well constrained (Öberg et al. 2017), and distinguishing between grain-surface and gas-phase pathways is not trivial, as one can continue to contribute also at the onset of the other (Walsh et al. 2014; Loomis et al. 2015; Öberg et al. 2017).

Interestingly, the H₂CO excitation temperature we obtained (between 20 and 30 K) is lower than the H₂CO freeze out temperature of ~70 K (Noble et al. 2012; Fedoseev et al. 2015; van’t Hoff et al. 2020), as it has been also shown for other sources (Terwisscha van Scheltinga et al. 2021; Temmink et al. 2023): this implies that H₂CO needs to be efficiently produced, either through gas-phase reactions or nonthermally released from icy grains, with faster timescales with respect to the freeze-out one.

Fig. 9

H₂CO emitting layer height as a function of radius (solid black line), and the related uncertainty (black ribbons), compared to the 2D temperature structure in the background (Law et al. 2024). The two solid green and blue lines refer to the ¹²CO and ¹³CO emission surfaces, with related uncertainties, respectively (Law et al. 2024). The dashed red line refers to the NIR scale height as a function of radius, from Keppler et al. (2018), while the red diamond is the height of the NIR emitting layer at the NIR ring ~ $0_{\cdot}^{''} 48$ $\[0^{\prime\prime}_\cdot 48\]$ (Law et al. 2024). Dashed black lines indicate fixed z/r layers.

5.5 Low emitting layer of c-C₃H₂

c-C₃H₂ has been previously detected in a few other disks (Bergin et al. 2016; Ilee et al. 2021; Booth et al. 2023b). In particular, Ilee et al. (2021) presented the detection of c-C₃H₂ in four disks in the MAPS sample (GM Aur, AS 209, HD 163296, and MWC 480). They extracted a disk average column density of c-C₃H₂ between 1.1 × 10¹² cm⁻² and 7.0 × 10¹² cm⁻² across the sample, which is consistent with the disk average column density of ~6.8 × 10¹² cm⁻² we extracted from the integrated rotational diagram analysis (see the second panel, first row, of Fig. 5). In particular, we applied a disk averaged rotational diagram, to fit both for the excitation temperature and for the molecular column density, as we previously outlined in Sec. 4.4.

The excitation temperature we obtained is (15 ± 4) K, which is low if we consider that c-C₃H₂ is expected to originate in the gas phase in the disk atmosphere (Loison et al. 2017), while T_ex ~ 15 K in the PDS 70 disk implies that the emitting layer is below z/r = 0.2, under the assumption of LTE, and of a geometrically thin emitting layer. Nevertheless, the low excitation temperature could also result from the drastically different emitting radii of the two c-C₃H₂ transitions, while the rotational diagram assumes a uniform temperature. Another explanation could be that the c-C₃H₂ is actually originating from a higher layer, but the LTE approximation is not valid anymore and the resulting excitation temperature is lower. On the other hand, both high and low emitting layers for the c-C₃H₂ molecule have been presented in other works: for example, Cleeves et al. (2021) extracted 0.25 < z/r < 0.4 for the TW Hya disk, while Paneque-Carreño et al. (2023) obtained z/r ~ 0.1 for the HD 163296 disk. More observations and detailed chemical modeling are needed to shed light on the emitting layers of these molecules to constrain their formation pathways.

Fig. 10

Channel maps around 4 km s⁻¹ showing the spatial distribution of the emission of three different molecules: SO (stacked), ¹²CO and H₂CO, from left to right. The white contour shows the 3σ level of the SO emission, which is reported in the other channel maps for reference. The pink contours indicate the bright ring in the band 7 continuum emission (Isella et al. 2019), while the two white dots mark the position of the two protoplanets (Wang et al. 2021). The white ellipse in the bottom left side of each panel is the beam of each imaged cube.

5.6 No accretion signatures detected in SO

The sulfur monoxide has been previously linked to a variety of processes occurring in planet forming disks (Sakai et al. 2014; Tabone et al. 2017), including planetary accretion (Booth et al. 2023a). Nevertheless, the sulfur monoxide has not been detected in many class II disks, except for the disk around AB Aurigae (Pacheco-Vázquez et al. 2016), Oph-IRS 48 (Booth et al. 2021), HD 100546 (Booth et al. 2023a), and HD 169142 (Booth et al. 2023b; Law et al. 2023).

In the PDS 70 disk, we did not detect any signal from the SO molecule either localized onto the two actively accreting protoplanets or as an extended emission across the disk, after applying boosting techniques both in the uv and in the image plane. We extracted an upper limit on the disk averaged column density of SO, from the uv-stacked data cube, obtaining N_t < 2.7 × 10¹² cm⁻², assuming T_ex = 30 K. Even if we rescale for the different disk masses (as N_t ∝ disk mass) of the sources where SO has been detected and we consider the higher excitation temperatures they find (around 100 K), the upper limit we extracted for the SO total column density is lower than the SO column density obtained for the sources listed above.

Additionally, we highlight the marginal detection of a localized 3σ emission (white contour in the first panel of Fig. 10) in the channel map of the uv-stacked SO data cube, located at ~ $1_{\cdot}^{''} 2$ $\[1^{\prime\prime}_\cdot 2\]$ outside the mm dust ring (pink contours), which is spatially and spectrally colocated with a kink in the ¹²CO Keplerian pattern (see the second panel) and with a gap in the H₂CO emission (see the third panel). Even if this detection is tentative, it would be interesting to further investigate a possible physical origin of this feature in high energy lines. Interestingly, the radial location of the observed SO marginal signature and of the H₂CO gap corresponds to the radial position of the peak intensity of CS and C₂H (see Fig. D.2), and to a dip in the brightness temperature profile of ¹²CO (Law et al. 2024). In this context, Jiang et al. (2023) showed that a giant planet forming inside a gap and outside the methane snowline can result in an enhanced gas-phase C/O ratio inside the gap.

5.7 C/O ratio

The molecular line emission of the observed transitions is dominated by C-rich molecules, including various simple organic molecules, such as HCN, C₂H, CS, and c-C₃H₂, while we did not detect CH₃OH and SO, which may suggest a gas phase O-depleted outer disk. This evidence is fundamental in the context of two important astrochemical open questions related to the disk evolution, which are the understanding of the sulfur chemistry (main carriers, phase and evolution, Kama et al. 2019; Öberg & Bergin 2021) and the fact that we observe a general depletion of volatile oxygen with respect to carbon in more evolved disks (Miotello et al. 2019, 2023).

In this context, there are several pieces of evidence pointing to an O-poor gas reservoir in the PDS 70 outer disk, and possibly to a higher C/O ratio, including:

bright emission of C₂H, which is enhanced in volatile O-depleted regions, preventing the destruction of the hydrocarbon (Bergin et al. 2016; Miotello et al. 2019; Bosman et al. 2021; Alarcón et al. 2021), and the detection of c-C₃H₂;
bright CS emission outside the submillimeter continuum ring, and the non-detection of SO, which has been previously linked to high C/O ratio (Guilloteau et al. 2016; Semenov et al. 2018; Le Gal et al. 2021; Booth et al. 2021; Keyte et al. 2023);
non-detection of CH₃OH.

Cridland et al. (2023) compared the observed radial profiles of ¹²CO, C¹⁸O, and C₂H with the modeled ones obtained trough the thermo-chemical code DALI (Bruderer et al. 2012; Bruderer 2013), in order to extract the volatile C/O ratio for the PDS 70 disk (we note that they inferred a star C/O ratio for PDS 70 of ~0.44). The preferred model they obtained assumes a volatile C/O ratio marginally above unity, even if they suggest that radial variations of the C/O ratio across the disk may better reproduce both the inner peak and the outer emission shoulder in the observed radial profile of the C₂H integrated intensity.

We can compare this result with the lower limit on the CS/SO ratio, which has been proposed to be a good direct tracer of the elemental C/O ratio (Guilloteau et al. 2016; Semenov et al. 2018). In particular, we extracted a lower limit on the CS/SO ratio, by taking the ratio of the CS average column density within a circle of radius 2″ and the upper limit on the SO column density we obtained from the uv-stacked data cube. By varying the excitation temperature between 20 and 50 K, we obtained a lower limit on the average CS/SO ratio between 0.8 and 1.1. Semenov et al. (2018) explored the effect of the C/O ratio on the CS/SO one, through thermo-chemical models, obtaining that an increase in the volatile C/O ratio from 0.46 to 1.0 would result in raising the CS/SO ratio from ~0.026 to ~13. The inferred lower limit on the disk-averaged CS/SO ratio suggests that the gas phase C/O ratio in the outer disk may become greater than 1. With the cavity wall presenting a C/O ratio of 1 (Cridland et al. 2023), this could be potentially indicating a radial gradient in the C/O ratio of the gas phase. In the future, a forward model extracting the gas C/O ratio from the very inner disk from the JWST data (Perotti et al. 2023) to the outer disk from the ALMA data presented in this paper should have high priority.

6 Conclusions

In this work, we have presented high-resolution and deep-sensitivity ALMA band 6 line emission observations of the planet-hosting disk around PDS 70. In particular, we have presented an astrochemical overview of the molecular complexity of the PDS 70 disk, and searched for weak lines through imageplane and uv-plane boosting techniques. We summarize our conclusions below:

We robustly detected ring-shaped ¹²CO, ¹³CO, C¹⁸O, H¹³CN, HC¹⁵N, DCN, H₂CO, C₂H, CS, and H¹³CO⁺, with the peak of the emission not always colocated with the millimeter continuum bright ring. Most of the molecules show a peak in the emission inside the millimeter dust ring, which we suggest is linked to the efficient irradiation of the cavity wall outside the planets’ orbit on the disk chemistry;
We present the first detection of the c-C₃H₂ 218 GHz lines in the PDS 70 disk, leveraging boosting techniques;
We suggest that the outer disk around PDS 70 has a O-depleted gas reservoir, supported by the non-detection of SO and CH₃OH, the detection of C-rich molecules (CS, C₂H, HCN, c-C₃H₂), and a CS/SO ≳ 1 ratio;
We extracted a radial profile of the excitation temperature and column density of H₂CO, from a radially resolved rotational diagram analysis. We propose a gas-phase formation pathway for the observed H₂CO emission, even if the grain-surface formation route cannot be completely ruled out;
We have found evidence of an azimuthal asymmetry in the higher-energy H₂CO (3_2,1−2_2,0) line emission, which could be linked to accretion heating near PDS 70b. Moreover, we marginally detected a spatially and spectrally localized 3σ-signature in the SO emission, from the outer disk, colocated with a kink in the ¹²CO and a gap in the H₂CO peak intensity.

Detailed chemical modeling is needed to quantify the C/O ratio across the dust cavity, shedding light on the composition of the gas accreted by the two giant planets, but also to test the suggested role of the irradiated cavity wall in shaping the molecular emission.

Acknowledgements

We thank the referee for an exceptionally thorough review of the paper and very helpful comments. We are grateful to Ryan Loomis and Ian Czekala for the useful suggestions. This paper makes use of the following ALMA data: 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), 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. L.R., S.F., and M.L. are funded by the European Union (ERC, UNVEIL, 101076613). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. S.F. also acknowledges financial contribution from PRIN-MUR 2022YP5ACE. Support for C.J.L. was provided by NASA through the NASA Hubble Fellowship grant No. HST-HF2-51535.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. The PI acknowledges the use of computing resources from the Italian node of the European ALMA Regional Center, hosted by INAF-Istituto di Radioastronomia.

Appendix A Disk-integrated spectra

Fig. A.1

Effect of the Keplerian shifting on the line spectrum. Top row: Disk-integrated spectrum of ¹²CO (J=2–1), C₂H (3_7/2−2_5/2) and H₂CO (3_2,1 − 2_2,0), respectively, on a circular de-projected region of radius 3″ for the ¹²CO and of 2″ for the other ones. The blue line is the shifted spectrum, obtained by correcting for the Doppler shift due to the Keplerian rotation of the disk, while the gray line corresponds to the un-shifted spectrum. Bottom row: Spectra of ¹²CO (J=2–1), C₂H (3_7/2−2_5/2) and H₂CO (3_2,1−2_2,0) as a function of radius, obtained after the Keplerian shifting: the spectra are expressed in brightness temperature, using the Rayleigh-Jeans approximation. The hyperfine components of the C₂H transition are visible in the middle column.

Fig. A.2

Spectra as a function of radius for the individual undetected lines SO 6₇−5₆ and 6₆−5₅, compared to the uv-stacked one. The white contours refer to the 3K level from the ¹²CO (J=2–1).

We present the stacked spectra obtained leveraging the Keplerian rotation of the disk in Fig. A.1, as we described in Sect. 4.3 (Teague et al. 2016). The top row in Fig. A.1 shows an example of three spectra for ¹²CO (2–1), C₂H (3_7/2−2_5/2) and H₂CO (3_2,1−2_2,0) transitions respectively: the shifted spectrum (in blue) visibly increases the S/N with respect to the original one (in gray), which shows the typical double horn effect due to the Doppler shift. The second column in Fig. A.1 shows the spectrally resolved hyperfine components of the C₂H transition, with the double horn of the hyperfine components being blended in the Doppler-shifted spectrum. We highlight the spatial distribution of the line emission by extracting the spatially integrated spectra as a function of radius, averaging over $0_{\cdot}^{''} 1$ $\[0^{\prime\prime}_\cdot 1\]$ wide annuli, after applying the Keplerian shifting: the result is presented in the bottom row of Fig. A.1, for the ¹²CO (2–1), C₂H (3_7/2−2_5/2), and H₂CO (3_2,1 − 2_2,0) lines. Figure A.2 shows the Keplerian shifted spectra as a function of radius, for SO 6₇−5₆, 6₆−5₅, and for the uv-stacked cube. The white contours refer to the 3K level from the ¹²CO (2–1), which highlights the absence of a clear Keplerian emission.

Appendix B Channel maps of the faint lines

In Fig. B.1 we present the channel maps of the detected c-C₃H₂ (218 GHz) and C-C₃H₂ (244 GHz) lines, and of the marginally detected SO, after the uv-stacking. The marginal SO signature in the outer disk, which we discussed in Sec. 5.6, is highlighted by the 3σ white contour and by the 4σ dot, at a velocity channel of 4.10 km s⁻¹, in the last row. The channel maps of the faint c-C₃H₂ lines (first and middle rows) show 3σ features at different velocity channels. The marginal detection of the signal from these faint lines throughout the channel maps highlights the importance of using boosting techniques, such as the ones we presented in this work, in order to asses the robustness of the detection.

Fig. B.1

Channel maps of c-C₃H₂ (218 GHz, top row), c-C₃H₂ (244 GHz, middle row) and SO stacked lines (bottom row), with the white contours showing the 3σ levels and the pink contours which indicate the extent of the bright ring in the band 7 continuum emission (Isella et al. 2019). The white dots correspond to the position of the two protoplanets (Wang et al. 2021).

Appendix C Test on uv-stacking

We tested the effect of the uv-stacking, we presented in Sec. 3.3, on the two individually detected lines of H₂CO: the result is presented in Fig. C.1. As expected from the effect of the CASA task concat (see Sec. 3.3 for more details), the integrated intensity profile of the stacked cube (pink line) lies between the higher and the lower profiles, representing the integrated intensity profiles of the bright (orange line) and the faint (blue line) transitions, respectively.

Fig. C.1

Test of the uv-stacking on the two individually detected lines of H₂CO, showing the radial profiles of the integrated intensity from the stacked cube (pink), which lies between the bright (orange) and faint (blue) ones, as expected.

Appendix D Gaussian fits to the radial profiles

As we described in Sec. 5.1, we performed a multiple-Gaussian fit (triple for¹²CO and H₂CO, double for the other molecules) of the radial profiles of the molecular integrated intensity, in order to retrieve the radial location of the peak of the emission, and related uncertainty. We performed the fit procedure using scipy.optimize.curve_fit, and the result is plotted in Fig. D.1 (dashed black line), together with the single Gaussians (dashed orange, pink, and red lines). The fit is in very good agreement with the radial profiles extracted from observations (solid blue lines).

We also tested the multiple-Gaussian fitting to the peak intensity radial profiles (see Fig. D.2), and we listed the peak positions with associated uncertainties in Table D.1. In Fig. D.3 we compare the peak positions extracted from the integrated intensity (marked in white) and from the peak intensity radial profiles (marked in red). There is good agreement between the two methods, except for the C₂H and ¹³CO molecules.

Fig. D.1

Radial profiles of the integrated intensity for the analyzed molecular transitions (blue lines), and corresponding multiple-Gaussian fit (dashed black line). Dashed orange, pink, and red lines show single Gaussian, which summed together result in the dashed black profiles. The vertical dashed blue line corresponds to the peak position directly inferred from the solid blue profile.

Fig. D.2

Same as Fig. D.1 but for the peak intensity.

Table D.1

Radial position of the peak of the molecular emission.

Fig. D.3

Same as Fig. 7, with a comparison between the radial position of the peak of the emission extracted from the radial profiles of the integrated intensity (marked in white), and from the radial profiles of the peak intensity (marked in red). Peaks and associated uncertainties were obtained from multiple-Gaussian fitting (see Fig. D.1 and D.2).

Appendix E Peak intensity maps

We present the peak intensity maps of the same lines as in Fig. 1, using the quadratic method implemented in the bettermoments tool (Teague & Foreman-Mackey 2019). We extracted the maps from the JvM-corrected and continuum-subtracted cubes imaged with natural weighting. The result is presented in Fig. E.1. We present the result in K, assuming the Rayleigh-Jeans approximation.

The structured emission morphology is visible in the peak intensity maps, as in the integrated intensity maps, with the peak of the emission not always colocated with the submillimeter dust ring. Most of the molecules show axisymmetric emission, except for the 2.5σ asymmetry along the major axis in H₂CO (3_2,1−2_2,0). The outer shoulders of emission of HC¹⁵N, DCN, H₂CO and C₂H, outside the outer edge of the submillimeter ring, are more evident in the peak intensity maps, with respect to the integrated intensity maps, and they can be clearly seen in the corresponding radial profiles, in Fig. D.2. The optically thick ¹²CO emission shows an arced feature along the major axis in the peak intensity map, highlighting the elevated emitting layer.

The radial dark features are artifacts due to sparse channelization: for reference, these are particularly visible in the maps of CS and H₂CO (3_0,3−2_0,2). The azimuthal arced features could also be linked to physical and chemical processes affecting the molecular emission: these are particularly evident in the maps of C₂H and H₂CO (3_0,3−2_0,2), along the major axis, outside the bright continuum ring.

Fig. E.1

Peak intensity maps of the detected lines, with the white contours showing the bright ring in the submillimeter continuum emission (Isella et al. 2019) and the white dots marking the position of the two forming planets (Wang et al. 2021). The brightness temperature is derived in the Rayleigh-Jeans regime. The white ellipse on the bottom left of each panel is the beam.

References

Ahrens, V., & Winnewisser, G. 1999, Zeitschrift für Naturforschung A, 54, 131 [CrossRef] [Google Scholar]
Alarcón, F., Bosman, A. D., Bergin, E. A., et al. 2021, ApJS, 257, 8 [CrossRef] [Google Scholar]
Amano, T., & Hirota, E. 1974, J. Mol. Spectrosc., 53, 346 [CrossRef] [Google Scholar]
Andrews, S. M. 2020, ARA&A, 58, 483 [Google Scholar]
Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41 [NASA ADS] [CrossRef] [Google Scholar]
Aoyama, Y., & Ikoma, M. 2019, ApJ, 885, L29 [NASA ADS] [CrossRef] [Google Scholar]
Atkinson, R., Baulch, D. L., Cox, R. A., et al. 2006, Atmos. Chem. Phys., 6, 3625 [NASA ADS] [CrossRef] [Google Scholar]
Bae, J., Teague, R., Andrews, S. M., et al. 2022, ApJ, 934, L20 [NASA ADS] [CrossRef] [Google Scholar]
Bae, J., Isella, A., Zhu, Z., et al. 2023, ASP Conf. Ser., 534, 423 [NASA ADS] [Google Scholar]
Benisty, M., Bae, J., Facchini, S., et al. 2021, ApJ, 916, L2 [NASA ADS] [CrossRef] [Google Scholar]
Benisty, M., Dominik, C., Follette, K., et al. 2023, ASP Conf. Ser., 534, 605 [NASA ADS] [Google Scholar]
Bergin, E. A., Du, F., Cleeves, L. I., et al. 2016, ApJ, 831, 101 [Google Scholar]
Bogey, M., Demuynck, C., & Destombes, J. 1981, Chem. Phys. Lett., 81, 256 [NASA ADS] [CrossRef] [Google Scholar]
Bogey, M., Demuynck, C., & Destombes, J. 1986, Chem. Phys. Lett., 125, 383 [NASA ADS] [CrossRef] [Google Scholar]
Booth, A. S., van der Marel, N., Leemker, M., van Dishoeck, E. F., & Ohashi, S. 2021, A&A, 651, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Booth, A. S., Ilee, J. D., Walsh, C., et al. 2023a, A&A, 669, A53 [Google Scholar]
Booth, A. S., Law, C. J., Temmink, M., Leemker, M., & Macías, E. 2023b, A&A, 678, A146 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bosman, A. D., Alarcón, F., Bergin, E. A., et al. 2021, ApJS, 257, 7 [NASA ADS] [CrossRef] [Google Scholar]
Bruderer, S. 2013, A&A, 559, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Bruderer, S., van Dishoeck, E. F., Doty, S. D., & Herczeg, G. J. 2012, A&A, 541, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Brünken, S., Fuchs, U., Lewen, F., et al. 2004, J. Mol. Spectrosc., 225, 152 [CrossRef] [Google Scholar]
Carney, M. T., Hogerheijde, M. R., Loomis, R. A., et al. 2017, A&A, 605, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Casassus, S., & Cárcamo, M. 2022, MNRAS, 513, 5790 [NASA ADS] [CrossRef] [Google Scholar]
Casassus, S., & Pérez, S. 2019, ApJ, 885, L27 [NASA ADS] [CrossRef] [Google Scholar]
Cataldi, G., Yamato, Y., Aikawa, Y., et al. 2021, ApJS, 257, 10 [NASA ADS] [CrossRef] [Google Scholar]
Cazzoletti, P., van Dishoeck, E. F., Visser, R., Facchini, S., & Bruderer, S. 2018, A&A, 609, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Cazzoli, G., & Puzzarini, C. 2005, J. Mol. Spectrosc., 233, 280 [NASA ADS] [CrossRef] [Google Scholar]
Cazzoli, G., Puzzarini, C., & Lapinov, A. V. 2004, ApJ, 611, 615 [NASA ADS] [CrossRef] [Google Scholar]
Chardon, J., & Guichon, D. 1973, J. Phys., 34, 791 [CrossRef] [EDP Sciences] [Google Scholar]
Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102 [Google Scholar]
Christiaens, V., Casassus, S., Absil, O., et al. 2019, MNRAS, 486, 5819 [NASA ADS] [CrossRef] [Google Scholar]
Cieza, L. A., González-Ruilova, C., Hales, A. S., et al. 2021, MNRAS, 501, 2934 [Google Scholar]
Clark, W. W., & De Lucia, F. C. 1976, J. Mol. Spectrosc., 60, 332 [CrossRef] [Google Scholar]
Cleeves, L. I., Bergin, E. A., & Harries, T. J. 2015, ApJ, 807, 2 [NASA ADS] [CrossRef] [Google Scholar]
Cleeves, L. I., Loomis, R. A., Teague, R., et al. 2021, ApJ, 911, 29 [NASA ADS] [CrossRef] [Google Scholar]
Cridland, A. J., Facchini, S., van Dishoeck, E. F., & Benisty, M. 2023, A&A, 674, A211 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Currie, T., Biller, B., Lagrange, A., et al. 2023, ASP Conf. Ser., 534, 799 [NASA ADS] [Google Scholar]
Czekala, I., Loomis, R. A., Teague, R., et al. 2021, ApJS, 257, 2 [NASA ADS] [CrossRef] [Google Scholar]
Facchini, S., Pinilla, P., van Dishoeck, E. F., & de Juan Ovelar, M. 2018, A&A, 612, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Facchini, S., van Dishoeck, E. F., Manara, C. F., et al. 2019, A&A, 626, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Facchini, S., Teague, R., Bae, J., et al. 2021, AJ, 162, 99 [NASA ADS] [CrossRef] [Google Scholar]
Fedoseev, G., Cuppen, H. M., Ioppolo, S., Lamberts, T., & Linnartz, H. 2015, MNRAS, 448, 1288 [NASA ADS] [CrossRef] [Google Scholar]
Fockenberg, C., & Preses, J. M. 2002, J. Phys. Chem. A, 106, 2924 [NASA ADS] [CrossRef] [Google Scholar]
Foreman-Mackey, D., Conley, A., Meierjurgen Farr, W., et al. 2013, Astrophysics Source Code Library [record ascl:1303.002] [Google Scholar]
Fuchs, U., Bruenken, S., Fuchs, G. W., et al. 2004, Zeitschrift Naturforschung Teil A, 59, 861 [NASA ADS] [CrossRef] [Google Scholar]
Fuchs, G. W., Cuppen, H. M., Ioppolo, S., et al. 2009, A&A, 505, 629 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Fulton, B. J., Rosenthal, L. J., Hirsch, L. A., et al. 2021, ApJS, 255, 14 [NASA ADS] [CrossRef] [Google Scholar]
Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Garg, H., Pinte, C., Hammond, I., et al. 2022, MNRAS, 517, 5942 [CrossRef] [Google Scholar]
Goldsmith, P. F., & Langer, W. D. 1999, ApJ, 517, 209 [Google Scholar]
Greene, T. P., Line, M. R., Montero, C., et al. 2016, ApJ, 817, 17 [NASA ADS] [CrossRef] [Google Scholar]
Gregersen, E. M., & Evans, Neal J., I. 2001, ApJ, 553, 1042 [NASA ADS] [CrossRef] [Google Scholar]
Guillot, T., Fletcher, L. N., Helled, R., et al. 2023, ASP Conf. Ser., 534, 947 [Google Scholar]
Guilloteau, S., Reboussin, L., Dutrey, A., et al. 2016, A&A, 592, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Guzmán, V. V., Bergner, J. B., Law, C. J., et al. 2021, ApJS, 257, 6 [CrossRef] [Google Scholar]
Haffert, S. Y., Bohn, A. J., de Boer, J., et al. 2019, Nat. Astron., 3, 749 [Google Scholar]
Hashimoto, J., Tsukagoshi, T., Brown, J. M., et al. 2015, ApJ, 799, 43 [CrossRef] [Google Scholar]
Heays, A. N., Visser, R., Gredel, R., et al. 2014, A&A, 562, A61 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Herbst, E., Messer, J., De Lucia, F. C., & Helminger, P. 1984, J. Mol. Spectrosc., 108, 42 [NASA ADS] [CrossRef] [Google Scholar]
Hidaka, H., Watanabe, N., Shiraki, T., Nagaoka, A., & Kouchi, A. 2004, ApJ, 614, 1124 [Google Scholar]
Hily-Blant, P., Magalhaes de Souza, V., Kastner, J., & Forveille, T. 2019, A&A, 632, L12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Hinkley, S., Carter, A., Ray, S., et al. 2023, AAS Meeting Abstracts, 55, 247.05 [Google Scholar]
Hiraoka, K., Ohashi, N., Kihara, Y., et al. 1994, Chem. Phys. Lett., 229, 408 [Google Scholar]
Hiraoka, K., Sato, T., Sato, S., et al. 2002, ApJ, 577, 265 [Google Scholar]
Högbom, J. A. 1974, A&AS, 15, 417 [Google Scholar]
Huang, J., Öberg, K. I., Qi, C., et al. 2017, ApJ, 835, 231 [Google Scholar]
Ilee, J. D., Walsh, C., Booth, A. S., et al. 2021, ApJS, 257, 9 [NASA ADS] [CrossRef] [Google Scholar]
Isella, A., Benisty, M., Teague, R., et al. 2019, ApJ, 879, L25 [Google Scholar]
Izquierdo, A. F., Facchini, S., Rosotti, G. P., van Dishoeck, E. F., & Testi, L. 2022, ApJ, 928, 2 [NASA ADS] [CrossRef] [Google Scholar]
Izquierdo, A. F., Testi, L., Facchini, S., et al. 2023, A&A, 674, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Jiang, H., Wang, Y., Ormel, C. W., Krijt, S., & Dong, R. 2023, A&A, 678, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Jorsater, S., & van Moorsel, G. A. 1995, AJ, 110, 2037 [Google Scholar]
Kama, M., Shorttle, O., Jermyn, A. S., et al. 2019, ApJ, 885, 114 [Google Scholar]
Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Keppler, M., Teague, R., Bae, J., et al. 2019, A&A, 625, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Keyte, L., Kama, M., Booth, A. S., et al. 2023, Nat. Astron., 7, 684 [NASA ADS] [CrossRef] [Google Scholar]
Klapper, G., Lewen, F., Gendriesch, R., Belov, S., & Winnewisser, G. 2000, J. Mol. Spectrosc., 201, 124 [NASA ADS] [CrossRef] [Google Scholar]
Lattanzi, V., Walters, A., Drouin, B. J., & Pearson, J. C. 2007, ApJ, 662, 771 [NASA ADS] [CrossRef] [Google Scholar]
Law, C. J., Loomis, R. A., Teague, R., et al. 2021a, ApJS, 257, 3 [NASA ADS] [CrossRef] [Google Scholar]
Law, C. J., Teague, R., Loomis, R. A., et al. 2021b, ApJS, 257, 4 [NASA ADS] [CrossRef] [Google Scholar]
Law, C. J., Booth, A. S., & Öberg, K. I. 2023, ApJ, 952, L19 [NASA ADS] [CrossRef] [Google Scholar]
Law, C. J., Benisty, M., Facchini, S., et al. 2024, ApJ, 964, 190 [NASA ADS] [CrossRef] [Google Scholar]
Le Gal, R., Öberg, K. I., Teague, R., et al. 2021, ApJS, 257, 12 [NASA ADS] [CrossRef] [Google Scholar]
Leemker, M., Booth, A. S., van Dishoeck, E. F., et al. 2022, A&A, 663, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Lees, R., & Baker, J. 1968, J. Chem. Phys., 48, 5299 [NASA ADS] [CrossRef] [Google Scholar]
Lissauer, J. J., Batalha, N. M., & Borucki, W. J. 2023, ASP Conf. Ser., 534, 839 [Google Scholar]
Loison, J.-C., Agúndez, M., Wakelam, V., et al. 2017, MNRAS, 470, 4075 [Google Scholar]
Long, F., Pinilla, P., Herczeg, G. J., et al. 2018a, ApJ, 869, 17 [Google Scholar]
Long, Z. C., Akiyama, E., Sitko, M., et al. 2018b, ApJ, 858, 112 [NASA ADS] [CrossRef] [Google Scholar]
Loomis, R. A., Cleeves, L. I., Öberg, K. I., Guzman, V. V., & Andrews, S. M. 2015, ApJ, 809, L25 [Google Scholar]
Loomis, R. A., Oberg, K. I., Andrews, S. M., et al. 2018, Astrophysics Source Code Library [record ascl:1802.006] [Google Scholar]
Lovas, F. J., Suenram, R., Ogata, T., & Yamamoto, S. 1992, ApJ, 399, 325 [CrossRef] [Google Scholar]
Maiwald, F., Lewen, F., Ahrens, V., et al. 2000, J. Mol. Spectrosc., 1, 166 [NASA ADS] [CrossRef] [Google Scholar]
McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, ASP Conf. Ser., 376, 127 [Google Scholar]
Mesa, D., Keppler, M., Cantalloube, F., et al. 2019, A&A, 632, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Millar, T. J., Bennett, A., & Herbst, E. 1989, ApJ, 340, 906 [NASA ADS] [CrossRef] [Google Scholar]
Miotello, A., Facchini, S., van Dishoeck, E. F., et al. 2019, A&A, 631, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Miotello, A., Kamp, I., Birnstiel, T., Cleeves, L. C., & Kataoka, A. 2023, ASP Conf. Ser., 534, 501 [NASA ADS] [Google Scholar]
Mollière, P., van Boekel, R., Dullemond, C., Henning, T., & Mordasini, C. 2015, ApJ, 813, 47 [Google Scholar]
Muñoz-Romero, C. E., Öberg, K. I., Law, C. J., et al. 2023, ApJ, 943, 35 [CrossRef] [Google Scholar]
Müller, H. S. P., & Lewen, F. 2017, J. Mol. Spectrosc., 331, 28 [CrossRef] [Google Scholar]
Müller, H. S. P., Klaus, T., & Winnewisser, G. 2000, A&A, 357, L65 [NASA ADS] [Google Scholar]
Müller, H. S. P., Thorwirth, S., Roth, D. A., & Winnewisser, G. 2001, A&A, 370, L49 [Google Scholar]
Müller, H. S., Schlöder, F., Stutzki, J., & Winnewisser, G. 2005, J. Mol. Struct., 742, 215 [Google Scholar]
Noble, J. A., Theule, P., Mispelaer, F., et al. 2012, A&A, 543, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Nomura, H., Tsukagoshi, T., Kawabe, R., et al. 2021, ApJ, 914, 113 [NASA ADS] [CrossRef] [Google Scholar]
Nomura, H., Furuya, K., Cordiner, M. A., et al. 2023, ASP Conf. Ser., 534, 1075 [NASA ADS] [Google Scholar]
Öberg, K. I., & Bergin, E. A. 2021, Phys. Rep., 893, 1 [Google Scholar]
Öberg, K. I., Guzmán, V. V., Merchantz, C. J., et al. 2017, ApJ, 839, 43 [Google Scholar]
Öberg, K. I., Cleeves, L. I., Bergner, J. B., et al. 2021a, AJ, 161, 38 [Google Scholar]
Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021b, ApJS, 257, 1 [CrossRef] [Google Scholar]
Öberg, K. I., Facchini, S., & Anderson, D. E. 2023, ARA&A, 61, 287 [CrossRef] [Google Scholar]
Pacheco-Vázquez, S., Fuente, A., Baruteau, C., et al. 2016, A&A, 589, A60 [Google Scholar]
Padovani, M., Walmsley, C., Tafalla, M., Galli, D., & Müller, H. 2009, A&A, 505, 1199 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Paneque-Carreño, T., Miotello, A., van Dishoeck, E. F., et al. 2023, A&A, 669, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
Perez, S., Casassus, S., Ménard, F., et al. 2015a, ApJ, 798, 85 [Google Scholar]
Perez, S., Dunhill, A., Casassus, S., et al. 2015b, ApJ, 811, L5 [NASA ADS] [CrossRef] [Google Scholar]
Perotti, G., Christiaens, V., Henning, T., et al. 2023, Nature, 620, 516 [NASA ADS] [CrossRef] [Google Scholar]
Pickett, H., Cohen, E., Brinza, D., & Schaefer, M. 1981, J. Mol. Spectrosc., 89, 542 [NASA ADS] [CrossRef] [Google Scholar]
Pickett, H., Poynter, R., Cohen, E., et al. 1998, J. Quant. Spectrosc. Rad. Trans., 60, 883 [CrossRef] [Google Scholar]
Pinilla, P., Benisty, M., & Birnstiel, T. 2012, A&A, 545, A81 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Pinte, C., Price, D. J., Ménard, F., et al. 2018, ApJ, 860, L13 [Google Scholar]
Pinte, C., van der Plas, G., Ménard, F., et al. 2019, Nat. Astron., 3, 1109 [NASA ADS] [CrossRef] [Google Scholar]
Pinte, C., Teague, R., Flaherty, K., et al. 2023, ASP Conf. Ser., 534, 645 [NASA ADS] [Google Scholar]
Portilla-Revelo, B., Kamp, I., Facchini, S., et al. 2023, A&A, 677, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Rosotti, G. P., Ilee, J. D., Facchini, S., et al. 2021, MNRAS, 501, 3427 [Google Scholar]
Roueff, E., Loison, J. C., & Hickson, K. M. 2015, A&A, 576, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Sakai, N., Sakai, T., Hirota, T., et al. 2014, Nature, 507, 78 [Google Scholar]
Sanchis, E., Testi, L., Natta, A., et al. 2021, A&A, 649, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Sastry, K., Helminger, P., Charo, A., Herbst, E., & De Lucia, F. 1981, ApJ, 251, L119 [NASA ADS] [CrossRef] [Google Scholar]
Sastry, K., Lees, R., & De Lucia, F. 1984, J. Mol. Spectrosc., 103, 486 [NASA ADS] [CrossRef] [Google Scholar]
Semenov, D., Favre, C., Fedele, D., et al. 2018, A&A, 617, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Shirley, Y. L. 2015, PASP, 127, 299 [Google Scholar]
Sierra, A., Pérez, L. M., Zhang, K., et al. 2021, ApJS, 257, 14 [NASA ADS] [CrossRef] [Google Scholar]
Stadler, J., Benisty, M., Izquierdo, A., et al. 2023, A&A, 670, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Suzuki, D., Bennett, D. P., Sumi, T., et al. 2016, ApJ, 833, 145 [Google Scholar]
Tabone, B., Cabrit, S., Bianchi, E., et al. 2017, A&A, 607, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Teague, R. 2019, J. Open Source Softw., 4, 1632 [NASA ADS] [CrossRef] [Google Scholar]
Teague, R. 2020, https://doi.org/10.5281/zenodo.4321137 [Google Scholar]
Teague, R., & Foreman-Mackey, D. 2019, Astrophysics Source Code Library [record ascl:1901.009] [Google Scholar]
Teague, R., Guilloteau, S., Semenov, D., et al. 2016, A&A, 592, A49 [CrossRef] [EDP Sciences] [Google Scholar]
Teague, R., Bae, J., & Bergin, E. A. 2019, Nature, 574, 378 [NASA ADS] [CrossRef] [Google Scholar]
Temmink, M., Booth, A. S., van der Marel, N., & van Dishoeck, E. F. 2023, A&A, 675, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Terwisscha van Scheltinga, J., Hogerheijde, M. R., Cleeves, L. I., et al. 2021, ApJ, 906, 111 [Google Scholar]
Thompson, S. E., Coughlin, J. L., Hoffman, K., et al. 2018, ApJS, 235, 38 [NASA ADS] [CrossRef] [Google Scholar]
Tiemann, E. 1982, J. Mol. Spectrosc., 91, 60 [CrossRef] [Google Scholar]
van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2016, A&A, 585, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
van’t Hoff, M. L. R., Harsono, D., Tobin, J. J., et al. 2020, ApJ, 901, 166 [Google Scholar]
Vigan, A., Fontanive, C., Meyer, M., et al. 2021, A&A, 651, A72 [EDP Sciences] [Google Scholar]
Visser, R., Bruderer, S., Cazzoletti, P., et al. 2018, A&A, 615, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Vrtilek, J., Gottlieb, C., & Thaddeus, P. 1987, ApJ, 314, 716 [CrossRef] [Google Scholar]
Wagner, K., Follete, K. B., Close, L. M., et al. 2018, ApJ, 863, L8 [Google Scholar]
Walsh, C., Millar, T. J., Nomura, H., et al. 2014, A&A, 563, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Walsh, C., Loomis, R. A., Öberg, K. I., et al. 2016, ApJ, 823, L10 [NASA ADS] [CrossRef] [Google Scholar]
Wang, J. J., Vigan, A., Lacour, S., et al. 2021, AJ, 161, 148 [Google Scholar]
Watanabe, N., & Kouchi, A. 2002, ApJ, 571, L173 [Google Scholar]
Watanabe, N., Nagaoka, A., Shiraki, T., & Kouchi, A. 2004, ApJ, 616, 638 [Google Scholar]
Weaver, E., Isella, A., & Boehler, Y. 2018, ApJ, 853, 113 [Google Scholar]
Winnewisser, M., Winnewisser, B. P., & Winnewisser, G. 1985, ASI Ser. C, 157, 375 [Google Scholar]
Winnewisser, G., Belov, S., Klaus, T., & Schieder, R. 1997, J. Mol. Spectrosc., 184, 468 [NASA ADS] [CrossRef] [Google Scholar]
Xu, L.-H., Fisher, J., Lees, R., et al. 2008, J. Mol. Spectrosc., 251, 305 [NASA ADS] [CrossRef] [Google Scholar]
Yen, H.-W., Koch, P. M., Liu, H. B., et al. 2016, ApJ, 832, 204 [NASA ADS] [CrossRef] [Google Scholar]
Zhang, Y., Snellen, I. A. G., Bohn, A. J., et al. 2021, Nature, 595, 370 [NASA ADS] [CrossRef] [Google Scholar]
Zhou, Y., Bowler, B. P., Wagner, K. R., et al. 2021, AJ, 161, 244 [NASA ADS] [CrossRef] [Google Scholar]

All Tables