Ice-Coated Pebble Drift as a Possible Explanation for Peculiar Cometary CO/$\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ Ratios

To model the dynamics of the gas bulk (defined as the bulk hydrogen gas, which experiences no source terms), we follow Lynden-Bell & Pringle [25], which is based on the $\alpha$-disk model of Shakura & Sunyaev [35]. Thus, we have the partial differential equation

in the absence of sources and sinks, where $\Sigma_{\mathrm{gas}}\equiv\int\rho_{\mathrm{gas}}\mathop{}\!\mathrm{d}z$ is the surface density of gas, $\nu$ is the viscosity, $R$ is the distance from the star in the $x$-$y$ plane, and $t$ is time. Viscosity is, in turn, given by

where $\alpha$ is a small parameter of our choosing, typically set to $10^{-4}$ to $10^{-2}$; $c_{s}$ is the local sound speed, given by

where $k_{B}$ is Boltzmann’s constant, $T$ is the local temperature, $\mu$ is the mean molecular weight, and $m_{p}$ is the proton mass; and $\Omega$ is the Keplerian angular frequency,

with $G$ the gravitational constant and $M_{\star}$ the central stellar mass. Equations 1, 2, 3, and 4 completely define the model of the gas bulk given parameters $\mu$, $M_{\star}$, and $\alpha$; the local temperature field $T=T\!\left(R,t\right)$ (see Section II.4; and the initial condition $\Sigma_{\mathrm{gas}}\!\left(R,t=0\right)$.

For the initial condition, we first define the self-similar solution,

with $\Sigma_{c}=20~\text{\text{g}}~\text{\text{cm}\textsuperscript{$-2$}}$, $R_{c}=20~\text{\text{au}}$, and $\gamma=0.5$; for reference, Andrews et al. [6] uses $-1\leq\gamma\leq 1$. Here, $\Sigma_{c}$ is the surface density at radius $R_{c}$, and $\gamma$ determines the slope of the power law part of the solution. Unfortunately, when $R\ll R_{c}$, this solution begins to blow up, which makes it computationally difficult to handle. We use a smooth interpolation between the self-similar solution and a flat, constant surface density profile, given by

as our initial condition. We take $p=5$ and $R_{\mathrm{trans}}=1~\text{\text{au}}$ so that the transition occurs close to the interior of the domain and the transition from the self-similar to the flat profile is not too sharp.

Though we have chosen to work in one dimension, some quantities depend on the local density $\rho$ rather than the surface density $\Sigma$. In these cases, we assume a vertical Gaussian distribution of material,

where the scale height $h_{\mathrm{gas}}=c_{s}/\Omega$.

We consider two solid populations in our model: a small “dust” population with radius $0.1$ µm and a “pebble” population with radius $1$ mm, with mass ratios $90\%$ and $10\%$, respectively. Following Birnstiel et al. [10], we define the surface density evolution for each population by the partial differential equation

in the absence of sources and sinks, where $\Sigma_{\mathrm{solid}}$ is the solid surface density for a single population and $F_{\mathrm{tot}}\equiv F_{\mathrm{adv}}+F_{\mathrm{diff}}$ is the total flux, with contributions from an advective and diffusive part. The fluxes are given by

and

In the above equations, the Stokes number is given by

in the Epstein regime, with $a_{\mathrm{gr}}$ the radius of a single (pebble or dust) grain and $\rho_{\mathrm{gr}}$ the density of the solid material (i.e., silicate). The radial velocity of the solids is given by

where

is the gas velocity and

is the velocity due to the gas pressure gradient. $E_{d}$ is a drift efficiency parameter and $p_{\mathrm{gas}}=\rho_{\mathrm{gas}}c_{s}^{2}$ is the gas pressure. Birnstiel et al. [10] gives more detail on these equations.

Again, we must make some assumption about the vertical distribution of solids to determine $\rho_{\mathrm{solid}}$. We make the same vertical Gaussian assumption as for the gas, but, to simulate settling, we allow the scale height of the pebbles to be a fraction $\xi_{\mathrm{pebbles}}$ of the gas scale height, so $h_{\mathrm{pebbles}}=\xi_{\mathrm{pebbles}}h_{\mathrm{gas}}$. We use $\xi_{\mathrm{dust}}=1$ such that the dust is not settled. For the pebbles, we take $\xi_{\mathrm{pebbles}}=0.1$.

Finally, adsorption, the process by which atoms and molecules stick to a solid surface, and desorption, in which the atoms and molecules leave the surface, must be included as source terms. Hollenbach et al. [24] gives the adsorption timescale as

where $n_{\mathrm{solid}}$ is the local number density of solids, $\sigma_{\mathrm{gr}}=\pi a_{\mathrm{gr}}^{2}$ is the cross-sectional area of a single grain, and $v_{\mathrm{therm}}=\sqrt{8k_{B}T/\pi m}$ is the thermal velocity of the atom or molecule of interest. Inverting the timescale, we find the adsorption rate

per atom or molecule.

For desorption, Hollenbach et al. [24] gives the rate per molecule of ice

where $\nu_{\mathrm{att}}$ is the attempt frequency — the vibrational frequency of the atoms and molecules on the surface — of order $10^{12}$ s^$-1$, and $T_{\mathrm{bind}}$ is the binding energy of the species of interest (4800 K for $\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ and 960 K for CO, Aikawa et al. 2).

Combining Equations 16 and 17, we find the volumetric source terms

To find the appropriate source term for the surface density equations above, we must integrate $s_{\mathrm{gas}}$ vertically and multiply by the species’ mass $m$. We find

and

(Equation 19 is derived in Appendix A.1.) These source terms both have units of g cm^$-2$ s^$-1$ and represent the rates at which the surface density changes due to adsorption and desorption processes, respectively.

Thus, the source terms for surface density equations are given by

and

These source terms encode the rate at which the surface densities of gas and solid species are changing due to the adsorption and desorption chemistry in our model.

The temperature field presents a challenge by itself. Temperature appears in Equation 1 through the viscosity term, and so it contributes to the gas dynamics. Yet the gas dynamics play a role in determining the dust dynamics, which, through radiative transfer from the central star, determine the temperature. In addition, the intrinsic luminosity of the star is expected to change significantly over the timescales presented here [36]. One way to solve this circular problem is through iteration, as in Price et al. [33]. However, that procedure would be more computationally costly when coupled to the code we have described here.

Instead of seeking a self-consistent solution, as in Price et al. [33], we follow a simpler procedure to capture the approximate temperature structure. Noting that the bulk surface density, and therefore dust grain surface density, does not change significantly over time, we use RADMC-3D version 0.41 [20] to compute a temperature structure with a self-similar dust initial condition (i.e., same form as Equation 5), assuming a dust-to-gas ratio of $0.01$ and using the interpolated DSHARP opacities [11]. We note that this procedure is an approximation, since we are not taking into account the evolving dust and pebble surface density, but it provides sufficient accuracy for our proof-of-concept purposes.

Next, we fit a power law $T\propto R^{-\beta}$ to the output from RADMC-3D, limited to the region between 2 au and 20 au to avoid edge effects and unphysical behavior far from the star. Though we run RADMC-3D with two dust populations, the temperatures are virtually equal, so we assume a power law slope of $-0.41$ and appropriate intercept parameter, which reasonably captures the behavior of both populations, and use that same power law for both when solving the differential equations.

To take into account a changing stellar luminosity over time, we use the Siess et al. [36] web server to compute stellar radii $R_{\star}$ and effective temperatures $T_{\mathrm{eff}}$ over the lifetime of the disk. Then, inspired by Chiang & Goldreich [16], Equation 12, we see that the disk temperature scales by a factor $f\propto R_{\star}^{1/2}T_{\star}$. We compute this factor from the isochrons and scale it by the initial value such that $f\leq 1$ at all times, i.e., the disk temperature is decreasing over time, primarily due to radial contraction decreasing the bolometric luminosity of the central star.

Finally, we perform a fit to the two regimes we observe in $f$ — a flat, early-time regime and a sloped, late-time regime — and join the two regimes by smooth interpolation. This interpolation takes the same form as Equation 6, but with a parameter $p=100$ that is more appropriate for this data. See Figure 1 for the parameters in each regime and the final interpolation. Figure 2 shows the resulting temperature that is used in the fiducial model alongside two fixed-temperature models representing the beginning and end state.

To solve Equations 1 and 8 with source terms given by Equations 21 and 22, we require approximations of first and second derivatives in radius. We use a logarithmically-spaced mesh in $R$ and second-order accurate finite difference derivatives estimated with Equations A3 and A4. Where appropriate, we switch to first-order accurate upwind finite difference derivatives.

To advance the solution in time, we use the backward differentiation formula (BDF) implementation in the Portable, Extensible Toolkit for Scientific Computation (PETSc) [7, 8, 9] time stepping (TS) [1] module. We use the PETSc internal colored finite difference Jacobian and solve the resulting linear system with the MUltifrontal Massively Parallel sparse direct Solver (MUMPS) [5, 4].

The system of partial differential equations we finally solve is in nine quantities. The bulk gas, pebble, and dust densities are treated according to Equations 1 and 8 with no source terms. Then, we consider $\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ and CO in gas, as ice on pebbles, and as ice on dust grains by adding the appropriate source term to the right-hand sides of Equations 1 and 8. We evolve the equations to $1$ Myr on the spatial domain $\left[0.5~\text{\text{au}},500~\text{\text{au}}\right]$.

In Figure 3 (first column), we show the behavior of the bulk gas, pebbles, and dust over time and radius in the fiducial model. While the gas behavior shows simple viscous spreading, the pebbles and dust show more interesting behavior. The 1 mm pebbles form a shallow gap-like structure at about 30–100 au. This position coincides with the radius where the Stokes number goes to unity, and thus where the pebbles move fastest. As a result, at smaller radii, the pebbles move inward, and, at larger radii, the pebbles move outward, resulting in a pebble deficit at about 100 au. The 0.1 µm dust forms a similar structure at larger radii.

Figure 4 shows that the ice-coated solids do not universally follow the same trends as the bulk. While $\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ on pebbles and dust forms the gap-like structure near 100 au, there is a second pebble and dust deficit at 1 au, the $\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ snowline, where there is also a rapid increase in $\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ vapor surface density. The behavior of CO is significantly different from that of $\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$. The gap at $30$–$100$ au is much shallower, and only clearly visible at late times. Analogous to $\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$, there is a rapid drop in CO dust and pebble surface density at the CO snowline. Figure 4 already shows a clear change in the CO/$\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ surface density ratio in the outer disk.

In Figure 5, we present the main results of this paper both for the fiducial model and for a small parameter study (see next section). Each panel in the figure shows the evolution of the CO/$\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ ratio in two ways: On the left, we show the variation over time and space on the vertical and horizontal axes, respectively. On the right, in a smaller panel, we integrate over radius in the region shown and show the evolution of the ice mass — total and where CO/$\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ $\geq 1$ — over time. In Figure 5a, we show the predicted ratio of CO/$\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ for our fiducial model, and we find that a maximum ratio near unity is achieved by 1 Myr in the region between about $20$ and $200$ au, and that this feature takes the shape of a funnel when observed in the space-time plane. This enhanced material accounts for an average of $40\%$ of the disk mass. See Table 2 for similar measurements of each model case that follows.

While our fiducial model results are encouraging in explaining anomalous, CO-enhanced comets, we also seek to understand the robustness of this result to changes in disk parameters, relatively unconstrained by observations or detailed simulations. The first parameters of interest are the viscosity parameter and drift efficiency. The viscosity parameter $\alpha$ ultimately sets the diffusion coefficient, which directly controls the gas’s diffusion and the diffusive flux of the solids. Figure 3 (second column) shows that reducing this parameter only has some minor effects on the dust and pebble evolution. The drift efficiency $E_{d}$ influences the coupling of solids to gas but only appears in the dust velocity, and therefore leaves gas motion unchanged. Figure 3 (third and fourth columns) show that increasing and decreasing this parameter dramatically changes the drift and therefore depletion of solids in the outer disk regions.

Figure 5b shows the enhancement of the CO/$\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ ratio in ice for a model with $\alpha$ reduced by a factor of ten compared to the fiducial model. Reducing $\alpha$ makes the viscosity smaller everywhere, which, in turn, amounts to making the diffusion coefficient smaller. Thus, we would expect that disk material would experience less viscous spreading in this case, and, indeed, we see that the characteristic “funnel” shape of the CO-enhanced region in time and space is truncated and does not reach 100 au, while drift still carries material inward towards the star. The amount of CO-enhanced ice is modest, but it is certainly present.

In Figures 5c and 5d, we show the enhancement of the CO/$\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ ratio in ice for the low and high drift models, respectively; for these test cases, we fixed $E_{d}$ at $0.01$ and $0.9$, changing the efficiency of the coupling to the gas pressure derivative. The models achieve about the same maximum CO/$\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ ratio, but with very different fractions of CO-enhanced ice; i.e., the low-drift enhancement feature is visibly smaller in the space-time plane. We immediately see, then, that the efficiency of the radial drift of pebbles and dust is very important for predicting the amount of mass available for making comets like 2I/Borisov and C/2016 R2 (PanSTARRS). We return to this in the discussion section.

The third parameter of interest is the initial CO/$\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ ratio. We test two possibilities in addition to the fiducial model: A high CO/$\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ value of $100\%$ and a low CO/$\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ value of $1\%$. We choose these end-member cases because, while typical comets have CO/$\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ of about $4\%$ [13], they may have $1\%$ or less [29], while the interstellar medium has up to $100\%$ with a large errorbar [30]. Note that the amount of CO has no effect on the bulk dynamics; it only affects the chemical evolution of the disk. In Figures 5e and 5f, we show the chemical evolution of the disk in these two cases. We find that the low CO model is not able to reach a CO/$\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ ratio of unity; see Figure 5e. On the other hand, the high CO model easily reaches values of CO/$\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ $\gtrsim 10$ by about $1$ Myr, as shown in Figure 5f, and this material accounts for a large fraction of the total disk mass.

The effect of a static, low temperature profile and a static, high temperature profile on our disk model is shown in Figures 5g and 5h, respectively. For these cases, we artificially fixed the temperature at its final or initial value, as appropriate (recall that the temperature strictly decreases with time, and see Figure 2 for the radially-dependent structures we adopted). These models reach roughly the same level of CO/$\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ ice enhancement, but the radii where the enhancement occurs are shifted. In the low temperature model, the onset of the enhancement is delayed in time. The high temperature model’s enhanced region is shifted to larger radii because the disk is warmer everywhere, and so ice will desorb off the grains in this model farther out than in the fiducial model. Most importantly, the details of the temperature structure and evolution are not critical for the formation of a substantial amount of CO ice; both static temperature models achieve at least 30% CO-enhanced ice (see Table 2).

Finally, we summarize our results numerically in Table 2 in terms of maximum CO/$\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ ratio, total number of CO-enhanced Halley-mass comets, and the mass fraction of the ice in the disk that is CO-enhanced by the end of the simulation. Most models achieve a CO/$\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ ratio above unity (only the low CO model achieves a lower maximum ratio). The lowest ratio (low initial CO model) is just over $0.1$, and the largest ratio (high initial CO model) is greater than $10$, revealing a monotonic dependence on CO initial abundances. The low drift model produces the next-least amount of CO-enhanced ice. The fraction of ice in the region $\left[5~\text{\text{au}},200~\text{\text{au}}\right]$ that is CO-enhanced is, on average, about $40\%$, but in the high-drift model it is all of $87\%$, indicating that most water ice has been lost from the system due to pebble drift. The maximum number of CO-enhanced Halley-like comets that could be formed in the disks is between $10^{8}$ and $10^{10}$, though this assumes a formation efficiency of 100% from the dust and pebbles and no additional mixing, trapping, or drift. While a large range of values are possible, we emphasize that there is almost always a region where there is some CO ice enhancement relative to $\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ ice.