A 3 D particle Monte Carlo approach to studying nucleation

Article history: Received 12 October 2017 Received in revised form 18 January 2018 Accepted 15 February 2018 Available online 24 February 2018


Introduction
Nucleation of aerosols is the fundamental process by which gas condenses to form stable clusters. These clusters can potentially grow all the way to sizes where they can serve as cloud condensation nuclei (CCN), typically 50-100 nm. The nucleation phenomenon has been observed all around the globe [23] and is also considered to contribute to the formation of clouds on brown dwarfs and exoplanets [12,25]. About half of all CCN are estimated to originate from nucleated aerosols [29], making nucleation a relevant topic not only for aerosol research but also for its implications for cloud formation. Additionally both aerosols and clouds are relevant for e.g. climate change due to their large forcing effects [3]. The key molecule for aerosol nucleation has long been thought to be sulphuric acid (with water or other stabilizing molecules) due to its ability to form strong bonds [5]. Recently it has been shown that highly oxygenated organic molecules are also able to nucleate at high altitudes [2]. Traditionally nucleation has been described by classical thermodynamic nucleation theory [11], kinetic numerical models [32,27,39], or parametrisations based on either nucleation theory (e.g., [37]) or experimental data (e.g., [7]). The parametrisations and numerical models have the advantage that they can be adapted for use in global modelling [35,41,31] due to being computationally quick. Kazil and Lovejoy [17] used a semi-analytical approach to add aerosols to a global model.
More recently an Atmospheric Cluster Dynamics Code (ACDC) has been developed solving the so-called birth-death equations, ordinary differential equations (ODEs) describing the temporal evolution of cluster densities of a given size [28]. The novelty of such a model is the automatic generation of ODEs for a given cluster size and its implementation into the solver whenever needed. This approach aims to reduce typographical errors when implementing ODEs manually. Such traditional numerical models can provide information with regards to particle size distributions and may reflect the physics satisfactory using actual data for condensation and evaporation. The information provided by these models is, however, focused on developments in time and not in space. If each molecule could be tracked in a three dimensional space it might be possible to achieve new insights into the process of nucleation.
Tracing individual particles in space and time is the main advantage of particle Monte Carlo codes. They are for example widely used to simulate the properties of lightning discharges [4,19,20] by tracing individual electrons and photons or to study the nanostructure growth of atoms on surfaces in electrochemical models [14,15,9]. Recently, molecular dynamic models have also been used to study vapour-to-liquid nucleation [6] or homogeneous water nucleation [1].
Monte Carlo models give the opportunity to include all relevant microphysical processes as well as the interaction amongst all involved particles. In contrast to kinetic models or pure parametrisations which are based on averaged quantities, such as the density or concentration of particles or the mean energy, these models are able to capture rare events initiated by single particles [33,13] as for example the production of gamma rays or positrons in the vicinity of lightning discharges [18].
The disadvantage of Monte Carlo codes is their runtime. Depending on the size of the problem Monte Carlo simulations can take up to several weeks whereas models based on averaged quantities take several hours to days [26]. For small systems, however, the time difference is not significant and Monte Carlo models offer a much better approach to the discrete nature of particles. We here present a particle Monte Carlo code to study the nucleation of sulphuric acid clusters in 3D, which to our knowledge has not been done before. We here emphasize, however, that we not only present a new model to simulate the nucleation of sulphuric acid clusters, but also that the implemented physics is sufficient enough to study the nucleation of neutral sulphuric acid clusters in order to benchmark our model. Previously Monte Carlo studies have been used on the sulphuric-acid water system to study cluster parameters such as cluster shape, conformation and dissociation [22] as well as cluster free energies [16]. The paper is organized as follows. In section 2 we introduce the Monte Carlo particle model for the study of nucleation processes and discuss all the ingredients of this model: the implementation of single particles, diffusion coefficients, the collision of particles, the evaporation coefficients as well as the choice of the time step. In section 3 we discuss the spatial and size distribution of single sulphuric acid clusters and compare nucleation rates calculated with the present model with values from literature. We finally conclude in section 4.

