CONVERGENCE OF A GENERALIZED WEIGHTED FLOW ALGORITHM FOR STOCHASTIC PARTICLE COAGULATION

We introduce a general family of Weighted Flow Algorithms for simulating particle coagulation, generate a method to optimally tune these methods, and prove their consistency and convergence under general assumptions. These methods are especially effective when the size distribution of the particle population spans many orders of magnitude, or in cases where the concentration of those particles that significantly drive the population evolution is small relative to the background density. We also present a family of simulations demonstrating the efficacy of the method. 1. Weighted methods for coagulation equations. Predicting the evolution of particle size distributions that are multi-dimensional in composition space and undergo coagulation is a challenge in many areas in science and engineering. Examples include the modeling of aerosol particles in the Earth’s atmosphere [27, 33, 38, 46], and of aerosol processes in industrial applications [10]. Recently, particle-resolved models have been introduced to tackle this problem [9, 32, 39]. These stochastic models fully resolve the multi-dimensional particle distribution as they simulate a representative group of particles distributed in composition space. Coagulation, condensation and evaporation, and other important processes can then be treated on an individual particle level. The first rigorous application of such a Monte Carlo approach was suggested by Gillespie, who developed the exact Stochastic Simulation Algorithm [16, 17, 18, 19] to treat the stochastic collision-coalescence process in clouds. A straight-forward implementation of these methods assumes that each computational particle represents a physical particle. However, this naive approach can 2010 Mathematics Subject Classification. Primary: 65C35, 60K35; Secondary: 60J28.

