Atmospheric Cluster Dynamics Code: a flexible method for solution of the birth-death equations

Abstract. The Atmospheric Cluster Dynamics Code (ACDC) is presented and explored. This program was created to study the first steps of atmospheric new particle formation by examining the formation of molecular clusters from atmospherically relevant molecules. The program models the cluster kinetics by explicit solution of the birth–death equations, using an efficient computer script for their generation and the MATLAB ode15s routine for their solution. Through the use of evaporation rate coefficients derived from formation free energies calculated by quantum chemical methods for clusters containing dimethylamine or ammonia and sulphuric acid, we have explored the effect of changing various parameters at atmospherically relevant monomer concentrations. We have included in our model clusters with 0–4 base molecules and 0–4 sulfuric acid molecules for which we have commensurable quantum chemical data. The tests demonstrate that large effects can be seen for even small changes in different parameters, due to the non-linearity of the system. In particular, changing the temperature had a significant impact on the steady-state concentrations of all clusters, while the boundary effects (allowing clusters to grow to sizes beyond the largest cluster that the code keeps track of, or forbidding such processes), coagulation sink terms, non-monomer collisions, sticking probabilities and monomer concentrations did not show as large effects under the conditions studied. Removal of coagulation sink terms prevented the system from reaching the steady state when all the initial cluster concentrations were set to the default value of 1 m−3, which is probably an effect caused by studying only relatively small cluster sizes.


Introduction
Atmospheric aerosol particles are known to have significant effects on both the global climate and human health (Pöschl, 2005). Secondary aerosols form as a result of a series of events, starting with the clustering of individual molecules and progressing ies indicate that atmospheric particle formation involves sulfuric acid and water (see, for example, Brus et al. (2011) and references therein), but these two components are not enough to explain all the observed particle formation events. Possible candidates to enhance sulfuric acid-water based particle formation are ammonia and dimethylamine (Kurtén et al., 2008). 10 Theory provides a useful tool to explore the nanometer-sized molecular clusters that are difficult to study experimentally. This paper gives an account of a method to solve the birth-death equations (BDE, the equations which describe the creation and destruction of molecular clusters by condensation and evaporation) as related to atmospheric clusters. Courtney (1962) explicitly solved the BDE for clusters up to one hundred 15 molecules of water. Nishioka and Fujita (1994) looked at the binary water/sulfuric acid system, solving the equations by using Euler's method. Wyslouzil and Wilemski (1995) examined six different systems using the Bulirsch-Stoer method. Vehkamäki et al. (1994) looked at two binary systems and solved the concentrations for the steady state using a matrix method. McGraw (1995) did something similar for the water-acid sys- 20 tem. All of the previous studies made an assumption that only monomers can collide and evaporate from the clusters, except for McGraw (1995), who reasoned that the acid hydrate is the most stable cluster (and therefore will be present in the highest concentration), so collisions with the hydrate are more important than with the bare acid monomer.

