Mapping the Vertical Gas Structure of the Planet-hosting PDS 70 Disk

PDS 70 hosts two massive, still-accreting planets and the inclined orientation of its protoplanetary disk presents a unique opportunity to directly probe the vertical gas structure of a planet-hosting disk. Here, we use high-spatial-resolution (≈0.″1; 10 au) observations in a set of CO isotopologue lines and HCO+ J = 4−3 to map the full 2D (r, z) disk structure from the disk atmosphere, as traced by 12CO, to closer to the midplane, as probed by less abundant isotopologues and HCO+. In the PDS 70 disk, 12CO traces a height of z/r ≈ 0.3, 13CO is found at z/r ≈ 0.1, and C18O originates at, or near, the midplane. The HCO+ surface arises from z/r ≈ 0.2 and is one of the few non-CO emission surfaces constrained with high-fidelity in disks to date. In the 12CO J = 3−2 line, we resolve a vertical dip and steep rise in height at the cavity wall, making PDS 70 the first transition disk where this effect is directly seen in line-emitting heights. In the outer disk, the CO emission heights of PDS 70 appear typical for its stellar mass and disk size and are not substantially altered by the two inner embedded planets. By combining CO isotopologue and HCO+ lines, we derive the 2D gas temperature structure and estimate a midplane CO snowline of ≈ 56–85 au. This implies that both PDS 70b and 70c are located interior to the CO snowline and are likely accreting gas with a high C/O ratio of ≈ 1.0, which provides context for future planetary atmospheric measurements from, e.g., JWST, and for properly modeling their formation histories.


INTRODUCTION
Planets assemble and acquire their compositions from gas and dust in their natal protoplanetary disks, while the planet formation process is also expected to simultaneously alter the disk physical and chemical structure (e.g., Cleeves et al. 2015;Facchini et al. 2018;Favre et al. 2019).In particular, vertical gas flows have been identified in multiple disks with suspected planets, which suggests that perturbations in the vertical gas structure of disks may be common planetary signposts (e.g., Teague * NASA Hubble Fellowship Program Sagan Fellow et al. 2019a; Yu et al. 2021;Galloway-Sprietsma et al. 2023;Izquierdo et al. 2023).However, the direct detection of planets embedded in disks remains difficult and only a handful of such systems have thus far been robustly identified (Keppler et al. 2018;Haffert et al. 2019;Currie et al. 2022;Hammond et al. 2023;Wagner et al. 2023), making it difficult to conclusively assess the potential influence of giant planets on disk structure.Of these systems, PDS 70 is the only source whose protoplanetary disk is at a favorable inclination (51 • .7;Keppler et al. 2019) to allow for a direct view of its vertical gas distribution and thus provides a unique opportunity Law et al. to understand if and how forming planets alter disk vertical structure.
The relative proximity and inclined orientation of the PDS 70 disk makes it possible to spatially resolve emission arising from elevated regions above and below the disk midplane using high-angular-resolution observations from the Atacama Large Millimeter/submillimeter Array (ALMA) (e.g., Rosenfeld et al. 2013).This then allows for the direct extraction of emission surfaces of vertically extended molecular lines using similar techniques that have now been applied to a substantial number of mid-inclination disks (e.g., Pinte et al. 2018;Teague et al. 2019b;Rich et al. 2021;Izquierdo et al. 2021;Law et al. 2021aLaw et al. , 2022;;Paneque-Carreño et al. 2022, 2023;Izquierdo et al. 2023;Stapper et al. 2023;Law et al. 2023).Using these methods for the PDS 70 disk, we can map out the line emitting heights and conduct a detailed search for vertical perturbations, such as those driven by the embedded planets PDS 70b and 70c.
By combining multiple optically-thick molecular lines, which trace different heights in the disk, we can also empirically derive the full two-dimensional (2D) temperature distribution (e.g., Dartois et al. 2003;Rosenfeld et al. 2013;Pinte et al. 2018;Law et al. 2021a;Leemker et al. 2022;Law et al. 2023).Having a well-constrained thermal structure is crucial as it allows for the determination of important volatile snowline locations, which have a direct impact on the expected atmospheric composition of embedded planets (e.g., Öberg et al. 2011;Mordasini et al. 2016;Mollière et al. 2022).Moreover, disk gas temperatures are vital inputs to numerical sim-ulations and thermochemical models (e.g., Bae & Zhu 2018;Zhang et al. 2021;Calahan et al. 2021) and are required to infer the radial and vertical locations of molecular species, e.g., large organic molecules, that are otherwise too faint to directly map out (e.g., Ilee et al. 2021).This is especially critical to determine what types of chemical reservoirs are accessible to be accreted by PDS 70b and 70c and to more generally connect planetary atmospheric composition with that of the disk environment.
Here, we extract line emission surfaces in the PDS 70 disk using high-spatial-resolution observations of CO isotopologues and HCO + .In Section 2, we describe the observations, imaging, and emission surface extraction process.We present the derived emission surfaces along with radial and vertical temperature profiles in Section 3. In Section 4, we explore the origins of the observed disk vertical structure and estimate the CO snowline location based on the derived thermal structure.We summarize our conclusions in Section 5.

Observational Details
We present long-baseline (LB) observations associated with the ALMA program (2019.1.01619.S; PI: S. Facchini) and Table 1 lists a summary of observation dates and details.Here, we restrict our analysis to a subset of the lines presented in the chemical inventory of Facchini et al. (2021) based on the short-baseline (SB) data from the same program, namely, the J=2-1 line of the 12 CO, 13 CO, and C 18 O isotopologues, which are sufficiently bright and vertically extended to allow for extraction of their vertical structure.Full details about these SB data, including the spectral set-up, are presented in Facchini et al. (2021) and the LB data for the remaining molecules will be published in forthcoming papers.We also make use of existing archival observations of the 12 CO J=3-2 and HCO + J=4-3 lines in the PDS 70 disk compiled from ALMA programs 2015.1.00888.S (PI: E. Akiyama) and 2017.A.00006.S (PI: M. Keppler), which were presented in Long et al. (2018); Keppler et al. (2019); Facchini et al. (2021).We also include unpublished LB data of the same lines from 2018.A.00030.S (PI: M. Benisty).See Benisty et al. (2021) for additional details about each of these programs.We include HCO + in our analysis, due to its vertically-elevated emission and the availability of existing high-angular-resolution line data.Since HCO + typically shows bright, opticallythick emission in disks (e.g., Booth et al. 2019;Huang et al. 2020), it is thus also useful for mapping the thermal structure of PDS 70 in combination with CO isotopologues.