1. Weighted methods for coagulation equations.Predicting the evolution of particle size distributions that are multi-dimensional in composition space and undergo coagulation is a challenge in many areas in science and engineering.Examples include the modeling of aerosol particles in the Earth's atmosphere [27,33,38,46], and of aerosol processes in industrial applications [10].
Recently, particle-resolved models have been introduced to tackle this problem [9,32,39].These stochastic models fully resolve the multi-dimensional particle distribution as they simulate a representative group of particles distributed in composition space.Coagulation, condensation and evaporation, and other important processes can then be treated on an individual particle level.The first rigorous application of such a Monte Carlo approach was suggested by Gillespie, who developed the exact Stochastic Simulation Algorithm [16,17,18,19] to treat the stochastic collision-coalescence process in clouds.
A straight-forward implementation of these methods assumes that each computational particle represents a physical particle.However, this naive approach can quickly become computationally infeasible, the reason being the diversity of the underlying particle population-either with respect to size or to composition-and the relative abundance of the sub-populations.There are two distinct conceptual scenarios in which these problems arise, which we term the "size multiscaling problem" and the "compositional multiscaling problem".
For the "size multiscaling problem", consider the example of initially monodisperse aerosol particles of sub-micrometer size coagulating due to Brownian diffusion in the atmosphere [26,Chapter 16.4.2].After some time of coagulation, the range of particle sizes present in the population can span many orders of magnitude, from tens of nanometers to several microns in diameter.Resolving the full particle size range is important since the small particles dominate the particle number concentration whereas the large particles dominate the particle mass concentration, and one typically aims to keep track of both.Moreover, the most likely coagulation events involve interactions of small and large particles, since the interactions are more likely to take place between particles with different diffusion properties [40,Chapter 13.3.1.3].However, the number concentrations of particles in the submicron size range are typically several orders of magnitude larger than the particles in the super-micron size range [40,Chapter 8.2].Therefore, to resolve the entire size spectrum accurately, the total number of computational particles in the simulation must be very large, otherwise the population of large particles will be under-resolved and the accuracy of the simulation significantly compromised.We will address this very problem in Section 6, where we simulate a cloud droplet population undergoing collision-coalescence, producing very few, but for the dynamics of the system crucially important, large drizzle drops [29].
A conceptually different problem is the "compositional multiscaling problem"; this arises when a certain type of particle is of special interest, and particles of this type are much rarer compared to the rest of the particles in the given size range.For example, the ability of an atmospheric aerosol population to absorb solar radiation is directly related to the aerosol impact on climate and hence is an important quantity [40,Chapter 24.2].One of the few absorbing aerosol types is soot, a by-product of the incomplete combustion of carbon containing material [36].Thus, to accurately predict how much solar radiation the aerosol population absorbs, the soot sub-population must be adequately captured.The size distribution of soot aerosol peaks in the accumulation mode range [22], a size range where many other particle types are present that do not contain soot [21,37].If the conditions are such that the number concentration of the soot sub-population is small compared to the non-soot particles, either a large number of computational particles needs to be expended to resolve the soot particles or the prediction of absorptivity will be inaccurate.
What this has in common with the previous example, however, is that there is an important sub-population of particles that have concentrations that are several orders of magnitude lower than other sub-populations.A naive representation of the particle population either leads to under-resolving the small but important subpopulation, or to requiring very large numbers of computational particles in the simulation.One strategy to circumvent these difficulties is to use the concept of "weighted" computational particles, where a single computational particle corresponds to some number of real particles.This idea has been pursued by several groups in the past; the methods can be broadly divided into two categories, either the weightings are attached to the particles or they are functions of the particle constituents.Variants of both particle-attached [35,30,47,24,25,48,49] and particle-function weightings [1,6,11,12,13,14,28,41,44] have been widely studied.
Building on this literature, and particularly on the work of Kolodko and Sabelfeld [28], we presented the development and application of the Weighted Flow Algorithms (WFA) for particle-resolved methods in DeVille et al. [7], as a generalized way of weighting computational particles to improve computational efficiency.Our analysis was able to demonstrate general conditions for consistency and optimality of this class of weighted methods.Applying this to the simulation of atmospheric aerosol particles, we used weighting functions that were power laws in the total particle mass, and showed that this can greatly increase the accuracy and performance of the simulation.One particular feature of these simulations was their tunability that allowed them to be optimized for different observables, e.g., one could choose parameters to make a simulation accurate in measuring particle number, and then change these parameters to have another simulation optimize for particle mass.Thus we have a family of simulations of such systems, each of which optimizes different observables to a different degree.To obtain a full simulation that takes into account several different observables (e.g., if we would like to track both particle number and mass well) it was necessary to perform several independent simulations, each optimized for their particular observable, and then combine them a posteriori.While this was done in DeVille et al. [7], a better method would be one that combines these different weightings during the running of the simulation.
We can now describe the content of this paper as follows: we have developed a method that solves both the size and compositional multiscaling problems at once, and also allows for a run-time mixing of multiple schemes.To do so, we have developed a particle method that allows for several different sub-populations that can interact; these sub-populations can correspond to populations at different length scales (speaking to the size multiscaling problem) or of different types (speaking to the compositional multiscaling problem).Each of these sub-populations can be weighted differently as the problem requires.
We also prove convergence of such methods under very weak assumptions.A more precise statement is included in the theorems below, but in essence we show that there are a family of "consistency conditions" that, when satisfied, will guarantee that a scheme developed with this method will converge to the correct solution in the limit of large numbers of particles.These conditions are easy to verify for a particular scheme.Moreover, we show how to develop an optimal method in the set of consistent methods.
Finally, we present a series of numerical simulations demonstrating the practical efficacy of the multiscaling methods developed in this paper.
2. Problem formulation and summary of results.For aerosol dynamics simulations, one of the most expensive portions of the calculations can be the parts that deal with coagulation of particles [15,16].The equations that govern coagulation are the Smoluchowski equations.
where K(µ, ν) is the coagulation rate between particles of masses µ, ν ∈ R n and we write [0, µ] to denote the set of points that are elementwise between 0 and µ.Now choose any partition of [0, ∞), which we denote {I n } ∞ n=0 .Then if we define it is not hard to see that (2.1) is asymptotically equivalent to where the asymptotic correction goes to zero as max j |I j | → 0. It could be debated whether we should consider (2.1) or (2.2) as the fundamental object to consider.
In this paper, we consider (2.2) as the fundamental object and consider how to simulate it stochastically.
2.2.Summary of results of the paper.We first develop a method of particle simulation for the Smoulchowski equation with the added feature that we are allowed to consider multiple, and possibly interacting, particle sub-populations.It is not surprising that one must choose a method so that the mean of the stochastic method corresponds to the evolution we are trying to approximate; however, we show below in Theorem 1 that this is also sufficient to obtain the correct mean for trajectories.
We then develop a family of consistency conditions in Section 4 that give the correct mean during the simulation.It turns out that, due to the many degrees of freedom in the parameters of the stochastic simulation, there are many choices that give consistency, so we also develop the conditions that give the optimal consistent method.Finally, we state and prove a family of theorems that guarantee that these methods are accurate in the limit of large population size.
3. Formal overview of algorithm.See Table 1 for a summary of the notation used throughout the paper.
3.1.Computational coagulation.We want to represent computational particles by both a mass and a tag representing their "class".We denote M as the set of classes and assume that masses are integer-valued.Thus the state of the system is the function Q a i (t) = "the number of computational particles of class a and mass i at time t." In this context, we use the word class very broadly: all it means is that we have made an a priori division of the particles.As can be seen from (3.1), each particle belongs to exactly one class, so the class can be regarded as an arbitrary additional particle property.The class is not a function of the other particle properties.The criterion for assigning particles to classes will depend on the application.As can be seen in Section 6, examples include: classes delineated by mass, with one class being large particles and the other small particles; classes delineated by composition, with one class containing the more common chemical constituents and the other class Variables Meaning Range a, b, c, d 1. Notation and variable ranges.containing the rare chemical constituents for examples of classes.It may even be that otherwise identical particles are divided between two classes.
We want to consider a splitting method to simulate (2.2) where we use a two-step process to determine which event to perform next.The fundamental computation is the propensity of the pair of particles (a, i) and (b, j) to coagulate and become the particle (c, i + j).The two steps are: 1. Choose a triplet of particles ((a, i), (b, j), (c, i+j)) with some probability, where (a, i) and (b, j) are represented by some positive number of computational particles, i.e., q a i (t) > 0, q b j (t) > 0. The rate at which we choose this triplet will be denoted T abc ij .2. Next, we decide what to do with this triplet; we can either remove an (a, i) particle or not, remove a (b, j) particle or not, and create a (c, i+j) particle or not.We thus define the eight probabilities p abc ij (α, β, γ), where α, β, γ ∈ {0, 1} denote whether the particle of type a, b, or c exists after the coagulation.We require, for all a, b, c ∈ M and i, j In the standard (what we might call "non-split") particle method, we would only allow α = β = 0 and γ = 1, meaning that when the triplet is chosen, we always destroy the first two particles and create the third.We could then interpret T abc ij to be the rate of coagulating (a, i) and (b, j) and obtaining (c, i + j).In this split method, we can see that the propensity to create the particle of type (c, i+j) in this coagulation is now given by the product T abc ij α,β∈{0,1} p abc ij (α, β, 1), and similarly for the destruction of the other two particles.
One note here about "split" and "non-split" methods.The reader has no doubt noticed that we are allowing for transitions that, e.g., violate mass conservation, because we can create the (c, i + j) particle without destroying the (a, i) and (b, j) particles.If we were interested in simulating a finite population of particles, this would be a problem.However, this algorithm is intended to generate a finite particle estimator of an infinite particle limit, namely the Smoluchowski equations (2.1).This gives additional design flexibility in choosing the particle dynamics, and as we show in Theorems 1 and 2, all of the choices we make below are asymptotically correct.In fact, as we show in several examples in Section 6, using "fictitious" jump laws can actually approximate the infinite particle limit better than those choices which are "realistic".
Thinking of these transitions as taking place in the context of a continuous time Markov chain (CTMC), we would like to define the rates of the CTMC as and the event jumps by where e a i is the unit vector corresponding to one particle of size i and class a.More concretely, this means that if we choose a family of independent Poisson processes of rate one, N abc ijαβγ (t), then we have Notice that if α = 0 (resp.β = 0), we will remove one (a, i) (resp.(b, j)) particle, and if α = 1 (resp.β = 1), we will do nothing to the (a, i) (resp.(b, j)) population.Conversely, if γ = 1 we will create a particle of type (c, i + j) and if γ = 0 we will do nothing to the (c, i + j) population.As we said above, the choice (α, β, γ) = (0, 0, 1) corresponds to the "direct" coagulation method where we always remove the two smaller particles and always create the larger one.Every other possibility exists in the split scheme, even (in theory) (α, β, γ) = (1, 1, 0), which would have no effect on Q (of course, we would like this to happen infrequently during an actual simulation, and this is a question dealt with in Section 4.3 below).
Of course, once we specify T abc ij and p abc ij (α, β, γ), this specifies a particular numerical scheme.The question remains as to whether this corresponds in any sense to an approximation for (2.2).Let us define a computational volume V , and then define the computational particle concentration Y a i (t), the physical particle concentration X a i (t), and the binned Smoluchowski variable X i (t) as follows: where w a i is a "weighting function" which can be interpreted as follows: each computational particle of class a and size i corresponds to w a i physical particles of class a and size i.The variable X i (t) corresponds to the concentration of all physical particles of any class with mass i, and the fundamental constraint is that X i (t) should somehow simulate the dynamical system (2.2).Of course, X i (t) is a random quantity, and x i (t) is deterministic, so some care is needed here.
We can deduce from (3.6) that or ) which is just a time-discretization for a continuous-time Markov chain [5,8,34].If it turns out that the quantity V −1 λ(V Y ; . . . ) is well-defined, and the variance of Y (t + h) − Y (t) goes to zero fast enough as V → ∞, then (3.11) limits to a deterministic ordinary differential equation (ODE) for y = E[Y ], and one can ask whether or not this ODE corresponds to (2.2).
In fact, we show below in Theorem 1 that if we define and we scale T abc ij so that this quantity is well-defined in the limit V → ∞, then in this limit, Y (t) is close to the deterministic quantity y(t), where (We specify precisely what we mean by "close" below in Theorem 1.) We want this system to correspond to that given in (2.2), so if we define the ensemble averages then (2.2) corresponds to a deterministic ODE for x i (t), and then, ultimately, y a i (t).We want these equations to match, and thus we define: Definition 1 (Consistency).Specifying T abc ij and p abc ij (α, β, γ) determines (3.13).Similarly, (2.2) gives an ODE for x i (t) and thus gives evolution equations for y a i (t) (q.v.(4.3) below).If these choices give us the same ODE, we say they are consistent.
It turns out that to check whether or not a choice of method is consistent a priori is a rather complicated computation.We perform this computation in the general case in Section 4 below, but show that a simplifying Ansatz allows us to develop a much simpler set of conditions, given in (4.13).
One of the two main theoretical results of this paper is the following theorem.
Theorem 1. Assume that K ij has compact support, i.e., there exists an N > 0 such that K ij = 0 whenever i > N or j > N .Assume that T abc ij and p abc ij (α, β, γ) are chosen to satisfy the consistency conditions defined in (4.13).Let Q be a solution to the CTMC whose rates and transitions are given by (3.4,3.12) and define Y = Q/V and X as in (3.9).Let x solve (2.2).Then X converges to x in the following sense: for any T > 0, h > 0, there exist Remark 1.The constraint that K ij have compact support seems restrictive, especially since many realistic kernels are typically formally written as polynomials, examples including the linear or multiplicative kernels, K ij = i + j, K ij = ij.However, this restriction is not significant in practice, as we prove below.Since there is a conserved quantity (the total mass) in this system, it is not possible to ever have a particle larger than this mass in (2.2).Any stochastic simulation of this equation that gives a good approximation is unlikely to ever generate such a particle, and if it does not then truncation makes no difference.(We make this statement precise in Theorem 2 below.)4. Computations of consistency and optimality.In this section, we first compute explicitly the consistency conditions on T abc ij and p abc ij (α, β, γ), i.e., the conditions on these functions which make (3.13) as a function of y evolve the same as (2.2) as a function of x. 4.1.Consistency.Combining the definitions from (3.14, 3.7, 3.8, 3.9), we obtain Substituting this into (2.2),we obtain or We now expand (3.13) using (3.1, 3.4, 3.5, 3.12).
Writing out the specific coordinates of y(t), we have Making the change of dummy variables in the first sum of b → a, j → i, β → α, and switching these same variables in the third sum, we obtain If we symmetrize Thus for consistency we must choose d∈M,α,β∈{0,1} c∈M,α,γ∈{0,1} We call these the consistency conditions on T, p as per Definition 1.These are still quite complicated as written, so we will make a simplifying Ansatz in the next section.

