CO Depletion in Protoplanetary Disks: A Unified Picture Combining Physical Sequestration and Chemical Processing

In this work, we focus on the region exterior to r = 20 au, as we are interested in the behavior around and outside the CO snowline. When setting the disk gas 2D density and temperature structure (both assumed constant in time), we use the same approach as detailed in Krijt et al. (2018, Section 2.1), with minor changes to some of the global disk parameters. Specifically, focusing on a disk around a Sun-like star, we have set r_c = 80 au, ${M}_{\mathrm{disk}}/{M}_{\star }=0.1$, and the midplane temperature ${T}_{1\mathrm{au}}=135\,{\rm{K}}$, making the disk slightly smaller, twice as massive, and slightly warmer than the disk in Krijt et al. (2018). The gas density and temperature structure is illustrated in Figure 1.

Zoom In Zoom Out Reset image size

At t = 0, we assume that all grains are $0.1\,\mu {\rm{m}}$ in size and are present everywhere in the disk at a dust-to-gas mass ratio of $({\rho }_{{\rm{d}}}/{\rho }_{{\rm{g}}})={10}^{-2}$. For this setup, $\approx 225\,{M}_{\oplus }$ of dust (excluding ices) is located beyond r = 20 au at the beginning of the simulations.

Our focus is on understanding five major carbon- and oxygen-carrying molecules: H₂O, CO, CO₂, CH₄, and CH₃OH. These species have been chosen because together, they represent the dominant carbon- and oxygen-bearing species in the colder regions of the outer disk.

For these five molecules, we solve a set of (semi-) analytical differential equations, which include the following CO processing pathways:

• 1.  
  Successive hydrogenation of CO ice, ultimately leading to the formation of CH₃OH;
• 2.  
  CO reacting with OH on grain surfaces, leading to the formation of CO₂ ice;
• 3.  
  Gas-phase CO reacting with He⁺, leading to the formation of CH₄.

These pathways, and the details of how they are approximated here, are discussed in more detail in Bosman et al. (2018b, Section 2.3) and in the Appendix.

The initial abundances for the five species are listed in Table 1. These initial conditions are based on Bosman et al. (2018b, Table B.1), and correspond to the majority of the carbon being in CO, while the oxygen is initially split between CO and H₂O (see also the Appendix). For a discussion of the validity of these assumptions, see Bosman et al. (2018b, Section B.4.2).

The timescales on which the chemical reactions take place is set by the cosmic-ray ionization rate, which we set to ${\zeta }_{\mathrm{CR}}={10}^{-17}\,{{\rm{s}}}^{-1}$ for our standard model. The surface densities we are interested in do not exceed $100\,{\rm{g}}\,{\mathrm{cm}}^{-2}$, justifying the use of a constant ζ_CR (Umebayashi & Nakano 1981; Padovani et al. 2018). We do not include X-rays as a source of ionization, the presence of which could potentially speed up CO destruction in very low-mass ($\leqslant 0.003{M}_{\odot }$) disks (Schwarz et al. 2018).

Desorption energies for the five major species are listed in Table 1. The balance between freeze-out and desorption on small dust grains is assumed to reach steady state for all five major molecular species during each timestep (as in Bosman et al. 2018b). This is a valid assumption in most regions of the disk, but can potentially lead to the models underpredicting the amount of gas-phase CO in the cold outer disk regions, especially when these become depleted in dust (Semenov et al. 2006; Krijt et al. 2018). However, for the relatively weak turbulence strength used in this work, these effects will be minor. For pebble-size particles that drift through the midplane snowline, time-dependent sublimation is included, which will result in pebbles releasing CO molecules over an extended radial region (see also Piso et al. 2015).

For the initial abundances, the midplane CO snowline (defined as the location where 50% of the present CO is frozen out) is located at r ≈ 34 au, where the midplane temperature equals T(r) ≈ 23 K.

The evolution of solids is handled in a way that is very similar to the approach described in Krijt et al. (2018), with small grains being treated as a fluid on a grid, while the dynamics of pebble-size particles are followed using a swarm of representative tracer particles.

For the small grains, the ices present on small grains, as well as trace gas-phase vapor species, diffusive transport is included using the approach of Krijt et al. (2018; see also Ciesla 2009). The strength of the turbulence is controlled via the dimensionless α-viscosity parameter (Shakura & Sunyaev 1973), which is assumed to be connected to the gas diffusivity via ${D}_{{\rm{g}}}=\alpha {c}_{s}H$. Motivated by Flaherty et al. (2015) and Teague et al. (2016), we set α = 10⁻⁴, an order of magnitude lower than in the fiducial runs of Krijt et al. (2018).

The formation and evolution of pebble-size particles is handled in a way that is very similar to the approach described in Krijt et al. (2018). For increased flexibility and more control over the pebble behavior, however, we simplify the approach of Krijt et al. (2018, Section 3) in two ways:

• 1.  
  We assume that all pebbles that form in a single simulation have the same size s_p, independent of their formation location.