Self Calibration and Imaging
The three J=2-1 CO isotopologue lines are included in one spectral setup of the Facchini et al. (2021) program.For these Band 6 data, we followed the same selfcalibration procedure, in line with Andrews et al. (2018); Benisty et al. (2021).The self-calibration was performed with CASA v5.8 (McMullin et al. 2007;CASA Team et al. 2022).We first re-aligned the data of the individual Execution Blocks (EBs) to a common disk center with the fixvis and fixplanets tasks.The center was computed by fitting an ellipse to the outer ring (also clearly discernible in the SB data, see Facchini et al. 2021).This approach had previously performed well for the self-calibration of the Band 7 data of the PDS 70 disk by Benisty et al. (2021).After flagging the lines listed in Facchini et al. (2021) within ±15 km s −1 from the line center, we spectrally averaged the data into 250 MHz channels to create pseudo-continuum ms tables.Individual EBs were rescaled in amplitude to provide the same disk flux density, after verifying that phase de-coherence was not artificially reducing it.The SB data were self-calibrated first, and then concatenated Law et al. to the LB data before running a new self-calibration.The tclean model was used as the source model for each step, where tclean was performed over an elliptical mask with a 1. ′′ 5 radius and inclination and position angle as in Facchini et al. (2021).Four rounds of phase-only and one round of amplitude self-calibration were performed on the SB data, and three rounds of phase-only and one of amplitude on the LB data, respectively.Before amplitude calibration, we tcleaned the data with a deep 1σ threshold.Using Briggs weighting with robust=0.5, the synthesized beam was 0. ′′ 15×0.′′ 12 (PA=−82.9• ).The resulting 233 GHz continuum image had a flux density of 58.3 mJy within the CLEAN mask (in agreement with Facchini et al. 2021), with an rms of 12.5 µJy beam −1 .We finally applied the gain solutions to the full spectral data.
We obtained archival Band 7 observations of the 12 CO J=3-2 and HCO + J=4-3 lines from the datasets presented in Benisty et al. (2021), for which we applied the same gain solutions to the line spectral windows.For both the Band 6 and Band 7 lines, we continuumsubtracted the data with the contsub task using a firstorder polynomial.
We imaged the five lines discussed in this paper with CASA v6.2.The channel maps were restored using tclean with a threshold of 3.5σ.Table 2 lists the imaging parameters of all cubes.We used Keplerian masks for the channels generated from the keplerian mask code (Teague 2020), where we assumed the disk to be extended 3 ′′ in radius, and we conservatively assumed that all emission is as elevated as 12 CO.Given the focus of the paper, extracting accurate brightness temperatures is imperative.We thus applied the so-called JvM correction to the final images (Jorsater & van Moorsel 1995;Czekala et al. 2021) to more accurately recover the low surface brightness intensities in the disk outer regions, and of weaker lines (e.g., C 18 O J=2-1).The effective rms, which can be underestimated when applying the JvM correction (Casassus & Cárcamo 2022), is not relevant for the analysis of this paper.The full image cube channel maps are availabe in Appendix A. We also imaged the non-continuum-subtracted line data by adopting the same imaging parameters as the line only emission image cubes.The non-continuum-subtracted image cubes are required for the calculation of gas temperatures (Section 3.4).Overall, the ability to extract line emitting heights depends on having sufficiently high angular resolution, spectral resolution, and line sensitivities.Considering this, we selected the imaging parameters that best suited surface extraction for this work.Table 2 lists the properties of the image cubes selected for analysis of the vertical line emission structure.
Figure 1 shows an overview the CO isotopologue and HCO + line emission velocity-integrated intensity, or "zeroth moment," maps, and the 344 GHz continuum emission.The continuum image is taken from Isella et al. (2019).We generated zeroth moment maps of line emission from the image cubes using bettermoments (Teague & Foreman-Mackey 2018) with the same Keplerian masks employed during CLEANing and with no flux threshold for pixel inclusion to ensure accurate flux recovery.

Emission Surface Extraction
We derived vertical emission heights from the line emission image cubes using the disksurf (Teague et al. 2021) Python code, which is based on the methodology originally presented in Pinte et al. (2018).We closely followed the methods outlined in Law et al. (2021a), which we briefly summarize below.3. The solid lines show the radial range used in the fitting, while the dashed lines are extrapolations.Lines of constant z/r from 0.1 to 0.5 are shown in gray.Vertical substructures and the locations of PDS 70b and 70c are marked in each panel (Wang et al. 2021).The lack of 12 CO J=2-1 heights in the inner 50 au (square box) is due to insufficient angular resolution rather than the absence of emission at these radii.The FWHM of the major axis of the synthesized beam is shown in the bottom right corner of each panel.The emission surfaces shown in this figure are available as Data behind the Figure.
We first adopted a PA=160 • .4 and inc=51 • .7 for the PDS 70 disk (Keppler et al. 2019;Facchini et al. 2021) and then extracted a deprojected radius r, emission height z, surface brightness I ν , and channel velocity v for each pixel associated with the emitting surface.We then employed an iterative fitting approach, which further refines the extracted surfaces using between one to five subsequent iterations depending on the line.To ensure the robustness of the extracted points, we also filtered those pixels with unphysical z/r values based on the expected disk structure and removed those points with low surface brightnesses (2-3×rms).At each step, we visually confirmed the fidelity of the derived emission surfaces.To further reduce scatter in the extracted surfaces, we also performed each of the same two types of binning (radially binned and moving averages) as in Law et al. (2021a) with half-beam spacings.
We then fitted parametric models to all line emission surfaces using an exponentially-tapered power law, which captures the inner flared surfaces and the plateau and turnover regions at larger radii (e.g., Teague et al. 2019b).We adopted the following functional form: This is a similar form employed by Law et al. (2021a) but here we include an additional term r cavity to describe the inner gas cavity of the PDS 70 disk.
We used the Monte Carlo Markov Chain (MCMC) sampler implemented in emcee (Foreman-Mackey et al. 2013) to estimate the posterior probability distributions for: z 0 , ϕ, r taper , r cavity , and ψ.Each ensemble used 64 walkers with 1000 burn-in steps and an additional 500 steps to sample the posterior distribution function.The posteriors were approximately Gaussian with no significant degeneracies between parameters.Table 3 shows the median values of the posterior distribution, with uncertainties given as the 16th and 84th percentiles, as well as the radial range, r fit, max , considered for each surface.a Characteristic z/r is computed as the mean and 16th to 84th percentile range within 80% of the r cutoff as in Law et al. (2022Law et al. ( , 2023)).

