© The Authors 2023

1 Introduction

The link between the gas and ice chemistry in protoplanetary disks and the resulting chemical structure in the atmospheres of (giant) exoplanets has become an important tool in probing the underlying physics of planet formation. The argument follows that chemical gradients in the disk gas and ice – particularly in the carbon and oxygen tracers – are encoded into the planet’s atmosphere as it accretes its atmosphere-forming gas. As a result, if accurate measurements of the atmospheric number ratio of carbon and oxygen atoms can be made, then we can infer the location in the disk from which the planet has accreted its material. Restricting the location from which known planets form helps to constrain planet formation models both in terms of their rates – growth versus migration – and their initial conditions.

Young planets, with ages younger than ~10 Myr, offer a unique view of the connection between planet formation and the chemistry of planetary atmospheres. First, they are still self-luminous because of their latent heat of formation, which makes direct spectroscopy (and imaging) more feasible than for billion-year-old planets. Second, and in particularly for the ideal case of PDS 70, their natal protoplanetary disk remains sufficiently massive that it can be detected and chemically characterised. In this way, we can directly measure chemical properties of the gas that is feeding the planetary atmospheres and use this knowledge to infer the chemical properties of the gas that had fed the planet during its growth phase. Unfortunately, the presence of the disk makes the first point more complex because it contaminates attempts at direct spectroscopy (discussed in more detail below). We focus more on the second point in this work. Finally, these young, warm planets did not exist sufficiently long for internal processes such as mixing and sedimentation to drastically alter the observable chemical properties of its atmospheres. We discuss this point in more detail in Sect. 5 below.

The carbon-to-oxygen number ratio, or more commonly, the carbon-to-oxygen ratio (C/O), is the total number of carbon atoms divided by the total number of oxygen atoms in either the refractory (rock) and/or volatile (gas and ice) components. It has become the most popular elemental ratio for tracing the formation of planets since its introduction by Öberg et al. (2011). The reasons for its popularity are threefold: first, carbon and oxygen are the most abundant elements heavier than hydrogen and helium in the Universe, and thus their molecular forms H₂O, CO, and CO₂ tend to be the most abundant molecules after H₂ in astrophysical settings. Second, these same molecular species have strong molecular features in the near- and mid-infrared that can be observed in spectroscopic studies through either direct observations or during a stellar occultation (a transit; see the review by Madhusudhan 2019). Finally, as a result of their differences in condensation temperature, there is the aforementioned C/O gradient as a function of radius through the protoplanetary disk. Different formation models that predict different starting locations for a giant planet should predict (observable) chemical differences in the atmospheres of gas giant planets.

Measuring the planet’s atmospheric C/O offers a link to the local C/O in the protoplanetary disk where the giant planet drew down the bulk of its gas. Because the later stages of gas accretion are typically very rapid, with timescales that are much shorter than the migration timescales (Pollack 1984; Pollack et al. 1996; Cridland et al. 2016; Emsenhuber et al. 2021), the atmospheric C/O is determined by a narrow range of radii in the protoplanetary disk. To zeroth order, the local C/O in the gas is determined by the freeze-out of the most abundant volatiles at their ice lines. The most famous ice lines are the water-ice line near a gas temperature of ~150 K and the CO ice line near a temperature of ~20 K. The former is famous for its importance in setting the water abundance of planetesimals in the Solar System (Ciesla & Cuzzi 2006, and references therein), while the latter is relevant because of its indirect detection in protoplanetary disks (e.g. in TW Hya; Qi et al. 2013b). However, other ice lines, such as the CO₂ ice line, will also contribute to the overall evolution of carbon- and oxygen-bearing species in planet-forming regions (see for ex. Eistrup et al. 2018).

The local C/O ultimately depends on chemical processes impacting the arrangement of carbon and oxygen atoms between the gas and ice phases as well as the total number of carbon and oxygen atoms in the system. The latter can be represented by a global C/O and be thought of as the C/O of the gas and ice of the primordial molecular cloud of the stellar system. Because the disk and star accreted from the same molecular cloud material, it is often assumed that the stellar C/O represents this global C/O. Having a measure of the global C/O via the carbon and oxygen abundances in the host star is crucial for understanding a planet’s atmospheric C/O in the context of its planetary system and formation processes.

The paper is organised as follows: in Sect. 2 we outline the important features of the PDS 70 system, and in Sect. 3 we outline our numerical and analytic methods. The results are summarised in Sect. 4 and are discussed, along with the model and assumption caveats, in Sect. 5. We conclude our results and offer direction for future study in Sect. 6.

2 PDS 70 system

2.1 Star

The system PDS 70 was first identified as a T Tauri star in the Pico dos Dias Survey (PDS) by Gregorio-Hetem et al. (1992). It was included in the fourth (Kordopatis et al. 2013) and sixth (Steinmetz et al. 2020) date releases of the Radial Velocity Experiment (RAVE), where it was measured at a medium resolution (R ~ 7500) to spectroscopically derive its stellar atmospheric properties such as effective temperature, gravitational acceleration, and metallicity. It has an effective temperature of 4237 ± 134 K, a log ɡ = 4.82 ± 0.2, and [Fe/H] = −0.11 ± 0.1. It is at a distance of 113 pc (Gaia Collaboration 2020).

For the purpose of this paper, it would be useful to have an estimate of the photospheric abundance of carbon and oxygen of the star. Barring this measurement, we can estimate the stellar C/O based on its inferred metallicity and the relation between these values in stars found by Suárez-Andrés et al. (2018). In their sample, they found a linear relation between the star’s metallicity and its C/O, in particular, between stars that host giant planets (as PDS 70 does). Using their fitted slope (0.41 ± 0.02) and a normalisation of C/O = 0.48 ± 0.15 at [Fe/H] = 0, which is slightly lower than the commonly accepted solar value of 0.54, we find a stellar C/O = 0.44 ± 0.19 for PDS 70.

2.2 Disk

The infrared excess identified in the IRAS catalogue by Gregorio-Hetem et al. (1992) classified PDS 70 as a T Tauri system. Based on our current knowledge of T Tauri systems, this implies the presence of a dusty disk. The excess was confirmed by Metchev et al. (2004), and the disk was soon after detected in coronagraphic images of its scattered light by Riaud et al. (2006). The inner gap and transition-disk nature of the PDS 70 disk was confirmed through H-band polarmetric and L′-band imaging by Hashimoto et al. (2012), and follow-up observations at 1.3 mm with the Sub-Millimeter Array (SMA) additionally found evidence for a compact disk inside the dust gap (Hashimoto et al. 2015). The disk was studied at a higher spatial resolution with the Atacama Large (sub)Millimetre Array (ALMA) at 0.87 mm continuum and CO J = 3−2 and HCO⁺ J = 4−3 line emission by Long et al. (2018). Similar to the dust, the gas emission shows a gap near the location of the two planets (Keppler et al. 2019, and see below), it also extends much farther (~200 AU) from the host star than the dust (~80 AU).

High-resolution observations of the system with VLT/SPHERE¹ and VLT/NaCo have mapped the small dust distribution around the host star in more detail and have confirmed the presence of the inner dust disk (Keppler et al. 2018). The size of the dust disk was constrained to less than 17 AU based on the polarised intensity of the dust emission, but could not be further constrained because of a degeneracy between the outer radius and dust depletion factor in their inner disk model. Along with the more detailed look at the physical properties of the PDS 70 disk, Keppler et al. (2018) also reported the first detection of an embedded planetary companion to the PDS 70 star. Now called PDS 70b, the direct detection of this embedded planet is the first of its kind. This detection was later confirmed by Haffert et al. (2019), who also detect a second embedded planet, PDS 70c, using Hα emission.

2.3 Planets

The detection of PDS 70b (Keppler et al. 2018) and PDS 70c (Haffert et al. 2019) represents a unique opportunity in the study of planet formation. At the time of writing, they are the two only confirmed embedded planets that have been directly detected² by the emission of their young atmospheres (Wang et al. 2021) through the Hα emission of currently accreting gas (Haffert et al. 2019) and through continuum emission coming from the circumplanetary disk (Isella et al. 2019; Benisty et al. 2021). Because of their coexistence with their natal protoplanetary disk and the inferred age of the system (~5 Myr), they are likely the youngest exoplanets ever detected, in a stage of evolution consistent with the final stages of planet formation. At this final state, the planets have opened a gap in the disk because its gravitational influence has exceeded the viscous and gas pressure force otherwise governing hydrodynamics. The planets have thus largely accreted the bulk of their gas, with any remaining accretion being fed by gas that approaches and then surpasses the gap edge through meridional flows (Morbidelli et al. 2014; Teague et al. 2019). We quote the planetary properties of the two embedded planets as derived by Wang et al. (2021) in Table 1, who inferred the orbital parameters based on the assumption that the system is dynamically stable over the lifetime of the system (~5 Myr).

Both planets were observed in the near-infrared at high spatial and spectral resolution using the GRAVITY instrument on the Very Large Telescope Interferometer (VLTI) by Wang et al. (2021). Their spectral energy distribution (SED) fitting included the Exo-REM models of Charnay et al. (2018), which allows an inference of the atmospheric C/O of both planets. Unfortunately, their adequate models, which have Bayes factors within a factor of 100 of the best fit, were unable to provide stringent restrictions on the atmospheric C/O other than a constraint of C/O> 0.4 in PDS 70b. For PDS 70c, none of the Exo-REM models were considered adequate, but the best-fit Exo-REM model, whose Bayes factor is 114 times smaller than that of the best-fit model, constrained the atmospheric C/O < 0.7. This constraint, however, lies very close to the upper limit in their model’s C/O prior and thus should be taken with a grain of salt.

Cugno et al. (2021) similarly struggled to place constraints on the chemical properties of the planets in the PDS 70 system. They used medium-resolution data from VLT/SINFONI and the molecular mapping method (Hoeijmakers et al. 2018) to try and find molecular emission from the planets. They predicted relatively high upper limits on the molecular abundances of CO and H₂O (10^−4.1 and 10^−4.0 relative to hydrogen, respectively), which disagrees with the results of Wang et al. (2021). A possible reason for deviations between the two observations might be extinction either by clouds or surrounding dust that impacts the overall efficiency of the molecular mapping technique.

2.4 Summary

For a chemical characterisation, there are sufficient data to model the global chemical structure in the PDS 70 system. The stellar C/O = 0.44 and the detection of carbon-rich molecular species such as C₂H implies that the disk has likely undergone some form of chemical evolution to remove gaseous oxygen from the emitting layer of the disk. This chemical processing of oxygen, mainly in CO and H₂O, out of the gas phase is known to result in unexpectedly low CO emission in protoplanetary disks (Bruderer et al. 2012; Favre et al. 2013; Du et al. 2015; Kama et al. 2016; Ansdell et al. 2016; Bosman et al. 2018b, 2021; Schwarz et al. 2018; Krijt et al. 2020; Miotello et al. 2022). We thus interpret the chemical difference between the stellar C/O and the expected high C/O of the disk (as discussed by Facchini et al. 2021) as the chemical processing of gaseous CO in the ice phase and the locking up of carbon and oxygen rich ices into pebbles that settle to the disk midplane.