• 2.  
  We approximate the dust-to-pebble conversion timescale as

  $$\begin{eqnarray}&&{\tau }_{{\rm{d}}\to {\rm{p}}}=\displaystyle \frac{{a}_{1}}{{\rm{\Omega }}}{\left(\displaystyle \frac{{\rho }_{{\rm{d}}}}{{\rho }_{{\rm{g}}}}\right)}^{-1},\end{eqnarray} \tag{ 1 }$$

  with Ω the local Keplerian frequency, and a₁ a constant.

When pebbles form, they inherit the ice composition of the small grains present in that particular grid cell. The local abundance of small grains is lowered accordingly, and the effect of the decreasing small-dust-to-gas ratio is included in the chemical network.

Vertical settling, radial drift, and turbulent transport are included following Krijt et al. (2018). We include the sublimation of volatile species present on pebble-size particles if and when they move to disk regions with a higher temperature (i.e., as they drift inward), but otherwise, pebbles are assumed to be chemically inert because of their low surface-to-mass ratio. In practice, this means that small dust and pebbles can have quite different ice compositions, even when they are present at the same location.

In the fiducial model, we set s_p = 1 mm throughout the disk, and as in Krijt et al. (2018), we imagine that further growth is frustrated by a combination of bouncing and radial drift that removes particles before they can gain more mass. We ignore collisional fragmentation, which is reasonable for the outer disk regions and the fairly low value of α = 10⁻⁴ (Birnstiel et al. 2012; Misener et al. 2019). In more complex models, the exact maximum particle size is a function of location, local dust-to-gas ratio (for the drift limit), and particle composition, but a constant s_p = 1 mm suffices for our purposes and reproduces the typical pebble size outside ∼20 au in more complete models to within a factor of a few (e.g., Stammler et al. 2017; Krijt et al. 2018). For our disk model, pebbles of this size have Stokes numbers of $\mathrm{St}=4\times {10}^{-3}$, 0.02, and 0.05 in the midplane at r = 20, 50, and 100 au, respectively.

Starting with all the dust present as submicron grains and the initial abundances listed in Table 1, the five molecular tracers (in either vapor and/or ice form), the small dust, and pebble tracers are evolved forward in time. During every iteration,

• 1.  
  Gas-phase and grain-surface chemistry is advanced in each cell using the simplified scheme outlined in the Appendix.
• 2.  
  The new sublimation and freeze-out balance is calculated, and the partitioning of volatile species between ice and vapor is updated in each cell.
• 3.  
  The exchange of small dust grains (carrying ices) and gas-phase species between neighboring cells is calculated.
• 4.  
  The positions of the pebble tracer particles are updated, taking into account vertical settling, radial drift, and turbulent diffusion.
• 5.  
  The sublimation of ices carried by pebbles that have moved to warmer disk regions is calculated.
• 6.  
  In each cell, a fraction of the locally available dust mass is stochastically converted into pebble particles. The ice on newly formed pebbles has the same composition as the dust in that cell.

Steps 1–6 are repeated until t = 3 Myr have passed, and the timestep for this operative splitting routine is set by the shortest diffusion timescale across a single cell (Krijt et al. 2018, Equation (17)). Similarly, we set a₁ = 10, which results in coagulation timescales similar to those found in 1D models (see Section 3.3).

The starting point is a model in which we only include chemical processing and freeze-out/sublimation of molecular species. Dust coagulation and material exchange between adjacent grid cells is not included (see Table 2). The results of this model are presented in Figure 2, in which the following locations, corresponding (roughly) to the temperature and density combinations shown in Figure A1, have been highlighted: (A) inside the CO snowline; (B) outside the CO snowline; (C) main freeze-out zone in the outer disk; (D) the warm molecular layer (e.g., Aikawa et al. 2002). The observed behavior can be directly compared to the chemically more sophisticated, but similarly static studies (e.g., Walsh et al. 2010; Henning & Semenov 2013; Bosman et al. 2018b; Schwarz et al. 2018, and others). Like Schwarz et al. (2018) and Bosman et al. (2018b), we find that the processing of CO in the WML (e.g., location D) is inefficient on timescales of some million years. In colder regions, where CO is either fully or partially frozen out on grain surfaces, more rapid processing can occur, which typically results in CO₂ ice (around T ∼ 20–25 K, see location A) and CH₃OH or CH₄ ice (where $T\lt 21{\rm{K}}$, see locations B and C) becoming the dominant carbon carrier (see also Section 2.2 and Bosman et al. 2018b, Figures 5 and 7). Similar findings across a variety of disk temperature and density profiles led Zhang et al. (2019) to conclude that chemical processing alone is not an efficient way of removing the warm CO gas that dominates the emission that is readily probed by ALMA.

Zoom In Zoom Out Reset image size

The next step is to include material transport through turbulent diffusion. Generally, diffusion will act to smear out existing spatial gradients in the concentrations of gas-phase and ice-phase species. In the absence of dust coagulation, and thus large variations in the dust-to-gas ratio, abundance gradients are set mainly by spatial variations in temperature, gas density, and the radiation field (e.g., Willacy et al. 2006; Ciesla 2009; Semenov & Wiebe 2011).