Emission Surfaces in the PDS 70 Disk
Figure 2 shows the derived emission surfaces in the PDS 70 disk.All of the CO isotopologue lines show characteristic emission surface profiles, i.e., a sharply-rising inner component, followed by subsequent plateau and turnover phases toward larger radii (e.g., Teague et al. 2019b;Law et al. 2021a).Due to the high sensitivity and angular resolution of the observations combined with the iterative surface fitting procedure of disksurf, we are also able to confidently track emission heights even at large disk radii.
Both the J=3-2 and J=2-1 lines of 12 CO show altitudes of z/r ≈ 0.3 and reach a maximum height of approximately 40 au at a radius of 150 au.The 12 CO J=3-2 heights are generally consistent with those derived in Keppler et al. (2019), who fit a single powerlaw surface profile to the Keplerian rotation map but did not extract individual heights.The 12 CO J=3-2 emission surface shows evidence for vertical variation due to the cavity wall at ≈45 au.From a radius of ≈15-40 au, the surface is quite flat (z/r ≈ 0.1) and located near the midplane, but then abruptly increases height at a radius of ≈45 au up to z/r ≈ 0.3.We are unable to extract emission heights for the 12 CO J=2-1 line in the inner 50 au (square box in Figure 2) despite 12 CO emission being clearly detected in the inner disk (Figure 1).This is likely due to the larger beam size of the J=2-1 observations (3 times larger in beam area) rather than an intrinsic difference in emitting heights, especially since the J=2-1 and J=3-2 12 CO lines have typically been observed to originate from the same heights in other disks (Law et al. 2023).The 13 CO J=2-1 line originates from a deeper layer in the disk (z/r ≈ 0.1).
We do not resolve the vertical structure of the C 18 O emission, which suggests that C 18 O arises at, or near, the disk midplane.In addition to CO isotopologues, we also measured the HCO + J=4-3 emitting surface, which lies at z/r ≈ 0.2.Beyond ≈100 au, the HCO + surface shows an extended plateau of flat, near-midplane emission (z/r ≲ 0.1), which is not observed in any other line.At <100 au, the HCO + arises from higher altitudes than 13 CO, but at >100au, HCO + is coming from lower altitudes than 13 CO.
In addition to the cavity wall seen in 12 CO J=3-2, each of the 12 CO and 13 CO surfaces shows a shallow vertical dip between ≈90-120 au.We confirmed that each of substructures are also detected in emission surfaces extracted using the non-continuum-subtracted image cubes and with no filtering of low surface brightness points, and thus, are not artifacts of the continuum subtraction or surface extraction process.To quantify the significance of these features, including the cavity wall, we follow the fitting procedure outlined in Law et al. (2021a).In brief, we remove a local baseline around each vertical substructure and then fit for its radial position, width, and depth.Table 4 shows the properties of each substructure, which are labeled with a "Z" to indicate that they are vertical variations followed by their radial location rounded to the nearest whole number in astronomical units.The cavity wall at Z42 is the deepest vertical feature, followed by the Z95 dip in 13 CO, while the dips at Z119 and Z113 in the 12 CO surfaces are the most shallow.
In Figure 3, we overlay the inferred emission surfaces on peak line intensity maps to better illustrate their 3D geometries.We generated all peak intensity maps with the 'quadratic' method of bettermoments using the full Planck function.

Radial Disk Chemical Substructure
While several molecular lines considered here have been previously observed toward the PDS 70 disk (e.g.,    Long et al. 2018;Facchini et al. 2021), these new data represent improvements in both angular resolution and sensitivity of more than a factor of two.Thus, in addition to mapping the vertical gas structure, we can also constrain the presence of small-scale radial substructures at the highest spatial resolution to date.To do so, we computed azimuthally-averaged radial profiles using the radial profile function in the GoFish Python package (Teague 2019) to deproject the zeroth moment maps.During deprojection, we incorporated the derived emission surfaces listed in Table 3.This is particularly important for highly elevated surfaces, e.g., 12 CO, to derive accurate radial locations of line emission substructures (e.g., Rosotti et al. 2021;Law et al. 2021b).Figure 4 shows the resultant radial profiles.
The location of line emission substructures are labeled according to their radial location rounded to the nearest au following standard nomenclature (e.g., Huang et al. 2018;Law et al. 2021b).The 12 CO lines show centrallypeaked profiles with a dip at 21 au, followed by an emission peak between 40-50 au and an extended plateau out to ≈250 au.The higher resolution 12 CO J=3-2 profile resolves an inner ring at ≈5 au, indicating the possible presence of small-scale inner disk substructure.The 13 CO and C 18 O J=2-1 lines instead show a central gap and broad ring at 67 au.The HCO + profile is more similar to that of the 12 CO lines with a dip at 25 au and ring at 50 au.However, HCO + is much more sharply centrally peaked with no indications of a central dip in the inner few au, despite the HCO + having the highest spatial resolution.Overall, the radial morphologies are consistent with existing lower resolution observations (Long et al. 2018;Facchini et al. 2021), namely that 12 CO and HCO + peak interior to the mm dust ring and 13 CO and C 18 O are co-located with the mm continuum peak.As discussed in Facchini et al. (2021), these differences are likely due to varying line optical depths and the presence of high temperatures and illumination at the edge of the cavity wall.
We also computed the disk size of each emission line by determining the radius in which 90% of the total flux is contained (e.g., Tripathi et al. 2017;Ansdell et al. 2018).Table 2 lists the derived values.The 12 CO lines have have sizes of ≈200 au, while 13 CO and HCO + have sizes of ≈165 au.The C 18 O is the most compact line at ≈140 au.The mm continuum edge of the PDS 70 disk is approximately 100 au.This implies gas-to-dust sizes of ≈1.4-2, which are typical, albeit on the lower end, for large, resolved disks (e.g., Law et al. 2021b;Sanchis et al. 2021;Long et al. 2022).

