FastChem Cond: Equilibrium chemistry with condensation and rainout for cool planetary and stellar environments

Cool astrophysical objects, such as (exo)planets, brown dwarfs, or asymptotic giant branch stars, can be strongly affected by condensation. Condensation does not only directly affect the chemical composition of the gas phase by removing elements but the condensed material also influences other chemical and physical processes in these objects. This includes, for example, the formation of clouds in planetary atmospheres and brown dwarfs or the dust-driven winds of evolved stars. In this study we introduce FastChem Cond, a new version of the FastChem equilibrium chemistry code that adds a treatment of equilibrium condensation. Determining the equilibrium composition under the impact of condensation is complicated by the fact that the number of condensates that can exist in equilibrium with the gas phase is limited by a phase rule. However, this phase rule does not directly provide information on which condensates are stable. As a major advantage of FastChem Cond is able to automatically select the set stable condensates satisfying the phase rule. Besides the normal equilibrium condensation, FastChem Cond can also be used with the rainout approximation that is commonly employed in atmospheres of brown dwarfs or (exo)planets. FastChem Cond is available as open-source code, released under the GPLv3 licence. In addition to the C++ code, FastChem Cond also offers a Python interface. Together with the code update we also add about 290 liquid and solid condensate species to FastChem.


INTRODUCTION
Condensation is an important process in astrophysical objects that strongly affect the chemical composition of their gas phases.Typically, it can become very important at temperatures lower than roughly 2000 K, when the first, highly refractory species usually start to condensate out of the gas phase.Not only are the abundances of molecules directly affected by this process (e.g.Visscher & Fegley 2005), but the presence of condensing species also play an important role in other physical and chemical processes in such objects.
For example, condensation in atmospheres of (exo)planets and brown dwarfs can lead to the formation of cloud particles that strongly impact the spectral appearance and temperature structures of these objects (e.g.Burrows & Sharp 1999;Ackerman & Marley 2001;Helling & Woitke 2006;Kitzmann et al. 2011;Morley et al. 2012;Lee et al. 2016;Gao et al. 2021).Moreover, the formation of dust grains in the atmospheres of asymptotic giant branch (AGB) stars, on the other hand, results in a massive stellar wind that finally leads to the loss of almost the entire stellar envelope to the interstellar medium over time (e.g.Hoefner et al. 1998;Winters et al. 2000;Gail & Sedlmayr 2013).Other environments where condensation is an important process are protostellar nebulae.For example, condensation sequences in the solar nebula are studied to determine the composition of planet building blocks and meteorites (e.g.Fegley ★ E-mail: daniel.kitzmann@unibe.ch& Palme 1985;Palme & Fegley 1990;Lodders & Fegley 1993).Thus, the treatment of condensation and its impact on the gas phase composition in an equilibrium chemistry code is a fundamental necessity to achieve a deeper understanding of cooler astrophysical environments.An extensive overview of condensation in astrophysical environments and a summary of modelling efforts is discussed by, for example, Ebel (2006).More detailed discussion and a historical overview of the different approaches to the condensation problem and corresponding chemistry codes is provided by Lodders & Fegley (1997) or Ebel (2021).
The chemical equilibrium approximation simplifies the problem of describing condensate formation significantly, albeit at the expense of assuming local thermodynamic equilibrium.Although this assumption is not flawless, it remains reasonable for various scenarios.Nevertheless, this equilibrium is not attained in regions with very low pressures and temperatures, such as the outer circumstellar shells of AGB stars.In such an environment, for example, the timescale of the hydrodynamic gas flow is much shorter than that of chemical reactions, resulting in a state where the chemical composition essentially becomes frozen (Gail & Sedlmayr 2013).
Yet, despite the inadequacy of chemical equilibrium condensation in certain environments, its computation remains a crucial necessity.Chemical equilibrium calculations establish a fundamental benchmark, serving as a foundation to extend computations to encompass kinetic effects on chemistry.Without knowledge of equilibrium conditions for reference, it becomes unfeasible to ascertain whether a system is in or out of equilibrium, or to determine the degree of deviation from equilibrium (Ebel 2021).
Equilibrium chemistry calculations are typically done either via the law of mass action in combination with the element conservation equations (Brinkley 1946(Brinkley , 1947) ) or by using a Gibbs free energy minimisation approach (White et al. 1957(White et al. , 1958)).Both formalisms yield the same solution but are numerically different (Zeleznik & Gordon 1960, 1968;Smith & Missen 1982;van Zeggeren & Storey 1970).The first forms a system of non-linear equations while the latter is a global optimisation problem.On the one hand, examples of computer codes that use the law-of-mass-action approach are CON-DOR (Lodders & Fegley 1993), GGChem (Woitke et al. 2018), or FastChem (Stock et al. 2018(Stock et al. , 2022)).The Gibbs free energy minimisation method, on the other hand, is implemented for example in the RAND method (White et al. 1958), the NASA-CEA code (Gordon & McBride 1994;McBride & Gordon 1996), SOLGAS & SOL-GASMIX (Eriksson 1971(Eriksson , 1975)), or STANJAN (Reynolds 1986).A more detailed overview and descriptions of the different numerical methods used to compute the chemical equilibrium composition can be found in, for example, Smith & Missen (1982) or van Zeggeren & Storey (1970).
Determining the equilibrium composition in a case where condensation occurs is quite more elaborate than a pure gas-phase calculation since only a limited number of condensates can co-exist in equilibrium with the gas phase at the same time according to an extension of Gibbs' famous phase rule (Gibbs 1876;Denbigh 1955).In particular, it is in general not a priori obvious which species are stable.The phase rule only limits the number of stable condensates but doesn't restrict the potential combinations of condensates that can be present.This is especially problematic at low temperatures when many different species could in principle condense and, thus, the number of possible combinations could be very large.
Equilibrium condensation is therefore tpyically treated as a twostep process that includes solving for the combined gas-condensed phase chemical composition as well as finding the set of stable condensates.The latter step is often done in an iterative way by adding and/or removing potential condensates from a list of stable condensate candidates until the phase rule is satisfied and element conservation by the gas and condensed phase species is achieved (e.g. Gordon & McBride 1994;Woitke et al. 2018).Such iterative methods, however, become very challenging once the number of potential condensates is large.In such a case, additional a priori information on condensation sequences might be required to limit the number of combinations.While an iterative way is an established approach to solve this problem in astrophysical contexts, alternative methods do also exist (e.g.Leal et al. 2016).
The first two version of our equilibrium chemistry code FastChem focused solely on the gas phase chemistry.In the original publication (Stock et al. 2018) the basic semi-analytical approach that our chemistry code uses was introduced but was strictly limited to chemical systems including hydrogen.In FastChem version 2 (Stock et al. 2022) we then, among other things, adapted the formalism to treat arbitrary element compositions.
In this study we adapt FastChem to include equilibrium condensation.With an internal version number of 3.0, we refer to this updated version as FastChem Cond.In addition to equilibrium condensation, FastChem Cond also employs the rainout approximation as an alternative condensation description that is commonly used in (exo)planetary science or for atmospheres of brown dwarfs (Marley et al. 2013).In Section 2 we first briefly summarise the theoretical concepts of equilibrium condensation and the related computational challenges.Section 3 describes the numerical framework implemented in FastChem Cond, while in Section 4 we test our new equilibrium chemistry code on a set of different scenarios from the literature that involve condensation.

