Dust condensation in evolving discs and the composition of planetary building blocks

ABSTRACT

Partial condensation of dust from the Solar nebula is likely responsible for the diverse chemical compositions of chondrites and rocky planets/planetesimals in the inner Solar system. We present a forward physical–chemical model of a protoplanetary disc to predict the chemical compositions of planetary building blocks that may form from such a disc. Our model includes the physical evolution of the disc and the condensation, partial advection, and decoupling of the dust within it. The chemical composition of the condensate changes with time and radius. We compare the results of two dust condensation models: one where an element condenses when the mid-plane temperature in the disc is lower than the 50 per cent condensation temperature (⁠|$\rm T_{50}$|⁠) of that element and the other where the condensation of the dust is calculated by a Gibbs free energy minimization technique assuming chemical equilibrium at local disc temperature and pressure. The results of two models are generally consistent with some systematic differences of ∼10 per cent depending upon the radial distance and an element’s condensation temperature. Both models predict compositions similar to CM, CO, and CV chondrites provided that the decoupling time-scale of the dust is of the order of the evolution time-scale of the disc or longer. If the decoupling time-scale is too short, the composition deviates significantly from the measured values. These models may contribute to our understanding of the chemical compositions of chondrites, and ultimately the terrestrial planets in the Solar system, and may constrain the potential chemical compositions of rocky exoplanets.

accretion, accretion discs, astrochemistry, solid state: refractory, solid state: volatile, stars: pre-main-sequence

1 INTRODUCTION

Our Solar system is an important testbed for studying the process of planet formation. The chemical composition of terrestrial planets and both differentiated and undifferentiated meteorites can constrain its formation and evolution. All terrestrial planets and many carbonaceous chondrites (CM, CO, and CV) are depleted in volatile elements and enriched in refractory elements relative to the composition of the Sun, as represented by CI chondrites (Asplund, Grevesse & Sauval 2005; Asplund et al. 2009). The degree of depletion of these elements is correlated with their volatility (Palme, Lodders & Jones 2014).

Fig. 1 shows the relative elemental abundances of terrestrial planets and a number of carbonaceous chondrites normalized to the Solar abundance and Si

$$\begin{eqnarray*} {\rm R_i=\frac{[\chi _{E}/\chi _{Si}]_P}{[\chi _{E}/\chi _{Si}]_{\odot }}}, \end{eqnarray*}$$

(1)

as a function of the ‘50 per cent condensation temperatures’, |$\rm T_{50}$|⁠, for each element as tabulated in Lodders (2003)(See also Table 1 here). The 50 per cent condensation temperature is the temperature where half of a given element is in the condensed phase under a total pressure of 10⁻⁴ bar with Solar composition. Here, χ_E and χ_Si stand for the concentrations of an element (E) and Si, respectively, and subscript P refers to a planet or chondrite. |$\rm R_i$| for terrestrial planets and CM, CO, and CV chondrites are smaller than unity for volatile elements while larger than unity for refractory elements. Determining causes of volatile element depletion in rocky planets, as well as the role that the physical and chemical environment plays in producing this depletion, is crucial for our understanding of the origin of our Solar system and is the subject of ongoing research in planetary science.

Figure 1.

Observed elemental abundances: Relative abundances, |$\rm R_i$|⁠, in terrestrial planets and CM, CO, and CV chondrites, normalized to the CI chondrites and Si, as a function of |$\rm T_{50}$| (Lodders 2003). Data sources: chondrites (Wasson & Kallemeyn 1988), bulk silicate Mercury (Morgan & Anders 1980), bulk silicate Venus (Morgan & Anders 1980), bulk silicate Earth (McDonough & Sun 1995), and bulk silicate Mars (Lodders & Fegley 1997). The horizontal line at 1 indicates where a body is neither enriched nor depleted in a given element relative to Si. The horizontal dashed lines show the median depletion in these objects at |$\rm T_{50}$| above and below 1310 K. This figure shows that all terrestrial planets and chondrites are depleted in volatiles, and are consequently enriched in refractory elements relative to the Sun. Highly siderophile elements in Earth are not plotted because they were affected by core formation. Elements in red font are modelled with the grains code (See Section 2.4).

Table 1.

The 50 per cent condensation temperatures from different models.

Element	50 per cent condensation temperature			grains T|$_2^{d}$|
	grainsa	Loddersb	Cassenc
H	...	...	...	...
He	...	...	...	...
C	...	40	...	...
N	...	123	...	...
O	...	180	...	...
Na	988	958	971	1069
Mg	1336	1336	...	1385
Al	1654	1653	...	1675
Si	1316	1310	1311	1478
P	1318	1229	1151	1358
S	669	664	652	711
Cl	440	948	...	446
K	906	1006	1002	1014
Ca	1500	1517	...	1669
Ti	1583	1582	...	1594
Cr	1304	1296	1276	1358
Mn	1154	1158	1191	1263
Fe	1335	1334	...	1360
Co	1350	1352	...	1361
Ni	1264	1353	1354	1326
Cu	1100	1037	1038	1232
Ga	1018	968	920	1159
Ge	911	883	826	1062
Mo	1548	1590	...	1644
Ru	1544	1551	...	1630
Pd	1336	1324	1334	1360
Hf	1730	1684	...	1749
W	1788	1789	...	1832
Re	1770	1821	...	1830
Os	1809	1812	...	1833
Ir	1584	1603	...	1691
Pt	1396	1408	...	1509
Au	1196	1060	1224	1322