4.2.
A simplifying Ansatz.We will now make a simplifying assumption on the probability distribution.As stated above, the probability distribution p abc ij (α, β, γ) has no constraints other than (3.3).
However, we now assume that the creation or destruction of each of the three particles involved in a coagulation are chosen independently, i.e., in a coagulation involving (a, i) and (b, j) coagulating into (c, i + j), we choose the three events "whether (a, i) is destroyed", "whether (b, j) is destroyed", and "whether (c, i + j) is created" independently.
More formally, we define numbers and and we assume that We see here that if we interpret p abc des,ij as the probability of destroying particle (α, i) during a coagulation event, and p abc cre,ij as the probability of creating particle (c, i + j) during this event, then (4.17) is simply asserting that these creations and destructions are done independently.
Writing out (4.13) under this Ansatz, we obtain It is more convenient to write these formulas if we switch (a, i) and (b, j) in the second equation and relabel j − i → j in the first equation, obtaining: Remark 1.In the case of just one type, where M = {1}, there is no sum on the RHS.This reduces to which is consistent with our earlier results in [7] (note the argument order for p death is reversed here).

4.3.
Optimality.Even under the independence Ansatz, there are still many choices we can make for T, p.In this section we discuss optimality, in the sense of having the least number of events generated by the method.Since they are probabilities, we require and this implies that we can ensure the consistency conditions if and only if c∈M where we use symmetry of T .For optimality we want to solve min subject to the above constraints and In general there will be many optimal solutions.One solution that is always optimal is to find and then set This expression is consistent for any choice of c * i+j .It is also optimal under the additional constraint that T abc ij = 0 for c = c * i+j .One nice practical possibility would be to require that c * is either a or b (as these are the two groups involved in the coagulation event already), and to pick c * to be the group with the lowest weight, thereby minimizing the statistical error.
Note that we can compute the creation and destruction probabilities by using the fact that T abc ij = 0 for c = c * , so giving If we define then we see that T can be written and Note again that the argument order convention for p des is reversed from [7].These optimal probabilities will be used in all numerical experiments in Section 6.
5. Theorem and proof.We will prove Theorem 1, which we restate here for completeness.
Theorem 1. Assume that K ij has compact support, i.e., there exists an N > 0 such that K ij = 0 for i > N or j > N .Assume that T abc ij and p abc ij (α, β, γ) are chosen to satisfy the consistency conditions defined in (4.13).Let Q be a solution to the CTMC whose rates and transitions are given by (3.4,3.12) and define Y = Q/V and X as in (3.9).Let x solve (2.2).Then X converges to x in the following sense: for any T > 0, h > 0, there exist C 1 , C 2 ∈ (0, ∞) such that Remark 2. Before we prove the theorem we give some motivation for the convergence in terms of generators.The generator of the stochastic process [34,42] is We will need to write this triple sum in many places below; whenever we write a sum without limits in the remainder of this section, we mean the triple sum which appears in (5.1).If we define and f (y) := f (V y), then we have a new operator Computing the first nonzero term in the Taylor series expansion on the right-hand side gives (5.7) Proof of Theorem 1.The main result that we will use in the proof is Kurtz's theorem [31], [42,Theorem 5.3].In fact, the statement of Theorem 1 is precisely the same as [42, Theorem 5.3] with the added assumption that the rescaled rates ρ abc ij (y; α, β, γ) be Lipschitz and uniformly bounded as functions of y a i (t), and that the vectors y a i (t) live in a finite-dimensional space.The functions ρ abc ij (y; α, β, γ) are locally Lipschitz, being polynomial.Moreover, under the assumption of compact support for K ij , they become globally Lipschitz.Similarly, since K ij = 0 if i, j ≥ N , we can also consider the finite-dimensional truncation of (2.2) without loss of generality, and note that this makes all of the sets of indices that we sum over above finite.
For the theorem to be true, we need a kernel of compact support.Assume now that we have a coagulation kernel K ij that has unbounded support.For any N > 0, we can consider the truncation of this kernel in the obvious way, by choosing The obvious question to ask is, if we choose two different truncations, , what is the probability that we get a different answer in the simulation scheme described in this paper?
First, define M t := i∈N1,a∈M iy a i (t), the total mass of the system at time t.Then notice that if |M t | < M , then y a i (t) = 0 for i > M .Moreover, sup so on the set y a i (t) < M , we have that y has finitely-many non-zero components and ρ is uniformly bounded and Lipschitz, as discussed in the proof above.Moreover, notice that M t is a martingale, i.e., that if Then we have the following theorem.
Theorem 2. Let N, N > 2M 0 and consider two truncated numerical methods, one with truncation parameter N and the other with truncation parameter N .More concretely, consider any numerical methods given by (4.13), but where we replace ; call these solutions X (N ) t and X (N ) t , respectively.Then for any set A and T > 0, there exist Proof.The first observation to be made here is that as long as M t < min(N, N ), then all of the statements in the theorem follow, because then the truncated kernels and thus the evolution (3.6) would be identical.So, in short, if we show that the probability that M t > 2M 0 is exponentially small, then the theorem follows.
Thus we consider either truncated process (we will choose X (N ) t for concreteness) and compute the probability of obtaining a realization with M t > 2M 0 .Denote the probability space of all outcomes of the stochastic process as Ω.Define Ω 2 := {ω ∈ Ω : M (X t ) > 2M 0 }, i.e., all outcomes in the probability space where the stochastic trajectory has twice the mass than it did at time zero.Notice that ).Then by Theorem 1, for any h, T > 0, there exist C 3 , C 4 > 0 such that and, by the same theorem, there exist C 5 , C 6 > 0 such that P(Ω 2 ) ≤ C 5 e −C6V .Then there are C 1 , C 2 > 0 with The theorem then follows by the triangle inequality.
6. Practical applications of method.Here we present a series of exhibitions of the practical effectiveness of the methods developed in this paper.In Section 6.1 we give an example of compositional multiscaling, in Section 6.2 we exhibit the usefulness of size-dependent weightings, and in Section 6.3 we give an example of both at once.
6.1.Independently weighted sub-populations.In this section we address the "compositional multiscaling problem", as introduced in Section 1.We use the formulation of Section 3 to weight particles in the rare sub-population differently to the rest of the particles.This allows more accurate statistics to be obtained for the sub-population while keeping the total particle number low and minimizing the computational cost.
Here we present this application in the context of two sub-populations: common sub-population 1 and rare sub-population 2. The extension to multiple subpopulations with varying frequency is straightforward.Take w a i to be independent of size i, so we have all particles in sub-population 1 weighted with w 1 and all particles in sub-population 2 weighted with w 2 .For coagulation events we will take the destination group c * to be whichever of the coagulating groups a or b has the lowest weight, as described in Section 4.3.That is, and we use the optimal expressions (4.37) and (4.38).
For a numerical test of coagulating independently weighted sub-populations, we use the Fuchs form of the Brownian coagulation kernel [26], where r i and r j are the radii of the coagulating particles, D p,i is the diffusion coefficient of particle i in air, δ i is the correction to account for particle i in the transition regime (see Eq. 15.34 in [26]), and vp,i is the thermal velocity of particle i in air.The exact expressions and constant values are specified in the file coag_kernel_brown.F90 and its dependencies in the PartMC code [45].Each sub-population begins at time t = 0 with a log-normal distribution with geometric mean diameter D g = 50 nm and geometric standard deviation log 10 σ g = 0.24.Sub-population 1 (blue in Figure 1) has an initial physical number concentration of N 1 = 10 5 cm −3 while sub-population 2 (red in Figure 1) has an initial physical number concentration of N 2 = 10 2 cm −3 .We assigned a density of 1000 kg m −3 and a molecular weight of 18 g mol −1 to the particles (pure water).The temperature was set to 288 K and the pressure was 10 5 Pa throughout the simulation.
Sub-population 1 particles have weight w 1 and sub-population 2 particles have weight w 2 .The resolution of each particle sub-population is related by the ratio Two special cases are r = 1, corresponding to equal weights for all particles (the equal weight case in Figures 1 and 2), and r = N 1 /N 2 = 10 3 , corresponding to equal numbers of computational particles used to initially represent each distribution on average (the equal number case in Figures 1 and 2).
The initial mean numbers of computational particles in each sub-population, denoted N p1 and N p2 , sum to the total number of computational particles, N p1 + N p2 = N p , and have weights corresponding to the physical number concentrations, so w 1 N p1 = N 1 and w 2 N p2 = N 2 .Using (6.3) gives ) In the case simulated here the physical concentration ratio is N 1 /N 2 = 10 3 , so when r = 1 we have N p1 ≈ N p and N p2 ≈ N p /10 3 , meaning that almost all computational particles are in the common sub-population 1 and almost none are in the rare sub-population 2. Correspondingly, when r = N 1 /N 2 = 10 3 we have N p1 = N p2 = N p /2, so the computational particles are evenly shared between the two sub-populations.
Figure 1 shows the number size distributions after t = 24 h of simulation for the equal weight and equal number cases, compared to very accurate finite-volume simulations.The error bars show the 95% spread for the particle simulations.The finite-volume simulations use a two-population scheme, with one set of finite-volume bins tracking each sub-population.Each combination of sub-populations coagulate using a high-resolution finite-volume scheme [3], with the coagulation products from 1-1 and 2-2 coagulations remaining in the respective sub-population, and 1-2 crosscoagulations going to the rare sub-population 2 following (6.1).
As we might expect, in Figure 1 the equal weights method substantially underresolves sub-population 2, while an equal-number weighting dramatically improves the resolution of sub-population 2 while only slightly degrading the sub-population 1 resolution.
The trade-off between sub-population 1 and sub-population 2 resolution as the ratio r is varied is quantified in Figure 2.Here ratios ranging from r = 1 to r = 10 6 are used for computational particle number varying from N p = 10 3 to N p = 10 5 .The special cases of equal weight (r = 1) and equal number (r = 10 3 ) are indicated by circles and squares, respectively.
For ratios between r = 1 and r = 10 3 we see a trade-off between sub-population 1 expected error and sub-population 2 expected error, with substantial improvements in sub-population 2 error accompanied by very small degradations in subpopulation 1 error.Equal-number weighting gives a reduction in sub-population 2 error by about 17 times from equal weights, at the expense of an increase of about 1.4 times in sub-population 1 error.
For ratios above r = 10 3 , the simulation has more computational particles resolving sub-population 2 than sub-population 1.This does not result in improved subpopulation 2 error, however, because the dynamics of sub-population 2 are highly dependent on sub-population 1 for coagulation partners.We thus see degrading performance in both sub-populations as the ratio increases and more computational particles are used for sub-population 2.
As shown in Theorem 1, we expect that the particle solutions will converge to the mean-field solutions generated by the very accurate finite-volume method as the number of particles increases, N p → ∞.The particle error versus N p is shown in Figure 3, and we see that the method is indeed converging as N p increases, as predicted, with rate N .Size distributions for two particle sub-populations using the scheme described in Section 6.1 with N p = 10 3 particles.Top: equal weightings for both distributions.Bottom: equal computational number for each distribution.The points are the mean of the particle process, while the error bars show 95% spread of realizations, not confidence intervals for the mean.The solid lines are very accurate finite-volume solutions.Using an equal-number rather than equal-weight weighting reduces sub-population 2 expected error by about 17 times, with an increase in sub-population 1 expected error of only about 1.4 times.
All particle simulations were implemented using the accelerated sampling scheme of Curtis et al. [4], so the computational cost is proportional to the number of coagulation events.In this subsection the number of coagulation events is almost independent of the weighting scheme chosen, because the two sub-populations have approximately equal size distributions and the coagulation kernel is only size-dependent, with no dependence on which sub-population contains the particle.Thus the computational cost is approximately constant for each value of N p in Figure 2, independent of the weight ratio r. , moving from upper-left to down and to the right on each line.circle points correspond to equal weights for both sub-populations (top panel in Figure 1), while the square points correspond to equal numbers of computational particles for each sub-population (bottom panel in Figure 1).Error bars are not shown as they are visually negligible.Convergence of the size distribution using the weighting scheme from Section 6.1 for different weight ratios r for the distribution of group a = 1 particles and group a = 2 particles.The baseline for computing the error is a very accurate finite-volume solution n fv .
6.2.Combined size-dependent weightings.Another practical problem that often arises in particle simulations concerns the relative scarcity of large particles formed by many coagulation events from small precursors, namely the "size multiscaling problem" discussed in Section 1.This occurs physically in raindrop formation, where approximately 10 6 small precursor droplets with diameter D ≈ 10 µm must coagulate to form a single, comparatively rare, rain droplet with diameter D ≈ 1 mm [29].
While it would in principle be possible to classify small and large particles into different sub-populations and use the differential weighting scheme discussed in Section 6.1, in practice it is often difficult to predict in advance just where the transition between these two sub-populations should be chosen, and how they should be independently weighted.In addition, it is undesirable to introduce an artificial cutoff in the accuracy of the simulation as a function of particle size, such as would arise from using discrete sub-populations to classify and weight particles.
An alternative approach, as developed in the Weighted Flow Algorithm (WFA) [7,28], weights particles with a smooth function of particle size.In DeVille et al. [7] it is shown that power-law weightings of the form w(D) ∝ D α are efficient for many physical simulation problems, where D is the particle diameter.This is equivalent to taking w a i ∝ i α/3 in our notation.While such exponential size-dependent weightings can be very efficient, they typically require a priori knowledge of the system evolution, so that the correct weighting function exponent can be selected.We can overcome this limitation and develop an a posteriori adaptive scheme by using the framework from Section 3.This extends the concept of composite weighting schemes [7,Section 6.7] to allow interactions between the component weighted particle sets.
For simplicity of exposition, we will only consider particles to be differentiated by size in this section, and we will discuss multiple sub-populations with size-dependent weightings in Section 6.3.We divide the particle population into two groups, so that each particle group can be doubled and halved independently to control total computational particle number, but the weight w a i will be independent of group number a.
To construct the weighting function, take precursor weight factors ν b i for b ∈ M and computational volumes V b .Taking a total computational volume V , we compute the weights w a i by the harmonic mean The harmonic mean is used as it corresponds to the minimum-variance estimator for a composite simulation [7, Section 6.7].In the numerical results for this section we use two precursor weight factors given by ν 1 i = 1 and ν 2 i = i −1 (corresponding to α = 0 and α = −3 diameter weighting exponents), and we use the optimal expressions (4.37) and (4.38).
For a numerical example we use a sedimentation kernel with an initial log-normal distribution of particles, having geometric mean diameter D g = 20 µm and geometric standard deviation log 10 σ g = 0.2.As in the example in Section 6.1, we assigned a density of 1000 kg m −3 and a molecular weight of 18 g mol −1 to the particles.The temperature was set to 288 K and the pressure was 10 5 Pa. Figure 4 shows the number (left column) and mass (right column) size distributions at time t = 0 and after t = 10 min of simulation for N p = 10 3 computational particles and three different weighting schemes.We see that the α = 0 (flat) weighting (top row) completely fails to capture the development of the second size peak around D = 5 mm and so mis-predicts both the number and size distributions.The α = −3 weighting (center row) is able to capture the second size peak and resolves the mass distribution fairly well, but fails to resolve the number distribution adequately.In contrast, the combined weighting scheme (bottom row) resolves both the number and size distributions well and accurately captures the evolution of the system.
The errors for varying exponent α and number of particles N p are shown in Figure 5. Changing the exponent from α = 0 (upper points) through to α = −3 (lower points) tends to increase number error while reducing mass error.This is because larger particles have lower weight with negative exponents, leading to greater resolution at larger sizes.The filled circles show the errors for the combined α = 0 and α = −3 weightings, and we see that this delivers equal-best mass error with better-than-best number error for a given number of computational particles N p .Compared to the simple α = 0 weight, the combined weighting has at least a factor of 2 improvement in number error and more than 20 times improvement in mass error.Compared to the simple α = −3 weight, the combined weighting has about the same mass error and about 5 times small number error.
The benefits of the combined weighting are less substantial when compared to the simple α = −1 or α = −2 weights.The main advantage of the combined weighting, however, is that it delivers excellent error performance without needing to be tuned to the particular simulation.Looking at the red solid number curves (finite-volume simulations) in Figure 4, we see that the peak-to-peak straight line fit has a slope of about α = −1.8 on the log-scale, which is why the α = −2 weight works well.For different initial conditions or simulation parameters, however, this weight will not perform as well.In contrast, the combined weighting scheme adapts to the simulation and so delivers excellent resolution without needing to know the solution behavior in advance.Theorem 1 implies that the particle solutions will converge to the mean-field solutions as the number of particles N p increases.The particle error versus N p is shown in Figure 6 for the number and mass distributions, which shows that the method is in fact converging as N p → ∞.For the flat weighting (α = 0) the convergent regime does not begin until N p ≥ 10 4 , as is expected from Figure 4, where we see very poor behavior for α = 0 with only N p = 10 3 particles.
The computational cost of the simulations shown in Figure 4 is again linear in the number of coagulations that occur, as described in Section 6.1.The inversemass (α = −3) and combined (α = 0 and α = −3) weightings result in simulations that are more expensive than the flat (α = 0) weighting, because there are more coagulation events being resolved.This is because the sedimentation kernel (6.7) is largest between particles of different sizes, so the simulations that better resolve the large particles will simulate more coagulation events between small and large particles.
6.3.Size-dependent sub-population weightings.The two practical problems addressed by the schemes of Sections 6.1 and 6.2, namely the composition multiscaling problem and the size multiscaling problem, frequently occur simultaneously in application scenarios [23,43].In such cases, highly efficient simulations can be obtained by combining the independent sub-population weightings of Section 6.1 with combined size-dependent weightings of Section 6.2.
We will not present a comprehensive convergence study of the combined algorithm, but rather illustrate its behavior in Figure 7.This shows the Brownian kernel test-case from Section 6.1 with four different simulations, and both the number and mass output distributions.The four cases (top to bottom) in Figure 7 are: (1) a naive unmodified simulation, with equal weights for all particles; (2) sub-population weights using the Section 6.1 scheme, but no size-dependence; (3) size-dependent weights using the combined scheme of Section 6.2, but no sub-population differentiation; and (4) using the schemes of Sections 6.1 and 6.2 simultaneously.
The observed numerical behavior corresponds to that described in detail in Sections 6.1 and 6.2.The simulations (2 and 4) with differentially-weighted subpopulations have dramatically improved error in sub-population 2, as explained in Section 6.1.The simulations (3 and 4) with combined size-dependent weighting have excellent behavior for both number and mass distributions, rather than just the number distribution, as explained in Section 6.2.In summary, simulation 4 combines the advantages of both the previously outlined methods.
7. Conclusions.This work describes the method development of particle simulations for the Smoulchowski equation where multiple particle sub-populations exist and possibly interact.Our particular target problem considered sub-populations that may exhibit large differences in their abundances.This is a common feature in cloud microphysics and atmospherically relevant aerosol applications and poses a challenging modeling problem.We conceptualized this challenge as the "size multiscaling problem" and the "compositional multiscaling problem".
We proved convergence for these methods under very weak assumptions and developed a family of consistency conditions that result in the correct mean during the simulation.Due to the many degrees of freedom in the parameters of the stochastic simulation, we showed that many choices gave consistency, so we also developed the conditions for the optimal consistent method.We stated and proved a family of theorems to guarantee that these methods are accurate in the limit of large population size.Finally, we presented numerical simulations illustrating the applicability and efficacy for two atmospheric relevant aerosol systems, one representative of the size multiscaling problem using the sedimentation kernel, and the other representative of the compositional multiscaling problem, using the Brownian kernel.Since in practice the two problems often occur simultaneously, we also showed that the methods can be combined and result in excellent behavior for both number and mass distributions.

