**3. Error analyses for operator splitting**

constant) the particle conditioning method should be more efficient than the full condition‐

Within the OSMC, both the full and particle conditioning methods have been implemented to

**Figure 1** shows the flowchart of OSMC including all aerosol dynamics. Only the first‐order

At time = 0, 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 <sup>=</sup> /0, where 0 is the initial number density. If simulation starts from an empty case, i.e., 0 = 0, 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 = 1010, which renders a tiny particle number density to approximate the condition 0 = 0. 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 [37], although they are usually not necessary in Monte Carlo simulation of aerosol dynamics, since mostly the simulation starts

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

The coagulation simulation process (those steps grouped in the dashed bounding box in **Figure 1**) shows how to implement the stochastic algorithm of Gillespie [25]. Updating

 = + 1,…,). The random coagulation time is generated according to Eq. (13). The com‐ parison statement sum <sup>+</sup> <sup>&</sup>lt; is to judge whether the time for two particles to coagulate is still admissible within the discrete time step . Here sum is the accumulated coagula‐ tion 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 coagula‐ tion, 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

,

)/, ( = 1,…, 1,

process according to its growth rate (usually depends on its size).

coagulation kernel is to calculate the pair coagulation rate <sup>=</sup> (

ing method.

choose from easily.

10 Aerosols - Science and Case Studies

from an empty case.

error.

**2.3. Full algorithm for the OSMC**

operator splitting is sketched in the figure.

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 transfor‐ mation 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

$$M(t) = \int\_0^\wp F(v) n(v, t) \, d\upsilon\_\prime \tag{14}$$

then from the Smoluchowski's equation (9), it can be derived [1]

$$\frac{\mathrm{d}M(t)}{\mathrm{d}t} = \frac{1}{2} \int\_{0}^{v} \int\_{0}^{v} [F(v+u) - F(v) - F(u)] \beta(v,u) n(v,t) n(u,t) \mathrm{d}u \mathrm{d}v. \tag{15}$$

Specifically, when () = , ( 0), () is the th order moment, and Eq. (15) transforms to the dynamic equation for the moment ()

$$\frac{\mathrm{d}M\_k(t)}{\mathrm{d}t} = \frac{1}{2} \int\_0^\upsilon \int\_0^\upsilon [(\upsilon + \mathsf{u})^k - \upsilon^k - \mathsf{u}^k] \beta(\upsilon, \mathsf{u}) n(\upsilon, t) n(\mathsf{u}, t) d\mathsf{u} d\mathsf{v}.\tag{16}$$

One interesting but somewhat overlooked conclusion [1] is that the right hand side of Eq. (16) is negative if < 1, and it is positive if > 1, which mean that will grow with time for  > 1, and it will decrease with time for < 1. Especially, 1() is constant with time. Actually, 1 is the volume fraction, which is conserved during the pure coagulation process. 0 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 [21] or with the Taylor expansion closure [22].

In the case of constant rate nucleation and constant kernel coagulation, the number density (0) satisfies the following simple model equation

$$\frac{\mathbf{d}M\_0}{\mathbf{d}t} = \mathbf{J} + \mathbf{K}M\_0^2. \tag{17}$$

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 <sup>=</sup> <sup>−</sup> <sup>1</sup> 20, where 0 is a constant coagulation kernel. In fact, in the moment dynamic Eq. (16) if the coagulation kernel is the constant 0 and let = 0 (to obtain the number density equation); we immediately have <sup>−</sup> <sup>1</sup> <sup>2</sup>00 <sup>2</sup> on the right‐hand side. The notation <sup>=</sup> <sup>−</sup> <sup>1</sup> <sup>2</sup>0 (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 0 that satisfies

$$
\beta\_0 \beta\_0 M\_0^2 = \beta\_0 \int\_0^\alpha \int\_0^\upsilon \mathfrak{n}(\upsilon, t) \mathfrak{n}(\mu, t) \mathrm{d}\mathfrak{n} d\upsilon = \int\_0^\alpha \int\_0^\upsilon \beta(\upsilon, \mu) \mathfrak{n}(\upsilon, t) \mathfrak{n}(\mu, t) \mathrm{d}\mathfrak{n} d\upsilon. \tag{18}
$$

Of course, the right 0 depends not only on the specific function, but also on the distribution . Here, we assume that a good estimate of 0 is available and will not delve into this issue further.

Given the initial condition 0( = 0) = 0 \*, the model Eq. (17) has analytical solution

$$M\_{\rm o}(t) = \frac{\sqrt{|K|}}{K} \tan \left[ t \sqrt{|K|} + \arctan \left( \frac{M\_{\rm o}^{\star} K}{\sqrt{|K|}} \right) \right]. \tag{19}$$

In fact, a solution of the same Eq. (17) is found as [7]

 > 1, and it will decrease with time for < 1. Especially, 1() is constant with time. Actually, 1 is the volume fraction, which is conserved during the pure coagulation process. 0 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 [21] or

In the case of constant rate nucleation and constant kernel coagulation, the number density

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 <sup>=</sup> <sup>−</sup> <sup>1</sup>

where 0 is a constant coagulation kernel. In fact, in the moment dynamic Eq. (16) if the coagulation kernel is the constant 0 and let = 0 (to obtain the number density equation); we

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 0

> *M nvtnut u v vunvtnut u v* <sup>2</sup> 00 0 0 0 0 0 = ( , ) ( , )d d = ( , ) ( , ) ( , )d d . ¥ ¥ ¥ ¥

> > *JK M K M t t JK <sup>K</sup> JK*

<sup>0</sup> ( ) = tan arctan .

Of course, the right 0 depends not only on the specific function, but also on the distribution . Here, we assume that a good estimate of 0 is available and will not delve into this issue

b

é ù æ ö ê ú <sup>+</sup> ç ÷ ç ÷ ë û è ø

ò ò ò ò

negative sign) is found to be able to greatly simplify the following analyses.

<sup>2</sup> on the right‐hand side. The notation <sup>=</sup> <sup>−</sup> <sup>1</sup>

0

+ (17)

(18)

\*, the model Eq. (17) has analytical solution

\* 0 20,

(19)

<sup>2</sup>0 (be aware of the

*<sup>M</sup> J KM <sup>t</sup>* 0 2

<sup>d</sup> = . <sup>d</sup>

with the Taylor expansion closure [22].

12 Aerosols - Science and Case Studies

immediately have <sup>−</sup> <sup>1</sup>

that satisfies

further.

b

 b

Given the initial condition 0( = 0) = 0

(0) satisfies the following simple model equation

<sup>2</sup>00

$$M\_0(t) = \frac{\sqrt{-fK}}{K} \tanh\left[t\sqrt{-fK} + \text{arctanh}\left(\frac{M\_0^\*K}{\sqrt{-fK}}\right)\right].\tag{20}$$

Since, <sup>=</sup> <sup>−</sup> <sup>1</sup> <sup>2</sup>0 < 0 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 <sup>−</sup> <sup>1</sup> 20) is more natural, the form in Eq. (19) proves more handy in the following analyses. Here, we expand Eq. (19) at <sup>=</sup> as Taylor series for later usage

$$M\_0(\boldsymbol{\delta}\_t) = M\_0^\ast + (\boldsymbol{\mathcal{K}} \boldsymbol{M}\_0^{\ast 2} + \boldsymbol{I})\boldsymbol{\delta}\_t + \boldsymbol{\mathcal{K}}(\boldsymbol{\mathcal{K}} \boldsymbol{M}\_0^{\ast 2} + \boldsymbol{I})\boldsymbol{M}\_0^\ast \boldsymbol{\delta}\_t^2 + \frac{1}{3}\boldsymbol{\mathcal{K}}(\boldsymbol{\mathcal{G}} \boldsymbol{\mathcal{K}}^2 \boldsymbol{M}\_0^{\ast 4} + \boldsymbol{\mathcal{G}}\boldsymbol{\mathcal{K}} \boldsymbol{M}\_0^{\ast 2} + \boldsymbol{I}^2)\boldsymbol{\delta}\_t^3 + \mathcal{O}(\boldsymbol{\delta}\_t^4). \tag{21}$$

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) nuc coag, (ii) coag nuc, (iii) 1 2 nuc coag 1 2 nuc, and (iv) <sup>1</sup> <sup>2</sup>coag 1 2coag.

nuc In the operator splitting method, the following two joint equations are solved

$$\frac{\text{d}M\_{\text{o}}^{(1)}}{\text{d}t} = \text{J},\tag{22a}$$

$$\frac{\text{d}M\_0^{(2)}}{\text{d}t} = \text{K}M\_0^2. \tag{22b}$$

(i) Lie scheme nuc coag (Denoted as Lie1)

To seek the solution at <sup>=</sup> , Eq. (22a) is solved first with the given initial condition at = 0, then Eq. (22b) is solved with the just obtained solution at <sup>=</sup> as the initial condition.

Eq. (22a) has solution:

$$M\_0^{(1)}(t) = ft + M\_0^"\,. \tag{23}$$

Hence, 0 (1)( ) = <sup>+</sup> 0 \*. Take 0 (1)( ) as the initial condition and solve Eq. (22b), then the final solution is found as

$$M\_{\rm o}^{(\rm OP)}(t) = -\frac{J\delta\_t + M\_{\rm o}^{"\prime}}{J\mathbf{K}\delta\_t t + K M\_{\rm o}^{"\prime}t - 1}.\tag{24}$$

Let <sup>=</sup> in the above equation, and expand it into Taylor series

$$M\_0^{(\text{OP})} (\delta\_t) = M\_0^"+ (\text{KM}\_0^2 + I)\delta\_t + \text{K}(\text{KM}\_0^{"\*2} + 2I)M\_0^"\delta\_t^2 + \mathcal{O}(\delta\_t^3). \tag{25}$$

