© ESO 2021

1 Introduction

Cloud formation has emerged as one of the major obstacles to fully benefitting from observations of exoplanetary and brown dwarf atmospheres and efforts are therefore gearing up to understand fundamental processes in enough detail to allow the development of efficient model descriptions (Helling 2020). Exoplanet cloud modelling approaches include a variety of parametrizations from imposing cloud particle properties as an opacity source (e.g. Roman et al. 2021; Parmentier et al. 2020) and adding gravitational settling and mixing (Lines et al. 2019), to kinetic modelling (Lee et al. 2016; Lines et al. 2018)as part of radiation-hydrodynamics atmosphere simulations. The data driven retrieval is also hampered by the need fora fast representation of clouds as opacity sources (e.g. Barstow 2020; Colón et al. 2020; Lacy & Burrows 2020) and it is utilised routinely for atmospheric abundance measurements (e.g. Nikolov et al. 2018; Spake et al. 2020).

Modelling the formation of cloud condensation nuclei (CCNs) is key to modelling cloud formation in brown dwarf and exoplanet atmospheres because it determines when and where cloud formation sets in. Given that the formation of exoplanet CCNs is largely unexplored, Ormel & Min (2019) introduced the nucleation rate (CCN formation rate) as a parameter into their physical model for cloud formation and transport, and they demonstrate that the nucleation rate is rather ill constrained by present observations for giant gas planets. On Earth, the CCN rate can be adjusted according to detailed observations of sand storms, bush fires, salt spray, and other factors. Nevertheless, the effect of CCN on cloud formation is the largest source of uncertainty in climate modelling on Earth, as described by the Intergovernmental Panel on Climate Change (IPCC) (Stocker et al. 2013).

The James Webb Space Telescope (JWST) guaranteed-time observer (GTO) target list contains a substantial number of gas giants, including HD 189733b, WASP-43b, and WASP-121b, in addition to the rocky planets of the TRAPPIST-1 system, for example. In order for cloud formation to occur in gas giants, brown dwarfs, and hot rocky planets, cloud condensation nuclei need to form from the gas phase. The quest for the primary condensate is not new and leads back to the Asymptotic Giant Branch (AGB) star dust formation modelling (Sedlmayr 1994). Withing the research of AGB star mass loss, this idea has recently been re-addressed by Boulangier et al. (2019). Using the modified classical nucleation theory with surface tension values for a wide variety of materials, Lee et al. (2018) suggest that a considerable range of CCN materials may be available to kick-off cloud formation in exoplanet and brown dwarf atmospheres. However, little information exists about their formation mechanisms in astrophysical environments. Works by Patzer et al. (1999, 2005), Lam et al. (2015), Decin et al. (2017) and Gobrecht et al. (2018) explore the formation of corundum (Al₂O₃ solid) through the study of Al-O clusters, and Boulangier et al. (2019) propose the Al₂O₃ solid as the most efficient condensation seed in AGB stars. Chang et al. (2013) explore the formation of Fe–C–H ring molecules (FeC₂H_n and FeC₃H_n for n = 0, 2, 4) while Patzer et al. (2014) explore Ti–C clusters (Ti_x C_y, x, y = 1…4). Chang et al. (2005), considering the neutral and charged inorganic cage molecules Kr@Y₁₂@Z (Y = Ni, Pd; Z = As, Sb, Bi; q = 0,−1,−3), report a stable system where a krypton atom is enclosed by a fullerene-like inorganic double cage. Only the exploration of the Ti–O clusters by Jeong et al. (2000) followed by Lee et al. (2015) led to the application of the cluster data to exoplanet cloud formation studies as part of kinetic cloud formation modelling by Woitke & Helling (2004), Helling & Woitke (2006), Helling et al. (2019) and Samra et al. (2020). The most exhaustive study of the formation process of CCN has been conducted for the Earth’s atmosphere, where simulations can be paired with laboratory works and in-situ measurements for the H₂SO₄-complex (Svensmark et al. 2013, 2016, 2017, 2020; Dunne et al. 2016). Recently, a 3D particle Monte Carlo (3D MC) approach has been developed to model the early stages of CCN formation for the Earth’s atmosphere (Köhn et al. 2018, 2020). This model, which tracks single molecules and clusters in a 3D space, captures fluctuations that are not included in most models but are essential for the formation of stable clusters (Olenius et al. 2018).

This paper focuses on the study of seed formation in exoplanet and brown dwarf atmospheres, and in particular on the formation of TiO₂ nucleation. TiO₂ has long been debated as a primary seed forming species (Sedlmayr 1994; Helling & Fomins 2013) and we aim to explore the formation of TiO₂ -CCN by two different methods and compare the results. Here we couple methods developed within astrophysics with models developed for terrestrial atmospheric science in order to study the growth of TiO₂ clusters. The two methods are the 3D particle Monte Carlo approach where we expand the approach presented by Köhn et al. (2018), and the kinetic approach that solves a set of rate equations that require thermodynamic data for the individual clusters (Sedlmayr 1994). We use the Monte Carlo approach to gain more in-depth understanding about the onset of the nucleation process, the transition from a time-dependent to a stationary process, as well as the effect of reaction efficiencies. We are particularly interested in the time evolution of individual and cumulative cluster nucleation rates in 3D space.

Section 2 outlines our approach and the required theoretical background for both methods applied. Section 3 presents the results on the TiO₂-cluster formation, their individual abundances, size distributions, and the effect of the reaction efficiency α. Section 4 presents the comparison of the results from both methods, MC and rate equations approach. We discuss our results in Sect. 5 and conclude with Sect. 6.

2 Approach

We approach the modelling of the formation of TiO₂ cloud condensation nuclei (CCNs) by applying two different methods, a 3D MC approach (Sect. 2.2) and a steady-stateapproach utilizing thermodynamic cluster data (Sect. 2.3) in order to understand TiO₂ cluster formation as a necessary step to the formation of TiO₂ nucleation seeds. The MC approach enables the study of the formation dynamics of the cluster growth, which is not accessible by a steady-state approach. Each method makes assumptions and requires certain input data that we outline in the following.

2.1 Terminology

As this work combines expertise from astrophysics and Earth science, there can be a confusion due to the difference in how various terms are applied. Here we explain how we use certain terms. A “cluster” is a collection of monomers sticking together through physical or chemical interaction. “Nucleation” is the formation of a thermodynamically stable cluster (i.e. crossing the critical size), the genesis of clusters of all other sizes is referred to as formation. “Cloud condensation nuclei” (CCN) are clusters large enough that they can be activated by another gas with a given supersaturation to form cloud drops (on Earth this gas is water). “Growth efficiency” and “sticking probability” is used interchangeably to indicate the probability that a collision between clusters leads to growth. “Reaction” refers to all processes changing the cluster size.

2.2 Three-dimensional Monte Carlo simulation of TiO₂ seed formation

We use the particle 3D MC code introduced by Köhn et al. (2018). This code tracks individual particles in a 3D space and in time, and hence allows a detailed study of the nucleation process, including rare events and fluctuations for individual cluster sizes. The computational time is large and it scales with the number of particles considered.