Notes. ^a|$\rm T_{50}$|s calculated from the grains code. If |$\rm T_{50}$| for an element is less than 300 K, we do not provide it here.

^b|$\rm T_{50}$|s showed in Lodders (2003).

^c|$\rm T_{50}$|s used in Cassen (1996).

^dThe condensation temperatures calculated from the grains code. When the amount of an element in the condensed phase is larger than 2 per cent of the total amount of the element, we define the corresponding temperature as the condensation temperature. Here the pressure is 10⁻⁴ bar.

Since planetesimals and planets form in protoplanetary discs, their composition should be a function of both the evolution of the gaseous disc and the condensation of the elements in it. Since the early 1970s, chondrites and chondritic components such as CAIs (Calcium-Aluminium-rich Inclusions) are hypothesized to form from the hot solar nebula during the early stage of the Solar system formation. To understand how this formation proceeds, one must model both the physical and chemical processes involved.

1.1 Previous work

Partial condensation models are favoured to explain the volatile depletion in rocky planets and many carbonaceous chondrites (CM, CO, and CV). Grossman (1972) was among the first to use a thermodynamic approach to study the chemical evolution of the Solar nebula and the condensation of dust. He suggested that the CAIs are direct condensates from the nebula. Other early calculations were performed by Urey (1955), Lord (1965), and Larimer (1967). Lewis (1974) applied the same concept to understand the chemical composition of the terrestrial planets.

An important milestone for our purposes is Cassen (1996), who modelled the chemical evolution of an evolving disc to explain the volatile element depletion in CM, CO, and CV chondrites. He assumed that an element begins to condense once the temperature at a given radius in the disc drops below |$\rm T_{50}$| (estimated at a pressure of 10⁻⁴ bar). Despite this simplifying assumption, he was able to broadly reproduce the volatile element depletions in carbonaceous chondrites by forming them between 1 and 3 AU. Cassen (2001) used a similar model to calculate the relative abundance of the elements in the large rocky materials at different radii. The main differences in this case were the inclusion of 2D temperature distributions within the disc, external illumination, and the migration of the coagulated solids towards the Sun.

Later, Ciesla (2008) investigated whether the cooling disc can explain the depletion of moderately volatile elements (with condensation temperature between ∼650 and 1350 K) in chondritic meteorites. Compared to Cassen (1996), he calculated the viscous evolution of the disc and included the effects of gas drag on migrating pebbles. He found that the observed depletion can be explained by his model within a small range of parameters. In order to reproduce the depletion pattern observed in the asteroid belt, the initial temperature of the disc needed to be higher than roughly 1350 K beyond 2 AU and the formation time-scale of planetesimals should be less than 10⁵ yr. From this result it remains challenging to explain the observations that chondrule formation began a few million years after CAI formation (Kita et al. 2000; Amelin et al. 2002; Connelly et al. 2012).

Other groups considered the chemical equilibrium of the system. Using the commercial software package hsc chemistry, Bond, Lauretta & O’Brien (2010a) calculated the equilibrium partitioning of 16 elements among 78 gaseous and 33 solid species at seven specific times in the evolving disc. They then modelled the formation of the terrestrial planets using N-body accretion simulations (O’Brien, Morbidelli & Levison 2006). Using the same method but different initial chemical compositions of host stars, Bond, O’Brien & Lauretta (2010b) also modelled the formation of extrasolar planets – showing that the resulting compositions of exoplanets are diverse, depending upon elemental abundances of the host star, ranging from Earth-like to carbide-rich planets.

Using the hsc chemistry software, three different disc models, and elemental abundances of the present-day Solar photosphere (Asplund et al. 2005), Elser, Meyer & Moore (2012) modelled the composition of the terrestrial planets. They adopted a self-consistent approach to determine the time in the evolving disc for calculating radial variations in disc chemistry assuming chemical equilibrium. They showed that disc models with high temperatures and pressures produce a large variety in the composition of the resulting planets. Pignatale et al. (2016) adopted a 2D disc model and used the same software to model 170 gas and 317 solid compounds that form from 15 elements. They used a snapshot of the disc at 1 Myr as a proxy for typical disc conditions (temperature and pressure in the mid-plane of the disc) and modelled chemical and mineralogical gradients in the disc.

Given this background, none of these previous simulations consider simultaneous physical and chemical disc evolution or how differential transport of dust and gas affect dust condensation calculation. As the disc evolves with time, the surface density, temperature, and pressure in the disc change, which affects the condensation of chemical elements. Therefore, a time-dependent disc model is needed which considers the effects of changing pressure and surface density on the condensation temperatures as a function of time and radial distance of the evolving disc.

In this work, we present composition results for a dynamically evolving disc that incorporates real-time dust condensation. We use two different models. The first is the same method as Cassen (1996) but with updated T₅₀s. The second is a modified version of the evolving disc from Cassen (1996) and to use grains code from Petaev (2009) to model the chemical equilibrium.