The operator splitting error can be found by comparing the two Taylor expansion series (25) and (21)

$$E r\_{0i\text{j}} = M\_0^{\text{(CP)}}(\mathcal{S}\_i) - M\_0(\mathcal{S}\_i) = M\_0^"\!/\text{K} \mathcal{S}\_i^2 + \mathcal{O}(\mathcal{S}\_i^3) = -\frac{1}{2} M\_0^"\!/\mathcal{J}\_0 \mathcal{S}\_i^2 + \mathcal{O}(\mathcal{S}\_i^3). \tag{26}$$

(ii) Lie scheme coag nuc (Denoted as Lie2)

The analysis is quite similar to case (i), except that the solution procedure for the two joint Eqs. (22a) and (22b) is reversed. The final solution is

$$M\_0^{(\text{OP})} (\delta\_t) = J \delta\_t - \frac{M\_0^{\ast}}{K M\_0^{\ast} \delta\_t - 1} \, \text{} \tag{27}$$

and its Taylor expansion is

$$M\_0^{(\text{OP})} (\mathcal{S}\_t) \equiv M\_0^"+ (\text{K}M\_0^2 + I)\mathcal{S}\_t + \text{K}^2 M\_0^{"\text{3}} \mathcal{S}\_t^2 + \mathcal{O}(\mathcal{S}\_t^3). \tag{28}$$