Monte Carlo particles: The particles that are considered in the 3D MC simulations are characterized by size, composition, and their position in Cartesian coordinates. TiO₂ monomers and clusters can be traced. The largest cluster size achieved in the present 3D MC simulations is N > 10⁵ monomer units per cluster for a simulation covering t = 10⁹s. The TiO₂ monomer mass and radius are kg and nm, respectively. The definition of cluster radii for sizes larger than the monomer (N > 1) is challenging, and we note that larger clusters can be more compact than smaller clusters with less elongated structures. This is explored in Fig. 1 where the results from different methods to calculate a cluster diameter are shown. To allow for comparison, the interacting volume of each cluster is calculated with the respective method and the radius of a sphere with the equivalent volume is plotted. The interacting volume of each cluster is larger than the geometrical volume assuming spherical particles because of van der Waals forces determining the actual chemical interaction. The interacting volume is calculated by finding the axes of the smallest box that contains the cluster of size N and adding the Van der Waals interaction length of 1.6 Å (Koch & Manzhos 2017) to each of the axes. The volume of that extended box is the interacting volume of the cluster. For a different method, plotted in orange, the volume was calculated using Delaunay triangulation and the ConvexHull algorithm as described by Barber et al. (1996). This algorithm uses polygons to approximate the surface of the cluster and therefore gives a smaller volume than the approximation through a box. After determining the sphere equivalent radius from this volume, the Van der Waals interaction length is added to account for the interaction distance that is larger than the physical dimensions spanned by the atomic cores of the cluster. In our 3D MC simulations, the cluster collisions are treated as a hit-and-stick coagulation process between clusters, hence contractions, compactification, or expansion are omitted.

Monte Carlo cluster growth: The frequency with which cluster collisions may occur depends on the ambient density. Figure 2 demonstrates that exoplanet atmospheres that fall into the thermodynamic regime of giant gas planets shown here are mainly in the hydrodynamic regime of a free molecular flow. Hence, the collision frequencies of the clusters in our 3D MC model will be determined by their mean free path. We note that diffusion is commonly assumed to determine the collision frequencies in molecular-dynamic simulations for Earth atmosphere nucleation studies (see e.g. Köhn et al. 2018, 2020), which is supported by a Knudsen number of 10⁻¹⁵ for Earth’s atmosphere (and hence as an example for rocky planets).

The calculation of the time interval, Δt, between twopossibly constructive cluster collisions is therefore guided by an analysis of the Knudsen number (Fig. 2). The Knudsen number (Kn = l_MFP∕(2a); l_MFP – mean free path, a – particle radius) is derived for results from DRIFT-Phoenix model atmosphere simulations (Witte et al. 2009, 2011). Results include the hydrostatic (T_gas, p_gas)-structures, which are calculated consistently with cloud formation, chemical equilibrium, and gas and cloud radiative transfer. As a result, the cloud structure can be described by mean particle size, a, material volume fraction, and element depletion. The maximum MC time step is therefore defined as (1)

with the local thermal velocity v_rel(T), and the mean free path l_MFP within the gas (see e.g. Eq. (10) in Woitke & Helling 2003). Hence for temperatures between 500 K and 2000 K, max(Δt) lies on the order of 10⁷ s; for better numerical resolution, we have chosen Δt = 100 s. The clusters move with their local thermal velocity for a time Δt. If their location overlaps after Δt, they are considered to merge. The merging can occur between clusters of different sizes and a monomer only (polymer-monomer) or between all available sizes (polymer-polymer). Such a merging of two MC particles depends on their collision probability but also on the stickiness (i.e. reactiveness) of the two interacting clusters. This is parameterised as a sticking probability, α. The effect of the sticking probability on the cluster growth is tested by simulating different sticking probabilities α ∈ [0, 1] and checking for r ≤ α with a random number r ∈ [0, 1) once two particles overlap. Only if this condition is fulfilled are two MC particles merged (coagulated in Earth science terms) by adding their volumes and their masses. After coagulation, a new monomer is added if the total number of individual particles is less than 10³ such that d f_tot∕dt ~ 0 for the total particle density f_tot. By adding monomers to the simulation domain, we simulate a constant influx of monomers whilst we keep the total particle density constant. As such, we take into account that cloud formation requires the transport, diffusion, and mixing of particles.

Besides the short-term interaction due to the collision of TiO₂ clusters, long-term interactions are caused by the fact that TiO₂ clusters are dipoles. TiO₂ has a dipole moment of approximately 7D (Brünken et al. 2008) where D ≈ 3.33 × 10⁻³⁰ C m is the Debye unit of the electric dipole moment. Small TiO₂ clusters show dipole moments with similar values to the TiO₂ monomer (e.g. (TiO₂)₃ and (TiO₂)₄), or no dipole moment like the symmetric dimer (TiO₂)₂. However, the dipole contribution to the velocity is on the order of v_dip ~ 2pΔt∕(r³ ⋅ m) ≈ 10⁻²⁶ m s⁻¹ for r = l_MFP,  Δt = 10⁷ s and , which is much smaller than the thermal velocity that is on the order of several hundreds of m s⁻¹; hence the dipole contribution is negligible.

Monte Carlo cluster evaporation: Clusters not only undergo constructive growth processes but can also evaporate if they are thermally unstable. The most unstable cluster is called the critical cluster, N_*, and all clusters larger than N_* will preferably grow and therefore be thermally stable. After every time step Δt, we check the probability of evaporation,γ ⋅ Δt, for particles that consist of more than one monomer. If for a random number r ∈ [0, 1), evaporation is simulated by emitting a single monomer randomly. We note that clusters are not arranged as a collection of smaller clusters but as a collection of monomers; hence the probability of a polymer escaping is quite small and we only model the emission of single monomers from large clusters.

The evaporation frequency, γ, for clusters of size R(N) and mass m(N) is described according to Yu (2005) as (2)

where is the phase-equilibrium vapour density for TiO₂ with respect to the solid TiO₂[s] phase (Woitke et al. 2018), the surface tension σ (Lee et al. 2015), and the density ϱ(N, T) of the clusters of size N are (3) (4) (5)

with a monomer radius of nm and a monomer mass of kg. Boltzmann’s constant is k_B ≈ 1.38 × 10⁻²³ J K⁻¹, T is the temperature, M₁ ≈ 80 g mol⁻¹ is the molar mass of TiO₂, and J (mol k)⁻¹ is the universal gas constant. Thus, the formation of clusters is determined by the relation between the collision and evaporation frequencies, which are both increasing functions of temperature. Using the evaporation rate (2) is the standard approach to include evaporation (Seinfeld & Pandis 2006; Yu 2005) as experimental evaporation rates are hard to obtain. For H₂SO₄–H₂O clusters, we have validated this approach against experimental results (Köhn et al. 2018).

Figure 1 (top panel) compares the values used in the 3D MC simulation (green line, Eq. (5)) to other representations of cluster diameters. Numerical checks (Fig. 1, lower panel) of the effect of these cluster radius differences showed little effect of these uncertainties on the results of the present study; hence the assumption of spherical particles works sufficiently well although small particles deviate from a spherical form.

By applying the surface tension we assume our clusters to have a non-planar surface. The Kelvin effect (i.e. the curvature) is taken into account through the exponential term in Eq. (2) (Kelvin equation), which models the effect of the cluster surface curvature. The same effect is accounted for in the rate equation approach (Sect. 2.3) when the surface energy of cluster N is expressed as Eq. (4.14) in Helling & Fomins (2013).

Figure 3 shows the evaporation frequency (Eq. (2)) for TiO₂ clusters as a function of constituting monomers N for different temperatures. The strong temperature dependence suggests that evaporation will be negligible for T ≤ 1000 K whilst evaporation becomes dominant for T ≥ 2000 K.

Monte Carlo cluster formation rates: The time evolution of the formation rate of clusters larger than size N is derived from our 3DMC results as (6)