Comparison with NIR Scattering Surface
The PDS 70 disk has been extensively observed at NIR wavelengths, which revealed an elliptical ring, inner disk component, and asymmetrical feature to the northwest of the star due to either a double ring or planet-disk interactions (e.g., Hashimoto et al. 2012;Keppler et al. 2018;Müller et al. 2018;Mesa et al. 2019;Juillard et al. 2022).Here, we only focus on the scattering surface in the outer disk, which is tracing the small dust grain population, and compare this to the derived line emission surfaces.Only a handful of disks have independent measurements of both NIR and gas heights, and as a result, the general relationship between the vertical distribution of gas and small dust in disks is not yet well established.Figure 5 shows the power-law NIR surface derived in Keppler et al. (2018) versus the extracted molecular gas heights.We also directly measured the height of the scattered light ring at ≈54 au in the PDS 70 disk by fitting an ellipse to the peak flux using the 2016 Keppler et al. (2018).To do so, we followed the same approach outlined in Rich et al. (2021) and described in detail in Appendix C of Law et al. (2023).We derived a height of 4.5 ± 0.1 au at a radius of 53.7 ± 0.1 au, which agrees well with the overall power-law trend from Keppler et al. (2018) and is shown as a diamond marker in Figure 5.The NIR surface is located at the same vertical height as the 13 CO surface over the entire radial range of the NIR disk.Existing measurements in other disks show that the NIR surface typically lies between the 12 CO and 13 CO layers, but with considerable diversity in the relative heights of the small dust grains (e.g., Rich et al. 2021;Law et al. 2022Law et al. , 2023;;Paneque-Carreño et al. 2023).In particular, the wellstudied IM Lup and HD 163296 disks show comparable heights in 12 CO gas and small dust grains within 100 au, but significantly more elevated 12 CO emitting heights at larger radii (Rich et al. 2021;Paneque-Carreño et al. 2023).Given that NIR emission in PDS 70 only extends to ≈115 au (Keppler et al. 2018), it is thus possible that the PDS 70 disk shows slightly lower than expected small dust heights compared to that of 12 CO.However, we caution that only a small number of disks have had their CO isotopologue and small dust heights jointly mapped and additional disk observations are required to better place PDS 70 in context.

Disk Thermal Structure
We extracted radial (Section 3.4.1)and vertical (Section 3.4.2) gas temperatures for the PDS 70 disk using 12 CO, 13 CO, and HCO + .To do so, we followed the procedures of Law et al. (2021a).We assume that all lines are optically thick and in LTE, which is expected to be the case at the typical densities and temperatures of protoplanetary disks (e.g., Weaver et al. 2018).We also assume that the line emission fills the beam due to the high angular resolution of the observations.We then repeated the surface extractions as in Section 2.3 on the non-continuum-subtracted image cubes, which ensures that we do not underestimate the line intensity along lines of sight containing strong continuum emission (e.g., Boehler et al. 2017).We then converted the peak surface brightness I ν of each extracted pixel to a brightness temperature using the full Planck function and assume the resulting brightness temperature is equal to the local gas temperature.
All subsequent radial and 2D gas temperature distributions are taken directly from individual surface measurements, rather than from mapping peak brightness temperatures back onto derived emission surfaces (i.e., Figure 3) or radially-deprojecting peak intensity maps (see Appendix B).Hence, we only consider the brightness temperature of those pixels with derived emission heights.

Radial Temperature Profiles
Figure 6 shows the radial temperature distributions along the emission surfaces.The 12 CO lines trace the warmest temperatures, followed by HCO + and then 13 CO.For all lines, the temperatures generally increase toward the central star, except in the inner ≈60-70 au where the temperatures of the CO lines plateau or decrease.This is unsurprising since the CO lines all either show a central cavity or deep gap in the inner disk, which means that at these radii, the lines are likely no longer optically thick and thus do not trace the true gas temperature.We also confirmed that the radial temperatures extracted directly from the emission surfaces in Figure 6 show excellent agreement with the peak intensity profiles generated from Figure 3, which are shown in Appendix B.
The J=2-1 and J=3-2 lines of 12 CO are generally comparable in temperature, but the J=3-2 line has a peak temperature of ≳40 K that is a few K warmer than that of the J=2-1 line, which is likely due to beam dilution at smaller radii, since the J=2-1 beam is more than twice as large than that of the J=3-2 line.While the 13 CO and HCO + temperatures profiles appear smooth, both 12 CO profiles show the presence of a dip at ≈120 au.This dip occurs exterior to the continuum edge and corresponds to a small vertical dip in molecular emission height at the same radius (Section 3.1).This temperature decrease could either be due to locally reduced gas surface density, in which the emission surface traces a deeper and thus cooler layer, or due to a change in heating or radiation properties related to the continuum edge.To better quantify these temperature dips, we use the same method as for the vertical substructures in Section 3.1.The dip in temperature for 12 CO J=2-1 occurs at 115 au, while in 12 CO J=3-2, it is at 127 au.Both features have comparable widths (≈20-30 au) and depths (≲10%).A full listing of temperature substructure properties is provided in Appendix B.
To better characterize the radial temperature profiles, we fitted each line with a power law profile as: where q is the slope and T 100 is the brightness temperature at 100 au.Table 5 shows the fitted parameters.
Overall, the profiles have slopes of q ≈ 0.4-0.8, which is generally consistent with other disks (e.g., Law et al. 2022), with 12 CO J=3-2 having the steepest radial temperature profile.Law et al.

2D Temperature Profiles
Figure 7 shows the thermal structure of the CO isotopologues and HCO + emitting layers as a function of (r, z).By combining all lines with derived surfaces, which trace different heights in the disk, we can construct a 2D model of the temperature distribution.Following Law et al. (2021a), we fit a two-layer temperature model in which disk midplane and atmosphere temperature are power laws smoothly connected with a hyperbolic tangent function (Dartois et al. 2003;Dullemond et al. 2020): T mid (r) = T mid,0 (r/100 au) where z q (r) = z 0 (r/100 au) β .The α parameter defines the height at which the transition in the tanh vertical temperature profile occurs, while β describes how the transition height varies with radius.We used the MCMC sampler in emcee (Foreman-Mackey et al. 2013) to estimate the following seven parameters: T atm,0 , q atm , T mid,0 , q mid , α, z 0 , and β.We used 256 walkers, which take 500 steps to burn in and then an additional 5000 steps to sample the posterior distribution function.Figure 8 shows the fitted 2D temperature profiles and Table 6 lists the fitted parameter values and uncertainties, which are given as the 50th, 16th, and 84th percentiles from the marginalized posterior distributions, respectively.The presence of embedded massive protoplanets may be expected to either locally or perhaps globally alter the vertical gas structure of protoplanetary disks, which in turn, would affect the properties of the material available to be accreted by those planets.The PDS 70 system provides an ideal environment to search for such effects.Here, we first discuss the potential presence of any local perturbations due to PDS 70b and 70c and then compare the overall vertical structure of the PDS 70 disk with other disks in the literature.

