Chemistry in an Evolving Protoplanetary Disk

Chemistry in an Evolving Protoplanetary Disk: Effects on Terrestrial Planet Composition


The composition of planets is largely determined by the chemical and dynamical evolution of the disk during planetesimal formation and growth. To predict the diversity of exoplanet compositions, previous works modeled planetesimal composition as the equilibrium chemical composition of a protoplanetary disk at a single time. However, planetesimals form over an extended period of time, during which, elements sequentially condense out of the gas as the disk cools and are accreted onto planetesimals. To account for the evolution of the disk during planetesimal formation, we couple models of disk chemistry and dynamics with a prescription for planetesimal formation. We then follow the growth of these planetesimals into terrestrial planets with N-body simulations of late stage planet formation to evaluate the effect of sequential condensation on the bulk composition of planets. We find that our model produces results similar to those of earlier models for disks with C/O ratios close to the solar value (0.54). However, in disks with C/O ratios greater than 0.8, carbon rich planetesimals form throughout a much larger radial range of the disk. Furthermore, our model produces carbon rich planetesimals in disks with C/O ratios as low as 0.65, which is not possible in the static equilibrium chemistry case. These results suggest that (1) there may be a large population of short period carbon rich planets around moderately carbon enhanced stars (0.65 C/O 0.8) and (2) carbon rich planets can form throughout the terrestrial planet region around carbon rich stars (C/O 0.8).

Subject headings:
planets and satellites: composition, planets and satellites: formation, planets and satellites: terrestrial planets, protoplanetary disks

1. Introduction

With the detection of over 700 confirmed exoplanets and over 3,000 exoplanet candidates (; Wright et al., 2011) in a wide variety of planetary system architectures, it has become abundantly clear that there are a diversity of outcomes to the planet formation process. Furthermore, planet mass and radius data coupled with interior structure models, as well as transmission and emission spectroscopy, indicate that exoplanets may be as diverse in their bulk compositions as they are in their orbital properties. Of the small but growing sample of exoplanets for which mass and radius and/or spectroscopic data exist, a range of bulk compositions have already been inferred, including: Mercury-like compositions (e.g. CoRoT-7b and Kepler-10b) with iron cores composing as much as 65 wt-% of the planet (Wagner et al., 2012) to water-worlds (e.g. GJ1214b; Berta et al., 2012). There are even indications of planets with carbon rich atmospheres (e.g. WASP-12b; Madhusudhan et al., 2010) and interiors (55 Cancri e; Madhusudhan et al., 2012).

The potential diversity of terrestrial exoplanet compositions has been addressed from a theoretical standpoint in several publications (Kuchner & Seager, 2005; Bond et al., 2010b; Carter-Bond et al., 2012a, b). The most recent three of these works sought to predict the range of exoplanet compositions that should exist based on the range of elemental abundances that have been observed in planet hosting stars. It was assumed that the stellar abundances reflect those of the initial protoplanetary disk and consequently there should be a connection between the composition of planets that formed from this disk and their host star. The simulations consisted of models of protoplanetary disk chemistry and late stage planet formation, which were coupled together by tagging the initial planetesimals of the planet formation simulation with the composition of the disk’s solid material at their initial location. The planets formed in these simulations had compositions ranging from Earth-like to almost entirely C and SiC. These carbon rich planets form in disks with C/O 0.8, where carbon readily condenses to form solids (Larimer, 1975).

The composition of the disk’s solid material in these simulations was determined with equilibrium chemistry calculations assuming the mid-plane pressure and temperature profile of the disk. Because the mid-plane pressure and temperature profiles of a disk change as the disk ages, it is not obvious which disk age to use when calculating equilibrium abundances. In their work, they adopted an age from a prior study (Bond et al., 2010a) that best reproduced the composition of the Solar system’s terrestrial planets. Elser et al. (2012) attempted a more self-consistent approach to determining the disk age for the equilibrium chemistry calculation by imposing a transition condition: the surface density of solid material predicted by the disk model must match the initial surface density of the dynamical simulations. With this condition they had difficulty reproducing the abundances of the Solar system terrestrial planets and found that their results varied depending on the disk model used.

