An algorithm for the numerical solution of the multivariate master equation for stochastic coalescence

In cloud modeling studies, the time evolution of droplet size distributions due to collision–coalescence events is usually modeled with the Smoluchowski coagulation equation, also known as the kinetic collection equation (KCE). However, the KCE is a deterministic equation with no stochastic fluctuations or correlations. Therefore, the full stochastic description of cloud droplet growth in a coalescing system must be obtained from the solution of the multivariate master equation, which models the evolution of the state vector for the number of droplets of a given mass. Unfortunately, due to its complexity, only limited results were obtained for certain types of kernels and monodisperse initial conditions. In this work, a novel numerical algorithm for the solution of the multivariate master equation for stochastic coalescence that works for any type of kernels, multivariate initial conditions and small system sizes is introduced. The performance of the method was seen by comparing the numerically calculated particle mass spectrum with analytical solutions of the master equation obtained for the constant and sum kernels. Correlation coefficients were calculated for the turbulent hydrodynamic kernel, and true stochastic averages were compared with numerical solutions of the kinetic collection equation for that case. The results for collection kernels depending on droplet mass demonstrates that the magnitudes of correlations are significant and must be taken into account when modeling the evolution of a finite volume coalescing system.


Introduction
The evolution of the size distribution of coalescing particles has often been described by the kinetic collection (hereafter KCE) or Smoluchowski coagulation equation, known under a number of names ("stochastic collection", "coalescence").The discrete form of this equation has the form (Pruppacher and Klett, 1997) where N (i, t) is the average number of droplets with mass x i , and K(i, j ) is the collection kernel related to the probability of coalescence of two droplets of masses x i and x j .In Eq. ( 1), the time rate of change of the average number of droplets with mass x i is determined as the difference between two terms: the first term describes the average rate of production of droplets of mass x i due to coalescence between pairs of drops whose masses add up to mass x i , and the second term describes the average rate of depletion of droplets with mass x i due to their collisions and coalescence with other droplets.
Within the kinetic approach (Eq.1), it is assumed that fluctuations are negligibly small.This assumption can only be correct if the volume and the number of particles are infinitely large.An alternative approach considers the coalescence process in a system of finite number of particles, with fluctuations that are no longer negligible.This finitevolume description is intrinsically stochastic and has been pioneered by Marcus (1968) and Bayewitz et al. (1974) and studied in detailed by Lushnikov (1978Lushnikov ( , 2004) ) and Tanaka and Nakazawa (1993).
Within the finite volume description a system of particles whose total mass is M T is considered.The mass distribution of the particles is described by giving the number n i of particles with mass i, i.e., n 1 , n 2 , n 3 , . . ., n N .Then, the state of the mass distribution of the particle system is described by the N-dimensional state vector n = (n 1 , n 2 , . .., n N ).The time evolution of the joint probability P (n 1 , n 2 , . .., n N ; t) that the system is in state n = (n 1 , n 2 , . .., n N ) at time t is calculated according to the equation (Tanaka and Nakazawa, 1993) K(i, j )(n i + 1)(n j + 1) × P (. .., n i + 1, . .., n j + 1, . .., n i+j − 1, . ..; t) × P (. .., n i + 2, . .., n 2i − 1, . ..; t) K(i, j )n i n j P ( n; t) The master Eq. ( 2) is a gain-loss equation for the probability of each state n = (n 1 , n 2 , . .., n N ).The sum of the first two terms is the gain due to transition from other states, and the sum of the last two terms is the loss due to transitions into other states.The gain terms show that the system may be reached from any state with an i-mer and a j -mer more, and one (i+j )-mer less.In Eq. ( 2) K(i, j ) is the collection kernel and the transition rates are K(i, j )(n i + 1)(n j + 1) if i = j and K(i, i)(n i +1)(n i +2) if i = j .From conservation of the total probability, P ( n; t) must satisfy the relation where the sum is taken over all states.Moreover, the total mass M T of the system must be conserved, and the particle number n i should be non-negative for any mass x i : Exact solutions of Eq. ( 2) are only known for a limited number of cases (constant, sum and product kernels) and for monodisperse initial conditions.For these special cases the master equation has been solved by Lushnikov (1978Lushnikov ( , 2004) ) and Tanaka and Nakazawa (1993) in terms of the generating function of P ( n; t).For general, multidisperse initial conditions, the solution of Eq. ( 2) is not known.Additionally, for stochastic coagulation, approximate solutions were calculated by using Van Kampen's system size expansion or expansion (Van Dongen and Ernst, 1987;Van Dongen, 1987) which permits finding solutions of Eq. (2) valid in the limit of a large system.However, the system size expansion gives less reliable results when applied to systems with a low number of particles or small volumes.
Then, in order to obtain solutions for more realistic kernels (Brownian motion, differential sedimentation, etc.), a small number of particles and general multidisperse initial conditions, it has to be solved numerically.In this paper, we present an algorithm that can be applied to obtain the solution of Eq. ( 2) for any type of kernel and initial conditions.By applying this method, numerical solutions of the master equation were obtained for realistic kernels relevant to cloud physics, along with calculation of the correlations for the number of droplets for different sizes.
It is worth mentioning that the stochastic simulation algorithm (SSA) developed by Gillespie (1975) also accurately reproduces the master equation.In Gillespie's method, the master equation is not solved directly, but a statistically correct trajectory (possible solution) of the master equation is generated.At any time, expected values at each droplet size can be obtained by averaging over many runs.However, a large number of realizations are necessary in order to obtain the desired accuracy at the large end of the droplet size distribution.A detailed comparison between the two methods will be made in Sect.3.
The problem of calculating correlation coefficients was also addressed by Wang et al. (2006), who derived what they called the "true stochastic collection equation" (TSCE), which is a mean field equation at the first order and contains correlations among instantaneous droplets of different sizes.The problem with this equation and similar ones is that the rate of change of moments of order n depends on moments of order (n + 1), as was remarked by Marcus (1968).
In our work, we overcome this drawback by calculating the true stochastic averages directly from the solution of the master equation.The main idea is to reduce the dimensionality by restricting the state space only to those states which have a finite probability of being accessed.It turns out that this provides a considerable improvement in numerical efficiency.
The paper is organized as follows: in Sect.2, the numerical algorithm is explained in detail.Numerical solutions for the sum and constant kernels with a comparison with analytical solutions and with the method of Gillespie (1975) are presented in Sect.3. The numerical results for mass-dependent kernels along with calculation of correlations for different droplet sizes are presented in Sect. 4. Finally, in Sect. 5 we briefly discuss the results and the possible applications of the numerical algorithm.