The timescale for vertical mixing to act on a scale Δx can be estimated as ${t}_{z}={({\rm{\Delta }}x)}^{2}/{D}_{{\rm{g}}}$. We first focus on vertical transport. For Δx = H and assuming a diffusion coefficient based on a turbulence strength described by α = 10⁻⁴ (Section 2.1), we obtain ${t}_{z}\approx 0.25$ and 1.6 Myr at 30 and 100 au, respectively. These timescales are shorter than or comparable to the chemical processing in the disk midplane (Figures 2 and A1), suggesting that turbulent transport and chemical processing of CO are intertwined and that the CO abundance is sensitive to transport processes.

The results of the chemistry and diffusion (DIFF) model are shown in Figure 3. When we compare the behavior to the static CHEM model (Figure 2), the effect of turbulent transport is indeed evident. First, the ice species produced near the midplane are effectively transported to the disk's upper layers, resulting in CH₃OH and CH₄ being present at high z/H after 3 Myr. Second, the narrow band of CO₂ that is seen hugging the T ∼ 20–25 K region in Figure 2 has disappeared, with CO₂ ice only being the dominant carbon carrier for r < 50 au. In this region, however, fast vertical mixing has resulted in CO₂ being dominant throughout the entire vertical column. Qualitatively similar behavior has been found in earlier works by Semenov & Wiebe (2011) and Furuya & Aikawa (2014).

Zoom In Zoom Out Reset image size

Finally, Figure 3 reveals a decrease in the gas-phase CO abundance in the WML of about an order of magnitude after 3 Myr. This can be understood in terms of vertical mixing: While local processing of CO in the WML is not efficient (see Figure 2), the addition of diffusion allows gas-phase CO molecules to be cycled through the colder midplane, where (1) they rapidly stick to grain surfaces, and (2) can be processed to form less volatile species such as CH₃OH. These products, and the portion of CO ice that escaped processing, will eventually return to the WML, but overall the process is asymmetric and results in the loss of gas-phase CO from the upper regions on timescales of some million years.

In the seminal work of Semenov & Wiebe (2011), the depletion of gas-phase CO in the outer disk was found to be inefficient (see also Semenov et al. 2006), its column density changing by no more than a factor ∼3–5 over $5\,\mathrm{Myr}$ when turbulent diffusion was included. While a direct comparison is difficult because of several differences in the setup (disk density and temperature structure, etc.), we attribute this difference to two reasons: (1) Semenov & Wiebe (2011) did not include H and H₂ tunneling. This results in a slower conversion of CO into CH₃OH on the grain surface, while making it easier to re-form CO as the incorporation of O into H₂O is also suppressed (see also Bosman et al. 2018b); and (2) Semenov & Wiebe (2011) used a ratio of diffusion to binding energy ${E}_{\mathrm{diff}}/{E}_{\mathrm{bind}}=0.77$ that is substantially higher than our value of 0.3, also reducing the efficiency of grain surface chemistry.

For completeness, we include here a model with dust evolution and transport only, without taking into account chemical processing (PEBB, Figure 4). This setup is essentially identical to the one used in Krijt et al. (2018), with the caveat that we have made several minor changes to the way the pebble size and coagulation timescale are calculated (e.g., Section 2.5).

Zoom In Zoom Out Reset image size

As the result of pebble formation and vertical settling, the local dust-to-gas ratio decreases in the disk surface layer and increases in the midplane (e.g., Krijt et al. 2018, Figure 6). At the same time, the inward drift of pebble-size particles results in a radial mass flux of solids throughout the disk. Typically, at several 10 s of au, this flux builds up rapidly, peaks around ${10}^{5}\,\mathrm{yr}$, and then decreases slowly on $\sim {10}^{6}\,\mathrm{yr}$ timescales (Krijt et al. 2018, Figure 12), although the detailed behavior varies with location and depends on the (initial) disk properties and details of the dust coagulation process (e.g., Birnstiel et al. 2012; Lenz et al. 2019).

Figure 7 shows the time evolution of the total mass of various pebble reservoirs present in the PEBB simulation as well as the cumulative pebble mass that has drifted through the inner boundary. Over the course of the PEBB simulation, $\sim 220{M}_{\oplus }$ worth of pebbles can be seen to drift inward of r = 20 au, with about 10% of them doing so in the first $0.2\,\mathrm{Myr}$. Pebble fluxes of similar magnitudes play an important role in modern planet formation theories, providing a means for planetary embryos to accumulate mass quickly through pebble accretion (Johansen & Lambrechts 2017; Ormel 2017). For example, roughly $\sim 150{M}_{\oplus }$ of pebbles are required to grow the cores of the solar system's giant planets within 3 Myr (Lambrechts & Johansen 2014).

The consequences of pebble formation and dynamics on gas-phase CO abundances have been discussed in Krijt et al. (2018), and we recover the same behavior: (1) the CO abundance in the WML drops by a factor of a few as ices are sequestered on pebble-size particles in the midplane; and (2) the sublimation of rapidly drifting and ice-covered pebbles results in a plume of gas-phase CO interior to the CO ice-line. The appearance of this plume will vary with time as the radial pebble flux evolves (see Figure 7), while its shape and peak CO abundance depend sensitively on the strength of the turbulence and its relation to the diffusion coefficient (Stammler et al. 2017, Figures 7 and 8).