Modelling
For the study of the growth of H 2 SO 4 -H 2 O clusters, we introduce a particle Monte Carlo code tracing individual H 2 SO 4 -H 2 O clusters. In order to benchmark this model, we perform simulations with initial sulphuric acid molecule densities of n = 10 7 cm −3 and of 10 8 cm −3 at 200 K and 300 K as well as of n = 2.5 · 10 7 cm −3 at 238 K with a relative humidity of approximately 50%.
In the simulations we do not distinguish between molecules and clusters, we simply refer to them as particles. Each particle is described as a sphere characterized by its position r = (x, y, z) in Cartesian coordinates. The particle's radius R consisting of k sulphuric acid molecules is determined through [38] implicitly including water content for a relative humidity of approx. 50%. The mass and the density are given through [37] with the Avogadro constant N A ≈ 6.022 · 10 23 mol −1 and the mole fraction which is a powerlaw fit to Fig. 2b of [38] for a temperature of approx. 300 K for a humidity of 50%. According to the same figure, we multiply (4) with a factor of ≈ 0.8 for temperatures of 200 K and 238 K.
After every time step t, the position is updated through 2 depends on the particle size R and on the diffusion coefficient [8] for molecules where R(1) = 0.329 nm [21] and m(1) = 2.0033 · 10 −25 kg are the initial size and the mass of single sulphuric acid molecules without any attached water molecules. k B ≈ 1.38 · 10 −23 J/K is the Boltzmann constant and P = 1 bar the ambient pressure. For T = 300 K, T = 238 K and T = 200 K the diffusion coefficients are D 0 ≈ 1.65 · 10 −6 m 2 /s, D 0 ≈ 1.16 · 10 −6 m 2 /s and D 0 ≈ 8.96 · 10 −7 m 2 /s.
The time step used is related to the time it takes a particle to diffuse the average separation between particles. The average separation between particles is 3 √ 1/n and the related diffusion length is Hence it is t ∼ 0.653 ms for n = 10 7 cm −3 and 300 K, t ∼ 0.504 ms for n = 2.5 · 10 7 cm −3 and 238 K, t ∼ 1.202 ms for n = 10 7 cm −3 and 200 K as well as t ∼ 0.141 ms for n = 10 8 cm −3 and 300 K. Note that the actual time step should be t and that t limits the maximal distance which a particle travels. Hence, we have chosen t = 100 μs for n = 10 7 cm −3 and 300 K, t = 250 μs for n = 2.5 · 10 7 cm 7 and 238 K and t = 1 ms for n = 10 7 cm 7 and 200 K. However, since the diffusion term in (5) depends on a Gaussian random number and since for polymers it is D(R) < D 0 , the actual travelled distance is smaller than the maximally possible value. With t = 0.1 ms or t = 1 ms, it is √ 2D 0 t ≈ 18 μm for 300 K and √ 2D 0 t ≈ 42 μm for 200 K. Since the mean free path of sulphuric acid molecules in air is approximately 10-60 nm [34], hence significantly smaller than 1 μm, the justification of the diffusion approach in this study is justified.
After every time step we check whether two particles i and j with sizes R i and R j overlap by evaluating the condition where r i are r j are the particles' positions. If this condition is fulfilled, the two particles are merged with new mass and radius according to (1) and (2); the new position is determined through Note that we add a new H 2 SO 4 -H 2 O cluster to the simulation domain at a random position after merging two particles if the particle number becomes smaller than the initial particle number; thus, the density of particles does not drop below the initial particle density. Vice versa we also check after every time step whether a cluster with radius R and mass m evaporates by emitting one sulphuric acid molecule added to the simulation domain leaving a cluster with reduced mass and volume behind. For that purpose, we check after every time step and for every cluster whether r ≤ P eva (10) for a random number r ∈ [0, 1) and for the evaporation probability P eva = 1 − exp(−γ t). Here the evaporation frequency γ for 300 K is adopted from [38] where the concentration n ∞ a,sol of H 2 SO 4 vapour molecules in the equilibrium vapour above a flat source and the surface tension are determined through [30,37]    Because of the lack of physical data, e.g. for the equilibrium vapour pressure, we multiply Eq. (11) with a factor of ≈ 10 −4 for 238 K and of ≈ 10 −5 for 200 K based on Fig. 4 of [38]. For better accuracy, we use Eq. (11) for k > 5 only; for k ≤ 5, we use distinct values extracted from [38]. Fig. 1 shows the evaporation coefficient γ as a function of k for 200 K, 238 K and 300 K. For all considered cases evaporation is most efficient for small cluster sizes and becomes less probable with increasing size. For 200 K and 238 K γ is several orders of magnitude smaller than for 300 K, thus compared to 300 K evaporation is much less likely and nearly negligible for 200 K. Since the emitted molecule is added to the simulation domain, the particle density slightly increases; however, we have observed that it does not exceed 1.002 · 10 7 cm −3 for an initial density of 10 7 cm −3 and 1.009 · 10 8 cm −3 for an initial density of 10 8 cm −3 .
In order to optimize the runtime of our simulations, we have chosen the volume of the simulation domain to be 10 −4 cm 3 for a density of 10 7 cm −3 , 4 · 10 −5 cm 3 for a density of 2.5 · 10 7 cm −3 , and 10 −5 cm 3 for a density of 10 8 cm −3 .
Consequently we initiate all simulations by placing 1000 individual H 2 SO 4 molecules at random positions into the simulation domain which scales with the concentration.