where f_≥N(t) is the number density of clusters above a given size N as a function of time t within the volume V (1 cm⁻³ for the present simulations). Thus, the number of clusters of size N is f_{= N}(t) ⋅ V (see Fig. 4).

Fig. 1 Top: different cluster radius definitions. Blue line: radius based on the maximal inter-atomic distance from DFT (density functional theory) simulations (Berardo et al. 2014). Orange line: Radius based on the ConvexHull volume approximation of cluster geometry from Berardo et al. (2014) with added Van der Waals interaction distance. Green line: (Eq. (5)). Purple line: radius based on volume of smallest box containing cluster with added Van der Walls interaction distance. Bottom: Cluster number densities from MC code using different radius definitions. Solid lines: (Eq. (5)); dotted lines: box approximation/van-der-Waals radius (purple line in upper panel).

Fig. 2 Knudsen number, Kn = lMFP∕(2a) (lMFP – mean free path, a – particle radius) for different exoplanet and brown dwarf DRIFT-Phoenix model atmosphere results, except for Jupiter and Venus. The number Kn indicates whether the cluster interactions occur through diffusion (Kn < 1) or in a free molecular flow (Kn > 1). Exoplanet and brown dwarf atmospheres are predominately in the regimes of a free molecular flow (Kn > 1). Knudsen numbers for Earth’s atmosphere and typical cloud particle sizes are ≈ 10−5 and therefore fall in the turbulent, viscous case.

Fig. 3 Evaporation frequency γ(N) (Eq. (2)) for TiO2 clusters as a function of their size N for different gas temperatures.

2.3 Kinetic seed formation through rate equations

The kinetic, steady-state approach describes the temporal evolution of the size distribution function as (Helling & Fomins 2013) (7)

where f(N) is the number density of a molecular cluster containing N i-mers contributing to the growth of the particle and is the effective flux (or transition rate) for the growth of the particle of size N − i to size N due to all associating or dissociating reactions involving an i-mer of the condensing species. This stationary flux through cluster space is (8)

summing over all chemical reactions r_i in which an i-mer is involved. Equation (8) is equivalent to Eq. (15) in Köhn et al. (2018), both giving the number of seed particles cm⁻³ s⁻¹. The growth time τ_gr(r_i, N − i, t) is the time by reaction r_i leading from cluster size N − i to cluster size N. The evaporation time τ_ev(r_i, N, t) is the time leading from size N to size N − i, (9) (10)

The densities n_f(r_i) and n_r(r_i) are the number density of the molecule of the growth (forward) process and of the evaporation (reverse) process for reaction r_i, respectively (Patzer et al. 1998). The surface A(N) is the surface of the cluster of size N. The average relative velocity between the growing or evaporating TiO₂ molecule and the cluster is v_rel, and is defined as the thermal velocity (11)

The reaction efficiencies for growth and evaporation via reaction r_i are given through α(r_i, N − i) and β(r_i, N). We note that α can also be referred to as the sticking probability. Equation (10) corresponds to Eq. (2) in the 3D MC approach, except that the thermal stability does not enter Eq. (10). Instead, we use the thermodynamic properties of the clusters according to Eq. (12). The challenge is that the growth and evaporation efficiency coefficients, α(r_i, N) and β(r_i, N), are often unknown for the different cluster sizes N, which is why we perform 3D MC simulations in addition to the kinetic approach. As Eq. (10) corresponds to Eq. (2), hence , we can estimate β = γ∕(A(N)v_rel(n_r(r_i), N) n_r(r_i)). For T = 1000 K, β is on the order of ≈ 10⁻¹–10⁰ for all considered sizes. Hence, we primarily use β = 1 when comparing the MC approach and the kinetic model.

For further considerations, we need to introduce a reference equilibrium state. In that, we follow Patzer et al. (1998) who show in their Appendix A that if the temperatures of all components are equal, the supersaturation ratio of a cluster of size N with respect to the bulk is . Hence, phase equilibrium between monomers and between the clusters and the bulk solid, plus simultaneous chemical equilibrium in the gas phase, plus thermal equilibrium (i.e. all components have the same temperature) characterise this equilibrium state. In such local thermodynamic equilibrium (LTE) between the gas phase and the clusters, where is the equilibrium number density, the principle of detailed balance holds for a single microscopic growth process and its respective reverse, that is, the evaporation process. This implies that under the condition of detailed balance, , which allows us to express the evaporation rate by the growth rate. The equilibrium number densities for the clusters and the monomers are , , , and . The law of mass action links these equilibrium particle densities to Gibbs free energies, (12)

The energy (kJ mol⁻¹) is the Gibbs free energy of formation and it can be calculated from the standard molar Gibbs free energy of formation of all reaction participants at the temperature T_d(N) = T_gas, which is (13)

Assuming LTE, the equilibrium cluster size distribution is therefore expressed by a Boltzmann distribution, (14)

with being the equilibrium density of the monomer (here TiO₂).

The difference in energy ΔG(N) is the free energy change due to the formation of a cluster of size N from the saturated vapour. It is related to the standard molar Gibbs free energy of formation of the N-cluster by (15)

The energy is the standard molar Gibbs free energy of the formation of the solid phase, and p^{° −} is the pressure of the standard state. Most often, p^{° −} is the atmospheric pressure on the Earth at which p_sat and were measured. The right hand side of Eq. (15) now contains quantities that can be determined from lab experiments or quantum-chemical calculations.

No assumptions were necessary with respect to size-independent surface tensions as required in classical nucleation theory, which is the droplet model applied in Sect. 2.2. We note that only a classical nucleation approach (droplet model) requires the consideration of the so-called Kelvin effect. An exponential (Eq. (2), also Eq. (18) in Powell et al. 2018) corrects the saturation vapour pressure, p_vap, data that is derived for flat surfaces. The Kelvin effect is therefore used to include some representation of the curvature of the clusters for the evaporation process. The use of thermodynamic cluster data in combination with the assumption of detailed balance conveniently side-steps calculating evaporation such that the concept of surface tension is not required to describe the effect of surface curvature on thermal stability.

3 Three-dimensional Monte Carlo dynamics: (TiO₂)_N cluster formation results

The results from the 3D MC simulations for TiO₂ cluster formation are presented for a set of local gas temperatures of T = 500 K, 1000 K, 1500 K, and 2000 K and a local total gas density of f_{tot, gas} = 10³ cm⁻³ of these TiO₂ clusters. These thermodynamic values are representative of the region where cloud particle seed formation may occur in exoplanet and brown dwarf atmospheres (e.g. Lee et al. 2015; Helling et al. 2017). The 3D MC simulations are run for a fixed volume of V = 1 cm³ as the computational domain in which the total number of particles remains constant for each run, namely 10³ particles per cm³. By assuring f_{tot, gas} = const, each simulation is run under isobaric conditions such that T = const. Our studies are therefore not affected by changing thermodynamic conditions and all 3D MC results presented here are caused by particle interaction alone.

For each temperature, we explore two cases for cluster growth: a) cluster–monomer growth (i.e. polymer–monomer), and b) cluster–cluster growth (i.e. polymer–polymer) where cluster growth can proceed by inter-cluster coagulations of all sizes. The growth tree is visualised in Fig. A.1 to help understand the presented results. In order to address the effect of parameter uncertainty, different growth efficiencies (sticking probabilities) α = 1.0, 0.1 and a random value rand were tested.

3.1 Cluster size distribution in a homogeneous TiO₂ gas: the effect of temperature and reaction efficiency

