Non-equilibrium condensation of the first Solar System solids

KineCond is a code developed to calculate the time-dependent condensation processes of minerals in the solar nebula, which is dominated by hydrogen. It is described in detail in Supplementary Section 3 with several validation tests, and here we give a summarized version. The elements considered in the system are H, He,O, Mg, Si, Fe, Al, Na, K, Ni, Ca, Cr, S and C. Here, the term element designates an atomic species. The system consists of a gas in interaction with 39 minerals. Most of the condensation codes published in the past and that allowed the investigation of the Equilibrium Condensation Sequence^3,4,5,7 rely on the Gibbs free-energy minimization (GFEM). GFEM calculates—for a given pressure and temperature—the most stable combination of minerals and gas; however, like any equilibrium calculation, it provides no information if the equilibrium state is established in a reasonable time, nor does it give the list of reactions through which the equilibrium state is realized (even though some additional physical arguments may help to determine those reactions, especially at high temperature when the number of minerals present is small). However, each reaction has its own kinetics, which depends on pressure, temperature and the local abundances of all elements. Thus, to design a time-dependent kinetic code, we must adopt a different strategy. We must explicitly specify the list of all reactions of interest, and advance each of them individually in a fully coupled way.

KineCond proceeds as follows: at each time step, we compute the number of atoms in the gas and the number of atoms in every mineral, keeping the total number of atoms constant. We assume that the pressure is constant and that only the temperature varies with time. The evolving variables of the system are: the number of moles of each molecule i in the gas (N_i^gas) and the number of moles of each mineral j in the system (N_j^min; Supplementary Table 2). We assume that gas–gas reactions are much faster than gas–mineral reactions and condensation reactions, so the gas molecular composition is always close to chemical equilibrium (then the gas molecular composition only depends on the current values of temperature and pressure, as well as N_i^gas). As we focus on gas–mineral processes, reactions are divided into two broad categories: (1) condensation or evaporation reactions; and (2) gas–mineral reactions. The rates of these reactions dictate the evolution of N_i^gas and N_j^min. The temperature varies linearly with time, dropping from 2,000 K to 130 K on the timescale T_c (ranging from 0.01 to 2000 years). The system evolves as follows: at each time t we first calculate the molecular composition of the gas (excluding mineral condensation) by computing the chemical equilibrium of the gas at T(t) and P and with elemental abundances N_i^gas(t). This is performed using the iconic Chemical Equilibrium with Application distributed by NASA (CEA-NASA) code⁴⁷, which includes about 1,500 gas species in total. The instantaneous gas molecular composition is then used to compute the different condensation and gas–mineral reactions and the rate at which they proceed. The different steps of the calculation are summarized in a flow chart provided in Supplementary Fig. 5. We detail these calculations below.

Condensation or evaporation reactions We follow a condensation and evaporation theory developed for forsterite evaporation in a H₂ gas^48,49, and we generalize it to many minerals. The net formation rate of a mineral j is the difference between an evaporation flux (J_j^e) and a condensation flux (J_j^c). Each of them must be computed explicitly. In the gas, the flux of any element E across a unit surface (in moles m^–2 s^–1) is calculated using the kinetic theory of gases^48,49:

where m is any gas molecule; \(\nu _m^E\) is the stoichiometric coefficient of element E in molecule m; P_m is the partial pressure of molecule m; µ_m the molar mass; and R is the ideal gas constant. The partial pressures of gas molecules (P_m) are obtained by running the CEA-NASA code. We now consider a mineral j with formula \(\\alpha _j^EE\\), where \(\alpha _j^E\) is the stoichiometric coefficient of element E in mineral j. The condensation flux of the mineral j (J_j^c) is determined by the smallest flux (over all elements E entering in its composition) so that:

γ_E,j is the sticking efficiency of atom E on mineral j. It ranges from 0 to 1 and is poorly constrained.

Here we set γ_E,j = γ = 0.1, as a standard value, for all minerals and all atomic elements.

The evaporation flux of the mineral j (J_j^c) is classically⁵⁰:

where P_j^sat is the saturating vapour pressure of mineral j at temperature T; P_j^sat is easily defined in vacuum and for a mineral that has an equivalent gas form with same formula (like H₂O or Fe).