In both of these approaches, the composition of the solids in the disk are determined assuming they all form at the same time. That is, equilibrium composition is calculated at a single age in the disk’s lifetime and all solids are assigned that composition. Age dating of meteoritic material suggests that solid material condensed out of the Solar Nebula over the course of about 2.5 Myr or more (Amelin et al., 2002). During this time, the Solar Nebula would have cooled significantly resulting in solids being formed over a range of temperatures.

As a protoplanetary disk cools, it will sequentially condense out each of the elements starting with the most refractory and progressing towards the most volatile. At the same time, the disk is also losing mass. Elements that condense at later times/cooler temperatures will be less abundant and therefore depleted compared to the more refractory elements. Cassen (1996, 2001) showed that taking into account the cooling of the Solar Nebula during planetesimal formation led to moderately volatile element ratios that were consistent with meteoritic abundances. Ciesla (2008), re-examined this result, this time incorporating the inward migration of small planetesimals, and found that the moderately volatile depletion patterns in the asteroid belt could only be produced in models with a narrow range of parameters that were inconsistent with the  3 million year timescale of planetesimal formation. They did, however, find that such depletion patterns could occur closer to the Sun, potentially affecting the composition of some of the terrestrial planets.

In this paper, we take an approach similar to Cassen (1996) and couple models of protoplanetary disk evolution, equilibrium chemistry and planetesimal growth in order to predict the chemical composition of planetesimals and the planets that they form. We apply this method to systems with solar-like compositions and systems with high C/O ratios and show that considering the growth of planetesimals in an evolving disk affects the chemical composition of the resulting planets. Lastly, we consider the effect that the depletion of oxygen in the outer disk has on planetesimal compositions in the inner disk.

2. Methods

We simulate the formation of terrestrial planets with a coupled chemical and dynamical model. A theoretical model of the protoplanetary disk provides the radial temperature, pressure and density structure of the disk as a function of time. We assume that the elemental abundances of the host star reflect the original elemental abundances of the disk. Combining the disk model and the stellar abundances, we calculate the chemical equilibrium composition of the disk as a function of radius. This provides a chemical model of the disk over time which is combined with a prescription for the planetesimal formation rate to determine the composition of planetesimals that form in the disk. We then follow the dynamical evolution of a population of planetesimals as they collide and build planets. The final composition of these planets is the combination of each of the planetesimals that they accreted.

2.1. Disk Model

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




and and are the initial mass and outer edge of the disk. The temperature of the disk in this region is given by




is the opacity and is taken to be a constant 3 , is the viscosity parameter, is the mean molecular weight, is the adiabatic index, and is the mass of hydrogen. The viscous timescale is:


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),


where is the temperature at which the majority of solids evaporate, is the opacity at and is the power law index taken to be -14. The surface density in the inner region of the disk is,




The temperature in the inner region is given by


The transition radius between these regions is


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






The temperature in the outer region is


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


For all calculations in this work we use , , and . This model corresponds to the second example in Stepinski (1998) which is consistent with a planetesimal forming disk.

2.2. Equilibrium Chemistry

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 reflect 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 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, 2007, and references therein).

CO Fe Mg Ni SN TiO


AlN Fe MgSiO(OH)
Ca(PO) FeC NaAlSiO
CaAlSiO FeO Ni
CaAlO FeS Si
Cr MgS TiN

Table 1Chemical species used in equilibrium calculations.

2.3. Sequential Condensation Chemistry

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 therin) 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:


where 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 effect 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 effect less than 15% of the material.

2.4. Dynamical Simulations

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 which 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 bimodel 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


where 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 bimodel 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 and about 1000 planetesimals of mass 0.0024 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 was given a random value between and 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.

2.5. System Compositions

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.

