Robust metric for quantifying the importance of stochastic effects on nanoparticle growth

Comprehensive representation of nanoparticle dynamics is necessary for understanding nucleation and growth phenomena. This is critical in atmospheric physics, as airborne particles formed from vapors have significant but highly uncertain effects on climate. While the vapor–particle mass exchange driving particle growth can be described by a macroscopic, continuous substance for large enough particles, the growth dynamics of the smallest nanoparticles involve stochastic fluctuations in particle size due to discrete molecular collision and decay processes. To date, there have been no generalizable methods for quantifying the particle size regime where the discrete effects become negligible and condensation models can be applied. By discrete simulations of sub-10 nm particle populations, we demonstrate the importance of stochastic effects in the nanometer size range. We derive a novel, theory-based, simple and robust metric for identifying the exact sizes where these effects cannot be omitted for arbitrary molecular systems. The presented metric, based on examining the second- and first-order derivatives of the particle size distribution function, is directly applicable to experimental size distribution data. This tool enables quantifying the onset of condensational growth without prior information on the properties of the vapors and particles, thus allowing robust experimental resolving of nanoparticle formation physics.


Simulated molecules and particles
The simulated molecules are representative of atmospheric low-volatile and extremely-low-volatile organic compounds (LVOC and ELVOC, respectively). The properties of these species, listed in Supplementary Table S 1, are based on ref. 1, and result in apparent growth rates GRapp similar to those observed. We simulate both a one-compound LVOC system and a two-compound LVOC-ELVOC system. The simulated particle size range includes all discrete particle compositions from single vapor molecules to particles of a mobility diameter of dp = 10 nm in case of LVOC, and dp = 4.5 nm in case of LVOC-ELVOC mixture. The size range is smaller for the two-compound system, as the number of all possible individual molecular compositions in a given size range, and consequently the number of coupled differential equations in the discrete GDE (Eq. (1)), increases drastically for systems of multiple chemical compounds. For LVOC-ELVOC, the size range up to 4.5 nm contains already 3730 individual particle compositions, and including larger sizes rapidly becomes computationally unaffordable. For LVOC, the size range up to 10 nm contains 1350 discrete particle sizes. In order to be able to study wider size ranges, we focus on the LVOC system, but report results also for the LVOC-ELVOC mixture to demonstrate that the presented analysis tools are valid also for multi-compound systems. The saturation vapor pressure of 10 −10 Pa used for ELVOC is somewhat higher than the lowest reported estimates of ca. 10 −13 -10 −12 Pa 1 . The higher value was used as (1) it has been proposed that these types of compounds may not have as low saturation vapor pressures as often assumed 2 , and (2) using lower values in the simulations results in very high particle concentrations and strong particle self-coagulation, making condensational growth less important in the sub-5 nm size range.
The evaporation rate of compound k from particle of size dp is in most simulations approximated with the Kelvin formula where βk is the rate constant of the corresponding molecular collision, and psat,k and χk are the pure liquid saturation vapor pressure and the particle-phase mole fraction of k, respectively. σ is the particle surface tension, kB is the Boltzmann constant, and T is the temperature. In addition, we perform simulations using qualitatively and quantitatively different evaporation rates to show that the obtained conclusions are independent of these rates.
The default simulation set-up corresponds to a laboratory chamber experiment (Supplementary Table S 2). The initial concentrations of vapors and particles are set to zero, a constant vapor source is turned on, and the simulation is run until the concentrations of particles in the simulated size range have reached a steady state. The walls of the chamber and the dilution of the chamber air act as sinks for the vapors and particles (Supplementary Table S  Sink corresponding to the walls of the laboratory chamber CLOUD S (dp) = 10 −3 nm s −1 / dp + 9.6•10 −5 s −1 , where dp is the mobility diameter of a particle or vapor molecule, and the size-independent term corresponds to dilution of the chamber air 17 Sink corresponding to the surface of larger aerosol particles in the ambient atmosphere S (dp) = Sref • (dp / dp,ref where Sref and dp,ref are a reference sink and size corresponding to LVOC vapor 15 . Sref was set to either 10 −3 or 5•10 −4 s −1 .

Additional simulation set-ups
To examine the effect of external conditions on the apparent growth of the particle population, additional simulations were performed using a different size dependence and magnitude for the sink rate, as well as using no sinks and a time-independent vapor concentration. The alternative sink rate (Supplementary Table S 2) corresponds to scavenging by a background population of larger aerosol particles. In addition, test simulations were performed using the default chamber set-up but different properties for the vapor and particles. The mass, density, surface tension and saturation vapor pressure were set to m = 98.08 amu, ρ = 1830 kg m −3 , σ = 5•10 −2 N m −1 , and psat = 10 −9 Pa, respectively, corresponding to quasi-unary, stabilized sulfuric acid 3 .
Finally, test simulations were conducted using different particle evaporation rate profiles. While accurate assessments of evaporation rates of small clusters remain unavailable, the rates are expected to exhibit discrete changes as a function of cluster size and composition instead of a smooth behavior such as given by Eq. (S1) 4,5 . This affects the details of the particle distribution, and thus we performed test simulations using different types of non-smooth evaporation profiles. These rate profiles correspond to the qualitative behavior and order of magnitude of evaporation rates calculated with different quantum chemical methods [4][5][6][7][8] . The rate profiles were obtained as follows: Rates calculated by the Kelvin approach using the LVOC properties were taken as a starting point, and random noise was introduced to the rates, with the magnitude of the noise decreasing as a function of particle size as In Eq. (S2), γmod and γKelvin are the modified and the Kelvin-formula-based evaporation rate constants, respectively, ΔOoM is the maximum change in the order of magnitude of the rate, and R{0…1} is a randomly generated number between 0 and 1. The purpose of the exponential factor is to decrease the absolute modification as the number of molecules N in the particle increases, making γmod approach the macroscopic rate γKelvin. Two sets of modified rates were generated using maximum changes of ΔOoM = 2 and ΔOoM = 3. These tests do not address the chemical accuracy of the properties assigned to the compounds; instead, the aim is to use evaporation rates of a realistic order, and study the effect of the shape of the rate profile on the applicability of the presented analysis tools.
In addition, test simulations were performed for representative sulfuric acid-ammonia (H2SO4-NH3) and sulfuric acid-dimethylamine (H2SO4-DMA) systems using available quantum chemical data for the smallest clusters (Supplementary Table S  (1) The rates were computed from the standard Kelvin formula (Eq. (S1)) with psat set to result in realistic vapor and particle concentrations.
(2) The rates were set to include a prefactor Γ, corresponding to an activity coefficient for non-ideal solutions, in order to mimic the strong acid-base interaction originating from hydrogen bonding and salt formation in the particle phase.  Evaporation rates γ for clusters consisting of up to 5 H2SO4 and 5 NH3 molecules from QC 6,18 , and for larger particles from Eq. (S1) with psat,H2SO4 = 10 −9 Pa, psat,NH3 = 10 −8 Pa, ρH2SO4 = 1830 kg m −3 , ρNH3 = 696 kg m −3 , and σ = 5•10 −2 N m −1 2 Otherwise same as H2SO4-NH3 case 1, but with psat,H2SO4 = 10 −8 Pa, psat,NH3 = 10 −7 Pa, and γ from Eq. (S1) modified to include an activity coefficient where χk is the mole fraction of the evaporating compound 3 Quasi-unary H2SO4-NH3 system with the molecular properties of H2SO4, and γ for the smallest clusters from QC along the main cluster growth pathway 18  Otherwise same as H2SO4-DMA case 1, but with psat,H2SO4 = 10 −8 Pa, psat,DMA = 10 −8 Pa, and γ from Eq. (S1) modified to include an activity coefficient where χk is the mole fraction of the evaporating compound

Numerical solution of the discrete GDE
The discrete GDE (Eq. (1)) was solved by a standard Euler method using an adaptive time step. The step size is adjusted as follows: At every integration step, the initial estimate for the step size Δt is determined based on the previous integration step and the maximum allowed relative change δCmax = 0.01 in the concentration of an individual particle i. The former, referred to as Δtnew, is given by Eq. (S4) below, and the latter, referred to as Δtmax, is obtained from the initial concentrations Ci (t0) and their derivatives dCi / dt as The time step Δt is set to be the smaller one of Δtnew and Δtmax, and particle concentrations after Δt are evaluated using a single integration step. The time interval Δt is then integrated again using two steps of length Δt / 2, re-evaluating the derivatives dCi / dt at the midpoint. The two solutions for each i are compared, and if their relative difference ε does not exceed a set tolerance of εmax = 10 −5 , Δt is accepted.
Otherwise, Δt is decreased as 9 and a new attempt is made until the criterion regarding εmax is satisfied. After successful integration, an estimate Δtnew for the next time step is determined according to Eq. (S4), and the procedure continues as described above. The criteria εmax and δCmax apply to all particles that have a concentration of at least 10 −12 cm −3 . The Eulerian integration was tested against the Fortran ODE solver VODE 10 using small simulation systems, ensuring that the method is reliable.
When a collision leads to a particle outside of the simulated size range, the particle is placed in an additional size bin, which represents all larger particles grown out of the simulation range. Particles in the additional bin have the same size corresponding to the size of particles at the boundary of the simulation range. These outgrown particles, the purpose of which is to account for the additional coagulational loss caused by larger particles, can coagulate with the simulated particles and each other and be lost to the sink, but they do not evaporate or grow in size.

Evaluation of the methods for determining
In Eq. (S5), the terms of the continuous condensation equation (Eq. (4)) are marked with "cont.", and the rest of the terms must be negligible for the continuous assumption to hold. The first terms on the righthand side of Eq. (S5) demonstrate that ∂ 2 c/∂i 2 needs to be minor compared to ∂c/∂i in order to describe vapor condensation and evaporation neglecting stochastic widening, also when other dynamic processes affect the size distribution.

Two-compound mixtures
The mixture of LVOC and ELVOC compounds was modeled assuming conditions similar to the LVOC simulations, but splitting the vapor source in a 90:10 ratio between LVOC and ELVOC according to the relative abundances of these types of compounds assessed by laboratory measurements 1 . A test simulation was performed using a 50:50 split. The behavior of GRapp and GRcond is similar to the LVOC system (Supplementary Figure S 5a), but the absolute values are increased by the lower-volatility compound. As for the LVOC simulations, the metric ∂ 2 :∂ is able to predict the size range where GRapp and GRcond converge (Supplementary Figure S 5b; Figure 2b). Also all the H2SO4-base systems simulated applying quantum chemical data (Supplementary Table S 6)) for a set of different model cases. Obviously, the stochastic processes, which lead to increased GRapp, are quantified differently by TREND: at the beginning of the formation event, GRTREND is high at the small sizes similar to GRapp, but the values decrease close to zero at later moments in time (Supplementary Figure S 7). This is because TREND gives an average growth rate based on the shift of specific subsections within the size distribution. At the beginning, the initial sizes build up, forming a "moving front", and TREND sees a positive shift towards larger sizes. After the initial build-up of the distribution, there exist larger particles that evaporate back towards the smallest sizes, and the forward and backward particle fluxes together result in a minor shift in TREND. In any case, it is important to note that GRTREND is determined under the assumption that particles of a given size and at a given point in time grow at the same rate. As a result a particle that is larger than another particle is assumed to remain largerno matter if there is growth or shrinkage. This is generally not the case in the size range governed by stochastic growth, and similarly to GRapp, GRTREND cannot be interpreted through condensation at these sizes (see also Figure 2a). However, the benefit of TREND is that it excludes the effects of particle sinks and coagulation, and gives also the time dependence of the growth rates. Supplementary Figure S 7 shows the size-and time-dependent GRTREND together with GRcond,vapor for representative LVOC simulations in which the growth is either driven by vapor and small clusters (panel (a)), or particle coagulation has a significant role (panel (b)). To compare the TREND results to those obtained using GRapp, Figure 3 presents GRTREND at the particle appearance times for the LVOC and also for the LVOC-ELVOC system at different vapor source levels. For high vapor sources which generate strong particle formation and coagulation, GRTREND and GRcond,vapor show stronger size-and time-dependence due to varying vapor concentration levels (Supplementary Figure S 7b), caused by the large number of formed particles becoming a significant sink for vapor. For these cases, GDE-based methods such as TREND are needed to extract the condensational growth rate, as GRapp cannot be interpreted solely through condensation even for the larger sizes ( Figure 3).

Effect of the definition of tapp on GRapp vs. GRcond
In this work, the appearance time tapp of a given particle size is by default defined according to the 50 % criterion, and GRcond is accordingly calculated and compared to GRapp at this moment. This definition is however a matter of choice, and different criteria, including a 5 % criterion, have also been used 11 . Supplementary Figure S 8a shows the size-and time-dependent region where the continuous condensation model is expected to be valid, based on ∂ 2 :∂ decreasing below five percent, together with tapp defined according to different criteria. The graph demonstrates that the choice of tapp affects the uncertainties related to the apparent growth rates: DGR is larger for GRapp defined according to the 5 % criterion (Supplementary Figure S 8b). At the beginning of the formation event, the size-dependent particle concentration is steeper as the distribution at larger sizes has not yet built up, and the higher gradient affects the diffusion term in Eq. (3). Based on this, it is recommended that tapp is defined such that GRapp is not determined until the distribution around the size of interest has taken a less steep shape, for example according to the 50 % criterion instead of the 5 % criterion.

Remarks on the interpretation of the apparent growth rate GRapp
GRapp determined from the simulation data exhibit a trend similar to previously reported experimental observations 1,11-13 : GRapp is low at a couple of nanometers, and shows a rather steep increase to higher values as the particle size increases. This type of behavior is often interpreted as a thermodynamic barrier hindering the growth at the smallest sizes. While evaporation definitely has a decreasing effect on particle growth rate, the time evolution of the size distribution at the smallest sizes, and consequently the apparent growth are not determined solely by deterministic condensation and evaporation processes. In addition to stochastics, due to their high mobility the smallest sizes are more strongly affected by sinks and coagulation processes compared to larger particles. Moreover, in the presence of a vapor source, the first particle sizes may appear when the vapor concentration is still increasing, further affecting the apparent size-dependent growth rate.
Supplementary Figure  demonstrates that the behavior of GRapp at the very first sizes is largely determined by particle sinks and sources: a test simulation performed using a constant vapor concentration and omitting the external sink gives a different behavior for GRapp at the smallest sizes compared to the default laboratory chamber setup including a vapor source and a particle sink onto the chamber walls (see also ref. 3). On the other hand, a simulation performed at chamber conditions with no Kelvin factor (the exponential factor in Eq. (S1)) results in a similar, increasing GRapp at the small sizes as simulations where the Kelvin factor is included (Supplementary Figure S 9b). This indicates that the existence of thermodynamic barriers in the growth of very small nanoparticles cannot be directly resolved based on the size dependence of the apparent growth rate. In general, vapor and particle sinks affect the apparent growth in different ways. First, sinks may affect ∂ 2 :∂ and DGR by increasing the gradient of the particle size distribution function ( Figure 2b). Second, sinks decrease particle life time, which shortens the time scale of steady-state equilibration and may lead to faster particle appearance and thus higher apparent growth 14 . Third, sinks affect the condensation and evaporation fluxes: scavenging decreases the concentrations of vapors and small clusters, decreasing the growth rate of a given particle. On the other hand, scavenging may also have an increasing effect on GRapp due to reduced concentrations of particles larger than the given size, as this reduction decreases the backward evaporation fluxes to the given particle 3 .

Nanoparticle growth rates in atmospheric aerosol formation modeling
In regional and global climate and air quality models, the initial steps of aerosol particle formation from atmospheric vapors are commonly approximated by (1) determining the initial formation rate J (dp,1) of nanoparticles of ca. dp,1 = 1-1.5 nm based on nucleation theories, and (2) scaling the initial formation rate to a larger size, most often ca. dp,2 = 3 nm, based on the estimated growth rate GR and particle scavenging sink S at sizes between dp,1 and dp,2 using the continuous GDE as 15,16 ( p,2 ) = ( p,1 )exp (− ∫ GR p,2 p,1 d p ). (S6) The growth rate GR is typically assumed to be size-independent, and the loss rate S is approximated to decrease with increasing particle size as S (dp) = Sref • (dp / dp,ref) −m , where Sref is the loss rate at reference size dp,ref, and m is a parameter normally set to −1.6 (Supplementary Table S While the continuous condensation approach is a source of uncertainty in Eqs. (S6) and (S7), the calculated J may be more severely distorted by the assumed growth rate GR. An overestimation of GR results in an overestimation of the formation rate J (dp,2), with the quantitative error depending on the sink Sref as well as the other parameters in the exponent of Eq. (S7). Assuming GR = 1 nm h −1 , using a typical boundary layer sink of Sref = 10 −3 s −1 , and scaling J from 1 to 3 nm, an overprediction of a factor of 2 or 5 in GR leads to an increase of a factor of 2.3 or 3.8 in J , respectively. For 5•10 −3 s −1 , the corresponding values are 66 and 810. Generally, assuming GR between 1 and 5 nm h −1 and Sref between 10 −4 and 10 −2 s −1 leads to an increase in J of a factor ranging from slightly more than one to over a hundred. Conditions where this error becomes particularly significant include relatively unpolluted environments with slow or average nanoparticle growth but non-negligible scavenging sink, as demonstrated by Supplementary Figure S 10.
Supplementary Figure S 10: Factor by which the particle formation rate scaled from dp,1 = 1 nm to dp,2 = 3 nm by Eq. (S7) increases when the growth rate GR is overestimated by a factor of two for different values of GR and Sref.