1 Figure 1
Figure1.Size distributions for two particle sub-populations using the scheme described in Section 6.1 with N p = 10 3 particles.Top: equal weightings for both distributions.Bottom: equal computational number for each distribution.The points are the mean of the particle process, while the error bars show 95% spread of realizations, not confidence intervals for the mean.The solid lines are very accurate finite-volume solutions.Using an equal-number rather than equal-weight weighting reduces sub-population 2 expected error by about 17 times, with an increase in sub-population 1 expected error of only about 1.4 times.

3 N p = 10 4 NFigure 2 .
Figure2.Errors for two particle sub-populations for varying number of particles N p .The weight ratio r = w 1 /w 2 is 1, 10, 10 2 , 10 3 , 10 4 , 10 5 , and 10 6 , moving from upper-left to down and to the right on each line.circle points correspond to equal weights for both sub-populations (top panel in Figure1), while the square points correspond to equal numbers of computational particles for each sub-population (bottom panel in Figure1).Error bars are not shown as they are visually negligible.

a = 1 , r = 10 6 a = 2 , r = 10 6 a = 1 , r = 10 3 a = 2 , r = 10 3 a = 1 , r = 1 a = 2 , r = 1 Figure 3 .
Figure3.Convergence of the size distribution using the weighting scheme from Section 6.1 for different weight ratios r for the distribution of group a = 1 particles and group a = 2 particles.The baseline for computing the error is a very accurate finite-volume solution n fv .

fv 2 ]Figure 5 .
Figure 5. Number and mass errors for the sedimentation simulations described in Section 6.2 shown in Figure 4 for varying number of particles N p .The solid lines show size-weighted simulations with varying weight exponent α ∈ {0, −1, −2, −3} (ordered from top to bottom).The filled circles are combined flat-and-inverse-mass weighted simulations.

Figure 6 .
Figure 6.Convergence of the number and mass size distribution using the weighting scheme from Section 6.2 for different weight exponents α and for the combined weighting.The baseline for computing the error is a very accurate finite-volume solution with number distribution n fv and mass distribution m fv .The exponents α = −1 and α = −2 have similar behavior (not plotted).