Finally, we present the results of the FULL model including chemistry, diffusive transport, coupled with dust coagulation and pebble dynamics in Figure 5.

Zoom In Zoom Out Reset image size

Because the processing of CO takes place on timescales of some million years, the picture at early times (e.g., at 0.3 Myr) is comparable to the one presented in Figure 4, with a similar plume of gas-phase CO developing centered around r ∼ 30 au, and some CO disappearing from the WML around r ∼ 50 au as it becomes sequestered on pebbles in the midplane. As in the PEBB model without chemistry, the plume inside the CO ice-line persists for 3 Myr, but the magnitude at the end of the simulation is significantly reduced. This can be seen clearly in Figure 6, where we compare radial profiles of the gas-phase CO abundance (in the midplane and in the WML) for the four models described in Sections 3.1–3.4. There are two reasons for this reduction. First, CO can be processed locally, mainly forming CO₂ at these pressures and temperatures (see Figure 2). Second, as a result of chemical processing on grain surfaces, pebbles that form and arrive at the CO ice-line late in the simulation will carry a significant fraction of their carbon in the form of CH₃OH and CH₄ rather than CO (we return to this point in Section 4.6). These species have their ice-lines further in, and will not be released to the gas-phase as these pebbles cross the CO ice-line.

Zoom In Zoom Out Reset image size

In the region inside and around the CO ice-line, i.e., for r < 50 au, the inward migration and evaporation of CO-ice-coated pebbles results in gas-phase CO always being the dominant carbon carrier. The situation in the outer disk, however, especially for $t\gtrsim 1\,\mathrm{Myr}$, more closely resembles the picture presented in Figure 3, with CH₄ and CH₃OH becoming the dominant carbon carriers. Focusing on the disk surface layer exterior to r ∼ 50 au, we see that the two CO removal mechanisms encountered earlier (sequestration on pebble surfaces, and turbulent transport followed by chemical processing in the midplane) seem to exacerbate the depletion of gas-phase CO, with local abundances dropping by up to 2 orders of magnitude over the modeled 3 Myr. Such depletion factors are reached exclusively in the FULL model.

In the framework presented here, nothing stops the pebbles from drifting inward in the models that include dust coagulation (Figures 4 and 5). Thus, the only way to prevent a resurgence of gas-phase CO just inside the snowline from occurring is to have chemistry outpace pebble drift, resulting in pebbles carrying carbon in a less volatile form (i.e., in CH₄ and CH₃OH in our calculations). This is illustrated in the left panels of Figure 8, where we have varied ζ_CR from ${10}^{-18}\,{{\rm{s}}}^{-1}$ (FULL-CR18), ${10}^{-17}\,{{\rm{s}}}^{-1}$ (FULL, the fiducial model), and ${10}^{-16}\,{{\rm{s}}}^{-1}$ (FULL-CR16). We note that the overall ionization structure of disk systems as a function of time is highly uncertain (e.g., Cleeves et al. 2013, 2015; Padovani et al. 2018) and the values used here represent a reasonable range for exploration. The pebble behavior has not been altered in these runs, so the flux of solids crossing the CO snowline is identical⁸ to the one depicted in Figure 7. In the high cosmic-ray model, the region with enhanced CO abundance around 20–40 au persists for less than a million years, with gas-phase CO being virtually absent after 3 Myr. The disappearance of the peak has two causes; (1) gas-phase CO is processed locally, mostly converted into CO₂ (see Section 3.1), and (2) chemical processing in the outer disk midplane means that pebbles arriving at times $t\gtrsim {10}^{5}\,\mathrm{yr}$ contain very little CO ice (see also Section 4.6 and Figure 12).

Zoom In Zoom Out Reset image size

Alternatively, the balance between chemical processing and radial drift can shift when pebble migration becomes less efficient or is completely halted. For example, the formation of larger bodies (i.e., planetesimals, see Section 4.6) or the trapping of pebbles in radial pressure bumps should reduce the radial pebble flux (e.g., Pinilla et al. 2020). In disks where massive planets have already formed, deep gaps in the gas surface density profiles can trap pebble-size particles outside the planet's orbit with a near 100% efficiency (Pinilla et al. 2012; Bae et al. 2019). Indeed, pebble accumulations observed with ALMA (e.g., by the DSHARP survey) are regularly located exterior to the inferred CO ice-line (Huang et al. 2018; Long et al. 2018; Andrews 2020), clearly demonstrating the viability of trapping pebbles and the ices they carry in the outer regions of protoplanetary disks.

