CHEMISTRY IN AN EVOLVING PROTOPLANETARY DISK: EFFECTS ON TERRESTRIAL PLANET COMPOSITION

The pressure and temperature structure of the disk, necessary for the equilibrium chemistry calculations, is calculated using a theoretical model derived in Chambers (2009), which describes a viscously heated and irradiated disk.

The Chambers model divides the disk into three regions. In the intermediate region the heating of the disk is dominated by viscous dissipation. The surface density is given by

$$\begin{equation} \Sigma (r,t) = \Sigma _{{\rm vis}}\left(\frac{r}{s_0}\right)^{-3/5}\left(1+\frac{t}{\tau _{{\rm vis}}}\right)^{-57/80}, \end{equation} \tag{ 1 }$$

where

$$\begin{equation} \Sigma _{{\rm vis}}=\frac{7M_0}{10\pi s_0^2}, \end{equation} \tag{ 2 }$$

and M₀ and s₀ are the initial mass and outer edge of the disk, respectively. The temperature of the disk in this region is given by

$$\begin{equation} T(r,t)=T_{{\rm vis}}\left(\frac{r}{s_0}\right)^{-9/10}\left(1+\frac{t}{\tau _{{\rm vis}}}\right)^{-19/40}, \end{equation} \tag{ 3 }$$

where

$$\begin{equation} T_{{\rm vis}}=\left(\frac{27\kappa _0}{64\sigma }\right)^{1/3} \left(\frac{\alpha \gamma k_B}{\mu m_H}\right)^{1/3} \left(\frac{7M_0}{10\pi s_0^2}\right)^{2/3} \left(\frac{GM_*}{s_0^2}\right)^{1/6}. \end{equation} \tag{ 4 }$$

κ₀ is the opacity and is taken to be a constant 3 cm² g⁻¹, α = 0.01 is the viscosity parameter, μ = 2.4 is the mean molecular weight, γ = 1.7 is the adiabatic index, and m_H is the mass of hydrogen. The viscous timescale is

$$\begin{equation} \tau _{{\rm vis}} = \frac{1}{16\pi }\frac{\mu m_H}{\alpha \gamma k_B}\frac{\Omega _0M_0}{\Sigma _{{\rm vis}}T_{{\rm vis}}}. \end{equation} \tag{ 5 }$$

In the inner region of the disk, heating is still dominated by viscous dissipation, but the temperature is too high for solids to condense, and consequently the opacity is expected to be several orders of magnitude lower than elsewhere in the disk. The opacity is assumed to follow the power-law form given in Stepinski (1998),

$$\begin{equation} \kappa =\kappa _0\left(\frac{T}{T_{e}}\right)^n, \end{equation} \tag{ 6 }$$

where T_e = 1380 K is the temperature at which the majority of solids evaporate, κ₀ is the opacity at T < T_e, and n is the power-law index, taken to be −14. The surface density in the inner region of the disk is

$$\begin{equation} \Sigma (r,t) = \Sigma _{{\rm evap}}\left(\frac{r}{s_0}\right)^{-24/19}\left(1+\frac{t}{\tau _{{\rm vis}}}\right)^{-17/16}, \end{equation} \tag{ 7 }$$

where

$$\begin{equation} \Sigma _{{\rm evap}} = \Sigma _{{\rm vis}}\left(\frac{T_{{\rm vis}}}{T_e}\right)^{14/19}. \end{equation} \tag{ 8 }$$

The temperature in the inner region is given by

$$\begin{equation} T(r,t)=T_{{\rm vis}}^{5/19}T_{e}^{14/19}\left(\frac{r}{s_0}\right)^{-9/38}\left(1+\frac{t}{\tau _{{\rm vis}}}\right)^{-1/8}. \end{equation} \tag{ 9 }$$

The transition radius between these regions is

$$\begin{equation} r_e(t) = s_0\left(\frac{\Sigma _{{\rm evap}}}{\Sigma _{{\rm vis}}}\right)^{95/63}\left(1+\frac{t}{\tau _{{\rm vis}}}\right)^{-19/36}. \end{equation} \tag{ 10 }$$

In the outer regions of the disk, viscous dissipation does not deposit significant amounts of energy, and so stellar irradiation dominates the heating of the disk. The surface density equation is