Results
In the following we discuss the spatial distribution and the size distribution of all particles as well as the nucleation rates of clusters with a radius of 0.85 nm. As supplementary material we have added movies showing the temporal evolution of the particle position and size for n = 10 7 cm −3 and 200 K as well as for n = 10 7 cm −3 and 300 K.    Fig. 3 shows the spatial distribution of all particles projected onto the xy plane after 50 s. For better visibility we have multiplied the real size of each particle by a factor of 10 6 . Fig. 4 illustrates the capability of the present model for full three dimensional simulations. It shows the spatial distribution after the same time steps and for the same densities and temperatures as in Fig. 3 a) and c). Fig. 3 and 4 demonstrate that in all cases the particles are distributed randomly within the simulation domain. Panels a) and b) of Fig. 3 show that for T = 300 K, there is no growth at all; the average radius of all particles after 50 s is approximately 0.329 nm. In contrast to a temperature of 300 K where the simulation is evaporation  Fig. 3 c) and d). The right y-axis in the right column shows the difference between the particle numbers after 80 s and after 50 s. dominated, there is a significant growth of clusters for a temperature of 200 K and 238 K where evaporation is less probable to negligible. The average size is 0.48 nm for 200 K and 0.46 nm for 238 K. After 80 s, the size distributions at 300 K have not changed at all. Panel a, however, shows that for a temperature of 200 K, the number of monomers has decreased enormously since there is only nucleation, but no evaporation. The right y-axis of panels b) and d) show that there is an increase of polymers with sizes above 0.8 nm; this effect is more significant for 200 K than for 238 K since evaporation is more likely for the latter case.

Nucleation rate
For 300 K, the simulation is evaporation dominated. Hence we do not observe any polymers and conclude that the nucleation rate is 0.0. In the case of 238 K and 200 K, there is an interplay between the nucleation through the approach of particles and the evaporation. As Fig. 3 and 5 show, there is a significant growth of the particle sizes. After every time step t we therefore count the number N(R, t) of particles with size R. Subsequently the nucleation rate of particles with size R as a function of time t is given through where V is the volume of the simulation domain.  Table 1 The nucleation rates ν sim calculated from our simulations, the nucleation rates ν par [7] and the nucleation rates by [43]. As an example, Fig. 6 shows the nucleation rate ν(R = 0.85 nm, t) of particles with a radius of 0.85 nm as a function of time for n = 10 7 cm −3 and T = 200 K. For comparison with literature, e.g. with [7], we have chosen a radius of 0.85 nm; however, this critical size is rather arbitrary. One of the advantages of the introduced model here is the possibility to determine the nucleation rate of particles of any size. For 200 K, evaporation is negligible and thus particles keep growing and consequently the nucleation rate increases as a function of time. In the beginning, the nucleation rate grows slowly since it takes time for the first clusters to reach the nucleation size. However, the slope of the nucleation rate flattens because of the decreased diffusion coefficient for larger clusters which diffuse more slowly than small clusters and single molecules and consequently collide less with surrounding particles.
We here define the mean nucleation rate of particles with a radius of 0.85 nm and above within time interval T as Dunne et al. [7] simulated the formation of atmospheric aerosol particles in extensive laboratory experiments. They determined the nucleation rates for aerosols at different temperatures and for different compounds and derived the parametrisation 1000.0 (16) for the neutral nucleation rate of sulphuric acid-water clusters with a radius of 0.85 nm as a function of the density n in units of 10 6 molecules/cm 3 and of the temperature T in Kelvin. Additionally, Yu [43] set up a website which calculates nucleation rates ν yu based on [40] and [42]. Table 1 compares the nucleation rates ν sim calculated from our simulations with the nucleation rates ν par ν yu . In all cases ν sim is comparable to the values obtained by Dunne et al. [7] or [43], hence there is a good agreement within the error bars. This might be due to underestimating the evaporation coefficient. However, from the available data, it is difficult to get the exact values for the evaporation coefficients. However, we here emphasize that this paper is not about finding the best input data, but to introduce a code that -provided sufficiently accurate input data -models nucleation on a particle level.

