© ESO 2018

1 Introduction

Deriving the chemical composition of gases in the atmospheres of stars, brown dwarfs and planets is one of the most fundamental modelling tasks in astronomy. All hydrodynamical, evolutionary and radiative processes depend on the thermodynamics of the gas. In particular, the interpretation of spectral observations requires detailed knowledge about the gas composition. Only then, can the inverse problem of inferring the internal physico-chemical structure of an object in space from spectroscopic observations be successfully solved.

Instrument development has progressed tremendously in recent years, allowing us to use a large range of wavelengths to analyse for example the chemical composition of brown dwarfs (see Allard et al. 1997; Helling & Casewell 2014 and references in these reviews) and exoplanets (e.g. Désert et al. 2008; Grillmair et al. 2008; de Kok et al. 2013; Kreidberg et al. 2014; Fraine et al. 2014; Birkby et al. 2017). We can also revisit the long-standing question of what drives the mass loss of oxygen-rich asymptotic giant branch (AGB) stars (Gail & Sedlmayr 1998; Woitke 2006; Höfner & Andersen 2007; Golriz et al. 2014; Gail et al. 2016; Decin et al. 2017), including the problem of seed particle formation and the search for the first condensate in space (see e.g. Patzer et al. 1995; Patzer 2007). High-resolution spectroscopy of molecular lines from CO and H₂O are also used, for example, to derive winds properties of extrasolar planets (Brogi et al. 2016). Infrared instruments on board of the James-Webb Space Telescope are expected to revolutionise our understanding of cool objects such as exoplanets, brown dwarfs and cool stars (e.g. Marley & Leggett 2009; Beichman et al. 2014). High-resolution ground-based spectroscopy, for example with CARMENES (e.g. Sanchez–Lopez et al. 2017), aim at a better understanding of late-type main-sequence stars and their low-mass planets. In the more distant future, PLATO (Rauer et al. 2016) will combine astroseismology for solar-like stars with the search for habitable, Earth-like planets. Exploring such high-end quality data requires substantial modelling efforts for which deriving the chemical gas composition is one of the key issues. Such modelling needs to be fast and sufficiently accurate over a wide range of temperatures and pressures. The predictions of the molecular composition of the observed gases depends critically on the availability and precision of thermo-chemical data.

In the first approximation, the gas is assumed to be in local thermodynamical equilibrium (LTE), so the gas composition can be calculated from the thermodynamical principle of minimisation of the system Gibbs free energy (e.g. White et al. 1958; Eriksson 1971). Thermo-chemical equilibrium is a part of the set of LTE assumptions used in particular to compute the abundances of molecules. Pioneering work in this area was done for the atmospheres of cool stars where molecules are the main opacity carriers (e.g. Tsuji 1965; Auer & Mihalas 1968; Gustafsson 1971). To date, applications of thermo-chemical equilibrium are widespread in astronomy, for example in gas giant planet atmospheres, brown dwarfs and low-mass dwarf stars (Allard & Hauschildt 1995; Tsuji et al. 1996; Lodders & Fegley 2002; Visscher et al. 2010; Marley & Robinson 2015).

Deviations from thermo-chemical equilibrium are generally expected to occur in diluted gases, in particular due to UV-fields, X-rays or cosmic rays as, for example, in protoplanetary discs (e.g. Semenov & Wiebe 2011; Woitke et al. 2016), in the upper layers of planetary atmospheres (e.g. Moses 2014; Rimmer & Helling 2016), or in the case of fast hydrodynamical processes where the dynamical timescale becomes smaller than the chemical timescale, for example in interstellar shocks (e.g. Hollenbach & McKee 1989).

Griffith & Yelle (1999) and Cooper & Showman (2006) argue for strong chemical dis-equilibrium effects and higher CO abundances in BD and exoplanet atmospheres. Concerning the atmospheres of cool giant planets, Moses (2014) have shown that thermo-chemical equilibrium only prevails in the deepest layers (p > a few bar), whereas transport-induced quenching and photochemistry are likely to be important in the upper layers. However, recent kinetic chemical models (Visscher et al. 2006, Visscher et al. 2010; Zahnle et al. 2009; Line et al. 2010; Moses et al. 2011; Kopparapu et al. 2012; Venot et al. 2012) have shown that in hot exoplanet atmospheres (T ≳ 1200 K), the chemical timescales are in fact quite short, and hence thermo-chemical equilibrium prevails. For the hottest planets, Line & Yung (2013) and Oreshenko et al. (2017) have shown that at a pressure of 100 mbar CH₄, CO, H₂O and H₂ should be close to chemical equilibrium even when fast atmospheric mixing is assumed. Line & Yung (2013) have pointed out that spectral observations probe layers in exoplanet atmospheres where quenching by atmospheric mixing can be relevant, but not photochemistry.

At low temperatures (≲2000 K), the formation of solids and liquids becomes an essential part of the problem of finding the equilibrium composition. The condensates selectively consume most of certain elements (e.g. Al, Si, Mg, Fe), whereas others remain in the gas phase down to substantially lower temperatures (e.g. S, Cl). The abundant elements H, He, C, N, O remain mostly gaseous before even they, eventually, will condense for example in form of water, ammonia and methane ices like on Titan. The remaining gas phase element abundances can differ by orders of magnitude from the initial, undepleted (e.g. solar) abundances before condensation (e.g. Woitke & Helling 2004; Juncher et al. 2017), which is essential to understand the spectroscopic appearance of cool and cold objects like AGB stars and their winds (Gail et al. 1984; Sedlmayr & Dominik 1995; Helling et al. 1996; Woitke 2001; Jeong et al. 2003), or the formation of clouds in brown dwarfs and planetary atmospheres (Lunine et al. 1986; Tsuji et al. 1996; Marley et al. 1999; Allard et al. 2001; Burrows et al. 2002; Tsuji 2005; Marley et al. 2002; Helling et al. 2008).

Since dust grains can easily de-couple from the gas phase (e.g. rain-out in planetary atmospheres, gravitational settling in protoplanetary discs, acceleration by radiation pressure in AGB star winds), the condensed elements may be carried away, gather in other places (e.g. oceans) or re-evaporate in other places (e.g. comet tails). These effects are very likely to cause most peculiar local element abundances, with large effects on gas phase abundances and spectroscopic appearance. Such effects are rarely studied in detail (see, however, Booth et al. 2017). Instead, researchers tend to rely on known sets of element abundances, in particular the solar element abundances (e.g. Asplund et al. 2009) and those derived from Earth rock analysis.

Deviations of element abundances from these known sets have mostly been treated by reducing all metal abundances by the same factor, that is, by changing the metallicity [M/H] (e.g. Helling & Lucas 2009; Witte et al. 2009; Crossfield et al. 2013; Hu & Seager 2014; Morley et al. 2017) or by changing the carbon-to-oxygen ratio (e.g. Madhusudhan et al. 2011; Madhusudhan 2012; Helling et al. 2017; Van Eck et al. 2017). Gaidos (2000) was one of the earliest works suggesting the importance of C/O in an exoplanet context. Mollière et al. (2015) have studied the impact of C/O on the (p, T)-structure and emergent spectra of irradiated planets. Planetary atmospheres with very different element abundances have been studied by Mbarek & Kempton (2016) and Mahapatra et al. (2017), considering an atmosphere made of evaporated rock.

In order to derive the composition of a gas in thermo-chemical equilibrium, different approaches are used by different groups working on stellar, brown dwarf and planetary atmospheres. Three major schools are: (i) use of law of mass action and molecular equilibrium constants based on Gibbs free energy data (Tsuji 1973; Stock 2008; Bilger et al. 2013), (ii) use of law of mass action and molecular equilibrium constants based on partition functions (Allard & Hauschildt 1995; Gustafsson et al. 2008; Barklem & Collet 2016; Van Eck et al. 2017), and (iii) minimisation of the system Gibbs free energy (Lodders 2003; Ackerman & Marley 2001; Miller–Ricci et al. 2009; Mbarek & Kempton 2016; Blecic et al. 2016). Helling et al. (2008) provide a summary of approaches used in brown dwarf and planetary atmospheres.

Approaches (i) and (ii) usually lead to a system of N algebraic equations for N unknowns, which can be solved by any root-finding algorithms, for example by the Newton–Raphson method. In contrast, the Gibbs free energy minimisation technique (iii) requires special numerical minimisation algorithms (Spang 1962), for example the Dantzig-simplex method, the Hessian-conjugate gradient method, or the Lagrangian steepest-descent method (Blecic et al. 2016). The Gibbs free energy minimisation approach can easily be generalised to include equilibrium condensation, which is the reason why this approach was introduced originally. However, finding the root of a coupled non-linear equation system F(x) = 0 can be done numerically in much more efficient ways than finding the minimum of min, because the vector F has much more information than the scalar F. Therefore, Gibbs free energy minimisation codes tend to be considerably slower.

Beyond thermo-chemical equilibrium models for gas phase and condensates, hybrid methods have been developed where the condensation is treated time-dependently, but the concentrations of gaseous molecules are calculated in thermo-chemical equilibrium (e.g. Gail et al. 1984; Sedlmayr & Dominik 1995; Woitke & Helling 2004; Woitke 2006; Höfner & Andersen 2007; Höfner et al. 2016; Helling & Fomins 2013), following the idea that dust formation is by far the slowest process, causing the first deviations from LTE.

The probably most common approach used in the retrieval community, however, is not to solve any equations at all, but to treat the molecular concentrations as free parameters which are then adjusted to match observations (see e.g. Table 1 in Helling et al. 2008). This is partly done for reasons of feasibility, as the retrieval technique requires to run tens of thousands of models, but also to allow the method to retrieve non-equilibrium effects whereas enforcing equilibrium would rather be considered as a limitation. Benneke (2015), Line et al. (2016) and Oreshenko et al. (2017) have shown that it is feasible to include equilibrium chemistry inexoplanet retrieval methods.

This paper summarises in Sect. 2 the theoretical background of our approach to compute the abundances of molecules and condensates in thermo-chemical equilibrium. Section 3 gives an overview of the thermo-chemical data available for molecules and condensates. We assess the level of (dis-)agreement between various data sources. Section 4 introduces our thermo-chemical equilibrium code. Section 5 present our results concerning the spectroscopically most active molecules, a benchmark test, and the condensation sequence of elements. We demonstratethat equilibrium condensation leads to a significant increase of the C/O ratio in the gas phase where phyllosilicates play a particular role. Section 6 re-addresses the question of the first condensate in space, which is of particular interest for dust formation in AGB star winds. Section 7 contains our summary and conclusions.

2 Chemical and phase equilibrium

2.1 Molecular equilibrium

Let us consider a molecule A_aB_bC_c made of three elements A, B, C, where a, b, c are the stoichiometric factors. Guldberg’s law of mass action (e.g. Berline & Bricker 1969) is given by (1)

where the p_i = n_ikT are the partial pressures [dyn cm⁻²], n_i the particle densities [cm⁻³] and p^⊖ is a standard pressure. is the Gibbs free energy of formation [J mol⁻¹] of the molecule at standard pressure from neutral atoms at the same temperature (2)