Previous work by Moriarty, Madhusudhan & Fischer (2014) did include a model that took into account both the evolution of disc and the disc chemistry, using the same chemical equilibrium calculation as Bond et al. (2010a, b) and modelling the disc and dust evolution in a fashion similar to our approach. Their motivation was to seed the initial composition of a planetesimal disc to model terrestrial planet formation in various extrasolar systems with different chemical compositions, rather than to compare with Solar system chondrites. They found a good agreement with previous studies for Sun-like stars with C/O ∼ 0.54 and predicted a number of carbon-rich planets around stars with C/O > 0.65. We build upon their approach by including a larger list of chemical elements and compounds (33 versus 16 elements and 762 versus 111 species, respectively). In addition, we examine the effects of changing the assumptions on the disc evolution and dust condensation time-scales.

2 MODEL AND METHODS

2.1 Algorithm overview

In our calculations, the elements are partitioned among three phases: one gaseous and two condensed (advected dust and decoupled dust) phases. The partitioning of elements is based upon the temperature of the disc mid-plane – and is either found directly from the full chemical equilibrium calculation, or estimated based on |$\rm T_{50}$|⁠. The dust phase decouples from the gas on a time-scale proportional to the local orbital period (see equations 10 and 11). The portion of the dust that decouples from the gas remains at the location where it condensed and no longer chemically interacts with other solid or gaseous material. The portion of the dust that remains coupled to the gas, the advected dust, follows the gas as the disc evolves. As the advected dust flows into new regions, we assume it is brought into chemical equilibrium with the ambient material and use the new composition to calculate new abundances for the condensed matter at that location. Fig. 2 shows how our algorithm tracks the chemical composition of the decoupled dust in the disc mid-plane.

Figure 2.

Algorithm flowchart: A chart showing the calculations done at each step in the evolution of the disc. We calculate the chemical equilibrium of the gas and advected dust (green) at each radial location. A portion of the dust then decouples (black) using the decoupling time-scale given by equations (10) and (11). The remaining dust advects to neighbouring regions and a new chemical equilibrium is calculated. The process is repeated as shown by the blue arrow. Over time the composition of the dust changes as it advects, mixes with the gaseous elements in its new environment, and decouples to settle in the disc mid-plane.

2.2 Gaseous disc evolution model

For this work, we test our code by adopting the model used by Cassen (1996) to calculate the evolution of the gaseous disc (H₂ and He). The surface density at semimajor axis r and time t is given by

$$\begin{eqnarray*} \Sigma \left(r,t\right) = \Sigma _{0}\left(t\right)\mathrm{exp}\lbrace -\left[{r}/{r_0}\left(t\right)\right]^2\rbrace . \end{eqnarray*}$$

(2)

Here Σ₀(t) is the surface density near the central star, which is

$$\begin{eqnarray*} \Sigma _{0}\left(t\right) = \frac{M_\mathrm{d}}{{\rm \pi} {r_0}^2}, \end{eqnarray*}$$

(3)

and r₀(t) is the characteristic radius, which is given by

$$\begin{eqnarray*} {r_0}\left(t\right) = \frac{1}{GM_{\ast }}\left(\frac{J}{M_\mathrm{d}\Gamma \left({5/4}\right)}\right)^2, \end{eqnarray*}$$

(4)

where M_d is the mass of the disc, G is the gravitational constant, M_* is the mass of the central star (in this paper, we use 1 M_⊙), J is the angular momentum of the disc (we adopt J = 3 × 10⁵² g cm² s⁻¹), and Γ is the Gamma function (i.e. Γ(5/4) = 0.9064).

The evolution of the disc mass is governed by

$$\begin{eqnarray*} {M_\mathrm{d}}\left(t\right) = {M_\mathrm{d0}}\left(1+\frac{t}{t_\mathrm{e}}\right)^{-0.5}, \end{eqnarray*}$$

(5)

where M_d0 (⁠|$0.21 \ \rm M_{\odot }$|⁠) is the mass of the disc at t = 0, and t_e = 2.625 × 10⁴ yr is the characteristic evolution time-scale of the disc.

The mid-plane of the disc is heated by the accretion of the gas, so the mid-plane temperature, T(r, t), can be calculated from (Cassen 1994)

$$\begin{eqnarray*} T^4 = \frac{3 G \tau M_{\ast } \dot{M}_{\rm I}}{64 {\rm \pi} \sigma _{\rm SB} r^3}, \end{eqnarray*}$$

(6)

where |$\dot{M}_{\rm I}$| is the mass accretion rate to the central star, σ_SB is the Stefan–Boltzmann constant, and τ = κΣ/2 is the optical depth (κ = 4 cm² g⁻¹ is the opacity of gas with Solar composition). For simplicity, we use a constant opacity throughout the disc.

The mid-plane pressure is

$$\begin{eqnarray*} P = \frac{\rho \mathcal {R} T}{\mu }, \end{eqnarray*}$$

(7)

where |$\mathcal {R}$| is the gas constant, μ = 2.34 is the mean molecular weight in atomic mass units, |$\rho =\Sigma /\sqrt{2 {\rm \pi} }H$| is the mid-plane density, and |$H=\sqrt{\mathcal {R}T/\mu }/\Omega$| is the scale height of the gaseous disc.

We note that our disc is not in a steady-state. The mass accretion rate of the gas is given by

$$\begin{eqnarray*} \dot{M}(r) = -1.5\Sigma \nu \Big(1-4\Big(\frac{r}{r_0}\Big)^2\Big). \end{eqnarray*}$$