Hence, the operator splitting error is

$$\operatorname{Err}\_{\text{ouj}} = -M\_0^\dagger \operatorname{IK} \delta\_t^2 + \mathcal{O}(\delta\_t^3) = \frac{1}{2} M\_0^\ast \operatorname{J} \beta\_0 \delta\_t^2 + \mathcal{O}(\delta\_t^3). \tag{29}$$

(iii) Strang scheme 1 2 nuc coag 1 2 nuc (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 0(0) = 0 \* to get the solution 0 (1)(). Second, 0 (1)( 1 <sup>2</sup> ) is used as the initial condition to solve Eq. (22b) to get the solution 0 (2)(). Finally, 0 (2)( ) is used as the initial condition at <sup>=</sup> <sup>1</sup> <sup>2</sup> to solve Eq. (22a) again. The final solution at <sup>=</sup> is found as

$$M\_0^{(\text{CP})}(\delta\_i) = \frac{1}{2}I\delta\_i - \frac{I\delta\_i + 2M\_0^{\ast}}{\text{J}\mathcal{S}\delta\_i^2 + 2KM\_0^{\ast}\delta\_i - 2},\tag{30}$$

and its Taylor expansion is

Hence, 0

Let <sup>=</sup>

and (21)

(ii) Lie scheme

and its Taylor expansion is

(iii) Strang scheme 1

Hence, the operator splitting error is

2

nuc coag

(1)(

14 Aerosols - Science and Case Studies

final solution is found as

) = <sup>+</sup> 0

d

d

Eqs. (22a) and (22b) is reversed. The final solution is

d

 d

coag nuc (Denoted as Lie2)

\*. Take 0

(1)(

*J M M t*

in the above equation, and expand it into Taylor series

(OP) 0 0 \*

*M M KM J K KM J M <sup>t</sup> <sup>t</sup> t t* (OP) \* 2 \* 2 \*2 3 0 00 0 0 ( ) = ( ) ( 2 ) ( ).

*t t tt t t Er M M M JK M J* (OP) \*2 3 \*2 3 0i) 0 0 0 0 0

*t t*

d

*<sup>M</sup> M J KM*

(OP) 0 0 \*

 d

*M M KM J K M <sup>t</sup> t tt* (OP) \* 2 2 \*3 2 3

*tt tt Er M JK M J* \*2 3 \* 2 3 0ii) 0 0 0

 d

<sup>1</sup> <sup>=</sup> ( )= ( ). <sup>2</sup> -+ +

 + ++ + d

0 00 0 ( )= ( )

d

1 2  d

<sup>1</sup> = ( ) ( )= ( )= ( ). <sup>2</sup>

The analysis is quite similar to case (i), except that the solution procedure for the two joint

( )= , <sup>1</sup> - -

*t*

 dd

 bd

nuc (Denoted as Strang1)

 dO O (29)

d

0

\*

 d

 + ++ + + d

*t t*

<sup>+</sup> - + d

( )= .

d

\*

0

1

d

 bd


*JK t KM t*

The operator splitting error can be found by comparing the two Taylor expansion series (25)

) as the initial condition and solve Eq. (22b), then the

(24)

 d

> d

(27)

O( ). (28)

O (25)

$$M\_0^{\text{(CP)}}(\mathcal{S}\_t) = M\_0^\* + (\text{K}M\_0^{\*2} + f)\mathcal{S}\_t + \text{K}(\text{K}M\_0^{\*2} + f)M\_0^\*\mathcal{S}\_t^2 + \frac{1}{4}\text{K}(4\text{K}^2M\_0^{\*4} + 6\text{K}|M\_0^{\*2} + f^2)\delta\_t^3 + \mathcal{O}(\delta\_t^4). \tag{31}$$

Hence, the operator splitting error is

$$E r\_{\text{out}} = \frac{I\mathcal{K}}{12} (2M\_0^{\*2}\mathcal{K} - f)\delta\_\iota^3 + \mathcal{O}(\delta\_\iota^4) = \frac{I\beta\_0}{24} (M\_0^{\*2}\beta\_0 + f)\delta\_\iota^3 + \mathcal{O}(\delta\_\iota^4). \tag{32}$$

(iv) Strang scheme 1 <sup>2</sup>coag nuc 1 2coag (Denoted as Strang2)

Similar to case (iii), the final solution is found as

$$M\_0^{(\text{CP})}(\delta\_t) = 2 \frac{\text{J} \text{K} \text{M}\_0^\* \text{S}\_t^2 - 2 \text{J} \delta\_t - 2 \text{M}\_0^\*}{-\text{J} \text{K}^2 \text{M}\_0^\* \delta\_t^{\text{S}} + 2 \text{J} \text{K} \delta\_t^{\text{S}} + 4 \text{K} \text{M}\_0^\* \delta\_t - 4},\tag{33}$$