Gas phase treatment
Below we very briefly summarise the equations describing chemical equilibrium in the gas phase using the law of mass action.For a more detailed description of the notation and formalisms used in the following we refer to Stock et al. (2018) and Stock et al. (2022).
The number density of gas phase species   is given by the law of mass action: where   is the number density of elements  ∈ E in their atomic form,   the temperature-dependent equilibrium constant of species , and    their stoichiometric coefficients.The set E contains all elements, while the set S includes all gas phase species.
With the relative, normalised element abundances ε  , the total number density of atomic nuclei of a given element  ∈ E is given by where the total number density of all atomic nuclei is defined by We can then write the corresponding element conservation equation for element  ∈ E: Finally, assuming an ideal gas, the gas pressure  g is given by where  B denotes Boltzmann's constant and  the gas temperature.
Equations ( 1), (4), and (5) form a system of linear and non-linear algebraic equations for the number densities of all gas phase species.The semi-analytical computational procedure used to solve this system is described in detail by Stock et al. (2022) (FastChem 2).

Stoichiometric condensates
In the following we only consider neutral, stoichiometric condensates1 .Other cases, such as solid solutions for example, are not considered here.
Let C = { 1 , . . .,  | C | } be the set of all considered condensates in the chemistry model, where |C| denotes the total number of condensates.
Unlike gas phase species described in the previous subsection, not all condensates can stably exist in equilibrium with the gas phase, but only a subset C  ⊂ C. Depending on the chosen values for the temperature , the gas pressure  g , and the element abundances ε  , C  can also be an empty set, i.e. no condensate exists.This is the case for sufficiently high temperatures, for example.Furthermore, we introduce the subset of elements E  ⊂ E that are contained in the stable condensates C  .
The stability criterion of a condensate  can be expressed by using its activity   , given by the corresponding law of mass action (Gail & Sedlmayr 2013): where   is the equilibrium constant for the condensate.
For numerical computations, Eq. ( 6) can also be more conveniently expressed in its logarithmic form: In the following, we will use this logarithmic form of the activity equation unless stated otherwise.

Coupling of gas phase and condensates
When condensates are taken into account, the element conservation equation is extended by a fictitious number density of condensed molecules   : where in analogy to Eq. ( 2)   is defined as with This fictitious number density   should not be confused with an actual dust particle number density.On the contrary, the   are merely used in Eq. ( 8) to keep track of the elements bound in the condensed material.It is also important to note that the fictitious number densities of the condensates do not contribute to the pressure of the gas phase in Eq. ( 5).
An alternative way to describe the depletion of elements from the gas phase by condensates is the degree of condensation (see Gail & Sedlmayr 2013, for example), given by: Using this degree of condensation, the effective element abundance   in the gas phase can be expressed via: where instead of   ,   is then used to calculate the normalised element abundances ε  .
With the element abundances without considering condensation, ε  (   = 0), and the total number density of atomic nuclei  ⟨g⟩ , we can also define the largest possible fictitious number density  ,max of a condensate: assuming that the least common element in a condensate is completely condensed.

Stability criterion for condensates
The activity equation ( 7) provides a stability criterion for condensates.For those that exist in equilibrium with the gas phase, the activity has to be   = 1 and, thus, All other (unstable) condensates have an activity of   < 1, yielding: In chemical equilibrium, therefore, only the following two different situations are permitted for all condensates:

The phase rule
Equation ( 14) can be seen as constraints for the densities   of the elements E  that are contained in stable condensates.Since the gas phase contains |E | unique elements, there is a maximum number of |E | linearly-independent equations ( 14).This a priori limits the number of stable condensates to |C s |≤|E |, which is usually a much smaller number than the total number of condensates |C| considered in a chemical system.
As described in Sect.2.1, the   are the result of a system of nonlinear equations involving the law of mass action for the gas phase species, the element conservation equations, as well as the ideal gas law.In a case that includes stable condensates, however, some   are also fixed by the solution of the Eq. ( 14).Consequently, they are therefore no longer free variables for the gas phase itself.
Thus, if all elements in the system are contained in stable condensates at the same time, their densities   in the gas phase are all determined by Eqs. ( 14).The densities of all molecules are then directly given by their corresponding law of mass action with the fixed   following Eq.( 1).This, however, creates the problem that the resulting   and   likely do not yield the correct pressure of the gas phase in Eq. ( 5) anymore.
Therefore, at least one incondensable element in the gas phase is needed, whose   is determined by the correct value of  g .Consequently, the total number of stable condensates is effectively limited to This is an extension of the famous, so-called Gibbs' phase rule (Gibbs 1876(Gibbs , 1878;;Denbigh 1955).While the phase rule limits the number of stable condensates, it unfortunately does not directly provide information on the condensates contained in the set C s , allowing for a potentially large number of possible combinations.

Solution strategies
One fundamental difference between the law of mass action for molecules (Eq.( 1)) in comparison to its formulation for condensates in Eq. ( 6) is that the former directly yields the number density of a molecule if the   are known, while the latter provides only a stability criterion.The number densities   appear only in the element conservation equation ( 8) and the degree of condensation (11).Thus, there is no direct way of obtaining   .
This non-linear problem can be solved in a variety of different ways, two of which are briefly summarised in the following.For a broader overview of other methods, we refer the reader to van Zeggeren & Storey (1970).We note, however, that the two schemes discussed next both assume that the set of stable condensates C  is already known when the non-linear system is solved.

Iterative solution with Newton's method
One possible solution for the coupled problem is to use an iteration based on Newton's method.This, for example, is done in the GGChem code by Woitke et al. (2018).Their model uses the effective element abundances ε  and the fictitious condensate number densities   as free variables.
Using Eq. ( 14) and the corresponding conservation equations for a considered set of stable condensates, GGChem assembles a linear system of equations to solve for the changes in the number densities   for given changes in the ε  .This is coupled to the gas phase calculations with ε  in a Newton's method scheme until   = 1 is achieved for all considered condensates.In this case, the Jacobian matrix is not given analytically but is evaluated numerically by repeatedly calling the gas-phase-only calculations.
A Newton's method is also used in the so-called NASA-method approach (Huff et al. 1951) 2 .In contrast to GGChem, however, the NASA method solves the gas and the condensed phase in a combined Newton's method rather than separating the compution into two distinct parts.

Direct solution
An alternative way is used in, for example, the Condor code by Lodders & Fegley (1993).Here, Eq. ( 14) for all considered, stable condensates is solved directly for the element densities   .Using these   in the gas phase calculations and comparing the results with those of the condensate-free case allows to obtain the degree of condensation   for all elements and, subsequently, the number densities   .No separate Newton's method is required here since the requirement of   = 1 is automatically satisfied by directly solving Eq. ( 14) and keeping the   fixed in the gas phase calculations.