Element Solar System 55 Cnc HD19994 HD 213240
Al 2.34 8.7 6.2 4.7
C 2.45 7.4 8.9 5.4
Ca 2.04 2.8 3.4 2.5
Cr 4.37 7.8 7.4 5.2
Fe 2.82 6.3 5.1 4.4
H 1.00 1.0 1.0 1.0
He 8.51 8.5 8.5 8.5
Mg 3.39 1.2 6.2 6.5
N 6.03 1.8 2.2 1.3
Na 1.48 3.9 6.5 3.6
Ni 1.70 3.6 3.3 2.4
O 4.57 7.4 7.1 1.2
P 2.29 8.5 6.0 4.6
S 1.38 2.1 1.4 1.3
Si 3.24 6.9 6.0 4.4
Ti 7.94 2.2 1.5 1.3

Table 2Elemental abundances for systems considered.

3. Results

3.1. Planetesimal Surface Density

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 effect 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:


This produced a distribution of planetesimal masses that is consistent with the 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 inwards 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 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 affects of different initial compositions without having to account for the stochastic nature of late stage planet formation.

3.2. Disk Chemistry

Figure 1.— Elemental abundances as a function of radius () in the disk for the four systems considered. The left column shows results from the sequential condensation model and the right column shows results from the equilibrium chemistry model.

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 years. We chose a disk age of years 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. 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 effected 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 1300K and 1330K). 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 (?), 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 600K and 1000K, and silicon carbide at temperatures greater than 1000K. It is also abundant in carbon ices below 78K (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 78K and 600K), very little carbon will condense out of the disk. In contrast, oxygen is the dominant component of rocky material (present at temperatures below 1300K) 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.

Figure 2.— The wt-% of carbon contained in planetesimals as a function of radius () in the disk. Each line represents a disk with a different initial C/O ratio ranging from 0.54 (solar) to 1. For comparison, the equilibrium chemistry results (dashed lines) are shown in addition to the sequential condensation results (solid lines). The sequential condensation model produces carbon rich planetesimals throughout a larger radial range of the disk and in disks with lower initial C/O ratios than the equilibrium chemistry model does.

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 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.

Figure 3.— The wt-% of carbon contained in planetesimals at different radii () in the disk for different planetesimal accretion rates. The planetesimal accretion rate used for each simulation is the nominal value needed to obtain a final planetesimal surface density similar to the MMSN scaled by the value in the legend.

3.3. N-Body Simulations

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-1.3M with a median mass of 0.7M. These results are consistent with the CJS simulations of O’Brien et al. (2006) from which they were based off of. 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.

Figure 4.— Material source regions for simulated planets. Each pie plot corresponds to a simulated planet and each row to one simulation. Each planet’s final orbital position is indicated by its locations along the x-axis and its mass is proportional to the radius of the pie plot cubed. Each colored slice represents material that originated from the region of the disk indicated in the color bar.

3.4. Planet Bulk Composition

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 over-abundant 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.

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.

Figure 5.— Elemental abundances of simulated planets normalized to abundances of corresponding Solar System planets. Results from our sequential condensation (SC) model and our equilibrium chemistry (EC) model are shown in addition to those of Bond et al. (2010a)

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 years 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.6AU.

Figure 6.— Composition of simulated planets in the HD19994 system. Each pie chart shows the fraction of a planet’s mass in each of the most abundant elements and its semi-major axis (along the x-axis). Results for the equilibrium chemistry model are shown in the top row and those for the sequential condensation model are shown in the bottom row.

4. Discussion

4.1. Why Should the Sequential Condensation Model be Used?

This and previous works both 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 regards 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 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, etc). Bond et al. (2010a) were able to find a model that produced planets with the correct O abundance, but these planets were over abundant in Na and S compared to the Solar System planets. Nevertheless, their model accurately reproduced the abundances of ten out of twelve 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.

4.2. Formation and Growth of Planetesimals

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 over abundance 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.