Methods
The Atmospheric Cluster Dynamics Code (ACDC) is a dynamical model to study the time development of molecular cluster distributions by explicit solution of the birth-death equations. The birth-death equations can be written as 5 where i is the cluster whose concentration is given by this equation, j is another cluster in the system, c i is the number density of cluster i , β i j is the collision coefficient between clusters i and j , γ i →j is the evaporation coefficient of a cluster i into two smaller clusters (one of which is j ), Q i is an outside source term of i , and S i includes other possible loss mechanisms for cluster i . The terms on the right hand side can be described 10 physically as the generation (birth) of clusters of type i through collisions of smaller clusters, the generation of clusters of type i through evaporation of larger clusters, the destruction (death) of clusters of type i through collisions with other clusters, the destruction of cluster i by fragmentation into smaller clusters, other creation mechanisms, and other destruction mechanisms, respectively. These other mechanisms depend on 15 the system being studied. We shall enumerate various such terms here for both neutral systems and those containing ions. All of the terms in the neutral system are also present in the ionic case, although the reverse is not true. In ACDC, the BDE are solved with the ode15s solver in the MATLAB software suite, which is effective in solving systems of stiff differential equations  ichelt, 1997). The novelty of ACDC is the generation of the equations that are fed into the MATLAB solver. Generation of the equations is essentially a series of logical checks over all possible cluster combinations to see which evaporations and collisions can create/destroy a given cluster. This work is tedious and prone to typographical errors if done by a human, but it is ideally suited to a computer code. In this case, we have 25 chosen to use the Perl scripting language to generate the equations, due to the relative 25267 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ease by which Perl handles string manipulation. This code was originally written to study the atmospherically-relevant system of sulfuric acid, ammonia, and bisulfate ion, but it has been designed to be flexible and can easily be extended to different systems (see Sect. 2.3 and Vehkamäki et al. (2011) for other examples). The equations can be modified through both an input file (which specifies the number and composition 5 of clusters) and command-line arguments (which are used to add and remove various source/sink terms in the equations themselves).
It should be noted here that Eq. (1) contains all possible evaporations and collisions in the system. Many previous studies have limited themselves to the case of only monomer collisions and evaporations. This is a valid assumption in many cases where the monomer concentrations are much higher than those of other clusters, and it greatly lowers the complexity of the resulting equations. However, given that the generation of the equations in ACDC is done by a Perl script, there is not a lot of extra effort required to include these clusters. The existence of small stable clusters in the atmosphere also suggests that non-monomer collisions may be important . 15 In addition, non-monomer evaporations might also become significant as the cluster size grows, as breaking into small stable clusters can be more energetically favorable than the evaporation of a monomer . In any event, like most of the features in ACDC, non-monomer collisions and evaporations can be included or ignored by a simple command-line argument.

Neutral clusters
The collision coefficients between clusters are taken from the kinetic gas theory, and they describe how often two spherical clusters following a Maxwell-Boltzmann velocity distribution collide with each other. For clusters i and j the collision coefficient, β i j in m 3 s −1 is given as where k b is the Boltzmann constant, T is the temperature, and m i and V i are the mass and volume of cluster i , respectively. In ACDC, the volume of a cluster is determined in a simple fashion by assuming that the compounds in the cluster are in liquid form. The molecular volume of the pure liquid is calculated from the molecular mass and saturated liquid density (at an atmospherically relevant temperature), which are input 5 parameters. The volume of the cluster is taken to be the sum of the molecular volume multiplied by the number of molecules of a given type for all molecular kinds in the cluster. While this is only an approximation, using more sophisticated methods (such as using the density of a solution with the same composition as the given cluster, or taking into account the temperature dependence of the liquid density) does not result in 10 large differences to the steady-state cluster concentrations. Given the extra complexity of including such functions for every possible cluster composition, the difficulty of doing this for charged clusters, and the relative insensitivity of the final result, ACDC retains the simplest approach to determining cluster volume. According to the condition of detailed balance, the evaporation coefficients can be 15 calculated from the collision coefficients and the free energies of formation of the mother and daughter clusters: where i and j are the daughter clusters, β i j is the collision coefficient between i and j , c e i is the equilibrium concentration of cluster i , ∆G i is the free energy of formation of 20 cluster i from the constituent monomers (which implies that all monomer free energies are zero), and c ref is the monomer concentration of the reference vapor for which the free energies were calculated. From ACDC's point of view, the origin of the free energies is not relevant, i.e. they can be calculated by any method from the classical liquid drop model to quantum chemistry. ACDC currently allows for one to enter the enthalpy often-used assumption that ∆H and ∆S remain constant over the temperature range of interest). Two major loss terms which need to be included in any realistic system are the loss of particles on the walls of the chamber and the loss of particles by collision with large aerosols not explicitly included in the simulation, i.e. coagulation. The former term 5 is only applicable in the case of laboratory experiments, while the latter is important when trying to predict the concentrations under atmospheric conditions where there are usually high populations of large pre-existing aerosols. This work assumes atmospheric conditions, and consequently the only loss term included in the simulations is the loss by coagulation. This term is derived from experimental data measured in Hyytiälä, Finland (Dal Maso et al., 2008). We use a constant coagulation sink coefficient of 2.6×10 −3 s −1 , which is the median condensation sink coefficient of sulfuric acid vapor on pre-existing aerosol particles. Using the parametrized formula from  for a cluster size dependent coagulation sink coefficient did not have a significant effect on the steady-state cluster concentrations. 15 One question that plagues the users of any simulation package in a finite system is that of boundary conditions. Boundary issues can manifest themselves in a code such as ACDC quite readily when the size of the largest tracked clusters is relatively small. This is because the code allows collisions to form large clusters which are not explicitly tracked in the system. As soon as these clusters form, they are "lost" (the material Introduction In addition to the steady-state cluster concentrations, one would also like to compute the cluster flux between different sizes and, by extension, the particle formation rate. The formation rate in our study is defined as the flux of clusters outside the system; since these clusters are not allowed to re-enter, it is as if they have become stable particles (a valid assumption if the clusters on the boundaries are large enough to 5 have negligible evaporation coefficients). However, in smaller systems, it is not clear that this should be the case. Indeed, as Ortega et al. (2011) point out, for systems with stable pre-critical clusters even clusters larger than the critical nucleus can rapidly decay through non-monomer evaporation. This issue is explored more fully in Section 3 for our particular system of interest. ACDC automatically tracks the particle formation 10 rate, and an extra command-line flag can be passed to keep track of the fluxes between clusters and the formation rates of all clusters, which makes for easy analysis after the simulation has been run.