The numerical algorithm
To solve Eq. (2) by brute force, the joint probability P (n 1 , n 2 , . .., n N ; t) must be discretized into a multidimensional array.The main drawback of this approach is its susceptibility to the curse of dimensionality (Bellman, 1961), i.e., the exponential growth in memory and computational requirements in the number of problem dimensions.

Calculation of all possible states
Instead of the brute force discretization of the multidimensional joint probability distribution, the solution for this problem lies on the generation of all possible states from an initial configuration, and the posterior calculation of the time evolution of the probability P ( n; t) for each generated configuration by using the master equation.From an arbitrary initial condition P (n 01 , n 02 , . .., n 0N ; 0) = 1 all possible states can be generated numerically.This can be performed by taking into account that the only transitions allowed are of the form n(+) 1 → n1 if i = j and n(+) 2 → n2 if i = j , where n(+) 1 , n1 and n(+) 2 , n2 are the state vectors: n1 = (n 1 , . .., n i , . .., n j , . .., n i+j , . .., n N ), n2 = (n 1 , . .., n i , . .., n 2i , . .., n N ). (5d) For a system consisting of N monomers at t = 0, R(N ) states (or N-dimensional vectors) can be realized, where R(N ) is the number of solutions in integers n of the Eq.(4) for conservation of mass.The number of possible configurations can be approximated from the equation (Hall, 1967) Note that, although R(N ) increases very quickly with N (for example, R(50) = 217 590 and R(100) = 190 569 232), a number of states that is manageable with an average computer is obtained (compare with the 50-dimensional array with 1.34 × 10 16 elements required for N = 50).Although Eq. ( 6) slightly overestimates the number of states, it gives estimates that can be used in order to check the performance of the algorithm.For N = 6, 10, 20, and 30 we obtained 11, 42, 627, and 5604 with the numerical algorithm, and 13, 48, 692, and 6078 by using Eq.(6).As an example, the 11 possible configurations generated from the initial state (6, 0, 0, 0, 0, 0) are displayed in Fig. 1.