However, most minerals do not have a corresponding gas form (like Mg₂SiO₄), making equation (3) not directly applicable to most minerals. Furthermore, the presence of the surrounding gas (mainly hydrogen and helium) modifies the chemical equilibrium and must be taken into account when computing the saturation vapour pressure. So, equation (3) is not directly applicable to compute the evaporation flux of all minerals immersed in a surrounding gas rich in H. To compute J_j^e from the kinetic theory of gases, we follow a strategy by which the evaporating flux is directly computed without the need to compute the saturating vapour pressure^48,49. When a mineral is in equilibrium with its surrounding gas, the evaporating and condensation fluxes exactly balance (equations (2) and (3)).

So, determining the evaporation flux of the mineral j is equivalent to determining the condensation flux of the mineral j in equilibrium with the surrounding gas at pressure P and temperature T. Following previous work⁴⁸ focusing on forsterite equilibrium with H₂ we have calculated the vapour composition at equilibrium with a mineral immersed in H₂ gas and derived the evaporative flux using equation (2). This calculation was performed using the CEA code for 39 minerals (Supplementary Table 2). The pressure of the H₂ gas varied in the 10⁻¹⁰ bar < P < 0.1 bar range, the temperature T was varied in the 130 K < T < 2,200 K range, and fO₂ was varied by ten orders of magnitude. All of these data are compiled into lookup tables and are interpolated at current (T, P, fO₂). So in KineCond the ideal evaporative flux (F_j) from the mineral j is obtained by reading a precalculated lookup table and interpolated between tabulated values at current values of ln(P), T and ln(fO₂). Examples of computed evaporative fluxes are shown in Supplementary Fig. 10. The result is multiplied by γ = 0.1 so the effective evaporative flux is:

The time derivative of number of moles of any mineral j due to competing condensation and evaporation processes is then:

Where S is the surface of contact of the grain with the gas; J_j^c and J_j^e are given by equations (2) and (4), respectively. Our code implements a first-order time solver (Euler), and during a time step dt, the number of moles of every mineral evolves according to equation (5). To conserve the total number of moles of each element, the atoms released and removed from the gas (N_i^gas) are counted according to the stoichiometry of every mineral. We adopt an operator-splitting approach in which condensation reactions during the dt time step are treated first (Section 3.2 of Supplementary Information), and gas-mineral reactions are treated in a second step (Section 3.3 of Supplementary Information).

The surface of contact of minerals with the gas (S) depends on the radius of the grain (r) and the number of grains (N) so that S ≈ N4πr². The self-consistent calculation of N and r necessitates the computation of the time-dependent nucleation process and that both the sticking and fragmentation processes of minerals during their settling and growth in the turbulent solar nebula are taken into account. So far, such a coupling (mineral condensation, coagulation and fragmentation) has never been realized, and some models take some processes into account (coagulation or fragmentation^51,52, or metal condensation²²). Models coupling dust settling, coagulation and fragmentation show that dust size distribution reaches rapidly steady-state inward 10 AU in a few orbital periods^51,53, and that which determines dust size is the equilibrium between coagulation and fragmentation in turbulence, rather than mineral growth. Following these lines, and to make the calculation tractable, we have used a characteristic dust size, which gives a gas–mineral surface of contact, rather than computing self-consistently a mineral growth model, and we limit ourselves to an order-of-magnitude calculation. In other terms, we assume that all minerals (despite their mass) have a surface in contact with the gas equivalent to a sphere with radius r = 10 µm (where r is a free parameter of the model). This size is typical of the minerals observed in chondrites. The number of grains in the system (N) should be controlled by the number of nucleation sites that first appear in the gas. As aluminium is the most refractory atom in our model, we approximated N ≈ M_Al/(4/3πρ_Alr³) where M_Al is the total mass of aluminium in our system and ρ_Al is the density of aluminium, so S = 3M_Al/(ρ_Alr). This does not mean that all minerals are 10 μm in radius, but rather that on average the total surface of contact of a mineral is equivalent to a population of minerals with an average surface of a sphere with a 10-µm radius (it could be fractal in shape). Of course, our calculation may not be accurate at the beginning of coagulation when the particle size is close to the monomer size, but changing the monomer size by a factor of 1,000 only changes by 4% the growth time and does not change the final size²². Particles with a radius of 70-μm were formed in about ten weeks⁸ (with a cooling rate of about 100 K per year at 10⁻⁴ bar; a pressure typical of the 0.1–1 AU region), which is comparable to the orbital period at the distance of Mercury. Other works find the formation of 10-μm grains in a few orbital periods at 1 AU (refs. ^51,53) or at 5 AU (ref. ⁵²). Of course, larger minerals can be formed, but 10 μm is typical of what is observed in meteorites.