These ices are thus invisible to detection via line emission observations but might be available for the accretion into the atmosphere of the young planets in principle. Because CO is also an abundant carrier of carbon, we would expect a similar depletion of carbon as of oxygen. As a result, many works that explored disk chemistry in situations in which C/O > 1 find that a depletion of both carbon and oxygen is needed to explain the observed fluxes (Bergin et al. 2016; Miotello et al. 2019; Öberg et al. 2021; Bosman et al. 2021). We find here that such a carbon depletion is not consistent with the disk line emission observations of Facchini et al. (2021) and the models presented here, and we discuss this below.

In order to constrain the chemical properties of the disk, we used the line emission survey of the PDS 70 system of Facchini et al. (2021), which included the three most abundant CO iso-topologues and many bright hydrocarbon lines. The detection of bright lines from C₂H J = 7/2−5/2, c-C₃H₂ J = 3₂₁−2₁₂, and H¹³CN J = 3−2 is suggestive of a high carbon abundance, relative to oxygen, in the outer disk. In this work, we model the volatile chemistry and line emission for a subset of the detected species, ¹²CO, C¹⁸O, and C₂H, which broadly constrain the temperature, density, and chemical properties of the PDS 70 disk, respectively.

Our strategy proceeds as follows: We model the line emission observations of Facchini et al. (2021) to estimate the current physical and chemical structure of the PDS 70 disk. We stay close to the derived physical structure of Keppler et al. (2019), who studied the dust continuum and CO line emission of the disk. We use the derived physical structure to estimate the physical structure of the disk back when the planets presumably began forming and compute the chemical structure of this younger disk. Finally, we compute the formation of the two planets over the lifetime of the PDS 70 system (~5 Myr) to estimate their current-day atmospheric C/O. We expand on these steps in more detail below.

3 Methods

This work combines a number of numerical and theoretical work to understand and model the physical and chemical properties of the PDS 70 system. The chemical modelling was perfromed using the code Dust and Lines (DALI; Bruderer et al. 2012; Bruderer 2013), which computes the chemical and thermal evolution of the protoplanetary disk gas and ice self-consistently. The DALI code also contains a radiative transfer module for computing the radiative heating and cooling (for the thermal evolution), but is also used to produce synthetic observations of the disk models.

3.1 Current disk model

Our first step was to determine the current physical and chemical structure of the PDS 70 disk at its current age of 5 Myr based on the line observations of Facchini et al. (2021). The gas-disk model is based on a simple parametric description of a transition disk proposed by Andrews et al. (2011) and is discussed in detail in Bruderer (2013). The model follows the standard self-similar solution of a viscously evolving protoplanetary disk assuming that the gas viscosity is constant with time and radially varies as a power law (Lynden-Bell & Pringle 1974). Solutions of this form follow $${\text{Σ}}_{\text{gas}}\left(R\right)={\text{Σ}}_{\text{c}}{\left(\frac{R}{R_{\text{c}}}\right)}^{-\gamma }\exp \left\{-{\left(\frac{R}{R_{\text{c}}}\right)}^{2-\gamma }\right\},$$(1)

where Σ_c is the critical surface density at a critical radius of R_c. The capital R denotes the midplane radius. We allowed the gas and dust density to extend out to a maximum radius of R_disk.

Transition disks are characterised by the presence of a large dust gap caused by the presence of giant planets (at least in the case of PDS 70). This gap can extend all the way to the host star (in which case it is colloquially called a cavity), but transition disks can also contain a smaller inner disk of gas and dust. Thus, rather than a radially constant dust-to-gas ratio Δ_dtg, the model proposed by Andrews et al. (2011) and used in DALI reduces Δ_dtg to arbitrarily low values between the outer radius of the inner disk (R_gap) and the inner radius of the outer disk (R_cav).