(8)

Over time, the value of r₀ increases from 3.2 AU at the beginning to 19 AU at 5t_e. The radius at which the accretion rate equals zero moves outwards, and the behaviour of the disc in the inner region is increasingly similar to a steady state disc. However, we overestimate the temperature in the inner region at early evolution of the disc, which affects slightly the location where the elements begin to condense and where the composition of the decoupled materials match that in terrestrial planets and chondrites.

Fig. 3 shows how the disc evolves with time. Fig. 3(a) shows the evolution of the normalized surface density of the disc. As time goes on, the surface density decreases in the inner region (radii less than several AU) while it increases in the outer region of the disc (i.e. the disc expands). The general trends of temperature (Fig. 3b) and pressure (Fig. 3c) with time are similar to the trend of surface density. Note that within 4 AU, the temperature decreases with time.

Figure 3.

Disc evolution: (a) Normalized surface density (Σ_0I is the initial value of Σ₀), (b) Mid-plane temperature, and (c) Mid-plane pressure for the disc at various points in time, (d) Radial velocity of the gas.

2.3 Motion and decoupling of condensable compounds

We assume that all elements in the gas phase have the same velocity as the H₂ and He (Jacquet, Gounelle & Fromang 2012). Their evolution is governed by

$$\begin{eqnarray*} \frac{\partial \Sigma _i}{\partial t}+\frac{1}{r}\frac{\partial \left(r V_r \Sigma _i\right)}{\partial r}=-F_{i}, \end{eqnarray*}$$

(9)

where Σ_i is the surface density of element i in the gaseous and advected dust phases, V_r is the radial velocity of H₂ and He, and F_i is the depletion rate due to the condensation. In this work, we only consider the advection of particles. We do not consider diffusion from turbulence, which may play a role in turbulent discs.

In our models, once an element begins to condense it also begins to decouple from the gas phase. The amount of the decoupled dust is determined by a decoupling rate that depends upon the surface density of the advected dust,

$$\begin{eqnarray*} \frac{{\rm d}\sigma _{i}}{{\rm d} t}=F_{i} =\epsilon \Omega \Sigma _{i,\rm adv}, \end{eqnarray*}$$

(10)

where σ_i is the surface density of element i in the decoupled dust, ϵ is the decoupling efficiency, Ω is the angular velocity around the central star, and |$\Sigma _{i,\rm adv}$| is surface density of element i in the advected dust. For the method that uses |$\rm T_{50}$|⁠, |$\Sigma _{i,\rm adv}$| is a step function around |$\rm T_{50}$| from 0 to Σ_i. For the chemical equilibrium method |$\Sigma _{i,\rm adv}$| is determined by chemical equilibrium calculation of the system. Initially, all the elements are of Solar abundance. After getting |$\Sigma _{i,\rm adv}$| and V_r at each time and each radius, we integrate Σ_i using the Euler method.

We assume the decoupling efficiency is uniform across the entire disc. The decoupling time-scale at 1 AU is given by

$$\begin{eqnarray*} t_{\rm dec}=\frac{1}{\epsilon \Omega \left(r=1 \rm AU\right)}. \end{eqnarray*}$$

(11)

At a given radial distance, the time for the dust to decouple depends upon the local orbital speed and we normalize the decoupling time to the orbital speed at 1 AU. A decoupling time-scale of zero corresponds to immediate decoupling, while an infinite time-scale implies no decoupling at all. Generally, the time-scale for dust to coagulate into planetesimals is 10³ to 10⁵ yr (Weidenschilling & Cuzzi 1993). For our fiducial models, we adopt a disc evolution time-scale of t_e = 2.625 × 10⁴ yr and a dust decoupling time-scale of t_dec = 1.5 × 10⁴ yr, which is the same as that used in Cassen (1996). All these models, and the parameters we used for them are shown in Table 2 where the historical model Cassen (1996) is labelled as MC, and the two new models are labelled as M1 and M2 (described below).

Table 2.

Different models and disc evolution and dust decoupling time-scales.

Time-scales	Historical model	Models with new |$\rm T_{50}$|		Time-dependent chemical equilibrium model
	(MC)	(M1)		(M2)
	Typical values	Typical values	Group 3	Typical values	Group 2	Group 3	Group 4
	(MC)	(M1G1)	(M1G3)	(M2G1)	(M2G2)	(M2G3)	(M2G4)
tdec yr−1	1.5e4	1.5e4	1.5e3	1.5e4	1.5e5	1.5e3	1.5e6
te yr−1	2.625e4	2.625e4	2.625e4	2.625e4	2.625e4	2.625e4	2.625e6

2.4 Model 1 and recalculating the 50 per cent condensation temperatures with grains

In our first model (M1), an element begins to condense out of the gas once the mid-plane temperature drops below |$\rm T_{50}$| (as was done in Cassen 1996). In order to compare our results with previous work, we first recalculate |$\rm T_{50}$| using the grains code (Petaev 2009). The elements used in grains are listed in Table 1. In this case, H and He account for 98.65 per cent (H: 71.34 per cent, He: 27.31 per cent) of the total mass with the remaining elements having Solar abundances (Lodders 2003). grains distributes 33 elements among 242 gaseous and 520 solid species by minimizing the Gibbs free energy of the system to identify the chemical equilibrium composition.