The equilibrium constants k_p are introduced as (3)

where (4)

As we are using cgs-units in this paper, the k_p have units . T is the gas temperature [K], R the ideal gas constant [J mol⁻¹ K⁻¹] and k the Boltzmann constant [erg K⁻¹].

To determine the chemical composition of the gas, we solve the element and charge conservation equations as (5)

where ϵ_k are the element abundances normalised to hydrogen (ϵ_H = 1), n_{⟨H ⟩} is the total hydrogen nuclei density n_⟨H⟩ = ∑_is_i,H n_i. n_i denote all gas particle densities including free electrons, neutral and charged atoms, and neutral and charged molecules. s_i,k is the stoichiometric factor of element k in gas particle i. We assume charge neutrality by the inclusion of the charge as an additional element “el” with zero abundance (ϵ_el = 0) where s_i,el = 0 for neutrals, s_i,el = +1 for ions and s_i,el = −1 for cations. The free electron has s_el,el = −1. The gas density ρ [g cm⁻³] is given by (6)

where the second part of Eq. (6) follows from Eq. (5), and the total gas pressure is given by (7)

m_i are the gas particle masses and m_k are the masses of the elements. The total particle density n = ∑_in_i and the mean molecular weight μ = ∑_im_i n_i∕n are results of the computations and hence depend on density and temperature. The total hydrogen nuclei particle density n_{⟨H ⟩} is always proportional to the mass density ρ, whereas n is not, causing the gas to deviate from an ideal gas.

After elimination of all molecular particle densities from Eq. (5) by using the k_p (T)-data (see Eq. (3)), Eq. (5) becomes a system of non-linear algebraic equations with K unknowns, namely the atomic partial pressures and the electron partial pressure, for example {p_k | k = H, …, W, el} where H is hydrogen and W is tungsten. All results depend solely on n_⟨H⟩, T and ϵ_k.

If a solution is requested for a given pressure and temperature, an iteration is performed where the mean molecular weight μ is initially guessed, Eqs. (7) and (6) are used to compute n_{⟨H ⟩} from p and T, the equilibrium chemistry is solved with n_⟨H⟩ and T, and finally μ is re-calculated. We find this iteration to converge after 1–5 calls of the equilibrium chemistry.

For high temperatures, T ≳ 1000 K, the solution is easy to find numerically as all particle concentrations stay within 10^{−a few ten}. However, for T → 500 K, the results become increasingly more extreme, for example the electron concentration approaches 10⁻³⁰, and for T → 100 K the problem turns into a numerical nightmare with some k_p > 10^{+several 1000} and some atom and electron concentrations <10^{−several hundred}. We overcome these problems by switching to quadruple precision arithmetics at low temperatures and by applying an iterative procedure according to the hierarchical order of elements to provide extremely good initial guesses of the unknowns for the final Newton–Raphson iteration. We explain the details of the numerical approach of the new GGCHEM-code in Appendix A.

2.2 Condensed phases

The supersaturation ratio of a condensate j is given by (8)

where p_j^vap(T) is the vapour pressure (9)

In phase equilibrium we have (10)

whereas S_j > 1 (supersaturation) is not possible in phase equilibrium. Equation (10) states that the partial pressure of any molecule is limited in phase equilibrium, if a corresponding condensed phase exists. Once p_j exceeds that threshold (its vapour pressure), that molecule condenses¹ and leaves the gas phase, which reduces a few gas phase element abundances ϵ_k. This reduction causes not only p_j to drop, but also affects on all other partial pressures in the gas (usually reduces them). This process continues until phase equilibrium is established, see Appendix B in Helling et al. (2008).

For condensates which have no corresponding molecule in the gas phase (which is true for most minerals), Eq. (8) is not applicable. However, we can consider the fictive nominal molecule of that condensate, for example j = A_aB_bC_c with unknown G^⊖ (A_aB_bC_c, T), and apply Eqs. (8) and (9). Using also Eqs. (3) and (4), we find the generalised supersaturation ratio to be (11)

where (12)

is the Gibbs free energy of formation of the condensed phase A_a B_b C_c at standard pressure from atoms in the gas phase at temperature T. The unknown Gibbs free energy G^⊖(A_aB_bC_c, T) of the fictive molecule cancels out.

The addition of Eq. (10) to chemical equilibrium models (→equilibrium condensation or phase equilibrium models) leads to a considerable complication of the mathematical structure of the problem to solve. The case differentiation in Eq. (10) means that we do not know a-priori which condensates to consider as only those which finally result to be present contribute to the number of equations to be solved (in form of S_j = 1), whereas the supersaturation ratios of all other condensates do not matter as long as they stay <1. Clearly, at very high temperatures, where no condensates are stable at all, the problem falls back to the pure gas phase equilibrium discussed in the previous section.

Our numerical method in GGCHEM is based on a Newton–Raphson iteration with nested calls of the equilibrium chemistry with reduced gas phase element abundances ϵ_k due to condensation, see Appendix B. This reduction becomes extremely large at low temperatures – up to hundreds of orders of magnitudes. The additional unknowns to be solved for are the reduced gas phase element abundances as affected by the condensation of the selected solid and liquid phases. This selection may change during the iteration. Albeit being somewhat sophisticated, we find that this method is fast and very accurate down to 100 K when using quadruple precision arithmetics.

3 Thermo-chemical data

3.1 Molecular equilibrium constants

The k_p-data of the various ions and molecules must usually be fitted as function of temperature in some way, before they can be used in models. An alternative approach was recently proposed by Blecic et al. (2016), who use all available NIST-JANAF data-points directly by means of internal spline-fits. Various fit-functions have been proposed as summarised below. From Eq. (4) we have the general relation between and k_p as (13)

where n = a + b + c is the sum of stoichiometric coefficients in the molecule, for example n = 3 for H and n = 1 for H.

Tsuji (1973) has used the following fit-function (14)

for 335 molecules, where θ = 5040∕T and a_i are the fit coefficients. All k_p formulae listed here are valid in cgs units.

Gail & Sedlmayr (1986) have used (15)

which is the functional form that was used in the old GGCHEM code until recently (e.g. Helling et al. 2006; Bilger et al. 2013; Helling et al. 2017). The old data collection had 205 molecules.

Sharp & Huebner (1990) have used (16)

fitting directly the [cal mol⁻¹] of 184 molecules. Here, the standard pressure is p^⊖ = 1 atm and R_cal = 1.987 cal mol⁻¹ K⁻¹.

A new approach was presented by Stock (2008) (17)

who fitted the dimensionless quantity for 924 molecules, to be used with p^⊖ = 1 bar. This functional form has a very smooth and stable behaviour towards low temperatures. It also gives automatically more weight to the data points at low temperatures, so we have adopted this functional form and most of their data in the new GGCHEM code.

Barklem & Collet (2016) have recently published partition functions for 291 diatomic, and charged diatomic molecules, from which k_p is derived, but without providing a fit formula. We note that their data is in SI units, and the definition of their k_p is different from ours, in particular for the charged molecules. We have done all necessary conversions (see Appendix C) and fitted the Barklem & Collet data with a Stock-function.

Figure 1 compares the fitted data for two example molecules found in all 5 datasets. The deviations are as large as 10 kJ mol⁻¹ for HS and about 5 kJ mol⁻¹ for CN. At T = 300 K, an uncertainty of ±6 kJ mol⁻¹ (0.06 eV) translates into an uncertainty in k_p of about one order of magnitude, same for ±60 kJ mol⁻¹ (0.6 eV) at T = 3000 K. The complete comparison catalogue (Worters et al. 2017) includes 2782 individual datasets for 1155 molecules, where similar plots as shown in Fig. 1 can be found for all molecules for which we found at least two different data sources. We observe astonishingly large deviations in some cases. The reasons for these deviations are not entirely clear, but some factors could be

• different primary data sources,

• extrapolation errors,

• the “art of fitting”, and human errors.

Gail & Sedlmayr (1986), Sharp & Huebner (1990) and Stock (2008) have all based their fits on the JANAF database, but using different editions. Gail & Sedlmayr have used the 2nd edition (Stull & Prophet 1971), Sharp & Huebner have used the 3rd edition (Chase et al. 1982; Chase 1986), and Stock has used the 4th edition (Chase 1998). All three works include some Tsuji-data for molecules missing in JANAF. Stock (2008) has refitted the Tsuji-data for T > 1000 K with his fit-function that is more reliable concerning extrapolation towards lower temperatures. Despite using more or less the same original data sources, some molecules show considerable deviations even at medium temperatures (see Worters et al. 2017). We conclude that errors of about half an order of magnitude for k_p can easily occur, especially at low temperatures, depending on which fit-formula is chosen and how the fit was obtained in detail. Tsuji (1973) and Barklem & Collet (2016) have used other, independent data sources. Some obvious differences occur in particular when different molecular dissociation energies are assumed. We also find a few outliers, where one data source is off by many orders of magnitude with respect to all others. Although we cannot exclude human errors on our side, our comparisons have been made in a highly automated way, so these outliers are quite puzzling.

Tsuji (1973) and Sharp & Huebner (1990) recommend to use their fits only for T ≳ 1000 K, where the Sharp & Huebner fit-function behaves somewhat more robustly in extrapolations towards low temperatures. This is because is an approximately linear function of T as shown in Fig. 1. The fit term causes the extrapolation problems at low T which is avoided in the Stock (2008)-fits, which are therefore much more reliable for extrapolations towards low T. The old GGCHEM fits (Gail & Sedlmayr 1986) should not be applied for T≲500 K, and behave in a similar way as the Tsuji fits, namely with strong extrapolation errors towards low temperatures, because of a term.

Appendix C highlights a few more details of this comparison and Table D.1 lists our selected set of k_p(T) functions for 568 molecules composed of 24 elements, which can be safely applied to temperatures 100–6000 K. We could use more elements and molecules, but the data listed in Table D.1 have been carefully checked and benchmarked against the TEA code (Blecic et al. 2016), see Sect. 5.1.

Fig. 1 Comparison of kp(T) data for molecules HS (mercapto) and CN (cyanogen) used by several groups. The second column of smaller plots shows the deviations of log10kp from the mean values (computed without the Tsuji data). On the right side, the data has been converted to [kJ mol−1] (see Eq. (13)), before computing the deviations from the mean values. The large deviations at small temperatures for the Tsuji (1973) and Sharp & Huebner (1990) data are due to extrapolation errors beyond the valid fit-range. HS is classified as belonging to the group “data disagrees”, and CN to “data disagrees at low temperatures” in our comparison catalogue (Worters et al. 2017) since even without the Tsuji data, there are relatively large differences to the new (Barklem & Collet 2016) data.

3.2 Condensed phase thermo-chemical data

We have used two main sources for thermo-chemical data of condensed phases in this paper, the NIST-JANAF database (Chase et al. 1982; Chase 1986)², and the geophysical SUPCRTBL database (Zimmer et al. 2016; Johnson et al. 1992). When extracting Gibbs free energy data and/or vapour pressuredata p^vap(T) from these sources, it is important to understand the reference states of the elements with respect to which the data is presented, see more details in Appendix D. NIST-JANAF uses temperature-dependent reference states. When subtracting the data, the reference state Gibbs free energies cancel out, so we can use (18)