Time evolution of the probabilities P (n; t) for each state
At t 0 = 0 for the initial state P (n 01 , n 02 , n 03 , n 04 , . .., ; t 0 ) = 1, and the probabilities for the rest of the states are set equal to 0. The probabilities of all generated configurations are up- dated according to the first-order finite difference scheme: It is clear from Eq. ( 7) that the state probabilities P ( n; t 0 + t) at t = t 0 + t will increase if the states from which transitions are allowed have a non-zero probability at t = t 0 (second and third terms in the right-hand side of Eq. 6) and will decrease due to collisions of particles from the same state at t = t 0 (fourth and fifth terms in the right-hand side of Eq. 7) if P ( n; t 0 ) is positive.The finite difference equation for P (1, 0, 0, 0, 1, 0) was written to illustrate the method.
P (1, 0, 0, 0, 1, 0; t 0 + t) is calculated from the equation In the second term on the right-hand side of Eq. ( 8), n 2 + 1 and n 3 + 1 are set equal to 1, as they are the number of particles in the second and third bins in the configuration (1, 1, 1, 0, 0, 0) at t = t 0 .In the third term, n 1 + 1 = 2 and n 4 + 1 = 1 as they are defined from the state (2, 0, 0, 1, 0, 0) and, finally, n 1 = n 5 = 1 in the fourth and last term.As an exercise, the time evolution of each state probability was calculated for the coalescence kernel K(i, j ) = (i 1/2 +j 1/2 )/40 from Marcus (1968).The results for 5 of the 11 possible configurations are displayed in Fig. 2.

Calculation of the expected values of the number of particles for each particle mass
The number of particles for a given mass n 1 , n 2 , . . ., n N are discrete random variables whose probability distributions can be obtained from Usually, the numerical implementation of Eq. ( 9) would involve calculating the sum of all elements of a multidimensional array, which is computationally very expensive.
Our approach is simpler: once the probabilities of all possible states are determined for all times, P (n, m; t) can be calculated just by summing over all states that have n m = n: P (n, m; t) = all states with n m = n, The expected values n m for the number of particles of mass m are then calculated from the equation As an example, for the system from Fig. 1, the probability distribution P (n, 1; t) of having n particles with mass m = 1 is displayed in Table 1.
3 Comparison with analytical solutions and the stochastic simulation algorithm (SSA) of Gillespie