Cassen (1996) and Lodders (2003) estimate their |$\rm T_{50}$|s based upon a single pressure of 10⁻⁴ bar. However, the |$\rm T_{50}$| depends upon both the pressure and the starting composition, which, in turn, affects the condensation process throughout the disc. Here, we recalculate |$\rm T_{50}$| for each element by fixing the value of pressure (from 10⁻¹² to 1 bar) and lowering the temperature from 2500 K to 300 K. We note that this calculation does not involve the disc evolution – so there are only the gas and dust phases.

The condensation histories of several elements as functions of temperature are shown in Fig. 4. In Fig. 4(a), we see that |$\rm T_{50}$| for Mg, Si, and Fe are 1336, 1316, and 1335 K, respectively. Due to the difference between the temperatures where an element begins to condense and where it is completely condensed, the element with lower |$\rm T_{50}$| may have higher relative abundance in the dust phase at certain temperatures.

Figure 4.

Condensation curves: Proportions of elements (Na, Mg, Al, Si, P, K, Ca, Ti, and Fe) in the condensed phases as a function of the mid-plane temperature. The relative abundances of each element match the solar composition and the pressure used in these calculations is 10⁻⁴ bar.

Fig. 5 shows the effect on the |$\rm T_{50}$|s when changing the pressure from 10⁻¹² to 1 bar. When the pressure increases, the |$\rm T_{50}$|s also increases, and the relative values of |$\rm T_{50}$| may change. Also shown in that figure are the |$\rm T_{50}$|s used in Cassen (1996) and in Lodders (2003) with pressure 10⁻⁴ bar and Solar composition. While the values are generally consistent with each other, several elements have significant differences (e.g. K, Cu, Ni, and P). Woitke et al. (2018) calculated chemical equilibrium for solar abundance at 1 bar. The |$\rm T_{50}$|s at 1 bar here are consistent with the values from fig. 5 in their paper. With the new |$\rm T_{50}$|s shown above, we use equation (10) to calculate the decoupling of the elements as in Cassen (1996).

Figure 5.

50 per cent condensation temperatures: |$\rm T_{50}$| as a function of elements for various pressures (a) and pressure for various elements (b). In panel (a), we also add |$\rm T_{50}$| used in Cassen (1996) (violet circle) and Lodders (2003) (purple cross) with pressure 10⁻⁴ bar and Solar composition.

2.5 Model 2: Time-dependent chemical equilibrium model

In our second model (M2), we use the grains code to calculate the chemical equilibrium dynamically throughout the evolution of the disc rather than assuming an immediate condensation at |$\rm T_{50}$| followed by the decoupling process. In this model, |$\Sigma _{i, \rm adv}$| is calculated from the chemical equilibrium of the system, and grows smoothly as the temperature falls below the point where the elements begin to condense (see Table 1 and Fig. 4). Note that the temperature where an element begins to condense (T|$\rm _2$|⁠, in Table 1) is higher than |$\rm T_{50}$| (where half of the element has already condensed, see Table 1).

To save computational time, we do not determine the chemical equilibrium at every time-step of the disc evolution. Rather, we calculate the chemical equilibrium everywhere in the disc whenever the mid-plane temperature at 1.8 AU changes by 1 K. These temperature steps at this radius are small enough to keep the size of the temperature steps at all other relevant radii small enough that the chemical equilibrium calculation is minimally affected. We verified that these choices do not affect our results by testing the convergence of the condensation results with several simulations where we fix the temperature step at 1.8 AU to be several values in the range between 1 and 100 K. Fig. 6 shows that if the temperature step is less than or equal to 10 K, the differences between the results are less than 1 per cent. Consequently, we need not use a temperature step smaller than 1 K – which yields errors at a level near 10⁻³.

Figure 6.

Relative values of |$\rm R_i$| for different temperature steps: Relative values of |$\rm R_i$|⁠, for a sample of elements, as a function of the change in temperature that triggers the evaluation of the chemical equilibrium. |$\Delta \rm Ri$| is the difference between |$\rm R_i$| for temperature steps ΔT and that for 1 K. |$\rm R_i$| in y-axis is the value for 1 K. Evaluating the chemical equilibrium when the temperature at a distance of 1.8 AU drops by 10 K produces differences at the level of 1 per cent. The abundance differences shown in this plot are evaluated at 3 AU.

3 RESULTS

Our first new model (M1) produces results which differ only slightly from the model from Cassen (1996) (MC). Thus, we do not discuss the details of that model. Rather, we first discuss the evolution of the decoupled dust from the second model after which we compare the results of the two models (M1 and M2).

3.1 Evolution of the surface density of decoupled dust

To give an example of how a system evolves, Fig. 7(a) shows how Si and Ca in the decoupled dust evolve with time. Si first condenses at around 4 AU. At smaller radii, the temperature of the disc is too high for Si to condense, while at larger radii the gas surface density is low – yielding a surface density for the decoupled materials that is similarly low. As time goes on, the condensation region expands both inward (as the temperature of the inner region decreases) and outward (as the gas surface density increases due to viscosity and the outward flow of disc material).

Figure 7.