Conclusions
We have presented a particle Monte Carlo model tracing individual H 2 SO 4 -H 2 O clusters for different initial densities, n, of sulphuric acid molecules and for different temperatures T taking the growth by collision of particles and evaporation by single H 2 SO 4 molecules into account.
Four cases were considered, n = 10 7 cm −3 with temperature either T = 300 K or T = 200 K, one case with n = 2.5 · 10 7 cm −3 and T = 238 K and one case with n = 10 8 cm −3 and T = 300 K. In the two cases at 300 K, there is no significant growth between 50 s and 80 s, in accordance with a low nucleation rate. However, for 200 K and 238 K, growth of clusters between 50 s and 80 s is observed.
The simulations make it possible to calculate the size distribution based on individual particle sizes. For 300 K, all particles are single molecules. This behaviour is rather independent of the initial density. In contrast at 200 K and 238 K, evaporation becomes improbable or negligible and thus the size distribution consists of fewer monomers compared to the case with T = 300 K which increases the number of clusters above 0.329 nm.
Finally we have calculated the nucleation rates as a function of time and the mean nucleation rates averaged over time. Since for 300 K there is an interplay between nucleation and evaporation, the nucleation rate is 0. For lower temperatures, evaporation is less likely and as such the nucleation rate increases with time. However, it tends to an upper limit since the probability of two particles colliding with each other is reduced as an effect of the decreased diffusion coefficient.
We compared the nucleation rates with values experimentally obtained by Dunne et al. [7] and F. Yu [43]. Within the given error bars we see a good agreement between our simulation results and experimental values which serves as a benchmark for our Monte Carlo code. We therefore conclude that the physics implemented in the present Monte Carlo model is appropriate to simulate the growth of sulphuric acid clusters. Its main advantage is that it traces individual particles and therefore reflects their distinct nature and as such reality much better than for example kinetic numerical models or pure parametrisations. Subsequently, the 3D particle approach allows us to study new phenomena. In future work, we intend to investigate how the temporal evolution of the nucleation rate changes for fixed particle densities and different volumes, thus different particle numbers. Such a task is ideal for a model tracing individual particles and cannot be studied with kinetic numerical models or pure parametrisations.
The model presented in this paper can easily be adjusted to study nucleation phenomena of other species in different temperatures and pressures. Lee et al. [25] studied the formation of TiO 2 and SiO clusters on brown dwarfs and extra-solar planets. They use a classical nucleation theory where they allow only monomers to attach to other monomers or clusters. In future work, we will adjust our model to include cluster-cluster collisions and to study their effects on the nucleation of clusters in the atmospheres of exoplanets.
In a forthcoming paper, we will finally present a more sophisticated model where we include HSO − 4 ions and investigate their influence on the nucleation rate. Not only is this of relevance for aerosol formation [10] and growth [36] in Earth's atmosphere, but this model can be adjusted to study the influence of ions in exoplanets' atmospheres which has not been considered before.