Figure 4 shows the (TiO₂)_N cluster size distribution for different temperatures and different reaction efficiencies (sticking probabilities, α) in a pure TiO₂ gas. The top panel shows the distributions for N = 1–50 and α = 1 for monomer–cluster (denoted as polymer–monomer) and cluster–cluster (polymer–polymer) growth scenarios at three different times. The middle panel shows the full size distribution for all six scenarios at 1000 K (three different values of α for monomer–cluster and cluster–cluster growth) at the same time (t = 4 × 10⁸ s), while the bottom panel shows cluster size distributions at temperatures of 500, 1000, and 1500 K for both monomer–cluster and cluster–cluster scenarios forα = 1 and α = rand at the same time (t = 1.7 × 10⁷ s). These times are chosen as the highest possible time that is shared between all the shown simulations.

In the top panel, we see that the monomer–cluster scenario and the cluster–cluster scenario develop in two fundamentally different ways. The former has a Gaussian-like shape with a well-defined mean cluster diameter while the latter looks more like a power-law that reaches larger sizes than the monomer–cluster cases but also has more free monomers. When only growth via monomers is allowed, an abundance of dimers is formed as can be seen for t = 3 × 10⁷ s. Due to the replenishment mechanism for monomers, a new monomer is created for each collision (as long as the total number of clusters is below the starting amount of 1000). So, the formation of a dimer effectively consumes (at least) one monomer. This depletes the amount of monomers and makes the resulting distribution look more like the growth of an initial burst of newly formed aerosols. For the cluster–cluster cases, more large clusters are formed and these are effective at scavenging the monomers that are then replenished (nearly) 1:1. Due to the replenishing mechanism, the total cluster mass increases during the simulations at different rates for the two scenarios. In the α = 1 cases for 1000 K after 4 × 10⁸ s seconds, there is a total amount of monomers (free + bound in clusters) of 55 406 for the cluster–cluster simulation and of 8372 for the monomer–cluster simulation. We note that these are monomers and there are never more than 1000 clusters in the simulation; however, since we add monomers to the simulation domain, whenever the total particle number (monomers and clusters) declines below 1000, the total number of free and bound monomers can reach values above 1000 monomers.

The middle panel allows a comparison of the different sticking probabilities. At the shown time, (t = 4 × 10⁸ s), the distributions for all three cluster–cluster scenarios are very similar, indicating that they have reached steady-state conditions except for the very largest clusters that still continue to grow. For the monomer–cluster scenarios, the α = 0.1 scenario is less developed than the two others, which is to be expected if a steady state has not yet been reached. Since only growth via monomers is expected, it makes sense that this scenario develops more slowly than the cluster–cluster cases.

In the bottom panel, we see the effect of temperature. For α = 1 and monomer-polymer collisions very little has happened for 500 K whereas the 1000 K distribution has grown to larger sizes. Higher temperatures cause more collisions, which produces larger clusters. The trend of increasing concentrations going from 500 K to 1000 K is similar for α = rand with cluster–cluster collisions allowed. If the temperature is increased further to 1500 K, the distributions shift towards lower sizes again in both shown scenarios. This is because the evaporation increases drastically (see Fig. 3); for T = 1000 K the evaporation probability is 10⁻⁹ to 10⁻¹⁰ s⁻¹ for small clusters while for T = 1500 K it has increased to 1 s⁻¹, yielding ≈ 1000 K as the most efficient temperature for TiO₂ nucleation (Boulangier et al. 2019). No size distributions for T = 2000 K are shown since very little happens at that temperature because of the very large evaporation frequency.

Thus the formation of clusters is determined by the relation between the collision and evaporation frequencies, which are both increasing functions of temperature. If the temperature is too low nothing happens because there are too few collisions, but if it becomes too high then evaporations dominate. This is consistent with results found by Lee et al. (2015).

Figure 5 elucidates the interplay between collisions and evaporations by showing the timescale for both processes as a function of cluster size (N) for T = 1000 K, α = 1, with and without cluster–cluster growth. These are the timescales observed in the actual simulations which, for evaporation, can be compared with the evaporation frequency calculated in Fig. 3. The waythe characteristic times are calculated is by looking for how long it takes for a given size N to increase or decrease and averaging this over a time interval (see the figures for the exact times). Therefore, any increase in for example N = 2 is counted as growth, but it might as well be evaporation from N = 3.

We see that the times are very close to each other for the monomer–cluster and cluster–cluster scenarios, respectively.The upper panel shows the growth and evaporation times averaged over the time until N = 20 peaks in the cluster–cluster case (cf. Fig. 6, bottom right). It is obvious here that the larger sizes develop much faster for the cluster–cluster case. In the bottom panel, the cluster-cluster data is the same, but the monomer–cluster data is extended until the time where N = 20 peaks for that scenario. We note here, like for the cluster size distributions, that the peak of the action is not at the largest nor smallest sizes but around N = 5 … 10.

Fig. 4 TiO2 cluster size distribution in a homogeneous TiO2 bath gas under exoplanet atmosphere conditions at ftot, gas = 103 cm−3. Top: cluster size distributions for T = 1000 K at different times for a constant α = 1.0 for poly–mono and poly–poly. Middle: cluster size distributions for T = 1000 K and various α at t = 4 × 108 s (rand – random choice of the sticking probability). Bottom: comparison of size distributions for different temperatures after ≈ 1.7 × 107 s.

Fig. 5 Average (net) TiO2 cluster growth and evaporation times for α = 1, T = 1000 K, and ftot, gas = 103 cm−3. pp: polymer–polymer; pm: polymer–monomer. Top: the values are evaluated until tmax = 4.12 × 107 s, which corresponds to the JN≥20 peak time for polymer–polymer cluster growth (see Fig. 6, bottom right). Bottom: the values are evaluated until different total run times t to capture the polymer–monomer peaks at t = 0.99 × 109 s (see Fig. 6, bottom left).

3.2 Temporal evolution of formation rates

Figure 6 shows the formation rate (Eq. (6)) of clusters of different sizes in pure TiO₂. The top panels show the early stages of formation rates of clusters larger than size N = 2, 5, and 20 for times up to 1 × 10⁷ s for T = 1000 K and 1500 K and for monomer–cluster scenarios (left) and cluster–cluster (right). The middle panel shows the same but for T = 500 K and 1000 K and extended up to 1 × 10⁹ s. The bottom panel focuses on the T = 1000 K case and includes clusters larger than 20, with times up to 2 × 10⁹ s. From the top panel, we notice that at 1500 K the formation rate of both cluster sizes is initially larger than for 1000 K even though we previously saw that the size distribution was further developed for 1000 K. This shows how the collision dynamics determines the cluster growth. The early part of the temporal evolution is dominated by collisions and an abundance of monomers. As time goes on fewer monomers are available, the evaporations take over, and we notice how the formation rate for N ≥ 2 becomes larger for the 1000 K case. In general there is a much larger overshoot of the formation rates for the 1500 K cases before they find a more stable level.

In the middle panels where the time frame is extended to t = 10⁹ s, we see that the 1000 K cases actually also overshoot but are slower not only to increase but also to relax. For clusters of N ≥ 2, the formation rates for 500 K and 1000 K eventually become nearly identical while for N ≥ 5 the value for 1000 K remains significantly larger than that for 500 K in the monomer–cluster case, while in the cluster–cluster case the values find the same final level.