Local Perturbations and Cavity Wall
In the 12 CO J=3-2 emission surface, we identify relatively flat (z/r ≈ 0.1) heights that overlap with the planet locations, i.e., within the central ≈45 au.This vertically-flat region extends to the edge of the cavity wall, where the heights sharply rise again to z/r ≈ 0.3 (Figure 2).This is the same region that shows a deep drop in brightness temperature, reaching as low as ≈15 K (Figure 6).
The flat emitting heights observed within ≈45 au are consistent with expectations of a transition disk with a deep gas-and dust-depleted central cavity, which in the case of the PDS 70 disk has been carved out by PDS 70b and 70c (Bae et al. 2019;Portilla-Revelo et al. 2023).Line emission surfaces are tracing regions where optical depths reach unity.Thus, in this region near the planets, the overall line optical depth decreases due to the reduced gas surface density and the surface is pushed deeper into the disk, i.e., close to the midplane, which naturally explains both the dip in vertical heights and the lower brightness temperatures.This is the first transition disk where this effect is directly observed in line emitting heights (see Section 4.2 for more details).For all other lines besides 12 CO J=3-2, no estimates of vertical heights within the central cavity were possible due to insufficient angular resolution or intrinsically flat emitting heights likely due to the low gas surface densities.

Outer Disk Vertical Structure of PDS 70 in Context
Although we do not see any clear local perturbations around PDS 70b or 70c beyond that of the cavity wall, the overall structure of the gas disk may still be influenced by the presence of both planets in the form of, e.g., dynamical or planet-disk interactions.In particular, here we focus on the outer disk vertical gas structure.One way to assess this is to compare the PDS 70 disk with others in the literature with similar constraints on their line emitting heights.
To do this, we computed the characteristic z/r heights of each line in the PDS 70 disk following the procedure outlined in Law et al. (2022), i.e., mean z/r within the steeply-rising portion of the surfaces, for consistent comparison with literature values.Table 3 lists the computed values.We first consider the CO isotopologues and then HCO + : CO Isotopologues: Figure 9 shows the characteristic z/r emission heights of 12 CO and 13 CO versus the stellar mass and disk size compiled from consistently derived heights (Law et al. 2023).We also include several expected scaling relations, i.e., assuming line emis-   6.The 12 CO and 13 CO J=2-1 (left) and 12 CO J=3-2 and HCO + J=4-3 (right) are shown in separate panels for visual clarity but the 2D temperature was fit to both sets of lines simultaneously.The same color scale is used for the data and fitted model in each panel.Points excluded from the fits are show as hollow markers and contours show constant temperatures.The uncertainty of the temperature measurements, which is not shown here, can be found in Figure 7.   2019).Several scaling relations, i.e., line emission heights scale with gas pressure scale height, and a previously-identified positive 12 CO height trend and disk size are labeled (see Law et al. 2022Law et al. , 2023, for more details).
sion heights scale with gas pressure scale heights, and a previously-identified positive 12 CO-R12 CO trend (for more details, see Law et al. 2022Law et al. , 2023)).With respect to the literature sample, the PDS 70 disk appears quite typical in terms of its 12 CO and 13 CO emitting heights compared to disks with similar host stellar masses and disk sizes.This is true for all disks, or if we only consider the sub-sample of transition disks.Thus, the two planets in the PDS 70 disk do not seem to be substantially altering the heights of the line emission surfaces via, e.g., planet-disk interactions.However, we note the caveat that many of the literature disks used as the comparison sample also have indirect evidence for the presence of embedded planets.Thus, the comparison of emitting heights in the PDS 70 disk versus those disks without planets is not necessarily as straightforward as Figure 9 suggests, but such a comparison nonetheless establishes that the PDS 70 disk has a vertical gas structure that is typical of the large, resolved disks thus far studied in detail with ALMA.The presence of embedded planets could, for instance, have resulted in flatter emitting surfaces due to the removal or redistribution of gas via vertical flows (e.g., Teague et al. 2019a;Yu et al. 2021;Galloway-Sprietsma et al. 2023), or alternatively, driven dynamical interactions which would inflate the disk gas vertical distribution (e.g., Montesinos et al. 2021;Kuo et al. 2022).However, the close-in radial locations of the PDS 70 planets may provide an explanation, in which the vertical disk gas distribution is decoupled from planet formation occurring in the inner few 10s of au.Several previous works (Paneque-Carreño et al. 2023;Izquierdo et al. 2023) have identified vertical variations in line emitting heights that are coherent across several tracers at the radial locations of suspected planetary companions or kinematic deviations.In the PDS 70 disk, we find vertical variations at radii between ≈95-120 au in both 12 CO and 13 CO.However, no evidence for kinematic deviations or additional embedded planets exists in the outer disk of PDS 70.Thus, it is more likely that these vertical substructures are instead related to the edge of the mm continuum disk due to, e.g., changing radiation fields at the dust edge.
HCO + : Only a few measurements of HCO + emitting heights exist in protoplanetary disks.Paneque-Carreño et al. (2023) derived heights of the J=1-0 line in the MAPS disks ( Öberg et al. 2021) and for the J=4-3 line in the WaOph 6 disk, while Huang et al. (2020) constrained the J=3-2 emitting surface in the GM Aur disk.
The J=1-0 lines showed typical heights of z/r ≲ 0.1, but cannot be directly compared to PDS 70, because the difference in excitation properties between the J=4-3 (E u ≈ 43 K) and J=1-0 (E u ≈ 4 K) lines mean that we do not necessarily expect both lines to trace the same disk vertical layers.The WaOph 6 measurement, although using the same transition as in PDS 70, suffered from high uncertainties with a possible z/r ranging from 0.1 to 0.5 due to the coarse beam size (≈0.′′ 3) and lower SNR of the observations.This range of heights is consistent with what we measure in PDS 70 but does not permit a detailed comparison of the emitting surfaces.The most similar previous measurement is in the transition disk GM Aur, where Huang et al. (2020) used high angular resolution observations to fit a single power law profile to the HCO + J=3-2 emission surface and found an approximate height of z/r ≳ 0.1.Although not identical, the J=3-2 (E u ≈26 K) line has more comparable excitation properties to that of the J=4-3 line.More-Law et al.
over, both PDS 70 and GM Aur are T Tauri stars of similar spectral types and host transition disks with central cavities of comparable size with the main difference being that GM Aur is at least twice as large in overall disk size in both its mm dust and 12 CO line emission.For consistency, we re-fit the HCO + J=3-2 data in the GM Aur disk from Huang et al. (2020) using an exponentially-tapered power-law with the same procedure as in Section 2.3.
HCO + is emitting from similar disk vertical regions (z/r ≈ 0.2) in both the GM Aur and PDS 70 disks, and thus, similar to the CO isotopologues, PDS 70 does not appear to show any significant differences, which could be ascribed to the influence of PDS 70b or 70c.As noted in Section 3.1 and shown in Figure 10, the relative heights of the 13 CO J=2-1 and HCO + J=4-3 surfaces cross over beyond ∼100 au.We find the same trend in the GM Aur disk, with the HCO + surface located at higher altitudes than 13 CO in the inner ∼200 au, while at large radii, the 13 CO is coming from a higher disk layer than HCO + .The cross-over point of these surfaces occurs directly exterior to bulk of the mm continuum emission in both disks, which suggests that this effect may be due to a change in illumination at the continuum edge.When combined with the lower densities expected at larger disk radii, this, in turn, would likely alter how far ionizing photons can penetrate in the disk and could thus naturally explain the lower heights of the HCO + surfaces in the outer disk.Although intriguing, with only two sources, it is difficult to make any additional conclusions about the vertical distribution of HCO + in protoplanetary disks more generally.We stress that high-angular resolution HCO + observations in additional disks, especially those with different properties than PDS 70 and GM Aur, are urgently needed to set stringent constraints on disk ionization structure.