or, respectively, (19)

with p^⊖ = 1 bar. The data is available in the 7th column of the NIST-JANAF files at temperature points mostly separated by 100 K. The Gibbs free energy data can be directly compared to Sharp & Huebner (1990), and the vapour pressure data can be compared to, for example, Yaws (1999), Weast (1971) and Ackerman & Marley (2001).

The geophysical SUPCRTBL database provides condensed phase Gibbs free energy data via formulas with tabulated coefficients which partly have a direct physical meaning. In this database, the Gibbs free energy of formation is defined in a different way, namely with respect to the standard states of the elements. These standard states are the most stable forms of the pure elements at standard pressure p^⊖ = 1 bar and room temperature T_ref = 298.15 K. These are atomic gases for the noble gases He, Ne, Ar, …; diatomic gases , , , , for H, N, O, F, Cl; liquids for Br and Hg; and crystalline solids for all other elements. In order to use these data for our purpose, we need to subtract the Gibbs free energy of formation of the free atoms at temperature T from the standard states at reference temperature T_ref. The respective thermo-chemical data of the atoms are not available in SUPCRTBL, however we can use the NIST-JANAF database to provide this link, see further details in Appendix D. To fit the resulting functions, we have generated 100 data points at log-equidistant temperature points between 100 and 2500 K.

While fitting the data points from either the NIST-JANAF or the SUPCRTBL database by smooth functions, we have tried various options and carefully selected those fits which perform best concerning precision and possible interpolation and extrapolation issues, see Appendix D. If the corresponding molecule of a condensate is known to exist as a free molecule in the gas, it is generally advantageous to use p^vap(T) fits, because the values to be fitted only correspond to a few 100 kJ mol⁻¹. However, if no such molecule exists (which is true for most minerals) we need to fit the directly, with values of several 1000 kJ mol⁻¹, where the relative precision of the fits becomes crucial. A particular challenge is to properly fit pairs of solid and liquid phase data, which are quite similar to each other, yet their intersection point should match the melting point, see Fig. D.1. Table D.2 lists our fit formula and fit coefficients for 97 condensates, mostly fitted by us to the NIST-JANAF database, with a few additions from other sources. The fits done by ourselves generally have precisions < 1 kJ mol⁻¹, much better in most cases. The melting points are matched to within < 20 K (or <1%), much better in most cases, see D.3. Table D.4 lists the data for 160 condensates extracted from the SUPCRTBL database.

Figure 2 shows a comparison between -fits for MgSiO₃ (enstatite) obtained from three sources. The reader is welcome to study an auxiliary document (Woitke et al. 2017) which contains additional figures like Fig. 2 for all 121 condensates extracted from the SUPCRTBL database. The different data generally agree well with each other, within about 5 kJ mol⁻¹, although there are a few exceptions, remarkably CaAl₂Si₂O₈ (anorthite) which seems off by about 100 kJ mol⁻¹ in Sharp & Huebner (1990). The Sharp & Huebner fits can safely be extrapolated to about 500 K, but may have extrapolation artefacts below. All new fits to the NIST-JANAF data can be used from about 3500 K down to 100 K. This does require some extrapolations for a small number of condensates where no such NIST-JANAF data exist, but we have carefully checked that our fits continue smoothly and at least do not produce any artefacts. The behaviour of the various fit functions is similar to the molecular k_p -fits as described in Sect. 3.1, the Stock-fit function is most robust towards low temperatures, but slightly less accurate at high temperatures.

Fig. 2 Comparison of condensed phase Gibbs free energy data from Sharp & Huebner (1990, orange), fits to the NIST-JANAF database done by the authors of this paper for GGCHEM (green), and fits to the SUPCRTBL data (Zimmer et al. 2016, blue). The SUPCRTBL data points have been generated for temperatures 100 to 2500 K according to the equations and parameters given in Zimmer et al. (black dashed curve). Lower plot: deviations [kJ mol−1] with respectto the mean of all three data sources.

4 The GGchem code

GGCHEM is an abbreviation for the German word “Gleich-Gewichts-Chemie” which means equilibrium chemistry. The code has been originally developed by H.-P. Gail in the 1970s (e.g. Gail & Sedlmayr 1986). Further works on the code have been done by C. Dominik, Ch. Helling and P. Woitke in Berlin until 2005. For this paper, we have completely re-written the code and given it a modern FORTRAN-90 code architecture. We have updated the thermo-chemical input data (Sect. 3.1), and have added a quadruple precision version and a new pre-iteration scheme (Appendix A) to improve code stability and range of applicability down to 100 K. We have added the equilibrium condensation part (Appendix B), and have added and revised the condensed phase thermo-chemical data (Sect. 3.2). The code is publicly available including all thermo-chemical data³.

GGCHEM allows the user to select elements, molecules and condensates via versatile input files. All elements from hydrogen to zirconium (atomic number 40) are supported, with additional options to include tungsten (atomic number 74) and charges. The element abundances are currently pre-compiled with options for solar abundances (Asplund et al. 2009), meteorite abundances or Earth crust abundances. The default choice is to select all molecules and condensates which are composed of the selected elements from provided data-files, but the user is welcome to develop their own input files, for example a customised small model optimised for speed, with just a few elements and hand-selected species. There are options to just solve the pure gas phase chemical equilibrium or the combined problem with equilibrium condensation. The gas mass density ρ or the total gas pressure p can be chosen as input parameter, besides the temperature T.

Fig. 3 Benchmark test against TEA (Blecic et al. 2016) at constant p = 1 bar and T = 100–6000 K, showing mixing ratios ni∕n. The coloured full and dashed curves are the GGCHEM results. The TEA results are overplotted with narrow grey lines. Only the most important molecules are shown, that is, those which reach a certain threshold concentration, depending on element. The TEA results stop at 400 K as TEA has difficulties converging at lower temperatures, even using 2000 iterations. The following molecules with notable concentrations are missing in TEA: SiH2, SiH3; CaH; TiS, TiH, TiN, TiC; CrH, CrS; MnH, MnS, MnCl, MnF, MnO; FeH; NiH, NiF, NiO, and the following molecules have been omitted by us in the TEA model to improve its convergence: Si(CH3)4, SiCH3 Cl3 and Ni(CO)4. Differences in the selection and availability of molecular species explain all visible deviations between GGCHEM and TEA.

5 Results

5.1 TEA benchmark test

To benchmark our results, we have used the public thermo-chemical equilibrium abundances (TEA) code (Blecic et al. 2016)⁴. TEA followsthe Gibbs free energy minimisation method as described by White et al. (1958) and Eriksson (1971), using a Lagrangian optimisation scheme with Lambda correction algorithm. The Gibbs free energies of the molecules are computed from the thermo-chemical data in the NIST-JANAF tables (Chase 1998), for most molecules given between 100 and 6000 K, separated between 100 K, provided by Dr. Thomas C. Allison in October 2012 through private communication. The program uses internal spline interpolations to compute the Gibbs free energy between these temperature points. The NIST-JANAF data is available for 600 gaseous molecular species (84 elements).

TEA has been developed in particular to determine the abundances of the major, most-abundant and spectroscopically most active gaseous species expected to be present in hot, giant exoplanetary atmospheres, as these species have a dominant influence on the planetary spectra in the optical and infrared (e.g. Seager & Deming 2010; Burrows 2014). TEA has been benchmarked against CEA (Chemical Equilibrium with Applications), see Gordon & McBride (1994) and McBride & Gordon (1996), as well as against several analytical codes (Burrows & Sharp 1999; Heng & Tsai 2016; Heng & Lyons 2016). These tests typically involve a few elements with a few dozens of molecules between 500 and 4000 K, where TEA is found to converge robustly within reasonable times (2–3 CPU-sec per (p, T)-point).

For our benchmark test between TEA and GGCHEM we have selected the following 24 elements of astrophysical interest: H, He, Li, C, N, O, F, Na, Mg, Al, Si, P, S, Cl, K, Ca, Ti, V, Cr, Mn, Fe, Ni, Zr and W with solar abundances (Asplund et al. 2009). The task was to compute the chemical composition of the gas at p = 1 bar from T = 6000 K down to 100 K (or as low as possible), with all available molecules.

Since the TEA-code can presently only handle neutral gas species, we have disregarded all charged species and switched off equilibrium condensation in GGCHEM for this test. The results are shown in Fig. 3. For practically all depicted molecules, the GGCHEM results (coloured lines) agree very well with the TEA results, to a precision better than the linewidth of about ~ 0.05 dex. All visible deviations can be explained by some molecules not available in NIST-JANAF, and are hence missing in TEA, or some very large molecules that we have de-selected on purpose in order to improve the convergence of TEA. The missing molecules are listed in the caption of Fig 3. GGCHEM uses 421 molecular species plus 24 neutral atoms for this test, whereas TEA uses altogether 400 species including the 24 neutral atoms. The agreement between TEA and GGCHEM is very reassuring, however, TEA failed to converge within 2000 iterations for temperatures below 400 K.

GGCHEM needs about 0.004 CPU-sec per point on a 2.8 GHz Linux Laptop. The time consumption is mainly caused by the internal matrix operations which scale as K³ where K is the number of elements. For the lower temperatures in this benchmark test, GGCHEM uses quadruple precision arithmetics. For temperatures >1000 K, the code switches to double precision which takes only 0.0013 CPU-sec per point for 24 elements. To compute the 100 temperature points requested for this benchmark test, the total computational time consumption of GGCHEM is 0.71 CPU-sec, where about half of it is initialisation. This is in sharp contrast to TEA which needs about 45 CPU-min per point with 2000 iterations, resulting in altogether 74 CPU-hours to complete this benchmark test. Clearly, we have pushed TEA to its limits in this benchmark test, TEA should normally be run with about 100 iterations where it performs as stated above. However, using TEA with 2000 iterations in this benchmark test to obtain results for as low as possible temperatures, we conclude that GGCHEM is more than 10⁵ times faster than TEA.

Fig. 3 continued.

Fig. 4 Concentrations of major molecules at 1 bar according to a simplified nine-molecule model in chemical equilibrium (full lines), in comparison to a model with 24 elements and 445 molecules (dashed) and a full model with equilibrium condensation (dotted).

5.2 Abundance of major molecules

In a first application, we have used the simplified setup of Heng & Tsai (2016) with 4 elements (H, C, N, O) and 9 molecules (H₂, H₂ O, CO, CO₂, N₂, NH₃, HCN, C₂ H₂, C₂ H₄) with solar abundances (Asplund et al. 2009). The results shown in Fig. 4 can be directly compared to the upper part of Fig. 1 in Heng & Tsai (2016). GGCHEM needs about 0.12 milli CPU-sec per point in this setup.