The lower panels allow for a detailed comparison of the different sticking probabilities at 1000 K. In the left panel (monomer–cluster cases), for N ≥ 2 there is very little differencebetween α = 1 and rand while α = 0.1 clearly lags behind for the early parts of the simulation. The difference between the sticking probabilities becomes clearer as we look at the larger sizes. For N ≥ 20, the largest probability produces the largest production rate for the full duration of the simulation. Looking at the right hand side panel (cluster–cluster collisions), all the production rates for different sticking probabilities become similar much faster. We now understand why the size distributions in Fig. 4 are almost the same for the cluster–cluster cases and only α = 0.1 is smaller than the others for monomer-cluster cases. While the different values of α affect the temporal evolution, the end result remains the same.

A general feature worth noting is that while the formation rates approach a stable value, there are fluctuations around this value. This is a feature of the MC simulations that is not captured in other models.

Fig. 6 Time evolution of the cluster formation rate, JN (cm−3 s−1), for different temperatures and sticking probabilities. Left column: for clusters larger than N = 2, 5, 20 from polymer–monomer cluster growth. Right column: for clusters larger than N = 2, 5, 20 from polymer–polymer cluster growth. Top: early time evolution until t = 107 s for T = 500 K, 1000 K; only J>2 and J>5 appear. Middle: intermediate time evolution until t = 109 s for T = 500 K, 1000 K, for J>2 and J>5. Bottom: long-term evolution until t = 2 × 109 s for J>2, J>5, and also J>20 for T = 1000 K only. No clusters >20 have yet formed for T = 500 K.

3.3 Spatial distribution

Monte Carlo simulations have shown that the cluster formation rates are affected by the stochastic nature of particle motion and particle collisions (Köhn et al. 2020). For α = 1 and T = 1000 K, Fig. 7 shows the spatial distribution of all clusters (top: projected onto the xy-plane; middle: whole 3D distribution) at 2.2 × 10⁹ s. It shows that most (TiO₂)_N clusters remain rather small but that there is one distinct cluster continuously growing (orange). Such a large cluster has a significantly large surface area, which makes it more likely to merge with other clusters in contrast to the less probable collision of two small clusters.

The bottom panel in Fig. 7 shows the mean and maximum radius of all clusters as a function of time. For the simulation allowing polymer–polymer interactions, it shows that the mean radius (orange line) hardly changes in time, limited below 2 nm until t = 2 × 10⁹ s, whilst the blue line depicts the creation of one big particle (orange in the upper two panels). This one distinct cluster forms very early in the nucleation process and subsequently becomes the dominantly growing cluster. The bursts in growth seen in the maximum size is due to scavenging of other large particles by the largest cluster.

The red and violet lines show the mean and maximum radius when only monomers are allowed to attach to clusters. In this case, R_max and R_mean are of comparable size and do not grow significantly until 2 × 10⁹ s, staying smaller than R_mean for polymer–polymer interactions.

Figure 8 shows the same spatial distribution for T = 1000 K and α = 1 as in Fig. 7, but for different realisations of random numbers throughout the whole simulation (e.g. checking for evaporation). It shows that the formation of one big cluster and of many small clusters is independent of the randomness of the collision process. Figure 8 also shows the maximum cluster size (left) and the mean cluster size (right) for three different simulations with different seeds of random numbers. For all simulations, the maximum and mean cluster size show a very similar temporal evolution. The maximum radius becomes ~13 nm after approx. 2 × 10⁹ s (the orange cluster in the panels above) whilst the mean radius for polymer–polymer and polymer–monomer interactions and the maximum radius for polymer–monomer interactions is limited to ≲ 2 nm.

Fig. 7 Spatial distribution of clusters formed by polymer–polymer collisions at t ≈ 2.2 × 109 s for α = 1 and T = 1000 K. The largest cluster (orange) is of size N = 281 755 with R = 13 nm (Eq. (5)). Each cluster radius is amplified by a factor of 107 for visualisation purposes in the plots. Top: projected onto the xy-plane. Middle: 3D distribution of all clusters that formed; the size is colour coded. Bottom: temporal evolution of maximum and mean cluster size.

4 Comparing Monte Carlo and kinetic results for TiO₂ cluster abundances

Figure 9 shows the cluster size distributions based on our kinetic rate equation approach (Eq. (12), Sect. 2.3) for the cluster sizes N = 1 … 10. These results utilise the thermodynamic cluster data as presented in Lee et al. (2015), which are based on the set of (TiO₂)_N cluster structures presented by Jeong et al. (2000). Combining Eqs. (7) and (8) results in (16) (17)

This system of ordinary differential equations (ODEs) allows us to consider the effect of cluster-monomer growth (for i = 1) and cluster-cluster nucleation, similar to Sect. 3. All variables have been defined in Sect. 2.3.

Equations (16) and (17) are solved for a temperature of T = 1000 K. The initial monomer density of f_{TiO_2} = f(1) = 10³ cm⁻³ is chosen in accordance with Fig. 4 in Lee et al. (2015), which is equal to the initial particle density f_tot used in our MC simulations. It is a representative value for the gaseous TiO₂ abundance in the atmospheric regions where T = 1000 …1500 K (and the gas pressure p ≈ 10⁻⁶ × 10⁻⁴ dyn cm⁻² for a model atmosphere of T_eff = 1600 K, log(g [cm s⁻²]) = 3) where TiO₂-seed formation can be typically expected in brown dwarfs and gas giant exoplanets. Higher TiO₂ number densities may occur for higher gas densities. As we focus here on the comparison of two models, the kinetic nucleation and the MC approach, and since MC simulations usually take a significant amount of time (Köhn et al. 2018), here we offer one example. In the future, we plan to extend our comparison to other model atmospheres for example f_tot = 10 or f_tot = 10⁸ cm⁻³.

Fig. 8 Top, middle: spatial distribution of clusters at t = 2.2 × 109 s for α = 1 and T = 1000 K as in Fig. 7, but for two different realisations of random numbers (left and right column). Bottom: maximum and mean cluster size as a function of time for simulations with different seeds of random numbers (solid, dashed, dotted). The black lines in the right hand panel and the red lines in both panels almost lie on top of each other.

Fig. 9 Cluster size distributions f(N) [cm−3 ] as a result of kinetic approach based on thermodynamic cluster data for clusters sizes N = 1 … 10 for α = 1 for T = 1000 K, β = 1.0 for the top two panels, and β = 10−5 for the bottom panel in order to visualise the relative importance of cluster evaporation. The differences between the cluster-monomer (top) and the cluster–cluster (bottom) results appear small for the small clusters, and are amplified if the evaporation efficiency is suppressed.

4.1 Results from thermodynamic cluster formation modelling

Analogously to the 3D MC simulations in Sect. 3, we explore cluster growth by cluster–monomer and by cluster–cluster process for the first 10 (TiO₂)_N clusters. ForT = 1000 K, the solution of Eqs. (16) and (17) is shown in Fig. 9. As previous work (Boulangier et al. 2019) and our MC approach have demonstrated for the considered cases, the growth of TiO₂ clusters is more dominant than the evaporation frequency at 1000 K, which is why we focus on this temperature in the following. Both cases look rather similar, except for the larger cluster sizes that appear earlier if cluster growth proceeds via cluster–cluster collisions. Inorder to understand the importance of cluster evaporation, we have solved the kinetic model for negligible evaporation (β = 10⁻⁵) such that the growth of clusters is the only relevant process. Large and size-dependent differences occur regarding the relaxation timescale to a steady state if the evaporation probability is decreased (i.e. suppressed) uniformly for all clusters (Fig. 9, bottom panel). The artificial suppression of cluster evaporation demonstrates that larger clusters become increasingly more stable against destruction compared to the smaller clusters, driving the cluster size distribution away from an equilibrium distribution. We conclude that the maximum cluster size of N = 10 monomers isstill too small for the divergence of the monomer–cluster and the cluster–cluster growth to appear as pronounced as in our 3D MC simulations.