Comparison to the LkCa 15 Transition Disk
Here, we compare the PDS 70 disk to that of LkCa 15, which is one of the few transition disks whose vertical emission structure has been mapped in detail with multiple molecular lines (Leemker et al. 2022;Law et al. 2023).LkCa 15 has a similar stellar spectral type (K5) and age (∼5 Myr) (Donati et al. 2019) and while the LkCa 15 disk is much larger than that of PDS 70 in both mm dust continuum and CO isotopologue emission by a factor of ≈2-3, both disks have nearly identical central cavity sizes (e.g., Piétu et al. 2006;Facchini et al. 2020).Moreover, both disks have comparable gas temperatures within 200 au and similar 12 CO and 13 CO emission layer z/r heights (Leemker et al. 2022;Law et al. 2023).Thus, they are ideal systems to compare how line emitting heights change across their central cavities.
Figure 11 shows line emission surfaces in each disk zoomed into the location of the central cavity, as indicated in the corresponding azimuthally-averaged mm continuum radial profile.The PDS 70 disk shows an increase in emitting heights in 12 CO J=3-2 associated with the edge of the cavity wall, while no such effect is seen in LkCa 15.While in Figure 11, we show the J=2-1 lines of 12 CO and 13 CO for the LkCa 15 disk, emitting height measurements also exist for 12 CO J=3-2 and 13 CO J=6-5 (Leemker et al. 2022), but neither of these surfaces show a similar height increase as in PDS 70.
However, the absence of this cavity wall effect in the LkCa 15 disk is likely an observational limitation, rather than indicating any fundamental differences in inner disk structure.For instance, in Figure 11, the 12 CO J=2-1 beam is ≈0.′′ 3, which at the distance of LkCa 15 (157 pc;Gaia Collaboration et al. 2021), corresponds to a physical scale of ≈50 au, which is approximately the size of the central cavity.In comparison, our 12 CO J=3-2 observations of PDS 70 trace physical scales of ≈11-15 au, a factor of several times better.We note that while the 13 CO J=6-5 surface in LkCa 15 from Leemker et al. (2022) was derived from observations of comparable angular resolution, emitting heights interior to ≈40 au could not be extracted.Moreover, 13 CO originates at lower elevations than that of 12 CO, which means that this height increase at the cavity wall would be intrinsically more difficult to observe using 13 CO.
Emission from 12 CO J=2-1, 3-2 and 13 CO J=6-5 have been detected in the central cavity of LkCa 15 (Jin et al. 2019;Leemker et al. 2022;Law et al. 2023), which suggests that emission height constraints here would be possible with sufficiently sensitive and highangular-resolution data.Additionally, Leemker et al. (2022) found that the gas column density in the dust cavity (≈45 au) of LkCa 15 drops by a factor of >2 compared to the outer disk, with an additional orderof-magnitude decrease inside 10 au.In comparison, the PDS 70 cavity shows a higher gas depletion (≳10×) over its full cavity (Portilla-Revelo et al. 2023).Thus, while we still expect flatter emission surfaces in the LkCa 15 cavity, this cavity wall feature may be less pronounced in the LkCa 15 disk due to its more modest gas density decrease.Nonetheless, observations from the ongoing exoALMA Large Program will probe similar physical resolutions as in the PDS 70 disk in several transition disks, including LkCa 15, and should be able to confirm if such changes in emitting heights at cavity walls are commonplace or if the PDS 70 disk structure is unique.