We illustrate the impact of these processes in the right-hand panels of Figure 8, where we compare the fiducial model to a case where the drift efficiency is reduced by 90% (FULL-DR90), and 99% (FULL-DR99), respectively. It is evident that the CO depletion pattern in the outer disk WML is identical, as this is set mainly by vertical transport. The decreasing importance of radial motions for pebbles has a clear impact on the CO peak interior to the snowline. For the FULL-DR90 model, a small peak still develops, but it takes several million yers to do so. In the FULL-DR99 model, where pebble drift is essentially stopped (as would be the case when pebbles are very efficiently converted into planetesimals at all radial locations), even the region interior to the CO snowline becomes depleted in CO.

Zoom In Zoom Out Reset image size

Recently, Zhang et al. (2019) derived radial profiles of the (gas-phase) CO depletion factor⁹ by combining spatially resolved CO isotopolog emission with spectral energy distribution (SED) modeling and thermo-chemical modeling for a handful of nearby disks (specifically, DM Tau, TW Hya, HD 163296, and IM Lup). While these disks vary in age, size, and temperature structure, we can qualitatively compare the findings of Zhang and collaborators to the results of our models.

Focusing on the WML region (solid curves) that is located beyond the CO ice-line, Figure 6 indicates that depletion factors between ∼10–100, as inferred for IM Lup and TW Hya, are exclusively reached in the FULL model, which combined chemical processing with turbulent mixing and the sequestration of CO ice in the disk midplane. Even with these processes combined, however, several million years of evolution were needed to reach such severe CO depletion, while IM Lup is believed to be younger, $\sim 1\,\mathrm{Myr}$ (Zhang et al. 2020a). Possible reasons for more rapid CO removal could include a higher ionization rate (Figure 8), higher turbulence in the disk surface layers (Xu et al. 2017; Krijt et al. 2018), or a colder disk in which the CO surface snowline lies farther away from the midplane (Zhang et al. 2019, Figure 4).

For HD 163296, the radial profile derived by Zhang et al. (2019) included the region just interior to the CO snowline, suggesting the presence of a plume of CO vapor, as predicted in our models that included pebble drift (Sections 3.3 and 3.4) and earlier works by Stammler et al. (2017) and Krijt et al. (2018). This local CO plume was further confirmed by the modeling of the ¹³C¹⁸O (2–1) line spectrum of HD 163296 by Zhang et al. (2020b), who found that a radial pebble flux of ∼15–60 M_⊕ Myr⁻¹ is needed to reproduce the elevated C/H ratio interior to the CO snowline. Comparable fluxes¹⁰ are seen in our models that include pebble evolution (see Section 3.3). The observed substructure in the dust continuum emission of HD 163296 (e.g., Isella et al. 2018) might suggest efficient trapping of pebbles; however, Rosotti et al. (2020) argue that even the largest grains currently present in the disk are relatively well coupled to the gas, finding $(\mathrm{St}/\alpha )\sim 20$ for the two rings firmly outside the CO snowline.¹¹ Such observational constraints on the radial pebble flux are very valuable because this quantity plays a key role in shaping the final masses and orbital architectures of both gas giant and more terrestrial planets (e.g., Johansen & Lambrechts 2017; Lambrechts et al. 2019). Last, the limited CO depletion in the outer disk for the relatively old HD 163296 is likely a result of the disk being fairly warm (Zhang et al. 2019; Dullemond et al. 2020), rendering CO processing less efficient (Bosman et al. 2018b; Schwarz et al. 2018).

For TW Hya, CO isotopolog emission reveals a significant amount of gas-phase CO between 5 and 20 au, with abundances significantly higher than in the outer WML, but still depleted relative to ISM levels (Schwarz et al. 2016; Zhang et al. 2017, 2019). Recent studies by Bosman & Banzatti (2019) and McClure et al. (2019) also argue that CO release inside 20 au is minimal, suggesting that carbon is either transported inward in a different form, or that pebble migration has ceased altogether (see Section 4.1). However, as discussed in Kama et al. (2016, Section 6.3), the C/Si ratio of gas accreting onto TW Hya is an order of magnitude above that of typical T Tauri stars, suggesting that some carbon makes it to the inner disk even as refractory elements like Si are locked up in planetesimals or even larger bodies.

We note that while the disk used in this work was quite massive (see Section 2.1), we expect CO depletion to proceed similarly in lower mass disks as the relevant timescales do not depend on disk gas mass directly, but rather on the dust-to-gas ratio (for pebble formation) and the turbulence strength (for vertical mixing). Nonetheless, the detailed temperature structure and ionization environment will set the timescale for chemical processing as well as the locations of snowlines. As such, in depth comparisons to observed protoplanetary disks warrant dedicated models for each individual object (e.g., Zhang et al. 2019).

In the absence of a direct gas mass tracer, gas-phase CO (isotopolog) line fluxes and dust continuum emission originating from the outer disk are frequently used to estimate disk gas masses. A common approach is then to assume ISM-like values for CO/H₂ ≈ 10⁻⁴ and/or ${M}_{\mathrm{dust}}/{M}_{\mathrm{gas}}\approx {10}^{-2}$ (e.g., Williams & Best 2014; Ansdell et al. 2016; Miotello et al. 2016; Pascucci et al. 2016, and many others). Uncertainties in these indirect approaches have recently been discussed in, e.g., Bergin & Williams (2017) and Kama et al. (2020). Indeed, using an alternative method of looking at disk dusk lines for several bright protoplanetary disks, Powell et al. (2019) obtain typical dust-to-gas ratios of $\sim {10}^{-3}$ in the outer disk regions, and total gas disk masses between ∼3–100× higher than those derived from CO.