Gas–mineral surface reactions (nebular reactions)

In addition to condensation and evaporation reactions, condensed minerals can also interact with the gas, leading to mineral transformation. Nebular reactions are those where the gas interacts with a pre-existing mineral (M1) and forms a new mineral (M2) and consuming M1. The number of such reactions is potentially infinite and for now KineCond implements reactions with the following generic form:

Where M1 and M2 represent two minerals; E1 and E2 represent any element in gas form; and α, β₁ and β₂ are stoichiometric coefficients (normalized so that the stoichiometric coefficient of M1 is 1). For example, Mg₂SiO₄(s) + 1 Si + 2O ⇒ MgSiO₃ (s).

Elements E1 and E2 can be implied in the reaction in any molecular form (silicon could be in the molecular form SiO, Si, SiO₂ and so on), but the radical of the molecule that is not used in the reaction is released into the gas and does not enter into the calculation of the reaction constant (due to the difference of the formation ∆G between the left and right sides of the equation; see Supplementary Sections 3.4 and 3.4.1 for a detailed calculation). For now, we are limited to reactions in which mineral M1 does not lose atoms to the gas. If M1 did, this could be an incongruent evaporation process, and this nebular reaction would be inconsistent with our hypothesis of treating condensation or evaporation as congruent processes using the Hertz–Knudsen formalism described above. Using combinatorial analysis, we found that 55 reactions of this type are possible with our mineral selections (Supplementary Table 2). After many tests, only 38 reactions were kept, the others playing a more minor or no role (Supplementary Table 3).

For a given temperature and for a given gas composition, we first compute whether any reaction listed in Supplementary Table 3 is kinetically possible, that is, if the mineral M2 is formed (see Supplementary Section 3.4 for a detailed calculation). If the reaction favours the formation of the mineral M2, then we calculate the rate at which M2 is formed and M1 disappears.

Rate of nebular reactions

We use the simple collision theory²³ model to compute the reaction rate of reaction 6, but modified the model to take into account improvements in our understanding of reaction rates on the surface of the mineral²⁴. We first determine the flux of incoming elements E1 and E2 by summing the fluxes of all molecules in the gas, which carry elements E1 or E2 (equation (1)) called J^c(E1) and J^c(E2), where c is condensation. We call them elementary fluxes. The smallest of the two fluxes (weighted by β₁ or β₂) controls the rate of progression of the reactions. So if reaction 6 is thermodynamically possible, then the rate at which mineral M1 appears (and mineral M2 disappears) is:

Equation (7) assumes that all collisions lead to the reaction and overestimates the reaction rate; it must therefore be corrected. Laboratory experiments show that two regimes of reaction rates exist^24,33,54,55: the linear and parabolic regimes. In the former, each molecular collision has a certain probability of producing a chemical reaction, parameterized by the activation energy \(E_a^l\). So the linear rate is:

We recover here the simple collision theory model²³. However, laboratory experiments show that after a first period during which the reaction rim grows linearly with time at the mineral’s surface^33,54, the reaction switches to a parabolic regime where atomic diffusion across the reactive rim limits the reaction rate. In this regime, the rim grows with the square root of time. In that case, the reactive layer with thickness H grows like H² = k(T)t where k(T) is a diffusion coefficient that depends on temperature (m² s^–1) and t is time³³. It is usual to write \(K(T)=Ce^-E_ap/RT\), where C and \(E_\rma^p\) represent a prefactor and the activation energy (in the parabolic regime), respectively. This process can be described by a Fick diffusion law where the reaction rate (in mol m⁻² s⁻¹) is

where c(x) is the concentration of mineral M2 at location x (x = 0 corresponds to the surface of mineral M1). To simplify the calculation, we assume that c(x) drops linearly within M1, so that dc/dx ≈ 1/H. So we solve:

where H is the thickness of the mineral M2 layer above the M1 surface and S the grain surface. In KineCond H is calculated as H = µ_M2N_M2/(Sρ_M2) where µ_M2 and ρ_M2 are the molar mass and density of mineral 2. To avoid non-physical high reaction rates when H is close to 0 we bound the parabolic rate to be always smaller than the rate of incoming atoms to the M1 mineral’s surface.

So, the reaction rate in the parabolic regime reads:

Unfortunately, there are only a handful of laboratory measurements²⁴ and most gas–mineral reactions are undocumented. For magnetite formation, laboratory experiments give \(E_\rma^p\,\approx \,90\,\mathrmKJ\,\mathrmmol^\mbox–1\) (ref. ⁵⁴), whereas for Troilite (FeS) the activation energy in the parabolic rate is reported as \(30\,\mathrmkJ\,\mathrmmol^\mbox–1 < E_\rma^p < 70\,\mathrmkJ\,\mathrmmol^\mbox–1\) and in the linear regime is \(28\,\mathrmkJ\,\mathrmmol^\mbox–1 < E_\rma^\rml < 94\,\mathrmkJ\,\mathrmmol^\mbox–1\) (ref. ³³). These measurements were done in the (T, P) range of stability of both minerals. In contrast, the rate of forsterite to Enstatite (Mg₂SiO₄ + 1Si + 2O → 2MgSiO₃) was measured at a temperature well above the stability of Enstatite or Forsterite (>1,700 K) and E_a^p ≈ 500 kJ mol^–1 was found⁵⁶. For magnetite, troilite and enstatite activation energies, we use the laboratory values for the parabolic regime. For the linear regime and for all other reactions due to many uncertainties and lack of data²⁴, we investigate different end-member scenarios defined below:

• FNR: (end member) E_a^l = 0 kJ mol^–1 and E_a^p = 0 kJ mol^–1

• MNR: E_a^l = 20 kJ mol^–1 and E_a^p = 20 kJ mol^–1

• SNR: E_a^l = 80 kJ mol^–1 and E_a^p = 500 kJ mol^–1.

The prefactor coefficient in k(T), C, is determined so that the linear and parabolic fluxes connect smoothly for a rim thickness of H = 1 µm. This also prevents unrealistically high reaction rates for small values of H (a well-known problem of Fick’s law). So C is equal to the linear rate divided by 1 µm, that is, in the range of transition rim thicknesses measured for troilite formation³³. For the enstatite formation experiment, the transition thickness is below 10 μm (ref. ⁵⁶).

Of course, the above procedure does not reflect the vast richness of gas–grain reactions and suffers from an important lack of experimental data. However, we have performed numerous tests, varying the values of \(E_\rma^l\) and E_a^p. We found that changing these values does not significantly impact our results as the processes investigated here are dominated by condensation processes, rather than gas-mineral surface interactions.

Putting all things together, the user first specifies the gas composition (solar in general), pressure and cooling time (in years). Temperature will decrease linearly from 2,000 K to 130 K in time T_c. The equivalent cooling rate is therefore (2,000 − 130)/T_c (K per year). Each time step is decomposed as follows.

1. 1.

   Compute the gas molecular composition at a given temperature and pressure (using only atoms present in gas form) using the CEA-NASA equilibrium code⁴⁷.

2. 2.

   Compute the resulting gas atomic fluxes (equation (1)).

3. 3.

   Condensation/evaporation phase: calculate the condensation and evaporation fluxes for every mineral (equation(5)), and evolve the number of minerals accordingly and the number of elements in the gas.

4. 4.

   Nebular reactions phase: Determine which nebular reactions occur in the gas (Supplementary Table 3) and compute the rate of production of mineral M2 and the rate of destruction of mineral M1 (equations (9) and (12)).

5. 5.

   Update the number of elements in the gas, enforcing mass conservation.

6. 6.

   Increment time and go back to the first step.

Several comparisons of KineCond outputs against various laboratory experiments and a comparison with the classical equilibrium condensation sequences obtained by GFEM are shown at the end of Supplementary Section 3. We find good, to very good, agreement for all these tests.