Calculation of the Ion-Ion Recombination Rate Coefficient via a Hybrid Continuum- Molecular Dynamics Approach
Tomoya Tamadate, Hidenori Higashi, Takafumi Seto, Christopher Hogan
11 Calculation of the Ion-Ion Recombination Rate Coefficient via a Hybrid Continuum- Molecular Dynamics Approach
Tomoya Tamadate , Hidenori Higashi , Takafumi Seto , Christopher J. Hogan Jr. Department of Mechanical Engineering, University of Minnesota, 111 Church St SE, Minneapolis, MN, 55455 Faculty of Natural System, Graduate School of Natural Science and Technology, Kanazawa University, Kanazawa, Japan Faculty of Frontier Engineering, Institute of Science and Engineering, Kanazawa University, Kanazawa, Japan
The following article has been published in the Journal of Chemical Physics as: Tamadate T., Higashi H., Seto T., & Hogan C. J. Calculation of the Ion-Ion Recombination Rate Coefficient via a Hybrid Continuum- Molecular Dynamics Approach. The Journal of Chemical Physics. 152: 094306. Doi: 10.1063/1.5144772. Please refer to The Journal of Chemical Physics (https://doi.org/10.1063/1.5144772) for the final published version and supporting information. * To whom correspondence should be addressed: E-mail: [email protected], T: 1-612-626-8312
ABSTRACT
Accurate calculation of the ion-ion recombination rate coefficient has been of long-standing interest, as it controls the ion concentration in gas phase systems and in aerosols. We describe the development of a hybrid continuum-molecular dynamics approach to determine the ion-ion recombination rate coefficient. The approach is based on the limiting sphere method classically used for transition regime collision phenomena in aerosols. When ions are sufficiently far from one another, ion-ion relative motion is described by diffusion equations while within a critical distance, molecular dynamics (MD) simulations are used to model ion-ion motion. MD simulations are parameterized using the AMBER force-field as well as by considering partial charges on atoms. Ion-neutral gas collisions are modeled in two mutually exclusive cubic domains composed of 10 gas atoms each, which remain centered on the recombining ions throughout calculations. Example calculations are reported for NH recombination with NO in He, across a pressure range from 10 kPa to 10,000 kPa. Excellent agreement is found in comparison of calculations to literature values for the 100 kPa recombination rate coefficient (1.0 x 10 -12 m s -1 ) in He. We also recover the experimentally observed increase in recombination rate coefficient with pressure at sub-atmospheric pressures, and the observed decrease in recombination rate coefficient in the high pressure continuum limit. We additionally find that non-dimensionalized forms of rate coefficients are consistent with recently developed equations for the dimensionless charged particle-ion collision rate coefficient based on Langevin dynamics simulations. Symbol Dictionary
Constants e
Electron charge 1.60 x 10 -19 [C] k b Boltzmann constant 1.38 x 10 -23 [J/K -1 ] T Temperature 298.15 [K] Permittivity of vacuum 8.85 x 10 -12 [Fm -1 ] Variables a ij Collision radius of ion i and j [m] b Collision parameter [m] b c Critical collision parameter [m] D i , D j Diffusion coefficient of ion i and j [m s -1 ] E Electric field [Vm -1 ] H Dimensionless collision kernel [-] I Ion flow [s -1 ] I Ion flow into limiting sphere surface [s -1 ] J Ion flux density [m -2 s -1 ] Kn D Diffusive Knudsen number [-] Kn Diffusive Knudsen number for limiting [-] sphere m i , m j Mass of ion i and j [kg] m ij Reduced masses of ion i and j [kg] = m i m j /( m i + m j ) n i , n j Ion i and j concentrations [m -3 ] n ∞ Ion i bulk concentrations [m -3 ] N col Number of collisions observed in the [-] simulation N obs Number of orbits observed in the [-] simulation N tot Total number of simulated trajectories [-] P Pressure [Pa] p Collision probability for ions entering [-] the limiting sphere p Collision probability for ions entering [-] the trapping sphere p F Fuchs’s collision probability [-] r Radial position (ion-ion distance) [m] r km Radial distance for atoms k and m [m] r m Minimum radial distance (ion-ion distance) [m] t Time [s] v Velocity at the limiting sphere surface [m s -1 ] v e Electrophoretic velocity at the limiting [m s -1 ] sphere surface v th Thermal velocity at the limiting sphere [m s -1 ] surface z k , z m Partial charge of specie k and m [-] ij he collision rate coefficient for type j ions [m s -1 ] with type i ions which have arrived at a ij he collision rate coefficient for type j ions [m s -1 ] with type i ions which have arrived at Trapping sphere radius [m] Limiting sphere radius [m] km Lenard-Jones parameter between atoms [J] k and m inter Intermolecular potential [J] LJ Lenard-Jones potential [J] E Electric potential [J] Enhancement factor [-] c , f Dimensionless continuum regime and free [-] molecular enhancement factors ij Ion-ion mean free path [m] Incident angle at limiting sphere surface [-] c Critical angle for collision [-] c, Critical angle to enter trapping sphere [-] km Lenard-Jones parameter for atoms k and m [m] E Magnitude of the Coulombic energy to [-] the thermal energy ratio at the limiting sphere surface
I. Introduction
Ion-ion recombination refers to instances where ions of opposing polarity collide with one another, and via charge transfer reactions the products are neutral vapor phase species.
Recombination is crucially important in many gas phase and aerosol systems; the rate of recombination relative to the rate of ion formation determines ion concentrations in the ambient
5, 6 , in high temperature combustion systems , and in intentionally irradiated systems utilized in modulating the charge on aerosol particles (i.e. in aerosol chargers
8, 9 ). Ion concentrations, in turn, can affect gas phase reaction and ambient particle formation rates , and via the charging of particles, particle-particle interactions and particle measurements.
14, 15
Gas phase and aerosol systems are thus highly sensitive to ions and ion-ion recombination. The recombination rate itself is the product of the number concentrations of the colliding ions, n i and n j for ion species i and j , respectively, and the recombination rate coefficient, β 𝑖𝑗 (which is typically denoted with the symbol , but which we denote as β 𝑖𝑗 for reasons noted subsequently). Recombination is typically considered to be a collision-controlled process, hence β 𝑖𝑗 is determined by the motion of two ions about one another until collision. While on the surface this appears to be a relatively simple process, calculation and measurement of the recombination rate has been a topic of scientific interest for more than a century,
2, 4, 16-18 and improved recombination rate coefficient calculations still remain of interest. Along these lines the purpose of this work is to develop a new method of calculating the ion-ion recombination rate coefficient, accounting precisely for the influences of ion structures, ion-ion potential interactions, and ion-neutral background gas collisions on the recombination process, at variable temperatures and pressures.
To better justify development of a new calculation approach, we first summarize prior ion-ion recombination rate coefficient measurements and calculations, noting why specifically improved calculation methods are important. Under tropospherically relevant conditions, measurements repeatedly suggest a recombination rate coefficient of order 1.0-3.0 x 10 -12 m s -1 .
4, 17, 19-21
In general, the recombination rate coefficient increases with increasing pressure at sub-atmospheric levels, and increases with decreasing temperature.
Further observations are dependencies on the ion chemical composition and background gas relative humidity (which also affects ion composition
5, 25 ), as well as a decrease in the rate with increasing pressure at elevated, super-atmospheric pressures.
4, 20, 26
Calculation approaches that can accurately predict how the recombination rate coefficient changes above or below atmospheric pressure, at variable temperature, gas composition, and ion composition are still lacking. The earliest calculation approaches relevant to atmospheric ions were performed by Thomson. Discussed in detail by Loeb & Marshall, Thomson modeled recombination such that even at low pressure, ion recombination is influenced by three-body trapping, i.e. as two oppositely charged ions migrate about one another, there is a non-negligible probability that one of them collides with a neutral gas molecule, and this ion-neutral collision alters the ion trajectory such that an ion-ion collision subsequently occurs when otherwise it would not. Three-body trapping can be contrasted with the two-body collision approach, where ion-neutral gas collisions are negligible, and which is widely utilized in modeling high energy electron collisions with positive ions and positively charged particles and probes in low pressure plasmas.
16, 27
It can also be contrasted with the continuum (Langevin) approach, where ion motion is resisted by multiple neutral gas collisions, as would be the case in a high pressure gas. Calculation approaches to bridge the combined two-body collision, three-body trapping low pressure regime and the continuum, high pressure regime were later developed by Natanson, Brueckner, and Bates & Flannery, all treating ions as structureless point entities. Subsequent calculation efforts were devoted to modeling this process using Monte-Carlo simulations
30, 31 and to scrutinizing implementation of combined three-body trapping- two-body collision models near atmospheric conditions. Overall, while agreement between theoretical calculations, Monte-Carlo simulations, and measurements can be achieved (often with some tuning of input values), no single calculation approach has proven applicable over a wide range of ion chemical compositions, temperatures, pressures, and gas compositions (see these references
5, 24, 25, 32 for examples of disagreement between theory and measurements). The issues with prior calculation approaches arise largely because it is difficult first to analytically model ion-ion relative motion across a wide pressure range, and second to model ion-neutral gas collisions without precisely accounting for ion-neutral interactions and instead treating ions as structureless point charges. To accurately capture ion-neutral gas collisional influences, molecular dynamics (MD) simulations, or at least gas molecule scattering calculations
33, 34 need to be carried out. Such calculations need to consider realistic ion and gas molecule structures (all-atom models), accurate ion-neutral potential interactions, the appropriate gas molecule kinetic energy (velocity) distribution for the temperature of interest, and a procedure to accurately describe gas molecule impingement and rebound during collision with an ion. While the solution would seem to be to simply apply detailed MD simulations to model ion-ion recombination, even with modern computational power, it is generally not feasible to carry out complete MD simulations modeling ion-ion motion in a domain large enough to circumscribe the region where ions Coulombically interact with one another. With this issue in mind, to implement MD simulations to precisely examine ion trajectories on close approach to one another and the possibility of multiple ion-neutral gas molecule collisions, here we develop and implement a hybrid continuum-molecular dynamics approach, wherein motion is modeled via continuum equations when ions are sufficiently far from one another, and via MD simulations as ions approach one another. Such approaches have been developed previously to investigate non-continuum fluid flow ; here the approach is simplified in that the continuum portion of the model can be treated analytically and the continuum model and the molecular dynamics model can be solved decoupled from one another. In combining spatially separated regions where ion migration is modeled via continuum and non-continuum approaches respectively, our approach builds upon the limiting sphere models of particle ionization (particle-ion collisions) of Fuchs and more directly, Filippov , who originally described a universal continuum-non-continuum limiting sphere model that is readily implementable with MD simulations. In utilizing an approach applicable to ion-ion collisions as well as particle-ion and particle-particle collisions, we also attempt to unify theories describing binary collision-limited reactions via a single calculation framework (hence utilizing 𝛽 𝑖𝑗 as the recombination coefficient, as it is typically used to denote the collision rate coefficient for all possible colliding partners in aerosols). In the section that follows, we describe the equations utilized in the hybrid model. A test case is presented of NH recombination with NO in a variable pressure Helium environment (10 kPa – 10 kPa) at 300 K. This test case is largely chosen for computational simplicity for initial model development (a monoatomic noble gas). Results are compared to the near atmospheric pressure recombination rate measurements of Lee & Johnsen ; this comparison shows excellent agreement near atmospheric pressure, where measurements were performed. II. Theory & Numerical Methods
We first present the limiting sphere approach of Fillipov, which, although not the original work presenting the limiting sphere approach, provides the most general form for this model, and forms the basis for the approach utilized in this study (section II.A). Prior to discussing MD calculations, we also provide a review of Fuchs’s assumptions in recombination rate coefficient and particle-ion collision rate coefficient calculation via the limiting sphere approach (section II.B), as well as a presentation of Hoppel & Frick’s modifications to Fuchs’s approach (section II.C.). We elect to review these two, alternative approaches because we believe full presentation of these theories and their limitations better motivates the development of a new calculation method. We additionally compare MD calculation results to predictions from these theories. Nonetheless, readers not concerned with the calculation details of prior theories may proceed to second II.D following section II.A., with minimal loss in scope. Following discussion of MD simulations in section II.D., we then discuss non-dimensionalization of the recombination rate coefficient (section II.E.). A. Limiting Sphere Theory
We consider a domain where ion j is positioned at the center, and the relative motion of type i ions is monitored about this central ion (Figure 1). The limiting sphere radius, , is defined such that at distances beyond it, relative motion can be described fully by continuum diffusion equations. Therefore, the spatial concentration evolution of ions of type i beyond obeys the equation: 𝜕𝑛 𝑖 (𝑟)𝜕𝑡 + 𝜕𝐽𝜕𝑟 = 0 (1) where: 𝐽 = −(𝐷 𝑖 + 𝐷 𝑗 ) 𝜕𝑛 𝑖 (𝑟)𝜕𝑟 + 𝑒𝑘 b 𝑇 (𝐷 𝑖 + 𝐷 𝑗 )𝐸(𝑟)𝑛 𝑖 (𝑟) (2) n i is the concentration of type i ions, t is time, r is radial position (ion-ion distance), e is the unit electron charge, D i is the diffusion coefficient of ion i , 𝑘 b is Boltzmann’s constant, T is temperature, and E is the electric field formed by the ions. In equation (2) both ions are assumed singly charged. In the present study, we also assume a simple Coulomb form for the electric field outside the limiting sphere radius, 𝐸(𝑟) = −𝑒4𝜋ε 𝑟 as short range potential interaction terms are negligible at long distances ( ε is the permittivity of free space). Combining equations (1) and (2) and assuming that the concentration profile rapidly reaches a steady-state yields: (𝐷 𝑖 + 𝐷 𝑗 ) 𝜕𝑛 𝑖 (𝑟)𝜕𝑟 + (𝐷 𝑖 + 𝐷 𝑗 ) 𝑒 𝑘 b 𝑇𝑟 𝑛 𝑖 (𝑟) = 𝐴 (3) where 𝐴 is a constant. The total flow of type i ions to any radial coordinate equal to or larger than the limiting sphere radius is then defined as: 𝐼 = 4𝜋𝑟 𝐴 = 4𝜋𝑟 (𝐷 𝑖 + 𝐷 𝑗 ) 𝜕𝑛 𝑖 (𝑟)𝜕𝑟 + (𝐷 𝑖 + 𝐷 𝑗 ) 𝑒 𝜀 𝑘 b 𝑇 𝑛 𝑖 (𝑟) (4) Equation (4) is a first order differential equation, whose solution for the boundary condition 𝑛 𝑖 → 𝑛 ∞ as 𝑟 → ∞ was determined by Fuchs to be: 𝐼 = 4𝜋(𝐷 𝑖 + 𝐷 𝑗 ) 𝑛 ∞ −𝑛 𝑖 (𝑟)𝑒𝑥𝑝( −𝑒24𝜋𝜀0𝑘b𝑇𝑟 )∫ 𝑒𝑥𝑝( −𝑒24𝜋𝜀0𝑘b𝑇𝑥 )𝑑𝑥 = (𝐷 𝑖 +𝐷 𝑗 )𝑒 𝜀 𝑘𝑇 𝑛 ∞ −𝑛 𝑖 (𝑟)𝑒𝑥𝑝( −𝑒24𝜋𝜀0𝑘b𝑇𝑟 )1−𝑒𝑥𝑝( −𝑒24𝜋𝜀0𝑘b𝑇𝑟 ) (5) Equation (5) would apply to the point of collision (i.e. r = a ij , where a ij is the point of ion-ion contact) if continuum transport equations were valid for all ion-ion separation distances. Therefore, equation (5) can be applied up to a point r = and then equated with the product of the type i ion concentration at and the collision rate coefficient for type j ions with type i ions which have arrived at β δ : 𝐼 δ = (𝐷 𝑖 +𝐷 𝑗 )𝑒 𝜀 𝑘𝑇 𝑛 ∞ −𝑛 𝑖 (δ)𝑒𝑥𝑝( −𝑒24𝜋𝜀0𝑘b𝑇δ )1−𝑒𝑥𝑝( −𝑒24𝜋𝜀0𝑘b𝑇δ ) = 𝑛 𝑖 (δ)β δ (6) Filippov showed that this limiting sphere collision rate coefficient can be calculated via the equation: β δ = 𝑝 δ (2−𝑝 δ ) ( b 𝑇𝑚 𝑖𝑗 ) (7) where 𝑚 𝑖𝑗 is the reduced mass of ions and 𝑝 δ is the probability that a type i ion entering the limiting sphere will collide with ion j . The recombination coefficient, β 𝑖𝑗 can then be defined by rearrangement of equation (6). Defining β 𝑖𝑗 = 𝐼 𝛿 𝑛 ∞ , and noting that 𝑛 𝑖 (δ) = 𝐼 δ β δ = β 𝑖𝑗 𝑛 ∞ β δ yields: β 𝑖𝑗 = (𝐷 𝑖 +𝐷 𝑗 )𝑒 𝜀 𝑘𝑇(1−𝑒𝑥𝑝( −𝑒24𝜋𝜀0𝑘b𝑇δ )) (1 + (𝐷 𝑖 +𝐷 𝑗 )𝑒 𝛽 δ 𝜀 𝑘 b 𝑇(𝑒𝑥𝑝( 𝑒24𝜋𝜀0𝑘b𝑇δ )−1) ) −1 (8a) β 𝑖𝑗 = 𝑖 +𝐷 𝑗 )δΨ δ (1−𝑒𝑥𝑝(−Ψ δ )) (1 + ( 𝜋2 ) δ )𝑝 δ 𝐾𝑛 δ Ψ δ (𝑒𝑥𝑝(Ψ δ )−1) ) −1 (8b) Ψ δ = 𝑒 𝑘 b 𝑇δ (8c) 𝐾𝑛 𝛿 = ( 𝑚 𝑖𝑗 𝑘 b 𝑇 ) 𝑖 +𝐷 𝑗 )δ (8d) In equation (8c), Ψ δ is the magnitude of the Coulombic energy to the thermal energy ratio at the limiting sphere surface, and in equation (8d), 𝐾𝑛 δ is a diffusive Knudsen number for the limiting sphere. Though equations (8b-d) are uniquely represented here (i.e. we do not believe these equations have been presented in this manner previously), their derivation follows the approach of Fillipov and as noted by Fillipov, equation (8a) is similar to that derived by Fuchs in developing the limiting sphere model. Fillipov additionally notes that equation (8) can be implemented for any selected value of , provided that beyond the continuum transport approximation is valid; in fact, the continuum transport approximation can remain valid within a portion of , provided that a method to accurately determine 𝑝 𝛿 for properly chosen values of is implemented. For this reason, equations (8a-d) are applied with the MD simulation approach discussed in section II.D. Figure 1.
A depiction of ion motion to the limiting sphere surface (upper left). Determination of the probability of collision for an ion initiated on the limiting sphere surface (entering ion) with an ion initiated at the center (center ion) of the sphere yields the ion-ion recombination rate via equation (8). Molecular Dynamics simulations (depicted in the upper right) can be used in lieu of traditional theories to compute the probability of collision accounting for ion-neutral gas encounters. Ion-ion motion yields either collision, non-collision, or orbit (lower line) and any theory defining the probability of collision must explicitly state how each event is defined and accounted for collision probability analysis. B.. Fuchs’s Assumptions in Probability Calculation
In developing and utilizing the limiting sphere approach to determine particle-ion collision rate coefficients (with extension to ion-ion recombination coefficients), Fuchs utilized equation (8a) with very specific assumptions. First, following the discussion of Wright, he selected the limiting sphere radius as: δ = 𝑎 𝑖𝑗3 𝜆 𝑖𝑗2 [ (1 + 𝜆 𝑖𝑗 𝑎 𝑖𝑗 ) − (1 + 𝜆 𝑖𝑗2 𝑎 𝑖𝑗2 ) (1 + 𝜆 𝑖𝑗 𝑎 𝑖𝑗 ) + (1 + 𝜆 𝑖𝑗2 𝑎 𝑖𝑗2 ) ] (9a) where 𝑎 𝑖𝑗 is the collision radius (the point at which the two ions or ion and particle collide), and λ 𝑖𝑗 , commonly referred to as the ion-ion mean free path, is given as: 𝜆 𝑖𝑗 = (𝐷 𝑖 + 𝐷 𝑗 ) ( 𝜋𝑚 𝑖𝑗 b 𝑇 ) (9b) Second, Fuchs assumed that ion relative motion inside the limiting sphere can be described completely by free molecular (collisionless) kinetics, hence can be expressed as: β δ = 𝜋δ ( b 𝑇𝜋𝑚 𝑖𝑗 ) 𝑝 F (10) Equation (10) is a free molecular collision rate coefficient (projected area and mean thermal speed product) multiplied by the probability ( 𝑝 F ) that an ion impinging upon limiting sphere does in fact collide with the central ion. To calculate 𝑝 F , Fuchs utilized simplified trajectory calculations for point ions (with a charged particle at the center), wherein ions were initiated on the limiting sphere surface with a single speed 𝑣 . Conservation of energy and conservation of angular momentum during ion migration yield: 𝑚 𝑖𝑗 𝑣 − 𝑒 δ = 𝑚 𝑖𝑗 [( 𝑑𝑟𝑑𝑡 ) + 𝑟 ( 𝑑𝜃𝑑𝑡 ) ] − 𝑒 𝑟 (11) 𝑏𝑚 𝑖𝑗 𝑣 = 𝑚 𝑖𝑗 𝑟 ( 𝑑𝜃𝑑𝑡 ) . (12) where 𝑟 and 𝜃 are the radial and angular coordinates of the incoming ion (with the central ion remaining fixed in the examined frame of reference), 𝑏 = δsin𝜃 is the initial impact parameter (Figure 1), 𝜃 is the initial entering angle, and 𝑑𝑟𝑑𝑡 and 𝑑𝜃𝑑𝑡 are the radial and angular velocities inside the limiting sphere, respectively. Combining these two equations yields the ion trajectory upon entering the limiting sphere as: ( 𝑑𝑟𝑑𝜃 ) = ( 𝑟 𝑏 ) {1 − ( 𝑏𝑟 ) + 𝑒 𝑖𝑗 𝑣 𝜀 ( − )} (13) The condition 𝑑𝑟𝑑𝜃 = 0 yields a relationship between the minimum radial distance ( 𝑟 m ) that the two species under examination approach one another and their initial impact parameter: 𝑏 = 𝑟 m √1 + 𝑒 𝑖𝑗 𝑣 𝜀 ( m − ) (14a) Equation (14a) can be then used to define a critical impact parameter ( 𝑏 c ) for collision by equating 𝑟 m with the ion-ion collision radius 𝑎 𝑖𝑗 : 𝑏 c = 𝑎 𝑖𝑗 √1 + 𝑒 𝑖𝑗 𝑣 𝜀 ( 𝑖𝑗 − ) (14b) For all impact parameters less than 𝑏 c , corresponding to 𝑟 m < 𝑎 𝑖𝑗 , collision would occur (given Fuchs’s assumptions). Additionally the assumption that ions enter the limiting sphere uniformly distributed in initial 𝜃 yields: 𝑝 F = ( 𝑏 𝑐 δ ) (15) From equations (10), (14b) and (15), a final form for β δ is obtained as: β δ = 𝜋𝑎 𝑖𝑗2 ( b 𝑇𝜋𝑚 𝑖𝑗 ) 𝛾(𝑎 𝑖𝑗 ) , (16a) γ(𝑎 𝑖𝑗 ) = 1 + 𝑒 𝑘 b 𝑇 ( 𝑖𝑗 − ) . (16b) In γ(𝑎 𝑖𝑗 ) , 𝑘𝑇 is used in place of 𝑚 𝑖𝑗 𝑣 . This substitution ensures that the equation (8a) recombination coefficient converges to the ballistic (free molecular) limit expression correctly. We note that Fuchs utilized 𝑘𝑇 instead of 𝑘𝑇 , leading to incorrect ballistic limit (low pressure) two-body calculations. Equations (16a-b), with equation (8a), yield predictions for the ion-ion recombination coefficient with little difficulty in computation. However, the development of Fuchs’s approach requires several assumptions rendering it invalid in ion-ion recombination coefficient prediction, and its predictions can differ by an order of magnitude from literature reported ion-ion recombination rates. First, while there is no issue with the definition of a critical radius beyond which a continuum transport approximation is valid, it is not correct to assume that within ion-neutral gas molecule collisions negligibly affect their trajectories, i.e. it is not correct to completely neglect three-body trapping. This assumption was a major concern of Hoppel & Frick in the development of a limiting sphere theory including three-body trapping, as well as the concern of others in expanding upon their work. Second, as noted by Gopalakrishnan & Hogan but hitherto unaddressed in limiting sphere based approaches, on the limiting sphere boundary, incoming ions are not uniformly distributed in 𝜃 and their initial speed distribution function at the limiting sphere boundary will deviate significantly from a thermal equilibrium distribution. The change in velocity and angle distributions is attributed to the influence of potential interactions on ion motion prior to entering the limiting sphere; ions are attracted to one another and hence accelerated to much higher velocities in the direction of their center-to-center vector. To demonstrate this latter point, we note that in the Cartesian coordinate system depicted in Figure 1, wherein the “y” direction is aligned with the ion center-to-center vector, at the initial relative velocity vector between ions is given as 𝑣 ⃗⃗⃗⃗ = (𝑣 x , 𝑣 y , 𝑣 z ) = (𝑣 th , 𝑣 th + 𝑣 e , 𝑣 th ) , where 𝑣 th is sampled from the Maxwell-Boltzmann 1D distribution function: 𝑓(𝑣 th ) = √ 𝑚 𝑖𝑗 b 𝑇 𝑒𝑥𝑝 (− 𝑚 𝑖𝑗 𝑣 th2 b 𝑇 ) (17a) and 𝑣 e is the average electrophoretic velocity at 𝑣 e = (𝐷 𝑖 +𝐷 𝑗 )𝑘 b 𝑇 𝑒 δ (17b) Example plots of the joint probability density function (pdf, denoted as 𝜕 𝑛 ∗ 𝜕𝑣 𝜕𝜃 ) for the initial ion speed ( 𝑣 = ‖𝑣 ⃗⃗⃗⃗ ‖ = √𝑣 x2 + 𝑣 y2 + 𝑣 z2 ) and the initial angle 𝜃 = cos −1 ( 𝑣 y 𝑣 ) for the collision of NO and NH in He background gas at 300 K and variable pressures are shown in Figure 2. These pdfs were determined via randomly sampling 10 velocities per plot using equation (17a), rejecting all trajectories which would not enter the limiting sphere. The properties of NO and NH under these conditions are given in Table 1; these properties are used in all reported calculations. δ in Figure 2 is calculated via equation (9a) using a collision radius of 𝑎 𝑖𝑗 =3.35Å . For example calculations, the diffusion coefficients were calculated using the equations of Fuller et al; while this technically applies to the neutral vapor molecule, this approximation does not have a substantial influence on the presented results as the test gas (helium) is negligibly polarizable at 300 K. The test conditions yield a relative mean thermal speed between the two ions of 700 m s -1 . At the lowest pressures, where the limiting sphere radius is largest, potential interactions do not strongly influence the initial speed and angle distribution. However, as pressure increases, the decreasing value of δ leads to large increases in the mean speed, shifting the joint pdf to the right, and leading to preferential motion in the y direction. For comparison, lines denoting θ c are displayed, where θ c is defined from equation (14b) via the relationship 𝑏 c = δsinθ c . Table 1.
A summary of the ion properties utilized in calculations.
Nitrogen dioxide ion, NO Ammonium ion, NH Nominal Molecular mass [Da] 46 18 Diffusion Coefficient*Pressure [m /s∙Pa] 7.02 7.42 Electrical mobility*Pressure [m /Vs∙Pa] 273 289 Approximated Radius [Å] 1.775 1.575 Figure 2.
Heat maps showing the joint probability distribution functions (pdfs, 𝜕 𝑛 ∗ 𝜕𝑣 𝜕θ ) for the initial angle ( 𝜃 , as per Figure 1) and initial relative speed for ions entering a Fuchs limiting sphere, considering NO (entering ion) recombining with NH (center ion). In contrast to these heap maps, traditional limiting sphere calculation approaches in aerosol particle charging assume either the entering ion speed is at the mean thermal speed (relative) or Maxwell-Boltzmann distributed in speed, with a uniform angle distribution. The true joint pdf varies with pressure (upper right of each heat map) because of how the limiting sphere radius varies with pressure. θ c denotes the critical angle in the Fuchs model for collision as a function of velocity; for θ < θ c , collision occurs. C. Hoppel & Frick’s Probability Modification
To account for three-body trapping in limiting sphere calculations, Hoppel & Frick invoked the concept of a trapping sphere of radius χ ; if the center-to-center distance of the entities under examination becomes less than χ , collision would occur with probability 𝑝 χ . Redefining the limiting sphere radius as: δ = λ 𝑖𝑗 + χ , they proposed that 𝛽 δ (equation 10) can be expressed as: β δ = 𝜋χ ( b 𝑇𝜋𝑚 𝑖𝑗 ) γ(χ)𝑝 χ (18) where γ(χ) is calculated via equation (16b), replacing 𝑎 𝑖𝑗 with χ . Equation (18) and the equations presented subsequently would apply in instances where χ > 𝑎 𝑖𝑗 ; otherwise, Hoppel & Frick suggest use of equation (16a) using their modified δ definition. Implementation of their modification hence requires calculation of 𝑝 χ and χ . For a given χ , it can be shown that: 𝑝 χ = 1 − λ 𝑖𝑗2 [1 − 𝑒𝑥𝑝 ( −2χ𝑐𝑜𝑠θ c,χ λ 𝑖𝑗 ) (1 + 𝑖𝑗 cosθ c,χ )] (19a) θ c,χ = sin −1 ( 𝑏 c,χ χ ) (19b) 𝑏 𝑐,χ = 𝑎 𝑖𝑗 √1 + 𝑒 b 𝑇ε ( 𝑖𝑗 − ) (19c) Unfortunately, Hoppel & Frick do not provide a means to determine χ . As they were specifically concerned with particle-ion collision rate coefficients, they proposed χ can be determined using the ion-ion recombination rate coefficient itself at a given pressure and temperature. With the recombination rate ( β input ) known, they proposed that the theoretical expression of Natanson can be employed to determine an ion-ion trapping sphere radius, from which a separate particle-ion trapping sphere radius can be determined. In the present study, the ion-ion recombination rate is the parameter of interest, hence χ would be the ion-ion trapping sphere radius itself. According to Natanson , the input ion-ion recombination rate is linked to this trapping distance via the expression: β input = 𝜋χ ( ) 𝑓(𝑔)[1+ 𝑒2λ𝑖𝑗4𝜋ε0𝑘b𝑇𝜒(χ+λ𝑖𝑗) ]𝑒𝑥𝑝[ 𝑒24𝜋ε0𝑘b𝑇(χ+λ𝑖𝑗) ]1+ 𝜋ε0𝑘b𝑇χ2𝑓(𝑔)𝑒2λ𝑖𝑗 [1+ 𝑒2λ𝑖𝑗4𝜋𝜀0𝑘b𝑇χ(χ+λ𝑖𝑗) ]{𝑒𝑥𝑝[ 𝑒24𝜋ε0𝑘b𝑇(χ+λ𝑖𝑗) ]−1} (20a) where function f ( g ) is: 𝑓(𝑔) = 2ω − ω (20b) ω = 1 + 2 | 𝑒𝑥𝑝(−𝑔)𝑔 + 𝑒𝑥𝑝(−𝑔)𝑔 − | (20c) g = 𝑖𝑗 (20d) While it first glance implementation of equations (20a-d) to infer χ and use of equation (18) to yield β δ may seem tractable, the end result is a circular argument; a recombination rate coefficient is required as an input to determine the recombination rate coefficient. Furthermore, Hoppel & Frick do not discuss the consequence of changing the system pressure on the trapping distance. Pressure changes most certainly affect the ion-ion recombination rate
20, 26 , the value of λ 𝑖𝑗 , and likely change the most appropriate estimate of χ . However, if χ must be recalculated for changes in temperature and pressure, the ion-ion recombination rate would need to be recalculated in the first place to again determine the ion-ion trapping distance. For this reason we believe this approach cannot truly be employed used to predict the ion-ion recombination rate a priori , and is best thought of as a fitting procedure wherein a single value of χ is determined for a single β input (at a single temperature and pressure), and subsequently this value is used to predict the recombination rate coefficient or particle-ion collision rate coefficient for variable pressure and temperature conditions. We also remark that neither should a complete theory predicting the ion-ion recombination rate (or particle-ion recombination rate) require the trapping distance as an input nor is the trapping distance a physically rigorous parameter. The continuum – molecular dynamics hybrid approach obviates the need to discuss a hypothetical three-body trapping distance. D. Molecular Dynamics based Probability Calculation
We now return to equations (8b-d), which simply require a means of determining 𝑝 δ for implementation. To utilize MD simulations to compute 𝑝 δ we position ion i (the incoming ion) at a Cartesian coordinate location (0, ,0) and position ion j (the central ion) at (0,0,0). The equations of motion are then solved for both ions. The initial center-of-mass velocity vector of ion i is sampled from the joint probability density functions in Figure 2 (i.e. using equations 17a & 17b), while ion j is modeled at rest initially. The individual atoms in each ion have initial velocities which are the sum of the initial center-of-mass velocity and an additional thermal term, determined via separate NVT equilibration simulations at 300 K. We carried out calculations both including and omitting the influences of equation (17b), i.e. with and without the influence of electrostatic forces on the incoming ion initial velocity. We first discuss MD simulations carried out in the absence of neutral gas. Potential interactions ( ϕ inter ) between atoms within different ions are modeled via the Lennard-Jones 6-12 potential with the long range electrostatic potential: ϕ inter = ϕ LJ + ϕ E = 4ϵ 𝑘𝑚 [( σ 𝑘𝑚 𝑟 𝑘𝑚 ) − ( σ 𝑘𝑚 𝑟 𝑘𝑚 ) ] + 𝑧 𝑘 𝑧 𝑚 𝑒4𝜋ε 𝑟 𝑘𝑚 (21) where ϵ 𝑘𝑚 and σ 𝑘𝑚 are the Lenard-Jones potential parameters between atom k and atom m , 𝑧 𝑘 is the partial charge of atom k , and 𝑟 𝑘𝑚 is the atom-atom center-to-center distance. Lennard-Jones parameters for NH and NO were selected based on the AMBER force-field , and the partial charges were determined via NIH (National Institute of Health) database values. Intramolecular bond and angle potentials are additionally included, and were also calculated here based on the AMBER force-field. The AMBER force-field was simply selected as it is well-studied force-field; however, we remark that the influence of force-field choice on collision-controlled gas phase reactions is less studied than for molecular dynamics in solution, and this choice may need to be more carefully scrutinized in future continuum-molecular dynamics implementations. As depicted in Figure 1, calculations in the absence of background gas proceeded until one of three outcomes occurred: (1) the ion-ion separation distance was less than the collision radius a ij (collision); (2) the ion-ion separation distance exceeded non-collision); or (3) the ions were observed to enter stable orbit about one another. To compute 𝑝 δ we made sure to monitor at least 100 case (1) collision events, and 𝑝 δ was determined both treating orbits as non-collision events ( 𝑝 δ = 𝑁 col 𝑁 tot ⁄ , where 𝑁 col = 100 is the number of collision events and 𝑁 tot is the total number of trajectory cases examined) and treating orbits as collision events ( 𝑝 δ = (𝑁 col + 𝑁 orb ) 𝑁 tot ⁄ , where 𝑁 orb is the number of orbiting trajectories observed). Calculations neglecting neutral gas need only consider a small number of atoms (< 10), hence they proceed efficiently irrespective of the size of selected. However, equation (9a) calculated limiting sphere radii of 10 nm and higher can lead to instances where the number of gas molecules which would occupy the domain bounded by the limiting sphere exceeds 10 (and approaches 10 ). As 𝑝 δ is anticipated to decrease to well below unity (by orders of magnitude) with decreasing pressure, a greater number of trajectories must also be modeled for accurate calculation at low pressure; combined this makes MD simulations considering all gas molecules in the limiting sphere domain computationally intractable across a wide pressure range. To circumvent this issue while accounting for neutral gas, we instead opt to surround each ion with a cubic domain containing 10 neutral He gas atoms and whose dimensions is determined by the background gas temperature and pressure (hence simulating only 2000 additional atoms). All He atoms are initially positioned at random locations within their assigned cubic domain, with random velocity directional vectors and random speeds sampled from the Maxwell-Boltzmann speed distribution at the input temperature. Within each domain, He atom motion is monitored via the velocity-Verlet algorithm, with a number of simplifications. First, He-He interactions are completely ignored. This assumption has no bearing on ion trajectories except at the highest pressures examined, where the gas mean free path approaches the ion size, but as 𝑝 δ → 1 in the high pressure limit (the continuum limit), even at high pressure this approximation has no bearing on MD calculations (though the diffusion coefficients of ions would need to be modified in equation (8)). Second, He atoms interact with all atoms within the ion in their assigned box via Lenard-Jones potential interactions (i.e. equation 17a omitting the Coulombic term). Helium-NH + and Helium-NO potential parameters were estimated from the Lorentz-Berthelot combination rule , with ϵ He−He =0.0203 kcal/mol σ He−He =2.556 Å . The cubic domains are moved throughout trajectory calculations such that they remain centered on their respective ions. Importantly, He atoms do not interact with atoms in the ion not in their assigned box; therefore as ions approach one another and their domains intersect, the effective gas density around each ion remains constant. When a He atom leaves its assigned domain, periodic boundary conditions are invoked; while this leads to a modest increase in the total kinetic energy in each simulation (as the ions accelerate one another) we find that including 10 atoms leads to this increase negligibly influencing results. Third, because gas atoms only influence ion motion via collision (close range interaction), gas atoms which are a distance of 5.0 nm or more from their assigned ion’s center-of-mass are moved with a variable time step shown in the supporting information, while those within this radial distance are moved with the ion’s time step of 1 fs. Lennard-Jones potential interactions were not calculated for atoms more than 1.5 nm from one another. Combined, the three aforementioned model assumptions yield tractable 𝑝 δ calculations for ion-ion recombination at pressures as low as 10 kPa (and likely lower, but at greater computational expense). They additionally enable detailed accounting of the effect of ion-neutral gas collisions, obviating the need to assume either specular or diffuse gas atom scattering upon collision, i.e. all degrees of freedom in ions and gas atoms are modeled
48, 49 . This is an important feature in ion-ion recombination calculations, as the model required to accurately depict gas atom or molecule scattering from an ion without modeling ion internal degrees of freedom (frozen ion models) has been found strongly dependent on ion chemical composition and size.
50, 51
An example trajectory calculation including He gas atoms is depicted in Figure 3, with variable orientation views provided. Videos depicting trajectories resulting in collision, non-collision, and orbit are also including as supporting information files (where gas atoms are not shown). The line trace in Figure 3 depicts the trajectory of the incoming ion relative to the central ion (repositioned to remain fixed at the domain center). Qualitatively, the trajectory displays features which cannot result in the absence of background gas; the reversals in trajectory direction are the result of ion-neutral collisions. Simulations with gas atoms are also carried out until 10 collision events are observed; in the presence of neutral gas atoms, stable orbits are not observed. Calculations do require input values of and a ij , yet the calculation result needs to be independent of these choices; with neutral gas considered, needs to be sufficiently large, while a ij needs to be sufficiently small such that ions cannot escape one another at this distance. We probe the influence of both of these parameters, and report on these sensitivity analyses in the Results and Discussion section (section III). The reported calculations were carried out using hardware both at the Research Center for Computational Science (Okazaki, Japan) and the Minnesota Supercomputing Institute (MSI) of the University of Minnesota. A custom written C-code was employed for all calculations using trajectory calculation algorithms based upon the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) source code, which is open-source and optimized for efficient trajectory calculations. Figure 3.
Depiction of relative ion trajectories during a single molecular dynamics simulation (100 kPa). Each ion (center ion: NH , entering ion: NO ) is surrounded by 10 He atoms which do not interact with one another but fully interact with the atoms in the ion they are assigned to. The displayed instance does not result in collision. E. Non-Dimensionalization
As a final point in recombination rate coefficient method development, we are also interested in examining collapsed, dimensionless forms of the recombination rate coefficient. A series of recent studies have shown that the dimensionless collision rate coefficient ( 𝐻) in particle-particle,
39, 53, 54 particle-vapor molecule, and particle-ion
44, 56, 57 collisions can be expressed, almost universally, as a function the diffusive Knudsen number ( 𝐾𝑛 D ) , with these two parameters defined as: 𝐾𝑛 D = √𝑚 𝑖𝑗 𝑘 b 𝑇η c 𝑓 𝑖𝑗 𝑎 𝑖𝑗 η f (22a) 𝐻 = β 𝑖𝑗 𝑚 𝑖𝑗 η c 𝑓 𝑖𝑗 𝑎 𝑖𝑗3 η f2 (22b) where 𝑓 𝑖𝑗 , defined as ( D i + D j )/ k b T for two ions, is the effective reduced friction factor. η c and η f are the dimensionless continuum regime and free molecular enhancement factors, respectively, which take different functional forms based upon the potential interactions between the two colliding species under consideration. Neglecting shorter range interactions, for attractive Coulomb interactions between single charged ions alone c and f are respectively expressed as: η c = Ψ E E ) (23a) η f = 1 + Ψ E (23b) Ψ E = 𝑒 𝑘 b 𝑇𝑎 𝑖𝑗 (23c) Gopalakrishnan & Hogan developed a regression equation to Langevin dynamics simulations for hard sphere collisions, showing that: 𝐻 HS (𝐾𝑛 D ) = D2 +25.836𝐾𝑛 D3 +√8𝜋11.211𝐾𝑛 D4 D +7.211𝐾𝑛 D2 +11.211𝐾𝑛 D3 (24) where the subscript “HS” denotes the hard-sphere curves. Interestingly, via appropriate calculation of η c and η f , this curve has been found exendable to short range, attractive potentials (van der Waals and image potentials) and to repulsive Coulomb interactions.
13, 57
However, equation (24) yields incorrect results for attractive Coulomb interactions outside the 𝐾𝑛 D → 0 (continuum) and 𝐾𝑛 D → ∞ (free molecular) limits. Recently, Chahl & Gapalakrishnan developed improved regressions for attractive Coulomb potentials over a wide E and Kn D range (0< E <60, Kn D <2000): 𝐻(𝐾𝑛 D , Ψ E ) = 𝑒 𝜇(𝐾𝑛 D ,Ψ E ) 𝐻 HS (𝐾𝑛 D ) (25a) μ(𝐾𝑛 D , Ψ E ) = 𝐶𝐴 (1 + 𝑘 𝑙𝑛(𝐾𝑛 D )−𝐵𝐴 ) − −1 𝑒𝑥𝑝 [− (1 + 𝑘 𝑙𝑛(𝐾𝑛 D )−𝐵𝐴 ) − ] (25b) 𝐴 = 2.5 (25c)
𝐵 = 4.528𝑒𝑥𝑝(−1.088Ψ E ) + 0.7091𝑙𝑛(1 + 1.537Ψ E ) (25d) 𝐶 = 11.36Ψ
E0.272 − 10.33 (25e) 𝑘 = −0.003533Ψ E + 0.05971 (25f) We compare the dimensionless expression provided by equations (25a-f) to continuum-molecular dynamics determined recombination rate coefficients, non-dimensionalizing calculation results using equations (22a-b). III. Results & Discussion
A. Limiting Sphere Radius Selection
Figure 4.
Resulting probability of collision (a) and recombination rate coefficient (b) for variable values of the limiting sphere radius, , normalized by the Fuchs limiting sphere radius from equation (9a). The recombination rate coefficient is normalized by the hard-sphere continuum value. In utilizing equation (9a), Fuchs was attempting to identify a center-to-center distance beyond which continuum equations were approximately valid and within which free molecular trajectory calculations could be applied. The latter restriction does not apply in Fillipov’s approach and does not apply in continuum-molecular dynamics calculations; nonetheless, still needs to be sufficiently large for proper calculation. Near atmospheric pressure (100 kPa), we varied the choice of to examine its influence on calculations considering neutral gas and with the initial incoming ion velocity sampled as 𝑣 ⃗⃗⃗⃗ = (𝑣 th , 𝑣 th + 𝑣 e , 𝑣 th ) . As a function of the ratio / F , where F denotes the equation (9a) prediction (with 𝑎 𝑖𝑗 = 3.35Å ), Figure 4a displays the resulting 𝑝 𝛿 , while Figure 4b displays β 𝑖𝑗∗ , the recombination rate coefficient calculated with equation (8b), normalized by 𝑖𝑗 (𝐷 𝑖 + 𝐷 𝑗 ) . With increasing / F beyond 0.5, we find that 𝑝 𝛿 decreases in manner leading to a near constant value of β 𝑖𝑗∗ , suggesting that values of larger than ½ F are sufficient for calculations at atmospheric pressure. For limiting sphere radii selected smaller than ½ F the assumption of purely continuum transport beyond the limiting sphere radius appears to break down, leading to large changes in β 𝑖𝑗∗ . For the remainder of the calculations reported, based upon this sensitivity analysis at atmospheric pressure, we employ from equation (9a) with 𝑎 𝑖𝑗 = 3.35Å . B. The Recombination Rate Coefficient without Neutral Gas
Table S1 of the supporting information summarizes all continuum-molecular dynamics calculations performed. Prior to discussing calculation results with ion-neutral gas atom collisions, we first examine calculations in the absence of neutral gas, as in instances where the incoming ion velocity is sampled omitting electrostatic influences, results should agree with Fuchs’s theory. Figure 5(a) and Figure 5(c) display plots of 𝑝 𝛿 for initial incoming ion velocities accounting for both thermal and electrostatic effects, and only accounting for thermal effects, respectively. Figures 5(b) and 5(d) show the resulting β 𝑖𝑗 calculations for 𝑝 𝛿 results provided in Figures 5(a) and 5(c), respectively. The probability plots contain dashed lines for the values in the continuum limit ( 𝑝 𝛿 → 1 ) and free molecular limit. The latter is given by the expression: 𝑝 F = ( 𝑎 𝑖𝑗 δ ) (1 + Ψ E ) (26) We display results both discounting and including orbiting trajectories as collisions. When discounting orbiting trajectories, 𝑝 𝛿 approaches the minimum of the continuum and free molecular limiting curves under all conditions, irrespective of whether the incoming ion velocity is sampled including or excluding electrostatic potential influences. This suggests that the initial velocity has little influence on non-orbiting ion-ion relative trajectories, presumably because the electrostatic potential terms in the equations of motion rapidly increase ion velocities to the appropriate values in the absence of neutral gas. Interestingly, Fuchs calculations are in excellent agreement with our calculations when excluding orbiting trajectories, suggesting that the assumptions made in the original limiting sphere theory development are consistent with omitting the possibility of orbit. Including orbiting trajectories as collisions drastically changes calculation results outside the continuum limit such that the deviations from Fuchs predictions are substantial. Using calculation values at 100 kPa discounting and including orbiting trajectories to determine (smaller values correspond to discounting orbiting trajectories), the Hoppel & Frick approach can be used to fit calculation results with reasonable agreement near the selected pressure used in fitting, for the orbiting case. However, because the Hoppel & Frick model does not converge to the two-body free molecular limit as pressure decreases, it cannot capture the behavior of the orbit-excluded models, and it is a poorer fit outside the continuum limit when the electrostatic velocity is included as an initial condition. As stable orbits cannot be established in the presence of neutral gas (the ions will either eventually be pushed towards each other or pushed away from one another), the disparity in results including and discounting orbiting trajectories highlights the need to carry out full MD simulations within the limiting sphere. Figure 5.
The probability of collision, (a) & (c) , and recombination rate coefficient, (b) & (d) , resulting from molecular dynamics calculations in the absence of neutral gas atoms. Dashed lines denote the continuum and free molecular limits in all plots. (a) & (b) correspond to ions initiated with initial velocities and angles accounting for both thermal and electrostatic effects, while (c) & (d) correspond to ions initiated with thermal effects only. Closed symbols denote orbiting cases considered as collisions, while open symbols denote results considering orbit ions as non-collision events. In (b) & (d) predictions of Fuchs (equations 8a, & 16a-b, green lines), and of Hoppel & Frick, (equations 8a, 16a-b, and 18, blue lines) are shown. The Hoppel & Frick predictions are based upon the recombination rate coefficient at 100 kPa from simulations as an input, leading to the values shown in each plot; smaller values correspond to calculations discounting orbits in collisions. C. The Recombination Rate with Ion-Neutral Collisions
We compute 𝑝 δ using MD simulations by randomly sampling initial conditions for ions and gas molecules; we found that this is the most computationally efficient method to determine the collision probability, and fixing the number of collision events to be monitored at 10 sets the counting statistic uncertainty in 𝑝 δ at 10% of 𝑝 δ itself. However, a more rigorous approach would be to compute 𝑝 δ (𝑣 , θ ) , i.e. the collision probability for specific initial incoming ion speeds and angles. The value of 𝑝 𝛿 would then be computed as: 𝑝 δ = ∫ ∫ 𝑝 δ∞0𝜋/20 (𝑣 , θ ) 𝜕 𝑛 ∗ 𝜕𝑣 𝜕θ 𝑑𝑣 𝑑θ (27) While we do not employ equation (27), examination of 𝑝 δ (𝑣 , 𝜃 ) and the product 𝑝 δ (𝑣 , θ ) 𝜕 𝑛 ∗ 𝜕𝑣 𝜕θ calculated with and without neutral gas is of interest; this highlights how ion-neutral collisions influence ion trajectories. Figures 6 (a), (b), and (c) display plots of 𝑝 δ (𝑣 , θ ) calculated without neutral gas excluding orbit, without neutral gas including orbit, and with neutral gas, respectively, at 100 kPa. Such plots are similar to the collision probability curves of Yang et al. and Halonen et al . in examining impact parameters and initial speeds of molecules and clusters leading to condensation in the free molecular regime. Guidelines denote θ c , the critical angle for collision in Fuchs theory (below this curve, collision occurs, above it, ions do not collide) and multiples of θ c . Excluding orbit yields 𝑝 δ (𝑣 , θ ) in line with Fuchs definition of θ c , while including orbit shifts the 𝑝 δ (𝑣 , θ ) heat map slightly on a log-log scale. However, the inclusion of neutral gas molecules completely changes the 𝑝 δ (𝑣 , θ ) surface (we note the streaks in Figure 6c appear to be the result of spline fitting to a finite number of tested 𝑣 , θ locations). Including neutral gas, collisions are neither certain for any initial speed or angle, nor are they prohibited at any initial speed-angle combination except large angle, large velocity initial trajectories (which almost immediately leave the limiting sphere after entry). This highlights the need to explicitly include gas atoms in detailed modeling of ion trajectories; complete exclusion of gas atoms and approximation of ion-neutral collision influences via a three-body trapping distance both cannot recover the 𝑝 𝛿 (𝑣 , θ ) surfaces obtained from precise examination of ion-neutral collisions. Figures 6 (d), (e), and (f) display heat maps of the product 𝑝 δ (𝑣 , θ ) 𝜕 𝑛 ∗ 𝜕𝑣 𝜕θ for the conditions corresponding to 6(a), 6(b), and 6(c), respectively. Heat maps highlight the most probable collision type for the prescribed calculation conditions. Through comparison of Figure 6(f) to 6(d) and 6(e), it is evident that inclusion of neutral gas shifts the initial angles leading to most collisions to larger values than are predicted to be possible using Fuchs theory or anticipated based on the exclusion of neutral gas. Stated differentially, at atmospheric pressure most ion-ion collisions occur with ions initially in grazing collision trajectories which would not collide in the absence of gas. Figure 6.
The probability of collision as a function of and v at 100 kPa, for simulations without neutral gas atoms and ignoring orbiting events (a) , without neutral gas atoms but counting orbiting as collision (b) , and including ion-neutral collisions (c) . For the three examined cases, the resulting products of the probability of collision and joint probability density function for the initial angle and velocity are shown in (d) , (e) , and (f) , respectively. Analogous to Figure 5, in Figure 7 we plot the probabilities of collision including neutral gas (Figure 7a) and the resulting recombination rate coefficient (in Figure 7b). Results in the absence of neutral gas are included for comparison, along with continuum and free molecular limiting lines. Gas inclusion increases the overall probability of collision, leading to recombination coefficients higher by an order of magnitude than the orbit excluded instances outside the continuum regime, and higher by a factor of ~2 than the orbit included instances outside the continuum regime. Including neutral gas near atmospheric pressure, we arrive at a recombination rate coefficient of 1.1 x 10 -12 m s -1 including electrostatic effects in the initial velocity. Lee & Johnsen report a recombination rate coefficient near 1.0 x 10 -12 m s -1 in He near atmospheric pressure, which is in outstanding agreement with calculations. Near a pressure of 30 kPa they report a recombination rate coefficient near 7 x 10 -13 m s -1 ; this is again in excellent agreement with our calculations (6.7 x 10 -13 m s -1 ) as pressure is reduced below atmospheric levels. An increase in recombination rate with increasing pressure at low pressure, and then its decrease with increasing pressure at high pressure has been documented in some of the earliest recombination rate experiments,
26, 60 and it is this unique behavior that served as driver for many of the initial theoretical investigations of ion-ion recombination. Our calculations capture this phenomenon; calculations outside the continuum limit reveal a monotonically increasing recombination rate coefficient with increasing pressure and at the continuum limit, a decreasing recombination rate coefficient with increasing pressure. In total, we find that the explicit inclusion of neutral gas in calculations has a substantial influence on the recombination rate coefficient calculation, and that without any fit parameters (i.e. using input for the test ions and neutral gas inferred from their tabulated properties, without adjustment), we can achieve excellent agreement between continuum-molecular dynamics hybrid calculations and experimental measurements. Figure 7.
The probability of collision (a) and recombination rate coefficient (b) , both considering neutral gas (triangles) and neglecting neutral gas (circles, squares).
D. Non-Dimensional Curves
A final pending question is the extent to which calculations agree with regression based dimensionless curves resulting from Langevin simulations. However, to non-dimensionalize results appropriately, the choice of the collision radius a ij needs to be carefully scrutinized, as neither ion is strictly spherical, and the collision radius, i.e. the center-to-center distance within which charge neutralization is certain, is not easily defined. For 100 kPa calculations including neutral gas and electrostatic potential effects on the entering ion velocity, in Figure 8 we plot the recombination rate coefficient as a function of the input collision radius. Alterations to the collision radius do not require new calculations; the same set of trajectories can be re-analyzed with a different assumed a ij value. The true collision radius will be the largest distance below which the collision radius assumed has no bearing on the inferred recombination rate coefficient. In Figure 8 the input value of is found to be far below the true collision radius of 12.8 Å. We thus elect to use 12.8 Å as the collision radius in non-dimensionlization. Figure 8.
The recombination rate coefficient as a function of the input collision radius (center-of-mass to center-of-mass distance considered to be collision) at 100 kPa. Figure 9 displays a plot of the dimensionless rate coefficient H from equation (22b) as a function of Kn D from equation (22a) from calculations including neutral gas and using the inferred collision radius of 12.8 Å in non-dimensionalization. Also plotted are the continuum ( 𝐻 = 4𝜋𝐾𝑛 D2 ) and free molecular ( 𝐻 = √8𝜋𝐾𝑛 D ) limit expressions, non-dimensionalized Fuchs (1963) predictions, the Chahl & Gopalakrishnan regression equation, and a non-dimensionalized Hoppel & Frick fit where = 21.6 nm is based upon the recombination rate coefficient at 100 kPa. In large part, both the Chahl & Gopalakrishnan regression equation and the Hoppel& Frick fit capture our calculations well at Kn D < 100. Deviations begin to manifest at Kn D > 100 and are significant at Kn D > 1000. This high Kn D limit is near the upper bound of the Chahl & Gopalakrishnan regression equation, and as discussed previously, the Hoppel & Frick model does not converge appropriately to the two body collision limit, which will eventually be reached at low enough pressure, i.e. sufficiently high diffusive Knudsen number. While dimensionlessly, our results apply only to a specific value of Ψ E = 43.4 , results do suggest that with improvements, the dimensionless collision rate coefficient curves developed for particle collisions can be extended to ion-ion recombination, which vastly simplifies recombination rate estimation. At the same time, when higher precision calculations are needed or the influences of ion structure are of interest, continuum-molecular dynamics hybrid calculations appear to be a tractable method to study ion-ion collision phenomena in gases and aerosols. Figure 9.
The dimensionless recombination rate coefficient as a function of the diffusion Knudsen number. All dimensionless ratios are defined using a collision distance of 12.8 Å. Predictions based on Fuchs , and Chahl & Gopalakrishnan (equation 25), and fitting using Hoppel & Frick are also plotted. IV. Conclusions
We develop a continuum-molecular dynamics hybrid calculation approach for the ion-ion recombination rate coefficient, which is based upon the limiting sphere approach of Fuchs and Fillipov (1993). Recombination rate coefficient calculations explicitly modeling ion-neutral gas atom (or molecule) collisions are made possible by including two distinct cubic simulation domains around ions whose dimensions are determined by the system temperature and pressure, while the number of gas atoms within each domain is set at 10 . Neglect of gas atom-gas atom interactions, variable time stepping, and ensuring that gas atoms only interact with their assigned ion all increase computational speed, which is necessary to simulate the required number of trajectories to determine collision probabilities and recombination rate coefficients. Including all degrees of freedom in ions and inclusion of the neutral gas enables accurate modeling of the outcomes of ion-neutral gas close encounters. We find that ion-neutral collisions fundamentally change the nature of the recombination process; recombining ions can approach one another at grazing angles, in comparison to traditional theories which largely predict “head-on” collisions. For this reason, we advocate adoption and implementation of continuum-molecular dynamics hybrid approaches in lieu of traditional limiting sphere theories where ion motion within the limiting sphere is imprecisely modeled. Continuum-molecular dynamics calculation predictions, without any fitting, are found to yield recombination rate coefficients in excellent agreement with prior measurements in He gas near atmospheric pressure. Though the test case examined was limited to Helium, by modeling N and O , the approach is also readily extendable to ion-ion recombination in atmospherically relevant conditions and for ion compositions relevant to aerosol charging. Such extensions will require inclusion of induced-dipole potentials between these diatomic gas molecules and ions. Condensable vapor molecules, i.e. water molecules, can also be added to calculations to examine their influence both on ion structure and on the recombination rate coefficient. Furthermore, continuum-molecular dynamics approaches are extendable to examine clustering reactions, condensation, and coagulation; for these processes they can be used for accurate rate calculations at variable temperature and pressure, provided the continuum transport properties of colliding species are known, and provided suitable potentials are developed (either all-atom or coarse-grained) to represent colliding species. As a finally note, beyond the collision rate coefficient, of interest in collisions in clusters is the structure and thermodynamic properties of the collision product, which can be out of thermal equilibrium with its surroundings and for this reason may be anomalously reactive until cooled by subsequent collisions with the background neutral gas. As Coulombically driven collisions occur at substantially elevated speeds, future implementations of continuum-molecular dynamics simulations may afford the opportunity to examine unique Coulombically facilitated chemical reactions for ions, clusters, and charged particles. Supporting Information.
Table S1 displays the Pressure, gas simulation conditions, resulting probability, recombination rate coefficient, diffusive Knudsen number, and dimensionless rate coefficient from all simulations performed. Information on the gas atom time steps in simulations are also provided, as are video depictions of ion-ion relative trajectories colliding, not colliding, and orbiting, in the absence of neutral gas. Acknowledgements.
C.J.H. acknowledges support from Department of Energy Award No. DE-SC0018202. Computations were performed using computing resources at the Research Center for Computational Science, Okazaki, Japan as well as the Minnesota Supercomputing Institute (MSI) at the University of Minnesota. T. T. was supported by the Hosokawa Powder Technology Foundation. H. H. and T. S. were supported by JST CREST Grant Number JPMJCR18H4, Japan. References K. A. Brueckner, The Journal of Chemical Physics (1964) 439. J. J. Thomson, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science (1924) 337. M. R. Flannery, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences (1982) 447. L. B. Loeb, and L. C. Marshall, Journal of the Franklin Institute (1929) 371. A. Franchin et al. , Atmos. Chem. Phys. (2015) 7203. D. Smith, and M. J. Church, Planetary and Space Science (1977) 433. H. F. Calcote, Symposium (International) on Combustion (1961) 184. B. Y. H. Liu, and D. Y. H. Pui, J Aerosol Sci (1974) 465. I. Ibarra, J. Rodríguez-Maroto, and M. Alonso, J Aerosol Sci (2020) 105479. R. P. Turco, J.-X. Zhao, and F. Yu, Geophysical Research Letters (1998) 635. Y. Wang et al. , Aerosol Science and Technology (2017) 833. G. Sharma et al. , J Aerosol Sci (2019) 70. H. Ouyang, R. Gopalakrishnan, and C. J. Hogan, The Journal of Chemical Physics (2012) 064316. A. Wiedensohler, J Aerosol Sci (1988) 387. L. Li, H. S. Chahl, and R. Gopalakrishnan, J Aerosol Sci (2020) 105481. H. M. Mott-Smith, and I. Langmuir, Physical Review (1926) 727. J. J. Thomson, and G. P. Thomson,
Conduction of Electricity through Gases (Cambridge University Press, London, 1928), P. Langevin, Annales de chimie et de physique (1903) 433. S. McGowan, Physics in Medicine and Biology (1965) 25. L. B. Loeb,
Fundamental Processes of Electrical Discharge in Gases (John Wikey & Sons, New York, 1939), M. E. Gardner, Physical Review (1938) 75. W. Gringel, K. H. Käselau, and R. Mühleisen, Pure and Applied Geophysics (1978) 1101. M. A. Biondi, Canadian Journal of Chemistry (1969) 1711. D. R. Bates, Planetary and Space Science (1982) 1275. H. S. Lee, and R. Johnsen, The Journal of Chemical Physics (1989) 6328. L. B. Loeb, Physical Review (1937) 1110. J. E. Allen, Phys. Scr. (1992) 497. G. Natanson, Soviet Physics Technical Physics (1960) 538. D. R. Bates, and M. R. Flannery, Journal of Physics B: Atomic and Molecular Physics (1969) 184. J. N. Bardsley, and J. M. Wadehara, Chemical Physics Letters (1980) 477. P. J. Feibelman, The Journal of Chemical Physics (1965) 2462. E. K. Parks, The Journal of Chemical Physics (1968) 1483. C. Larriba, and C. J. Hogan, The Journal of Physical Chemistry A (2013) 3887. C. Larriba, and C. J. Hogan, Journal of Computational Physics (2013) 344. X. B. Nie et al. , Journal of Fluid Mechanics (2004) 55. N. A. Fuchs, Geofis. Pura Appl. (1963) 185. A. V. Filippov, J Aerosol Sci (1993) 423. W. A. Hoppel, and G. M. Frick, Aerosol Science and Technology (1986) 1. R. Gopalakrishnan, and C. J. Hogan, Aerosol Science and Technology (2011) 1499. P. G. Wright, Discussions of the Faraday Society (1960) 100. X. López-Yglesias, and R. C. Flagan, Aerosol Science and Technology (2013) 688. X. López-Yglesias, and R. C. Flagan, Aerosol Science and Technology (2015) 196. Y. G. Stommel, and U. Riebel, Aerosol Science and Technology (2007) 840. R. Gopalakrishnan, and C. J. Hogan, Physical Review E (2012) 026410. E. N. Fuller, P. D. Schettler, and J. C. Giddings, Industrial & Engineering Chemistry (1966) 18. W. D. Cornell et al. , Journal of the American Chemical Society (1995) 5179. T. A. Halgren, Journal of the American Chemical Society (1992) 7827. T. Tamadate et al. , Aerosol Science and Technology (2019) 260. R. Lai, E. D. Dodds, and H. Li, The Journal of Chemical Physics (2018) 064109. C. Larriba-Andaluz et al. , Physical Chemistry Chemical Physics (2015) 15019. H. Ouyang et al. , Journal of the American Society for Mass Spectrometry (2013) 1833. S. Plimpton, Journal of Computational Physics (1995) 1. T. Thajudeen, R. Gopalakrishnan, and C. J. Hogan, Aerosol Science and Technology (2012) 1174. A. M. Boies et al. , Small (2019) 1900520. R. Gopalakrishnan, T. Thajudeen, and C. J. Hogan, Journal of Chemical Physics (2011) 054302. H. S. Chahl, and R. Gopalakrishnan, Aerosol Science and Technology (2019) 933. R. Gopalakrishnan et al. , J Aerosol Sci (2013) 60. H. Yang, E. Goudeli, and C. J. Hogan, The Journal of Chemical Physics (2018) 164304. R. Halonen et al. , Atmos. Chem. Phys. (2019) 13355. W. Mächler, Zeitschrift für Physik (1937) 1. A. Maißer et al. , J Aerosol Sci (2015) 36. L. Ahonen et al. , The Journal of Physical Chemistry Letters (2019) 1935. H. Yang, Y. Drossinos, and C. J. Hogan, The Journal of Chemical Physics (2019) 224304.(2019) 224304.