Condensed surface density and relative abundance over time: (a) Evolution of Si and Ca in the decoupled dust for the model (M2G1). The solid lines denote the evolution of Si while dashed lines indicate Ca. Over time the condensation fronts move inward as the temperature decreases and outward as the surface density rises.

The general trend of the evolution of Ca is similar to that of Si except the maximum surface density and the inner and outer boundaries of the condensed region change. The different maximum surface densities and outer boundaries are determined by the different gas surface densities of the elements and their condensation temperatures. The inner boundaries are correlated with the |$\rm T_{50}$| for the elements.

3.2 Comparison of the two models

In Fig. 8, we show the evolution of |$\rm R_i$| for Al, Si, S, Os, and Au at 3 AU for our two models M1 and M2. Al and Os are both refractory elements with |$\rm T_{50}$| for Al (1654 K) and Os (1812 K) both higher than that for Si (1316 K). Thus, their |$\rm R_i$| decrease over time as more Si condenses. The final ratios calculated from the chemical equilibrium are slightly larger than those calculated using |$\rm T_{50}$|⁠. Both S and Au are more volatile with their |$\rm T_{50}$|s being lower than that for Si: S (669 K) and Au (1196 K). Thus, their |$\rm R_i$| increases over time from 0 to their final values.

Figure 8.

Abundances versus time: |$\rm R_i$| as a function of time for Al, Si, S, Os, and Au at 3 AU. Al and Os are refractory elements, which begin condensing early while S and Au are more volatile and condense later. Both chemical equilibrium calculation (M2G1, solid lines) and fixed |$\rm T_{50}$| (M1G1, dotted lines) of the elements in the evolving disc are plotted (the models used in Cassen 1996 and in Fig. 9). The ratios of the final values of |$\rm R_i$| between the two models for Al, S, Os, and Au are 0.979, 1.03, 0.978, and 1.01, respectively, with refractory elements being underestimated and volatile elements being overestimated by the 50 per cent condensation temperature approximation. The corresponding temperature (grey line) is also plotted with the temperature scale along the right-hand side.

Since the dust begins to condense earlier in the chemical equilibrium calculation (M2 models) than it does in the |$\rm T_{50}$| calculations (M1 models) and because the surface density of the disc decrease quickly at early time in the inner region, the approximation from |$\rm T_{50}$| systematically overestimates the volatile abundance (3 per cent for S and 1 per cent for Au) while it systematically underestimates the refractory abundances (2 per cent each for Al and Os). The dotted lines in Fig. 9 show |$\rm R_i$| by using updated |$\rm T_{50}$|s calculated from the grains code (M1G1). There are more refractory elements in Fig. 9 than are shown in the results of Cassen (1996) because the grains code considers more elements.

Figure 9.

Equilibrium versus 50 per cent condensation models (M1G1 versus M2G1): Comparison of |$\rm R_i$| calculated from the two new models at the end of the simulation (5t_e, where t_e = 2.625 × 10⁴ yr). The chemical equilibrium model (M2G1, solid lines) and the fixed |$\rm T_{50}$| model (M1G1, short dotted lines) for the elements in the evolving disc are shown.

The general trends at 1 and 2 AU are similar to those from the historical model. At 3 AU, |$\rm R_i$| for elements with |$\rm T_{50}$| higher than Ca (1500 K) are almost the same, because the initial temperature at this radius is near |$\rm T_{50}$| for Ca and these elements begin condensing simultaneously. Because the temperature interior to 4 AU decreases with time, the value of |$\rm R_i$| for an element depends upon the initial value of the temperature at a given radius. If the initial temperature at a point is lower than |$\rm T_{50}$| for Si (1316 K), |$\rm R_i$| will be 1. If it is higher than |$\rm T_{50}$|⁠, |$\rm R_i$| will decrease with decreasing |$\rm T_{50}$|⁠. At 4 AU, the |$\rm R_i$| of the two models are largely the same for elements with |$\rm T_{50}$| higher than that for Cu (1100 K).

4 QUANTIFYING THE EFFECTS OF DECOUPLING AND DISC EVOLUTION TIME-SCALES

4.1 Changes to the decoupling time-scale

To illustrate the conditions under which our simulations can reproduce the chondrite abundances, we look at the effects of changes to the disc evolution time-scale and decoupling time-scale. We first change the decoupling time-scale, t_dec, to be an order-of-magnitude longer than the nominal case (1.5 × 10⁵ yr – M2G2; Fig. 10a). We find that these results are almost the same as those for the nominal case. Because the decoupling time-scale is very long, the material flows with the gas and leaves behind the same relative abundances of decoupled dust. Small ϵ does not affect the relative |$\Sigma _{i,\rm adv}$| and, consequently, not the σ_i. Thus, while less material decouples overall, the relative amount of the elements in the decoupled dust stays almost the same – a result that persists with even longer decoupling time-scales.

Figure 10.

