Langevin Dynamics Calculation of Brownian Coagulation Coefficient for Spherical Equal-size Aerosol Particles in Transient Regime

Coagulation coefficient of aerosol particles due to Brownian motion is an important issue to describe change in particle size distribution. Motion of aerosol particles is diffusive in continuous region (small Knudsen number; Kn), or like free molecular motion of gaseous molecular in free molecular region (large Kn). Fuchs (1964) presented an expression of coagulation coefficient in transition regime by a so-called “Flux Matching” method. In his method, transportation of particles inside of the “limiting sphere” is assumed to be like free molecular, or diffusive outside of the sphere. These days, some researchers presented coagulation coefficient of aerosol particles by direct calculation of motion of aerosol particles. They employed Langevin dynamics equation to represent the stochastic motion of aerosol particles. In this study, we developed new model to calculate the coagulation coefficient. Our model employed spherical calculation space in which one scavenging particle is in the center of it: the calculation sphere moves together with the motion of the scavenging particle. The coagulation coefficient can be calculated from the mean time between collisions and the concentration of collision particles. By using the above numerical model, we have calculated the coagulation coefficient of spherical particles of from 4 nm to 100 nm in diameter.


Introduction
Condensation of vapor to particles and coagulation between particles are major process leading to change in particle size distribution. In a free molecular regime ( ∞ → Kn ) the coagulation rate can be expressed by cross section of particle and relative thermal speed between particles, in contrast by diffusion flux of particles in continuum regime ( 0 → Kn ). Fuchs (1964) proposed an expression of coagulation coefficient β by means of "Flux-matching" method in which particle diffusive transport and particle motion like free molecular were combined. Fuchs' expression remains the most widely used, based on reasonable agreement with experimental measurements (Szymanski et al., 1989;Kim et al., 2003).
In recent years, owing to the progress of computer, kinetic approach has been becoming realizable (Gopalakrishnan and Hogan Jr., 2011;Polovnikov et al., 2016), however the capability of computers is not enough to reproduce the motion of aerosol particles in large scale. In this study, a novel kinetic numerical model focusing on both diffusive nature of particles far from other particle and like free molecular nature was developed. Fuchs (1964) estimated the coagulation coefficient of spherical particles as follows. Imaging one of the particles of Dpi in diameter to be stationary, while other particles of Dpj collide with the stationary particles due to Brownian motion. Since two spherical particles come into contact when the distance between their centres is equal to (Dpi+Dpj)/2, the stationary particle can be replaced by an absorbing sphere of radius (Dpi+Dpj)/2, and the other particles by mass points. Assuming the continuum regime ( 0 D → Kn ), concentration distribution of particle of Dpj around the absorbing sphere can be described by diffusion equation in spherical coordinate by using the effective diffusion coefficient De=Di+Dj, which is expressed as C(r)=A1/r+A2 (A1 and A2 are constants). Otherwise, the existence of so-called "limiting-sphere" is assumed: it is concentric with the absorbing sphere, and Dpi+ Dpj+2gij in diameter, and the motion of particle inside the sphere cannot be described by diffusion equation. Thickness of limiting sphere gij is roughly estimated as mean distance of mean free path of particle from the surface of absorbing sphere in all directions. In the limiting sphere, the particles move as if they were in a vacuum with a mean relative speed of cij.