Ionic clusters
If one is interested in ion-induced particle formation, the situation becomes a little more 15 complex compared to the neutral case. All of the above terms are still included in the system, but additional questions must now be addressed. These questions include how the ions are introduced to the system (ion source terms), how the ions disappear from the system (recombination with ions of opposite charge), and how ions collide with neutral aerosols. These terms are included here for completeness, despite that 20 the focus of this manuscript is on the neutral clusters.
A simple way to introduce ions into the system is to add a constant source term to the equation describing the concentration of the ionic form of a monomer of interest, similar to the way that other compounds (such as sulphuric acid and ammonia) are introduced to the system. This results in all ions being added to the system as 25 monomers, which can then grow by combination with neutral clusters. However, from an experimental/observational standpoint, this is not very realistic. Ionization in the atmosphere is caused by cosmic rays or radiation from radon decay (  . The concentration of sulfuric acid is much lower than the concentration of other air molecules, and consequently it is unlikely that the ion sources will be ionization of the sulfuric acid monomer (this is made even more unlikely by the fact that atmospheric sulfuric acid is mainly present as hydrates or attached to a base ). Instead of adding a source term to the ionic form of e.g. the sulfuric acid monomer, 5 ACDC has the option to introduce a new equation to the system. This equation keeps track of the concentration of a generic negative ion (it currently has the mass and molecular volume of an oxygen molecule, as the oxygen concentration in the atmosphere is many orders of magnitude above the sulfuric acid concentration).
In Eq. (4), Q ion is the ionization rate of the air (oxygen molecules), α rec is the recombination rate of positive and negative ions (taken to be the usual literature value of Israël, 1970;Bates, 1982), and the summation is over every neutral cluster in the system. In addition to this equation, every neutral cluster has a loss term and every ionic cluster has a source term corresponding to the third term on the 15 right in Eq. (4). All ionic clusters also have loss terms corresponding to the recombination term, while all neutral clusters have a similar source term (it is assumed for simplicity that recombination produces a single neutral cluster, different to the mother ionic cluster only by conversion of the bisulfate ion to a neutral sulfuric acid). It is well-known that the collision rate coefficient between ionic and neutral clusters 20 is higher than between two neutral clusters (Tammet and Kulmala, 2005). In sulphuric acid containing systems, this is due to the fact that the ion interacts strongly with the permenant dipole moment of the acid molecule, resulting in more attractive forces and a larger collision cross-section. Consequently, Eq.
(2) needs to be multiplied by an enhancement factor in the case that one of the clusters contains an ion (if both of the clusters contain ions of the same polarity, electrostatic repulsion will prevent then from colliding, so such collisions are not allowed in ACDC). The exact form of the enhancement factor is not well known, and several formulae exist. The first one that is present in ACDC is to simply multiply every ion-neutral collision (and, because of detailed balance, every evaporation of an ionic cluster) by a constant factor (taken to be equal to ten). A sensitivity analysis revealed that the exact value of the enhancement factor can have a significant impact on the results, so more realistic descriptions were 5 examined. One of these is given by Hoppel and Frick (1986), and depends on the size of the ionic cluster (using the rational that the more solvated the ion, the less impact it should have, so the value should tend towards unity as the cluster increases). A second (and more complicated) description is given by Lovejoy et al. (2004), and it depends also on the nature of the colliding cluster. Preliminary results have shown that 10 neither of these more realistic descriptions give results that are outside of the estimated uncertainty range of the results obtained from using a constant factor. Similarly, in the current model collisions of ionic clusters with the wall are enhanced by a factor of two.
In future work containing ionic clusters, this effect will be examined in more detail.

15
As with any other computational method, ACDC must be validated against a known system before the results can be analyzed in detail. That is, it must be shown to give the same answer as other methods for the same problem, or the differences must be thoroughly explained. The ordinary differential equation solvers in MATLAB (in particular, the ode15s routine used here) have already been well-tested and found to be 20 robust (Shampine and Reichelt, 1997); therefore, what remains to be validated here is the method of generating the birth-death equations, i.e. the Perl script. Validation of this script has been performed in three ways. The first way is the most straightforward: visual examination of the resulting equations. If the expected terms appear error-free in their entirety, the code is obviously computer codes. While continual inspection of this nature is certainly useful, it cannot be relied upon by itself to confirm the accuracy of the method. Consequently, two other tests have been employed. The second way is by comparison to classical nucleation theory. The equations of classical nucleation theory (CNT) give expressions for the steady-state nucleation (particle formation) rate and cluster concentrations in multi-component systems which is an accurate approximation to the solution of the birth-death equations when one considers only monomer collisions and evaporations. ACDC has been tested against CNT for the case of homogeneous water vapor-liquid nucleation and found to give differences of far less than 1 % .
A final test is the nucleation of pure aluminum at high temperatures. Li et al. (2007) have determined very accurate monomer addition free energies using Monte Carlo simulation of the equilibrium constants and several high level corrections for the aluminum dimer through the 60-mer. Girshick et al. (2009) explored the kinetics of this system by using these addition free energies and modified collision rate coefficients (Li and Truhlar, 2008) to predict the nucleation rate and steady-state concentrations of the vapor-liquid nucleation of aluminum. Even without using the modified collision rate coefficients (i.e., using Eq. (2) above), ACDC provides results in very close agreement with what is seen by Girshick et al. (2009). This is perhaps unsurprising, since Girshick et al. (2009) use CNT as a basis for their dynamics calculations, which ACDC has al-20 ready been validated against through the second test in this section, but it illustrates the flexibility of ACDC; the change from test two to test three required minor changes to the input file and the cluster free energies.

Results
When running ACDC, the first choice that needs to be made is the source of the thermodynamic data used to calculate the evaporation coefficients by Eq. (3) be either dimethylamine (DMA) or ammonia) clusters explored by Ortega et al. (2011). We wanted to include in our model only clusters for which we have commensurable high level quantum chemical data, thus restricting ourselves to clusters with a maximum of four acid and four base molecules. In the future, we will be able to add larger clusters when quantum chemical data becomes available with increased computer power, or 5 through the use of liquid drop model properties for larger clusters. Water is not included in the system, because sufficient quantum chemical data for clusters containing water, sulfuric acid, and a base are not yet available. It should be noted that Nadykto et al. (2011) have recently published results on a similar system, which could also be used for this test; however, Ortega et al. (2011) also report results for larger clusters, which reduces the boundary effects.  has also pointed out that the exchange/correlation density functional used by Nadykto et al. (2011) can significantly underestimate the stability of DMA/sulfuric acid clusters. As noted above, the operative definition for the particle formation rate in smaller systems is somewhat unclear. If it is defined as the flux of all material out of the 15 system, the rates produced in this system might be artificially inflated. This happens because some collisions occur which result in clusters not in the system, but which should (at least partially) evaporate back into the system. For example, the collision of the DMA monomer (which is present in a high concentration) with a cluster containing four DMA molecules (and any number of sulfuric acid molecules) forms a cluster that 20 contains five amine molecules and is therefore not included in the system, but might not be very stable  observed that clusters with more bases than acids are generally not stable). Nonetheless, this flux leaves the system and does not return. Because of this, we have decided to only allow clusters to leave the system by collisions where both colliding clusters contain sulfuric acid. All of the values of J 25 reported here were computed by: where the indices i and j refer to the number of sulphuric acid molecules in the first and second cluster, and k and l refer to the number of base molecules, subject to the constraint that i + j > 4 and/or k + l > 4, so that the resulting cluster's concentration is not explicitly tracked in the simulation. It should be noted here that the largest clusters in our system are approximately 1 nm in diameter. As the particular formation rates 5 given by Eq. (5) may be artificially overestimated, caution must therefore be stressed when comparing the rates presented here to true particle formation rates. Figure 1 shows the particle formation rate for sulfuric acid-ammonia and sulfuric aciddimethylamine systems as a function of acid and base concentrations in T = 298.15 K. The ranges of the sulfuric acid and ammonia monomer concentrations were chosen 10 based on typical atmospheric concentrations measured in Hyytiälä (Janson et al., 2001;Petäjä et al., 2009). Because the monomer concentration in our system is not fixed, the source terms of each molecule type need to be set so they result in the desired steady-state monomer concentrations. This was done in the following way: first, the concentrations of the sulfuric acid and base monomers were set to the target constant 15 values, i.e. the differential equations (Eq. 1) for the monomers were set to zero. After the simulation finished, the source terms were calculated from the monomer equations in the steady state, and the simulation was re-run without fixing the monomer concentrations using these calculated source terms. Finally, it was checked that the monomer concentrations have the desired values in the accuracy of 3 %. 20 When solving a series of differential equations, it is important to make sure that the simulation has reached the steady state. All the simulations were initially run for 50 000 s, and the concentrations of species at several conditions were examined to ensure they reached the steady state. In addition, all of the rates in the right graph of Fig. 1 were run for 5 000 000 s (results not shown here); the rates differed only in 25 the tenth significant figure or so, providing further evidence that 50 000 s is enough to reach the steady state.
Several conclusions can be drawn from Fig. 1. The first is that the rate is higher in the DMA system than in the ammonia system at atmospherically relevant monomer 25276 the rate by around three orders of magnitude at low acid/high DMA concentrations). At high acid concentrations, the DMA concentration does not have a large effect at least for the ranges studied here.
From Fig. 1, we can determine the conditions for the rest of the tests presented here. A typical atmospheric sulfuric acid concentration is 10 12 m −3 . If we take the DMA 10 concentration to be equal to 5×10 13 m −3 , ACDC predicts a particle formation rate of around 10 6 m −3 s −1 , which is a typical new particle formation rate seen in Hyytiälä (Dal Maso et al., 2005;Manninen et al., 2009). This DMA concentration is on the high range (but approximately the same order of magnitude) as what is measured in rural boreal areas (Ge et al., 2011). Therefore, we will use these concentrations as our "standard" 15 concentrations for examining the effect of the other parameters below. Figure 2 shows how long the concentrations in the system take to stabilize. This number is calculated by examing the concentrations throughout the entirety of the simulation, determining the smallest time at which the concentration of every species is within 0.1 % of its final concentration. As this number is computed by looking at the fi-probably not observed in the atmosphere for this system. Figure 3 shows the effect of only allowing monomer collisions and evaporations in the system. While this assumption (that non-monomer collisions and evaporations are not important) is valid under certain conditions (Arstila, 1997), it clearly will not be valid when there are stable pre-nucleation clusters . From Fig. 3, we 5 can see a rather large difference in the rates by excluding non-monomer interactions (several orders of magnitude under certain conditions). This suggests that there are stable pre-critical clusters in this system; indeed, e.g. Ortega et al. (2011) have noticed the relative stability of clusters consisting of two acids and two DMA molecules. Figure 4 shows the steady-state concentrations of all the clusters when the monomer 10 concentrations are set to the standard values described above. In Fig. 4, clusters with the highest concentrations are located around the diagonal of the acid and base number matrix. This seems sound, since clusters consisting mainly only of acid or base molecules are not expected to be stable . Of particular interest is that the concentrations stay relatively high everywhere on the diagonal, i.e. 15 larger clusters are also present in relatively high concentrations. A plot of the total flux (collisions minus evaporations) is shown in Fig. 5. These pathways were determined by the following algorithm: (1) Every flux above a certain cutoff size (10 3 m −3 s −1 ) was recorded, (2) For each flux, the path back to the monomer using the highest flux option was traced. This means that not all flux pathways are shown. 20 For example, the concentration of one DMA/one acid clusters is high enough that selfcollision to form the two amine/two acid cluster is significant; however, that collision is not included in Fig. 5 because the addition of one acid to form the two acid/one DMA cluster, followed by addition of a single DMA is a more favorable pathway. It appears that the dominant pathway involves the formation of the one DMA/one acid cluster, 25 which is very rapid (due to the high concentrations of the acid and DMA monomers). This cluster is stable enough to be present at a fairly high concentration, and could be a platform for growth into the larger sizes. It is interesting to note from this graph that the main flux out of the system is along the diagonal (equal numbers of acids and 25278 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | bases), and that significant fluxes are observed through non-monomer collisions. The aforementioned boundary effects were examined with the option that disables collisions resulting in clusters outside of the system. Figure 6 shows the differences in concentrations with and without these collisions. It has to be noted here that the upper and lower limits of the color scale are set to be the same in all the figures describing 5 differences when a certain parameter is changed. This allows for better comparison between the figures; however, it has to be kept in mind that values outside of the limits are now scaled to be equal to the limit values. A bit surprisingly, disallowing collisions from leaving the system did not have much effect on the concentrations in the chosen standard conditions (left panel of Fig. 6). This is probably due to the pathways 10 out of the system by cluster collisions being much smaller in magnitude than those going to the coagulation sink or passing between clusters inside the system (which can somewhat be seen in Fig. 5). Therefore, this option was also tested at high monomer conditions ( system. At higher concentrations, the flux out of the system from a certain cluster is now comparable to the flux to and from neighboring clusters, which means that eliminating this pathway results in a surplus of material to rearrange around the system. This, of course, results in higher cluster concentrations near the boundary. An essential parameter in studying the dynamics of a molecular system is the tem-20 perature. Figure 7 shows the differences in concentrations when the temperature is decreased from T 1 = 298 K to T 2 = 248 K. This range was chosen because it covers a broad range of tropospheric temperatures. As can be seen quite clearly in Fig. 7, decreasing the temperature increased the concentrations of almost all the non-monomer clusters. The increase in the concentrations can be explained by comparing the colli-25 sion and evaporation coefficients at both temperatures. As the temperature decreases, so do both the coefficients, but the decrease in the evaporation coefficients (where the temperature dependence is exponential) is much higher than that in the collision coefficients (where the dependence is in the square root of the temperature). Therefore, although the velocity of the particles decreases and they collide less frequently, they also stick together much more tightly, leading to more stable clusters. This is seen most dramatically in the clusters which are very unstable (the clusters which do not have around equal numbers of acids and bases). One of the major particle sinks in the system is the condensation of clusters on large 5 pre-existing aerosols (a similar effect to the condensation of clusters onto the chamber walls in laboratory experiments). The effects of increasing this sink or turning it off are shown in Fig. 8. It can be seen that the concentrations of the largest clusters are reduced by increasing the losses due to this sink. This is expected, since the flux will be attenuated with each step towards a larger cluster. It must be noted here that removing 10 the coagulation sink terms resulted in the system not reaching the steady state if the initial monomer concentrations are set to the default value of 1 m −3 . In this case, running for a simulation length of 5 × 10 4 s produced different cluster concentrations than running for 1 × 10 7 s. However, when the initial monomer concentrations are set to the desired steady state values, the system does find a steady state, which is presented in Fig. 8. In this case, the concentrations of the largest clusters are increased, which also is an expected result. Since collisions between clusters play a major part in the birth-death equations, the question of sticking probabilities naturally arises. While many kinetic codes (including ACDC) increase the sticking probabilities for ion-neutral cluster collisions (as described 20 above), to our knowledge the reverse case has not been explored in detail, i.e. when two neutral clusters collide, they will always stick together. However, there is no guarantee that this is always the case, and therefore the effect of sticking probabilities less than unity was examined. It has to be noted here that the sticking probability can also be thought to be taken into account in the evaporation (Kulmala and Wagner, 2001). 25 Figure 9 shows the difference in concentrations when the sticking probability is reduced to 0.1 for collisions involving clusters that have the highest concentrations at the standard conditions. These clusters include the monomers and clusters consisting of (1) one acid and one DMA, (2) two acids and one DMA and (3) two acids and two DMA cluster concentrations. This implies that the sticking probabilities in collisions involving the most numerous clusters can have a significant effect on the cluster distribution, and the issue of sticking probabilities should be examined in more detail in the future. In particular, the exact form of the probability should be explored (it will certainly be a function of the cluster size, for example). From Fig. 9 it can be seen that decreasing 10 the sticking probabilities reduces the concentrations of all the clusters, the effects being more pronounced for larger clusters. What is interesting here is not that the concentrations decrease if one allows for non-sticking clusters; rather, the fact that an order of magnitude decrease in the sticking probabilities of only a few clusters can reduce the concentration of other clusters (and consequently, also the particle formation rate) by 15 many orders of magnitude. This illustrates once again the highly non-linear behavior of these systems and the difficulty of predicting the extent to which small changes can manifest themselves.

Conclusions
We have presented a new program for modeling the kinetics of clusters by explicit 20 solution of the birth-death equations, using an efficient computer script for their generation and the MATLAB ode15s routine for their solution. This script, referred to as the Atmospheric Cluster Dynamics Code (ACDC), can easily use cluster free energies calculated by any method, from liquid drop theory to quantum chemical calculations. By using the recent formation free energies computed by Ortega et al. (2011) for sulfuric acid and dimethylamine containing clusters, we have examined the effect of changing various system parameters at atmospherically relevant monomer concentrations, 25281 non-monomer collisions, and monomer concentrations can all have significant effects under certain conditions. Removal of coagulation sink terms prevented the system from reaching the steady state when the initial monomer concentrations were set to 1 m −3 , which probably results from the small system size. If the starting concentrations for monomers were set to the wanted steady state values, the system was able to find Introduction