Deriving the set of stable condensates
The two computational approaches discussed in the previous subsections assume that the set of stable condensates is already known.Usually, however, this is rarely the case and the solution of this problem therefore has to be embedded in some form of iterative scheme.
A typical calculation would, thus, start with determining the gas phase composition in the absence of any condensates.This provides initial guesses for the   .Using these initial solutions, one can test for the potential presence of stable condensates by calculating their activities.If all   are smaller than unity, then no condensate can be present and the calculation is finished.
Else, if some condensates have activities larger than unity, one then has to find a solution in terms of the element number density   and the fictitious density   , such that the activities of all stable condensates are unity and the element conservation is fulfilled.The solution also has to agree with the phase rule and potentially stable condensates considered in the solution have to yield a linearly independent system of equations.
Especially at low temperatures many different condensates can initially have activities   > 1 and, thus, are potential members of the set C  .This problem is commonly solved by constantly adding and removing condensates from C  , trying to solve the non-linear system with the methods described above (e.g. Gordon & McBride 1994;Woitke et al. 2018).This iteration has to continue until all condensates in C  have an activity of   = 1, while all others have   < 1.
During the course of the iteration it might become necessary to remove condensates again from the set C  , if, for example, they yield negative densities   .Considering, however, that many condensates compete for the same elements it is not always immediately obvious, which condensates need to be removed to stabilise the system.
Given the large number of possible combinations of stable condensates at low temperatures, finding the final set C  is, therefore, extremely challenging without any prior knowledge.One possibility is to start at high temperatures and slowly cool down the system, using the results of the higher temperatures as initial guesses for the lower temperature calculations (e.g.Burrows & Sharp 1999).GGChem, for example, additionally creates a lookup table of the effective element abundances as a function of pressure and temperature that allows to obtain the results in later calculations more easily.This lookup table, however, obviously depends directly on the individual element abundances and needs to be created for every set of metallicity values.

CONDENSATE TREATMENT IN FASTCHEM
As illustrated in the previous section, solving the combined condensate-gas phase problem is challenging, especially under low temperature conditions when many different condensates could potentially exist.The solution to this problem is made complicated by the fact that initially the final set of stable condensates C  is not known.The phase rule and the requirement that the equations ( 14) for all considered condensates need to be linearly independent offers some guidance.However, even with this supplementary information, an iterative approach of adding and removing condensates as described above does often not converge properly without further a priori information or good initial values for e.g.Newton's method (Burrows & Sharp 1999).
The necessity of selecting an initial set C  is caused by the fact that Eq. ( 14) is only valid for stable condensates.Consequently, it is not possible to include all potential condensates at once in the system to derive the set C  by solving the combined problem as evidently some of them might require a solution with   < 1 in case they are not stable.
In this new FastChem version, we choose an alternative approach to this problem.This new approach is based on the basic ideas presented by Kulik et al. (2013) and Leal et al. (2016), the latter of which can be seen as a variant of the NASA method (Huff et al. 1951;van Zeggeren & Storey 1970).In the following, we adapt their equations and numerical methods to the general computational framework used in FastChem.

Modified condensate equations
The most important modification of the standard approach is to adapt the activity equation in chemical equilibrium (Eq.( 14)), so that it is valid for all condensate species C instead of only for species in the (initially unknown) set C  .This is achieved by introducing a correction factor   ≥ 0, such that Obviously, for stable condensates ( ∈ C  ), the correction factor is   = 0, while for the unstable ones ( ∈ C \   ), we require   = − ln   in equilibrium.As discussed by van Zeggeren & Storey (1970) and Leal et al. (2016), the new correction coefficients   are directly connected to the Lagrange multipliers of the Gibbs free energy minimisation approach (e.g.White et al. 1957White et al. , 1958)).
Since a new variable   is introduced, an additional equation is required to solve the combined system of equations.This auxiliary equation combines   with the condensate number density   : This is supplemented by the element conservation equation, Eq. ( 8), involving all condensates.
The combined system ensures that the two possible states for equilibrium condensation are a natural outcome of the problem, i.e.
This approach, however, comes at the price of a larger system of equations, i.e. its dimension is in general For a numerical treatment, Eq. ( 17) needs to be adapted slightly.The constraint that either   or   need to be exactly equal to zero is numerically too strong.As suggested by Leal et al. (2016) the equation is, therefore, modified with a constant   , such that The constant   is a very small number, chosen such that it does not impact the outcome of calculation to a large degree.For   we use the expression where  is a small number with a default value of 10 −15 .The quantity   effectively controls the minimum non-zero number density of unstable condensates.The use of   also implies that for stable condensates, ln   is not exactly zero but has a very small value of the order of .
Since both   and   can become very small numbers, we use the logarithmic form of Eq. ( 18): In summary, we have to solve the following new, non-linear system of equations: for the unknown quantities   ,   , and   .Additionally, we require that all of these quantities have values larger than zero.We note that the phase rule is never used explicitly in the construction of the system.Yet, the solution of the system of equations and, thus, the final set of stable condensates will automatically satisfy the phase rule.Furthermore, we also have not to pay any special attention to which condensates can be included to avoid linearlydependent equations.

Newton's method
The combined system Eqs.( 21) -( 23) introduced in the last subsection is solved by using Newton's method with an analytic Jacobian matrix.Instead of solving for the variables   ,   , and   , we replace them by their logarithmic counterparts ln   , ln   , and ln   .This provides a better numerical stability and also guarantees that these quantities cannot become negative.
For a Newton step  we perform first-order Taylor series expansions of the functions    ,    , and    around the initial values , and  ( )  .For    this yields: which can be re-written as with For the activity functions    we perform the Taylor series expansion with respect to ln   and ln   : where ln ).This can also be written as with A similar expansion for the element conservation equations    yields: Here we have used the law of mass action (Eq.( 1)) to analytically calculate the derivatives of the species densities   with respect to those of the elements   and to express the number densities as a function of the initial element densities  ( ) .This can be rearranged to: with . (32) The final system of equations then has the general form where  and 0 are the identity and zero matrix, respectively.The other components of the Jacobian matrix are given by the |E | × |E | matrix the as well as the |C| × |C| diagonal matrix The components of the vector of unknowns are: while those of the right-hand-side vector are given by Eqs. ( 26), ( 29) & (32).The total number of unknowns for this system is 2 |C| + |E |.
Technically speaking, the system of equations is usually smaller than that because elements not forming any potentially stable condensate don't need to be included in the system.With the computed solution for the quantities Δ ln   , Δ ln   , and Δ ln   , the corresponding quantities are then corrected following ln where  is either   ,   , or   .As mentioned earlier, the use of logarithmic quantities prevents these values from becoming negative or zero by construction.

Reduction of the system of equations
The system of equations ( 33) can be become rather large if the number of potential condensates is very high, which typically happens at low temperatures.Its size, however, can be reduced significantly by removing both Δ ln   and Δ ln   from the system of equations.
To remove Δ ln   we first solve Eq. ( 25) for Δ ln   : In Eq. ( 28) we can now replace the quantity ln Δ  , which yields b This reduces the linear system of equations to the form for the |C| + |E | unknowns Δ ln   and Δ ln   .With the Δ ln   from the solution of this system, we can then calculate the corresponding corrections Δ ln   from Eq. ( 40).
The same approach can be used for the unknowns Δ ln   (Leal et al. 2016).Using Eq. ( 41), we can derive an equation for Δ ln   : and use it to eliminate Δ ln   in Eq. ( 31).This yields a system of equations where the Δ ln   are the only unknowns: with the right-hand-side vector components b With the resulting Δ ln   we can calculate the corrections Δ ln   from Eq. ( 44) and subsequently again Δ ln   from Eq. ( 40).Instead of the initial |E | +2|C| dimensions, the system of equations is now reduced to just |E | equations for the elements.Given that usually |E | ≪ |C|, this provides a huge increase in computational performance.
Unfortunately, using this replacement for all condensates can yield a numerically unstable system.The reduced Jacobian matrix in Eq. ( 45) has terms of the form   /  .For stable condensates,   approaches values of close to zero.In fact, with the original equation ( 17) we would effectively perform a division by zero here in case of stable condensates with   ≫ 1.This generates very large, offdiagonal entries in the Jacobian, which can potentially make the matrix essentially numerically singular.
For unstable condensates we have   ≈ 0 and   ≈ − ln   .In such a case, the problematic terms vanish from the Jacobian and the matrix remains well-conditioned.Thus, while stable condensates have to be included in the Jacobian matrix for numerical stability, unstable ones can be safely removed.The system of equations would, therefore, have a dimension of |C  | + |E | once all unstable condensates have been identified and accordingly removed from the Jacobian matrix.