$$\begin{equation} \Sigma (r,t) = \Sigma _{{\rm rad}}\left(\frac{r}{s_0}\right)^{-15/14}\left(1+\frac{t}{\tau _{{\rm vis}}}\right)^{-19/16}, \end{equation} \tag{ 11 }$$

where

$$\begin{equation} \Sigma _{{\rm rad}}=\Sigma _{{\rm vis}}\left(\frac{T_{{\rm vis}}}{T_{{\rm rad}}}\right), \end{equation} \tag{ 12 }$$

and

$$\begin{equation} T_{{\rm rad}}=\left(\frac{4}{7}\right)^{1/4}\left(\frac{T_*kR_*}{GM_*\mu m_H}\right)^{1/7}\left(\frac{R_*}{s_0}\right)^{3/7}T_*. \end{equation} \tag{ 13 }$$

The temperature in the outer region is

$$\begin{equation} T(r,t)=T_{{\rm rad}}\left(\frac{r}{s_0}\right)^{-3/7}. \end{equation} \tag{ 14 }$$

The transition radius between the outer irradiated region and the intermediate viscous region is given by

$$\begin{equation} r_t(t) = s_0\left(\frac{\Sigma _{{\rm rad}}}{\Sigma _{{\rm vis}}}\right)^{70/33}\left(1+\frac{t}{\tau _{{\rm vis}}}\right)^{-133/132}. \end{equation} \tag{ 15 }$$

For all calculations in this work we use M₀ = 0.1M_☉, s₀ = 33 AU, T_* = 4200 K, and R_* = 3R_☉. This model corresponds to the second example in Stepinski (1998), which is consistent with a planetesimal-forming disk.

Earlier works (e.g., Bond et al. 2010a; Elser et al. 2012) have used equilibrium chemistry calculations in the protoplanetary disk to predict planetesimal compositions. Each calculation used the pressure and temperature profile of a protoplanetary disk at one point in its evolution, and so the composition of the planetesimals reflects the composition of the disk only at that epoch. The planets built from these planetesimals had compositions that, to first order, were consistent with the bulk compositions of the terrestrial planets of the solar system. As a comparison with these works, we perform similar equilibrium calculations.

Like the aforementioned works, we use the software package HSC Chemistry (version 7.1). Given the pressure, temperature, and elemental abundances, HSC Chemistry will calculate the equilibrium abundances for a list of species by minimizing the total Gibbs free energy of the system using the GIBBS solver (White et al. 1958). The species considered are the same as those used in Bond et al. (2010a) and are listed in Table 1. Each of the solid species was considered a pure substance (i.e., no solid solutions). Calculations were performed between 0.3 and 4 AU in radial extent to match the region considered in the dynamical simulations.

The accuracy of our chemical model is limited by two main factors: (1) The completeness of the list of chemical species considered. This list is limited by both the computational complexity of calculating equilibrium for a large number of species and HSC Chemistry's database. A noticeable shortfall is the lack of carbon in our models throughout the asteroid belt (Figure 1), where carbonaceous meteorites can contain as much as a few percent by mass in carbon. This is expected, though, because most carbon in carbonaceous meteorites is in the form of organic macromolecules (Pizzarello et al. 2006), which are not included in our study. (2) The assumption of equilibrium conditions. It is unknown exactly how much of the solar nebula experienced conditions where chemical equilibrium could occur. However, the material that makes up the asteroids and terrestrial planets of the solar system is believed to have formed under near-equilibrium conditions (Davis & Richter 2005 and references therein).

Zoom In Zoom Out Reset image size

In a more realistic scenario, planetesimals are formed over the course of a few million years. Because the temperature profile of a protoplanetary disk changes significantly during this period, it is necessary to consider the changes in the composition of solid material available for building planetesimals as the disk ages.

Taking a similar approach to Cassen (1996), we start with an initially hot protoplanetary disk and evolve it forward in time while keeping track of the chemical composition of dust, gas, and planetesimals in the disk. In our model, we break the disk up into two components: the gas/dust disk and the planetesimal disk. The planetesimal disk is assumed to not react chemically with the gas/dust disk as only the thin outer shell of a planetesimal is in contact with the gas and dust of the disk. The planetesimal disk is initially empty but builds up over time as planetesimals form from the available solid material. To calculate the composition of these planetesimals, we iterate through the lifetime of the disk. At each time step and in each radial zone of the disk, we perform the following steps.