In our simulations, the amounts of warm CO, dust, and pebbles vary significantly with time as a result of chemical processing and dust evolution. We illustrate the effects that these processes can have on disk mass estimates in Figure 9, which shows the time evolution of the mass-weighted CO/H₂ and solids-to-gas conversion factors for the various models shown in Figure 8. The left panel shows warm CO interior to the midplane CO snowline, while the panel on the right illustrates the warm molecular layer in the outer disk. In both panels, the horizontal axis represents the total (i.e., dust+pebble) solid-to-gas mass ratio integrated over the entire disk. The background colors represent factors of $3\times$ reductions in the CO abundance (blue) or solids-to-gas ratio (red).

Zoom In Zoom Out Reset image size

The trajectories depicted in Figure 9 illustrate that the accuracy of using either CO emission (while assuming an ISM-like abundance) or dust emission (while assuming an ISM-like dust-to-gas ratio) as a total disk mass tracer depends sensitively on the details of the ongoing chemical and dust evolution. For example, in the FULL-DR90 and FULL-DR99 models, the reduced drift efficiency means that the the solid mass hardly changes, while the warm CO in the outer disk drops by a factor $\sim {10}^{2}$ over the modeled 3 Myr. In disks evolving in this fashion, the dust mass would be a more reliable gas mass tracer, as found to be the case in the Lupus star-forming region by Manara et al. (2016). Conversely, in the FULL-CR18 model (see Section 4.1), slow processing of CO and fairly rapid dust drift result in CO emission from the outer regions being a more reliable gas mass tracer (especially at times $\lesssim \mathrm{Myr}$), even if the CO/H₂ conversion factor has to be adjusted by a factor of a few. For the canonical FULL model, both the total solid and warm CO mass in the outer disk are reduced by over an order of magnitude in 3 Myr.

Generally, however, both the warm CO (at least in the outer disk) and total solid masses are decreasing functions of time, implying that studies looking for evolutionary trends in gas disk masses by comparing star-forming regions of different ages have to be careful when using CO/H₂ or dust-to-gas conversion factors that are constant in time, especially if a significant fraction of the disk mass is located beyond the CO snowline.

The processes described described in Section 3 will alter key elemental ratios in different ways. Recognizing that the effect of radial drift on gas and solid elemental abundances in the midplane has been discussed in detail by multiple authors (e.g., Cuzzi & Zahnle 2004; Ciesla & Cuzzi 2006; Estrada et al. 2016; Öberg & Bergin 2016; Booth et al. 2017; Booth & Ilee 2019), we focus here primarily on the warm molecular layer.

4.4.1. Gas-phase C/H and C/O

Figure 10 shows the gas-phase (C/H) and (C/O) ratios after 3 Myr for CHEM, DIFF, and FULL simulations. The behavior seen for (C/H)_gas largely follows the CO depletion pattern, as described for each simulation in Section 3 in more detail.

Zoom In Zoom Out Reset image size

The story for (C/O)_gas is somewhat more involved. Initially, CO dominates the carbon and oxygen reservoirs in the gas, resulting in (C/O)_gas ≈ 1 at t = 0. As a result of chemical processing, this ratio tends to increase with time for the pressures and temperatures we are considering. The reason is that of the chemical products considered here, CH₄ is the most volatile (after CO, see Table 1), its gas abundance pushing (C/O)_gas > 1 as the more oxygen-rich products (e.g., H₂O, CO₂) readily freeze out. This effect is similar to the one described in Schwarz et al. (2019, Section 4.3.3), who reported (C/O)_gas as high as $\sim {10}^{8}$ in the midplane at r = 12 au after $6\,\mathrm{Myr}$, when virtually all CO had been processed.

In our simulations, the magnitude of the increase in (C/O)_gas and the extent of the region where it is seen depend on the transport processes that are included. In the static CHEM model (left panels of Figure 10), only a relatively small region shows (C/O)_gas > 1 because a significant fraction of the initial CO is still present after 3 Myr (see Figure 2).

For the DIFF model, the combination of more efficient CO removal combined with the upward transport of CH₄ formed closer to the midplane results in a large portion of the WML having (C/O)_gas ∼ 10–30. There is a small region in the upper inner part of the disk that has (C/O)_gas < 1, due to sublimation of CO₂ that formed at lower temperatures and was subsequently transported there.

The FULL model also results in carbon-rich gas in the WML, but the increase is not as extreme as for the DIFF model without pebble formation and evolution. There are two reasons for this. First, in the inner disk, the enhancement in CO vapor interior to the CO snowline (see Section 3.4) dominates all other effects, forcing (C/O)_gas ≈ 1. Thus, if radial drift is taking place, we would not expect to see the extremely high (C/O) ratios in the inner disk midplane gas as found by Schwarz et al. (2019). Farther out, e.g., in the WML between r ∼ 50–100 au, the sequestration of a substantial amount of CH₄ ice on settled pebbles likely plays a role in the (C/O)_gas ratio being slightly lower than in the DIFF model, but still significantly ∼3–10.