Comparison with analytical solutions
The expected values for each particle mass calculated with the numerical algorithm were tested against the analytical solutions of the master equation reported in Tanaka and Nakazawa (1993) for the constant (Eq.12) and sum (Eq.13) kernels (K(i, j ) = A, K(i, j ) = B(x i + x j )) obtained for the monodisperse initial condition P (N 0 , 0, 0, . .., 0; 0) = 1.They are In Eqs. ( 12) and ( 13), N 0 is the initial number of particles, C N 0 m is the binomial coefficient and n m values are the true stochastic averages for each particle mass m at time t.In Eq. ( 12) τ = AN 0 t, where A = 1.2×10 −4 cm 3 s −1 is the constant collection kernel.Finally, in Eq. ( 13), T = BN 0 v 0 t, where v 0 is the initial volume of droplets and B = 8.82 × 10 2 cm 3 g −1 s −1 .Turning to a concrete numerical example, the evolution of a cloud system with an initial monodisperse droplet size distribution of N 0 = 10 droplets of 10 µm in radius (droplet mass 4.189 × 10 −9 g) at t 0 , and a volume of 1 cm 3 was calculated with the numerical algorithm.The time Table 1.Probability distribution P (n, 1; t) of finding n particles of size m = 1 at time t, for a system with the initial condition P (6, 0, 0, 0, 0, 0; 0) = 1.
step was set equal to t = 0.1 s.Comparisons between the numerical and analytical results for both the sum and constant kernels at t = 1200 s are shown in Figs. 3 and 4 with an excellent agreement between the two approaches.

Comparison with the SSA of Gillespie
As was mentioned in the introduction, the algorithm of Gillespie generates a statistically correct trajectory of the stochastic master equation.It was presented in Gillespie (1975), and popularized in Gillespie (1977) were it was used to simulate chemical systems.As we know, in Gillespie's SSA, the ensemble mean for the number of droplets at each droplet mass is calculated from the expression (Gillespie, 1975) where N r is the number of realizations of the stochastic algorithm, N i (m; t) is the number of droplets of mass m in the i realization at time t, and N (m; t) is the ensemble mean.From expression (14) it is clear that in order to obtain the correct expected values (N (m; t)) at the large end of the droplet size distribution, we will need a large number of realizations of the SSA.
To further investigate this question, the evolution of a cloud system with an initial monodisperse droplet size distribution of N 0 = 30 droplets of 14 µm in radius (droplet mass 1.1494 × 10 −8 g) at t 0 , and a volume of 1 cm 3 was calculated with both the numerical algorithm and Gillespie's SSA for the sum kernel (K(i, j ) = B(x i + x j ), with B = 8.82 × 10 2 cm 3 g −1 s −1 ).The results obtained by the two methods were then compared with the analytical solution of the master equation (Eq.13) obtained by Tanaka and Nakazawa (1993) for the same conditions.The averages calculated from Gillespie's method for N r = 10 3 realizations and the analytical solution at t = 1200 are displayed in Fig. 5.As can be observed, both the Monte Carlo averages and the analytical solution are closely coincident for the small end of the droplet size distribution.However, due to the small number of realizations, the SSA fails to reproduce the distribution for the expected values at the large end (see Table 1).
For a more detailed analysis, the expected number of particles for each droplet size calculated from the analytical solution, the numerical algorithm and the SSA of Gillespie (for 1000 and 10 000 realizations) are displayed in Table 1.As can be seen in the table, the size distributions are almost identical for the small end.However, they differ substantially at the large end since the SSA produces no particles larger than 12 v 0 and 16 v 0 for 1000 and 10 000 realizations, respectively (v 0 = 1.1494 × 10 −8 g, mass of a 14 µm droplet).
For 1000 realizations, the Monte Carlo averages differ from the analytical solution for bin numbers larger than 8.For 10 000 realizations we have the same situation for bin numbers larger than 13.
As expected, for 1000 and 10 000 realizations, no states with droplets 30 times larger than monomer-sized ones were realized.The numerical algorithm described in this paper performed very well at the large end, with expected values that are very close to the analytical solution (see Fig. 5 and Table 3).
It can be concluded that our method will be suitable if we need to accurately calculate the large end of the droplet spectrum for small systems (with < 50 monomer droplets in the initial state).As the SSA requires a large number of realizations, it will be computationally very expensive.Then, for a small number of particles, our algorithm will be a good alternative, as it provides the desired accuracy to detect the possible small differences between different numerical approaches.It can also work as a benchmark for different Monte Carlo methods for the collision-coalescence process.

