Aerosol dynamics are described by the population balance equation (PBE). In principle, three typical methods (i.e., direct discretization, method of moments, and stochastic method) have been widely used to solve the PBE. Stochastic method is the most flexible among the three methods. However, stochastic method is computationally expensive. Recently, an operator splitting Monte Carlo (OSMC) method has been developed so as to improve the numerical efficiency while preserving the flexibility of the stochastic method. Within the OSMC, nucleation and surface growth are handled with deterministic means, while coagulation is simulated with a stochastic method (the Marcus‐Lushnikov stochastic process). The stochastic and deterministic treatments of various aerosol dynamic processes are synthesized under the framework of operator splitting. Here, the operator splitting errors of various schemes are analyzed rigorously, combined with concrete numerical examples. The analyses not only provide sound theoretical bases for selecting the most efficient operator splitting scheme for the usage of the OSMC, but also shed lights on how to adopt operator splitting in other PBE solving methods, i.e., direct discretization, method of moments, etc.
- Monte Carlo simulation
- aerosol dynamics
- operator splitting
Aerosols are dispersed systems of solid or liquid particles suspended in air. Aerosols are created during various natural processes and industrial operations , such as clouds, smoke particles from forest fires, sand and dust storms, volcanic dust, spores and seeds from plant life, fluidized catalysts, spray drying of fluids, atomization of liquid fuels, fly ash from stacks, soot in combustion flames, etc. Aerosol evolution is very complex, involving various physical and chemical processes. The most popular description of aerosol evolution is using the population balance equation (PBE) , which is also called the general dynamic equation in aerosol field . The PBE is a convection‐diffusion type equation with source terms taking into account the particle formation and depletion from various aerosol dynamics, mainly nucleation, surface growth (condensation), and coagulation.
Analytical solutions to the PBE are available for only a few specific cases [4–7]. Mostly, numerical means are resorted to solve the PBE. There are generally three typical methods for solving the PBE, i.e., direct discretization, method of moments, and stochastic methods. Direct discretization may be in the particle size space, such as discretization by sections  or finite difference/finite element  and may also be in the functional space, such as those in the collocation method , Galerkin method , and weighted residual method . Among the various discretizations, the sectional method pioneered by Gelbard et al.  is the most popular. There are quite a few sectional methods developed aiming to alleviate the nonpositiveness and nonconservation problems found in earlier works , and to obtain better efficiency and flexibility, such as the fixed/moving pivot method [14, 15] and cell average method [16, 17].
Other than solving the discretized PBE, a widely used alternative approach is the method of moments, which solves a group of moments equations derived from the PBE. These moments (generally several low order moments) provide important information on the particle size distribution (PSD) function, i.e., number density, volume fraction, polydispersity, etc. However, owing to the nonlinearity of the PBE, the governing equations for lower order moments generally contain higher order moments, which are not closed after cutting off at a certain order of moment. Various closure methods have been developed, such as unimodal log‐normal approximation , quadrature/direct quadrature [19, 20], interpolative closure , Taylor expansion , etc. Although the method of moments is computationally highly efficient, the complexity of the closure makes it difficult to accommodate complex physical models for aerosol dynamics with high flexibility. Another issue closely tied to method of moments is the realizability problem, when the moments sequence fails to find its corresponding distribution function [23, 24]. Among all the drawbacks, the inability to provide the PSD is the most severe, when the PSD is highly needed.
Stochastic simulation (or Monte Carlo simulation) does not solve the PBE, but to mimic the evolution of aerosol particles through a stochastic particle system. Based on the Marcus‐Lushnikov process modeling of coagulation , various stochastic methods have been developed to simulate the general aerosol evolution. It is very difficult to categorize precisely so many diversified methods. According to the evolution of time, there is time‐driven or event‐driven Monte Carlo [30, 31]. According to the different ways to weight the numerical particles, there are the mass flow algorithms [32, 33], differentially weighted algorithm , and the more general weighted flow algorithms . Also, there are various types of methods to improve the simulation efficiency, such as the majorant kernel method  and the leaping method .
Among the three typical methods, stochastic simulation is the most flexible for accommodating complex physical models for various aerosol dynamics, and is easy to implement without special numerical issues, which are often encountered in the other two methods, e.g., nonpositiveness and nonconservation problems in direct discretization, realizability problem in method of moments. However, stochastic simulation generally has convergence rate proportional to , where is the number of total numerical particles used. The simulation cost (CPU time and memory) is usually proportional to . Hence, stochastic simulation is not numerically efficient when pursuing a high precision solution. Recently, Zhou et al.  developed an operator splitting Monte Carlo (OSMC) method, which combines stochastic simulation for coagulation process and deterministic integration for nucleation and surface growth (condensation). Stochastic simulation for coagulation is much more efficient than directly solving the Smoluchowski integro‐differential equation (the coagulation model equation) with traditional discretization . On the other hand, deterministic integration for nucleation and surface growth processes is far more efficient than stochastic simulations. The stochastic simulation and the deterministic integration in OSMC are synthesized under the framework of operator splitting. The accuracy and efficiency of the OSMC have been demonstrated through various testing cases . Especially, the OSMC has been implemented with Message Passing Interface (MPI) to achieve very promising parallel simulation efficiency . Here, the impact of various operator splitting schemes on the accuracy and efficiency of the OSMC is addressed.
2.1. Governing equations
Aerosol particles are generally described statistically in a mesoscopic scale through defining the PSD, , and denotes the number density of aerosol particles at location at time with aerosol properties between and . Here, is usually called an internal variable, which can be volume, surface, constituent, etc. For simplicity, it is assumed to be the particle size (volume or diameter specified by the context). It is worth pointing out that the OSMC method can also deal with the cases with a vector type , i.e., combined properties to describe particles. In fact, the current OSMC code treats the diameter, surface, and volume of an aerosol particle independently, which can easily model particles of various morphologies, such as particles with fractal dimensions. The evolution of the PSD function satisfies the PBE (also called the general dynamic equation) 
where is the gas flow velocity, and is the diffusion coefficient. The source terms of Inuc, Gcond , and Ccoag on the right hand side denote nucleation, surface growth (condensation), and coagulation, respectively. The determination of the convection velocity requires solving the fluid dynamics control equations, and the coupling of the fluid dynamics with the aerosol dynamics is another important research field. Here, the gas velocity is assumed known, using some techniques (often computational fluid dynamics). It is worth pointing out that the coupling between the fluid dynamics and the aerosol dynamics is generally handled by the operator splitting method [11, 28–30], i.e., to solve the fluid dynamics and the aerosol dynamics one by one with the feedback from the other dynamics temporarily freezing, and to iterate between these two dynamics to obtain a converged solution.
Another simplification in this work is to neglect the diffusion term in the PBE. For aerosol particles, the mass diffusion due to the number concentration gradient () is usually very small compared to the air viscous diffusion, which can be described by the Schmidt number (i.e., the ratio of the viscous diffusion to the mass diffusion). For example, under standard atmospheric conditions, a spherical particle with a diameter of 10 nm has a Schmidt number equal to 290, and a 100 nm particle has Schmidt number . So it is legitimate to neglect the aerosol mass diffusion under very general conditions.
With the assumptions that the gas velocity is known and the mass diffusion term is negligible, the PBE can be simplified as
The convection term on the left hand side of the original PBE (1) has been implicitly solved with the introduction of the Lagrangian time
The transformation converts the Eulerian point of view in Eq. (1) to the Lagrangian point of view in Eq. (2), with the help of the assumed known velocity field . There are many practically interested conditions, where the flow field can be considered as steady and one‐dimensional, e.g., burner stabilized flames and counter flow diffusion flames. Then it is very easy to integrate Eq. (3) to build up a one‐to‐one mapping between the time and the spatial coordinate. Hence, the convection can be readily assimilated in the model Eq. (2). For general cases of aerosol evolution in three dimensional flows, it is possible to use the Lagrangian particle scheme to simulate the aerosol evolution along Lagrangian trajectories [31, 32]. So the model Eq. (2) and the subsequent OSMC method may have very broad applications. Later on, is simply denoted as , since the Lagrangian point of view is not explicitly referred any longer.
The aerosol dynamic processes on the right hand side describe the interactions among molecules (or cluster of molecules) and particles. Nucleation is the process that dozens or hundreds of molecules form a stable critical nucleus. Surface growth includes physical condensation and surface chemical reactions, which involves interactions between gas phase molecules and an aerosol particle. Surface chemical reactions will not be discussed here due to their vastness. Coagulation is the process that two particles coalesce to form a new one. More details about the aerosol dynamic processes are presented in Section 2.2.
2.2. Operator splitting Monte Carlo method
Operator splitting is an efficient and powerful method to solve such evolution equations  as Eq. (2). In an operator splitting method, the time evolution of is still handled by the classical methods, such as the forward Euler method or Runge‐Kutta method. Other than to integrate all aerosol dynamic processes together in one time step , operator splitting separates the integration into multiple fractional steps. There are various ways to construct the multiple fractional steps, such as or , where 1/2 denotes a half integration step. The flowchart (Figure 1) shows an example how the operator splitting scheme is carried out. Details of the various models for the aerosol dynamics are described in the following subsections, and the error analyses for various operator splitting schemes are given in Section 3.
Nucleation is the gas‐to‐particle transition to generate nascent nanoparticle. There are various nucleation models, largely based on the classical nucleation theory developed in the 1920–1940s . In the classical nucleation theory, there are two key parameters for the nucleation modeling, i.e., the nucleation rate (how many particles are nucleated per unit volume, unit time), and the critical size, which is the minimum size of nucleated particles with stable nucleus. According to the classical Becker‐Döring theory [3, 35], the nucleation rate and critical diameter are
Here, is the vapor pressure, is the mole fraction of the vapor, is the molecular mass of the nucleation vapor, is the particle density, is the temperature, is the vapor saturation ratio, is the Boltzmann constant, is the surface tension, and is the molecular effective volume of the nucleation vapor. Hence, the nucleation term in Eq. (2) can be expressed with the help of the Dirac delta function as
which means that only particles with the critical size are nucleated.
Within the OSMC, nucleation model is not limited to the classical nucleation theory. Any other nucleation model can be readily included, should only the nucleation rate and the critical size be provided.
2.2.2. Surface growth (condensation)
The functional form of the particle volume growth rate due to condensation depends on the Knudsen number , i.e., the ratio between the air mean free path and the particle radius . In the free molecular regime () and the continuum regime (), the condensational growth rates can be derived from the gas kinetic theory and the gradient diffusion model, respectively. In the transition regime (), interpolation formulas are usually used, such as the harmonic mean of the free molecular and continuum regimes. The expressions for volume growth rate and the condensation term are 
Here, is the equilibrium vapor pressure for particles of size , which is usually assumed to be the saturation pressure of the condensation vapor. For particles close to the critical size, the Kelvin correction is also required .
Within the OSMC, any other surface growth model can also be accommodated, should only the particle volume growth rate be provided.
2.2.3. Stochastic simulation of coagulation
Coagulation is the process that two particles coalesce/aggregate to form a new particle. The coagulation dynamics is modeled by the well‐known Smoluchowski's equation 
The coagulation kernel function describes the rate at which particles of size coagulate with particles of size . The first term on the right‐hand side denotes the formation rate of particles with volume due to coagulation between particles and . (Volume is conserved before and after coagulation.) The half factor () is to correct the double‐counting effect. The second term on the right‐hand side describes the depletion rate of particles with volume due to their coagulation with any other particles. Practically, the coagulation kernel function is very hard to find, except for simplified conditions. In the free molecule regime (), has the form 
In the continuum regime (), takes the form 
where is the dynamic viscosity of the gas‐phase. In the transition regime (), interpolation formula taking into account the former two regimes are generally used.
Although the Smoluchowski's equation was developed one century ago, and it has received broad attention in diverse fields, no analytical solutions are known for the kernels discussed above. Analytical solutions known so far are only limited to the following three simple kernels (or their linear combination) , i.e., (i) , (ii) , and (iii) .
Other than to solve the Smoluchowski's equation directly through various numerical means, the coagulation process is simulated by the Marcus‐Lushnikov stochastic process in the OSMC. The Marcus‐Lushnikov process is the natural stochastic analogue of the finite discrete Smoluchowski coagulation equation. The well‐known review paper  provides a wide‐ranging survey on the Marcus‐Lushnikov process and the Smoluchowski's equation. Here, we follow the work of Gillespie , and give an intuitive introduction to the stochastic simulation of coagulation.
In a virtual box of volume , there are “well‐mixed” particles. Any two particles may coalesce/aggregate to form a new one at a random time. The propensity of coalescence between two particles and is determined by the coagulation kernel . Let denote the probability that the next coalescence will occur in the time interval between and (), then it is found 
In essence, the simulation is to randomly pick up two particles and to coalesce (forming a new particle) at the random time according to the joint probability density function . Repeat the former process as time advances. During the simulation, the number of particles drops, and the particle number density drops also.
The key step in the simulation is to generate a random point according to the density function . The random coagulation time satisfies a Poisson process (the exponential distribution), and it can be generated as
where , and is a random number picked from the uniform distribution in the unit interval. Eq. (13) gives the classical way to generate a random variate from the exponential distribution with parameter .
After selecting the time at which two particles coagulate, the next step is to randomly select two particles according to the probability function . Loosely speaking, is the total coagulation probability between any two particles, and is the fraction of the coagulation probability over . Gillespie  introduced two methods to choose and , the so‐called full conditioning and partial conditioning methods. In the full conditioning method, the selection of and is accomplished in two steps. During the first step, calculate the partial sum , , then combine the partial sums in sequential to form a range with the index from 1 to denoting the part for each partial sum, and then it is straightforward to build up a linear mapping between the range and the unit interval . Draw a random number from the uniform distribution, and check the random number landing on which corresponding part in the range, and then select particle by the part index. During the second step, combine in sequential to form another range, and build up the linear mapping between the new range with the unit interval, then draw another random number to determine the particle in a similar way as in the first step. In the partial conditioning method, it is actually an acceptance‐rejection method to generate a bivariate random number from the function . In a general acceptance‐rejection algorithm to generate a vector random variable with density , it needs to choose a density function and a constant such that . The algorithm proceeds as follows :
Generate from density .
Generate from the uniform distribution on the interval .
If , then set ; else, return to 1.
The density function can be chosen freely. However, the optimal choice is to make the interval in 2 as small as possible, so that there is a better chance that so as to accept as the required without too often rejections and repetitions of the whole process. In the acceptance‐rejection method  to generate and according to , is chosen as a uniform function, and is an upper bound of all (preferably the smallest upper bound). It is argued  that when the values of are relatively comparable with each other (ideally constant) the particle conditioning method should be more efficient than the full conditioning method.
Within the OSMC, both the full and particle conditioning methods have been implemented to choose from easily.
2.3. Full algorithm for the OSMC
Figure 1 shows the flowchart of OSMC including all aerosol dynamics. Only the first‐order operator splitting is sketched in the figure.
At time , the simulator is initialized with given parameters, i.e., number of Monte Carlo (MC) particles , simulator size , and the initial PSD. The simulator size is determined as , where is the initial number density. If simulation starts from an empty case, i.e., , then the simulator is initialized with particles of uniform size (the real value of the size is immaterial to the simulation results), and the simulator size is set to a huge value, say , which renders a tiny particle number density to approximate the condition . If simulation starts from a case with a specific PSD, then the size of the simulation particles is assigned by a stochastic process to satisfy the initial PSD. There are various convenient ways to generate random number to satisfy a given distribution , although they are usually not necessary in Monte Carlo simulation of aerosol dynamics, since mostly the simulation starts from an empty case.
In the particles nucleation step, the simulator size is adjusted to reflect the change of particle number density due to nucleation. Nascent nucleated particles are added to the pool of MC particles. If the total number of MC particles exceeds the maximum allowable threshold value (which is set beforehand to balance the stochastic error and numerical cost), then down‐sampling is performed, i.e., exceeding MC particles are randomly removed from the pool to satisfy the number constraint. And then every MC particle undergoes the surface growth process according to its growth rate (usually depends on its size).
The coagulation simulation process (those steps grouped in the dashed bounding box in Figure 1) shows how to implement the stochastic algorithm of Gillespie . Updating coagulation kernel is to calculate the pair coagulation rate . The random coagulation time is generated according to Eq. (13). The comparison statement is to judge whether the time for two particles to coagulate is still admissible within the discrete time step . Here is the accumulated coagulation time in the coagulation step, which is initialized to zero before a coagulation step begins. If the comparison statement is true, two particles are selected to perform coagulation, and the number of MC particles decreases by one. The subsequent up‐sampling step is to keep the number of MC particles above a given threshold to avoid severe stochastic error.
3. Error analyses for operator splitting
Aerosol dynamics usually contain nucleation, surface growth, and coagulation processes. These processes have very different effects on the evolution of the total number density and volume fraction, which are the two most important statistics of aerosol particles. For the number density, it increases due to nucleation while decreases due to coagulation, and it does not change in the surface growth process (mainly condensation). For the volume fraction, it is conserved in the coagulation process. The volume fraction is generally dominated by surface growth process, and nucleation has limited contribution to the volume fraction directly, because nucleated particles are very small in size.
Nucleation is generally a function of the temperature and the nucleated vapor concentration, and it does not depend on the aerosol status (number density, particle size, etc). However, both condensation and coagulation are closely related to the number density, which is largely determined by the nucleation process. On the other hand, condensation and coagulation are usually loosely coupled through particle size and number density. Therefore, it is convenient and meaningful to consider two isolated cases: (i) nucleation + coagulation, which determine the number density, and (ii) nucleation + condensation, which determine the volume fraction.
3.1. Case of nucleation + coagulation
Before discussing the case of nucleation and coagulation, we introduce an integral transformation of the Smoluchowski's equation (9), which is very useful to understand the evolution of the moments of the PSD. Suppose is an arbitrary function of , and let
Specifically, when , is the order moment, and Eq. (15) transforms to the dynamic equation for the moment
One interesting but somewhat overlooked conclusion  is that the right hand side of Eq. (16) is negative if , and it is positive if , which mean that will grow with time for , and it will decrease with time for . Especially, is constant with time. Actually, is the volume fraction, which is conserved during the pure coagulation process. is the total number density, which does decay due to coagulation. It is worth pointing out that although fractional moments (with being a fraction number) do not possess obvious physical meaning, they are very useful in the method of moments with the interpolative closure  or with the Taylor expansion closure .
In the case of constant rate nucleation and constant kernel coagulation, the number density () satisfies the following simple model equation
One the right‐hand side, the first term models the contribution to the number density from nucleation, and represents the nucleation rate (number of particles nucleated per unit volume per unit time). The second term models the contribution from coagulation, and , where is a constant coagulation kernel. In fact, in the moment dynamic Eq. (16) if the coagulation kernel is the constant and let (to obtain the number density equation); we immediately have on the right‐hand side. The notation (be aware of the negative sign) is found to be able to greatly simplify the following analyses.
The assumption of a constant nucleation rate and a constant coagulation kernel does not severely limit the applicability of the following discussion on the operator splitting errors for the more general case of nucleation and coagulation. The analyses focus on the local error in a generally small single time step. Hence, even a time‐dependent nucleation rate can be readily assumed constant within the single time step. On the other hand, the coagulation kernel is a bivariate function of colliding particles sizes. Mathematically, it is possible to find a that satisfies
Of course, the right depends not only on the specific function, but also on the distribution . Here, we assume that a good estimate of is available and will not delve into this issue further.
Given the initial condition , the model Eq. (17) has analytical solution
Since, the solution Eq. (19) involves imaginary numbers. Using their relation to the complex exponential function, it is easy to show that Eqs. (19) and (20) are actually equivalent. Although the solution Eq. (20) (with replaced by ) is more natural, the form in Eq. (19) proves more handy in the following analyses. Here, we expand Eq. (19) at as Taylor series for later usage
On the other hand, we can try to use the operator splitting method to solve the same Eq. (17), and compare the obtained solutions with the true solution (19) to quantify the operator splitting errors. There are various operator splitting schemes. Here, we will focus on two first‐order Lie schemes and two second‐order Strang schemes, i.e., (i) , (ii) , (iii) , and (iv) .
In the operator splitting method, the following two joint equations are solved
(i) Lie scheme (Denoted as Lie1)
Eq. (22a) has solution:
Hence, . Take as the initial condition and solve Eq. (22b), then the final solution is found as
Let in the above equation, and expand it into Taylor series
The operator splitting error can be found by comparing the two Taylor expansion series (25) and (21)
(ii) Lie scheme (Denoted as Lie2)
and its Taylor expansion is
Hence, the operator splitting error is
(iii) Strang scheme (Denoted as Strang1)
Here, the solution to the joint Eqs. (22) is carried out in three steps. First, Eq. (22a) is solved with the initial condition to get the solution . Second, is used as the initial condition to solve Eq. (22b) to get the solution . Finally, is used as the initial condition at to solve Eq. (22a) again. The final solution at is found as
and its Taylor expansion is
Hence, the operator splitting error is
(iv) Strang scheme (Denoted as Strang2)
Similar to case (iii), the final solution is found as
and its Taylor expansion is
Hence, the operator splitting error is
To summarize, the two Lie schemes have the same magnitude in splitting error, both are of second order in within a time step. So the global integration error during time is of order , since the error will accumulate in steps. The Lie1 scheme underestimates the number density, while the Lie2 scheme overestimates that. This is because coagulation rate closely depends on the number density, the higher the number density, the faster the coagulation rate. Coagulation rate would be overpredicted if nucleation is integrated before coagulation, hence the number density would be underpredicted as in Lie1, and vice versa.
The splitting errors in the two Strang schemes are of third order in within a time step, which renders them a second‐order scheme globally. The Strang1 scheme is found to overestimate the number density, while the two Strang2 scheme to underestimate that. The magnitude of the absolute errors in these two schemes can be compared through investigating the coefficients of in Eqs. (32) and (35)
If , i.e., coagulation rate is smaller than nucleation rate, the absolute splitting error in Eq. (32) is smaller than that in Eq. (35). When , the error in Eq. (33) approaches a half of that in Eq. (36). On the other hand, when , the error in Eq. (32) nearly doubles that in Eq. (35). Hence, in the Strang schemes it is preferable to split the higher rate process into symmetric substeps to reduce the splitting error.
Only the most common first and second‐order operator splitting schemes are discussed so far. They are systematical ways to construct even higher order schemes. Especially, a third‐order scheme  can be obtained through the combination of the two Lie and two Strang schemes. Here, we are not interested to discuss the even higher order schemes other than the Lie and Strang schemes, since higher order schemes are complicate to implement and may not be worthy for the solving the aerosol dynamics. One interesting conclusion  is that the scheme is of third order, which can be used to check the accuracy of the above derivations. Keep the third‐order term in the error functions (26), (29), (32), and (29), then calculate the error in the third‐order scheme constructed as mentioned shortly before, it is found that only terms involving and higher orders survive. This is an indirect proof that all the derivations of the error functions are correct.
Following, we will briefly discuss the splitting errors for higher order moments other than . Since the total volume is conserved during coagulation, nucleation, and coagulation are independent of each other with respect to (volume fraction). Hence, there are no splitting errors for predicting for , the model equation for nucleation and coagulation can be easily derived from Eqs. (6) and (16) as
and is determined by nucleation only, which can be integrated as
Here, is the critical volume of nucleated particle, and is the initial value at . The solution to Eq. (37) is
Similar to the above error analyses, the four splitting errors for solving Eq. (37) can be readily obtained as
The splitting errors for the two second‐order Strang schemes are independent of the initial conditions.
3.2. Case of nucleation + condensation
Here, a case with constant nucleation rate and constant condensation growth in diameter is considered. A constant growth rate in diameter is not arbitrary at all. Eq. (7) gives the volume growth rate due to condensation as a piecewise function. The particle diameter growth rate function can be obtained easily through the sphere volume formula . The detailed function can also be found in reference . In fact, the particle diameter growth rate in the free molecular regime is independent of its size, which can be seen as a constant determined by environmental parameters (vapor pressure and temperature). In the following discussion, the particle diameter growth rate is assumed as . In this subsection, the moments are defined with respect to the particle diameter other than the volume . No attempt is sought to distinguish the symbols for moments defined differently, with the hope to reduce new symbols and to bring no confusion.
The dynamic equations for moments due to pure nucleation (to generate particles with diameter at rate ) can be easily obtained with the help of the Delta function (6), which are
The dynamic equations for moments due to condensation (at a constant diameter growth rate ) are
Hence, the dynamic moments equations for nucleation and condensation are
Assume the initial condition at is , then Eq. (43) has solution
(i) Lie scheme (Denoted as Lie1)
(ii) Lie scheme (Denoted as Lie2)
(iii) Strang scheme (Denoted as Strang1)
(iv) Strang scheme (Denoted as Strang2)
To summarize, there is no splitting error for the number density (), which is solely determined by nucleation. For to , all Lie schemes are of first order (globally), and the leading errors are of opposite sign between Lie1 and Lie2. For Strang schemes, the splitting error for is also zero, and the errors for and are of second order. In Lie1 and Strang1 scheme, the splitting solutions overpredict the high order moments, since nucleation in the prestep produces more than real particles to grow in condensation. On the other hand, Lie2 and Strang2 underpredict the high order moments. Since the splitting errors are in closed form of the time step , it is possible to construct higher order splitting scheme with zero splitting error for all moments.
Furthermore, the leading order (with respect to ) splitting errors for any two subsequent moments differ by a factor of , which renders the relative errors (splitting errors divided by the corresponding moments values) comparable among various moments, since two subsequent moments also differ by a factor of the particles average diameter according to the moment definition. The splitting errors are independent of the initial conditions, i.e., , , , and , which seems strange at first glance. In fact, the splitting errors are anchored to the moment's evolution. Take the splitting error in the Lie1 scheme as example, the error contains the factor , while the leading order in the analytical solution for in Eq. (58) also contains exactly the same factor . Then only the time step matters when calculating the relative error. Generally speaking, the relative errors due to the operator splitting are small, almost independent of the nucleation rate and the condensational growth rate.
4. Numerical examples
4.1. Case of nucleation + coagulation
Corresponding to the error analyses in Section 3.1, the case of constant nucleation and constant coagulation are simulated with the OSMC. The coagulation kernel and nucleation rate are set as unit , and . Initially particles are mono dispersed with unit volume size, and the number density is . Nucleated particles also have unit volume size.
Figure 2 shows the number density () simulated by the four schemes, compared against the analytical solution (19). All the simulations use fixed time step , with Monte Carlo particles. The OSMC results are averaged over dozens of independent simulations to eliminate the randomness. Lie1 scheme is found to systematically underpredict the number density, while Lie2 to overpredict. Strang1 gives result negligibly higher than the analytical solution. Strang2 slightly underpredicts the true solution. All these findings agree extremely well with the analyses in Section 3.1. It is also found that Strang1 shows better precision than Strang2, that is because nucleation rate is higher than the coagulation rate (), and the higher rate process nucleation is split in Strang1, as explained in Section 3.1 also.
Figure 3 compares among the four splitting schemes and the analytical solution (39). Lie1 is found to overestimate a little bit, while Lie2 to underestimate. The second‐order Strang schemes give result almost indistinguishable from the analytical solution. Although second‐order Strang schemes give better results than the Lie schemes, the discrepancy is not so significant as that for in Figure 2.
Both the error analyses and the numerical examples above show that the second‐order Strang schemes have much higher precision than the first‐order Lie schemes. Within the OSMC, coagulation is the most computationally intensive part, and nucleation only makes a tiny fraction of the total time cost. The Strang1 scheme has much higher precision than the Lie schemes, and only needs to integrate nucleation twice often as in the Lie schemes, where the computational overhead is negligible. Strang2 needs nearly double computational cost as Strang1, and Strang2 outperforms Strang1 in precision only if . The best performance of Strang2 over Strang1 is to nearly double the precision when . Hence, Strang1 is the optimal splitting scheme taking both the numerical precision and cost into account.
Another general conclusion is that only needs to be concerned when considering the operator splitting errors. Operator splitting generally produces less errors for higher order moments. The splitting errors for the Strang schemes contain the factor of the critical volume of nucleated particles, which is generally very small. Taking a more practical numerical example of dibutyl phthalate (DBP) aerosol nucleation and coagulation in a turbulent mixing layer , the critical volume is of order , which renders very small splitting errors for high order moments. However, the splitting error for may be still significant, and needs proper splitting scheme combined with small enough time step to satisfy the given error criteria.
4.2. Case of nucleation + condensation
Corresponding to the error analysis in Section 3.2, the case of constant nucleation and constant condensational growth in diameter is simulated with the OSMC. Both the nucleation rate and growth rate are set to unit and . Initially particles are monodispersed with unit volume size, and the number density is . Nucleated particles also have unit volume size. The time step is .
Figures 4 and 5 show the evolution of moments and , respectively, obtained by the four splitting schemes and the analytical solutions (44d). While directly represents the particle volume fraction (lacking a factor of ), do not have simple physical meaning, which may be used to define an average diameter through . For both and , the characteristics of splitting errors are quite similar. Lie1 is found to overpredict, while Lie2 to underpredict. The splitting errors for the two Strang schemes are very small, almost indistinguishable from the analytical solutions. According to the error analyses (49c) and (51c), the splitting errors for Strang schemes are zero. Close scrutiny shows in Figure 4 that the numerical error is not actually zero, which actually arises from the ordinary differential equation (ODE) integration error, since a relatively large time step is used. For a more appropriate smaller time step, both the ODE integration error and splitting error are much smaller. Comparing Figures 4 and 5, the relatively errors are comparable for both and here. All the findings agree very well with the analyses in Section 3.2.
The operator splitting Monte Carlo (OSMC) method has been developed recently  for solving the population balance equation for aerosol dynamics. Within the OSMC, nucleation and surface growth are handled with deterministic means, while coagulation is simulated with a stochastic method (the Marcus‐Lushnikov stochastic process). The deterministic and stochastic approaches in the algorithm are synthesized under the framework of operator splitting. The ultimate goal of the OSMC is to greatly improve the numerical efficiency and to preserve the extraordinary flexibility and applicability of a stochastic method for aerosol dynamics simulation.
The OSMC has been validated through quite a few testing cases , and has also been used to simulate aerosol evolution under various conditions [32, 39]. However, the operator splitting errors in the OSMC have not been systematically investigated. Here, focused on two representative cases, i.e., constant nucleation and coagulation, and constant nucleation and condensation, the splitting errors for four splitting schemes (two Lie schemes and two Strang schemes) are analyzed rigorously, combined with concrete numerical examples.
The number density, one of the most important statistics for aerosol particles, is found to be underestimated for the Lie1 scheme (), while overestimated for the Lie2 scheme (). This is because coagulation rate would be overpredicted if nucleation is integrated before coagulation, hence the number density would be underpredicted as in Lie1, and vice versa. The second‐order Strang schemes exhibit much better precision that the first‐order Lie schemes. If coagulation rate is higher than the nucleation rate, Strang2 () has better precision than Strang1 (). Otherwise, Strang1 has better precision. However, for the OSMC, Strang1 is always the most preferable splitting scheme after taking the numerical cost and precision into account.
The operator splitting errors for the case of nucleation and condensation show that splitting errors for different moments (except for ) are comparable with each other. Splitting between nucleation and condensation has marginal impact on accurately predicting the moments. This is in great contrast to the prediction of , where splitting between nucleation and coagulation has significant impact.
To construct a good operator splitting scheme for the general case of nucleation, coagulation, and condensation the first priority is to consider the splitting error for . When a good scheme for splitting between nucleation and coagulation is constructed, splitting among other aerosol dynamic processes can be constructed at more freedom. The overall splitting errors are pretty much determined by the splitting error for .
The analyses not only provide sound theoretical bases for selecting the most efficient operator splitting scheme for the usage of the OSMC, but also shed lights on how to adopt operator splitting in other PBE solving methods, i.e., direct discretization, method of moments etc.
The financial support from the National Natural Science Foundation of China (No. 11402179 and 11572274) is greatly acknowledged. KZ would like to thank Dr. F. Bisetti at KAUST, Saudi Arabia, for his insightful discussion and valuable help. The financial support from the General Research Fund, Research Grants Council of the Hong Kong SAR, China (No. PolyU 152663/16E) is also greatly appreciated for allowing the authors’ further development and extension of this OSCM foundation method.