Numerical solution
To allow for more flexibility in the numerical calculations, FastChem Cond contains all three different system of equations introduced so far.Before any of these systems of equations can be solved, the corresponding Jacobian matrix requires some rescaling to ensure numerical stability.The different matrix components of the Jacobian imply that entries can potentially differ by many orders of magnitudes based on the current values of   ,   , or   .This can result in an unstable numerical behaviour.Before solving the system, each row of the matrix is therefore scaled by the inverse of its largest entry.Mathematically, this corresponds to the multiplication of the Jacobian with a non-singular, diagonal scaling matrix (cf.Deuflhard 2011) and usually leads to a well-conditioned Jacobian matrix.
We solve the system of equations directly by using an LU decomposition with partial pivoting from the Eigen library3 .This method provides a fast and stable way to obtain the solution.In principle, we could also exploit the fact that the Jacobian has a very sparse structure, since it contains more zero entries than non-zero ones.Thus, we could use iterative methods that have been especially designed for such systems with a sparse coefficient matrix.However, the overall dimension of our system is normally still quite small and these iterative methods, therefore, tend not to perform better than the direct solution via the LU decomposition.
In some cases, the system of equations can still become illconditioned which makes a numerical solution challenging.Especially the Jacobians for the reduced systems ( 33) and ( 45) are prone to become numerically singular.If, through the course of the Newton iterations, both   and   for certain condensates become very small, the corresponding coupling terms in the submatrices J  C and J   for these condensates are so small that they yield an almost singular matrix.This, of course, creates problems when using the LU decomposition with partial pivoting because this method a priori assumes that the matrix is invertible.Thus, the numerical solution for such an ill-conditioned system can in general not be trusted.
The full system (33) is much less affected by these issues, since there are still coupling terms remaining in the Jacobian that stabilise the numerical solution.However, in rare cases even the full Jacobian matrix can become ill-conditioned.
The Eigen library also offers, for example, an LU decomposition with full pivoting.While this method is slower than the version with partial pivoting, it can detect cases with singular matrices.FastChem Cond is able to use both solvers but by default employs the one with partial pivoting as it provides a much higher computational speed.The user can optionally choose the LU decomposition with full pivoting.
If FastChem Cond uses the method with full pivoting and detects a non-invertible matrix, it will either use a singular value decomposition to solve the system of equations or perturbs the Jacobian to make it invertible (Dennis & Schnabel 1996).However, we note that both methods can lead to solutions that produce so minor corrections, such that the Jacobian matrix will remain non-invertible in the next iteration step.In that case, the iterative method in FastChem Cond will fail to converge.Consequently, in such a situation it is usually advantageous to take a dampened Newton step into any random direction provided by the unstable LU decomposition.Most of the time this produces a well-conditioned matrix again in the next iteration step.
We also note that the full system of equations ( 33) is much less affected by a potentially singular Jacobian matrix.Even in the case of very small values for   and   , there are still additional coupling terms left in the matrix to yield a stable matrix inversion.In fact, using the full system tends to require less iterations to converge than the reduced system of equations.However, the numerical cost per iteration step is much higher, which makes this approach overall often slower in terms of computational time.FastChem Cond will switch to this method if the iteration with the reduced system does not converge.The user can also optionally force FastChem Cond to always use the full Jacobian matrix in its calculations, which might help to overcome convergence issues.

Coupling with the gas phase chemistry
The condensation part of FastChem Cond described above, needs to be iteratively coupled to the gas phase chemistry calculations.At the start of each coupled calculation, FastChem will first determine the gas phase composition, neglecting any condensates.With that solution, the activities of all condensates are obtained.In case no activity is larger than unity, no stable condensates are possible and the calculation is finished.
If FastChem finds some species, with initial   larger than one, it selects these condensates as potential candidates for the condensed phase and proceeds with the solution of the system described in the previous section.In principle, the mathematical description allows to include all condensates available in the chemical model, even those with activities smaller than unity.Since, however, those are not stable by definition, only condensates that can potentially be present are included in the calculation to limit the computational cost of solving the system of equations.
Having determined the set stable condensates and their corresponding densities   , the depletion of the elements from the gas phase is evaluated using the degree of condensation from Eq. ( 11) and the effective element abundances   ,  ∈ E  .In principle, one could now use these   and   values to re-calculate the gas phase for all elements.However, at low temperatures some of the   can become numerically exactly unity, which yields an effective element abundance of   = 0. Consequently, this can create numerical issues in the gas phase calculations for these elements.
Instead, we fix the number densities of the elements in their atomic form,   , in the gas-phase calculations to the solution we obtain from the calculation of the condensed phase.As mentioned earlier, these   are essentially fixed by the solution of the activity equations and, therefore, should actually not be considered anymore as free variables in the gas phase.
An obvious advantage of using the second approach is that the calculation of the gas phase is simplified, since some elements do not need to be calculated at all as their   are already determined.This is especially helpful at lower temperatures when potentially many elements are taking part in condensation.
FastChem Cond will iteratively perform the calculations of the condensed and gas phase chemistry until the coupled system is converged.Normally, the first condensation iteration takes the longest time since this is the step where the selection of the stable and unstable condensates is made.Later calls of the condensation system typically need much fewer Newton steps because the sets of stable and unstable condensates usually don't change significantly.After each calculation of the gas phase FastChem Cond checks for additional, potential condensates that need to be added to the system.The calculation proceeds until all condensates fulfill their stability criterion from Sect.(2.2.2) and all elements are conserved.

Rainout approximation
Besides the standard equilibrium condensation, FastChem Cond also supports calculations using the so-called rainout approximation.Sharp & Huebner (1990).The plots show the partial pressures of gas-phase species as a function of temperature and roughly mirror the results shown by Sharp & Huebner (1990) with some species moved to other panels for better visualisation.We note that that our calculations extend to lower temperatures than those presented by Sharp & Huebner (1990).Following the corresponding figures by Sharp & Huebner (1990), solid lines refer to the chemistry calculations with equilibrium condensation while dashed lines refer to those without.
In the former case, each temperature and pressure point is treated individually.However, as pointed out by, for example, Lodders & Fegley (2002) or Marley et al. (2013) this pure equilibrium condensation is not suitable for calculating the chemical composition of planetary atmospheres or brown dwarfs, because rainout affects the element distribution as function of altitude.
As discussed by Lodders & Fegley (2002), condensates in such high-gravity environments tend form directly from the gas phase (primary condensates) that then settle into cloud layers.In a lowgravity case, such as a protoplanetary disk or a pre-stellar nebula, these primary condensates do not settle into distinct layers but remain dispersed in the gas phase.This allows for the formation of secondary condensates via chemical reactions between the gas phase and (primary) condensates when the temperature decreases further.
To simulate the rainout, FastChem Cond will first start the calculation at the lower boundary (highest pressure) for a given temperatepressure profile  ( ) and work its way upwards towards lower pressures.If condensation is encountered at some pressure , FastChem Cond will solve the coupled condensation and gas phase system as previously described.The converged system then yields the effective element abundances   of the condensing elements left in the gas phase (see Eq. ( 12)).Using these effective abundances, FastChem will then change the actual element abundances   to the computed   at all pressures below the level where condensation occurred.This simulates the rainout of these elements into condensate layers and leads usually to a rapid decrease in the abundances of condensing elements in the upper parts of an atmosphere.
An example to illustrate the differences between equilibrium condensation and rainout is shown in Sect.4.3 for an atmosphere of a brown dwarf.