CO Snowline Location
The elemental carbon-to-oxygen (C/O) ratio is of vital importance to understand the accretion history of protoplanets and connect disk chemistry with planetary organic compositions and atmospheres (e.g., Madhusudhan 2019;Öberg & Bergin 2021).While chemical evolution may alter the C/O ratio (e.g., Eistrup et al. 2018;Cridland et al. 2019;Krijt et al. 2020), the condensation of major volatile species across the disk has the most dramatic influence on the elemental ratio of gas-phase material (e.g., Öberg et al. 2011).These condensation fronts, often referred to as snowlines, occur at specific temperatures and thus, with detailed knowledge of the disk temperature structure, it is possible to estimate the location of these snowlines directly.Here, we focus on the CO snowline, which is readily observationallyaccessible at spatial scales probed by ALMA.
Figure 12 shows the midplane gas temperature as inferred from the 2D empirical fit (Figure 8).It is important to note that the midplane temperature from this fit is an extrapolation from higher disk heights since we lack lines that directly trace the disk midplane interior to a radius of 100 au.As a result, we also show the C 18 O J=2-1 peak brightness profile temperature in Figure 12, which we expect to be emitting at, or near, the disk midplane (Section 3.1).Beyond the central cavity, the C 18 O brightness temperature is generally consistent with our estimated midplane temperature and is ≈5-10 K colder at larger radii (≳100 au).As C 18 O is not expected to be fully optically thick, especially at large radii, it only provides a lower limit on the gas temperature, which is consistent with our warmer midplane temperature estimates.If we adopt a reasonable range of possible CO freezeout temperatures from 18-22 K (e.g., Collings et al. 2004;Öberg et al. 2005;Martín-Doménech et al. 2014;Facchini et al. 2017), we estimate an approximate snowline radius of ≈56-85 au.Despite this range, we can confidently conclude that both PDS 70b and 70c are located interior the CO snowline.Thus, the nearby gas from which both planets are accreting likely has a superstellar1 C/O ratio.Depending on the exact location of the CO 2 snowline, this gas likely has a C/O of ≲1.0 ( Öberg et al. 2011), with the exact value depending on the details of radial transport and chemical processing (e.g., Krijt et al. 2018Krijt et al. , 2020)).From our midplane temperature estimate, we would estimate a CO 2 snowline location of ≈21-26 au (assuming freeze-out temperatures of 42-52 K), which is consistent with C/O ≈ 1, with the caveat that these are extrapolations to radii within the inner cavity where we lack direct temperature measurements.
Overall, this is consistent with the findings of Facchini et al. (2021), who found that the PDS 70 molecular layer hosts a high C/O ratio inferred from bright C 2 H line fluxes and their lower limit on the CS/SO column density ratio (e.g., Bergin et al. 2016;Semenov et al. 2018;Miotello et al. 2019;Le Gal et al. 2021).The fact that both planets in the PDS 70 disk are likely accreting gas with a high C/O ratio provides important context for interpreting future atmospheric planetary measurements from, e.g., JWST, and for properly modeling their formation histories within the PDS 70 disk (e.g., Cridland et al. 2023).
Further observations of temperature tracers, such as additional J lines of CO isotopologues (e.g., Leemker et al. 2022), especially those closer to the midplane and in the inner disk would allow for a better constrained temperature model and lead to a more accurate prediction of the snowline location.Complementary to this, observations of N 2 H + , which has been shown to be an accurate observational tracer of the CO snowline (Qi et al. 2019), would provide an independent estimate of the CO snowline location.If combined with our temperature-based estimate, this would provide a unique opportunity to measure the CO binding energy in an astrophysical context and provide important and independent context to laboratory-based estimates.

CONCLUSIONS
We present observations of a suite of CO isotopologue and HCO + lines toward the PDS 70 disk at high angular resolution (≈0.′′ 1).We extracted line emission surfaces and conclude the following: 1. 12 CO emission surfaces originate in elevated disk regions (z/r ≈ 0.3), while the less abundant 13 CO emitting heights trace deeper disk layers (z/r ≈ 0.1).We also derived emitting heights for the HCO + J=4-3 line, which arises from z/r ≈ 0.2 in the inner 100 au and shows an extended, flat (z/r ≲ 0.1) component at larger radii.
2. In the 12 CO J=3-2 line, we clearly resolve the vertical dip and steep rise in emitting heights at ≈15-45 au due to the cavity wall in the PDS 70 disk.This is the first transition disk where this effect has been directly seen in emitting heights.
3. The emitting heights of the CO isotopologues in the outer disk of PDS 70 appear typical for its stellar mass and disk size and do not appear to be substantially altered by the presence of its two embedded planets.Overall, this suggests that the outer disk is largely decoupled from the planet formation occurring in the inner few 10s of au.
4. By combining CO isotopologue and HCO + lines, we derive an empirical 2D temperature structure for the PDS 70 disk.Using this, we estimate an We extracted radial profiles from the peak intensity maps (Figure 3) by adopting the inferred emission surfaces and following the same procedure as in Section 3.2.In Figure 14, we compare them to the radial temperatures inferred directly from the emission surfaces.Overall, we find excellent agreement between the two curves.The most notable differences occur in the inner 100 au of the 13 CO J=2-1 profile, where the peak intensities have temperatures which are at most ≲5 K larger than the surface-extracted ones.We also note that the dips identified at ≈120 au in both 12 CO J=2-1 and J=3-2 from Figure 6 are also present in the peak intensities.
We also find a small dip in the peak intensity profiles of 12 CO J=3-2 at ∼50 au that is not apparent in the radial brightness temperature profiles.This dip was identified previously in integrated intensity by Keppler et al. (2019) and was attributed to continuum absorption of the back side of the disk, since it appears coincident with the mm continuum gap.However, given its appearance in the peak intensity profile, which should be less affected by continuum absorption, this suggests that it is a real feature and may hint at a true gas and dust gap.While we do not see this feature in any of the other CO lines, this is perhaps not surprising given the small width of this feature and the larger beam sizes of the other CO lines.Additional observations at high angular resolution of other CO isotopologues are necessary to confirm the nature of this feature.
As noted in Section 3.4.1,we also determined the properties of the temperature dips seen in Figures 6 and 14 using the same method as for the vertical substructures, i.e., via removing a local baseline (Law et al. 2021a).Table 7 shows the properties of each of these features, including the 50 au dip in 12 CO J=3-2 identified in the peak intensity radial profiles.

Figure 1 .
Figure 1.CO isotopologue and HCO + zeroth moment maps and continuum images (from top left to bottom right) of the PDS 70 disk.Panels for each disk have the same field of view.The locations of the planets PDS 70b and 70c are marked in each panel (Wang et al. 2021).Color stretches were individually optimized and applied to each panel to increase the visibility of outer disk structure.The synthesized beam and a scale bar indicating 50 au is shown in the lower left and right corner, respectively, of each panel.The 344 GHz continuum is from Isella et al. (2019).Table 2 lists details about each image.