Figure 4 compares these results to two other GGCHEM models. The dashed lines show the results from the full model with 24 elements and 445 molecules as specified in Sect. 5.1, and the dotted lines show the results from a full model with equilibrium condensation as specified in Sect. 5.3. The comparison shows that, although the simplified nine-molecule model results in very reasonable approximations for the concentrations of the major molecules between 500 and 3000 K, there are some notable deviations around 1000–2000 K which are likely to be relevant for the retrieval of atmospheric properties from spectroscopic observations. These deviations are due to the consumption of oxygen by (i) the formation of SiO, Mg(OH)₂ and Fe(OH)₂ molecules, and (ii) silicate and phyllosilicate condensation. Both effects reduce the concentrations of H₂ O and CO₂, and increase C₂ H₂ and C₂ H₄ concentrations, by up to 0.5 dex around 1500 K. This effect is discussed further in Sect. 5.6.

5.3 The condensation of the elements

In the following, we study the results of a model with equilibrium condensation switched on, at constant pressure p = 1 bar, for 22 elements (H, He, Li, C, N, O, Na, Mg, Al, Si, S, Cl, K, Ca, Ti, V, Cr, Mn, Fe, Ni, Zr, W) and charged species, with solar abundances. We have omitted fluorine and phosphorus in this model because of our possibly incomplete data for fluorine and phosphorus condensates. Altogether 388 gas phase and 190 condensed species are taken into account in this model. Gas phase species include 365 molecules, 22 free atoms and the free electron. Condensed species include 38 liquids (Table D.3). The thermo-chemical data of the condensed species are mainly taken from the SUPCRTBL database (Table D.4), completed/replaced by our own fits to NIST-JANAF data of 92 condensed species (Table D.2), plus five further condensed species (CaTiO₃[s], MnS[s], NiS[s], NiS₂[s], Ni₃ S₂[s]) taken from Sharp & Huebner, which are neither available in NIST-JANAF nor in SUPCRTBL.

Figure 5 shows the dust to gas mass ratio as function of temperature in phase equilibrium, calculated as (20)

where j[cond] is the index of a condensed species, n_j[cond] its particle density (i.e. number of condensed units per volume) and m_j[cond] its mass. Tungsten establishes a dust to gas ratio of ~10⁻⁹ at T≲2250 K, zirconium dioxide brings it to ~10^−7.5 at 2000 K, and aluminium increases it to ~10⁻⁴ at 1900 K. Subsequently, calcium, iron, silicon and magnesium compounds build up a dust to gas ratio of ~ 10^−2.38 ≈ 0.004 at 1500 K, which then remains about constant towards much lower temperatures. The condensation of titanium, nickel, vanadium, chromium, manganese, sodium and potassium compounds make only minor contributions to dust/gas(T). The value of dust∕gas ≈ 0.004 is significantly lower than the standard value of dust∕gas ≈ 1∕100 (e.g. Beckwith et al. 1990; Chiang & Goldreich 1997; D’Alessio et al. 1998; Dullemond et al. 2002). At around 650 K, sulphur starts to condense which increases dust/gas to about 0.0045 with further minor contributions by chlorine and lithium. Just below 500 K, several minerals start to become hydrated by the incorporation of either water or hydroxyl groups, forming so-called phyllosilicates. This additional intake of oxygen atoms increases the dust to gas ratio to about 0.0052. Around 240 and 140 K, water ice and ammonia condense, respectively, which finally brings the dust to gas ratio to a value of almost exactly 1/100.

The lower part of Fig. 5 shows the depletion of the elements from the gas phase due to condensation. At 1000 K, for example, the gas consists mainly of S and Cl besides H, He, C, N and O, whereas all other elements, which are usually more abundant like Fe, Si and Mg, are already depleted by more than 8 orders of magnitude. The intake of oxygen in silicates and phyllosilicates leads to a substantial increase of the C/O ratio in the gas phase, from C∕O = 0.55 at T = 2500 K to C∕O = 0.71 at T = 1000 K (silicates) and C∕O = 0.83 at T = 300 K (phyllosilicates). This increase of C/O may well have observational consequences in brown dwarf and giant gas planet spectra. Once water condenses, the C/O ratio becomes very large, for example C∕O > 10⁶ at 150 K, and the gas consists mainly of CH₄ and NH₃ (beside H₂ and noble gases), similar to the atmosphere of Titan (Fulchignoni et al. 2005)⁵. Remarkably, at p = 1 bar, not a single carbon atom has been incorporated into any condensed state in phase equilibrium down to 100 K for solar abundances.

The onset of condensation: Figure 6 visualises the first stages of condensation in phase equilibrium at p = 1 bar. W[s] (crystalline tungsten) is found to be the first stable condensate at 2216 K, followed by ZrO₂[s] (baddeleyite) at 2027 K. Further phase transitions are as follows: Al-gas → Al₂ O₃[s] (corundum) at 1957 K, Ti-gas → CaTiO₃[s] (perovskite) at 1913 K, Ca-gas → Ca₂ Al₂SiO₇[s] (gehlenite) at 1880 K, Fe-gas → Fe[l] (liquid iron) at 1841 K. These are all type-1 (gas → condensed) phase transitions, featured by a smooth, rounded-off increase of the concentration of the new condensate as function of temperature. This gradual increase is mirrored by a gradual decrease of one element abundance in the gas phase, an element which was not affected by condensation before, to keep the partial pressure equal to the vapour pressure p^vap(T) of the new condensate.

At 1820 and 1777 K, the first type-2 transitions (condensed → condensed) occur, namely Fe[l] → Fe[s] (solid iron) and Al₂O₃[s] → MgAl₂O₄[s] (spinel), respectively. These transitions transform one combination of condensates into another, without substantial changes of any gas phase abundances. Type-2 transitions occur suddenly, they do not change the number of elements affected by condensation, hence the number of independent variables in the equation system to be solved, and thus the total number of present condensates (see Appendix B). Type-2 transitions have a box-like appearance in Fig. 6.

The appearance of the major condensates: The condensation of the elements continues in Fig. 6 with a number of further type-1 transitions: Si-gas → SiO[s] (silicon monoxide) at 1729 K at 1 bar, V-gas → VO[s] (vanadium monoxide) at 1710 K, Ni-gas → Ni[s] (crystalline nickel) at 1690 K, Mg-gas → Mg₂ SiO₄[s] (fosterite) at 1661 K and Cr-gas → Cr[s] (crystalline chromium) at 1513 K. At this stage, the build-up of the dust to gas ratio to a value of about 0.004 is mostly completed and there are 11 elements which are strongly reduced by condensation (see lower part of Fig. 5), which sets the number of simultaneously present condensates to 11 (as explained in Appendix B).

Increasing complexity 1: A number of further type-2 transitions occur in Fig. 6 towards lower temperatures, namely Ca₂ Al₂SiO₇[s] → Ca₂ MgSi₂O₇[s] (akermanite) at 1685 K, MgAl₂O₄[s] → CaAl₂Si₂O₈[s] (anorthite) at 1440 K, CaTiO₃[s] → Ti₄ O₇[s] (titanium oxide) at 1420 K, Ca₂MgSi₂O₇[s] → CaMgSi₂O₆[s] (diopside) at 1420 K, Cr[s] → MgCr₂O₄[s] (picrochromite) at 1336 K, ZrO₂[s] → ZrSiO₄[s] (zircon) at 1334 K, SiO[s] → MgSiO₃[s] (enstatite) at 1292 K, Ti₄O₇[s] → MnTiO₃[s] (pyrophanite) at 1211 K, and VO[s] → V₂ O₃[s] (karelianite) at 1194 K. These type-2 transitions can have some side-effects on other more abundant condensates, so strictly speaking, each type-2 phase transition affects a combination of solids where one is formed, one is consumed and a number of other, more abundant condensates adjust their concentrations. This applies in particular to MgSiO₃[s] and Mg₂ SiO₄[s] which are the main reservoir for condensed silicon and magnesium, which are affected by most type-2 transitions.

On the left side of Fig. 6, a few more type-1 transitions take place, namely Mn-gas → MnS[s] (alabandite) at T = 1290 K, Na-gas → NaAlSi₃O₈[s] (albite) at T = 1231 K and K-gas → KAlSi₃O₈[s] (microcline) at T = 1154 K, increasing the number of simultaneously present condensates to 14. At 1000 K, all three major components of feldspar are present (KAlSi₃O₈[s]– NaAlSi₃O₈[s]– CaAl₂Si₂O₈[s]) which is the most common mineral on Earth making up about 41% of its continental crust (Anderson & Anderson 2010). Feldspar is the major constituent of basalt, which is a key component of oceanic crust as well as the principal volcanic rock on Earth.

Increasing complexity 2: Figure 7 continues to show the model results down to 460 K. Most prominently, we see three more type-1 transitions: S-gas → FeS[s] (troilite) at 678 K, Cl-gas → NaCl[s] (halite) at 613 K and Li-gas → LiCl[s] (lithium-chloride) at 573 K, increasing the number of simultaneously present condensates to 17. Further type-2 transitions occurs as well, in particular MnTiO₃[s] → CaTiSiO₅[s] (sphene) at 814 K, MnS[s] → Mn₃Al₂Si₃O₁₂[s] (spessartine) at 670 K, KAlSi₃O₈[s] → KMg₃AlSi₃O₁₂H₂[s] (phlogopite) at 520 K, CaAl₂Si₂O₈[s] → MgAl₂O₄[s] (spinel) at 514 K, NaAlSi₃O₈[s] → NaMg₃AlSi₃O₁₂H₂[s] (sodaphlogopite) at 502 K, and CaTiSiO₅[s] → FeTiO₃[s] (ilmenite) at 486 K. KMg₃AlSi₃O₁₂H₂ and NaMg₃AlSi₃O₁₂H₂ are the first two phyllosilicates, with additional intake of hydrogen and oxygen.

The formation of phyllosilicates: Figure 8 shows the results down to 100 K. Again, a number of complex type-2 transformations occur, where most of the condensates are replaced by phyllosilicates, step by step. The most significant phyllosilicate is Mg₃ Si₂O₉H₄[s] (lizardite) which replaces Mg₂ SiO₄[s] at 345 K. At 270 K, the main condensates are Mg₃Si₂O₉H₄[s], FeS[s], Fe₃O₄[s] (magnetite), NaMg₃AlSi₃O₁₂H₂[s], Fe₃ Al₂Si₃O₁₂[s] (almandine), Ca₃ Al₂Si₃O₁₂[s] (grossular) and Ni₃ S₂[s] (heazlewoodite), whereas the more simple condensates commonly known from higher temperatures (fosterite, enstatite, solid Fe, etc.) are all gone.

The condensation of water and ammonia: The condensation of the elements is completed by two major type-1 transitions: O-gas → H₂ O[s] (water ice) at 247 K and N-gas → NH₃[s] (ammonia ice) at 141 K. We do not observe the condensation of CO[s] or CO₂[s] (dry ice) in our equilibrium models with solar abundances, although at 1 bar, CO₂ is well-known to deposit to a solid at − 79 °C = 195 K on Earth. This is because the Earth atmosphere is oxygen-rich and relatively hydrogen-poor. For solar abundances, all carbon is locked up in gaseous CH₄. At 100 K, all elements except hydrogen, carbon and noble gases have condensed into solid phases, and the composition of the gas phase is extremely simple – just H₂, He, CH₄, and noble gases. The next major phase transition to occur would be C-gas → CH₄[s] (methane ice), but this transition occurs below 100 K and our models do not reach these temperatures. We note that we are currently lacking organic liquids/solids in our collection of condensed species, which might change some of these conclusions.