4.3. Migration of Planetesimals

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 a knowledge of the distribution of planetesimal sizes and masses. However, it is possible that migration may significantly effect the composition of forming planetesimals and merits future study.

4.4. C/O Enhancements

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 has been studied in the context of the distribution of water in the Solar Nebula (Ciesla & Cuzzi, 2006), but its 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/HO flux ratio and disk masses for T Tauri stars. They suggest that this trend could be explained if a high HCN/HO 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/HO 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-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: (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 inwards. 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 is undoubtedly connected to the formation and growth of planetesimals. In our planetesimal formation prescription this is true because the planetesimal formation rate is a 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.

5. Summary and Conclusions

We have implemented a model of protoplanetary disk chemistry that accounts for the time evolution of the disk and the formation and accretion of planetesimals. Previous models performed equilibrium chemistry calculations at a single point in the disk’s lifetime. We have compared the results of the sequential condensation model to those of the equilibrium chemistry model and find three main results:

  • For a disk of solar composition, the sequential condensation model produces planets with elemental abundances very similar to those produced by the equilibrium chemistry model. Of the elements considered, the abundances match the true values of the Solar System terrestrial planets remarkably well with the exception of sodium and sulfur. Whereas the equilibrium chemistry model was already optimized to reproduce the compositions of the Solar System planets, the full parameter space of the sequential condensation model has yet to be explored. Tuning the sequential condensation model should produce better fits to the abundances of the more volatile elements in the Solar System planets.

  • In carbon rich disks, the sequential condensation model predicts planetesimals rich in carbon over a wider range of orbital radii than the equilibrium chemistry model. When compared to the Solar System optimized case, this results in carbon enriched planets at larger semi-major axes. Overall, the sequential condensation model produces planetesimals with compositions that vary less over the radial range of the terrestrial planet forming region.

  • The sequential condensation model predicts that carbon enriched planetesimals can form in initially oxygen dominated protoplanetary disks. These planetesimals are only produced in the innermost regions of the disk using the current model parameters. However, it is possible that an exploration of parameter space may find that the sequential condensation model can produce carbon rich planetesimals at greater distances from the star. Additionally, if a significant fraction of known close in super-Earths formed in situ, the number of carbon enriched exoplanets may be significantly more than previously indicated.

We thank the referee, John Chambers, for helpful comments used to improved this paper. NM acknowledges support from the Yale Center for Astronomy and Astrophysics (YCAA) at Yale University through the YCAA postdoctoral prize fellowship.


  1. Amelin, Y., Krot, A. N., Hutcheon, I. D., & Ulyanov, A. A. 2002, Science, 297, 1678, 1678
  2. Asplund, M., Grevesse, N., & Sauval, A. J. 2005, in ASP Conf. Ser. 336, Cosmic Abundances as Records of Stellar Evolution and Nucleosynthesis, ed. T. G. Barnes & F. N. Bash, 25
  3. Beirão, P., Santos, N. C., Israelian, G., & Mayor, M. 2005, A&A, 438, 251, 251
  4. Berta, Z. K., Charbonneau, D., Désert, J.-M., et al. 2012, ApJ, 747, 35, 35
  5. Blum, J., & Wurm, G. 2008, ARA&A, 46, 21, 21
  6. Bond, J. C., Lauretta, D. S., & O’Brien, D. P. 2010a, Icarus, 205, 321, 321
  7. Bond, J. C., O’Brien, D. P., & Lauretta, D. S. 2010b, ApJ, 715, 1050, 1050
  8. Carr, J. S., & Najita, J. R. 2011, ApJ, 733, 102, 102
  9. Carter-Bond, J. C., O’Brien, D. P., Delgado Mena, E., et al. 2012a, ApJL, 747, L2, L2
  10. Carter-Bond, J. C., O’Brien, D. P., & Raymond, S. N. 2012b, ApJ, 760, 44, 44
  11. Cassen, P. 1996, M&PS, 31, 793, 793
  12. —. 2001, M&PS, 36, 671, 671
  13. Chambers, J. E. 1999, MNRAS, 304, 793, 793
  14. —. 2001, Icarus, 152, 205, 205
  15. —. 2009, ApJ, 705, 1206, 1206
  16. Chiang, E., & Laughlin, G. 2013, MNRAS, 431, 3444, 3444
  17. Ciesla, F. 2008, M&PS, 43, 639, 639
  18. Ciesla, F., & Cuzzi, J. 2006, Icarus, 181, 178, 178
  19. Cuzzi, J. N., Hogan, R. C., & Shariff, K. 2008, ApJ, 687, 1432, 1432
  20. Davis, A. M., & Richter, F. M. 2007, in Treatise on Geochemistry, ed. H. D. Holland & K. K. Turekian (Oxford: Pergamon), 1–31
  21. Ecuvillon, A., Israelian, G., Santos, N. C., et al. 2004, A&A, 426, 619, 619
  22. —. 2006, A&A, 445, 633, 633
  23. Elser, S., Meyer, M. R., & Moore, B. 2012, Icarus, 221, 859, 859
  24. Gilli, G., Israelian, G., Ecuvillon, A., Santos, N. C., & Mayor, M. 2006, A&A, 449, 723, 723
  25. Helling, C., Woitke, P., Rimmer, P., et al. 2014, Life, 4, 142, 142
  26. Johansen, A., Oishi, J. S., Low, M.-M. M., et al. 2007, Nature, 448, 1022, 1022
  27. Kargel, J. 1993, Icarus, 105, 1, 1
  28. Kokubo, E., & Ida, S. 2002, ApJ, 581, 666, 666
  29. Kuchner, M. J. 2004, ApJ, 612, 1147, 1147
  30. Kuchner, M. J., & Seager, S. 2005, arXiv, 0504214
  31. Larimer, J. W. 1975, Geochimica et Cosmochimica Acta, 39, 389, 389
  32. Lodders, K. 1997, Icarus, 126, 373, 373
  33. —. 2003, ApJ, 591, 1220, 1220
  34. Madhusudhan, N., Lee, K. K. M., & Mousis, O. 2012, ApJ, 759, L40, L40
  35. Madhusudhan, N., Harrington, J., Stevenson, K. B., et al. 2010, Nature, 469, 64, 64
  36. Morgan, J. W., & Anders, E. 1980, PNAS, 77, 6973, 6973
  37. Najita, J. R., Ádámkovics, M., & Glassgold, A. E. 2011, ApJ, 743, 147, 147
  38. Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, Nature, 743, L16, L16
  39. O’Brien, D. P., Morbidelli, A., & Levison, H. F. 2006, Icarus, 184, 39, 39
  40. Palme, H., Larimer, J. W., & Lipschutz, M. E. 1988, in Meteorites and the Early Solar System, ed. J. F. Kerridge & M. S. Matthew (Tucson, AZ: Univ. Arizona Press), 436–461
  41. Pascucci, I., Herczeg, G., Carr, J. S., & Bruderer, S. 2013, ApJ, 779, 178, 178
  42. Pizzarello, S., Cooper, G. W., & Flynn, G. J. 2006, Meteorites and the Early Solar System II, 625, 625
  43. Raymond, S. N., Barnes, R., & Mandell, A. M. 2008, MNRAS, 384, 663, 663
  44. Stepinski, T. 1998, Icarus, 132, 100, 100
  45. Wagner, F. W., Tosi, N., Sohl, F., Rauer, H., & Spohn, T. 2012, A&A, 541, A103, A103
  46. Weidenschilling, S. J. 1977, Ap&SS, 51, 153, 153
  47. White, W. B., Johnson, S. M., & Dantzig, G. B. 1958, JChPh, 28, 751, 751
  48. Wright, J. T., Fakhouri, O., Marcy, G. W., et al. 2011, PASP, 123, 412, 412
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description