4.2 Comparison of 3D MC and thermodynamic cluster formation modelling

Our study of (TiO₂)_N cluster formation as a precursor for nucleation seed formation in exoplanet and brown dwarf atmospheres, as well as in AGB star envelopes, points towards considerable timescales before a steady-state solution can be reached. The timescale shortens if all possible clusters participate into the growth process, instead of cluster growth by monomer addition only. This finding is consistently confirmed by both approaches, the molecular-dynamics 3D MC and thermodynamic cluster formation modelling as shown in Figs. 6 and 9. Köhn et al. (2020) have further demonstrated that a decreasing gas density doesnot considerably affect the timescale on which the steady-state formation rate is reached, and that the gas density doesnot affect the value of the steady-state formation rate. The study by Köhn et al. (2020) addresses H₂SO₄–H₂O cluster growth for Earth atmosphere conditions, however.

A comparison between our two methods can only be carried out for the cluster sizes N = 1 … 10 because of the presently available thermodynamic data for the kinetic approach. For T = 1000 K, we compare the cluster–cluster results only. Figure 10 separates the comparison for the small cluster (N = 1 … 5) and the intermediate cluster sizes (N = 6 … 10) because of their qualitative differences. Generally, the cluster number densities for N = 1 … 10 and their time evolution agree reasonably well between the 3D MC and the kinetic nucleation approach. Both reach the steady-state values after a considerable time, hence there is an individually constant number density for all clusters shown. Both approaches predict the small clusters to be the most abundant and the largest cluster sizes to be the least abundant. Figure 10 shows that the density f(N = 10) of 10-mers after one year ~3 × 10⁷ s is on the order 10¹ cm⁻³; using an improved nucleation theory, which abandons equilibrium assumptions, discards growth restrictions, and uses quantum mechanical properties, Boulangier et al. (2019) find f(10) ≈ 10²–10³ cm⁻³ after 1 yr for AGB wind conditions for an initial mass density ϱ_Ti ≈ 2.84 × 10⁻¹⁵ kg m⁻³. As their ϱ_Ti is approximately 35.5 times larger than ϱ_Ti ≈ 8 × 10⁻¹⁷ kg m⁻³, as used in this work, we see that their values agree well with the values obtained by our simulations. We also note that the nucleation rate of clusters does not depend linearly but rather superlinearly on the initial particle density (Dunne et al. 2016; Köhn et al. 2018).

Although the agreement between the MC and the kinetic approach for T = 1000 K may seem a simple observation to make, it is a very important result to emphasise as both methods approach the problem from rather different angles. The characteristic material properties are only described by an evaporation frequency in the 3D MC (Eqs. (2) to (5)) whereas the kinetic nucleation approach requires individual thermodynamic properties for each cluster size (and ideally also their isomers; Eq. (15)). However, the simplicity of the 3D MC approach is offset by the immense computational time demand for a system as small as 10³ particles in 1 cm³.

For the case of α = 1 and β = 1, the individual cluster number densities do not vary by more than a factor of 5 between 3D MC and the kinetic approach in the steady state. As Fig. 4 shows, there is no significant difference in the evolution of the size distributions and thus of the cluster number densities when using α = rand instead of α = 1; for α = 0.1, the temporal evolution of the number densities is delayed by a factor of ~ 10 compared to α = 1. All 3D MC cluster number densities are lower in steady state than the kinetic results, except for N = 1 and N = 2. The 3D MC monomer density (f(N = 1)) fluctuates around a nominal value that results from the condition of keeping the total number of particles within the 3D MC computational domain constant. The divergence of the 3D MC monomer density (solid purple line, top panel in Fig. 10) from the kinetic nucleation values (dashed purple line) occurs when other cluster sizes become increasingly more abundant. That is when the 3D MC code replenishes the monomers such that the total number density is kept constant to the value of 10³ cm⁻³.

All 3D MC cluster number densities in Fig. 10 are lower in steady than the kinetic nucleation results (expect for N = 1, 2) because the3D MC simulations offer reaction paths to larger cluster sizes N > 10 whereas thekinetic model as formulated in Eqs. (16) and (17) is limited by N = 10. Therefore, the growth of larger clusters of N > 10 will deplete the number of smaller clusters of N < 10 in the 3D MC results, but not in the kinetic nucleation results. All 3D MC cluster number densities in Fig. 10 for N > 1 overshoot the kinetic nucleation values during the beginning of the nucleation process at t < 10^7.5 s, meaning they form more efficiently at earlier times than with the kinetic approach. That clusters appear faster with the 3D MC approach than with the kinetic nucleation could be because the randomness of the MC code produces fluctuating collision rates that deviatefrom the average values used by the kinetic nucleation code. Since evaporation is so dominant for the smallest sizes, faster growth due to an above average number of collisions can result in faster growth rates for small clusters, as shown by Olenius et al. (2018) and Köhn et al. (2020).

Another reason could be 3D effects in the MC simulations due to local density inhomogeneities (as seen in Fig. 7), which are not seen in a kinetic modelling approach. However, this finding does suggest that processes that introduce an intermittentdensity distribution of the gas-phase may have the same effect on the cluster formation rate early on in the nucleation process. Such a process is turbulence and its amplifying effect for chemical reactions has been demonstrated in engineering as well as in cloud formation modelling. Helling et al. (2004) simulated turbulent cloud condensation by modelling driven turbulence on mesoscopic scales in brown dwarf and exoplanet atmospheres. It was demonstrated that the efficiency of seed formation increases in a turbulent compared to a laminar gas, and that the formation of CCNs may occur in an event-like fashion. The amplifying effect of a turbulent fluid field on chemical processes can generally be understood as a turbulent fluid field creating an increased reactive surface compared to a laminar flow, an effect that has been demonstrated for thermonuclear deflagration in supernovae (Schmidt et al. 2005). Another effect that can produce this kind of overshoot is a large injection of monomers, such as when a host star rises on a planet. An initial burst of small clusters is then dampened as larger clusters start to scavenge them or the monomer concentration is depleted.

Fig. 10 Comparison of the cluster number densities for N = 1… 5 (top) and N = 6… 10 (bottom) resulting from the 3D Monte Carlo (3D MC, solid) and from the kinetic approach (dotted) for α = 1, β = 1 and T = 1000 K. Deviations are largest for the larger clusters as the kinetic approach is limited to N = 10 as the maximum but the MC approach can grow to N > 10, hence depleting the smaller cluster numbers.

5 Discussion

We have explored the formation of (TiO₂)_N clusters as precursors of cloud formation using a 3D MC approach and a thermodynamic rate equation approach in order to validate both approaches and to test assumptions that are needed on both sides. The resulting cluster number densities are promisingly comparable within a factor of 5. This conclusion is limited to the cluster sizes for which thermodynamic data are available (N = 1…10).

Both approaches demonstrate that the time to approach a steady-state cluster size distribution is substantial for a monomer number density that is reasonable to expect for exoplanet and brown dwarf atmospheres inside their thermodynamic seed formation window. This suggests that the atmospheres of exoplanets and brown dwarfs, as well as the optically thin parts of AGB outflows, may contain mainly monomers in the dynamic phase of cluster formation, but that a more uniform cluster abundance can be a fingerprint for the steady-state situation.