Inward of R_cav, the gas and dust can also be further depleted relative to their expected density from Eq. (1). For the gas, this is justified by the assumption that gas accretion is slightly less efficient inside the dust cavity due to the gravitational influence of the embedded planets, hence ${{\sum }^{\prime }}_{\text{gas}}\left(R<R_{\text{cav}}\right)={\delta }_{\text{gas}}{\sum }_{\text{gas}}\left(R\right)$. The dust density has the form $${\text{Σ}}_{\text{dust}}=\{\begin{array}{ll}{\delta }_{\text{dust}}{\text{Δ}}_{\text{dtg}}{\text{Σ}}_{\text{gas}}\left(R\right)\hfill & R<R_{\text{gap}}\hfill \\ {10}^{-15}\hfill & R_{\text{gap}}<R<R_{\text{cav}}\hfill \\ {\text{Δ}}_{\text{dtg}}{\text{Σ}}_{\text{gas}}\left(R\right)\hfill & R>R_{\text{cav}}\hfill \end{array}.$$(2)

The final two important radii for the model are the sublimation radius $R_{\text{subl}}\sim 0.07\sqrt{L_{*}/L_{\odot }}\text{AU}$ inward of which the disk is sufficiently warm for the dust to evaporate (Bruderer 2013), and the outer radius R_out, which describes the largest radius in the simulation. For the stellar spectrum, we used a black body with the effective temperature noted above. The star had a bolometric luminosity of 0.35 L_⊙.

The vertical distribution of the gas followed a Gaussian profile with a physical scale height following H = H_c(R/R_c)^ϕ. The so-called scale-height angle h = H/R is thus h = h_c(R/R_c)^ϕ−1, where the subscript ‘c’ continues to denote the value of any variable at the critical radius R_c. The vertical distribution of the dust is more complex than that of the gas and is split into two populations. The small population represents grains with sizes ranging between 0.005 and 1 µm, and the large population represents grains 1–1000 µm. While the small grains tend to be well coupled to the gas and thus follow the vertical distribution more closely, the larger grains will tend to settle to the midplane. The vertical volume density of the dust in cylindrical coordinates follows $$\begin{array}{l}{\rho }_{\text{dust},\text{small}}=\frac{\left(1-f\right){\text{Σ}}_{\text{dust}}}{\sqrt{2\pi }Rh}\exp \left[-\frac{1}{2}{\left(\frac{0.5\pi -\theta }{h}\right)}^2\right]\text{\hspace{0.17em}and}\hfill \\ {\rho }_{\text{dust},\text{large}}=\frac{f{\text{Σ}}_{\text{dust}}}{\sqrt{2\pi }R\chi h}\exp \left[-\frac{1}{2}{\left(\frac{0.5\pi -\theta }{\chi h}\right)}^2\right],\hfill \end{array}$$(3)

where θ = arctan (R/z), χ < 1 parametrises the amount of settling that occurs for the large dust grains, and f sets their mass fraction. DALI runs in two modes. The first run takes Eqs. (1)–(3) with a list of prescribed parameters to produce the disk gas and dust distributions, and the second run allows users to prescribe their own choice of dust and gas radial and vertical distribution (discussed in more detail below). We show the adopted values for the discussed parameters in Table 2. When a single value is shown in the second column of the table, it coincides with the derived disk model of Keppler et al. (2018) based on the dust distribution.

The quantity of gas within the orbital range of the planets and the dust cavity is still an open question. Numerical simulations (e.g. Lubow & D’Angelo 2006) have shown that gas accretion can occur across the planet-induced gap, but results in a lower accretion rate onto the host star than the accretion rate outward of the gap. Theoretical studies have shown (for example, see Crida et al. 2006) that a giant planet modifies the local gas density surrounding it as its gravitational torque competes with the viscous torques of the surrounding protoplanetary disk. It is thus expected that the gas surface density should be lower around the embedded planets of the PDS 70 system.

Keppler et al. (2019) used ¹²CO J = 3−2 line emission data to constrain the gas properties within the dust cavity. They found evidence of a gas gap around the PDS 70b planet consistent with a mass of 5 M_Jup. No gas gap was inferred for PDS 70c, but the numerical simulations of Bae et al. (2019) showed that the two planets would share the gas gap. In this work, we can include the impact on the line emission of CO and C₂H from the presence of a planet-driven gap. In order to add the impact of the gas gap, we further depleted the gas density within the dust cavity by a factor of δ_planet (see Fig. 1).

In Fig. 1, we show an example of how the different forms of depletion impact the structure of the midplane volume (number) densities. while the solid blue and green lines represent the fiducial distribution of the gas and dust (respectively) in a transition disk model, the dashed lines show the effect of reducing the inner disk densities by δ_gas = δ_dust = 0.01. The orange line shows the impact of imposing a reduction in the gas density by a factor of δ_planet = 10^−2.5 within the dust cavity, which approximately equates to the 5 M_Jup model of Keppler et al. (2019).

In Fig. 2, we show an example of the vertical distribution of the gas and dust in a fiducial model at a radius of ~80 AU, outward of the cavity outer edge. The vertical distribution of the gas (blue curve) is described by a single-Gaussian profile, and the vertical profile of the gas (green line) is described by a pair of Gaussian profiles (i.e. Eq. (3)). The dashed orange lines show each of the Gaussian profiles from Eq. (3), converted into number density rather than mass density, to demonstrate the fact that the dust profile is effectively the sum of the two Gaussian profiles.

Fig. 1 Examples of possible models of the gas and dust distributions to be fed into DALI. We show the number volume density defined as ${\text{n}}_{\text{gas}}={\sum }_{\text{gas}}/\left(\sqrt{2\pi }H\mu m_{\text{H}}\right)$ for the gas. For the dust, ndust is defined to ensure the correct dust-to-gas mass ratio. The solid and dashed curves show the effect of depleting the inner disk by constant factors of δgas and δdust, and the orange curve shows the impact of adding a planet-induced gas gap with a prescribed depth of δplanet.

Fig. 2 Vertical distribution of the gas and dust at a radius of ~80 AU outward of the outer edge of the dust cavity. The dashed orange line represents the different distributions of the large and small dust (Eq. (3)). The large dust grains are distributed closer to the midplane because they tend to settle with respect to the smaller grains. As a result, the dust-to-gas ratio is slightly higher near the midplane than far above the midplane.

3.2 Disk chemistry

As previously mentioned, the chemistry was handled using the code DALI (Bruderer et al. 2012; Bruderer 2013). DALI simultaneously computes the gas and dust temperature via radiative transfer (both continuum and line) as well as the chemical evolution of the volatiles molecular species. As a result, the code computes an accurate picture of the temperature structure and molecular abundance in the upper atmosphere of the disk while simultaneously predicting the relative molecular line strengths and line emission from a wide range of species.

In particular, we focused on three molecular lines: ¹²CO J = 2−1, C¹⁸O J = 2−1, and C₂H J = 7/2−5/2. These molecular lines are sensitive (mainly) to the temperature, density structure, and chemistry of the disk, respectively. These lines, among others, were observed by Facchini et al. (2021) using the Atacama Large (sub)Millimeter Array (ALMA) at a resolution of ~0.4”, and they are shown in Fig. 3. One of the main features found in this observational work is that at this resolution, the integrated intensity is not centrally peaked. It peaks within the dust cavity instead. In addition, the strong C₂H emission indicates that a high C/O, likely above unity, is required in order for C₂H to have a significant abundance in the disk.

We focused on constraining the global C/O in the disk as a first step in understanding the chemical environment in which the planets formed. Our chemical model, discussed below, is sensitive to chemistry-driven local changes to the volatile (gas/ice) C/O such as the freeze-out, but we neglect local changes to the C/O that may be driven by the radial transport of volatiles (as seen in for ex. Booth et al. 2017; Bosman et al. 2018a; Krijt et al. 2020). Given the presence of the two planets, the strong dust trap outside of their orbit, and the evidence of tenuous gas in the inner disk (Keppler et al. 2019), an exploration of the radial dependence of C/O would be an interesting focus for future work.

We employed the chemical network used by Miotello et al. (2019) to model the global C₂H emission from an ALMA survey of protoplanetary disks. The chemical network is based on the network originally developed for DALI (Bruderer et al. 2012) and pulls reactions and rates from the UMIST Database for Astrochemistry (UDfA; Woodall et al. 2007; McElroy et al. 2013). It is optimised for the abundances of CO and related species in the gas phase and includes UV cross sections and photodissociation rates of Heays et al. (2017). Visser et al. (2018) added to the network to include nitrogen chemistry, which often acts to compete with hydrocarbons (e.g. C₂H) for available carbon.

While chemical networks that include isotope-specific reactions exist (e.g. Miotello et al. 2016), we did not include these reactions here. We instead assumed that species such as C¹⁸O scale with their more abundant isotopologues by the standard isotopic ratios in the interstellar medium (ISM). The isotopic ratios of carbon and oxygen in the ISM that are relevant to our work are [¹²C]/[¹³C] = 77 and [¹⁶O]/[¹⁸O] = 560 (Wilson & Rood 1994). By simply scaling the abundance of species such as C¹⁸O to the abundance of the standard isotope (¹²C¹⁶O), we ignored possible chemical effects such as isotopic selective photodissociation, which causes the less abundant isotopologues to be further depleted because their self-shielding is less efficient. As a result, the less abundant isotopologue becomes optically thick lower in the disk atmosphere than would be assumed from simple scaling models. Miotello et al. (2016) showed that this will slightly weaken the line intensity of these less abundant iso-topologues because the emission comes from deeper in the disk atmosphere, from a slightly cooler region of the disk. Furthermore, the excess photodissociation results in more free atomic oxygen isotopes that can be incorporated into other species.

We began the chemical simulations with the majority of chemical elements in the atomic phase (i.e. the reset scenario of Eistrup et al. 2016), except for hydrogen, which began in its molecular form, and a fraction of nitrogen that began as frozen NH₃ and gaseous molecular N₂. The initial abundances followed those selected by Miotello et al. (2019), except in the case of carbon and oxygen, which were varied. The initial abundances are summarised in Table 3. Our fiducial model used the carbon and oxygen abundances shown in the table. These values resulted in a marginally carbon-rich disk (C/O= 1.01), which is what is expected from the presence of C₂H emission. We also explored a range of C/H and O/H (and thus C/O) in order to confirm this expectation.

Fig. 3 Radial profile observed by Facchini et al. (2021) using ALMA. The grey band shows the approximate location of the (mm) dust ring according to Keppler et al. (2019). The line fluxes peak in slightly different positions, with the optically thin tracers peaking near the edge of the dust cavity, and the optically thick tracers peaking inside the cavity.

3.3 Disk model at the time of planet formation

More than likely, the physical properties of the PDS 70 disk were different when planets b and c started forming ~4−5 Myr ago. The particular starting time for planet formation is a remarkably difficult feature to constrain, particularly for exoplanet systems, for which we lack the isotopic data that are used in Solar System studies on the subject. We assumed that the cores of both planets formed within the first 1 Myr of the lifetime of the protoplanetary disk, and began our gas accretion simulation from there.

We used analytic models of protoplanetary disk evolution to determine the density structure of the gas at a disk age of 1 Myr. These models, derived by Tabone et al. (2022), assume that the disk evolves through either viscous evolution or magnetohydrodynamic (MHD) disk winds. It is currently unclear which mechanism is the main driver of angular momentum transport and thus protoplanetary disk evolution, so that modelling both separately accesses the two extremes of disk evolution. The main physical difference between the two mechanisms is the direction of their angular momentum transport. Viscous evolution transports the bulk of the angular momentum radially outward along with a small amount of the disk mass, while the bulk of the disk mass moves inward. In other words, viscous disks evolve by spreading outward to larger radii with a small fraction of their mass, while the remainder moves inward toward the host star.

Protoplanetary disks evolving due to MHD disk winds move the bulk of their angular momentum to larger heights by launching gas along the magnetic field lines emerging from the top of the disk. Unlike in viscously evolving disks, wind-driven disks do not spread their mass outward and are thus expected to remain at roughly the same size throughout their evolution (Tabone et al. 2022). For each of the above cases, Tabone et al. (2022) found that the gas surface density and disk size should evolve as $$\begin{array}{l}{\text{Σ}}_{\text{gas}}={\text{Σ}}_{\text{c}}\left(t\right){\left(\frac{R}{R_{\text{c}}\left(t\right)}\right)}^{-1+\xi }\exp \left(-\frac{R}{R_{\text{c}}\left(t\right)}\right)\text{\hspace{0.17em}where}\hfill \\ R_{\text{c}}\left(t\right)=R_{\text{c}}\left(t=0\right)\left(1+\frac{t}{\left(1+\psi \right)t_{\text{acc},0}}\right)\text{\hspace{0.17em}and}\hfill \\ {\text{Σ}}_{\text{c}}\left(t\right)={\text{Σ}}_{\text{c}}\left(t=0\right){\left(1+\frac{t}{\left(1+\psi \right)t_{\text{acc},0}}\right)}^{-\left(\frac{5}{2}+\xi +\frac{\psi }{2}\right)}.\hfill \end{array}$$(4)

The relevant timescale, t_acc,0, in these equations is the initial accretion timescale. It denotes the time required for a gas parcel to accrete from an initial position of R_c/2 to the inner region of the disk. It has the functional form of $$t_{\text{acc},0}=\frac{R_{\text{c}}\left(t=0\right)}{3h_c{c}_{s,c}\tilde{\alpha }\left(t=0\right)},$$(5)

where c_s,c is the sound speed at the (initial) critical radius, and $\tilde{\alpha }$ is the effective disk-α (Shakura & Sunyaev 1973) due to either viscous or MHD disk winds. The parameter ψ = α_DW/α_SS is the ratio of the effective α due to disk winds (DW) and viscosity (SS, the standard α from Shakura & Sunyaev 1973). In the limit of a purely viscous disk, ψ = 0, while in the opposite limit, ψ = ∞. In the limit that ψ → ∞, the time evolution of R_c and Σ_c in Eq. (4) becomes $$\begin{array}{l}{R}_{\text{c}}\left(t\right)=R_{\text{c}}\left(t=0\right)\text{\hspace{0.17em}and}\hfill \\ {\text{Σ}}_{\text{c}}\left(t\right)={\text{Σ}}_{\text{c}}\left(t=0\right)\exp \left(-\frac{t}{2t_{\text{acc},0}}\right).\hfill \end{array}$$(6)

The final parameter, ξ, denotes the mass ejection index and quantifies the effect on the mass surface density profile from the ejection of mass through the disk wind. The change in surface density profile follows from the significant fraction of the gas that can be lost to the disk wind as it extracts angular momentum. For the case of a pure wind (ψ = ∞), it has the form ξ = 1/[2(λ – 1)], where λ is the magnetic lever arm parameter and has a typical value of between 2 and 5. We show the disk parameters used in setting the disk properties at an age of 1 Myr for the PDS 70 system in Table 4.

3.4 Planet formation

We used a standard approach to the growth and migration of our planetary embryo. In the two disk models presented in Sect. 3.3, we initiated planetary embryos at a range of initial radii between 20 and 80 AU in steps of 2 AU. For completeness, we also computed the growth and evolution of planets originating from 1 to 20 AU in steps of 1 AU. The initial mass of the embryos was set by either the pebble isolation mass, using the prescription of Bitsch et al. (2018), or the total dust mass exterior to the initial orbits of the embryos, whichever was lower.

The gas accretion and planetary migration followed the work of Cridland et al. (2020a), and we summarise the important features of the model in Appendix A. For the purpose of this work, the chemical accounting of the carbon and oxygen-bearing species is important to determine the atmospheric C/O. We traced the quantity of mass flux onto the atmosphere through Eqs. (A.1), (A.2), and/or (A.5), as well as the rate that the growing planets move through their disk using Eq. (A.6). We followed the method of Cridland et al. (2020a) to compute the acquisition of carbon and oxygen into the atmospheres. We summarise the relevant information below.

3.5 Chemical acquisition

During the age of planet formation, the disk physical and chemical properties were kept constant to make the problem more easily tractable. As was done in Cridland et al. (2020a), this simplicity ignores the impact that an evolving disk can have on the overall evolution of the gas chemistry, but it allowed us to include the chemical impact of meridional flow on the chemistry of planetary atmospheres. We view this as an appropriate tradeoff for this study because ALMA line emission observations tend to probe the chemistry of the disk from which the meridional flow originates.

Before a gap in the disk was opened by the growing planet, we computed the average carbon and oxygen abundance within the planet’s feeding zone, which we equated to the proto-planetary Hill sphere. When the gap was open, we assumed that the gas flow towards the planet followed a meridional flow (Mordasini et al. 2014; Teague et al. 2019), and we assumed that the planet-feeding gas comes from between one and three gas scale heights. In this case, we averaged the carbon and oxygen abundance in the gas and ice at the edge of the gas gap over one and three gas scale heights.

We assumed that any ice frozen onto the population of small dust grains was entrained with the gas and was delivered into the planetary atmospheres. This process occurs both while the planet is embedded in the disk and after the gap is opened. The large-grain population is not as entrained as the small population, but is nevertheless expected to flow with the gas near the planet. Before the gas gap was opened, we assumed that ices frozen on the large dust grain population can accrete into the planetary atmosphere at a rate that was reduced relative to the rate of gas accretion by a factor of 1/(1 + St²), where St is the average Stokes number of the population of large dust grains.

The refractory component of the disk is known to include a significant fraction of the overall carbon (about 50%; Bergin et al. 2015). During the growth of the atmosphere, we did not include any excess carbon from the dust. The small grains dominate the dust mass above the midplane in a region that is likely to have been photochemically processed to remove this excess carbon (Anderson et al. 2017; Klarmann et al. 2018; Binkert & Birnstiel 2023). This could have contributed to the high C/O that is observed in the PDS 70 gas disk. The dust near the midplane is dominated by larger grains that were able to maintain their carbon. While we allowed for their accretion, they did not significantly contribute to the carbon abundance of the planetary atmospheres.

4 Results

4.1 Current disk model

Our first goal was to constrain the chemical and physical properties within the PDS 70 disk by comparing DALI forward models to the observed radial profile of the ¹²CO, C¹⁸O, and C₂H emission reported in Facchini et al. (2021). We began with our fiducial model using the parameters derived in Keppler et al. (2019) and neglected any deviations from the standard surface density profile due to depletion in the inner disk or the embedded planets (δ_gas = δ_dust = δ_planet = 1.0). In the following, we begin by constraining the density structure using the optically thin C¹⁸O line. With the constrained density structure, we then investigate the impact on the optically thick ¹²CO line and the chemical dependence of the optically thin C₂H line.

The synthetic images were produced using the radiative transfer package in DALI using the observational parameters of inclination = 52 degrees and distance =113 pc. The data were then convolved with the same beam as in Facchini et al. (2021), were deprojected, and were averaged over annuli with the same sizes as in Facchini et al. (2021). In addition, we used the same channel width and number of channels as in Facchini et al. (2021).

4.1.1 C¹⁸O

In Fig. 4a, we show the radial profile of C¹⁸O emission in the case of a disk that is nearly void of material inward of the outer dust ring edge, while using the remaining fiducial parameters presented in Table 2. In this case, we are exclusively sensitive to emission from the inner edge of the dust ring, and the outer disk. With fiducial parameters R_cav = 45 AU and R_disk = 120 AU, we find a good agreement between the peak location of the C¹⁸O emission in the observations (orange line) and our model (blue line). We overestimate the flux coming from the outer disk, however, suggesting that the outer disk is too extended or too diffuse. A smaller R_disk may better replicate the observed flux, but because we neglected isotopologue-selected photodissociation, the outer disk could be less abundant in C¹⁸O than is represented in our model (see for ex. Miotello et al. 2014).

In Fig. 4b, we show the C¹⁸O radial profile for the preferred model, which includes an inner disk. Adding a small inner disk into the model adds a new source of optical depth, which both cool the inner edge of the dust cavity and protect molecular species from photodissociation. The inner disk might also increase the flux of C¹⁸O at small separations, which requires that it is fairly tenuous. The presented model involves a gas depletion δ_gas = 0.01 and dust depletion δ_dust = 0.0001. The inner disk had a radius of 6 AU in this model.

As a final modification, we added an additional gas depletion within the planetary gap caused by the presence of the young pair of giant planets. In Fig. 4c, we add an additional depletion δ_planet = 10^−2.5 to the dust cavity (R_gap < R < R_cav). Reducing the gas density inside the gap leads to a slightly warmer dust wall and to a slight increase in the peak flux of C¹⁸O. In addition, the flux coming from the inner regions is lowered back to a similar level as in the model shown in Fig. 4a. Because it reproduces the peak flux as well as its location, the model in Fig. 4c is our preferred model.

Fig. 4 Radial profiles of C18O emission for three different disk models compared to the observational data of Facchini et al. (2021). The grey band shows the width of the major axis of the beam. In each case, the relevant parameters are listed in the top right corner of the figure. From left to right: in model a, the inner disk is almost completely removed, and there is effectively no gas or dust from Rcav inward. In model b, the inner disk is replaced, with the noted depletions of gas and dust. Finally, in model c, the same inner disk parameters are used, but between Rgap and Rcav, the gas is depleted by a further factor.

Fig. 5 Radial profile of the 12CO emission according to the model presented in Fig. 4c.

4.1.2 ¹²CO

In Fig. 5, we show the resulting ¹²CO emission radial profile from the disk model presented in Fig. 4c. Surprisingly, the peak emission is somewhat shifted with respect to the observed radial profile. This may be linked to the exact gas structure near the inner edge of the outer disk. In our model, we assumed that the gas density sharply cuts off at the inner edge of the outer disk (R_cav). However, its likely that the disk smoothly transitions from the density in the outer disk down to the lower density in the cavity and lower disk (van der Marel et al. 2016; Leemker et al. 2022). The additional width due to this smooth transition might easily contribute to the (optically thick) ¹²CO emission while only weakly contributing to the optically thin tracers.

The outer disk flux of ¹²CO also falls off more quickly than is seen in observations, which might be exaggerated if the outer disk is truncated, as was discussed for C¹⁸O. Our model overestimates the ¹²CO flux between R_cav and R_disk (= 120 AU or ~1”) before dropping below the observed flux. This suggests that there is some material in the PDS 70 disk that is outside R_disk, but is sufficiently tenuous that it does not contribute to the C¹⁸O flux. In addition, the slope and peak of the ¹²CO may suggest that the temperature profile of the emitting layer is not consistent between model and observations. This might be caused by our ignoring the heating contribution from the embedded planets, or because our dust distribution does not provide the proper amount of opacity. The shifted peak suggests that the tenuous gas in the gap is warmer than we computed in our model. We leave the modelling of these aspects of the data to future work. For the purpose of understanding the planet formation within the disk, the tenuous gas in the outer disk likely does not play an important role.

4.1.3 C₂H

In Fig. 6a, we show the resulting C₂H flux from the preferred model presented in Fig. 4c. Our model predicts a similar peak flux, roughly 30% higher than observed, with a similar peak location. Our model peak is shifted slightly inward, towards the inner disk, with respect to the observed peak location. This is due to bright C₂H emission that is found in the warm inner disk. This inner disk emission contributes to a bright C₂H line, and regardless of our choice of disk parameters, we find a similarly C₂H-enhanced inner disk.

Our modelling goal was to constrain the global carbon-to-oxygen in the PDS 70 disk. Our fiducial model was set such that the carbon abundance was that of the ISM, while the oxygen abundance was set lower so that C/O = 1.01. We expected a marginally carbon-rich disk based on the detection of C₂H, as well as a few other carbon-rich molecular species (Facchini et al. 2021). We tested whether this hypothesis was robust by computing the C₂H emission for a disk with an ISM-like abundance of carbon and oxygen.

In Fig. 6b, we showed the preferred model from Fig. 4c, but with an enhanced oxygen abundance such that the disk C/O = 0.4. As expected, the total C₂H flux is suppressed in this model compared to the marginally carbon-rich disk. Furthermore, the flux peak shifts to smaller radii, which suggests that the flux is dominated by the inner disk in this model. For the low C/O model, the carbon and oxygen abundances were set to their standard ISM abundances, while in the higher C/O models, the oxygen abundances were depleted to O/H = 1.0 × 10⁻⁴, and the carbon abundance was adjusted to set the C/O. Depleting the abundances relative to the ISM value is consistent with previous studies of the emission of C₂H in protoplanetary disks (Miotello et al. 2019; Bosman et al. 2021), but their the abundances are depleted further than is done here. We further discuss our choice of carbon and oxygen abundances in the discussion section below.

In Fig. 6c, we show the C₂H emission for a model with C/O = 1.1. The peak flux shifts farther out and is more consistent with the location of the observed flux peak. However, the model emission is about four times brighter than is observed, which heavily disfavours these carbon-rich models.

In all three panels of Fig. 6, we find that the inner disk is bright in C₂H. In order to assess the impact of inner disk C₂H, we artificially removed the C₂H from the inner 40 AU of each disk model and recalculated the radiative transfer of the C₂H emission. The results are shown in Fig. 7. We find that in both of the carbon-rich models, panels a and c, the removal of the inner disk C₂H emission suppressed the total C₂H flux, while in the carbon-poor model (panel b), the profile is unaffected.

Our C₂H modelling suggested that the inner and outer disk have different chemical properties, with the former showing evidence of having a lower elemental abundance than the outer disk due to its dominant C₂H flux in the preferred model. Future studies of PDS 70 may need to vary the elemental abundances as a function of radius, and higher spacial resolution studies of the gas emission will help to constrain this variation.

Fig. 6 Radial profiles of C2H emission compared to the work of Facchini et al. (2021). In panel a, we show the emission from the preferred disk model of Fig. 4c, which has C/O = 1.01. The peak flux is slight overestimated, but it reproduces the peak location better than models with lower C/O. In panel b we show the same disk model, but with O/H = 2.525 × 10−4, so that C/O = 0.4. The peak flux is suppressed and the peak moves inward. In this model, the inner disk dominates the C2H flux. In panel c, we set C/H = 1.1 × 10−4 and O/H = 1.0 × 10−4 such that C/O = 1.1. We find that the predicted flux greatly overestimates the observed C2H flux.

Fig. 7 Same as Fig. 6, but for the model in which C2H is artificially removed in the inner 40 AU. The C/O = 0.4 model b does not show a large change, but in both of the carbon-rich models, we find that the inner disk contributes strongly to the C2H emission. Both the C/O = 1.01 (a) and C/O = 1.1 (c) models show similar C2H flux profiles.

4.1.4 Summary of the chemical modelling

In this section, we have tested forward DALI models of disks using the derived outer disk parameters of Keppler et al. (2018, 2019) with a grid of models that varied the properties of the inner disk. These models have shown that the outer disk parameters reproduce the C¹⁸O peak flux and peak location adequately. In addition, the peak ¹²CO emission is well reproduced, but is shifted with respect to the observed peak. This shift is consistent with previous findings of transitions disks (see for example van der Marel et al. 2016), where the gas and dust cavity edges are shifted with respect to each other (although they are the same in our model), and might be linked to a more gradual decline in gas density at R_cav than we modelled here.

Using C₂H emission as a tracer of chemistry, we find that both the outer and inner disk contribute to its flux. For marginally carbon-rich models, the flux peak is shifted towards the outer disk, which is more consistent with the observed radial profile. Furthermore, we showed that for C/O = 1.1, the C₂H flux is far too strong to be consistent with the observed emission, suggesting that the carbon and oxygen abundances of the disk are nearly equal, but slightly favours carbon. In our preferred model, we find that the inner disk is too bright in C₂H, which may suggest a difference in the volatile C/O or C/H between the inner and outer disk. This chemical difference might be linked to the sequestering of either element into larger trapped dust grains, as argued in Sturm et al. (2022). Because of the high temperatures in the inner disk, the exact composition of the gas there can be directly probed by the James Webb Space Telescope (JWST).

In comparison to other studies of high carbon protoplanetary disk chemistry, we computed the integrated flux of C₂H and ¹³CO (whose azimuthally averaged distribution is shown in the appendix). We find that ${\log }_{10}\left(L_{{\text{C}}_2\text{H}}\right)=7.63$, in units of mJy km s⁻¹ pc², and ${\log }_{10}\left(F_{{13}_{\text{CO}}}/F_{{890}_{\text{μm}}}\right)=1.15$, which places our preferred model in a similar range of disks as was studied by Miotello et al. (2019). For the total gas mass of the preferred model (at an age of 5 Myr), 1.5 × 10⁻³ M_⊙, we find that our predicted C₂H flux is consistent with the C/O models of Miotello et al. (2019), which are above unity. The column density of C₂H is higher than 10¹⁴ cm⁻² across the whole disk on average, which is similarly consistent with the models of Bosman et al. (2021), in which C/O is above unity. Because of what we find in comparing the radial distribution of the line emission and the comparisons to other works studying the C₂H flux in protoplanetary disks, we favour a marginally carbon-rich model for the current properties of the PDS 70 disk.

Fig. 8 Midplane gas density profiles (solid lines) used for planet formation compared to the current (at 5 Myr) density profile (green). The wind (orange) and viscous (blue) models are related to the turnback models presented in Table 4 and Eq. (4) at 1 Myr. The dashed lines show the gas temperature through each of the models. We assume that the disk temperature is dominated by the direct irradiation of the host star.

4.2 Planet formation in the PDS 70 disk

4.2.1 Disk model at the time of planet formation

In Fig. 8, we show the two gas density profiles (solid lines) determined by computing the evolution of the PDS 70 disk through viscous (blue) and disk wind (orange) evolution. Viscous disks must transport material outward to carry angular momentum away from the host star. As a result, there is less gas in the outer region (R > 45 AU) of the disk at a younger age than there is now (green). Wind-driven disks, on the other hand, transport angular momentm vertically and do not spread (Tabone et al. 2022). As a result, its outer disk is more massive than both the outer gas disk today and the viscous model.

The dashed lines in Fig. 8 show the midplane gas temperature through each of the models. The temperature profiles are similar between each of the disk models, with the exception of the model representing the disk now. While disk models generally cool as they age, we neglected viscous heating in our chemical models, so that the disk now is slightly warmer because the gas density has dropped, allowing for stellar irradiation deeper into the disk. By 1 Myr, our assumed age of the viscous and wind models, most evolving analytic models find that heating is dominated by direct stellar irradiation outside of a few AU (Chambers 2009).

4.2.2 Growth of planets b and c

Figures 9a,b show the results of planet formation in the wind-and viscous-driven disk models, respectively. In each figure, the faded points denote synthetic planet formation tracks that fail to generate PDS 70-like planets, while the regular points are the successful simulations. The general evolution of synthetic planets is from the bottom and right of the figure towards the top left (shown with the black arrows; outward migration is not possible in our model), and thus the PDS70-like planets begin their evolution farther from PDS 70 than their current orbit.

The masses and orbital radii of the PDS 70-like synthetic planets are within twice the error of the measured values of PDS 70 b and c at a simulation age of 4 Myr. This simulation age corresponds to a system age of 5 Myr because we assumed that planet formation begins at a system age of 1 Myr, and it is consistent with the system properties inferred by Keppler et al. (2019). Because of our choice of proximity within twice the uncertainty of the currently known properties of the planets, in the viscous model (Fig. 9b), the PDS 70b-like planet starts near the current orbit of PDS 70b and migrates slightly inward. This synthetic planet, however, is the only planet that has a similar mass and orbital radius as the current planets of PDS 70.

The drop in planet mass at all time steps in the outer regions of the disk is related to our choice of initial mass for the embryos. As the initial radius of the embryos moves outward, the dust mass decreases (because there is less disk exterior to its initial radius), and thus it starts at a lower mass than embryos beginning closer to the host star. The most extreme example of this is that planets forming inward of ~30 AU in the wind model and ~10 AU in the viscous model have started sufficiently massive that their gas accretion already underwent a runaway by the first time step (60 kyr) shown. In the outer disk, the embryos are not sufficiently heavy to undergo runaway gas accretion by 60 kyr, but they do so later in the simulation.

As a result of the choice of system age, we find that the wind-driven model can reproduce the mass and orbital radii of both planets better than the viscous model. This is mainly linked to the fact that the viscous model requires that the density structure at 1 Myr is more centrally concentrated than in the wind model. As a result, there is less material (both gas and dust) available in the outer disk for the construction of both the core and the gas envelope in the viscous model. Regardless, we can identify a few of our synthetic planets in both models that contain enough mass and orbits near to the current orbit of PDS 70 planets to investigate their resulting atmospheric chemistry.

In the wind-driven case, we overestimate the mass of PDS 70b by a factor of about two. We note here that as in previous works, the planet formation of each synthetic planet was handled separately, and thus the formation of planet c does not have any bearing on the formation of planet b. In reality, the growth of planet c likely impacted the gas flow from the disk to planet b, unless it formed much earlier than planet c. Our model is ot sensitive to differences in the formation start time, but this indicates other interesting lines of research for future studies of this system.

Fig. 9 Evolution of synthetic planets in the disk model given by disk-wind evolution (a) and viscous evolution (b). The colour of each point denotes a particular timestep since the begining of the simulation (t0 = 1 Myr). The faded circles are synthetic planets that do not result in PDS 70-like planets at a time of 4 Myr after the beginning of the simulation. The squares show the successful PDS 70b synthetic planets, and the diamonds show the successful PDS 70c planets. The black circles and error bars are the values for PDS 70 b and c, as shown in Table 1. The grey lines connect the points from individual synthetic planets. The wind model can reproduce both planets b and c, although it over estimates the mass of planet b. The viscous model cannot reproduce planet c and underestimates the radial location of planet b. The inset of panel b shows the collection of successful synthetic planets for both the viscous and wind models. See the text for more details.

4.3 Chemistry in the atmospheres of PDS 70b and c

We computed the number of carbon- and oxygen-carrying molecules that are available to planets b and c planets. We investigated two possible scenarios that are involved in setting the global C/O in the disk. The first was that the processes that lead to a high C/O (discussed above) occurred very early in the lifetime of system PDS 70, and thus the global C/O is the same during the era of planet formation (~1 Myr). On the other hand, the processes that deplete oxygen and some carbon in the disk might have occurred later in the disk lifetime, possibly after the planets had accreted most of their gas. In this scenario, the global C/O would be the same as the stellar value ~0.4.

In Table 5, we show the average C/O in all synthetic planets whose mass and orbital radii correspond to that of the PDS 70 pair of planets (within twice the uncertainties). Where multiple synthetic planets contributed to the average (i.e. in the wind model), we show the range of C/O from each individual planet to an additional decimal place. The spread is small among separate planets. As expected, the planets resulting from the higher C/O models have higher C/O in their atmospheres. In the wind-driven model, where a satisfactory planet c can be found, the PDS 70b planet has a slightly lower C/O than planet c. Planet b results in a slightly different atmospheric C/O in the wind and viscous models. In both cases, small differences in the formation history and the underlying chemical model lead to their varied atmospheric C/O.

In Fig. 10, we show a representation of the CO ice line in our chemical models along the disk midplane. While the two-dimensional structure of the chemistry is important (see for example Cridland et al. 2020a), the midplane abundances can be illustrative for the purpose of the following discussion.

The traditional definition of an ice line is the point in the disk at which the temperature is such that the gas and ice phases of a species have equal abundances. Observationally, this is difficult to constrain, and hence molecular tracers are used to infer the location of ice lines. For CO, N₂H⁺ is a popular molecular tracer because its main destruction pathway involves reactions with CO (Qi et al. 2013a; van’t Hoff et al. 2017). In TW Hya, N₂H⁺ has been used to infer the location of the CO ice line (Qi et al. 2013b).

The CO ice line is located between 15 and 30 AU along the midplane. Additionally, the grey bars in the figure show the initial location of the embryos that led to the best PDS 70-like planets. They are, from left to right, PDS 70b in the viscous model, PDS 70b, and PDS 70c in the wind model. By coincidence, we find that planet b in the viscous model begins its formation in the range of radii where the abundance of CO is near its minimum. The planets forming in the wind model form in a region in which the CO abundances is higher. This higher abundance follows from higher interstellar radiation in the outer disk, which results in higher dust temperatures and thus less efficient freeze-out.

As argued in Cridland et al. (2020b), the gas tends to be more carbon rich than the ices (also see Öberg et al. 2011), and thus the relative amount of gas and solid abundance impacts the atmospheric C/O. In the outer disk, when CO begins to freeze out either as pure CO ice or more likely as other oxygen-rich ice species (see for example Eistrup et al. 2018), the bulk of the carbon and oxygen can become trapped in the ices. This is particularly important for the PDS 70b planet-forming in the viscous model and can help explain its lower C/O compared to the wind model planets.

For the two planets forming in the wind-driven model, we find that their final masses are similar, and thus they acquire very similar atmospheric C/O. The (very) small difference in atmospheric C/O (at the level of 1%) may be due to the slight difference in their starting location and the resulting differences in the available CO in the gas phase. More available gaseous CO will lead to slightly higher atmospheric C/O.

In the low C/O disk model, we find that both PDS 70b and c show nearly equal atmospheric C/O in the two disk models. Unlike in the high C/O models, planet formation in the outer disk appears to be less constrained by chemical gradients in the disk than in the planet formation history of the planets.

Fig. 10 Midplane distribution of molecular abundances for CO and its tracer, N2H+, for the viscous and wind models at 1 Myr. The CO abundance dips for both at the same position as a steep increase in N2H+, which is indicative of the CO ice line. The three grey bands show the range of initial embryo radii that result in the synthetic planets b and c in the viscous and wind models.

5 Discussion

5.1 Why study young systems?

There are several billion years between the point at which a typical exoplanet protoplanetary disk evaporates and the measurement of its atmospheric chemistry. Several dynamical processes, such as gravitational scattering between other planetary bodies (i.e. the Nice model; Gomes et al. 2005; or for hot Jupiters, Beaugé & Nesvorný 2012) and/or scattering by passing stars (Shara et al. 2016; Hamers & Tremaine 2017; Wang et al. 2020, 2022) might change the planet’s orbital radius. These processes, however, should have no impact on the atmospheric chemistry if the planet has already formed fully, and thus the aforementioned framework can still be used to interpret observations.

Over timescales of billion years, the mass loss driven by photoevaporation can potentially change the chemistry of the upper atmosphere (Yelle 2004; García Muñoz 2007; Murray-Clay et al. 2009; Owen & Jackson 2012). In the PDS 70 system, however, these processes have limited effect on the chemistry of the two embedded planets. Firstly, they orbit very far from their host stars compared to planets found in the well-known hot Neptune desert (Owen & Lai 2018), at orbital periods shorter than a few days, which are the main focus of photoevaporation studies. Furthermore, while heavy elements have been observed to be in the evaporating winds of exoplanets (Fossati et al. 2010; Sing et al. 2019), they usually remain coupled to the outflow (Koskinen et al. 2013), which would maintain their relative abundances in the remaining atmosphere (Hunten et al. 1987). Finally, the effect of photoevaporation on the chemical structure of giant planets likely does not play an important role in planets more massive than Saturn (Mordasini et al. 2016; Fossati et al. 2018).

A final source of long-timescale evolution that has the potential to change the observable chemical abundances of the atmospheres away from their primordial values are chemical reactions between the atmosphere and the planet core. One possible direction for chemical evolution is through the envelope-induced core erosion that could transfer heavy elements from the core into the atmosphere via convection (Stevenson 1982, 1985; Guillot et al. 2004; Soubiran et al. 2017). The study of this topic represents a rapidly developing field, but many studies have argued that the cores of gas giant planets consist of a diffuse outer and inner core with steep chemical gradients that suppress adiabatic convection in favour of less efficient heat and compositional transport from the core (Stevenson 1985; Chabrier & Baraffe 2007; Leconte & Chabrier 2012; 2015, 2016, 2018; Wahl et al. 2017; Moll et al. 2017).

In the other direction, differences in the average condensation temperatures of refractory (silicon, magnesium) and volatile-bearing (oxygen, carbon) species causes these species to rain out at different altitudes. This can thus cause the chemical gradients in the planetary atmospheres and their observable chemical abundances to change as the planet loses its accretion energy (Stevenson et al. 2022). These internal processes can lead to difficulties in interpreting the chemical structure of old planets in the context of planet formation, depending on which molecular tracers are used. Young planets, on the other hand, still have most of their accretion heat and have not had enough time for many of these mixing processes to proceed. They thus provide an excellent test bed for studies linking the chemical properties of exoplanet atmospheres to planet formation physics.

5.2 Overview of the PDS 70 system

The PDS 70 system offers a unique look at the planet formation process. While the star’s metallicity is similar to that of the Sun, and so we would expect a stellar C/O that is near solar, the detection of many carbon-rich molecular features suggests that the current disk C/O is higher, perhaps above unity. The chemical modelling presented here seems to support this, but at low C/O, the inner disk contributes strongly to the C₂H flux. This might be caused by two separate assumptions regarding the chemistry in the disk, which we discuss Sect. 5.3.

Computation of the classic core-accretion scenario in the younger version of the PDS 70 disk has shown that super-Jupiter-mass planets at the orbital radii of the known planets are possible, but this is very sensitive to the accretion history of the disk. Viscously accreting disks that have begun to become more centrally compact struggle to accumulate enough material in the outer disk to build the two giant planets. Here we assumed that the initial mass of the planetary embryos is the higher of the pebble isolation mass given by Bitsch et al. (2018) and the total quantity of solid mass exterior to the planet’s starting location, assuming a gas-to-dust ratio of 100. As a result, the initial core mass of our synthetic planets struggles to be high enough to draw down significant quantities of gas, particularly in the viscous model.

We considered two scenarios that might drive the disk C/O away from the assumed stellar value. The first scenario is that the responsible physical or chemical process occurs earlier in the disk life than the era of planet formation, and thus the planets form in the same chemical environment that we see today (i.e. high C/O). The second scenario imagines that the process leading to high C/O occurs late in the disk life, at least after the bulk of planet formation has already occurred, and thus the planets form in a low C/O environment. As expected, the resulting atmospheric C/O are sensitive to this choice: synthetic planets growing in the low C/O environment have a low atmospheric C/O (effectively stellar), while in the opposite case, the resulting C/O is super-stellar.

The measured C/O of PDS 70b and c are currently not well constrained, but the models of Wang et al. (2021) predict an atmospheric C/O for PDS 70b of between 0.54 and 0.7 depending on their choice of model. The sensitivity in the PDS 70c planet is lower than for PDS 70b, but using the same models, Wang et al. (2021) find a range of C/O = 0.49–0.65. Clearly, the current observations of the PDS 70 planets are not sufficiently sensitive to confidently distinguish between different models, and future more sensitive observations are required.

5.3 ISM abundances of carbon

In our chemical model, we assumed an ISM abundance of carbon and modified the abundance of oxygen to set our disk C/O. We note that previous work exploring the link between C₂H flux and volatile C/O, particularly, Miotello et al. (2019) and Bosman et al. (2021), found that they must deplete the carbon and oxygen abundance by a factor of 100 relative to the ISM in order to better match the observed and modelled column densities and fluxes. We find that when we depleted the carbon and oxygen abundance, we completely lost our C¹⁸O flux that comes from the inner edge of the outer dust ring (see Appendix B). We can therefore say that a global depletion of carbon and oxygen is not consistent with the observations of Facchini et al. (2021).

In Fig. B.3, we show the most abundant carbon carriers for the model in Fig. 4a and a model with a similar setup, but with carbon and oxygen abundances reduced by a factor of 100 (i.e. consistent with Miotello et al. 2019) and a gas mass enhanced by a factor of 100. When the carbon and oxygen abundances are reduced, the primary carbon carrier becomes CH₄, which freezes out onto the dust grains. This can be understood from a chemical kinetic perspective as follows: The production of CO scales with the abundance of both carbon and oxygen, and thus its production scales with n_C ṅ_O $n_{\text{X}}^2$, where X denotes the gas metallicity. Thus, when the abundances of carbon and oxygen are reduced by a factor of 100 each, the production of CO is slowed by a factor of 10 000. CH₄, on the other hand, is produced at a rate that scales with n_X because it only contains one heavy element. It can therefore become the dominant carbon carrier in the midplane of the disk at 5 Myr.

Based on this, it is possible that local depletions (or enhancements) of the carbon and oxygen occur in the PDS 70 disk. For example, there is a very obvious dust trap outwards of the orbital position of PDS 70c that might be locally enhanced in carbon and oxygen ices because of the ongoing flux of volatile-rich ices there. The rest of the outer disk might be subsequently depleted in heavy elements. Furthermore, the lack of many optically thin lines in the inner disk might suggest that it too is depleted in carbon and oxygen, particularly given that our chemical models consistently predict a high C₂H flux there, even at lower C/O. Further radiation-dependent studies of the chemistry in the PDS 70 disk are needed, along with upcoming high-resolution ALMA data (Facchini et al. in prep.) and new data from the JWST.

5.4 Choice of chemical network

We followed the work of Miotello et al. (2019), using their chemical network, which was first developed by Visser et al. (2018). The network contains 64 construction and destruction reactions for C₂H with a variety of reaction partners. The carbon chains are a maximum of two carbon atoms (i.e. C₂H, C₂H₂, and C₂H₃), which can result in erroneously large abundances of these species.

As an illustration, Wei et al. (2019) computed the chemistry in a protoplanetary disk for C/O greater than unity. They generally found that the carbon that was not incorporated into CO tended to inhabit long chain hydrocarbons and/or cyanides such as HCN in their inner disk (r < 5 AU). If these longer chain hydrocarbons were available in our chemical model, the inner disk chemistry might have shifted away from C₂H, but would have largely not affected the emission in the outer disk. We leave this investigation to future work.

5.5 Core accretion versus gravitational instability

We have used the standard picture of core accretion to build our synthetic planets. As discussed above, this planet formation scheme struggled at times to build planets similar to PDS 70c, at least in the case of a disk that evolved via viscous evolution. We neglected the other popular planet formation mechanism that often leads to the generation of massive planets: gravitational instability (GI; i.e. Boss 1997).

Very few astrochemical studies of GI currently exist, but Ilee et al. (2017) explored the chemical evolution of a pair of gravitationally unstable clumps in hydrodynamic simulations. They found that the clumps that survived throughout their whole simulation evolved to have the same atmospheric C/O as was initialised in their disk model. We might therefore expect that if PDS 70b and c formed through GI, then they would have the same C/O as the disk.

If this is true, then it would be difficult to differentiate between the different formation models for PDS 70c with C/O alone because its C/O seems to replicate the disk C/O here. In this case, we may need to include an additional tracer of planet formation such as the solid-to-volatile ratio (Schneider & Bitsch 2021), or the nitrogen abundance (Bosman et al. 2019; Bitsch et al. 2022). For PDS 70b, enough information may be encoded in the atmospheric C/O that can be used to differentiate the two formation mechanisms, at least if the disk C/O is high at the time of its formation.

In the future, it would be useful to compute a similar experiment as was done in Ilee et al. (2017) to reproduce the PDS 70 pair of planets. A numerical simulation like this would also help us to understand how the two planets impact their mutual growth and chemical evolution. In this work, we ignored any mutual interaction between the two planets in their growth, which certainly oversimplifies the physical system.

6 Conclusion

We have modelled the chemical and physical structure of the PDS 70 disk in order to understand the environment in which planets PDS 70b and PDS 70c formed. We find that the physical model of Keppler et al. (2018, 2019) adequately reproduces the C¹⁸O and ¹²CO flux, and there is evidence of a tenuous inner disk. The carbon-to-oxygen ratio in the volatiles of the PDS 70 disk is likely marginally above unity, based on the flux of C₂H emission. The inner disk (r < 10 AU), however, is too bright in C₂H, which might suggest that it is depleted in carbon and oxygen relative to the ISM-like abundances that we assumed in this work. The current outer disk (at 5 Myr) is consistent with high abundances of carbon and oxygen, and we showed that if we deplete these elemental species, we cannot reproduce the flux of C¹⁸O that comes from the inner edge of the outer dust ring.

To understand how the disk would have looked at the stage at which planet formation began (at 1 Myr), we used analytic prescriptions for the evolution of the surface density and critical radius under the assumption that the disk evolved through viscous or MHD disk winds. The disk models resulting from the different driving mechanisms result in slightly different surface density profiles at an age of 1 Myr. The main difference between the two disk models is the quantity of material in the outer disk, which affects the formation of the synthetic planets at approximately the location in which the PDS 70 planets currently orbit.

We used a simple prescription for the growth of the PDS 70 planets, initialising their mass with the higher of the pebble isolation masses and the remaining total solid density exterior to the initial orbital radius of the embryo. We accreted the gas and accounted for the collection of carbon and oxygen from the gas-and volatile-rich ices. Because it is unclear when the chemical processing of the volatiles occurred during the evolution of the disk, we tested two scenarios. If the disk is chemically similar to what it is now during the early era of planet formation, the planets have a high C/O, while in the opposite case, when the initial disk composition resembled that of the host star, the planets tended to have stellar C/O. The lower C/O models may be a better match to our current understanding of the planetary C/O, but the planetary C/O is only weakly constrained by observations, which makes comparisons difficult.

The PDS 70 system represents a fantastic environment for studying the link between planet formation and the natal disk of the planet. In particular, understanding the chemical properties of the disk offers a unique opportunity to understand how giant planets acquire higher-mass elements such as carbon and oxygen. This study will benefit greatly from improved observational programs of the disk at submillimeter wavelengths with ALMA and from the infrared with the JWST, along with better constraints on the atmospheric C/O within the planets. Fortunately, all of these studies are well on the way, and so PDS 70 will continue to be an excellent system for studying planet formation for years to come.

Acknowledgements

Thanks to the anonymous referee for their comments which greatly improved the clarity of this work. Thanks to the ExoGRAVITY team for stimulating discussion and observations that triggered this study. Astro-chemistry in Leiden is supported by the Netherlands Research School for Astronomy (NOVA). A.C. and M.B. acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (PROTOPLANETS, grant agreement no. 101002188). E.F.v.D. acknowledges support from the A-ERC, grant agreement no. 101019751 MOLDISK. S.F. is funded by the European Union under the European Union’s Horizon Europe Research & Innovation Programme 101076613 (UNVEIL). This project made use of the following software: Astropy (Astropy Collaboration 2022), SciPy (Virtanen et al. 2020), NumPy (van der Walt et al. 2011), and Matplotlib (Hunter 2007).

Appendix A Planet formation details

Appendix A.1 Gas accretion

In this work, we focused primarily on the overall chemistry of the atmospheres of planets PDS 70b and c. As a result, we ignored the initial build up of the planetary core under the assumption that the core does not contribute significantly to the bulk chemistry in the atmosphere. The gas accretion is limited by a number of different mechanisms depending on the current evolutionary stage of the planet. While the planet is still embedded within the protoplanetary disk (i.e. has not yet opened a gap) we assumed that the gas accretion rate was limited by either the Kelvin Helmholtz timescale (KH) or the Bondi timescale, which ever is longer.

The KH timescale is related to the rate at which the collapsing envelope can release its gravitational potential energy as heat. Its timescale has the functional form of (Alessi & Pudritz 2018) $$t_{\text{KH}}={10}^7\text{yr\hspace{0.17em}}{\left(\frac{{\text{M}}_{\text{plnt}}}{{\text{M}}_{\oplus }}\right)}^{-2},$$(A.1)

where the exponents are determined by comparing population synthesis models of planetesimal core formation to populations of known exoplanets. In principle, the KH timescale depends directly on the opacity of the collapsing envelope, but this is not well constrained. The best-fit exponents of Equation A.1 include variations in envelope opacity when compared to the known exoplanet population.

The Bondi radius describes the region of the disk in which the gas is likely to be captured by the planet. The kinetic energy of a gas parcel is lower than its gravitational potential energy relative to the planet. The planet can thus not accrete more gas than is available within the Bondi radius, or it can be resupplied by viscous processes in the disk. The accretion timescale associated with accreting gas through the Bondi radius is (D'Angelo et al. 2010) $$t_B={\left(C_B\text{Ω}\right)}^{-1}\left(\frac{M_{*}}{a^2\text{Σ}}\right){\left(\frac{a}{H}\right)}^{-7}{\left(\frac{M_{\text{plnt}}}{M_{*}}\right)}^{-2},$$(A.2)

where a is the orbital radius of the planet, Ω is the Kepler orbital frequency, and C_B ≃ 2.6 is a constant meant to match this simple prescription to full three dimensional hydrodynamic simulations. The growth of the planet during the embedded phase is thus Ṁ_plnt = M_plnt/max(t_KH, t_B).

When the planet becomes sufficiently massive, its gravitational influence begins to dominate the surrounding gas over the viscous torques and gas pressure forces. At this point, it opens a gap locally in the protoplanetary disk. The criterion for opening a gap in the disk is a planet of mass M_plnt such that (Crida et al. 2006) $$\frac{3}{4}\frac{H}{R_H}+\frac{50}{qℛ}\lesssim 1,$$(A.3)

where R_H = a(M_plnt/3M_*)^1/3 is the Hill radius of the planet, q = M_plnt/M_* is the planet-to-star mass ratio, ℛ = a²Ω/v is the Reynolds number, and v = αc_sH is the gas viscosity under the standard α-prescription of Shakura & Sunyaev (1973).

When the gap opens, the geometry of the accretion flow changes. This is mainly due to the fact that while the gravitational influence on the gas at the disk midplane is strong, the gravitational force is necessarily weaker straight above the mid-plane. As a result, the gas can leak across the edge of the gap and lose its pressure support. It effectively free-falls towards the mid-plane of the disk, and the planet. This gas motion, often called ‘meridional flow’ (for ex. in Morbidelli et al. 2014), was found in high-resolution studies of CO gas velocities in HD 163296 to coincide with the expected location of an embedded giant planet (Teague et al. 2019). Furthermore, in numerical studies of disk gas hydrodynamics around giant planets, meridional flows have been found to contribute a large fraction (up to ~ 90%; Szulágyi et al. 2014) of the gas flux into the planet's region of gravitational influence.

Morbidelli et al. (2014) outlined the rate of gas accretion onto a growing planet that had recently opened a gap. Broadly speaking, it is limited by the delivery of material to the outer edge of the gap, the disk accretion rate, but they described a cycling of material driven by gas falling through a meridional flow into a decreting circumplanetary disk and back to the outer edge of the planet-induced gap. This recycled gas returns to hydrostatic equilibrium with the rest of protoplanetary disk gas and can return to the meridional flow and the growing planet. This process evolves on the dynamic rather than the viscous timescale, and thus, their mass accretion rate into the gap follows (Morbidelli et al. 2014) $${\dot{M}}_{\text{gap}}=8\pi v\left(r/H\right){\text{Σ}}_{\text{gas}}.$$(A.4)

This accretion rate is faster than the disk equilibrium mass accretion by a factor of 8/3(r/H), which, depending on the disk scale aspect ratio, can be higher by between one and two orders of magnitude.

The efficiency with which the rate of mass accretion into the gap is transferred into an accretion rate onto the planet depends on the local flow around the planet. Batygin (2018) proposed that the magnetic field of the young planet acts to deflect incoming gas into the circumplanetary disk, reducing the accretion and growth efficiency. Cridland (2018) derived the connection between the magnetic field strength and the mass accretion efficiency, resulting in the scaling $$\frac{{\dot{M}}_{\text{plnt}}}{M_{\oplus }/yr}=\frac{4}{3^{3/4}}{\left(\frac{R_0}{R_H}\right)}^2{\left(\frac{M_{\text{plnt}}}{M_{\oplus }}\right)}^{-2/7}{\left(\frac{{\dot{M}}_{\text{gap}}}{M_{\oplus }/yr}\right)}^{3/7},$$(A.5)

where the constant $R_0={\left({\pi }^2/2{\mu }_0{ℳ}^4/G{M}_{\oplus }^2/yr\right)}^{1/7}$ has units of length and depends on the magnetic moment $ℳ=B{R}_{\text{plnt}}^3$ of the (assumed) magnetic dipole of the planet. We assumed that the young planet has a (constant) magnetic field strength of 500 Gauss, which is two orders of magnitude above that of Jupiter, but lower than the typical magnetic field strength of ~1000 K brown dwarf stars. Interpolating the magnetic field strength in this way is consistent with our general understanding of the geo dynamo and solar dynamo (Christensen et al. 2009).

We assumed a constant planetary radius of R_plnt = 2R_Jupiter during this phase because this is the nominal size of young, self-luminous planets with masses greater than that of Saturn, as reported by Mordasini et al. (2015). Cridland (2018) explored the impact of the planet size on the final planet mass due to gas accretion through the above mechanism and found that final masses only varied by a factor of a few, depending on the whether the planet was a cold-start (planet radius equal to the current radius of Jupiter) or a hot-start planet with a radius three times larger than the current Jupiter radius.

We employed standard planet formation methods to compute the growth of the synthetic PDS 70 planets. For simplicity, however, we kept the disk physical and chemical properties constant throughout the planet formation process. This ignores the fact that the disk generally cools and becomes less dense as a function of time. At later times, as the gas is accreted onto the host star, the surface density of the gas drops and the gas accretion rate onto the planet can be reduced through equation A.2. The drop in mass accretion rate will be a particular problem in the viscous model because it already predicts low-mass giant planets. The wind model, on the other hand, may better predict the masses of planets b and c if the accretion rate is lower at later times.

Appendix A.2 Planet migration

Our gas accretion simulations began with a planetary core mass that was set to the smaller of the pebble isolation masses (Bitsch et al. 2018) and the total dust mass exterior to the initial radii. These masses are a factor of a few lower than the gas-gap opening mass (Crida et al. 2006), and thus type I migration affects the overall results of our work only little. We kept the initial core stationary until it reached the gas-gap opening mass, at which point, we allowed it to migrate via type II migration (Lin & Papaloizou 1986). Our choice of ignoring the type I phase (i.e. Ward 1991) of planet migration is related to our ignoring the chemical impact that the planet core might have on the chemical properties of the atmosphere. The location in which the core is built is therefore less important to the overall chemical properties of the atmosphere, and thus, we ignored type I migration and the early build up of the core for simplicity.

When the young planet reached a mass that satisfied Equation A.3, we began to evolve its orbital radius through standard type II migration. In type II migration, the planet has opened a gap in the gas that effectively changes the way in which the protoplanetary disk can transport its angular momentum. Because of its connection to the angular momentum transport, which is typically assumed to be due to viscous evolution in gas, we typically assumed that type II migration evolves on the viscous timescale, $${\dot{r}}_{\text{plnt}}=r_{\text{plnt}}/t_v,$$(A.6)

where the viscous time $t_{\nu }=r_{\text{plnt}}^2/\nu$, and v = αc_sH follows the standard α-prescription of Shakura & Sunyaev (1973). The gas viscosity v depends on the gas scale height H and the gas sound speed c_s³. At a certain point, type II migration stalls when the planet exceeds a critical mass of M_crit = πr²Σ, the total mass of the disk gas inward of the current orbital radius of the planet. After the planet passes this mass, its migration time becomes t_v → t_v(1 + M_plnt/M_crit).

Updated hydrodynamic models of type II migration have shown that the inner and outer disks are not perfectly separated by the planet-induced gap, and some gas crosses the gap (Dürmann & Kley 2015). The gas crossing the gap can also be accreted by the growing planet (Dürmann & Kley 2017), which reduces the efficiency of the gap crossing. Robert et al. (2018) showed that even given these complications, the rate of type II migration continues to be proportional to the gas viscosity, which justifies the use of Equation A.6 in the face of more complex hydrodynamic processes. The proportionality is related to the fact that if the planet migrates on a timescale that is shorter than the viscous timescale, the gas ahead of the migrating planet will be compressed by the planetary torques, while the space behind the planet becomes evacuated. As a result, the inner and outer torque would be strengthened and weakened, slowing the migration rate.

Appendix B Depleting carbon and oxygen abundance

As mentioned in the main text, a common finding in studies focused on C₂H flux is that the carbon and oxygen abundances need to be depleted by up to a factor of 100 times the ISM value in order to match to CO isotopologue and H₂O observations. The depletion is expected to occur because of the freeze out of CO and H₂O onto dust grains that subsequently settle to the mid-plane, hiding the bulk of the volatile abundance of carbon and oxygen. Additional chemical reactions transform CO into other species such as CO₂, CH₃OH, and other hydrocarbons (Bosman et al. 2018b; Krijt et al. 2020). We have tested how this depletion effects the flux of both C¹⁸ O and C₂H in order to better understand our chemical picture.

Figure B.1 shows the C¹⁸O flux from a series of models that were meant to test whether depleting the carbon and oxygen abundances is consistent with observations. In panels (a) and (c), the abundances are depleted, but the density was kept the same as in the model presented in Figure 4(a). The C¹⁸ O flux at the edge of the dust ring is completely lost. In order to possibly recover the flux, we enhanced the gas density in the disk by the same factor as the carbon and oxygen abundances were depleted. These tests are shown in panels (b) and (d). They generally show that gas enhancements cannot recover the missing flux caused by depleting the carbon and oxygen abundances.

Fig. B.1 Test of depleting carbon and oxygen abundances in the PDS 70 disk using C18O emission. Each panel represents a different model, and their relevant parameters are noted. Panels (a) and (c) use the fiducial critical surface density as in the preferred model in the main text. In panels (b) and (d), the critical densities are enhanced by the same factor as the abundances of carbon and oxygen were depleted. A depleted abundance is clearly inconsistent with observations.

In figure B.2 we show the C₂H flux for same disk models as presented in figure B.1. The C₂H flux is not greatly impacted by a depletion of a factor of 10 of the carbon abundance compared to the fiducial disk model (panel a), and it is even weakened when the gas density is enhanced by a factor of 10 (panel b). The reduced flux is related to the equal increase in the dust density in the models in which we enhanced the critical surface density because we kept the dust-to-gas ratio constant when we enhanced the gas density. The additional dust acts as an additional source of opacity, particularly at the edge of the dust ring. This weakens the C₂H flux.

Fig. B.2 Same as Figure B.1, but for C2H emission. Unlike in the C18O case, the emission of C2H can remain strong even with depleted abundances. This is consistent with the finding of Miotello et al. (2019).

The observations are not consistent with a global carbon and oxygen abundance depletion by a factor of 100, which is the depletion reported by Miotello et al. (2019). Here we find a drop in C₂H flux by nearly an order of magnitude (panel c), and more so when the gas (and dust) density is enhanced (panel d). Because we focused on a model with effectively no inner disk because we used the disk model with δ_gas = δ_dust = 10⁻¹⁵ (Figure B.1(a)), the outer disk is inconsistent in general with a global depletion of carbon and oxygen.

A global depletion of the carbon and oxygen abundance is not consistent with the observed radial profile of Facchini et al. (2021), but it might help to explain the unusually high C₂H flux from the inner disk of our model (< 6 AU here). If the inner disk is locally depleted in carbon and oxygen, then it would not emit in C₂H as brightly as we predict. This depletion might be caused by the inefficient transport of material from the outer disk by the elements that are sequestered in the dust and are trapped, where observations seem to prefer ISM-like abundances of carbon and oxygen. A future study of the radial distribution of gas and dust, along with their volatile abundances, would be a useful way of furthering our understanding of the PDS 70 system.

To understand the sudden drop in C¹⁸ O flux when the elemental abundances of carbon and oxygen are reduced, we show the primary carbon carriers for a pair of models in figure B.3. The figure shows the relative number of the most abundant carbon carrier in the volatile and ice (labelled with a ‘J’) phases in both models. When the abundances of carbon and oxygen are reduced by a factor of 100 compared to the ISM (right panel), the most dominant carbon carrier in the disk near the dust ring (50-70 AU) is CH₄, which remains frozen onto the dust grains near the midplane. We understand this shift in the context of the formation rate of CO, which depends on the number density of both carbon and oxygen. CH₄ only depends on the number density of carbon alone, and when the abundances of carbon and oxygen are reduced, the reduction in production rate that follows is therefore less drastic for CH₄ than it is for CO.

Fig. B.3 Most abundant species in the model presented in figure 4(a), (left panel) and the model presented in figure B.1(d), (right panel) for a disk age of 5 Myr. Molecules with a ‘J’ denote species that are frozen onto dust grains. The abundances are integrated between 50 AU and 70 AU (i.e. near the millimeter dust ring) and from the midplane to z/r = 0.43. It therefore includes the main reservoir of CO that contributes to the C18O flux. It is clear that when the abundances of carbon and oxygen are reduced, the main reservoir of carbon is shifted towards hydrocarbons that are frozen onto the dust.

Appendix C Radial profile of¹³CO J = 2-1

In figures C.1 we show the radial profiles of the ¹³CO emission generated by the preferred model. This line was not used to select the preferred model, but we include it here for completeness. The location of the peak flux is slightly shifted outward, similarly to ¹²CO, while the peak flux is slightly overestimated. This suggests that its emission is slightly optically thick, but not to the same extent as the ¹²CO flux.

Fig. C.1 Comparison between the observed 13CO flux and that computed by the preferred model. 13CO appears to be partially optically thick because its peak flux is shifted outward relative to the observed flux, similarly to the 12CO flux shown in figure 5.

References

All Tables

All Figures

Fig. 1 Examples of possible models of the gas and dust distributions to be fed into DALI. We show the number volume density defined as ${\text{n}}_{\text{gas}}={\sum }_{\text{gas}}/\left(\sqrt{2\pi }H\mu m_{\text{H}}\right)$ for the gas. For the dust, ndust is defined to ensure the correct dust-to-gas mass ratio. The solid and dashed curves show the effect of depleting the inner disk by constant factors of δgas and δdust, and the orange curve shows the impact of adding a planet-induced gas gap with a prescribed depth of δplanet.

In the text

Fig. 2 Vertical distribution of the gas and dust at a radius of ~80 AU outward of the outer edge of the dust cavity. The dashed orange line represents the different distributions of the large and small dust (Eq. (3)). The large dust grains are distributed closer to the midplane because they tend to settle with respect to the smaller grains. As a result, the dust-to-gas ratio is slightly higher near the midplane than far above the midplane.

In the text

	Fig. 3 Radial profile observed by Facchini et al. (2021) using ALMA. The grey band shows the approximate location of the (mm) dust ring according to Keppler et al. (2019). The line fluxes peak in slightly different positions, with the optically thin tracers peaking near the edge of the dust cavity, and the optically thick tracers peaking inside the cavity.
In the text

Fig. 4 Radial profiles of C18O emission for three different disk models compared to the observational data of Facchini et al. (2021). The grey band shows the width of the major axis of the beam. In each case, the relevant parameters are listed in the top right corner of the figure. From left to right: in model a, the inner disk is almost completely removed, and there is effectively no gas or dust from Rcav inward. In model b, the inner disk is replaced, with the noted depletions of gas and dust. Finally, in model c, the same inner disk parameters are used, but between Rgap and Rcav, the gas is depleted by a further factor.

In the text

	Fig. 5 Radial profile of the 12CO emission according to the model presented in Fig. 4c.
In the text

Fig. 6 Radial profiles of C2H emission compared to the work of Facchini et al. (2021). In panel a, we show the emission from the preferred disk model of Fig. 4c, which has C/O = 1.01. The peak flux is slight overestimated, but it reproduces the peak location better than models with lower C/O. In panel b we show the same disk model, but with O/H = 2.525 × 10−4, so that C/O = 0.4. The peak flux is suppressed and the peak moves inward. In this model, the inner disk dominates the C2H flux. In panel c, we set C/H = 1.1 × 10−4 and O/H = 1.0 × 10−4 such that C/O = 1.1. We find that the predicted flux greatly overestimates the observed C2H flux.

In the text

	Fig. 7 Same as Fig. 6, but for the model in which C2H is artificially removed in the inner 40 AU. The C/O = 0.4 model b does not show a large change, but in both of the carbon-rich models, we find that the inner disk contributes strongly to the C2H emission. Both the C/O = 1.01 (a) and C/O = 1.1 (c) models show similar C2H flux profiles.
In the text

Fig. 8 Midplane gas density profiles (solid lines) used for planet formation compared to the current (at 5 Myr) density profile (green). The wind (orange) and viscous (blue) models are related to the turnback models presented in Table 4 and Eq. (4) at 1 Myr. The dashed lines show the gas temperature through each of the models. We assume that the disk temperature is dominated by the direct irradiation of the host star.

In the text

Fig. 9 Evolution of synthetic planets in the disk model given by disk-wind evolution (a) and viscous evolution (b). The colour of each point denotes a particular timestep since the begining of the simulation (t0 = 1 Myr). The faded circles are synthetic planets that do not result in PDS 70-like planets at a time of 4 Myr after the beginning of the simulation. The squares show the successful PDS 70b synthetic planets, and the diamonds show the successful PDS 70c planets. The black circles and error bars are the values for PDS 70 b and c, as shown in Table 1. The grey lines connect the points from individual synthetic planets. The wind model can reproduce both planets b and c, although it over estimates the mass of planet b. The viscous model cannot reproduce planet c and underestimates the radial location of planet b. The inset of panel b shows the collection of successful synthetic planets for both the viscous and wind models. See the text for more details.

In the text

	Fig. 10 Midplane distribution of molecular abundances for CO and its tracer, N2H+, for the viscous and wind models at 1 Myr. The CO abundance dips for both at the same position as a steep increase in N2H+, which is indicative of the CO ice line. The three grey bands show the range of initial embryo radii that result in the synthetic planets b and c in the viscous and wind models.
In the text

Fig. B.1 Test of depleting carbon and oxygen abundances in the PDS 70 disk using C18O emission. Each panel represents a different model, and their relevant parameters are noted. Panels (a) and (c) use the fiducial critical surface density as in the preferred model in the main text. In panels (b) and (d), the critical densities are enhanced by the same factor as the abundances of carbon and oxygen were depleted. A depleted abundance is clearly inconsistent with observations.

In the text

	Fig. B.2 Same as Figure B.1, but for C2H emission. Unlike in the C18O case, the emission of C2H can remain strong even with depleted abundances. This is consistent with the finding of Miotello et al. (2019).
In the text

Fig. B.3 Most abundant species in the model presented in figure 4(a), (left panel) and the model presented in figure B.1(d), (right panel) for a disk age of 5 Myr. Molecules with a ‘J’ denote species that are frozen onto dust grains. The abundances are integrated between 50 AU and 70 AU (i.e. near the millimeter dust ring) and from the midplane to z/r = 0.43. It therefore includes the main reservoir of CO that contributes to the C18O flux. It is clear that when the abundances of carbon and oxygen are reduced, the main reservoir of carbon is shifted towards hydrocarbons that are frozen onto the dust.

In the text

	Fig. C.1 Comparison between the observed 13CO flux and that computed by the preferred model. 13CO appears to be partially optically thick because its peak flux is shifted outward relative to the observed flux, similarly to the 12CO flux shown in figure 5.
In the text