Figure 2 .
Figure2.CO isotopologue and HCO + emission surfaces for the PDS 70 disk.Large gray points show radially-binned surfaces and small, light gray points represent individual measurements.The orange lines show the exponentially-tapered power law fits from Table3.The solid lines show the radial range used in the fitting, while the dashed lines are extrapolations.Lines of constant z/r from 0.1 to 0.5 are shown in gray.Vertical substructures and the locations of PDS 70b and 70c are marked in each panel(Wang et al. 2021).The lack of 12 CO J=2-1 heights in the inner 50 au (square box) is due to insufficient angular resolution rather than the absence of emission at these radii.The FWHM of the major axis of the synthesized beam is shown in the bottom right corner of each panel.The emission surfaces shown in this figure are available as Data behind the Figure.

Figure 3 .
Figure 3. Peak intensity maps with overlaid contours showing the fitted emission surfaces, as listed in Table3.No surface could be derived for the C 18 O J=2-1 line, which is consistent with emission originating from near the midplane.The synthesized beam and a scale bar indicating 50 au is shown in the lower left and right corner, respectively, of each panel.

Figure 4 .
Figure 4. Deprojected radial integrated intensity profiles lines and the 344 GHz continuum in the PDS 70 disk.Shaded regions show the 1σ uncertainty.Solid lines mark emission rings and dotted lines indicate gaps.The FWHM of the synthesized beam is shown by a horizontal bar in the upper right corner of each panel.

Figure 5 .
Figure 5. Emission surfaces of 12 CO and 13 CO J=2-1 (top) and 12 CO J=3-2 and HCO + J=4-3 (bottom) versus the NIR scattering surface from Keppler et al. (2018) in the PDS 70 disk.The red marker shows the individual height measurement of the PDS 70 NIR ring derived in this work.The lines show the moving average surfaces and gray shaded regions show the 1σ uncertainty.The FWHM of the major axis of the synthesized beam is shown in the bottom right corner of each panel.

Figure 6 .
Figure 6.CO, 13 CO, and HCO + radial brightness temperature profiles of the PDS 70 disk.These profiles represent the mean temperatures computed by radially binning the individual measurements, similar to the procedure used to compute the radiallybinned surfaces (see Section 3.4).Vertical lines show the 1σ uncertainty, given as the standard deviation of the individual measurements in each bin.The solid gray lines show the power fits from Table 5, while the dashed lines are extrapolations.The FWHM of the major axis of the synthesized beam is shown in the bottom right corner of each panel.

Figure 7 .
Figure 7. 2D temperature profiles of the 12 CO and 13 CO J=2-1 (left) and 12 CO J=3-2 and HCO + J=4-3 (right) lines in the PDS 70 disk.Points are those from the binned surfaces and error bars are the 1σ uncertainties in z.Temperature measurements with radii near the emission gap (Figure 6) are marked by hollow markers and not used in the 2D temperature fits.The 2D temperature profiles shown in this figure are available as Data behind the Figure.
Influence of Embedded Planets PDS 70b and 70c on Disk Vertical Structure

Figure 8 .
Figure 8.Comparison of the measured temperatures (points) with the fitted 2D temperature structures (background) for the PDS 70 disk, as listed in Table6.The 12 CO and 13 CO J=2-1 (left) and 12 CO J=3-2 and HCO + J=4-3 (right) are shown in separate panels for visual clarity but the 2D temperature was fit to both sets of lines simultaneously.The same color scale is used for the data and fitted model in each panel.Points excluded from the fits are show as hollow markers and contours show constant temperatures.The uncertainty of the temperature measurements, which is not shown here, can be found in Figure7.

Figure 9 .
Figure9.Comparison of the characteristic 12 CO (top row) and 13 CO (bottom row) emission heights of the PDS 70 disk (in pink) vs other disks in the literature(Law et al. 2023, and see references therein).Transition disks are shown as symbols with a hollow center.All stellar masses are dynamical, z/r values are extracted homogeneously via directly-extracted emission surfaces, and disk sizes represent the radius enclosing 90% of the total flux.The dynamical stellar mass of PDS 70 is taken fromKeppler et al. (2019).Several scaling relations, i.e., line emission heights scale with gas pressure scale height, and a previously-identified positive 12 CO height trend and disk size are labeled (seeLaw et al. 2022Law et al. , 2023, for more details).

Figure 10 .
Figure 10.Top: Comparison of HCO + surfaces in the PDS 70 (left) and GM Aur (right) disks.The GM Aur data were obtained from Huang et al. (2020) and we extracted the emission surface as in Section 2.3.The single power law fit from Huang et al. (2020) is also shown for reference.Otherwise, as in Figure2.Bottom: Zoom-in on the HCO + (orange) and 13 CO (gray) emission surfaces.The mm continuum disk size is marked by a vertical line.The 13 CO surface in the GM Aur disk is fromLaw et al. (2021a).

Figure 11 .
Figure 11.Zoomed-in emission surfaces of transition disks around PDS 70 and LkCa 15 (top row ) compared to azimuthallyaveraged millimeter continuum radial profiles (bottom row ).Lines of constant z/r from 0.1 to 0.5 are shown in gray.Emission surfaces and continuum profiles in LkCa 15 are from Leemker et al. (2022); Law et al. (2023) and Facchini et al. (2020), respectively.Only the PDS 70 disk shows an increase in line emitting heights at the location of the cavity wall.

Figure 12 .
Figure 12.Midplane gas temperature as a function of radius in the PDS 70 disk as inferred from the 2D empirical temperature fit.The blue horizontal region indicates a range of CO freezeout temperatures (18-22 K), which corresponds to a midplane snowline radius of ≈56-85 au.The peak brightness temperature profile of C 18 O J=2-1 extracted in a narrow ±5 • wedge along the disk major axis is shown in orange.

Figure 13 .Figure 14 .
Figure 13.Channel maps of the 12 CO J=2-1 emission of the PDS 70 disk.The synthesized beam is shown in the lower left corner of each panel and the LSRK velocity in km s −1 is printed in the upper right.

Table 1 .
Details of ALMA LB Observations of PDS 70

Table 2 .
Image Cube Properties

Table 4 .
Properties of Vertical Substructures

Table 5 .
Radial Temperature Profile Fits Local temperature dips (D) labeled according to their approximate radial location in au.See Appendix B for further details.

Table 6 .
Summary of 2D Temperature Structure Fits