Effects of changes to the decoupling time-scale (M2G2 and M2G3): These panels show the effect of changing the decoupling time-scale of the dust relative to the disc evolution time-scale. Values of |$\rm R_i$| at different radii and at time |$t=5\ t\rm {_e}$| (⁠|$t\rm {_e}=2.625\times 10^4\ yr$|⁠) calculated from the evolution of disc and using the grains code. The decoupling time-scales here are (a) 1.5 × 10⁵ yr (M2G2) and (b) 1.5 × 10³ yr (M2G3), which are of one order of magnitude longer or shorter than that in the typical case, respectively. The elements are ordered from high to low |$\rm T_{50}$| calculated from the grains code. In panel (b), we also plot the results calculated from |$\rm T_{50}$|s (M1G3, dotted lines). Note that with the longer decoupling time there is a consistent trend in the relative abundances with distance. With the short decoupling time we observe peaks in the abundances where the initial temperature in the disc correlates with the condensation temperature of the element at that particular radial distance. Also, with short condensation times, the innermost region (where the gas flows towards the star rather than away) the refractory elements are depleted somewhat as the incoming gas from more distant regions contains a higher proportion of volatile elements (see the text for details).

Fig. 10(b) shows the corresponding abundances at different radii when the decoupling time-scale is an order-of-magnitude shorter than the nominal time-scale (M2G3, 1.5 × 10³ yr). Compared to Fig. 9, the abundances at 3 AU peak at Ca, with a |$\rm T_{50}$| of 1500 K. The more refractory elements have lower |$\rm R_i$|⁠, because they condense and quickly decouple in the inner regions and do not flow outwards. The most refractory elements at 3 AU have the lowest |$\rm R_i$|⁠. Ca has the largest abundance because it condenses and decouples at the beginning of the simulation and additional Ca is advected into that region from the inner disc (where it was too hot for Ca to condense) while the more refractory elements remain behind. At 4 AU, a similar effect happens with Cu. The more refractory elements also have lower |$\rm R_i$| while the more volatile elements have higher |$\rm R_i$| for elements with |$\rm T_{50}$| larger than that for Ca (1100 K).

The overall effect of the relationship between the disc evolution time-scale and the dust decoupling time-scale (i.e. the planetesimal formation time-scale) is that there will be peaks in the radial distribution of different elements that depend upon their condensation temperatures. If planetesimals form quickly, then the more refractory elements will not flow with the gas, and regions fed by the warmer parts of the disc will be somewhat depleted in those elements. On the other hand, if planetesimals form more slowly (a long decoupling time-scale), then the advecting refractory materials will enrich neighbouring parts of the disc. Definitive statements about the consequences of this effect will require a more realistic disc model and lie beyond the scope of this work.

4.2 Changes to the disc evolution time-scale

To see how the disc evolution time-scale affects |$\rm R_i$|⁠, we change both the decoupling time-scale and the disc evolution time-scale to be two orders-of-magnitude longer than the typical case while keeping the ratio the same (M2G4). With the disc evolution time-scale two orders-of-magnitude longer than the typical case, the accretion rate from the disc to the central star is smaller by two orders-of-magnitude, and the temperature is approximately one third of the typical case (See equation 6).

Because the temperature across the disc is now much lower, and is even below |$\rm T_{50}$| of Si at 1 AU (1316 K), most of the elements condense immediately at the start of the simulation. The most volatile elements that are close to the central star condense later than the refractory elements and have lower |$\rm R_i$|⁠. The resulting abundances are shown in Fig. 11. The refractory elements at different radius all have the same |$\rm R_i$| – around unity. Comparing with the results of the typical case (Fig. 9) shows that the abundances can be strongly affected by the temperature history of the disc.

Figure 11.

Changing both disc and dust evolution time-scales (M2G4): Values of |$\rm R_i$| at different radii calculated from the evolution of disc and condensation from the grains code. The decoupling time-scale and disc evolution time-scales are two orders-of-magnitude longer than the typical case (⁠|$t\rm _{dec} = 1.5\times 10^6$| yr and |$t\rm {_e}=2.625\times 10^6$| yr, respectively). The lower temperature of the disc, except at the innermost location, causes most elements to condense simultaneously with Si. The elements are ordered from high to low |$\rm T_{50}$| calculated from the grains code.

5 DISCUSSION AND COMPARISON WITH OBSERVATIONS

To compare our results to the element compositions of the CM, CO, and CV chondrites we adopt the decoupling time-scale of 1.5 × 10⁴ yr and disc evolution time-scale |$t\rm {_e}=2.625\times 10^4$| yr (identical to Cassen 1996). MC model reproduces the general trends for volatile elements of the chondrites as do our two new models. However, the relative abundances are different at different radii. Our full condensation model most closely matches the observed abundances (CM, CO, and CV) at 3 AU, i.e. consistent with the values of the relative abundances of the chondrites.

The more refractory elements of CM, CO, and CV chondrites with |$\rm T_{50}$| higher than that for Ti (1583 K) have lower relative abundances (Fig. 9). With our nominal decoupling and disc evolution time-scales, the quick depletion of the refractory elements affects the transfer of elements to different radii. This effect can explain the low relative abundances of the most refractory elements. For example, the three most refractory elements Re, Os, and W have lower abundances compared to Ti in chondrites.

To the extent that the composition of the condensed materials reflects the composition of the terrestrial planets, our results also match the general trends for those planets. The average values of |$\rm R_i$| for refractory elements for Mercury, Venus, Earth, and Mars are 1.80, 1.17, 1.32, and 1.02, respectively. For volatile elements, the values are 0.0877, 0.274, 0.105, and 0.381 (See Fig. 1). Generally, the values for the refractories are high at small orbital distances (except Venus) and low for volatile elements – consistent with the results shown in Fig. 9.