Jeong et al. (2000) provide the wavenumbers for the vibrational transitions for their DFT Gaussian94 (TiO₂)_N cluster simulations (DFT/B3P86/6-31G(d)) for cluster sizes N = 1 … 6, which may serve as a guide for observational approaches to disentangle the state of CCN formation in exoplanets, brown dwarfs, or also for AGB stars similar to Decin et al. (2017). Jeong et al. (2000) point out that larger clusters of size N have higher intensities for IR-active vibrations than smaller clusters. Although lower energy isomers of (TiO₂)_N, N = 3–6, exist (Berardo et al. 2014; Lee et al. 2015), structural (geometric) arguments support our assumption that the cluster growth proceeds via the ‘gaiter’-shaped clusters derived by Jeong et al. (2000). The timescale for a metastable isomer to relax to its global minimum structure can be long, in particular if the related potential energy surfaces (PES) consist of more than one funnel (Doye & Wales 1996). The number of funnels for the relaxation of the (TiO₂)_n, n = 3–6, clusters is unknown, since we did not perform a (computationally very demanding) scan of the PES. However, based on their atomic coordination, a relaxation of the (TiO₂)_n, n = 3–6, clusters used in this study involves a re-arrangement (breaking and formation) of several bonds and is likely to proceed via several funnels.

The line opacity, A_k(N) ⋅ f(N) [10⁻⁸ cm⁻¹ s⁻¹], for the integrated absorption coefficients data A_k [10⁻⁸ cm² s⁻¹], provided by Jeong (2000) for N = 1 … 6, is shown in Fig. 11. The line opacity is evaluated at two different times: t = 10⁷ s (top panel) when the cluster formation invery dynamic and for t = 10⁹ s (bottom panel) when the cluster formation has reached a steady state (see Fig. 10). The line opacities are rather similar for both cases such that an observational differentiation between the two nucleation stages may be difficult. We note that cluster data are re-evaluated all the time; however, what ultimately determines the result is the number of clusters of a given size. The higher abundance of the smaller clusters outweighs the stronger IR line absorption cross sections of the larger clusters. It may however be feasible at some point to observe TiO₂ clusters and Fig. 11 indicatesthat clusters N = 4, 5, 6 would be the best candidates. We therefore reproduce the detailed vibrational absorption spectrum, A_k [10⁻⁸ cm² s⁻¹], for (TiO₂)₄, (TiO₂)₅, and (TiO₂)₆ in Fig. 12 (data from Appendix C in Jeong 2000).

For each of these line spectra, the bottom right panel shows the Lorentz folded line spectra, which is constructed as follows. For given N, we construct the Lorentzian for the wavelength λ_0,i of the ith peak, where we choose a fixed width of w = 1 μm and where C_i is chosen such that the maximum of L_N,i agrees with the maximum of the ithline. Then for each N, the enveloped Lorentz profile is calculated as . As the total line spectrum depends on the particle number for each individual N, we calculate the total spectrum as the sum L₄ + L₅ + L₆ (dashed line).

The (TiO₂)₄, (TiO₂)₅, and (TiO₂)₆ vibrational lines would be well covered by JWST/MIRI (Mid-Infrared Instrument) for λ < 45μm or by ELT (Extremely Large Telescope) instruments like MIDIR (mid-IR imager). However, it may be difficult to disentangle individual cluster lines but instead be more feasible to detect a subset of lines of different clusters for example at the morning terminator of ultra-hot Jupiters. The morning terminator of ultra-hot Jupiters can be expected to be cloud free, hence optically thin to pressures of p_gas ≈ 10⁻² × 10⁻⁴ bar where the atmospheric temperature reaches the 1000K-window. In the future, the cross-correlation technique (e.g. Brogi & Line 2019; Hood et al. 2020; Serindag et al. 2020) may be most suitable to explicitly search for cloud condensation clusters in exoplanet atmospheres if high-resolution opacity data for such clusters are available.

Fig. 11 Line opacity, Ak(N) ⋅ f(N) [10−8 cm−1 s−1], for (TiO2)N clusters N = 1 … 6 in the wavelength window relevant for potential JWST observations. The (TiO2)N integrated line absorption coefficient, Ak [10−8 cm2 s−1], data is from Appendix C in Jeong (2000). The 3D MC cluster number density, f(N) [cm−3 ], is used for t = 107 s (top panel) and t = 109 s (bottom panel). Both times show little difference in opacity despite their different cluster number densities (Fig. 10). (TiO2)4, (TiO2)5, and (TiO2)6 have the largest line opacity amongst the N = 1 … 6 clusters. Their detailed line spectra are shown in Fig. 12.

Fig. 12 Cluster line opacities. (TiO2)N clusters for N = 4, 5, 6 have the strongest IR vibrations line absorption coefficient, Ak [10−8 cm2 s−1], according to Jeong et al. (2000) (their Eq. (3)) amongst their TiO2 cluster ensemble of N = 1…6. The data is reproduced from Appendix C in Jeong (2000). The Lorentz profile folded lines of (TiO2)N N = 4, 5, 6 are shown in the lower right panel for t = 109s.

6 Conclusions

The study of (TiO₂)_N cluster formation presented here has used a 3D MC and a kinetic approach that is based on thermodynamic cluster data (kinetic nucleation). The 3D MC approach enables the following insights into (TiO₂)_N cluster formation as key precursors for the formation of CCNs in exoplanet and brown dwarf atmospheres, and in AGB star envelopes:

• Modelling the formation of (TiO₂)_N clusters by monomer–cluster or by cluster–cluster collisions results in two fundamentally different cluster size distributions. The monomer-cluster scenario has a Gaussian-like shape with a well-defined mean cluster diameter. The cluster size distribution appears power-law like if cluster–cluster growth is enabled;

• Larger clusters form faster by cluster–cluster collisions instead of monomer–cluster collisions only;

• A sweet-spot temperature occurs for most efficient cluster growth. This sweet-spot temperature value (1000 K for TiO₂, Boulangier et al. 2019) also depends on the gas density as both determine the collisional rates of the clusters. Therefore, cloud formation in exoplanet and brown dwarf atmospheres and the dust formation in AGB stars will be triggered within this temperature range.

• The cluster size distributions are different for different temperatures, which is driven by the interplay between growth and evaporation processes;

• We tested the effect of missing material data by studying different sticking probabilities. This material data insufficiency matters for the time evolution of the cluster size distributions but not for the final steady-state results;

• The onset of cloud (or dust) formation may be traced by observing (TiO₂)₄, (TiO₂)₅, and (TiO₂)₆ vibrational lines, which are present at temperatures around 1000 K and may be detectable with instruments such as JWST/MIRI (λ < 45 μm) or ELT/MIDIR. A dedicated search applying cross-correlation would be most desirable but requires a complete-as-possible, high-resolution line list for each cluster.

A comparison of results from the 3D MC method and the kinetic approach based on thermodynamic cluster data led to the following conclusions:

• The 3D MC approach enables the study of the time-resolved formation dynamics of the individual cluster growth, which is not accessible to the kinetic steady-state approach;

• For T = 1000 K, both methods agree well regarding cluster number densities for N = 1 … 10, the vivid onset of cluster formation, and the long transition into a steady state;

• The faster onset of cluster formation in the 3D MC compared to the kinetic results from the increased collisional rates in spots of increased gas density. A turbulent fluid field will cause similar effects;