• 1.  
  Calculate the equilibrium chemical composition of the dust/gas disk using HSC Chemistry.
• 2.  
  Remove a fraction of the solid material from the gas/dust disk and add it to the planetesimal disk. The amount of material that is converted to planetesimals is determined by the planetesimal formation rate, which is described shortly.
• 3.  
  Calculate the radial movement of the gas and dust according to the disk equations described in Section 2.1. In our model we assume that all dust is perfectly coupled to the gas and that planetesimals are large enough that they do not experience orbital migration due to gas drag and thus remain where they formed.
• 4.  
  Recalculate the chemical inventory throughout the disk as it has changed due to radial motion of the disk and formation of planetesimals.

We repeat each of these steps until the amount of material being converted into planetesimals is negligible (about 3 million years in our model). The compositions of the planetesimals at the end of this portion of the simulation are taken as the input composition for the next stage of our model (the N-body simulations of late-stage planet formation).

As solid material grows from dust to planetesimal size, it must cross a size threshold where it goes from being chemically reactive to mostly locked away in the interior of a planetesimal. We will call the rate at which material crosses this threshold the planetesimal formation rate. The planetesimal formation rate in protoplanetary disks is currently unknown. In fact, the formation of planetesimals, in general, is not well understood. A number of mechanisms have been proposed to explain their formation, which can be categorized into two general processes: (1) planetesimals formed from the pairwise collisions and sticking of larger and larger aggregates (Blum & Wurm 2008, and references therein), and (2) they formed from the collective self-gravity of a large mass of small particles (e.g., Johansen et al. 2007; Cuzzi et al. 2008).

Regardless of the mechanism of planetesimal formation, the planetesimal formation rate will vary in time and space due to differing conditions in the disk. For the sequential condensation chemistry model it is essential to know how the planetesimal formation rate varies. For lack of a consensus on the planetesimal formation mechanism, we will use the simple prescription for the planetesimal formation rate used in Cassen (1996), which is consistent with formation from coagulation

$$\begin{equation} \dot{\Sigma }_p \propto \Sigma _{{\rm solid}}\Omega, \end{equation} \tag{ 16 }$$

where Σ_solid is the surface density of the solid material available to build planetesimals at a point in the disk and Ω is the orbital angular speed at that point. The exact form of this expression may change as planetesimal formation models are revised. Nevertheless, a planetesimal formation rate that increases with increasing solid surface density and a decreasing dynamical timescale is a logical starting point.

Despite the fact that only the planetesimals that form interior to 4 AU are considered in further steps of the model (see Section 2.4), it is necessary to track the formation of planetesimals farther out in the disk as well. Planetesimals forming in the outer regions will deplete the dust and gas disk of any elements that are solid, which, when the disk material moves inward, can affect the chemistry in the terrestrial-planet-forming region. For this reason, our model extends to 10 AU, at which point the maximum depletion of any element is <10%.

As discussed previously, the assumption of equilibrium chemistry is likely valid for the inner regions of the disk. Our model does not take into account deviations from equilibrium chemistry that may occur in the regions of the disk exterior to 4 AU. Because the planetesimals formed in this region are not considered in the dynamical simulations, these deviations are only important to the extent that they change the composition of the inward-moving material. Material is depleted no more than 15% by the time it reaches 4 AU. Therefore, deviations from equilibrium chemistry will only affect less than 15% of the material.

The final step of our planet formation model is to determine the number and masses of planets that can form in a system and from what radial range they accreted their material. This step is important because it allows us to track the compositional mixing of planetesimals from different regions of the disk. To do this, we tag a population of planetesimals with the final composition of planetesimals in the sequential condensation chemistry model and follow their growth into planets.

Rather than simulating the entire growth of small planetesimals all the way to planets, we employ the results of earlier works and begin at the end of the oligarchic growth phase. This allows us to skip a computationally expensive step that is not necessary for our purposes. Kokubo & Ida (2002) find that accretion through the oligarchic growth phase proceeds locally assuming that the radial migration of small fragments produced in planetesimal collisions does not significantly alter the surface density profile. Consequently, any planetesimal or embryo surviving at the end of this phase should have a composition similar to the original planetesimals at that location. On the other hand, during the final stage of planet formation, significant radial mixing occurs as planetary embryos scatter and collide with each other and slowly sweep up the remaining planetesimals in the system.