Numerical calculation of correlation coefficients
The evolution equation for the expected values of the random variables can be obtained by multiplying Eq. ( 7) by n k and summing over all states (see Bayewitz et al., 1974): The KCE is obtained from Eq. ( 15) by assuming that n i n j = n i n j , i.e., that the correlation between the random variables is zero.A form of Eq. ( 15) was deduced in Tanaka and Nakazawa (1993) and in Wang et al. (2006) for a general type of kernel.Bayewitz et al. (1974) have quantified the deviation of the size distributions calculated with the KCE from the exact distribution obtained from the master equation for a constant kernel.From Eq. ( 15) it can be concluded that as long as the correlations remain appreciable, the results of the KCE will not match the true stochastic averages.The correlation (or correlation coefficient) between two random variables n i and n j denoted as ρ i,j is In Eq. ( 16), the covariance (cov(n i , n j )) is calculated according to Where E n i n j is the expected value of the product n i n j which, for the bivariate case, is In Eq. ( 18), f (n i , n j ) is the two-dimensional joint probability mass function (pmf) which was calculated similarly to how it was done in the univariate case (see Eq. 10): f (n, l) = P (n, i; l, j ; t) = all states with n i = n and n j = l, P n 1 , n 2 , . .., n i = n, . .., n j = l, . .., n N ; t .
In the former equation, P (n, i; l, j ; t) is the probability of having n droplets of mass i and l droplets of mass j .

Numerical results for the constant, sum and product kernels
Correlation coefficients (ρ 1,2 and ρ 2,3 ) were obtained by Wang et al. (2006) using the analytical solution obtained by Bayewitz et al. (1974) for a constant collection kernel.They found that, even for this case, the magnitude of correlations could be quite large.We will extend their analysis by calculating the time evolution of the correlation coefficients ρ 1,2 and ρ 2,3 for the constant, sum and product kernels (see Fig. 6).For each case, the simulations were conducted for two systems containing 10 and 40 droplets of 14 µm in radius, respectively, and a volume of 1 cm −3 .As can be observed in the figure, in all the cases we have non-zero correlations.From the evolution of ρ 1,2 for all the kernels, we can infer that the random variables n 1 and n 2 are, at the beginning of the simulation, strongly anticorrelated.This is due to the fact that in the initial stage of evolution of the system we have mainly collisions between size 1 droplets to form size 2 droplets.On the other hand, the random variables n  to collisions with size 1 droplets will increase the number of size 3 droplets (Wang et al., 2006).
At t = 1800 s, the true stochastic averages (see Eq. 11) obtained numerically from the master equation are displayed in Figure 9.Time evolution of the correlation coefficients ρ 1,3 and ρ 1,5 for a 1 cm 3 system modeled with the turbulent hydrodynamic kernel and containing initially 20 droplets of 14 µm and 10 droplets of 17 µm.
Table 2. Analytical size distributions of the kinetic collection equation (KCE) calculated with monodisperse initial conditions.
T = AN 0 t Note: parameters β, B and C are constants, x and y are the masses of the colliding drops.N 0 is the initial concentration and v 0 is the initial volume of droplets.The index i represents the bin size.Fig. 7, together with the mean values for each droplet mass calculated from the analytical solutions of the KCE (see Table 2).For the three cases, at the large end of the spectrum, results differ substantially.This is in agreement with the analytical study of Tanaka and Nakazawa (1994), who demonstrated that the true stochastic averages coincide well with those obtained from the kinetic collection Eq. ( 1) if the bin mass k satisfies the inequality k 2 M 0 , where M 0 is the total mass of the system.