Fig. 5 Upper part: dust to gas mass ratio in phase equilibrium as function of temperature at constant p = 1 bar for solar element abundances. The arrows indicate where the elements start forming condensed phases in significant amounts. The elements marked in red have a significant influence on the shape of the dust/gas(T)-curve. Lower part: remaining element abundances in the gas phase.

Fig. 6 Onset of condensation at p = 1 bar for solar abundances in phase equilibrium. The graph show the concentration of the various condensed species with respect to hydrogen nuclei.

Fig. 7 Concentration of condensed species at medium temperatures 460–1000 K at p = 1 bar for solar abundances. Please note the changed scaling of the y-axis in comparison to Fig. 6.

Fig. 8 Concentration of condensed species at temperatures 100–460 K at p = 1 bar for solar abundances.

Fig. 9 Impact of CO-blocking on the composition of the gas at constant pressure and temperature. In the upper plot, the full lines represent the results from a pure gas phase model without condensation (24 elements and 445 molecules, see Sect. 5.1), and the dashed lines show the results from the model including equilibrium condensation (Sect. 5.3). ϵO is a fixed value here, whereas ϵC varies.

5.4 CO-blocking

Figure 9 demonstrates the well-known impact of the carbon-to-oxygen ratio C/O = ϵ_C∕ϵ_O on the composition of the gas phase, and on the condensation. The figure is designed to match the results of Mollière et al. (2015), although here (in contrast to Mollière et al.) we are plotting particle concentrations. At temperatures between about 850 and 4000 K at 10 mbar, CO is the dominant gas species in the carbon-oxygen system with a dissociation energy of about 11.1 eV. Therefore, CO almost completely consumes all carbon when C∕O≲1 or all oxygen when C∕O ≳ 1, an effect known as CO-blocking (see e.g. Sedlmayr & Dominik 1995; Sedlmayr 1997; Beck et al. 1992). The effect is particularly strong on hydro-carbon molecules like C₂ H₂ (acetylene) which are suppressed by more than 6 orders of magnitude if C∕O < 1.

Closer inspection of our results (Fig. 9) shows that the transition actually takes place around C∕O ≈ 0.96, which agrees well with the Mollière et al. results, who have been using the CEA code (Gordon & McBride 1994; McBride & Gordon 1996). Figure 1 in Madhusudhan (2012) actually shows the same effect, although Madhusudhan has used a log-scaling of the C/O-axis where deviations of order 5% are simply not visible. The reason for this asymmetry is the formation of the SiO molecule with a dissociation energy of 8.2 eV. If C∕O≲1, the small amount of excess oxygen is first of all consumed to form SiO molecules, with similar consequences for the abundances of water and hydro-carbon molecules as used from the CO-blocking. We measure numerically that the abundances of H₂ O and CH₄ intersect at C∕O = 1.00 when ϵ_Si → 0, at 0.98 when ϵ_Si = 7.3, at 0.96 when ϵ_Si = 7.51 (solar), at 0.92 when ϵ_Si = 7.7, and at 0.82 when ϵ_Si = 8.

Figure 9 also shows the influence of the CO-blocking on the condensation. At 1500 K and 10 mbar, C[s] (graphite) can condense at C∕O ≳ 1.07, SiC[s] (silicon carbide) at C∕O ≳ 0.96 and TiC[s] (titanium carbide) at C∕O ≳ 0.93. At smaller C/O, the oxygen-containing solids can form, here in particular Al₂ O₃[s] (corundum), Ca₂ Al₂SiO₇[s] (gehlenite), MgAl₂O₄[s] (spinel), Ca₂ MgSi₂O₇[s] (akermanite), CaTiO₃[s] (perovskite) and Mg₂ SiO₄ (fosterite). Other condensates like W[s] and ZrO₂[s] (not shown) and Fe[s] are not affected by C/O, because they either do not contain oxygen or carbon, or have even higher dissociation energies than CO, like ZrO₂[s].

The feedback of condensation on the molecular abundances for C∕O > 1 is more substantial than for C∕O < 1. The molecular concentrations of C₂ H₂ and CH₄ are found to not increase any further once carbon dust is present, as any surplus carbon is simply converted into additional carbon dust, but not into additional hydro-carbon molecules in phase equilibrium. We have actually used the total (gas + dust) abundances of oxygen and carbon to create Fig. 9, because the gas-phase C/O only increases to about 1.1 even if much larger values of the total C/O are used.

At lower temperatures (T < 850 K at p = 10 mbar), however, CH₄ becomes the dominating carbon gas phase species, whereas CO becomes a trace species, and so the spell is broken. In a hydrogen-rich gas, CH₄ is thermodynamically more favourable than any solid carbon compound at low temperatures, which leads to high water concentrations independent of C/O, and no carbon dust in phase equilibrium. This confirms our earlier findings in kinetic cloud formation models for carbon-rich atmospheres (Fig. 6 in Helling et al. 2017).

5.5 Phase diagrams of the elements

Figure 10 shows some results of the equilibrium condensation model (22 elements, charges, 365 molecules and 190 condensed species, see Sect. 5.3) in the entire (p, T)-plane between 10⁻⁹ and 100 bar, and between 100 and 2500 K, for solar abundances. The figures indicate the most abundant species that contain certain elements, after multiplication by their stoichiometric factor, hence the main reservoirs of the elements as function of pressure and temperature.

We generally see increasing complexity from high to low temperature, from ions over atoms to simple molecules, to complex molecules, and then to condensates with increasingly complex stoichiometry and phyllosilicates at the lowest temperatures. In the first approximation, all molecular transitions, as well as the type-1 and type-2 phase transitions discussed in Sect. 5.3 simply occur at lower temperatures when the pressure is decreased, roughly a factor of 2–4 in temperature for 10 decades of pressure decrease. This leads to numerous slanted and slightly curved phase transition lines, similar to phase diagrams in thermodynamics. Most prominent molecular transitions are H → H₂ (at ~2200 K), CO → CH₄ (at ~650 K) and N₂ → NH₃ (at ~320 K), where the numbers in brackets refer to a pressure of p = 10⁻⁴ bar.

The details, however, are in fact much more complicated, in particular concerning the gas → solid transitions. The order ofcertain phase transitions as described in Sect. 5.3 may be different, and some condensates may be present only in certain parts of the (p, T)-plane. The main silicates can form below about 1700 to 1100 K, and eventually the elements Mg, Si, Na, Al, K and Mn can form thermodynamically stable phyllosilicates at temperatures below about 600 to 200 K.

We want to stress, however, that all these results depend on the chosen element abundances, as demonstrated for the case of varying C/O in Sect. 5.4. The results do also depend technically on the choice of elements, so if another element is taken into account in addition, it will interfere with all others which can potentially change the results shown in Fig. 10. Therefore, our “phase diagrams” are not as general as their name suggests. The results should not be taken out of context and be applied to, for example, the atmospheres of rocky planets. The proper way to do this is to run the GGCHEM-code for the expected pressures, temperatures and local element abundances in the atmosphere.

Fig. 10 Gas phase/condensed species which contain most of the elements in the (p, T)-plane. Solid species are marked by “[s]” and liquid species by “[l]”. At low temperatures, the elements Mg, Si, Na, Al and K are mostly present in form of phyllosilicates. The only further element found to form phyllosilicates in phase equilibrium is manganese.

Fig. 10 continued.

5.6 The C/O-ratio as affected by metal oxide molecules and condensation

Figure 11 (right side) shows the C/O ratio in the gas phase as affected by condensation in the (p, T)-plane. One of the most significant results of this paper again shows, namely that the solar value of about C/O = 0.55 only holds prior to the condensation of the main silicates. The consumption of oxygen by silicate formation is substantial, and increases the C/O ratio to about 0.71. Once the silicates become hydrated to form phyllosilicates, which happens to the most abundant silicates between about 200 to 400 K in the equilibrium model (depending on pressure), the C/O ratio jumps to 0.83, and then gets very large once water condenses⁶. The little blue area in the central top of Fig. 11 is not a computational error – it is a small region in the (p, T)-plane where we find solid carbon C[s] (graphite) to condense.

To clarify the effect of condensation on gas phase C/O, we compare our results to a pure gas phase model (without condensation) where we compute an “effective” C/O as (21)

where only atoms and molecules are counted here if they are made of the elements H, C, N and O, but we set the stoichiometric coefficient if a molecule i contains any other (metal) elements like Mg, Si or Fe. In a pure gas phase model, as plotted on the left side of Fig. 11, the proper C∕O ≈ 0.55 is of course constant everywhere in the (p, T)-plane by assumption, but will be different because of the oxygen and carbon that is bound in molecules which are neglected in Eq. (21). We do this because in fast spectral retrieval models (e.g. Lavie et al. 2017), the influence of these metals on the concentration of molecules like H₂ O and CO₂ is often neglected, that is, is assumed.

In Fig. 11 (left side) we see that may become substantially larger than C/O at lower temperatures. increases from the solar value 0.55 to about 0.58 due to the formation of SiO between about 2000 and 3000 K, and then increases further to values above 0.8 once Mg(OH)₂ and Fe(OH)₂ become abundant in the gas phase, which happens at temperatures between about 600 and 1100 K.

In both discussed cases, the large pure gas phase model and the equilibrium condensation model, the reason for the increase of the (effective) C/O ratio is the same – the consumption of oxygen atoms by the formation of metal compounds, either oxide and hydroxide molecules or condensates, which are then missing for the formation of gaseous molecules like H₂ O and CO₂. The effect is actually similar in both models, because each Mg, Si and Fe atom, on average, consumes about 1–2 oxygen atoms in both cases at low temperatures. The main difference is that in models including condensation, the establishment of large C/O ratios already occurs at higher temperatures, between about 1200 to 1600 K. We discuss this point further in Sect. 7.

Figure 12 shows the abundance of the three most significant Mg-silicates in the (p, T)-plane. As described in Sect. 5.3, the two major reservoirs for condensed Mg and Si at relatively high temperatures are MgSiO₃[s] (enstatite) and Mg₂ SiO₄[s] (fosterite). However, some more stable silicates exist, which contain Ca, Fe, Ti, Mn or Cr, for example. These can form in smaller amounts only, due to element abundance constraints, so they “steal” some Mg and Si from enstatite and fosterite, which affects the concentrations of enstatite and fosterite in complicated ways. Fosterite generally becomes stable at slightly higher temperatures than enstatite, but once formed, enstatite can become the most abundant condensate in certain parts of the (p, T)-plane, along with Fe[s], Fe₂ SiO₄[s], FeS[s] or Fe₂ O₃[s].

Below 200–400 K (depending on pressure) both Mg-silicates are combined and hydrated to form the phyllosilicate Mg₃ Si₂O₉H₄[s] (lizardite). Interestingly, there is no feedback of water condensation on this phyllosilicate, the concentration of lizardite stays constant towards very low temperatures, showing that the formation of lizardite is thermodynamically favoured over the formation ofwater ice.