Designed to be focused on the processing of CO, the reduced chemical network we have used here does not consider further processing of CH₄ and the formation of other, generally less volatile, hydrocarbons. As such, the abundance of CH₄ at the end of our calculations, and the increase of (C/O)_gas associated with its appearance, should be treated as an upper limit.

4.4.2. Total Volatile C/H and C/O

In models without dust coagulation and pebble dynamics, the total volatile (i.e., gas + ice) C/O ratio does not change with time, since there is either no transport between neighboring cells (for the static CHEM model) or because the dynamics of gas-phase species and small ice-carrying grains are essentially identical (for the DIFF calculation). For these models, the initial concentrations described in Table 1 result in ${({\rm{C}}/{\rm{O}})}_{\mathrm{gas}+\mathrm{ice}}\approx 0.5$ and ${({\rm{C}}/{\rm{H}})}_{\mathrm{gas}+\mathrm{ice}}\approx {10}^{-4}$ everywhere in the disk.

When dust evolution is included, however, the dynamics of carbon and oxygen can differ depending on how they are partitioned between gas and ice, leading to spatial and temporal variations in the total elemental ratios (e.g., Öberg & Bergin 2016; Booth & Ilee 2019). Figure 11 shows the total¹² volatile C/H and C/O ratios for the FULL-CR18, FULL, and FULL-CR16 models (described earlier in Section 4.1).

Zoom In Zoom Out Reset image size

The behavior shown by the (C/H)_gas+ice can be understood in terms of dust coagulation. As dust grows, especially in the regions exterior to the CO snowline, a significant amount of carbon becomes locked up in pebbles, reducing (C/H)_gas+ice on a timescale comparable to the pebble formation timescale. This process is slightly less effective in the disk surface layers, where the carbon in CO and CH₄ is in the gas-phase. Inside the CO snowline, sublimation of inward-drifting pebbles (e.g., Section 3.4) leads to a (C/H)_gas+ice ratio that is elevated by a factor ∼2–3 compared to the initial value. For the FULL-CR16 model, the values for (C/H)_gas+ice in the inner regions are lower, as the more rapid chemical processing makes it easier for pebbles to carry away carbon in the form of CO₂, and later CH₄, ice.

The variations in (C/O)_gas+ice are caused by a mix of chemical processing and dust evolution, and can be understood as follows. In the cold outer disk midplane, virtually all carbon and oxygen is frozen out and (C/O)_gas+ice ≈ (C/O)_ice. Even as pebbles form and leave this region, the remaining ice dominates the carbon and oxygen budget, leaving (C/O)_gas+ice relatively unaffected. In the surface layers, and inside and around the CO snowline, the picture is different. Here, a significant amount of C and O is initially found in the gas, and the ice is generally more oxygen-rich than the gas, with the more volatile species (CO and CH₄) having C/O > 1. As the oxygen-rich ice is removed over time, (C/O)_gas+ice slowly increases, and even becomes >1 as the abundance of CH₄ becomes comparable to that of CO. This occurs faster in models with a higher ζ_CR, explaining the increasing (C/O)_gas+ice values going from FULL-CR18, to FULL, and FULL-CR16.

Elevated values (C/O)_gas+ice have been inferred from observation of multiple disks, with (C/O)_gas+ice ranging between 0.8 and 2.0 (Bergin et al. 2016; Cleeves et al. 2018; Bergner et al. 2019; le Gal et al. 2019; Miotello et al. 2019). This is comparable to the range of values found in the FULL-CR18, FULL, and FULL-CR16 models. The models predict a strong correlation between CO abundance and (C/O)_gas+ice ratio in the surface layers, however, which so far has not been observed (Miotello et al. 2019). Again, the lack of further processing of CH₄ could lead us to overpredict the (C/O)_gas+ice in the top layers of the disk. Even so, the dynamical and chemical processes studied here are at least in part responsible for the elevated (C/O)_gas+ice inferred from observations.

Following the seminal work of Öberg et al. (2011), many studies have sought to connect chemical abundances in the gas and solids, as found in the disk midplane, to the compositions of giant planet atmospheres (e.g., Öberg & Bergin 2016; Booth et al. 2017; Booth & Ilee 2019; Cridland et al. 2019). Recently, motivated by the theoretical and observational evidence that massive planets accrete gas primarily from the disk surface layers (Morbidelli et al. 2014; Szulágyi et al. 2014; Teague et al. 2019), Cridland et al. (2020) explored how gas accretion from reservoirs farther up in the disk can alter the giant planet atmospheric composition, generally finding that the final C/O ratio was lower as the result of the accretion of oxygen-rich ices from between 1 and 3 scale-heights. The vertical profiles of the molecular abundances at the basis of the accretion model of Cridland et al. (2020) were based on otherwise static disk model. It would then be interesting to explore how the changes predicted here (e.g., Figure 11) would impact the make-up of giant planets forming in the outer disk. However, only a relatively small fraction of the planet's final mass was accreted from the disk surface layers, especially for planets that (start) forming outside ∼30 au (Cridland et al. 2020, Figure 10). A more significant change in giant planet C/O ratios could come from the accretion of ice-rich pebbles in the disk midplane (e.g., Johansen & Lambrechts 2017). Understanding the balance between radial drift and chemical processing timescales (e.g., Section 4.1) is then key, as pebbles that carry processed carbon in the form of CH₃OH or CO₂ can hold on to their carbon for much longer than pebbles that are covered in more volatile CO ice.