and its Taylor expansion is

$$M\_{0}^{(\text{CP})}(\delta\_{i}) = M\_{0}^{"\prime} + \{\text{KM}\_{0}^{"\prime 2} + f\}\delta\_{i} + \text{K}(\text{KM}\_{0}^{"\prime 2} + f\}M\_{0}^{"}\delta\_{i}^{2} + \frac{1}{4}\text{K}(4\text{K}^{2}\text{M}\_{0}^{"\prime 4} + 5\text{K}\text{[M}\_{0}^{"\prime 2} + 2f^{2}\text{)}\delta\_{i}^{3} + \mathcal{O}(\delta\_{i}^{4}).\tag{34}$$

Hence, the operator splitting error is

$$\operatorname{Err}\_{\mathrm{(0v)}} = \frac{I\mathcal{K}}{12} (-\boldsymbol{M}\_0^{\ast 2}\boldsymbol{K} + 2\boldsymbol{f})\boldsymbol{\delta}\_t^3 + \mathcal{O}(\boldsymbol{\delta}\_t^4) = -\frac{I\beta\_0}{48} (\boldsymbol{M}\_0^{\ast 2}\beta\_0 + 4\boldsymbol{f})\boldsymbol{\delta}\_t^3 + \mathcal{O}(\boldsymbol{\delta}\_t^4). \tag{35}$$

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 all is of order , since the error will accumulate in all/ 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 3 in Eqs. (32) and (35)

$$\mathcal{Z}(\boldsymbol{M}\_{\boldsymbol{0}}^{\ast 2}\boldsymbol{\beta}\_{0} + \boldsymbol{I}) - (\boldsymbol{M}\_{\boldsymbol{0}}^{\ast 2}\boldsymbol{\beta}\_{0} + \boldsymbol{4}\boldsymbol{I}) = \boldsymbol{M}\_{\boldsymbol{0}}^{\ast 2}\boldsymbol{\beta}\_{0} - \boldsymbol{2}\boldsymbol{I}.\tag{36}$$

If 1 <sup>2</sup>0 \*20 <sup>&</sup>lt; , i.e., coagulation rate is smaller than nucleation rate, the absolute splitting error in Eq. (32) is smaller than that in Eq. (35). When 1 <sup>2</sup>0 \*20 ≪ , the error in Eq. (33) approaches a half of that in Eq. (36). On the other hand, when 1 <sup>2</sup>0 \*20 ≫ , 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 [38] 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 [38] is that the scheme 2 3(Strang1 + Strang2) + <sup>1</sup> 4(Lie1 + Lie2) is of third order, which can be used to check the accuracy of the above derivations. Keep the third‐order term 3 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 4 and higher orders survive. This is an indirect proof

Following, we will briefly discuss the splitting errors for higher order moments other than 0. Since the total volume is conserved during coagulation, nucleation, and coagulation are independent of each other with respect to 1 (volume fraction). Hence, there are no splitting errors for predicting 1 for 2, the model equation for nucleation and coagulation can be easily derived from Eqs. (6) and (16) as

that all the derivations of the error functions are correct.

$$\frac{\text{d}M\_2}{\text{d}t} = \text{J}\upsilon^{\text{\textquotedblleft}2} + \frac{1}{2}\beta\_0 \int\_0^\upsilon \int\_0^\upsilon (2\upsilon u)\mathfrak{n}(\upsilon, t)\mathfrak{n}(u, t)\text{d}u\text{d}\upsilon = \text{J}\upsilon^{\text{\textquotedblleft}2} + M\_{1\prime}^2. \tag{37}$$

and 1 is determined by nucleation only, which can be integrated as

To summarize, the two Lie schemes have the same magnitude in splitting error, both are of

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 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

> *M J M JM J* \* 2 \* 2 \* 2 0 0 0 0 0 0 2( ) ( 4 ) = 2 .

\*20 <sup>&</sup>lt; , i.e., coagulation rate is smaller than nucleation rate, the absolute splitting

nearly doubles that in Eq. (35). Hence, in the Strang schemes it is preferable to split the higher

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 [38] 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 [38] is that the scheme

(29), then calculate the error in the third‐order scheme constructed as mentioned shortly before,

Following, we will briefly discuss the splitting errors for higher order moments other than 0. Since the total volume is conserved during coagulation, nucleation, and coagulation are independent of each other with respect to 1 (volume fraction). Hence, there are no splitting errors for predicting 1 for 2, the model equation for nucleation and coagulation can be easily

 b

+- + - (36)

<sup>2</sup>0

4(Lie1 + Lie2) is of third order, which can be used to check the accuracy

<sup>2</sup>0

\*20 ≪ , the error in Eq. (33)

3 in the error functions (26), (29), (32), and

4 and higher orders survive. This is an indirect proof

\*20 ≫ , the error in Eq. (32)

all/

the number density would be underpredicted as in Lie1, and vice versa.

bb

error in Eq. (32) is smaller than that in Eq. (35). When 1

of the above derivations. Keep the third‐order term

that all the derivations of the error functions are correct.

it is found that only terms involving

derived from Eqs. (6) and (16) as

approaches a half of that in Eq. (36). On the other hand, when 1

rate process into symmetric substeps to reduce the splitting error.

within a time step. So the global integration error during time

steps. The Lie1 scheme underestimates the number

all is of order

second order in

16 Aerosols - Science and Case Studies

, since the error will accumulate in

3 in Eqs. (32) and (35)

of

If 1 <sup>2</sup>0

2

3(Strang1 + Strang2) + <sup>1</sup>

$$M\_1(t) = f\upsilon^\ast t + M\_1^\ast. \tag{38}$$

Here, \* is the critical volume of nucleated particle, and 1 \* is the initial value at = 0. The solution to Eq. (37) is

$$M\_{2}(t) = \frac{1}{3} f^{2} v^{"2} \beta\_{0} t^{3} + f v^{"} M\_{1}^{"} \beta\_{0} t^{2} + (f v^{"2} + \beta\_{0} M\_{1}^{"}) t + M\_{2}^{"}.\tag{39}$$

Similar to the above error analyses, the four splitting errors for solving Eq. (37) can be readily obtained as

$$E r\_{2(l)} = \frac{1}{3} \beta\_0 l v^\ast (2Jv^\ast \delta\_t + 3M\_1^\ast) \delta\_{t\\_\cdot}^2 \tag{40a}$$

$$E r\_{2\left(\text{ii}\right)} = -\frac{1}{3} \beta\_0 l v^\dagger (f v^\dagger \delta\_\iota + 3M\_\text{1}^\ast) \delta\_\iota^2 \,, \tag{40b}$$

$$E r\_{2\text{(iii)}} = -\frac{1}{12} \beta\_0 I^2 v^{"2} \delta^3\_{\text{t}} \,\text{\textdegree}\tag{40c}$$

$$Err\_{2^{\langle \rm{iv} \rangle}} = \frac{1}{6} \mathcal{B}\_0 I^2 v^{"\-2} \delta\_t^3. \tag{40d}$$

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 <sup>=</sup> <sup>6</sup> 3. The detailed function can also be found in reference [31]. 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 <sup>=</sup> ∫() d. 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

$$\frac{\mathbf{d}M\_0^{(1)}}{\mathbf{d}t} = \mathbf{J},\tag{41a}$$

$$\frac{\text{d}M\_{\text{i}}^{(\text{l})}}{\text{d}t} = \text{J}d\_{p\text{\textquotedblleft}p\text{\textquotedblright}}^{\*}\tag{41b}$$

$$\frac{\mathbf{d}M\_2^{(1)}}{\mathbf{d}t} = \mathbf{J}d\_p^{\*2} \,. \tag{41c}$$

$$\frac{\mathbf{d}M\_3^{(1)}}{\mathbf{d}t} = \mathbf{J}d\_p^{\*3}.\tag{41d}$$

The dynamic equations for moments due to condensation (at a constant diameter growth rate ) are

$$\frac{\partial \mathbf{d}M\_0^{(2)}}{\partial t} = \int -\frac{\partial (\boldsymbol{\eta} \mathbf{G})}{\partial d\_p} \mathbf{d}d\_p = \mathbf{0},\tag{42a}$$

$$\frac{\text{d}M\_1^{(2)}}{\text{d}t} = \int -\frac{\partial(nG)}{\partial d\_p} d\_p \text{d}d\_p = G \Big[ n \text{d}d\_p = \text{G}M\_0^{(2)}\Big] \tag{42b}$$

$$\frac{\text{d}M\_2^{(2)}}{\text{d}t} = \int -\frac{\partial(nG)}{\partial d\_p} d\_p^2 \,\text{d}d\_p = 2\,\text{G} \int \text{nd}\_p \,\text{d}d\_p = 2\,\text{G}M\_1^{(2)}.\tag{42c}$$

Operator Splitting Monte Carlo Method for Aerosol Dynamics http://dx.doi.org/10.5772/65140 19

$$\frac{\mathbf{d}M\_{\mathrm{g}}^{(2)}}{\mathbf{d}t} = \int -\frac{\partial (nG)}{\partial d\_{\mathrm{p}}} d\_{\mathrm{p}}^{3} \mathbf{d}d\_{\mathrm{p}} = \Im G \Big\{ n d\_{\mathrm{p}}^{2} \mathbf{d}d\_{\mathrm{p}} = \Im G M\_{\mathrm{g}}^{(2)}.\tag{42d}$$

Hence, the dynamic moments equations for nucleation and condensation are

function can also be found in reference [31]. 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

is sought to distinguish the symbols for moments defined differently, with the hope to reduce

The dynamic equations for moments due to pure nucleation (to generate particles with

*p <sup>M</sup> Jd <sup>t</sup>* (1) <sup>d</sup> <sup>1</sup> \* = , <sup>d</sup>

*p <sup>M</sup> Jd <sup>t</sup>* (1) <sup>d</sup> <sup>2</sup> \* 2 = , <sup>d</sup>

*p <sup>M</sup> Jd <sup>t</sup>* (1) <sup>d</sup> <sup>3</sup> \* 3 = . <sup>d</sup>

The dynamic equations for moments due to condensation (at a constant diameter growth rate

*p*

*pp p*

*pp pp*

*<sup>M</sup> nG <sup>d</sup> t d* (2) <sup>d</sup> <sup>0</sup> ( ) = d = 0, <sup>d</sup>

*<sup>M</sup> nG d d G n d GM*

*<sup>M</sup> nG d d G nd d GM*

<sup>d</sup> ( ) <sup>=</sup> d =2 d =2 , <sup>d</sup>

2 2 (2)

<sup>d</sup> ( ) = d= d= , <sup>d</sup>

1 (2)

*p*

*p*

*t d*

*t d*

(2)

(2)

*p*

¶ - ¶ <sup>ò</sup> (42a)

0

¶ - ¶ ò ò (42c)

¶ - ¶ ò ò (42b)

1

*<sup>M</sup> <sup>J</sup> <sup>t</sup>* (1) <sup>d</sup> <sup>0</sup> = , <sup>d</sup>

\* at rate ) can be easily obtained with the help of the Delta function (6), which are

d. No attempt

(41a)

(41b)

(41c)

(41d)

with respect to the particle diameter other than the volume <sup>=</sup> ∫()

new symbols and to bring no confusion.

diameter

18 Aerosols - Science and Case Studies

) are

$$\frac{\mathbf{d}M\_0}{\mathbf{d}t} = f\_{\prime} \tag{43a}$$

$$\frac{\text{d}M\_1}{\text{d}t} = \text{J}d\_p^\circ + \text{G}M\_{0^\circ} \tag{43b}$$

$$\frac{\text{d}M\_2}{\text{d}t} = \text{J}d\_p^{\prime 2} + 2\text{GM}\_{1\prime} \tag{43c}$$

$$\frac{\text{d}M\_3}{\text{d}t} = \text{Id}\_p^{\ast 3} + \text{\mathcal{R}}\text{GM}\_2.\tag{43d}$$

Assume the initial condition at = 0 is [0 \*,1 \*,2 \*,3 \*], then Eq. (43) has solution

$$M\_{\mathbf{o}}(\mathbf{t}) = \mathbf{J}\mathbf{t} + M\_{\mathbf{o}}^{\*}\tag{44a}$$

$$M\_1(t) = \frac{1}{2} fGt^2 + (fd\_p^\* + GM\_0^\*)t + M\_{1'}^\* \tag{44b}$$

$$M\_{\rm 2}(t) = \frac{1}{3} \| \mathbf{G}^2 t^3 + (\mathbf{J}d\_p^\prime + \mathbf{G}M\_0^\prime) \mathbf{G} t^2 + (\mathbf{J}d\_p^{\prime 2} + 2\mathbf{G}M\_1^\prime) t + M\_{\rm 2}^\prime \tag{44c}$$

$$M\_3(t) = \frac{1}{4} f \text{G}^3 t^4 + (\text{Id}\_\rho^\* + \text{GM}\_0^\*) \text{G}^2 t^3 + \frac{3}{2} (\text{Id}\_\rho^{\*2} + 2 \text{GM}\_1^\*) \text{G} t^2 + (\text{Id}\_\rho^{\*3} + 3 \text{GM}\_2^\*) t + M\_3^\*. \tag{44d}$$

On the other hand, the operator splitting can be used to solve Eqs. (41) and (42). Parallel to the analyses in Section 3.1, four splitting schemes are considered as follows.

(i) Lie scheme nuc cond (Denoted as Lie1)

After solving Eqs. (41) and (42) sequentially, the final solution at <sup>=</sup> is found as

$$M\_0^{\text{CP}}(t) = ft + M\_{0\prime}^\* \tag{45a}$$

$$M\_1^{\text{CP}}(t) = fGt^2 + (fd\_p^\* + GM\_q^\*)t + M\_{1'}^\* \tag{45b}$$

$$\boldsymbol{M}\_{2}^{\rm{OP}}(t) = \boldsymbol{f}\boldsymbol{G}^{2}\boldsymbol{t}^{3} + \langle\boldsymbol{\mathcal{Z}}\boldsymbol{f}\boldsymbol{d}\_{\boldsymbol{p}}^{\prime} + \boldsymbol{\mathcal{G}}\boldsymbol{M}\_{0}^{\prime}\rangle\boldsymbol{G}\boldsymbol{t}^{2} + \langle\boldsymbol{\mathcal{J}}\boldsymbol{d}\_{\boldsymbol{p}}^{\prime2} + \boldsymbol{\mathcal{Z}}\boldsymbol{\mathcal{G}}\boldsymbol{M}\_{1}^{\dagger}\rangle\boldsymbol{t} + \boldsymbol{M}\_{2}^{\prime}\tag{45c}$$

$$\mathbf{M}\_{\text{3}}^{\text{CP}}\{\mathbf{t}\} = \mathbf{J}\mathbf{G}^{3}\mathbf{t}^{4} + \langle \mathbf{G}\mathbf{J}\mathbf{d}\_{\text{p}}^{\prime} + \mathbf{G}\mathbf{M}\_{\text{0}}^{\prime}\rangle \mathbf{G}^{2}\mathbf{t}^{3} + \langle \mathbf{G}\mathbf{J}\mathbf{d}\_{\text{p}}^{\prime 2} + \mathbf{G}\mathbf{G}\mathbf{M}\_{\text{1}}^{\prime}\rangle \mathbf{G}\mathbf{t}^{2} + \langle \mathbf{J}\mathbf{d}\_{\text{p}}^{\prime 3} + \mathbf{3}\mathbf{G}\mathbf{M}\_{\text{2}}^{\prime}\rangle \mathbf{t} + \mathbf{M}\_{\text{3}}^{\prime}.\tag{45d}$$

The splitting error is found by comparing the splitting solution (45) with the exact solution (44)

$$Er\_{0(i)} = 0,\tag{46a}$$

d2 1(i) 1 = , <sup>2</sup> *<sup>t</sup> Er JG* (46b)

$$\text{Er}\_{2(i)} = \text{Jd}\_p^\ast \text{G} \delta\_t^2 + \frac{2}{3} \text{J} \text{G}^2 \delta\_t^3 \text{.}\tag{46c}$$

$$E r\_{3(l)} = \frac{3}{2} J d\_p^{"\prime 2} G \delta\_t^2 + 2 J G^2 d\_p^{"\prime} \delta\_t^3 + \frac{3}{4} J G^3 \delta\_t^4. \tag{46d}$$

(ii) Lie scheme cond nuc (Denoted as Lie2)

After solving Eqs. (42) and (41) sequentially, the final solution at <sup>=</sup> is found as

$$M\_{\rm{o}}^{\rm{CP}}(t) = \text{Jt} + M\_{\rm{o}}^{\rm{\*}} \tag{47a}$$

$$M\_{1}^{\text{CP}}(t) = (\text{Id}\_{p}^{\text{}} + \text{GM}\_{0}^{\text{}})t + M\_{1}^{\text{'}} \tag{47b}$$

$$\mathbf{M}\_{2}^{\rm{OP}}(t) = \mathbf{G}^{2}M\_{\mathrm{o}}^{\*}t^{2} + (\mathbf{J}d\_{p}^{\*2} + 2\mathbf{G}M\_{1}^{\*})t + \mathbf{M}\_{2}^{\*} \tag{47c}$$

$$M\_3^{\rm OP}(t) \equiv G^3 M\_0^" t^3 + \Im G^2 M\_1^" t^2 + (\{d\_p^{\*3} + \Im G M\_2^"\} t + M\_3^". \tag{47d}$$

The splitting error is found by comparing the splitting solution (47) with the exact solution (44)

$$Er\_{0(u)} = 0,\tag{48a}$$

$$\mathcal{L}\boldsymbol{E}\boldsymbol{r}\_{\text{1(ii)}} = -\frac{1}{2}\mathcal{J}\boldsymbol{G}\boldsymbol{\delta}\_{t}^{2} \,. \tag{48b}$$

$$Err\_{2(0)} = -Jd\_p^"\mathcal{G}\delta\_t^2 - \frac{1}{3}J\mathcal{G}^2\delta\_t^3,\tag{48c}$$

$$\operatorname{Err}\_{\mathfrak{A}(\mathfrak{u})} = -\frac{3}{2} \operatorname{Id}\_{p}^{\ast 2} \mathbf{G} \boldsymbol{\delta}\_{t}^{2} - 2 \operatorname{Id}\_{p}^{\ast} \mathbf{G}^{2} \boldsymbol{\delta}\_{t}^{3} - \frac{1}{4} \operatorname{J} \mathbf{G}^{3} \boldsymbol{\delta}\_{t}^{4}. \tag{48d}$$

(iii) Strang scheme 1 2 nuc cond 1 2 nuc (Denoted as Strang1)

*M t Jt M* OP \*

*M t JGt Jd GM t M <sup>p</sup>* OP 2\* \* \*

*M t JG t Jd GM Gt Jd GM t M p p* OP 2 3 \* \* 2 \*2 \* \*

*M t JG t Jd GM G t Jd GM Gt Jd GM t M ppp* OP 3 4 \* \* 2 3 \*2 \* 2 \*3 \* \*

1(i)

d

2(i)

nuc (Denoted as Lie2)

After solving Eqs. (42) and (41) sequentially, the final solution at <sup>=</sup>

3(i)

(ii) Lie scheme cond

20 Aerosols - Science and Case Studies

0 0 ( )= , + (45a)

0(i) *Er* = 0, (46a)

= , <sup>2</sup> *<sup>t</sup> Er JG* (46b)

0 0 ( )= , + (47a)

<sup>1</sup> 0 1 ( )=( ) , + + (47b)

2 0 1 2 ( )= ( 2 ) , ++ + (47c)

3 01 2 3 ( )= 3 ( 3 ) . + ++ + (47d)

is found as

= , <sup>3</sup> *pt t Er Jd G JG* (46c)

<sup>1</sup> 0 1 ( )= ( ) , ++ + (45b)

<sup>2</sup> <sup>0</sup> 1 2 ( ) = (2 ) ( 2 ) , + + ++ + (45c)

<sup>3</sup> <sup>0</sup> <sup>1</sup> 2 3 ( ) = (3 ) (3 3 ) ( 3 ) . + + + + ++ + (45d)

The splitting error is found by comparing the splitting solution (45) with the exact solution (44)

d2

+ \* 2 23

2

+ + \*2 2 2 \* 3 3 4

3 3 =2 .

 d

dd

2 4 *p t p t <sup>t</sup> Er Jd G JG d JG* (46d)

1

d

*M t Jt M* OP \*

*M t Jd GM t M <sup>p</sup>* OP \*\*\*

*M t G M t Jd GM t M <sup>p</sup>* OP 2 \* 2 \* 2 \* \*

*M t G M t G M t Jd GM t M <sup>p</sup>* OP 3 \* 3 2 \* 2 \* 3 \* \* After sequentially solving Eq. (41) with half time step /2, then Eq. (42) with full time step, and then Eq. (41) with half step again, the final solution at <sup>=</sup> is found as

$$M\_{\rm o}^{\rm CP}(t) = \text{J}t + M\_{\rm o}^{\prime} \tag{49a}$$

$$M\_{\rm 1}^{\rm OP}(t) = \frac{1}{2}I\mathcal{G}t^2 + (\mathcal{J}d\_p^\circ + \mathcal{G}M\_\bullet^\circ)t + M\_{\rm 1}^\circ \tag{49b}$$

$$M\_2^{\rm OP}(t) = \frac{1}{2} f \text{G}^2 t^3 + (f \text{d}\_p^\circ + \text{G} \text{M}\_0^\circ) \text{G} t^2 + (f \text{d}\_p^{\prime 2} + 2 \text{G} \text{M}\_1^\circ) t + M\_{2\prime}^\circ \tag{49c}$$

$$M\_{\odot}^{\text{CP}}(t) = \frac{1}{2}J\mathbf{G}^{\text{3}}t^{4} + (\frac{3}{2}\text{Id}\_{\rho}^{\prime} + \text{GM}\_{0}^{\prime})\mathbf{G}^{2}t^{3} + (\frac{3}{2}\text{Id}\_{\rho}^{\prime 2} + \text{G}\mathbf{G}\mathbf{M}\_{1}^{\prime})\mathbf{G}t^{2} + (\text{Id}\_{\rho}^{\prime 3} + \text{G}\mathbf{G}\mathbf{M}\_{2}^{\prime})\mathbf{t} + \mathbf{M}\_{3}^{\prime}.\tag{49d}$$

The splitting error is found by comparing the splitting solution (49) with the exact solution (44)

$$Er\_{0(iii)} = 0,\tag{50a}$$

$$Er\_{1(iii)} = 0,\tag{50b}$$

$$Err\_{2^{\text{(iii)}}} = \frac{1}{6} fG^2 \delta\_t^3 \,, \tag{50c}$$

$$E r\_{\mathfrak{z}(\mathfrak{ii})} = \frac{1}{2} J G^2 d\_p^\* \delta\_t^3 + \frac{1}{4} J G^3 \delta\_t^4. \tag{50d}$$

(iv) Strang scheme 1 <sup>2</sup>cond nuc 1 2cond (Denoted as Strang2)

After sequentially solving Eq. (42) with half time step /2, then Eq. (41) with full time step, and then Eq. (42) with half step again, the final solution at <sup>=</sup> is found as

$$M\_{\rm o}^{\rm OP}(t) = \text{J}t + M\_{\rm o}^{\rm \, \prime} \tag{51a}$$

$$M\_{1}^{\text{OP}}(t) = \frac{1}{2} f \text{G} t^{2} + (\text{Id}\_{p}^{\*} + \text{GM}\_{0}^{\*})t + M\_{1}^{\*},\tag{51b}$$

$$M\_2^{\text{OP}}(t) = \frac{1}{4} f \text{G}^2 t^3 + (f d\_p^\ast + \text{G} M\_0^\ast) \text{G} t^2 + (f d\_p^{\ast 2} + 2 \text{G} M\_1^\ast) t + M\_{2\prime}^\ast \tag{51c}$$

$$M\_3^{\text{CP}}(t) = \frac{1}{8} \text{J} \text{G}^3 \text{t}^4 + (\frac{3}{4} \text{Id}\_\rho^\circ + \text{GM}\_0^\circ) \text{G}^2 \text{t}^3 + (\frac{3}{2} \text{Id}\_\rho^{\circ 2} + \text{3GM}\_1^\circ) \text{G} \text{t}^2 + (\text{Id}\_\rho^{\circ 3} + \text{3GM}\_2^\circ) \text{t} + M\_3^\circ. \tag{51d}$$

The splitting error is found by comparing the splitting solution (51) with the exact solution (44)

$$Er\_{0(lv)} = 0,\tag{52a}$$

$$Er\_{\rm I(iv)} = 0,\tag{52b}$$

$$Er\_{2(\nu)} = -\frac{1}{12} J G^2 \delta^3\_{\text{t}} \,\text{\textdegree}\tag{52c}$$

$$E r\_{\mathfrak{A}(\mathsf{v})} = -\frac{1}{4} J G^2 d\_p^\ast \delta\_t^3 - \frac{1}{8} J G^3 \delta\_t^4. \tag{52d}$$

To summarize, there is no splitting error for the number density (0), which is solely deter‐ mined by nucleation. For 1 to 3, 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 1 is also zero, and the errors for 2 and 3 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 subse‐ quent 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., 0 \*, 1 \*, 2 \*, and 3 \*, which seems strange at first glance. In fact, the splitting errors are anchored to the moment's evolution. Take the 1 splitting error in the Lie1 scheme as example, the error contains the factor , while the leading order in the analytical solution for 1 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.