We conclude from Fig. 10 that Mg, Si, Na, Al, K and Mn can form thermodynamically stable phyllosilicates at low temperatures. Beside the major phyllosilicate Mg₃ Si₂O₉H₄[s] (lizardite), we also find NaMg₃AlSi₃O₁₂H₂[s] (sodaphlogopite), CaAl₂Si₂O₁₀H₄[s] (lawsonite), Mg₃ Si₄O₁₂H₂[s] (talc), KMg₃AlSi₃O₁₂H₂[s] (phlogopite), KFe₃AlSi₃O₁₂H₂[s] (annite), KMn₃AlSi₃O₁₂H₂[s] (Mn-biotite) and MnAl₂SiO₇H₂[s] (Mn-chloritoid) is less amounts, due to element abundance constraints.

Fig. 11 C/O ratio in the gas phase affected by the formation of metal oxide and hydroxide molecules (left) and condensation of silicates, phyllosilicates, and water ice (right). On the left side, we plot the “effective” C/O that is contained in the gas in form of molecules made of H–C–N–O only, but do not count metal oxide andhydroxide molecules. The contour lines shown are 0.58 (after SiO formation) and 0.8 (after Mg(OH)2 and Fe(OH)2 formation). On the right side, we plot the resulting gas phase abundances after the condensation of solids and liquids. The increase of C∕O > 0.7 is due to the consumption of oxygen by the major silicates like MgSiO3[s] and Mg2SiO4[s], and the further increase to >0.8 is due to phyllosilicates. The sharp increase at low temperatures (red area) is due to water condensation.

Fig. 12 Concentration (xj = nj∕n⟨H⟩) of the three major magnesium-silicates in the (p, T)-plane in the equilibrium condensation model.

Fig. 13 Left: nucleation rate of tungsten per hydrogen nucleus J⋆∕n⟨H⟩ [s−1] as function of gas pressure and temperature. Centre: size of the critical cluster N⋆. Right: nucleation timescale τnuc [yr] according to Eq. (25). The white dashed contour line marks saturation S = 1, the gas is supersaturated to the left of this line.

6 The nucleation of tungsten

An interesting side-result of our models is that solid tungsten W[s] is found to be the first condensate, that is, the first solid material that becomes thermodynamically stable in any galactic/stellar/atmospheric cooling flow which is initially too hot to contain any condensed materials. The question arises whether tungsten could provide the first nucleation seeds required for the growth of all other, more abundant solid/liquid materials which only become stable at lower temperatures. See Gail & Sedlmayr (1998) for a more detailed introduction to the search for the first condensate.

6.1 Is there enough tungsten?

To estimate the total number of seed particles that can be provided by tungsten, let us assume that one seed particle consists of N_l ≈ 10–1000 tungsten atoms. Using the solar abundance of tungsten ϵ_W = 7.1 × 10⁻¹² (Asplund et al. 2009), the concentration of tungsten seed particles cannot exceed (22)

Next, we assume that the various silicates, iron and all other abundant condensates will grow on top of these seed particles later, to form a solid mantle of unique thickness. The total available volume of solid material per H nucleus is about (23)

where we have used solar abundances and where V₁ ≈ (1.2–3) × 10⁻²³ cm⁻³ is the dust volume per heavy atom (the sum of Si, Mg and Fe atoms in the condensate), estimated from the mass and solid densities of iron, MgSiO₃ and Mg₂ SiO₄. The radius of the dust particles forming this way would be given by , that is, (24)

which falls surprisingly well into the size range of observed dust particles in space. This proves nothing of course, but it shows that there is enough tungsten in a gas of solar abundances to provide the seed particles needed for astrophysical bulk condensation. In fact, from the simple estimate presented here we can conclude that the concentration of seed particles must not exceed about 10⁻¹⁵ to 10⁻¹³ per hydrogen nucleus, otherwise too many, too small dust particles would form by condensation.

In application to M-type AGB-star winds, Höfner et al. (2016) have shown that only silicate grains in a certain size range 0.1 to 1 μm can efficiently drive an outflow by radiation pressure, via their large scattering opacity around λ ~ 1 μm. In contrast, if their absorption coefficients dominate at those wavelengths (as is true when they contain some iron), the grains rather get heated in the strong IR radiation field and evaporate (Woitke 2006). In order to obtain these particular particle sizes, Höfner et al. (2016) needed to assume a certain seed particle concentration in their time-dependent hydrodynamical and dust formation models, n_seeds∕n_⟨H⟩≈ 3 × 10⁻¹⁴ to 3 × 10⁻¹⁵.

6.2 Timescales and supercooling

The nucleation rate of crystalline tungsten has been calculated according to classical nucleation theory in the formulation of Gail et al. (1984), with surface tension σ = 3340 erg cm⁻² (Tran et al. 2016) and parameter N_1∕2 = 10 (the cluster size where the surface tension reduces to half its bulk value). The supersaturation ratio S and the particle density of atomic tungsten n_W are calculated by a pure gas-phase GGCHEM model without condensation. The resulting nucleation rates per hydrogen nucleus are shown in Fig. 13, from which we estimate that a supercooling of about 300–400 K is needed to start the effective nucleation of tungsten. The time required to form these seed particles is (25)

From Fig. 13 we read off that the formation timescales of the tungsten seed particles would be rather long and would require high gas pressures, for example τ_nuc = 1 yr at 0.1 bar and 1600 K. At lower gas pressures, the tungsten nucleation timescale increases fast, exceeding the age of the universe around 10⁻⁶ bar. This fast increase is mainly due to an increase of the critical cluster size N_⋆ (central part of Fig. 13).

7 Summary and conclusions

A fast and versatile public code has been developed, called GGCHEM, to calculate the chemical composition of astrophysical gases in chemical equilibrium down to 100 K, with or without equilibrium condensation. We have collected and compared the thermo-chemical data for molecules and condensates from different sources, and have presented tables with temperature-fits of molecular equilibrium constants k_p (T) and Gibbs free energies which behave robustly towards low temperatures and can safely be applied down to 100 K. The data comprises 552 molecules and 257 condensates, including 38 liquids.

For pure gas-phase applications, GGCHEM has been thoroughly tested against the TEA-code (Blecic et al. 2016), revealing close to perfect agreement. Differences arise only when the selection of molecules differs. In GGCHEM, we are using a number of molecules with thermo-chemical data from Barklem & Collet (2016), which are not available in NIST-JANAF, in particular metal hydrides like CaH, TiH and FeH. GGCHEM is found to converge fast and robustly down to 100 K under all tested circumstances including very small and very large selections of elements and molecules, as well as for unusual element abundances.

Concerning very sparse models with just a few molecules and elements (for example H, C, N, O to quickly determine the abundances of the spectroscopically active molecules like H₂ O and CH₄), we want to emphasise that the combined effect of elements like Mg, Si and Fe is substantial and might lead to an overestimation of the C/O ratio if these elements are neglected. For solar abundances (C∕O = 0.55) the combined abundance of Si, Mg and Fe alone is about 10⁻⁴ which is larger than the nitrogen abundance (ϵ_N = 6.8 × 10⁻⁵) and can be compared to the carbon and oxygen abundances (ϵ_C = 2.7 × 10⁻⁴, ϵ_O = 4.9 × 10⁻⁴). Since molecules like SiO, and form at low temperatures, every Si, Mg and Fe atom takes away about x ≈ 1–2 oxygen atoms from the C–N–O system. If we consider the effective C/O ratio remaining in the C–N–O system, the result for solar element abundances is (26)

Thus, if C/O ratios are derived from observations based on sparse chemical models, where elements like Mg, Si and Fe are lacking, the results must be expected to be too large.

We have discussed the condensation sequence of the elements in phase equilibrium, generally confirming the results found earlier by Sharp & Huebner (1990) and Lodders (2003). Differences arise mostly from the different selection/availability of the thermo-chemical data for condensed species, but we did not find any evidence for systematic differences due to different numerical techniques. An important step here was to create a link to the geophysical database SUPCRTBL (Zimmer et al. 2016), from which we have extracted the thermo-chemical data for 121 condensed species, including phyllosilicates. The numerical methods rather affect the computational speed and stability of the codes.

A straightforward result from our models is that the dust to gas mass ratio in a solar composition gas is expected to be 0.0045 rather than 0.01. The latter (standard) value of 0.01 is only obtained in our models once water and ammonia ices condense at low temperatures. The ice formation increases the mass of the condensates by a factor of ~2.5 and the volume of the condensates by a factor of ~6 (from 2.7 × 10⁻²⁷ to 1.7 × 10⁻²⁶ cm³ per hydrogen nucleus to be precise).

Our models show that phyllosilicates (“wet silicates”) are the most abundant condensates at low temperatures in a solar composition gas in phase equilibrium. The most abundant phyllosilicate is identified to be Mg₃ Si₂O₉H₄[s] (lizardite) which contains the majority of Mg and Si below 200–400 K in phase equilibrium, depending on pressure. In a similar way, a few other phyllosilicates are found to consume most of Na, Al, K and Mn, some of which may form already at 500 K. The phyllosilicates have a notable influence on the dust to gas ratio and will increase the C/O ratio in the gas phase to about 0.83 because of the additional intake of oxygen.

Whether or not these phyllosilicates (and all other condensates) actually form in space is not discussed in this paper, as this is a paper just about chemical and phase equilibrium. Thi et al. (2017) have proposed a kinetic approach for the formation of phyllosilicates in protoplanetary discs, based on a warm surface chemistry model. Thi et al. have found that despite the high energy barriers involved in chemisorption and diffusion, water and hydroxyl can efficiently migrate into the bulk lattice within typical lifetimes of protoplanetary discs if the gas and dust temperatures are of order 80–700 K. The maximum intake of water and hydroxyl is a parameter in these models.

As a side result, we discussed whether tungsten could be the first condensate in space and could provide the seed particles for astrophysical dust formation. We argue that solid tungsten becomes stable already at very high temperatures of order 2000 K, that it could provide just about the right numbers of seed particle per H-atom to form micron-sized silicate grains, but that a strong super-cooling of order 300–400 K is required and that high pressures of the order of milli-bars are necessary to make tungsten nucleation an efficient process.

Acknowledgements

We thank Dr. Wing-Fai Thi for pointing us towards the SUPCRTBL database and for his comments on the manuscript. The research leading to these results has received funding from the European Union Seventh Framework Programme FP7-2011 under grant agreement no 284405. Jasmina Blecic is supported by NASA trough the NASA ROSES-2016/Exoplanets Research Program, grant NNX17AC03G. Christiane Helling highlights the hospitality of the Rijksuniversiteit Groningen and the Universiteit van Amsterdam with travel support from NWO and LKBF.

Appendix A Solving gas phase chemical equilibrium down to 100 K

Solving the system of non-linear equations (Eq. (5)), that is, finding the root vector of atomic partial pressures which satisfies the element conservation equations after elimination of all molecular partial pressures in chemical equilibrium according to Eq. (4), seems to be a simple numerical standard problem. Indeed, a simple Newton–Raphson iteration may lead to quick success at high temperatures. However, at low temperatures, the solution vector becomes more and more degenerate, with some atomic partial pressures approaching values < 10^{−several hundred} dyn cm⁻², for example carbon, whereas others like hydrogen stay of order 1. This leads to a number of numerical problems, in particular:

• (i)

  the conditional number of the Jacobi matrix becomes large, i.e. the resolution of the linearised equation system is not possible without substantial losses of precision;

• (ii)

  the Newton–Raphson method only works satisfactory when a good initial estimate of the solution vector is provided, which becomes increasingly difficult at low temperatures.

We have developed a stable algorithm that solves these problems under all tested circumstances, including unusual element abundances and pressures, down to 100 K. The first necessary step was to switch to quadruple precision arithmetics for temperatures T <1000 K, which solvesproblem (i).

Figure A.1 sketches how we solve problem (ii). The algorithm starts by sorting the elements according totheir abundances. This hierarchy can change substantially by condensation as certain elements will be removed almost entirely from the gas phase, for example Ca and Al, whereas other elements stay in the gas phase for longer, for example S and Cl. In the example shown in Fig. A.1, we consider K = 9 elements. There are nine pre-iteration steps performed in this case, before a final Newton–Raphson iteration is carried out to solve the element conservation equations (Eq. (5)) for all elements. The task of the pre-iterations is to provide good initial estimates of the atomic partial pressures for this final iteration, they are the key to obtain code stability.

In each pre-iteration, we only consider a subset of elements and molecules which are composed of these elements. Each pre-iteration starts with an estimate of the atomic partial pressure of a new element (short red bar) using the values of all other atomic partial pressures computed so far (long grey bar). The pre-iteration then continues by improving the atomic pressures of the last N elements (long red bars, N = 4 in the figure) by taking into account the mutual feedbacks, however without changing the partial pressures of the other, even more abundant elements on the left (short grey bars).

For example, pre-iteration 1 (“pre-1”) only considers hydrogen with species H and H₂, from which a first guess of p_H is obtained. Next comes oxygen with additional species O, OH, O₂ HO₂, H₂ O and O₃. During the first stage of pre-2, a first guess of p_O is obtained by solving the element conservation equation for O at given p_H. During the second stage of pre-2, the feedback between O on H is taken into account by solving the two coupled conservation equations for H and O by a Newton–Raphson method to improve p_H, p_O. The next element is carbon with additional species C, CH, CO, C₂, HCO, CH₂, H₂ CO, CH₃, CH₄, CO₂, C₂ H, C₂ H₂, C₂ H₄, C₂ H₄O, C₂ O, C₃, C₃ O₂, C₄, C₅ from which p_C is first estimated and then determined, including the mutual feedback on H and O, by solving three coupled non-linear algebraic equations, and so on. The algorithm continues this way until Mg, where we stop refining p_H in stage 2 of pre-5. We are using the hierarchical order of the elements here, that is, whatever the outcome of p_Mg will be, it can only have a minor feedback on p_H since ϵ_Mg ≪ ϵ_H. When sulphur enters the pre-iteration phase at pre-8, the atomic partial pressures of O, C and N are already frozen. The algorithm continues this way, before the full system of element conservation equations (Eq. (5)) is solved with the Newton–Raphson method, using the atomic pressures from the pre-iteration phase as initial estimates.

Fig. A.1 Algorithm used in GGCHEM to iteratively solve the chemical equilibrium problem for total number of elements K = 9, and number of elements taken into account during the pre-iterations N = 4, see text.

These pre-iterations in fact work so well, that the final full Newton-Raphson iteration only makes minor corrections < 5% in all our benchmark tests, at most, which produces no visible changes in Fig. 3, for example, only another boost in performance. Figure A.1 sketches the algorithm for N = 4, however we achieved best overall performance with N = 5. Further improvements of code stability and performance are achieved by dividing the resulting atomic partial pressures by their initial guesses after each Newton–Raphson iteration, then store these ratios into the memory, and use them as correction factors next time to obtain improved initial guesses. This simple idea was the key to achieve code convergence down to 100 K.

At T = 96 K, however, the equilibrium constant k_p of the molecule hits 10⁺⁴⁹³², the maximum number that can be represented in quadruple precision, and the code crashes. It might be technically possible to go deeper in temperature by eliminating a few big molecules, however these are just the molecules which dominate at those temperatures, see Fig. 3.

Appendix B Solving phase equilibrium

We formulate the element conservation equations in the presence of condensed species in the following way (B.1)

where are the total element abundances (or “initial” element abundances before condensation), ϵ_k are the element abundances in the gas phase and is the concentration of condensed species j per hydrogen nucleus. N is the current number of present condensates, and s_j,k are the stoichiometric factors of elements k in condensate j. Equation (B.1) is a generalised variant of Eq. (5). We note that all equations for the gas phase in Sect. 2.1remain valid when using ϵ_k.

An essential realisation is that Eq. (B.1) should never be used in computer codes to “subtract the dust”, as this would unavoidably lead to a substantial loss of numerical precision at low temperatures. Some ϵ_k are reduced by several hundreds of orders of magnitude at low temperatures, which will soon exceed the numerical resolution required to perform this subtraction precisely, even when using quadruple precision arithmetics.

We have developed an iterative algorithm which avoids this problem by applying corrections δϵ_k to ϵ_k to eventually solve (B.2)

at given temperature T and density n_⟨H⟩. Provided that we have selected the present condensates correctly (all non-selected condensates must have S < 1), this equation simplifies to (B.3)

The dependencies on ϵ_k are marked here because all gas phase equilibrium results only depend on (n_⟨H⟩, T, ϵ_k) (where n_{⟨H ⟩} and T are parameters to the problem), and so do the supersaturation ratios S_j. We choose N_ind element abundances as independent variables of the problem (B.4)

whereas the dependent variables are (B.5)

A_j,i is the conversion matrix and D the number of dependent elements with indices k(j). One iteration step is completed by

The concentrations of the condensed species c_j are simply abyproduct of these computations, but they are not used during the iterations. The element conservation is implemented in the conversion matrix A_j,i as explainedbelow, and we carefully check that Eq. (B.1) remains valid after each iteration step.

To explain these definitions, let us consider the condensation of Al₂ O₃, MgSiO₃ and Mg₂ SiO₄ as an example. There are three condensates is this case (N = 3) and four affected elements (N_ind + D = 4). In order to proceed, we must select three of them (N_ind = 3, D = 1) to obtain three equations (Eq. (B.3)) for three unknowns. As default rule, we pick the least abundant elements (here Al, Si and Mg) as independentvariables, i = 1 is Al, i = 2 is Si and i = 3 is Mg. The other elements (here O, D = 1) become dependent elements, k(4) is O. To find the conversion matrix A_j,i we write down the conservation equations for all elements affected by condensation as

By a stepwise elimination scheme, we find the conversion matrix A in this case to be

The conversion matrix A is a N × (N + D) matrix which describes how the dependent variables change when the independent element abundances are changed. This transformation allows us to solve Eq. (B.3) as system of N = 3 coupled non-linear algebraic equations for the N = 3 unknowns (δϵ_Al, δϵ_Si, δϵ_Mg), which succeeds easily by a Newton–Raphson iteration with several internal calls of the equilibrium chemistry to compute the Jacobi matrix at every iteration step. We have to make sure, however, to limit the Newton–Raphson step δx to prevent negative c_j or negative ϵ_k. This method is found to generally converge very quickly at all temperatures, within less than about 10 iteration steps, depending on the quality of the initial guesses of ϵ_k (here ϵ_Al, ϵ_Si, ϵ_Mg).

To solve the problem of bad initial guesses for ϵ_k, we store all successfully computed results (ϵ_k, c_j) as function of (n_⟨H⟩, T) into a database and take the initial guesses from the closest (n_⟨H⟩, T) database point next time. This way, the program automatically becomes more stable and faster over time. But initially, we must run GGCHEM with small temperature steps from warm to cold.

Some practical problems may still arise, because of the somewhat unclear selection of condensed species and independent elements. The selection of condensed species is guided by monitoring the supersaturation ratios of all condensates during the iterations. If a new condensate becomes supersaturated, it will be added to the set of selected condensates, or will replace one of them. We have implemented a number and hand-crafted criteria here. For example, one condensate cannot be present simultaneously as liquid and solid. Application of the algorithm explained above would fail already in the first step, as there is only one affected element but two condensates – only one can be present. Another example is CaTiO₃, MgTiO₃ and CaMgSi₂O₆, where the latter is a linear combination of the former two, in which case the elimination scheme to obtain the conversion matrix fails – only two of them can be present under any circumstances. This leads us to some quite interesting general insights into the nature of phase equilibrium:

• (a)

  N = N_ind must hold, that is, the number of simultaneously present condensates must equal the number of independent elements. N_ind is limited by the number of elements affected by condensation K (K = N_ind + D with D ≥ 0). For any given set of condensates, this means that we can simply count the number of elements from which they are made, K, and the number of present condensates in phase equilibrium N must not exceed that number: N ≤ K.

• (b)

  Stoichiometric linear combinations of condensates must not be present simultaneously in phase equilibrium. This applies in particular to the solid/liquid phase transitions.

• (c)

  The number of dependent elements D is given by the number of elements which cannot be completely converted into condensates, due to element conservation constraints, at temperature T. Towards low T, N increasesand D decreases monotonically. Eventually, D becomes very small, only counting elements like O, N, C and H which are affected by condensation but need very low temperatures to completely condense⁷. All other elements eventually vanish almost entirely from the gas phase at low temperatures.

For example, N increases by one around 600 K, as soon as FeS[s] becomes stable. Since , the stoichiometry of that condensate makes it possible to completely remove sulphur from the gas phase. Without such new opportunities to lock away elements, some sudden transformations might occur between several condensates (type-2 phase transitions), but the total number of them will remain constant. Practically, we observe that condition (b) is also valid for linear combinations of condensates which conserve all elements but oxygen, the abundance of which is only marginally affected anyway by condensation.

Appendix C The molecular k_p data

During the preparation phase for this paper, we have collected molecular k_p data from Tsuji (1973), Sharp & Huebner (1990), Lüttke (2002), Stock (2008), and Barklem & Collet (2016) and compared these to our own previous data collection going back to Gail & Sedlmayr (1986). The names of 1155 molecules have been homogenised and the different functional fits of k_p(T) have been overplotted, similar to Fig. 1, and characterised to assess the level of agreement.

Figure C.1 shows the results of the complete comparison study (Worters et al. 2017). Most of the molecules fall into the category “data agrees” (30%) with deviations better than 0.4 dex at 200 K and better than 0.1 dex at 3000 K, disregarding the Tsuji data. However, for the majority of the molecules, the disagreement is worse, with some cases of quite obvious disagreements.

Similar conclusions have recently been published by Barklem & Collet (2016), who have reviewed the thermo-chemical data of 291 diatomic and charged diatomic molecules. These deviations are mostly due to uncertainties in the dissociation energies. Barklem & Collet (2016) published their data in terms of partition functions and dissociation energies (C.1)