A more direct comparison can be made between the ice mantles in our simulations and the compositions of left-over planetesimals in the solar system and beyond. For example, Eistrup et al. (2019) recently compared constraints on cometary compositions to the time and location dependent ice abundances in their static chemical model, finding that the majority of cometary properties can be explained by a formation location in the vicinity of the CO ice-line.

The most direct connection is arguably possible for small bodies in the Cold Classical Kuiper Belt (CC-KB). This population of planetesimals, orbiting the Sun between ∼42 and 47 au on low-inclination, low-eccentricity orbits, is unique in that it contains some of the most pristine bodies in the solar system (Prialnik et al. 2020). Moreover, the orbits of these bodies are mostly undisturbed by the outward migration of Neptune (Nesvorný 2015), suggesting that their current location is very close to their formation location. Finally, there is compelling evidence that these bodies formed through the gentle collapse of a cloud of pebble-size particles concentrated via the streaming instability (e.g., Nesvorný et al. 2019; McKinnon et al. 2020), implying that their bulk composition at the time of formation closely resembles that of the pebble precursors.

The best constraints on CC-KB object compositions come form the NASA New Horizons fly-by Arrokoth¹³ early in 2019 (Stern et al. 2019). While a body as small as Arrokoth cannot hold on to hypervolatiles like CO (Brown 2012), a surprising finding was the clear detection of CH₃OH ice, coupled with a lack of strong H₂O ice features in the reflected near-IR spectrum (Grundy et al. 2020). While it is unclear to what extent these surface features represent the bulk composition (see the discussion in Grundy et al. 2020), we can attempt to compare these findings to results in our models.

Figure 12 shows the (mass-weighted) composition of all pebble particles present between 42 and 47 au (i.e., the location of the CC-KB) as a function of time for three models with different cosmic-ray ionization rates. For the model with ${\zeta }_{\mathrm{CR}}={10}^{-18}\,{{\rm{s}}}^{-1}$, the grain surface chemistry is slow and the ice mantles of the pebbles are always composed of a large amount of CO ice, with ≲5% CH₃OH ice. For the standard FULL model (middle panel), the highest concentration of methanol ice (relative to water ice) is reached after ∼2 Myr. In the model with a high cosmic-ray rate (right panel), a similar situation is reached already after ∼0.4 Myr, after which the contribution of CH₃OH begins to diminish.

Zoom In Zoom Out Reset image size

The current mass of the entire CC-KB population is estimated to be $\sim 3\times {10}^{-4}{M}_{\oplus }$ (Fraser et al. 2014), and the region probably only suffered a factor of ∼2 dynamical depletion during the outward migration of Neptune (Nesvorný 2015). Thus, even after 3 Myr in our simulations, the region between 42 and 47 au contains enough mass in pebbles to produce a CC-KB like population of planetesimals (see Figure 7).

Looking beyond the solar system, observations of second-generation CO gas (and its derivatives) in old (>10–100 Myr) and otherwise gas-poor debris disks are being employed to constrain the combined CO+CO₂ mass fraction of planetesimals at the top of the collisional cascade (Wyatt 2020, and references therein). For example, Matrà et al. (2017) constrain the combined CO+CO₂ ice mass fraction to lie between ∼0.3%–30% for the planetesimals in HD 181327 and Fomalhaut. As these constraints become tighter and more systems are added, comparing such constraints to models like the ones presented here, and looking for trends with, e.g., stellar luminosity or planetesimal belt location will undoubtedly provide valuable insights.

Last, the interstellar interloper 2I/Borisov appears to be very rich in CO ice compared to solar system comets (Bodewits et al. 2020; Cordiner et al. 2020), which has been taken as evidence for its formation exterior to the CO snowline. In the context of the models presented here, we would argue that if it is indeed CO-rich, 2I/Borisov and the exocomets in, e.g., Fomalhaut not only formed exterior to the CO snowline in their systems, but also early during disk evolution, before significant processing of CO ice on grain surfaces took place.

We end by stressing that in the simulations presented here, the local dust-to-gas ratio in the midplane is not readily enhanced sufficiently to result in planetesimal formation via the streaming instability (see also Krijt et al. 2016, 2018). It appears then that additional physics such as disk photoevaporation (Carrera et al. 2017) or the presence of (short-lived) pressure traps (Lenz et al. 2019) are needed to facilitate planetesimal formation outside of r ≳ 10–100 au.