We begin our late-stage planet formation simulations at the end of the oligarchic growth phase. As described in Kokubo & Ida (2002), at this point, planetesimals have a bimodal distribution in mass. Planetary embryos are much more massive and make up about 50% of the mass, while the rest of the mass is in smaller planetesimals. The embryos are typically separated by about 10 mutual Hill radii. Following Chambers (2001) and O'Brien et al. (2006), we start our simulations with such a distribution of planetesimals placed throughout the disk such that the surface density profile of the disk obeys the relation

$$\begin{equation} \Sigma = \Sigma _0\left(\frac{r}{1\,{\rm AU}}\right)^{-3/2}, \end{equation} \tag{ 17 }$$

where Σ₀ = 8 g cm⁻² and then falls off linearly between 0.7 and 0.3 AU. We note that a fully self-consistent model would use the surface density profile calculated from the sequential condensation chemistry model rather than that of Equation (17). However, this would require separate N-body simulations for each different system. Because late-stage accretion is a stochastic process, it would then be very difficult to disentangle the effects of the different initial planetesimal compositions from differences in the N-body simulations. As described in Section 3.1, these two surface density profiles are compatible with each other. For these simulations, we use the N-body integrator, Mercury (Chambers 1999), which allows for a set of larger bodies that interact gravitationally with each other, as well as a set of smaller bodies that interact gravitationally with the larger bodies but not each other. This is ideally suited for the initial bimodal distribution of planetesimals used. Planets then grow through the collisions of these bodies. The collisions are assumed to be totally inelastic (i.e., they conserve momentum and form a single body with a mass of the combined mass of the colliding planetesimals).

We performed a set of four simulations, similar to those of O'Brien et al. (2006), with an initial population of 26 embryos of mass 0.09 M_⊕ and about 1000 planetesimals of mass 0.0024 M_⊕ distributed between 0.3 and 4 AU around a solar-mass star such that they obeyed the surface density relation given in Equation (17). The eccentricity and inclination of each planetesimal were given a random value of 0–0.01 and 0°–0.5°, respectively, and the longitude of the ascending node, argument of periapsis, and mean anomaly were assigned random values between 0° and 360°. Each system was integrated for 250 Myr with a time step of 5 days.

We chose four systems with stellar abundances that covered the range of Mg/Si and C/O ratios observed in stars. These ratios have the biggest influence on the mineralogy of the disk (Bond et al. 2010b), and so covering the full range should highlight the major differences between the equilibrium chemistry and sequential condensation models. The stellar abundances we used are the same as those used in Bond et al. (2010b), which were taken from Ecuvillon et al. (2004), Ecuvillon et al. (2006), Beirão et al. (2005), and Gilli et al. (2006), and the solar abundances were taken from Asplund et al. (2005). These abundances are shown in Table 2. We note that we only use the elemental compositions of these stars in our simulations. We do not use other physical parameters of the systems (e.g., stellar mass and planetary companions) for the N-body simulations. We choose to do this so that the same N-body simulations can be used for each system.