• Our comparison also supports the 3D MC approach for H₂SO₄–H₂O cluster formation by Köhn et al. (2018, 2020) to model the early stages of CCN formation.

Acknowledgements

Ch.H and M.B.E. acknowledge funding from the European Union H2020-MSCA-ITN-2019 under grant agreement no. 860 470 (CHAMELEON), J.-P.S. acknowledges funding from the St Andrews St Leonard College international scholarship. D.G. acknowledges support from the ERC consolidator grant 646758 AEROSOL. We thank H. Svensmark for discussions on the paper’s topic.

Appendix A Supplementary figures

Fig. A.1 Decision tree of how six initial monomers can grow. Top: Without cluster-cluster growth. Bottom: With cluster-cluster growth allowed. In red we show those paths that do not appear in the polymer-monomer case. The slash notation in the bottom panel gives the individual number of N-mers ordered by their position in the slash line. For polymer-polymer growth, all possible paths end up with one single 6-mer.

References

All Figures

Fig. 1 Top: different cluster radius definitions. Blue line: radius based on the maximal inter-atomic distance from DFT (density functional theory) simulations (Berardo et al. 2014). Orange line: Radius based on the ConvexHull volume approximation of cluster geometry from Berardo et al. (2014) with added Van der Waals interaction distance. Green line: (Eq. (5)). Purple line: radius based on volume of smallest box containing cluster with added Van der Walls interaction distance. Bottom: Cluster number densities from MC code using different radius definitions. Solid lines: (Eq. (5)); dotted lines: box approximation/van-der-Waals radius (purple line in upper panel).

In the text

Fig. 2 Knudsen number, Kn = lMFP∕(2a) (lMFP – mean free path, a – particle radius) for different exoplanet and brown dwarf DRIFT-Phoenix model atmosphere results, except for Jupiter and Venus. The number Kn indicates whether the cluster interactions occur through diffusion (Kn < 1) or in a free molecular flow (Kn > 1). Exoplanet and brown dwarf atmospheres are predominately in the regimes of a free molecular flow (Kn > 1). Knudsen numbers for Earth’s atmosphere and typical cloud particle sizes are ≈ 10−5 and therefore fall in the turbulent, viscous case.

In the text

	Fig. 3 Evaporation frequency γ(N) (Eq. (2)) for TiO2 clusters as a function of their size N for different gas temperatures.
In the text

Fig. 4 TiO2 cluster size distribution in a homogeneous TiO2 bath gas under exoplanet atmosphere conditions at ftot, gas = 103 cm−3. Top: cluster size distributions for T = 1000 K at different times for a constant α = 1.0 for poly–mono and poly–poly. Middle: cluster size distributions for T = 1000 K and various α at t = 4 × 108 s (rand – random choice of the sticking probability). Bottom: comparison of size distributions for different temperatures after ≈ 1.7 × 107 s.

In the text

Fig. 5 Average (net) TiO2 cluster growth and evaporation times for α = 1, T = 1000 K, and ftot, gas = 103 cm−3. pp: polymer–polymer; pm: polymer–monomer. Top: the values are evaluated until tmax = 4.12 × 107 s, which corresponds to the JN≥20 peak time for polymer–polymer cluster growth (see Fig. 6, bottom right). Bottom: the values are evaluated until different total run times t to capture the polymer–monomer peaks at t = 0.99 × 109 s (see Fig. 6, bottom left).

In the text

Fig. 6 Time evolution of the cluster formation rate, JN (cm−3 s−1), for different temperatures and sticking probabilities. Left column: for clusters larger than N = 2, 5, 20 from polymer–monomer cluster growth. Right column: for clusters larger than N = 2, 5, 20 from polymer–polymer cluster growth. Top: early time evolution until t = 107 s for T = 500 K, 1000 K; only J>2 and J>5 appear. Middle: intermediate time evolution until t = 109 s for T = 500 K, 1000 K, for J>2 and J>5. Bottom: long-term evolution until t = 2 × 109 s for J>2, J>5, and also J>20 for T = 1000 K only. No clusters >20 have yet formed for T = 500 K.

In the text

Fig. 7 Spatial distribution of clusters formed by polymer–polymer collisions at t ≈ 2.2 × 109 s for α = 1 and T = 1000 K. The largest cluster (orange) is of size N = 281 755 with R = 13 nm (Eq. (5)). Each cluster radius is amplified by a factor of 107 for visualisation purposes in the plots. Top: projected onto the xy-plane. Middle: 3D distribution of all clusters that formed; the size is colour coded. Bottom: temporal evolution of maximum and mean cluster size.

In the text

Fig. 8 Top, middle: spatial distribution of clusters at t = 2.2 × 109 s for α = 1 and T = 1000 K as in Fig. 7, but for two different realisations of random numbers (left and right column). Bottom: maximum and mean cluster size as a function of time for simulations with different seeds of random numbers (solid, dashed, dotted). The black lines in the right hand panel and the red lines in both panels almost lie on top of each other.

In the text

Fig. 9 Cluster size distributions f(N) [cm−3 ] as a result of kinetic approach based on thermodynamic cluster data for clusters sizes N = 1 … 10 for α = 1 for T = 1000 K, β = 1.0 for the top two panels, and β = 10−5 for the bottom panel in order to visualise the relative importance of cluster evaporation. The differences between the cluster-monomer (top) and the cluster–cluster (bottom) results appear small for the small clusters, and are amplified if the evaporation efficiency is suppressed.

In the text

Fig. 10 Comparison of the cluster number densities for N = 1… 5 (top) and N = 6… 10 (bottom) resulting from the 3D Monte Carlo (3D MC, solid) and from the kinetic approach (dotted) for α = 1, β = 1 and T = 1000 K. Deviations are largest for the larger clusters as the kinetic approach is limited to N = 10 as the maximum but the MC approach can grow to N > 10, hence depleting the smaller cluster numbers.

In the text

Fig. 11 Line opacity, Ak(N) ⋅ f(N) [10−8 cm−1 s−1], for (TiO2)N clusters N = 1 … 6 in the wavelength window relevant for potential JWST observations. The (TiO2)N integrated line absorption coefficient, Ak [10−8 cm2 s−1], data is from Appendix C in Jeong (2000). The 3D MC cluster number density, f(N) [cm−3 ], is used for t = 107 s (top panel) and t = 109 s (bottom panel). Both times show little difference in opacity despite their different cluster number densities (Fig. 10). (TiO2)4, (TiO2)5, and (TiO2)6 have the largest line opacity amongst the N = 1 … 6 clusters. Their detailed line spectra are shown in Fig. 12.

In the text

Fig. 12 Cluster line opacities. (TiO2)N clusters for N = 4, 5, 6 have the strongest IR vibrations line absorption coefficient, Ak [10−8 cm2 s−1], according to Jeong et al. (2000) (their Eq. (3)) amongst their TiO2 cluster ensemble of N = 1…6. The data is reproduced from Appendix C in Jeong (2000). The Lorentz profile folded lines of (TiO2)N N = 4, 5, 6 are shown in the lower right panel for t = 109s.

In the text

Fig. A.1 Decision tree of how six initial monomers can grow. Top: Without cluster-cluster growth. Bottom: With cluster-cluster growth allowed. In red we show those paths that do not appear in the polymer-monomer case. The slash notation in the bottom panel gives the individual number of N-mers ordered by their position in the slash line. For polymer-polymer growth, all possible paths end up with one single 6-mer.

In the text