Fig. 10(b) shows that |$\rm R_i$| for refractory elements can decrease towards small radii while increasing slightly for volatile elements. A short decoupling time-scale may explain the slightly lower refractory abundances of Venus compared to the Earth. The surface density and the temperature in the inner region of the model are higher than what we would expect from a more realistic disc model (Li & Sui 2017). For example, at one disc evolution time, the mid-plane temperature at 2 AU is around 1000 K, while at the same time, but using the disc model from Li & Sui (2017), 1000 K occurs near 1 AU. Therefore, if we obtain the same relative abundance for the terrestrial planets in our disc, the semimajor axes from our model will be larger than the actual locations of the planets in a real system.

Observations suggest that chondrule formation began after CAI formation and persisted for a few million years (Kita et al. 2000; Amelin et al. 2002; Connelly et al. 2012). Our models yield chemical compositions of the condensed materials, reproducing the elemental depletion patterns in some chondrites, but they do not address the formation of either CAIs or chondrules. Accounting for the differences between the formation of CAIs and chondrites would require simulating the evolution of the disc for several million years to track how the environment in the disc changes with time – roughly two orders-of-magnitude beyond our simulation times. In addition, since the pattern of the abundance of elements in CM, CO, and CV chondrites are measured from cm-scale samples (Wood 2005; Ciesla 2008), one should include some details of the planetesimal formation process in order to make more detailed comparisons with Solar system planetesimals.

6 CONCLUSIONS

We present two condensation models that predict the chemical compositions of the chondrites (CM, CO, and CV) and terrestrial planets. The first model follows Cassen (1996) by using |$\rm T_{50}$| as a fiducial value below which the elements will begin to decouple. In the second model, we calculate the chemical equilibrium of 33 elements by minimizing the Gibbs free energy of the system at each radius and specific time-steps in the evolving disc.

Compared with previous works, our newest model improves upon them by (1) using a chemical equilibrium calculation to model the condensation of the elements rather than fixed |$\rm T_{50}$|⁠, (2) significantly increasing the number of elements used in the calculation, and (3) embedding the elements in a dynamically and thermally evolving protoplanetary disc rather than choosing one or several snapshots in the evolving discs.

From our results we draw the following conclusions:

• For each element, the condensation front expands inward and outward from its initial condensation region. The inner edge of this region advances because the temperature decreases with time. The outer edge of the region retreats because of the outward movement of the gas in the disc, which increases the disc surface density.

• The relative abundances of the condensed elements in the disc vary with time and distance, decreasing for refractory elements and increasing for volatile elements within ∼3 AU. The resulting condensed materials are enriched in refractory and depleted in volatile elements within ∼3 AU. The chemical compositions of chondrites (CM, CO, and CV) are most closely matched by the condensate formed at 3 AU in these models.

• Approximating the condensation process with a fixed |$\rm T_{50}$| is similar to the full chemical equilibrium model. However, there are systematic differences from refractory to volatile elements. The relative elemental abundances, |$\rm R_i$|⁠, calculated from the chemical equilibrium are higher for refractory and lower for volatile elements within 3 AU. That is, the 50 per cent condensation model underestimates refractories and overestimates volatiles. Nevertheless, under certain circumstances, the |$\rm T_{50}$| model may be sufficient. The systematic differences can be as large as 10 per cent and may become more important as disc models are further refined.

• The abundances do not change significantly when the decoupling time-scale is equal to or longer than the disc evolution time-scale. Rather, the surface density of the material decreases because advected dust drifts out of the region of interest.

• If the decoupling time is short, there will be peaks in the elemental distribution as a function of the condensation temperature. The relative abundance of refractory elements are diminished at small semimajor axis because refractory elements from more distant regions decouple while the volatile elements will flow towards the star, diluting the refractory elements. This may explain the estimated lower abundance of refractory elements in Venus compared to Earth.

• When the initial temperature in the disc is high, the relative abundances of refractory elements are above unity, while these are low for volatile elements. When the initial temperature is low, the relative abundances of refractory elements are around unity since they condense at the same time as Si. Thus, we expect the relative elemental abundances, |$\rm R_i$|⁠, to be affected by the temperature history of the disc.

Using more elements in the calculations allows a more detailed comparison of with chondrite composition. In addition, by including the effects of the transport of matter in the disc, our results can provide more information on the history of dust within the disc and may provide better estimates for the dust in planet forming discs generally. Our model offers one explanation for both the decreasing of |$\rm R_i$| with |$\rm T_{50}$| and the cases that do not follow these trends. We notice that the thermal history can have a significant effect on the results. By including this aspect of the evolution of the disc we can improve our understanding of chemical compositions of the meteorites and planets.

Some limitations to our current disc model can be improved upon with future work. The surface density and the temperature in the inner region of our disc are high compared to a more realistic model (e.g. Li & Sui 2017). An improved model will result in a smaller semimajor axes for the formation location of the planetesimals and planets and may affect the relative elemental abundances.

ACKNOWLEDGEMENTS

The authors thank the referee for constructive comments that improved this manuscript. JHS, ML, and SH are supported by the NASA grant NNX16AK08G and National Science Foundation grant AST-1910955. ZZ acknowledges support from the National Science Foundation under CAREER Grant Number AST-1753168 and Sloan Research Fellowship.