Flux matching method and resent studies 2.1 Flux matching method of Fuchs
From the above assumption, Fuchs expressed the coagulation coefficient by assuming that the particles enter the limiting sphere by diffusion colloid with absorbing sphere by free molecular like motion as ) ( Note that Fuchs supposed that the particles were diffused from infinity, and accepted the discontinuity of concentration distribution at the surface of limiting sphere.

Recent studies
Gopalakrishnan and Hogan Jr. (2011) presented a dimensionless coagulation coefficient as (2) Here, µij and fij are reduced mass (mimj/(mi+mj)) and reduced friction coefficient of two particles (fifj/(fi+fj)), respectively. The friction coefficient of single particle is expressed as Here, k and T are Boltzmann constant and absolute temperature, respectively. Supposing the relative motion of two particles, they employed the diffusive Knudsen number expressed as Here, C1, C2, C3, and C4 are tuning constant that describe the transition regime. Nondimensionalized coagulation coefficient H are plotted in Figure 1. Grey dashed line and dotted line are free molecular and continuum limiting expressions by Eqs. (5) and (6), respectively. Black solid line is Fuchs' coagulation coefficient (Eq. (1)) nondimensionalized by Eq.
(2). It approaches asymptotically to grey lines at both large and small KnD. Gopalakrishnan and Hogan Jr. (2011) presented the dimensionless coagulation coefficient in transition regime by computing "mean first passage time" by kinetic method. They calculated motion of one or two particles in a cubic space with a side length of 25Dp by Langevin dynamics method (LD method; Chandrasekhar, 1943).
They presented their tuning parameter of transition regime as C1=25. 836, C2=11.211, C3=3.502, and C4=7.211. Their results (denoted as "G&H" in the legend) are also shown in Figure 1 by blue line. The difference between Fuchs' expression and their results are not standing out in the log-log scale plot, however there is meaningful difference between them.

Numerical Method
In this study, we developed a kinitic method like Gopalakrishnan and Hogan Jr. (2011) to calculate coagulation coefficients focusing on the concentration distribution of particle. A schematic illustration is shown in Figure 2. We assume a "scavenging particle" of Dp1 in diameter and concentric "calculation sphere" of DCS=100Dp1 in diameter. The scavenging particle has absorbing surface: other particles which collide with the scavenging particle are vanished. The calculation sphere moves together with the motion of the scavenging particle. Multiple "collision particles" (N=10) of Dp2 are generated randomly in the calculation sphere in an initial stage of the calculation. The motion of both scavenging particle and collision particles are calculated by using Langevin dynamics method (Ermak and Buckholz, 1980) at every time steps. If one of the collision particles escapes from the calculation sphere (a), it re- entries from opposite position of the scavenging particle (b). At every time steps, distances from center of the scavenging particle to center of collision particles are calculated. They are inspected if the collision particle collides with the scavenging particle, or distance < (Dp1+Dp2)/2. If the collision is occurred (c), the elapsed time from a previous collision, c t ∆ is recorded, then the collision particle is re-generated at random position of the surface of calculation sphere (d). In order to evaluate concentration distribution of the collision particles around the scavenging particle, space of the calculation sphere is divided into 640 concentric spherical shells in inverse of radius from radius of collision sphere to radius of calculation sphere. During the continuous numerical procedure of Ntotal time steps, number of collision particles which centers are in the i-th spherical shell, Ni are accumulated in every time steps.
Number concentrations of collision particles in i-th spherical shell of Vi in volume are evaluated as Ni/(Vi Ntotal).
Since it is concerned that static concentration distribution is not established in an initial stage of above procedure, accumulation of concentration distribution and record of time between collision is started after 10 collisions. In this study, concentration distribution of collision particles and mean time between collision ( c t ∆ ; MTBC) are accumulated through more than 100,000 collision event, in order to reduce relative errors in Monte Carlo method which is estimated as number trial 3 . Concentration of collision particles at infinity is estimated by least square approximation from concentration distribution of collision particles far from scavenging particle in which horizontal axis is 1/r. Assuming that the collision between the scavenging particle of (1/V) in concentration and collision particles of (N/V) in concentration, change in concentration of collision particles is Here, assuming that change in number of collision particles (∆N) equal -1 during mean time between collisions ( c t ∆ ), coagulation coefficient βLD can be written as In this study, all calculations were carried out at the condition where T=298.15 K and particle density of 1000 kg m -3 . Diameter of collision particles Dp2 are equal to that of scavenging particle Dp1. Dp1 and Dp2 are varied from 4 nm to 100 nm.

Concentration distribution
Concentration distribution of the collision particles around the scavenging particle with a horizontal axis of r and 1/r are shown in Figures 3 and 4, respectively. Both the scavenging particle and collision particles are 4 nm in diameter, and the diffusive Knudsen number, KnD of 10.6. Vertical dotted line, dashed line, and alternate long and short dash line in grey indicate the position of collision sphere, limiting sphere predicted by Fuchs, and calculation sphere, respectively. Horizontal black  dashed line is over-all concentration of collision particle calculated from the number of collision particle N, and volume of calculation sphere V as N/V. Red solid line shows concentration distribution of the collision particles obtained by LD calculation which is accumulated during through the calculations. Thickness of limiting sphere gij was much larger than the particle diameter, because the calculations were carried out in near free molecular regime. The concentration of collision particles is almost constant at the outside of the limiting sphere, and decreased near the scavenging particle. From the 1/r plot shown in Figure 4, concentration distribution inside the limiting sphere of Fuchs are not on the straight line, because the transport of collision particles near scavenging particle are not governed by diffusion equation but motion of particles like free molecular. Concentration distribution of 1/r plot near surface of calculation sphere was approximated by straight line by mean least square method. Blue dashed line is an extrapolated line from surface of collision sphere to infinity. Intercept of the line in 1/r plot are the concentration of collision particles at infinity, C should be larger than N/V in order to counter balance the decrease in concentration near the scavenging particle, however the increase in ∞ 2 C was quite small due to very thin diffusion layer of collision particles at this condition as shown in Table 1.
Concentration distribution obtained by the calculations for particles of 20 nm in diameter, KnD of 0.9862 are shown in Figures 5 and 6, in the same manner. As shown in Figure 6, decrease in concentration of collision particles at outside of limiting sphere is clear, and the 1/r plot shows that the concentration distribution is straight wise but slightly off the straight line at near the surface of limiting sphere. This means that particle  transport far from the scavenging particle are governed by diffusion equation. Due to increase in the thickness of diffusion layer, relative increase of collision particle concentration ∞ 2 C at infinity is increased to +0.97%. Concentration distribution obtained by the calculation for particles of 100 nm in diameter, KnD of 0.110 are shown in Figures 7 and 8, in the same manner. At this condition, the thickness of diffusion layer is more increased, and the concentration of collision particles at the surface of limiting sphere is decreased. As shown in 1/r plot in Figure 8, concentration distribution outside the limiting sphere are almost in straight wise, but slightly off the straight line at near the surface of limiting sphere: this tendency is same to Figure 6. Due to significant increase in the thickness of diffusion layer, relative increase of collision particle concentration at infinity ∞ 2 C was increased to +2.57%. Considering the increase in concentration at infinity,

Coagulation coefficient
From the estimated concentration of collision particles at infinity (2) then we found tuning parameters of Eq. (7) Figure 9 by red line.
As described in Section 2, Fuchs derived his expression of the coagulation coefficient by dividing the space around the absorbing sphere by the surface of limiting sphere into the space in where particles are transported by diffusion and the space in where the particles move like in vacuum.
He allowed the discontinuity of concentration distribution at the surface of limiting sphere. As shown in Figures 3-8, we confirmed the space where the particle transport cannot be described by diffusion equation as Fuchs predicted. However, the concentration distribution resulted in our calculation are on continuous line. In other word, there are no discontinuity in concentration distribution. Continuous change in the nature of particle motion should cause the difference in the coagulation coefficient.
The mean first passage time method by Gopalakrishnan and Hogan Jr. (2011) is same to our calculation method in terms of the kinetic method, however our results of coagulation coefficients didn't agree with their results. For the calculation of large KnD condition, their cubic space used in their method was seems to be not enough to represent the particle motion like free molecular. For instance, mean free path of particle of 4nm in diameter is 64.79 nm. In our calculations, we employed the calculation sphere of 400 nm in diameter, but cubic space with a side length of 100 nm was employed in their calculation in their study. For the calculating the coagulation coefficient of small KnD number condition, we used the extrapolated concentration at infinity and N/V are increased as KnD decreased: the difference should be more enhanced when the calculation sphere is smaller.

Conclusions
In this study, we developed a new kinetic model to represent the Brownian coagulation of aerosol particles. Our numerical method is kinetic one, and similar to Gopalakrishnan and Hogan's model (2011), however, we employed larger space to reproduce the Brownian particle motion and continuous calculations were carried out. As the results the limiting sphere around a particle predicted by Fuchs (1964) was almost reproduced, and our results successfully elucidated that the concentration distribution of particles is continuous even in transition regime.