Numerical results for the turbulent hydrodynamic collection kernel
Collisions between droplets under pure gravity conditions are simulated with a collection kernel of the form The hydrodynamic kernel (Eq.19) does not take into account the turbulence effects and considers that droplets with different masses (x i and x j and corresponding radii r i and r j ) have different settling velocities.In Eq. ( 20), E(x i , x j ) are the collection efficiencies calculated according to Hall (1980).In turbulent air, the hydrodynamic kernel should be enhanced due to an increase in relative velocity between droplets (transport effect) and an increase in the collision efficiency (the drop hydrodynamic interaction).These effects were taken into account by implementing the turbulenceinduced collision enhancement factor P Turb (x i , x j ) calculated in Pinsky et al. (2008) for a cumulonimbus cloud with dissipation rate, ε = 0.1 m 2 s −3 , and Reynolds number, Re λ = 2 × 10 4 for cloud droplets with radii ≤ 21 µm.Then, the turbulent collection kernel has the form K Turb (x i , x j ) = P Turb (x i , x j )K g (x i , x j ). (21) In the simulation for turbulent air, a system corresponding to a cloud volume of 1 cm 3 and a bidisperse droplet distribution was considered: 20 droplets of 14 µm in radius and another 10 droplets of 17.64 µm in radius, corresponding to a liquid water content (LWC) of 0.436 g m −3 .For the turbulent collection kernel the true stochastic averages at t = 1200 and 1800 s are displayed in Fig. 8, and compared with the mean values for each droplet mass calculated numerically from the KCE with kernel (20).Also, for this case, at the large end of the spectrum, results obtained from the KCE differ substantially from the stochastic means.The time evolution of the correlation coefficients ρ 1,5 and ρ 1,3 displayed in Fig. 9 confirms the fact that correlations cannot be neglected.Finally, the time variations of n 1 , n 5 , n 15 and n 20 were calculated and compared with the time evolution of the averages calculated from the KCE with the same initial conditions and coalescence rate.We can see from Fig. 10 that for the small masses k = 1 and 5, both solutions are closely coincident up to 1800 s, and that for the larger masses k = 15 and 20, the results are different at all times.