where A and B are atoms, AB is a diatomic molecule, m_red = m_Am_B∕(m_A + m_B) is the reduced mass, m_el the electron mass, Q_A, Q_B and Q_AB are the partition functions, and is the dissociation energy. To compute the equilibrium constants k_p(T) according to our definition from this data, also for charged molecules, we need to combine Eq. (C.1) with Saha functions

to obtain

χ_A is the ionisation potential of the neutral atom and χ_A⁻ is the electron affinity of the negative ion. Tables for , χ_A, Q_A (T), Q_A⁺(T), Q_A⁻(T), Q_AB(T), Q_AB⁺(T) and Q_AB⁻(T) are given in the electronic appendices of Barklem & Collet (2016), on 42 temperatures points between 10⁻⁵ and 10⁴ K. The electron affinities are not included, however, so we adopt the following values⁸ (C.7)

Fig. C.1 Statistics of molecular kp(T) agreement between the functional fits provided by Sharp & Huebner (1990); Lüttke (2002); Stock (2008), and Barklem & Collet (2016).

We have fitted these data between 50 and 6000 K with a Stock-function (Eq. (17)).

Our final choice of k_p(T) data for GGCHEM is shown in Table D.1. We have given preference to the NIST-JANAF data fitted by (Stock 2008). For molecules available in (Barklem & Collet 2016), but not available in NIST-JANAF, we used the Barklem & Collet data. The TEA benchmark test revealed a few mismatches between the (Stock 2008) data and the TEA data. For each of those molecules, we have re-fitted the NIST-JANAF data ourselves. We do not use any of the (Tsuji 1973) or (Sharp & Huebner 1990) data directly, although some of the (Stock 2008) data we are using are based on the Tsuji-data.

Appendix D The condensed phase data

We have extracted condensed phase Gibbs free energy data from two sources. First, the NIST-JANAF database (Chase et al. 1982; Chase 1986)⁹ provides -data [kJ mol⁻¹] at p^⊖ = 1 bar in the 7th column of their data-files. In NIST-JANAF, these are the Gibbs free energies of formation from the respective NIST reference states of theelements, which are temperature-dependent (D.1)

Fig. D.1 Fitting the vapour pressures of solid and liquid iron from the NIST-JANAF database. The two fitted pvap (T) curves intersect at 1824 K which is close to the known melting point of 1809 K according to NIST-JANAF. The modelwill automatically give preference to the phase with the lower vapour pressure. Solid data exist between 100 and 2200 K, and liquid data between 298 and 4000 K. The solid curve shoots upward for T > 2500 K, but this does not matter as iron condenses as a liquid at those temperatures. When designing our solid Fe vapour pressure pvap (T)-fit, we just haveto make sure that is does not produce a second, spurious intersection point with the liquid curve.

To obtain the Gibbs free energy of formation of a condensate from free gas atoms, our as input for Eq. (12), we need to subtract as (D.2)

The required input for the vapour pressure (Eq. (9)) is obtained in a similar way by subtracting instead.

Second, the geophysical SUPCRTBL database (Zimmer et al. 2016; Johnson et al. 1992)¹⁰ provides condensed phase data. In this database, the Gibbs free energies of formation are defined in a different way, namely with respect to the reference state of the elements at reference temperature T_ref = 298.15 K (D.3)

The database actually provides according to a term in Eq. (2) of Zimmer et al. (2016). However, V ^⊖ is of order ofa few J bar⁻¹, and can safely be neglected at p  < a few bar as assumed in this paper. In order to arrive at the input for our Eq. (12), we need to subtract corrections as (D.4)

namely the Gibbs free energy differences between the free atoms at temperature T and the respective elements in their standard state at T_ref. We could not figure out a way how to consistently perform these conversions from the SUPCRTBL data, as most of the atoms are lacking. However, assuming that the standard states of the elements at T_ref are identical in both NIST and SUPCRTBL, we can utilise NIST-JANAF as

where , [] and are given in the 7th, 5th and 3rd columns of the NIST-JANAF data-files. We have computed high-precision functional fits of these conversions and list the fit coefficients and functional form in Table D.5.

By applying these atomic corrections to the SUPCRTBL data according to Eq. (D.4), we have converted that data to our reference states of neutral atoms needed in Eq. (12). We have extracted data this way for a large number of minerals known on Earth listed in Table D.4. The table lists functional fits for for 121 solid species, interestingly including some phyllosilicates which host OH or H₂ O in their lattice structure. We have ordered the minerals in Table D.4, somewhat arbitrarily, by , where N is the sum of stoichiometric factors, which gives a first expression of what could be the most thermodynamically stable solid materials (note that SUPCRTBL does not have tungsten). We did not include “aqueous” species (i.e. species solved in liquid water), and applied the following additional selection criteria:

• in case of multiple minerals with the same stoichiometric factors we only included the most stable compound at standard temperature,

• minerals with broken stoichiometric factors, or with any stoichiometric factor >16 are ignored,

• we did not include arsenic (As) or gallium (Ga) compounds.

References

All Tables

All Figures

Fig. 1 Comparison of kp(T) data for molecules HS (mercapto) and CN (cyanogen) used by several groups. The second column of smaller plots shows the deviations of log10kp from the mean values (computed without the Tsuji data). On the right side, the data has been converted to [kJ mol−1] (see Eq. (13)), before computing the deviations from the mean values. The large deviations at small temperatures for the Tsuji (1973) and Sharp & Huebner (1990) data are due to extrapolation errors beyond the valid fit-range. HS is classified as belonging to the group “data disagrees”, and CN to “data disagrees at low temperatures” in our comparison catalogue (Worters et al. 2017) since even without the Tsuji data, there are relatively large differences to the new (Barklem & Collet 2016) data.

In the text

Fig. 2 Comparison of condensed phase Gibbs free energy data from Sharp & Huebner (1990, orange), fits to the NIST-JANAF database done by the authors of this paper for GGCHEM (green), and fits to the SUPCRTBL data (Zimmer et al. 2016, blue). The SUPCRTBL data points have been generated for temperatures 100 to 2500 K according to the equations and parameters given in Zimmer et al. (black dashed curve). Lower plot: deviations [kJ mol−1] with respectto the mean of all three data sources.

In the text

Fig. 3 Benchmark test against TEA (Blecic et al. 2016) at constant p = 1 bar and T = 100–6000 K, showing mixing ratios ni∕n. The coloured full and dashed curves are the GGCHEM results. The TEA results are overplotted with narrow grey lines. Only the most important molecules are shown, that is, those which reach a certain threshold concentration, depending on element. The TEA results stop at 400 K as TEA has difficulties converging at lower temperatures, even using 2000 iterations. The following molecules with notable concentrations are missing in TEA: SiH2, SiH3; CaH; TiS, TiH, TiN, TiC; CrH, CrS; MnH, MnS, MnCl, MnF, MnO; FeH; NiH, NiF, NiO, and the following molecules have been omitted by us in the TEA model to improve its convergence: Si(CH3)4, SiCH3 Cl3 and Ni(CO)4. Differences in the selection and availability of molecular species explain all visible deviations between GGCHEM and TEA.

In the text

	Fig. 3 continued.
In the text

	Fig. 4 Concentrations of major molecules at 1 bar according to a simplified nine-molecule model in chemical equilibrium (full lines), in comparison to a model with 24 elements and 445 molecules (dashed) and a full model with equilibrium condensation (dotted).
In the text

	Fig. 5 Upper part: dust to gas mass ratio in phase equilibrium as function of temperature at constant p = 1 bar for solar element abundances. The arrows indicate where the elements start forming condensed phases in significant amounts. The elements marked in red have a significant influence on the shape of the dust/gas(T)-curve. Lower part: remaining element abundances in the gas phase.
In the text

	Fig. 6 Onset of condensation at p = 1 bar for solar abundances in phase equilibrium. The graph show the concentration of the various condensed species with respect to hydrogen nuclei.
In the text

	Fig. 7 Concentration of condensed species at medium temperatures 460–1000 K at p = 1 bar for solar abundances. Please note the changed scaling of the y-axis in comparison to Fig. 6.
In the text

	Fig. 8 Concentration of condensed species at temperatures 100–460 K at p = 1 bar for solar abundances.
In the text

	Fig. 9 Impact of CO-blocking on the composition of the gas at constant pressure and temperature. In the upper plot, the full lines represent the results from a pure gas phase model without condensation (24 elements and 445 molecules, see Sect. 5.1), and the dashed lines show the results from the model including equilibrium condensation (Sect. 5.3). ϵO is a fixed value here, whereas ϵC varies.
In the text

	Fig. 10 Gas phase/condensed species which contain most of the elements in the (p, T)-plane. Solid species are marked by “[s]” and liquid species by “[l]”. At low temperatures, the elements Mg, Si, Na, Al and K are mostly present in form of phyllosilicates. The only further element found to form phyllosilicates in phase equilibrium is manganese.
In the text

	Fig. 10 continued.
In the text

Fig. 11 C/O ratio in the gas phase affected by the formation of metal oxide and hydroxide molecules (left) and condensation of silicates, phyllosilicates, and water ice (right). On the left side, we plot the “effective” C/O that is contained in the gas in form of molecules made of H–C–N–O only, but do not count metal oxide andhydroxide molecules. The contour lines shown are 0.58 (after SiO formation) and 0.8 (after Mg(OH)2 and Fe(OH)2 formation). On the right side, we plot the resulting gas phase abundances after the condensation of solids and liquids. The increase of C∕O > 0.7 is due to the consumption of oxygen by the major silicates like MgSiO3[s] and Mg2SiO4[s], and the further increase to >0.8 is due to phyllosilicates. The sharp increase at low temperatures (red area) is due to water condensation.

In the text

	Fig. 12 Concentration (xj = nj∕n⟨H⟩) of the three major magnesium-silicates in the (p, T)-plane in the equilibrium condensation model.
In the text

	Fig. 13 Left: nucleation rate of tungsten per hydrogen nucleus J⋆∕n⟨H⟩ [s−1] as function of gas pressure and temperature. Centre: size of the critical cluster N⋆. Right: nucleation timescale τnuc [yr] according to Eq. (25). The white dashed contour line marks saturation S = 1, the gas is supersaturated to the left of this line.
In the text

	Fig. A.1 Algorithm used in GGCHEM to iteratively solve the chemical equilibrium problem for total number of elements K = 9, and number of elements taken into account during the pre-iterations N = 4, see text.
In the text

	Fig. C.1 Statistics of molecular kp(T) agreement between the functional fits provided by Sharp & Huebner (1990); Lüttke (2002); Stock (2008), and Barklem & Collet (2016).
In the text

Fig. D.1 Fitting the vapour pressures of solid and liquid iron from the NIST-JANAF database. The two fitted pvap (T) curves intersect at 1824 K which is close to the known melting point of 1809 K according to NIST-JANAF. The modelwill automatically give preference to the phase with the lower vapour pressure. Solid data exist between 100 and 2200 K, and liquid data between 298 and 4000 K. The solid curve shoots upward for T > 2500 K, but this does not matter as iron condenses as a liquid at those temperatures. When designing our solid Fe vapour pressure pvap (T)-fit, we just haveto make sure that is does not produce a second, spurious intersection point with the liquid curve.

In the text