Code update and availability
The numerical treatment for condensation and rainout described in the previous subsections has been added to the FastChem code.It is available as open source on GitHub (https://github.com/exoclime/FastChem).Like the previous versions of FastChem, the code is released under the GNU General Public License version 3 (gnu.org2007).This new FastChem release has a version number of 3.0 and is referred to as FastChem Cond.We note that for calculations not involving condensation, the results of FastChem Cond are identical to those calculated with the previous version FastChem 2.
As described by Stock et al. (2022), FastChem is written mainly in object-oriented C++ but additionally offers a Python interface (pyFastChem) that allows it to be imported as a normal Python module.Some example Python scripts are included in the repository that showcase the use of FastChem Cond for several different scenarios.In addition to the GitHub repository, pyFastChem is also available as a PyPI package that can easily be installed via pip.
The Python interface of FastChem Cond has been adapted to account for calculations involving condensates.New parameter options have also been included in the interface that allow the user to access special internal parameters.A detailed description of the C++ code and the Python interface is available in the manual, located in the GitHub repository.
Thermochemical data for condensates mostly based on the JANAF tables (Chase 1998) has been added to FastChem Cond.In particular, the current version of FastChem Cond includes about 290 solids and liquids for 27 elements more abundant than germanium for the solar element abundances by Asplund et al. (2009).Details on the condensate data can be found in Appendix A. We aim to add more elements, gas phase species, and condensates in future releases of FastChem.

TEST CALCULATIONS
In this section we test the new condensation approach in FastChem Cond for three different scenarios of astrophysical interest.In the first case we replicate the equilibrium condensation calculations from Sharp & Huebner (1990).In the second scenario we evaluate the numerical stability of our new code by calculating the equilibrium condensation in a protoplanetary disk.Since the temperatures within the disk drop to values of roughly 10 K, this scenario presents a numerically very challenging case.In the last scenario we show the differences in the chemical composition between a standard equilibrium condensation calculation and the rainout approximation for the atmosphere of a brown dwarf.All results shown in the following can be easily reproduced by the Python scripts contained in the FastChem repository.

Comparison with Sharp & Huebner (1990)
The study by Sharp & Huebner (1990) focuses on the chemical composition of the gas phase under the influence of condensation in a general astrophysical context.They adapted the SOLGASMIX code (Eriksson 1971;Besmann 1977) to include the condensation of species for solar element abundances.In contrast to FastChem Cond, SOLGASMIX used a Gibbs free energy minimisation approach to solve the equilibrium condensation problem.The study presents the resulting gas phase composition as a function of temperature for a fixed pressure of 0.5 bar.
For this test calculation, we replicate the temperature and pressure conditions used by Sharp & Huebner (1990).The resulting abundances of selected gas phase species are shown in the plots of Fig. 1, which roughly mirror those presented by Sharp & Huebner (1990).We note, however, that we extended the temperature range of the chemistry calculations to even lower temperatures than those used by Sharp & Huebner (1990) to illustrate the numerical capabilities of our new code.
Comparing our results with those of Sharp & Huebner (1990) suggests a very good agreement in terms of the resulting molecular abundances.Minor changes can be explained by the different element abundances.Sharp & Huebner (1990) employed element abundances based on Cameron (1973), wheres FastChem Cond uses the ones from Asplund et al. (2009).Other notable differences are the different sets of gas phase species and condensates as well as the underlying thermochemical data, the fitted equilibrium constants are based on.
In addition to the gas phase, we also show the degree of condensation for several elements in Fig. 2. The results indicate that most species condense out completely below a certain temperature.Because oxygen is much more abundant than the other more refractory elements, only about 20% of O is condensed out until the lowest considered temperature of 500 K.For a solar element abundance mixture, oxygen would condense out almost completely at even lower temperatures in the form of water ice.
Examples for the condensation sequences of titanium and magnesium-bearing condensates are shown in Fig. 3. Titanium first condenses into the high-temperature condensate perovskite (CaTiO 3 (s)).At lower temperatures, Ti then switches to titanium oxide condensates, namely Ti 3 O 5 , Ti 4 O 7 , and TiO 2 .The transitions between the different Ti condensates are very sharp, suggesting that titanium is only present in a single condensate at a time.
In contrast, magnesium can be found in several stable condensates simultaneously as suggested by the lower panel of Fig. 3.It first condenses into spinel (MgAl 2 O 4 ), followed shortly afterwards by diopside (CaMgSi 2 O 6 ) and forsterite (Mg 2 SiO 4 ).Lastly, magnesium also condenses into enstatite (MgSiO 3 ), which competes with forsterite for the most abundant magnesium-bearing condensate.In contrast to the titanium condensation sequence, none of the magnesium condensates actually disappear after they have formed, though their abundances can change as a function of temperature.

Protoplanetary disk
In this section, we test our updated FastChem code on a temperaturepressure structure of the midplane of a protoplanetary disk.The temperature and pressure profile is based on the disk model described in Emsenhuber et al. ( 2021 14 orders of magnitudes, from about 10 −3 bar to 10 −17 bar.The temperature starts at roughly 1700 K near the central star and drops to values of almost 10 K in the outermost part of the disk.The low temperatures and densities, especially at larger distances, make this scenario numerically very challenging.In the outer parts of the disk, essentially all condensates considered in FastChem Cond have initially activities larger than unity, which makes it obviously quite difficult to obtain the final set of stable condensates according to the phase rule.We note that this calculation is only done to test the numerical stability of the chemistry and condensation scheme and not to properly model the chemical composition of a protoplanetary disk.In reality the disk chemistry is rather complex, involving various additional chemical and physical processes, such as photoevaporation for example (see, e.g.Henning & Semenov 2013).For the calculations here, we also removed germanium from the list of considered elements, because currently Ge has no associated gas phase molecules or condensates due to the lack of corresponding data in the JANAF tables (Chase 1998).
In total, this calculation included 26 elements.Thus, following the phase rule a maximum of 25 different condensates should be able to co-exist at most.As the lower panel of Fig. 5 shows, the highest number of stable condensates in our calculation reaches a value of 22 in the outer disk (see Table 1), clearly satisfying the phase rule.The number of stable condensates should also not be higher than the number of elements contained in them.In our case, the converged system also satisfies this requirement.
The upper panel of Fig. 5 shows the degree of condensation of selected elements as a function of distance (upper panel).Just like in the previous case of Sharp & Huebner (1990), elements tend to completely condense out below a certain temperature.Iron is lost from the gas phase very close to the inner disk boundary.The same also applies to other elements contained in high-temperature condensates, such as Ca or Al (not shown).More volatile elements, such as Mn, Na, K, or Cl start to condense at distances between 0.3 and 1 au.Oxygen is also contained in many of the high-temperatures condensates.Due to the high element abundance of O compared to the other elements only about 20% of the element is lost from the gas phase at small distances, though.However, once the temperaturepressure profile crosses the water ice line of the protoplanetary disk, oxygen condenses out completely into H 2 O(s) at distances larger than roughly 3.5 au.In the outer part of the disk, only hydrogen and the noble gases argon, neon, and helium remain in the gas phase.Hydrogen is also contained in some condensates, most notably water, ammonia, and methane ice.However, due to its very high element abundance, its degree of condensation never exceeds 0.2%.
Phosphorus is the only element that has a distinctly non-monotonic degree of condensation.With increasing distance it condenses out close to 1 au, returns to the gas phase near 3 au and finally condenses out again at  > 200 au.
Even though many refractory elements are almost completely bound in condensates from the inner to the outer regions of the disk, the actual condensates formed by them change with distance.Since many of the elements are found in multiple different condensates, they all compete with each other for the available elements.This is illustrated for a sample of magnesium, aluminum, and titaniumbearing condensates in the upper panel of Fig. 6.
With increasing distance from the disk's centre magnesium first condenses out into enstatite, forsterite, and spinel.Up to distances of 3 au, enstatite and forsterite compete for being the highest-abundance magnesium-bearing condensate.At 3 au, Mg(OH) 2 forms and takes up most of the available magnesium.Consequently, forsterite is no longer able to exist, while enstatite is still able to co-exist as it requires smaller amounts of magnesium to form.
Aluminum is first found predominantly in the form of spinel (MgAl 2 O 4 ).Near 0.5 au, it competes with NaAlSi 3 O 8 for the element silicon and finally disappears close to 2 au in favour of corundum (Al 2 O 3 ).The latter, though, is only stable for a short range of distances and is then replaced by Al 2 SiO 5 .This species can easily form there because the disappearance of forsterite in favour of Mg(OH) 2 provides an abundant amount of silicon.As more silicon becomes available, the abundance of corundum is decreased again by the formation of NaAlSi 3 O 8 , but still remains stable.Titanium is predominantly found in only two condensates: TiO 2 and MgTiO 3 .Unlike magnesium or aluminium it does apparently not form several different condensates simultaneously.The innermost stable Ti-compound is TiO 2 .It is replaced by MgTiO 3 near 1 au.However, with the formation of Mg(OH) 2 not enough magnesium is left and titanium returns to form TiO 2 .
As shown in the bottom panel of Fig. 6, at lower temperatures the volatile species H 2 O, CH 4 , and NH 3 condense beyond 3 au, 4 au, and 70 au, respectively.These condensing species effectively remove the elements O, C, and N from the gas phase.The last element to leave the gas phase is phosphorus.As already noted above, P has a nonmonotonic degree of condensation.Near 1 au, it first condenses out into Mg 3 P 2 O 8 .With the formation of Mg(OH) 2 , though, magnesium is no longer available.Since also no other phosphorus-bearing condensate has a high enough activity, phosphorus returns back to the gas phase.Only at much lower temperatures in the outer part of the disk, some magnesium becomes available again to form Mg 3 P 2 O 8 , removing P from the gas phase once more.
As discussed in Sect.2.3, the phase rule and the requirement that the resulting system of activity equations for the set of stable condensates needs to be linearly independent limits the number of stable condensates.The latter requirement, however, also suggests that each condensed species should be linked with a related element that is associated with its stability.Often, this should be the least abundant element contained in a condensate normalised with the corresponding stoichiometric coefficient (see also Appending B).
To illustrate this, we focus on the condensates present at the outermost grid point of the protoplanetary disk.According to Fig. 5, this point has the highest possible number of elements contained in condensates.All other elements left in the gas phase are the noble gases (He, Ne, and Ar) that don't have any condensate species associated with them in FastChem Cond.In Table 1 we list all 22 condensate species that are stable in the outer part of the disk.We also list the associated element for each condensate and the fraction of the fictitious condensate density   and its largest possible value  ,max (see Eq. ( 13)).If   is smaller than  ,max , the linked element is contained in more than one condensate.The assignment of the condensates to their respective elements is done by the procedure described in Appendix B.
The table clearly suggests that each element is uniquely linked with one specific condensate.Even though some elements, such as Mg or Al, can form several different condensates simultaneously they are the associated species in only a single condensate.In many cases, the associated element is also the least abundant one in a specific condensate, often containing the total available inventory of that element.This also explains why, for example, the titanium-bearing condensates form very distinct condensation sequences (see Figs. 3 & 6).Species like TiO 2 (s,l), Ti 2 O 5 (s,l), or MgTiO 3 (s,l) cannot exist simultaneously as Ti would be the associated and least abundant element in all of them.Magnesium, on the other hand, can also be present in other condensate species besides Mg(OH) 2 (s), such as MgSiO 3 (s,l) in this example here, because it cannot fully condense out in a single condensate due to its higher abundance compared to the other elements in these condensates.Thus, elements with higher element abundances will usually tend to form multiple different condensates if those also contain elements with smaller abundances.

Equilibrium condensation vs. rainout approximation
Finally, we apply the FastChem code to an atmosphere of a typical brown dwarf.As discussed in Sect.3.6, for atmospheres of planets and brown dwarfs the rainout approximation is a more suitable approach to calculate the chemical composition than the pure equilibrium condensation used in the previous two examples.In the following, we will show the impact of these two different approaches on the forming condensates and the chemical species remaining in the gas phase For the atmospheric temperature-pressure structure of the brown dwarf we use output from the Sonora grid of brown dwarf models published by Marley et al. (2021), available on Zenodo4 .In particular, we choose a model with an effective temperature of 750 K and a surface gravity of  = 3.16 × 10 4 cm s −2 (log  = 4.5), resembling a typical T5 dwarf.The corresponding temperature-pressure profile is shown in Fig. 7. Based on the general understanding of brown dwarf atmospheres, we expect the lower atmosphere to be dominated by the condensation of iron and various silicates, while in the upper atmosphere species like Na 2 S, NaCl, or KCl should appear (see Marley et al. (2013) or Morley et al. (2012), for example).
The temperature-pressure profile is used as an input for two different chemistry calculations with FastChem Cond, the first one using the standard equilibrium condensation approach from the two previous test cases and a second case with the rainout approximation.As described in Sect.3.6, for the latter one we start the chemistry calculations at the bottom of the atmosphere and proceed upwards towards lower pressure, while calculating the gas phase composition and the condensed phase, respectively.Once condensation is encountered, we use the resulting fictitious number densities   of the condensates to calculate the effective abundances   of these elements left in the gas phase.These new, reduced element abundances are then used in the subsequent calculations of the upper atmosphere.
In the top panel of Fig. 8 we show some important iron, calcium, and magnesium-bearing condensates comparing both condensation descriptions.Iron is one of the first condensates to become stable at a pressure of roughly 70 bar and a temperature of about 2200 K. Above this pressure no condensate currently included in FastChem Cond can stably exist because of the high temperatures at the bottom of the atmosphere.For the equilibrium condensation, iron is found in the form of Fe(s,l) in the lower atmosphere, while in the upper atmosphere the dominant iron-bearing condensates are fayalite (Fe 2 SiO 4 (s)), magnetite (Fe 3 O 4 (s)), and FeS(s).
In contrast to that, the condensation into Fe(s,l) at pressures of about 10 bar already depletes the gas phase iron element abundance in the upper atmosphere for the rainout approximation.Consequently, no other Fe compounds can form there as essentially no iron is left in the gas phase at pressures lower than about 4 bar.The rainout of iron into an Fe(s,l) condensate layer and, consequently, its absence in the upper atmosphere agrees with results previously described in studies such as Marley et al. (2021) or Visscher et al. (2010).
The corresponding results for major calcium and magnesiumbearing condensates are depicted in the top middle panel of Fig. 8.Here we can clearly see that in the equilibrium condensation case, calcium is first condensed in CaTiO 3 (s) as well as Ca 2 SiO 4 (s).Around 20 bar Ca is briefly mostly contained in CaSiO 3 (s), while the upper atmosphere is then dominated by CaMgSi 2 O 6 (s).
The latter two condensates do not form in the rainout approach, however.Here, Ca is depleted in the lower atmosphere by condensation into CaTiO 3 (s) and Ca 2 SiO 4 (s).As a result, more complex, secondary condensates like CaMgSi 2 O 6 (s) cannot form in the upper atmosphere.
Magnesium is mostly contained in the two condensates forsterite (Mg 2 SiO 4 (s,l)) and enstatite (MgSiO 3 (s,l)).Condensation of forsterite takes place at pressures below the 20 bar level in the equilibrium condensation calculation, followed by MgSiO 3 (s,l) between 10 bar and 0.1 bar.Both condensates co-exist in this pressure range, competing both for the available magnesium.Below pressures of 0.1 bar, however, MgSiO 3 (s,l) stops forming, while forsterite becomes the dominating magnesium-bearing condensate.The second-most abundant magnesium condensate assuming the equilibrium condensation approach in the upper atmosphere is Mg 3 P 2 O 8 (s,l).
In the rainout approach, on the other hand, enstatite does not form.Condensation of magnesium into forsterite near 20 bar already depletes the gas phase of magnesium in the upper atmosphere, such .Impact of condensation on the chemistry of a brown dwarf atmosphere.The top panel shows abundance profiles for selected condensates: iron-bearing condensates (left), magnesium and calcium species (middle), sodium and potassium condensates (right).Bottom panel: abundance profiles for silicon-bearing condensates (left), impact of the two different condensation approaches on the two gas-phase species hydrogen sulfide (H 2 S) and phosphine (PH 3 ) (middle), the degree of condensation for selected elements (right).The solid lines refer to the equilibrium condensation calculations.Dashed lines indicate the results for the rainout approximation.
that MgSiO 3 (s,l), as well as other Mg-bearing species are unable to form.
The top right panel of Fig. 8 shows the condensation of sodium and potassium-bearing species.In the equilibrium condensation calculation, Na first condenses into albite (NaAlSi 3 O 8 (s)) near a pressure of 9 bar and into halite (NaCl(s)) below 1 bar.Around 0.2 bar, albite is replaced by sodium silicate (Na 2 SiO 3 (s,l)) as the most abundant sodium-bearing condensate.Potassium is only found in two stable condensates: microcline (KAlSi 3 O 8 (s)) in the lower atmosphere and KCl(s) (sylvite) in the upper part.
In the rainout approach model, albite is unable to form because neither silicon nor aluminum are available as they have already been depleted in the deeper atmosphere.Consequently, sodium stays in the gas phase at lower pressures compared to the equilibrium condensation case and finally condenses out into Na 2 S at about 3 bar and NaCl near 1 bar.At roughly the same pressure level, potassium rains out into a sylvite layer without forming an additional condensate species.
The condensation sequence of silicon-bearing species is depicted in the bottom left panel of Fig. 8.As the figure suggests, this sequence is rather complicated, with silicon being present in as much as five different condensates simultaneously.In the lower atmosphere, Si is mostly contained in Ca-bearing condensates and SiO(s) whereas in the upper atmosphere it is predominantly chemically bound in forsterite, fayalite, diopside, and Na 2 SiO 3 (s,l).The absence of quartz (SiO 2 (s)) as a stable condensate in the equilibrium condensation approach is consistent with studies such as Visscher et al. (2010), though it has been predicted to form under non-equilibrium conditions (Helling & Woitke 2006).In the rainout approximation, silicon rains out already in the deeper atmosphere into three distinct condensates: SiO(s), larnite (Ca 2 SiO 4 (s)), and forsterite.
The two different condensation treatments do not only have an impact on the condensing species, but also on the chemical composition of the gas phase.Examples for two gas-phase molecules are shown in the bottom middle panel of Fig. 8.The most well-known example is the impact on H 2 S, as already extensively discussed in the literature (see Lodders & Fegley 2002, for example).As already discussed, in the equilibrium condensation approach, iron is mostly contained in FeS(s,l) in the upper atmosphere, depleting the gas phase of both iron and sulphur.As a result, the H 2 S molecule is unable to efficiently form in the gas phase in the upper atmosphere.
In contrast to the equilibrium condensation approach, we can see that H 2 S is abundantly present in case of the rainout approximation.This is caused by the aforementioned depletion of iron into Fe(s,l) in the lower atmosphere.As a result, FeS(s,l) is unable to form in the rainout case and sulphur remains available to form H 2 S.
Another gas-phase species strongly affected by the differences of the two condensation treatments is phosphine (PH 3 ) as suggested by the results shown in Fig. 8.In the equilibrium condensation approach, phosphorus condenses into Mg 3 P 2 O 8 (s,l) in the upper atmosphere, as depicted in Fig. 8.However, this condensate cannot be formed in the calculation with the rainout approximation because Mg is removed from the gas phase in the deeper atmosphere.Consequently, enough phosphorus remains in the gas phase to form phosphine.
The impact of the two different condensation descriptions can also seen in the form of the element depletion, depicted in the bottom right panel of Fig. 8.While the degree of condensation for elements in high-temperature condensates, such as Fe, Si, or Mg, is not strongly affected by the condensation treatment, larger differences can be noted in the aforementioned elements S and P. For the equilibrium condensation approach, both condense out at 1 bar and about 0.1 bar, respectively.In the rainout scenario, however, only about 10% of sulphur actually condenses into Na 2 S(s,l), while the degree of condensation for phosphorous remains zero throughout the entire atmosphere.Other elements, such as sodium or potassium, condense out completely in both condensation approaches, albeit they leave the gas phase at considerable lower pressures for the rainout approximation.Conversely, chlorine already rains out at higher pressures than in the equilibrium condensation calculation.

SUMMARY
In this study we introduce a new version of FastChem.With a version number of 3.0, we refer to this updated code as FastChem Cond.While the two previous versions of FastChem only accounted for the gas phase composition, the new version now introduces the treatment of condensates to the code.In particular, FastChem Cond can compute scenarios involving equilibrium condensation as well as the rainout approximation often used in the modelling of atmospheres of (exo)planets and brown dwarfs.Together with the code update we also add about 290 liquid and solid condensate species to FastChem.
The numerical treatment of condensates is based on the basic ideas presented by Leal et al. (2016) and have been successfully adapted to the numerical formalism used in FastChem.A major advantage compared to other established equilibrium condensation codes in the field of astrophysics is the automatic selection of stable condensates satisfying the phase rule.
The updated version of FastChem is tested for several scenarios involving condensation.In particular, we reproduce the equilibrium condensation calculations by Sharp & Huebner (1990) and also test FastChem Cond for the numerically very challenging scenario of the midplane of a protoplanetary disk.Finally, we use a theoretical atmosphere model of a typical T5 brown dwarf to showcase the differences between equilibrium condensation and the rainout approximation.
FastChem Cond is programmed in object-oriented C++ and offers the additional Python module pyFastChem that allows it to be run directly from any Python script.The code is released as open source at GitHub (https://github.com/exoclime/FastChem)under the GNU General Public License version 3 (gnu.org2007).The repository also contains several Python example scripts that showcase the use of FastChem Cond under Python, as well as an extensive manual that provides detailed descriptions of the code, its required inputs and available, optional parameters.and usually need to be converted to the reference state used in our formalism.
If the number densities   are given in units of cm −3 , then in the activity equation ( 14) for the condensates   has units of It is related to the dimensionless mass action constant K () by ln   () = ln K () − ln where  − • is the pressure of the standard state of the thermochemical data.We note that this relation explicitly makes use of the ideal gas law.For the data used in the FastChem code,  − • is equal to 1 bar.As stated by Stock et al. (2018), we parameterise ln K as a function of the temperature  in form of the polynomial where   and   are the corresponding fit coefficients that are tabulated in the FastChem input file for each species.Following the approach of the gas-phase data in the previous versions of FastChem, we base our choice of species mainly on the JANAF tables (Chase 1998).In particular, we use all condensates for the elements contained in FastChem included in the JANAF tables.Species that have both liquid and solid-form data are fitted independently.Thus, some condensate species have two different sets of polynomial coefficients attached to them.FastChem Cond in principle allows a condensate with a certain stoichiometric formula to have even more sets of fit coefficients.It is, thus, possible to fit also different crystal phases independently if so desired by the user.A list of all species and the fit coefficients can be found in Table A1.
In addition to the JANAF tables, we also add selected condensates from other data sources.This includes condensate species from Sharp & Huebner (1990) not available in the JANAF tables, some volatile, low-temperature condensates such as water, ammonia, methane, or carbon dioxide, as well as silicon monoxide (SiO(s)) from Gail et al. (2013).
For species where the thermochemical data is given in terms of the saturation vapour pressure  vap, () rather than the Gibbs free energy, we use the relation , where  − •  () is the Gibbs free energy of the corresponding molecule.For the latter we again use the gas-phase data from the JANAF tables.
We note that we exchanged the data for SiO 2 (s,l) from the JANAF tables with that from Barin (1995).The JANAF data suggests that SiO 2 (l) would switch into solid SiO 2 , more specifically quartz, at a temperature of 1696 K.However, the standard phase diagram of SiO 2 as well as the data from Barin (1995) or the NASA Glenn coefficients (McBride et al. 2002) all agree that at standard pressure, SiO 2 (l) first transitions into the solid cristobalite at a temperature of 1996 K before its stable form finally changes to quartz at lower temperatures.Due to the discrepancies between the JANAF tables and the other literature, we choose to use the data from Barin (1995) instead.
Thermochemical data for condensates are often only available over a very restricted temperature range.Therefore, by default FastChem Cond avoids extrapolating the fit equation (A6) outside of the temperature range of the underlying data.A special configuration parameter is available in FastChem Cond, though, that allows the user to override this limitation.However, this parameter should only be used with great care as the extrapolation might result in an unphysical, non-monotonic behaviour of the ln   fit for some condensate species.

APPENDIX B: ELEMENT-CONDENSATE MAPPING
In order to determine the correlation between the condensed species with their associated elements even at the highest number of possible stable condensates according to the phase rule one can apply a simple scheme that is briefly discussed in the following.The scheme is applied to the outermost region of the protoplanetary disk from section 4.2 and its results are summarised in Table 1.
Figure B1 depicts the end-result of how a unique mapping between the chemical elements and their corresponding condensates according to the phase rule can be found.The right hand side of the figure defines a matrix, which entries are "1", if the condensate contains the respective element.The number in the outermost left column on the left hand side of the figure indicates how many condensates an element can be associated with.The elements are then sorted by this number, e.g.carbon receives a "1", because the only condensate containing carbon is CH 4 , nitrogen receives a "2", because NH 4 Cl and NH 3 include nitrogen, and so on.The mapping of the elements C, Ca, . . ., Zn is already determined, since they have only one condensate to associate with.These condensates are therefore not available to other elements.In the next step, the numbers on the left hand side in Fig. B1 are updated accordingly (second column).The procedure is repeated until no condensate is left.
In our example the final condensate is H 2 O which can be assigned to either hydrogen or oxygen, since the total number of elements exceeds the total number of condensates by one.If a bĳective mapping is desired, a secondary criterion must be utilised, such as for example the stoichiometrically normalised element abundance.In this case, we would assign oxygen, the less abundant element, to H 2 O. Note, that the domain of the mapping must then be reduced by hydrogen which is left over.The described algorithm works because the condensates are linear independent (see sections 2.3 and 3).
Table A1.An overview of condensates and the coefficients for the equilibrium constants included in FastChem.Species with both liquid and solid phases have two separate sets of coefficients.The first line always refers to the solid species, while the second line represents the liquid phase.It is important to note that the tabulated ln K use the elements in their monatomic form as reference state (see text for details).Moses et al. (1992), Prydz & Goodwin (1972), NIST Chemistry Webbook; (3) Goodwin (1985); (4) Yaws (1999) 10) Barin (1995);

Figure 1 .
Figure 1.Comparison with equilibrium condensation calculations ofSharp & Huebner (1990).The plots show the partial pressures of gas-phase species as a function of temperature and roughly mirror the results shown bySharp & Huebner (1990) with some species moved to other panels for better visualisation.We note that that our calculations extend to lower temperatures than those presented bySharp & Huebner (1990).Following the corresponding figures bySharp & Huebner (1990), solid lines refer to the chemistry calculations with equilibrium condensation while dashed lines refer to those without.

Figure 2 .
Figure2.Degree of condensation for selected elements for theSharp & Huebner (1990) scenario as a function of temperature and a gas pressure of  g = 0.5 bar.

Figure 3 .Figure 4 .
Figure 3. Examples for condensation sequences of selected elements for the Sharp & Huebner (1990) scenario.The upper panel shows titanium-bearing condensates, while the lower panel depicts the condensation sequence for magnesium.

Figure 5 .
Figure5.Degree of condensation of selected elements (upper panel) and the number of stable condensates (lower panel) as a function of the distance from the central star for the protoplanetary disk.The lower panel also additionally shows the total number of elements contained in stable condensates to verify that the results satisfy the phase rule.

Figure 6 .
Figure 6.Abundance profiles of selected condensates as a function of distance in the protoplanetary disk model.The upper panel focuses on Mg, Ti, and Al-bearing refractory condensates, while the lower panel shows the more volatile condensate species.

Figure 7 .
Figure 7. Temperature-pressure profile of a typical T5 brown dwarf used for the chemistry calculations.The profile is taken from the grid of Marley et al. (2021) for an effective temperature of 750 K and a surface gravity of  = 3.16 × 10 4 cm s −2 (log  = 4.5) .
Figure8.Impact of condensation on the chemistry of a brown dwarf atmosphere.The top panel shows abundance profiles for selected condensates: iron-bearing condensates (left), magnesium and calcium species (middle), sodium and potassium condensates (right).Bottom panel: abundance profiles for silicon-bearing condensates (left), impact of the two different condensation approaches on the two gas-phase species hydrogen sulfide (H 2 S) and phosphine (PH 3 ) (middle), the degree of condensation for selected elements (right).The solid lines refer to the equilibrium condensation calculations.Dashed lines indicate the results for the rainout approximation.

Table 1 .
List of stable condensates in the outer parts of the protoplanetary disk with the associated elements according to the phase rule.

Table A1 -
continued An overview of and the fit coefficients for the equilibrium constants included in FastChem.A1continuedAn overview of condensates and the fit coefficients for the equilibrium constants included in FastChem.