Discussion and conclusions
The full stochastic description of the growth of cloud droplets in a coalescing system is a challenging problem.For finite volume systems or in systems of small populations, statistical fluctuations become important and the mathematical description relies on the master equation which has analytical solutions for a limited number of cases.In an effort to solve this problem, we have introduced a new approach to numerically calculate the solution of the coalescence multivariate master equation that works for any type of kernel and initial conditions.
For the constant, sum and product kernels, the true stochastic averages calculated numerically were compared with analytical solutions of the master equation, with an excellent agreement between the two approaches.
A numerical procedure to calculate the correlation coefficients was implemented, which were calculated for mass-dependent kernels (sum, product, and kernels modified by turbulent processes).Also numerical solutions of the master equation for bivariate initial conditions and collection kernels modified by turbulent processes were obtained and compared with size distributions obtained from the numerical integration of the KCE.The two equations give different values at the large end of the droplet size distribution.It was also shown that, for small k, the true stochastic averages n k and the solution of the KCE are closely coincident up to 1800 s.For larger masses, the results are different at all times.
A topic of discussion can be the limits of applicability of the finite volume approach to problems of precipitation formation, since such small volumes would not remain undisturbed for a long time in a real cloud.However, in defense of the finite system approach, it might be argued that in the early stages of cloud development, due to small terminal velocities of the droplets, the coalescence process is a fairly localized process; i.e., two droplets in widely separated parts of the cloud are not going to be coalescing with each other.This was the approach followed by Bayewitz et al. (1974) (and endorsed in Gillespie, 1975).In their paper, for comparing the stochastic and kinetic approaches, they partitioned the cloud into many sub-volumes, with no collisions being permitted for two droplets of different sub-volumes.However, interactions between sub-volumes through sedimentation, diffusion or other physical processes were not considered.
For a constant collection kernel, a more complex model that uses the master equation formalism and introduces the interactions between the sub-volumes was developed by Merkulovich andStepanov (1990, 1991).This model is based on a scheme proposed by Nicolis and Prigogine (1977) for chemical reactions.Within this theory, the whole system is subdivided into sub-volumes (coalescence cells) that can be considered spatially homogeneous.Coalescence events are permitted only between droplets from the same subvolume, and interactions between neighbors occur through the diffusion process.That leads to a set of master equations for each sub-volume.Although very complex, it could be a 12325 Table 3. Expected values for each droplet mass obtained at t = 1200 s for the analytical solution, the numerical algorithm proposed in this work, and Gillespie's SSA (for N r = 1000, 10 000 realizations).Calculations were performed with the initial condition P (30, 0, 0, 0, . .., 0; 0) = 1, and the sum kernel K(i, j ) = B(x i + x j ) with B = 8.82 × 10 2 cm 3 g −1 s −1 .
Expected values for each droplet size: < n i >, t = 1200 s Bin number Analytical solution Numerical algorithm SSA SSA (N r = 1000) (N r = 10 000) starting point in order to consider the interactions between small coalescence volumes through sedimentation or other physical mechanisms.However, fluctuations will be also very important, if the collection kernel K(i, j ) increases sufficiently rapidly with i and j and a giant droplet with mass comparable to the total mass of the system is formed.In that case, the total mass predicted by the KCE starts to decrease.This is usually interpreted to mean that the system exhibits a phase transition (also called gelation).After this moment, the true averages calculated from the master equation will differ from the averages obtained from Eq. ( 1) and there is a transition from a system with a continuous droplet distribution to one with a continuous distribution plus a giant cluster (Alfonso et al., 2013).After the sol-gel transition the KCE breaks down: the second moment of the size distribution diverges at the gel point and, as was remarked, the first moment decays, i.e., mass is not conserved.
The limitation of the KCE equation arises from the fact that it is a deterministic equation with no fluctuations or correlations included.Then it describes an inherently stochastic process with a single metric, the mean cluster distribution (Matsoukas, 2015).Then, in order to model properly the system behavior after the giant cluster is formed, the role of fluctuations should be considered.
By using the finite volume approach, the expected values at the large end of the droplet size distribution can be obtained in the post-gel region (Lushnikov, 2004;Matsoukas, 2015) and be compared with the expected values obtained from the kinetic approach.As a result, it is expected to obtain broader droplet mass distributions by using the stochastic approach.A follow-up paper will be devoted to a more detailed analysis of all these problems.

) 12321 Figure 6 .
Figure 6.Time evolution of the correlation coefficients ρ 1,2 and ρ 2,3 for the constant, sum and product kernels (in panels a, b and c, respectively) for two systems with a volume of 1 cm −3 and containing 10 and 40 droplets of 14 µm.

Figure 7 .
Figure 7.Comparison of the size distributions obtained from the stochastic master equation (solid line) with that of the KCE (dashed line) at t = 1800 s for a 1 cm 3 system containing initially 40 droplets of 14 µm.The expectation values are shown for the constant, sum and product kernels (in panels a, b and c, respectively).For the small end the size distributions are closely coincident, for the large end the two equations give different values.

Figure 8 .
Figure8.For the turbulent hydrodynamic kernel, comparison of the size distributions obtained from the stochastic master equation (solid line) with that of the KCE (dashed line) at t = 1200 and 1800 s for a 1 cm 3 system containing initially 20 droplets of 14 µm and 10 droplets of 17 µm.For the small end the size distributions are closely coincident, for the large end the two equations give different values.

Figure 10 .
Figure 10.For the turbulent hydrodynamic kernel, comparison of the expected values n k obtained from the stochastic master equation (solid line) with that of the KCE (dashed line), for a 1 cm 3 system containing initially 20 droplets of 14 µm and 10 droplets of 17 µm.The time evolution of the expected values are shown for k = 1, 5, 15 and 20 (panels a, b, c and d, respectively).For the small masses k = 1 and 5, both solutions are closely coincident up to 1800 s.For the larger masses k = 15 and 20, the results are different at all times.