Before the sequential condensation model could be used to predict planetesimal composition, it was necessary to calibrate the planetesimal formation rate. Changing the constant of proportionality in the prescription for the planetesimal formation rate will increase or decrease the total mass in planetesimals that have formed by the end of the simulation and will also affect the chemical composition of the planetesimals. The planetesimal surface density that we model needs to match that which we believe existed during the early stages of planet formation. We can place constraints on this planetesimal surface density based on the masses of the terrestrial planets in the solar system (i.e., the minimum-mass solar nebula or MMSN; Weidenschilling 1977). Once planetesimals are formed, planet formation is a fairly efficient process, with planets accreting about 50% of the original mass (O'Brien et al. 2006). This means that although the MMSN is a lower limit, it is probably not far from the true amount of solid material that existed in the disk. The concept of the MMSN has been applied to extrasolar systems (e.g., Chiang & Laughlin 2013; Kuchner 2004) with a range of possible disk masses. In Section 3.2 we consider the effects of changing this normalization factor on our predicted planetesimal compositions.

In order to match the planetesimal surface density of the MMSN, we found that the planetesimal formation rate in our model should be

$$\begin{equation} \dot{\Sigma }_p = 9.5\times 10^{-6}\ {\rm yr}^{-1}\Sigma _{{\rm solid}}\frac{\Omega }{\Omega _{1\,{\rm AU}}}. \end{equation} \tag{ 18 }$$

This produced a distribution of planetesimal masses that is consistent with the r^−3/2 profile of the MMSN. N-body simulations of planet formation in the inner solar system (e.g., Chambers 2001; O'Brien et al. 2006) required surface density profiles that decreased inward of 0.7 AU in order to form low-mass Mercury analogs. However, the surface density profile produced in the sequential condensation model continues to increase interior to 0.7 AU. This discrepancy can be explained in a number of ways. For example, the prescription for the planetesimal formation/growth rate may be oversimplified and not appropriate for the innermost region of the disk, or planetesimals may become depleted in this region due to faster migration rates associated with higher gas densities. In our model, when we seed the disk with planetesimals for the dynamical simulations, we do not assume the surface density profile output by the sequential condensation model. Instead, we assume the surface density profile of Equation (17). The choice of surface density profile does not significantly affect the results because they are consistent with each other to begin with. However, we chose to use Equation (17) so that we could use the same set of N-body simulations for each of the different system compositions. This allows us to compare the effects of different initial compositions without having to account for the stochastic nature of late-stage planet formation.

Figure 1 shows the composition of planetesimals as a function of radius in the disk for each of the four different systems for both the equilibrium chemistry model and the sequential condensation model. In the equilibrium chemistry model, planetesimals were assumed to form instantaneously at 1.5 × 10⁵ yr. We chose a disk age of 1.5 × 10⁵ yr because it produces planets with compositions that matched those of the terrestrial planets of the solar system well. This figure highlights two major differences that occur by incorporating the time evolution of the disk into the chemistry model: (1) the radial range in which elements are present in planetesimals is extended in the sequential condensation model relative to the equilibrium chemistry model; and (2) the sequential condensation model leads to smoothly varying elemental abundances in planetesimals with increasing radius in the disk compared to the equilibrium chemistry model, where elements go from not present to fully condensed out over a small radial range.

The radial range where carbon makes up a significant fraction of the planetesimals' composition is much larger in the sequential condensation model, extending out to ∼3 AU, compared to the equilibrium chemistry model where carbon is only present in planetesimals interior to ∼1 AU. In carbon-rich systems, carbon is present in solid form at temperatures higher than 550 K. In the equilibrium chemistry model, carbon-rich planetesimals are limited to the region of the disk where these temperatures occur at the specific disk age used. On the other hand, planetesimals in the sequential condensation model are the accumulation of the solid material formed throughout the lifetime of the disk. This leads to carbon accumulating onto planetesimals farther out in the disk when it was younger and hotter, and closer to the star as the disk cools.

The relative amounts of the most abundant elements in systems with more solar-like chemistry are much less affected by the choice of chemistry model. Figure 1 shows that the ratios of Mg, Fe, and Si are relatively unchanged between models for these systems. This is due to the fact that these elements all condense out at about the same temperature (between 1300 K and 1330 K). The ratio of O to these three elements is also quite similar between the models interior to 2 AU. However, beyond 2 AU the sequential condensation model shows a much quicker increase in O abundance compared to the equilibrium chemistry model. This occurs because the ice line continues to move closer to the star as the disk cools, which is not captured in the equilibrium chemistry model. In addition to O, the relative abundances of some of the other more volatile elements can be seen to vary between the two models. Specifically, the amounts of S and Na show a gradual decrease closer to the star in the sequential condensation model compared to the abrupt change seen in the equilibrium chemistry model. Although this gradual decrease is qualitatively consistent with the volatile element fractionation patterns seen in solar system meteorites (see Palme et al. 1988 for a review), Ciesla (2008) found that this trend cannot be reproduced throughout the asteroid belt where it is observed in the solar system in sequential condensation-like models. This result suggests that other mechanisms, such as evaporation of volatiles (Davis & Richter 2005), may play an important role in determining the composition of planetesimals.

A more subtle effect of the sequential condensation model is that it depletes the disk of certain elements relative to others. As a parcel of gas and dust migrates through the disk, it deposits solid material in the form of planetesimals. The disk becomes depleted in the elements that composed this material, or in a relative sense, enriched in the elements that were not present. The exact amount of depletion seen in any given parcel will depend on where it originated from, the disk model that was used, and the planetesimal accretion rate.

This effect is extremely important when considering the C/O ratio of a disk. Carbon is abundant as a solid in the form of graphite at temperatures between ∼600 K and ∼1000 K, and silicon carbide at temperatures greater than ∼1000 K. It is also abundant in carbon ices below 78 K (Lodders 2003; we do not model these low-temperature ices because the disk does not get to low enough temperatures in the region and time of interest.) For a large portion of the disk (where temperatures are between 78 K and 600 K), very little carbon will condense out of the disk. In contrast, oxygen is the dominant component of rocky material (present at temperatures below 1300 K) and is even more abundant as a solid beyond the ice line. This means that throughout much of the disk, oxygen is being depleted much more than carbon, leading to an increase in the C/O ratio. In some cases the C/O ratio can reach high enough levels that once a parcel moves closer to the star and heats up, carbon will condense out of the disk, creating carbon-rich planetesimals. This can occur for systems that were not originally significantly enriched in C. Figure 2 shows the wt-% of carbon contained in planetesimals for disks with initial C/O ratios between solar and 1.00. The initial elemental abundances of the disk for each curve are all solar, except for the oxygen abundance, which has been decreased in order to obtain each C/O ratio. In contrast to the sequential condensation model, carbon-rich planetesimals can only be formed in the equilibrium chemistry model if the initial C/O ratio is greater than 0.8. Furthermore, the C/O enhancement of the sequential condensation model allows a much larger region of the disk to form carbon-rich planetesimals.

Zoom In Zoom Out Reset image size

The amount of carbon deposited in planetesimals will depend not only on the original C/O ratio of the disk but also on the disk model and the planetesimal accretion rate that is used. Exploring the possible range of disk models is beyond the scope of this paper, but their effect on protoplanetary disk chemistry merits further study. We increased the planetesimal accretion rate in the disk from the nominal value needed to achieve a disk mass comparable to the MMSN to three times that value. Figure 3 shows the effect of this increase on the carbon abundance of planetesimals in a disk with an original C/O ratio of 0.8. Higher accretion rates lead to faster depletion of oxygen and a correspondingly fast increase in the local C/O ratio. As a consequence, carbon condenses out of the disk gas at more distant radii, leading to carbon-enriched planetesimals throughout more of the disk.

Zoom In Zoom Out Reset image size

In each of the four N-body simulations, we form two or three planets within 2 AU (see Figure 4). They range in mass from 0.1 to 1.3 M_⊕ with a median mass of 0.7 M_⊕. These results are consistent with the Circular Jupiter and Saturn (CJS) simulations of O'Brien et al. (2006) off of which they were based. Figure 4 highlights both the importance and stochasticity of radial mixing of planetesimals in the final stages of terrestrial planet formation. The origin of material for any given planet can range from a small area around the planet's final location to practically the whole terrestrial-planet-forming zone. This stochasticity acts to increase the diversity of exoplanet compositions and also means that any one system is not representative of the others. The main importance of radial mixing is that it weakens any compositional gradients that existed at earlier stages.

Zoom In Zoom Out Reset image size

We first performed a benchmark test to see if our equilibrium chemistry simulations could reproduce the bulk compositions of the terrestrial planets of the solar system. As can be seen from Figure 5, the simulated compositions are in good agreement with the solar system values except for some of the more volatile elements, which are overabundant in the simulated planets. It is possible that non-equilibrium effects, such as evaporation during planetesimal collisions (Bond et al. 2010a), were important in the solar nebula. The elemental abundances for the solar system planets were obtained from Morgan & Anders (1980), Kargel (1993), and Lodders (1997) and are expected to have errors as high as 25% (Bond et al. 2010a). Our results are consistent with the findings of Bond et al. (2010a), which can also be seen in Figure 5. Small differences exist between our model and theirs, but these are expected due to the stochastic nature of planet formation and differences in the protoplanetary disk model used.

Zoom In Zoom Out Reset image size

Also shown in Figure 5 are the bulk elemental abundances of simulated planets using the sequential condensation chemistry model. In general, the results are very similar to those of the equilibrium chemistry. The only elements where we expected differences could occur based on the abundance curves in Figure 1 are O, Na, and S. These elements have intermediate condensation temperatures and therefore are not in the same proportions throughout the terrestrial planet formation zone. The differences seen in the abundance curves between the sequential condensation and equilibrium chemistry models were likely washed out due to the radial mixing of planetesimals during late-stage formation simulations, leading to the similarity in bulk compositions between the models.

In the simulations of carbon-rich systems there are significant differences between the equilibrium chemistry and sequential condensation models. Figure 6 shows the weight percent of the most abundant elements in the simulated planets for the HD 19994 system. The innermost two planets accreted most of their material from a region where carbon was solid in chemical equilibrium at 1.5 × 10⁵ yr and thus have a large fraction of their mass in carbon. This results in planets with similar bulk composition between the two models. However, the outermost planet of this system accreted little of its material from the carbon-rich region of the disk in the equilibrium chemistry model and consequently has a very small fraction of its mass in carbon. On the other hand, carbon is abundant in planetesimals beyond 2 AU in the sequential condensation model, thus forming a planet with a significant portion of carbon at 1.6 AU.

Zoom In Zoom Out Reset image size

This and previous works found that an equilibrium chemistry model of the solar nebula produces planets with remarkably similar abundances to the true terrestrial planet abundances. The sequential condensation model produces similar results but at the cost of a much more complex model. This begs the question, why use the sequential condensation model at all?

Firstly, the equilibrium chemistry model is not physically motivated in regard to its connection with planetesimal formation. The implicit assumption of the equilibrium chemistry model is that all planetesimals formed at the same time everywhere in the disk and on timescales much shorter than the timescales for disk evolution. For any given region of the disk, planetesimals should begin forming as soon as conditions are right. These conditions should occur at different times for different parts of the disk. Assuming that they all form simultaneously throughout the disk is unrealistic. Furthermore, we know from age dating of meteoritic material that planetesimal formation lasted for about 2.5 Myr (Amelin et al. 2002). The assumptions of the equilibrium chemistry model have been overlooked in the past because of the model's success in predicting the bulk composition of the solar system terrestrial planets.

The success of the equilibrium model can largely be attributed to the high condensation temperatures of most of the major rock-forming elements. For all but the hottest disks, these elements will be fully condensed out, and thus in solar proportions, throughout the entire terrestrial-planet-forming region. The refractory elements are of approximately solar proportions in the terrestrial planets (excluding Mercury), and so the models naturally match these. All that's really required then is a disk model that can reproduce the abundances of the more volatile elements (e.g., O, Na, S). Bond et al. (2010a) were able to find a model that produced planets with the correct O abundance, but these planets were overabundant in Na and S compared to the solar system planets. Nevertheless, their model accurately reproduced the abundances of 10 out of 12 of the most common elements in the terrestrial planets.

The situation is different for carbon-rich systems. In chemical equilibrium, carbon is abundant as a solid at temperatures higher than ∼550 K but is a gas at lower temperatures. For most disk ages, 550 K lies in the middle of the terrestrial-planet-forming region. Consequently, the predicted composition of planets will vary greatly in the equilibrium chemistry model depending on the choice of disk age. Because there are no carbon-rich systems in which we know the composition of the planets, it is not obvious which, if any, disk age will produce the correct results.

We suggest that the equilibrium chemistry model's ability to reproduce the abundances in the solar system is a natural consequence of the condensation properties of a disk of solar composition, but that it does not necessarily produce realistic results for carbon-rich systems. The sequential condensation model provides a more realistic picture of planetesimal formation in which planetesimals can form in any region of the disk at any time as long as the conditions are correct. This is consistent with the idea that planetesimals form over a few million years, and the model produces planets with elemental abundances in good agreement with those of the solar system's planets. Although the results for both models are similar in the case of the solar system, the sequential condensation model likely predicts more realistic compositions for planets in carbon-rich systems.

The most uncertain part of the sequential condensation model is the prescription for planetesimal formation and growth. This is necessarily so due to the current uncertainties in planetesimal formation theory. The prescription used in this work is likely oversimplified.

We showed that simply increasing the magnitude of the planetesimal accretion rate has a significant effect on the chemistry. Much larger effects would likely occur by changing its functional form. For example, the planetesimal formation rate could be a stronger function of the density of solids in the disk. In such a model, comparatively more material would form planetesimals earlier in the disk's lifetime, when it is denser and hotter, causing the planetesimal disk to be more enhanced in refractory elements and depleted in volatile elements. Such a scenario might be able to reduce the overabundance of the more volatile elements we see in our current solar system model. Applied to moderately carbon-rich disks, this could lead to enough depletion of oxygen early on that solids containing carbon could condense farther away from the star, thus increasing the amount of solid carbon in the terrestrial planet formation zone.

In the sequential condensation model, planetesimals form at a given distance from the star and remain there for the duration of the simulation. However, planetesimals are expected to migrate radially due to gas drag and should accrete material from different parts of the disk. Including planetesimal migration is beyond the scope of this work as it would require knowledge of the distribution of planetesimal sizes and masses. However, migration may significantly affect the composition of forming planetesimals and merits future study.

The sequestration of oxygen in planetesimals has been used as a mechanism to explain the high observed C/O ratios in some gas giant exoplanets (Öberg et al. 2011; Helling et al. 2014). Additionally, the sequestration of oxygen in the outer disk and resulting depletion in the inner disk have been studied in the context of the distribution of water in the solar nebula (Ciesla & Cuzzi 2006), but their influence on the abundances of other chemical species was not considered. Najita et al. (2011) pointed out that this depletion may influence the abundances of hydrocarbons in the inner disk and that these differences should be observable. Carr & Najita (2011) and Najita et al. (2011) find a correlation between the HCN/H₂O flux ratio and disk masses for T Tauri stars. They suggest that this trend could be explained if a high HCN/H₂O flux results from a high C/O ratio and that the C/O ratio is enhanced more in larger disks due to more efficient planetesimal formation. Similarly, Pascucci et al. (2013) find higher average HCN/H₂O fluxes in brown dwarf disks than T Tauri stars and attribute this to more efficient planetesimal formation in disks around brown dwarfs. They go on to point out that the C/O ratios needed to explain their observations would have important implications for the compositions of rocky planets forming in these carbon-rich regions of the disk.

Our models explore exactly this point. Indeed, we find that a portion of the disk can contain a substantial amount of solid carbon for disks with an initial C/O ratio as low as 0.65. This finding is unique compared to the equilibrium chemistry model, in which carbon-rich planetesimals can only form in disks with initial C/O ratios above 0.8.

The carbon-rich region in the sequential condensation model is generally confined to less than ∼0.5–0.6 AU for disks with C/O ratios between 0.65 and 0.8. As can be seen from Figure 4, this region accounts for a small fraction of a planet's material source region, creating planets with at most a few percent of their mass in carbon. The dynamical simulations in this work are limited to solar-system-like initial conditions (i.e., the planetesimal disk extends from 0.3 to 4 AU). However, it is quite possible that planets in other systems could form much closer to their host stars (e.g., Bond et al. 2010b; Raymond et al. 2008) and thus potentially accrete more carbon-rich material. Chiang & Laughlin (2013) discuss the possibility that most of the known close-in exoplanets formed in situ rather than migrating there. If this is the case or even partially the case, then we predict that carbon-rich planets may be significantly more common than previously indicated.

It is also possible that a further exploration of parameter space in our model will change the radial range where carbon-rich planetesimals can exist. There are a number of factors that could influence where and when oxygen-rich material is removed from the disk and turned into planetesimals, including the following. (1) The magnitude and functional form of the planetesimal accretion rate. We showed that by increasing the overall planetesimal accretion rate, the carbon-rich region of the disk moved inward. Additionally, the rate at which planetesimals in different parts of the disk form and grow may not be governed by Equation (16). Variation from this form would completely change the pattern of oxygen depletion and thus the C/O enhancement. (2) The disk temperature and density structure as a function of time. The structure and evolution of the disk are undoubtedly connected to the formation and growth of planetesimals. In our planetesimal formation prescription this is true because the planetesimal formation rate is dependent on the surface density of the disk. By changing the equations of disk evolution, we would also be changing the timing of planetesimal formation and thus their composition as well. (3) The overall amount of (C+O)/(Mg+Si). For a fixed C/O ratio, if oxygen (and carbon) are depleted relative to their solar values but Mg and Si are not, then the relative fraction of O that bonds with Mg and Si is higher, and thus the effective O depletion is larger and consequently the C/O enhancement is stronger.