The mechanics and dynamics of molecular collisions in air are investigated by thoroughly validated atomistic molecular dynamics (MD) simulations that treat oxygen and nitrogen as true diatomic molecules accounting for their non-spherical shape and, most importantly, force field. Due to their rotational motion and non-spherical shape, molecules follow complex trajectories at close enough separations experiencing a great variety of collision events. Most of the collisions are bimolecular. However, some can involve up to four molecules as pairs (or even triplets) of molecules that collide repeatedly are observed. Following their initial encounter, these molecules separate briefly, come back, and collide again and again creating even “orbiting” collisions, before they split apart to collide with other molecules. Identifying such rather spurious collisions and filtering them by hazard plot analysis was a key step in correctly determining collision densities and accumulating collision event statistics. By systematically recording the distribution of free paths (distances traveled by molecules between genuine collisions), the mean free path, λ, of air is determined as 38.5 ± 1 nm at 300 K and 1 atm. This is 43% smaller than the 67.3 nm widely accepted λ today at these conditions and quite robust to the employed MD force field as long as it accurately matches the experimentally determined macroscopic properties of air (density, viscosity, and diffusivity).

## I. INTRODUCTION

Molecular collisions determine the key transport properties of gases while their mean free path along with particle size largely control aerosol behavior in the environment and several industrial processes.^{1} From a theoretical point of view, an elegant treatment of molecular collisions and transport in gases is offered by the kinetic theory of gases,^{2} which leads to closed-form analytical expressions for their diffusivity, *D*, viscosity, *η*, and thermal conductivity, *κ*.^{3} For example, starting from the Maxwell–Boltzmann distribution of molecular velocities and the fundamental conservation laws for the total mass, energy (thus elastic collisions), and momentum assuming that the interaction between any two molecules depends only on their separation *r*, a closed-form expression for the angle of deflection *χ* is derived. This involves the inter-molecular potential *u*(*r*), the magnitude *ν _{r}* of the initial relative velocity of the two molecules, the distance of closest approach

*r*, and the impact parameter

_{m}*b*of the collision (or

*r*in the absence of molecular interactions).

_{m}^{3}

At a higher level of analysis,^{3} one solves the non-linear integrodifferential Boltzmann equation^{4} with the Chapman–Enskog method^{5,6} invoking a series expansion in terms of an ordering parameter.^{3} This helps express the transport coefficients in terms of weighted integrals (the so-called collision integrals) of the deflection angle *χ* carrying information about the dynamics of bimolecular collisions, hence about the intermolecular potential.^{3}

For the hard-sphere model, the collision integrals can be evaluated analytically. The same is true for extensions^{7–11} of the hard-sphere model, such as the variable hard-sphere (VHS) model^{7} and the variable soft-sphere (VSS) model,^{8,9} wherein a strong dependence of the total cross-section on the relative velocity of the collision pair is incorporated in order to predict the correct temperature dependence of the viscosity coefficient.^{12,13} For more sophisticated potential functions, however, such as the Lennard–Jones, the exp-6 (Ref. 14) and the Kihara potential,^{14} numerical integration is needed. Hirschfelder *et al.*^{15} have studied the transport coefficients of 6-12 Lennard–Jones (LJ) gases in great detail, with the two LJ parameters *σ* and *ε* obtained by fitting data for the viscosity and the second virial coefficient of the gas as a function of temperature. For spherical molecules described by a realistic potential (i.e., Lennard–Jones), Hirschfelder *et al.*^{15} were the first to examine the anatomy of their collisions. They showed that grazing, head-on, and even orbiting collisions can take place due to the intermolecular interaction through their Lennard–Jones potential.

Despite these insights, unfortunately, the assumption of structureless particles undergoing elastic collisions (rigid sphere theory) is widely invoked. This considerably limits the consistency and applicability of the estimated parameters. For example, LJ parameters obtained by fitting strictly viscosity data may not agree with LJ parameters obtained by fitting second virial coefficient data.^{15} Likewise, by keeping the LJ parameters *ε* and *σ* constant, it may not be possible to obtain satisfactory agreement of the transport properties over an extended temperature range.^{15}

To overcome these limitations, a more realistic treatment of the mechanics of collisions is needed where one accurately accounts for the diatomic nature of oxygen and nitrogen molecules and their interactions at the level of atoms (and not at the level of molecules).^{15} An excellent such alternative is offered by molecular dynamics (MD), which can provide an exact solution to the corresponding statistico-mechanical problem over approximate solutions. As MD is a deterministic methodology, it allows the validation of force field performance against macroscopic observables (i.e., for air, density, diffusivity, and viscosity) and track exact trajectories of molecules and their constituent atoms. Today, the MD method has evolved considerably due to a combination of advances such as development of more accurate parameterizations of experimentally verified force fields,^{16} efficient algorithms^{17} often combined with machine learning,^{18} and a steadily increasing computational power from hardware advances.^{19} As a result, MD is a powerful tool for unraveling basic phenomena that dictate the dynamics and evolving structure of nanoparticles and their intriguing performance that impacts several applications (see Sec. S1 in the supplementary material).^{20–28}

Here, MD is used to rigorously characterize molecular collisions of gases with a thoroughly-validated fully atomistic (FA) model incorporating bonded and non-bonded interactions and identify distinct patterns of collision dynamics. Also derived from the molecular collisions is the distance traveled by molecules between collisions, the so-called mean free path, *λ*.^{29,30} This is a fundamental quantity in a variety of fields, most notably aerosol science, where it is widely used to distinguish particle motion in the free molecular, transient, or continuum regimes.^{31} To compute *λ*, the simple kinetic theory (wherein molecules are taken to be spherical and rigid with a diameter *σ*) assumes that a collision between molecules *i* and *j* occurs whenever the impact parameter *b* is shorter than $ \sigma i j = ( \sigma i + \sigma j ) / 2$. In reality (and our MD simulations), however, gas molecules are diatomic and non-spherical, exhibiting several modes of interaction; thus, the above simple definition strictly does not apply.

To appreciate the results of the FA simulations in comparison with the current state of the art, two limiting cases are explored also: In the first, air molecules are approximated as hard spheres (HS model) undergoing strictly elastic collisions; this enables a direct comparison with classic kinetic theory for which analytical results for the collision densities and the mean free path are available. In the second, such spheres are empowered by a Lennard–Jones force field (LJ model), as Hirschfelder *et al.*^{15} first did, to bridge here the gap between the simplified HS model and the detailed FA one; this helps examine the importance of molecular shape (of relevance to polyatomic molecules) on top of that of the precise potential energy function. Last but not least, the mean free path is calculated also using various force fields and MD models to further assess the robustness and accuracy of the *λ* obtained here.

## II. METHODS

### A. Molecular dynamics (MD)

To precisely calculate molecular collisions and investigate their mechanics, an accurate force field is needed. In the literature, several atomistic models have been derived for air treating nitrogen and oxygen as diatomic molecules, including intramolecular bond stretching interactions (usually described by harmonic potentials) and intermolecular repulsive and attractive interactions (represented by LJ functions). The force field parameters are typically extracted by fitting simulation predictions to various thermodynamic observables (including, but not limited to, liquid density, Raman vibrational frequency, enthalpy of vaporization, solvation free energy, second virial coefficient, and viscosity). The force field of Zambrano *et al.*^{32} is most reliable for air as it combines the parameters of Laufer and Leroi^{33} for oxygen (based on solid-state thermodynamic data) and Murthy *et al.*^{34} for nitrogen [obtained by fitting both solid-state (lattice spacing and sublimation energy of the crystal at absolute zero) and second virial coefficient data]. Recently, Wang *et al.*^{35} also proposed a force field for air that reproduced well bond geometry, Raman vibration frequency, liquid density, and vaporization enthalpy. Both of these models are two-site potentials. Several three-site potentials, including electrostatic interactions, have also appeared.^{36–39} The third atom is massless, resides at the molecule's geometric center, and serves to generate a quadruple moment term in the potential energy function. Nevertheless, and as often noted in the literature,^{40–43} the quadruple moment term is essential to include in simulations of condensed phases or for addressing adsorption phenomena; otherwise, it can be safely neglected.

To check the sensitivity of the obtained simulation results on the precise value of the parameters describing non-bonded interactions, the transport properties and mean free path of air are compared from all the above-mentioned fully atomistic force fields, either 2- or 3-site ones,^{32,35–39} omitting, however, the phantom third atom of the three-site models.^{36–39}

It is worth noting that *ab initio* force fields^{44–47} address air properties at extremely high temperatures (above 1000 K) by accounting for chemical reactions (i.e., bond dissociation).^{48–51} Such force fields are described by complex interaction potentials, including many-body terms, and characterized by tens or even thousands of parameters. However, their examination goes beyond the scope of the present study as reactions are not considered here. On the other hand, there are two classic two-body potentials for pure nitrogen^{41,43} parameterized on the basis of such *ab initio* calculations.^{47,52} As above, the performance of the employed force field here,^{32} with respect to nitrogen density, viscosity, diffusivity, and *λ*, is also compared to these two potentials,^{41,43} as well as to available experimental measurements, to fully assess the robustness of its predictions.

All MD simulations have been performed with a force field^{32} incorporating bonded (i.e., bond stretching) and non-bonded interactions (Table S1 of the supplementary material). This is the so-called FA model. The nitrogen–nitrogen and oxygen–oxygen bonds are considered to be rigid with the corresponding holonomic constraints addressed by the SHAKE algorithm,^{53} while the rest of intermolecular contributions are described by a 12-6 Lennard–Jones (LJ) potential. The pairwise LJ interactions are truncated at a cutoff distance of 12 Å and to compensate for this, tail corrections are applied for both energy and pressure.^{54} The simulations are carried out primarily in the microcanonical (NVE) ensemble, while canonical (NVT) and isothermal-isobaric (NPT) ones are examined also with LAMMPS.^{55} The Nosé–Hoover thermostat^{56,57} and the Parrinello–Rahman barostat^{58} were used to keep *T* and *P* at their prescribed values with the corresponding relaxation constants at 50 and 400 fs, respectively. The equations of motion were integrated using the velocity-Verlet method^{59} with a time step of 1 fs. Large cubic simulation cells (of length 126.7 nm, subject to periodic boundary conditions in all three space directions) were used containing 5 × 10^{4} air molecules (consisting of 79 mol. % N_{2} and 21 mol. % O_{2} and corresponding to air concentration of about 0.041 M) to reduce statistical uncertainty and ensure the complete absence of system-size effects (see Sec. S5 of the supplementary material).

As already mentioned in the Introduction, additional simulations were performed in limiting cases wherein O_{2} and N_{2} were assumed to be hard spheres (HS model) or Lennard–Jones spheres (LJ model) to directly compare the results from the present FA model to those from the kinetic theory of gases but also assess the significance of each one of its key assumptions: elastic collisions and spherical molecule shape.

^{3}

^{,}

*r*denotes the separation distance. To avoid the discontinuity at $ r = \sigma i j$, the HS potential is approximated with a “softer” version for $ r \u2264 \sigma i j$,

^{60}

^{,}

*A*= 5000 kcal/mol, an energy prefactor.

^{60}The van der Waals diameters are

^{15}

*σ*

_{1}= 3.467 Å for oxygen and

*σ*

_{2}= 3.798 Å for nitrogen (

*i*= 1 and 2 denote oxygen and nitrogen molecules, respectively).

For the LJ model that “bridges” the HS and FA models, the following parameters were used:^{15}^{,} *ε*_{1} = 0.212 03 kcal/mol and the above *σ*_{1} for oxygen, as well as *ε*_{2} = 0.141 98 kcal/mol and the above *σ*_{2} for nitrogen. For cross interactions, the Lorentz–Berthelot mixing rules^{54} were applied for *ε*_{12} = 0.173 505 kcal/mol and *σ*_{12} = 3.6325 Å. All parameters of the FA force field are summarized in Table S1 of the supplementary material.

### B. Collision criterion

This is based on the calculated interatomic separation between any pair of molecules (or their constituent atoms) where the corresponding inter-molecular (or inter-atomic) force crosses from negative to positive values, meaning that one molecule (or atom) starts exerting a repulsive force to the other. For the HS model, this threshold distance *r _{m}* is equal to the van der Waals sum of radii of the two molecules,

*σ*. For the LJ and FA models, it is equal

_{ij}^{54}to 2

^{1/6}

*σ*, where

_{ij}*i*and

*j*refer to entire nitrogen or oxygen molecules (for the LJ model) or their atoms (for the FA model). The

*r*for the three models are shown in Table S2 of the supplementary material. To save computational time, first, candidate pairs of molecules or atoms for collision are identified by monitoring changes (even as small as a 1° angle for possible grazing collisions

_{m}^{15}) in the time autocorrelation function of the unit momentum vector for all air molecules. The next step is to check if the collision is true or spurious, i.e., if it involves successive collisions between the same molecules. To precisely distinguish and filter spurious collisions, we used hazard plot analysis,

^{61}as explained immediately below.

### C. Detection of spurious collisions

*λ*from its last collision experiences a new collision after having traveled length

*λ*+ d

*λ*. The cumulative probability (cumulative hazard) will be $ H ( \lambda ) = \u222b 0 \lambda h ( \lambda \u2032 ) d \lambda \u2032$. Then, the expectation of the cumulative hazard at the length of the

*m*th ordered collision (out of

*n*collisions) is

^{62}

^{,}

*ξ*), the cumulative hazard is a straight line

^{62}(with a slope equal to 1/

*ξ*). Events deviating from the straight line correspond, therefore, to correlated collisions or to collisions not obeying the exponential distribution.

### D. Mean free path, λ

To calculate *λ* from the NVE MD simulation trajectories (in the NPT simulations, the measured mean free path was not accurate due to the presence of the barostat that included contributions also from the volume fluctuations of the simulation cell; see Selection of Ensemble in Sec. S2 of the supplementary material), first, all collisions and all lengths traveled by molecules are recorded within a pre-specified time interval (the so-called observation time, *t*_{obs}), and the corresponding probability density distribution is computed. The *λ* is calculated by directly averaging over this distribution or by regressing the latter to an exponential function (as suggested by kinetic theory). Then, the observation time is increased and the calculations are repeated. This helps sample longer free paths, enrich their distribution, and thus update the value of *λ*. The procedure is repeated until convergence is attained in the sense that although the distribution of free paths continues to change slightly as longer and longer free paths are sampled and added, the *λ* approaches asymptotically a constant value. This is also calculated by regressing the *λ* simulation predictions with a hyperbolic function of the form $ \lambda ( t obs ) = ( b 1 + b 2 t obs ) / ( 1 + a \u2009 t obs )$ with *a*, *b*_{1}, and *b*_{2} numerical constants, and then taking the limit *t*_{obs} → ∞ to get $ \lambda = lim t obs \u2192 \u221e \lambda ( t obs ) = b 2 / a$. In the following, the so-extracted *λ* values are compared with results from theoretical expressions based on the kinetic theory of gasses (hard-sphere model) and extensions thereof such as the variable hard-sphere^{7} (VHS) and variable soft-sphere^{8,9} (VSS) models. A brief account of these theories, which all consider air molecules as perfect elastic spheres, is given in Sec. S3 of the supplementary material.^{12,13}

## ΙΙΙ. RESULTS AND DISCUSSION

### A. Validation

All simulations were carried out with rather large simulation cells at temperature *T* of 300 K and pressure *P* of 1 atm in the microcanonical (NVE) ensemble (see Selection of Ensemble in Sec. S2 of the supplementary material). The air velocity distributions were identical to those by the classic Maxwell–Boltzmann expression^{29} as shown in Fig. S2. The employed atomistic force field was validated by comparing the MD-predicted air density (1.180 ± 0.01 kg/m^{3}) to the experimentally measured one for dry air (1.177 kg/m^{3}).^{63} A second validation involved a comparison of the MD-predicted diffusivity, *D*, and viscosity, *η* of air with measurements.^{64–66} The *D* of air was predicted from the slope of the mean square displacement of O_{2} and N_{2} molecules in the Fickian (linear) regime.^{64} For the viscosity, *η*, the corresponding Green–Kubo equation^{54} was used involving the integral of the time autocorrelation function of the off diagonal components of the atomic stress tensor. The FA MD-predicted *D* and *η* were 0.203 ± 0.005 cm^{2}/s and 18.1 ± 0.2 *μ*Pas, respectively, both in excellent agreement with experimental^{65,66} or theoretical studies^{64,67} giving *D* = 0.203 cm^{2}/s and *η* =18.5 *μ*Pas, respectively.

The predictions for the *ρ*, *D*, and *η* of air from the two limiting cases of the FA model were, respectively, 1.17 ± 0.01 kg/m^{3}, 0.185 ± 0.005 cm^{2}/s, and 18.0 ± 0.5 *μ*Pas from the HS model, and 1.18 ± 0.01 kg/m^{3}, 0.2 ± 0.005 cm^{2}/s, and 18.1 ± 0.5 *μ*Pas from the LJ model. Overall, the three models provide accurate predictions of air density, diffusivity, and viscosity, except for the diffusivity from the HS model, which is underestimated by about 9%. For the FA model, these excellent predictions should have been expected. For the LJ model, this is the direct result of its parameterization on the basis of experimental data for the viscosity.

### B. Molecular collisions: Definition and distinction from spurious ones

The accumulated MD trajectories were analyzed to identify molecular collisions based on the distance between colliding molecules.^{3} A collision is declared when pairs of molecules or atoms have interatomic distance smaller than the collision distance *r _{m}* (Table S2 of the supplementary material). Figure 1 shows representative collision trajectories (bold black lines) from the three models in relative coordinates (i.e., one of the two molecules is kept at rest). Figure 1(a) (Multimedia view) shows the impact parameter

*b*and the angle

*θ*at the distance of the closest approach according to the HS model. Once a molecule hits another one, it changes direction and moves along a different straight line. The angle

_{m}*θ*is related to the deflection angle through $ \chi + 2 \theta m = \pi $; for the collision in Fig. 1(a),

_{m}*θ*= 38.4°. In general, when a collision occurs according to the HS model, the relative velocity along the line of impact is reversed, while that at right angle to this line remains unchanged.

_{m}^{3}

Figure 1(b) (Multimedia view) displays such a collision according to the LJ model. This is largely similar to that of the HS model except just before contact, where the LJ potential between molecules leads to a slight bending of the trajectory as these are attracted to each other, increasing the probability of their collision. Similarly, Fig. 1(c) (Multimedia view) displays a collision according to the FA model also depicting the constituent atoms of the molecules. The collision trajectory in Fig. 1(c) is not that different from that by the simpler LJ model, even though the shape difference (dumbbell vs spherical) enhances, even more, the collision probability. As A and B molecules approach each other at *θ _{m}* = 60.1° [Fig. 1(c)], they come even closer by attractive van der Waals forces, bending the trajectory of molecule B just before its atom B

_{1}hits atom A

_{2}of molecule A.

As just shown in Fig. 1, the attractive van der Waals forces that are explicitly accounted for by the LJ and FA models alter the straight molecular trajectories of classic kinetic theory revealing a longer and persistent interaction between colliding molecules.^{3} These interactions can lead to molecular trajectories that differ from those of classic kinetic theory as shown first by Hirschfelder *et al.*^{15} For example, molecules can split apart for a very short time (e.g., few picoseconds, ps) and after having traveled a very short distance (e.g., less than 0.2 Å), they can collide again and again before interacting with any other molecule. Such collisions are spurious that include “orbiting” ones.^{15} The origin of this motion is traced to the attractive component of the inter-atomic potential that is fully accounted for by the LJ and FA models.

Figure 2(a) (Multimedia view), for example, depicts snapshots from a pair of molecules experiencing seven spurious collisions within 17.1 ps according to the LJ model. The snapshots are in relative coordinates showing the trajectory of the moving molecule relative to the stationary. Thin red lines depict the trajectory of the actual spurious collisions while bold black lines the molecular trajectory before and after them. This is clearly a spurious collision event [Fig. 2(a)], that forms a so-called orbiting collision.

Similarly, Fig. 2(b) (Multimedia view) depicts the relative trajectory of a colliding molecule according to the FA model that experiences three successive collisions with the reference (stationary) molecule, moving on a highly curved trajectory before eventually splitting to continue its free motion along a new path. These events constitute a direct manifestation of the complex interaction between molecules due to their attraction, a feature that is absent in the HS model as in kinetic theory.

From far away, however, the mechanics of collisions seem similar between the three models: two molecules approach one another, interact strongly as the separation between them gets close to *r _{m}*, and finally they depart to different directions along straight paths before they are engaged in a new collision. Through hazard plot analysis,

^{61,62}spurious collisions are identified and filtered. All spurious collisions between two molecules are counted as one collision event in the calculations of collision density and

*λ*. For the HS model [Fig. S5(a) in the supplementary material], spurious collisions are totally absent as only regular collisions occur.

By tracking the motion of all molecules as closely as possible, almost every fs (i.e., every MD time integration step), multi-particle collisions can be identified also (see Sec. S7 of the supplementary material). By sorting out spurious collisions, the fraction of 2-, 3- [Fig. 3(a), Multimedia view], or even 4-body collisions [Fig. 3(b), Multimedia view] can be determined (Table S7 in the supplementary material). By far, the vast majority is 2-body or bimolecular collisions, as expected. For the HS model, this is always the case. However, for the LJ and FA ones, there also exists a small fraction of events involving more than two molecules. Three-body collisions makeup about 0.5% of all collisions according to the LJ model and about 0.8% according to the FA. Moreover, due to the attractive part of the potential, a tiny fraction of collisions (likelihood less than 5 × 10^{−3}% and 0.5 × 10^{−3}% for the FA and LJ models, respectively) involving four molecules was also identified. Such multi-body collisions should take place even at lower pressures than those employed here as there the mean free path becomes longer. Examples of how interatomic distances evolve in such 3- or 4-body collisions in the simulations by the FA model are shown in Fig. 3.

There, as the molecules move, they also rotate, resulting in wavy trajectories due to their non-spherical shape rather than the straight ones of kinetic theory. Note that the inter-atomic separation defines whether a collision takes place rather than the center-of-mass distance that can be longer than the collision threshold distance, *r*_{m}, between colliding molecules.

### C. Statistics of multi-particle collisions

Figure 3(a) shows that at *t* ≈ 2.5 ps, the inter-atomic separations *r*_{12}, *r*_{13}, and *r*_{23} of one O_{2} and two N_{2} molecules involved in a collision lie in pairs simultaneously below the threshold distance *r _{m}* for collision as indicated by the dashed orange and black lines (

*r*= 3.726 Å for two nitrogen atoms, and

_{m}*r*= 3.541 Å between nitrogen and oxygen atoms, see also Table S2 in the supplementary material), quantitatively identifying this as a 3-body collision. Similarly, in Fig. 3(b) at

_{m}*t*≈ 1.7 ps, three out of the six pairs of four molecules (all N

_{2}) that are involved in the event have inter-atomic separations below the shortest threshold distance for collision (

*r*= 3.726 Å for two nitrogen atoms), a direct manifestation of the 4-body character of this event.

_{m}### D. Mean free path

Table I shows the computed collision densities from the three models along with those from kinetic theory. The latter are in excellent agreement with those from the HS model, providing another strong validation of our MD simulations. Table I shows also that collision densities from the models accounting for attractive interactions (LJ and FA) are the highest as inter-molecular attraction enhances molecular proximity and, thus, collision frequency. For comparison, the corresponding collision densities accounting for spurious collisions by the LJ and FA models (Table S8 of the supplementary material) were more than an order of magnitude larger than those in their absence (Table I).

. | Collision densities (10^{33} m^{−3}s^{−1})
. | |||
---|---|---|---|---|

Oxygen–oxygen . | Nitrogen–nitrogen . | Nitrogen–oxygen . | Total . | |

Kinetic theory | 3.14 | 57.0 | 26.8 | 86.94 |

HS model | 3.2 ± 0.1 | 56.9 ± 1.6 | 26.7 ± 0.6 | 86.8 ± 1.7 |

LJ model | 5.3 ± 0.1 | 88.3 ± 2.4 | 43.3 ± 1.0 | 136.8 ± 2.6 |

FA model | 5.8 ± 0.2 | 96.1 ± 2.7 | 47.4 ± 1.2 | 149.3 ± 3.0 |

. | Collision densities (10^{33} m^{−3}s^{−1})
. | |||
---|---|---|---|---|

Oxygen–oxygen . | Nitrogen–nitrogen . | Nitrogen–oxygen . | Total . | |

Kinetic theory | 3.14 | 57.0 | 26.8 | 86.94 |

HS model | 3.2 ± 0.1 | 56.9 ± 1.6 | 26.7 ± 0.6 | 86.8 ± 1.7 |

LJ model | 5.3 ± 0.1 | 88.3 ± 2.4 | 43.3 ± 1.0 | 136.8 ± 2.6 |

FA model | 5.8 ± 0.2 | 96.1 ± 2.7 | 47.4 ± 1.2 | 149.3 ± 3.0 |

Established theoretical expressions^{1,3,29,68} for the mean free path [Eqs. (S7)–(S10) in Sec. S3 of the supplementary material] are based on the average time between collisions resulting in *λ* between 66.3 and 67.3 nm, for air at 300 K and 1 atm. Here, the *λ* is computed by recording and analyzing the actual distances traveled by molecules between collisions. A difficulty encountered with such a detailed approach is that, theoretically, these distances extend up to infinity rendering their quantification from the MD trajectories quite challenging. For example, molecules moving with exceptionally high speeds would result in free paths of exceptional length.^{29}

Consequently, the analysis had to be extended to sufficiently long times to capture rare occurrences of very long free paths and, of course, check if any collisions are spurious ones. To this, the observation time, *t*_{obs}, over which molecular paths are counted was progressively increased until convergence is attained (see Fig. S7 of the supplementary material). Then, at each *t*_{obs} the *λ* was calculated from the probability distribution of free paths accumulated up to that time.

Figure 4 shows how the three models result in slightly different distributions of *λ*′ for the same observation time, here *t*_{obs} = 6 ns. For the HS model, the corresponding regression with an exponential function is also shown. By far, the longer free paths are obtained by the HS model and the shortest ones by the FA model. This is not unexpected as attractive forces between molecules shorten such paths [Figs. 1(b) and 1(c)]. The free paths from the HS model follow largely an exponential distribution (in accord with kinetic theory) except for the very long ones which are characterized by the smallest probability of appearance.

Figure 5 shows that the asymptotic convergence of *λ* at very long observation times is faster for the FA and LJ models than for the HS. This is due to the fast convergence of shorter paths sampled with increasing observation time that outscores the appearance of longer ones. The latter are characterized, however, by an exponentially decreasing occurrence probability (Figs. S7 and S8 in the supplementary material). The so-derived *λ* is 42.3 ± 1.1 nm for the LJ model (circles) and 38.5 ± 1 nm for the FA model (triangles). For the HS model, *λ* is 70.6 ± 2.3 (solid line, squares) or 65.8 ± 1.9 nm (broken line) through direct averaging or averaging after performing the exponential regression (as suggested by kinetic theory), respectively. The latter (*λ* = 65.8 ± 1.9 nm) agrees very well with that obtained from Eqs. (S7)–(S9) of the supplementary material, *λ* = 66.4 nm.

The corresponding *λ* from direct averaging (70.6 ± 2.3 nm) is higher than this by about 6%. This level of agreement reveals that the distribution of free paths is predominantly (albeit not strictly) exponential, as has also been reported recently.^{69} Note that the “refined” value of 67.3 nm derived from Eq. (S10)^{70} of the supplementary material lies between these two values.

By far, however, the most important result is the significantly smaller *λ* obtained here than that in the literature^{1,68} when accounting for the force field and shape of air molecules by the FA model, *λ* = 38.5 ± 1 nm. This is almost 43% smaller than the widely used *λ* from kinetic theory (67.3 nm) but rather close to that from the LJ model (42.3 ± 1.1 nm).

That the true value of *λ* is smaller than that derived by the simple kinetic theory is consistent with the predictions of the kinetic theory variants. The calculation of the *λ* from these models is based on experimental measurements for the transport properties of air and all are summarized in Table S5 of the supplementary material. For example, the VHS model^{12} [Eq. (S12) in the supplementary material] gives *λ*_{VHS} = 53.6 nm. This is 20% smaller than that of the widely accepted value of *λ* = 67.3 nm, at T = 300 K and P = 1 atm. By making use of the VHS model expressions of Koura and Matsumoto,^{8} the mean free path can also be calculated from the air viscosity [Eq. (S15) of the supplementary material] or diffusivity [Eq. (S16) of the supplementary material] resulting in *λ _{η}* = 54.2 nm or

*λ*= 63.7 nm, respectively. Similarly, the VSS model

_{D}^{12,13}[Eqs. (S14) and (S24) of the supplementary material], a generalization of the VHS model, yields

*λ*

_{VSS}= 55 nm. All these

*λ*values are summarized in Table S5 of the supplementary material, including some additional ones referring solely to pure nitrogen gas. Although all these models consider air molecules as perfect elastic spheres, they suggest

*λ*that range from about 53.6 nm to about 67.3 nm. However, when accounting for the actual molecular nature of air (i.e., in terms of shape and interactions), the

*λ*is considerably smaller.

The *λ* from the FA model is fully consistent with the corresponding hazard plot analysis [Fig. S5(b) of the supplementary material], which shows that free paths between normal (non-spurious) collisions occur with a rate of 0.002 59 Å^{−1} suggesting a *λ* of 1/(0.002 59 Å^{−1}), i.e., of 38.6 nm. Similarly, for the LJ model [Fig. S5(a) of the supplementary material], the hazard rate for free paths between normal collisions is 0.002 39 Å^{−1}, suggesting a *λ* of 41.8 nm, again consistent with the asymptotic value *λ* = 42.3 ± 1.1 nm from the LJ model. For the simpler HS [Fig. S5(a) of the supplementary material], the calculated hazard rate is 0.001 45 Å^{−1}, corresponding to a *λ* of 69.0 nm, which also lies between the *λ* from the two ways of averaging the HS model results.

Very similar *λ* for the LJ and FA models are also obtained by just substituting into Eq. (S9) of the supplementary material the corresponding MD-predicted collision densities (see Table I): *λ* = 42.1 ± 0.9 nm for the LJ model, and *λ* = 38.6 ± 0.8 nm for the FA model. Overall, the agreement of the three different estimates (Table II) of the *λ* [(1) direct averaging and asymptotic analysis, (2) analysis of hazard rates, (3) kinetic theory based on collision densities—Eq. (S7) of the supplementary material] reflects the consistency of the above calculations.

. | Mean free path, λ (nm)
. | |||
---|---|---|---|---|

Direct averaging . | Hazard rate analysis . | Collision density . | ||

HS model | 70.6 ± 2.2 | 69.0 | 66.4 ± 1.8 | |

LJ model | 42.3 ± 1.1 | 41.8 | 42.1 ± 0.9 | |

FA model | 38.5 ± 1.0 | 38.6 | 38.6 ± 0.8 |

. | Mean free path, λ (nm)
. | |||
---|---|---|---|---|

Direct averaging . | Hazard rate analysis . | Collision density . | ||

HS model | 70.6 ± 2.2 | 69.0 | 66.4 ± 1.8 | |

LJ model | 42.3 ± 1.1 | 41.8 | 42.1 ± 0.9 | |

FA model | 38.5 ± 1.0 | 38.6 | 38.6 ± 0.8 |

The *λ* from the LJ and FA models in Table II have been computed for a collision distance equal to the characteristic LJ distance 2^{1/6}*σ _{ij}* at which the force between a pair of atoms (

*i*,

*j*) changes from attractive to repulsive as they approach each other. In Sec. S10 on the supplementary material, the

*λ*of neon gas is calculated using different collision criteria and compared to the corresponding results from Dongari

*et al.*

^{71}showing the significance of selecting consistently the collision distance with the employed force field.

### E. The impact of the selected force field on *λ*

To assess the sensitivity of the MD-derived *λ* on the parameterization of the employed FA model, the air density, viscosity, diffusivity, and *λ* were calculated with almost all two-body FA potentials available in the literature^{35–39} for air. The LJ parameters of these potentials are reported in Table S9 of the supplementary material. Then, in Table S10 of the supplementary material, the predicted density, diffusivity, and viscosity are reported from these force fields based on NVE simulations at *T* = 300 K and *P* = 1 atm. For three of these potentials (i.e., Bouanich,^{36} Vrabec *et al.*,^{37} and Zhang and Siepmann^{39}), these estimates agree with literature values and are fully consistent with those from the Zambrano *et al.* model^{32} employed here. On the other hand, the Wang *et al.* potential^{35} underestimates the density by ∼16% and overestimates the diffusivity and viscosity by ∼19% and ∼8%, respectively. Finally, the Vácha *et al.* potential^{38} yields the correct density but underestimates the diffusivity and viscosity by ∼50%.

The corresponding predictions of the FA models using above potentials for the *λ* of air are presented in Fig. S11 of the supplementary material. For all practical purposes, they coincide with the *λ* obtained from the Zambrano *et al.* model,^{32} thus confirming that the *λ* of 38.5 nm reported here is force field-independent. The only exception is the *λ* from the Vácha *et al.* potential^{38} which comes out to be considerably smaller, equal to *λ* = 20.5 nm; however, the large underprediction (up to 50%) of air transport properties by this model indicates that such a *λ* cannot be accepted.

As mentioned in Sec. II A, there are two “classic” potentials for nitrogen (one from Tokumatsu and Matsumoto^{41} and another from Kosyanchuk and Yakunchikov^{43}) that have been parametrized on the basis of *ab initio* simulations. Thus, to further check the sensitivity of the MD-derived *λ* on the atomistic force field parameterization, the *λ* for pure nitrogen was determined also with these two models and compared to that using the Zambrano *et al.* model^{32} employed here. The comparison is shown in Table S5 of the supplementary material. The corresponding predictions of the two models for the density and transport properties of nitrogen and their comparison against experimental data for density (obtained from the Engineering ToolBox^{63}), diffusivity^{72} and viscosity^{65,73} are shown in Table S11 of the supplementary material. The comparison is excellent for the density and diffusivity, as the largest discrepancy between simulation predictions and measured data are only 2%. For the viscosity the comparison is slightly noisier, with the largest deviation between the experimental data and the corresponding MD extracted values being equal to 9% for the Kosyanchuk and Yakunchikof potential.^{43}

Overall, all three force fields reproduce with remarkable accuracy both the transport properties and the density of nitrogen. Then, from Fig. S12 of the supplementary material, which shows the corresponding *λ* as a function of observation time, the following *λ* are obtained for pure nitrogen at 300 K and 1 atm: *λ* = 38.3 nm for the Tokumatsu and Matsumoto potential,^{41} and *λ* = 40.5 nm for the Kosyanchuk and Yakunchikov potential.^{43} The corresponding result from the Zambrano *et al.* model^{32} is 37.6 nm. The internal consistency of the three values is noteworthy given that intermolecular parameters differ markedly; for example, their *σ* differ by 7.3% and the *ε* by as much as 50%. Remarkably, the common value of *λ* derived from the MD simulations from the three models is about 40% smaller than what is suggested by Jennings' equation^{70} for pure nitrogen (*λ* = 66 nm). It appears that, not only for air but also for pure nitrogen, the *λ* predicted by atomistic MD simulations is force field-independent and, not surprisingly by now, much smaller than what it is accepted today from well-adopted theoretical correlations based on classic kinetic theory.

## IV. CONCLUSIONS

A detailed analysis of molecular collisions in dry air at ambient conditions was carried out by fully atomistic (FA) molecular dynamics models accounting for the force field and the diatomic state of O_{2} and N_{2}, for the first time, to the best of our knowledge. Closely following molecular trajectories, the occurrence of collisions was determined first. Then, by spotting and filtering out spurious collisions (i.e., counting all spurious collisions between two specific molecules as one collision event), the frequency and duration of genuine ones were computed. Two limiting cases assuming spherical molecules were examined also: (a) with (LJ model) and (b) without a force field resulting in purely elastic collisions (HS model) to directly compare the latter with classic kinetic theory.

Strong interactions between colliding molecules were revealed that resulted in various molecular trajectory types (including grazing, orbiting, etc. as first envisioned by Hirschfelder *et al.*^{15}) due to the rotating motion and non-spherical shape of air molecules. These interactions forced the molecules to spend considerable time together upon collision (“dragging”) before separating and moving to new directions. They also led to spurious collisions (multiple successive ones between the same molecules before separation) and a few 3- and even 4-body ones. All these are a direct result of the attractive component of the force field of gas molecules that is neglected by classic kinetic theory. A similar behavior was observed in the collisions with the LJ model. In contrast, no such collision dragging was observed with the simpler HS model, consistent with kinetic theory.

The collision densities between molecules as determined by the HS model were in agreement with kinetic theory, further validating the present analysis. However, the calculated collision densities by the LJ and, more importantly, by several FA models were much higher. This is attributed again to the attractive part of the non-bonded component of the inter-molecular potential in the LJ and FA models. Treating O_{2} and N_{2} as true diatomic molecules enhanced the relative population of multi-body collisions also. In essence, the analysis of intermolecular collisions presented here and the richness of collision events observed through FA simulations extends the theoretical study of Hirschfelder *et al.*^{15} for spherical LJ molecules to account for the exact diatomic nature of air molecules and their detailed interactions.

As a direct result of the enhanced interactions between oxygen and nitrogen when treated as diatomic molecules experiencing both repulsive and attractive interactions, their mean free path is considerably smaller than that given from the classic kinetic theory of gases and its variants.^{7,8,74} Although the latter still consider air molecules as elastic spheres, their predicted *λ* can be smaller than that of the classic kinetic theory (e.g., *λ*_{VHS} = 53.6 nm). However, when taking into account the actual molecular nature of air, its mean free path is 38.5 nm, 43% smaller than the widely used value of 67.3 nm at 300 K and 1 atm. Spurious collisions had been eliminated for the determination of the mean free path of air. The *λ* calculations were repeated using almost all classical force fields available for air today that can reliably reproduce its density, diffusivity, and viscosity. Practically identical *λ* were obtained using all these force fields in the FA model, proving that the estimated *λ* does not depend on the employed force field as long as the latter accurately estimates the experimentally determined viscosity, diffusivity, and density of air.

The above analysis points out that particle transport in the continuum regime (Navier–Stokes equations) can be safely employed down to smaller particles than used today. In contrast, much more detail is needed for process design in synthesis of much smaller particles (less than 5 nm). The interactions of such nanoparticles with surrounding gas are far more complex than those of larger particles that readily attain asymptotic size distributions (i.e., self-preserving) and structures (i.e., fractal-like) greatly facilitating process design and operation by computational fluid dynamics in manufacturing of carbon blacks, titania, and fumed silica.^{75} The present analysis has an immediate impact on the diffusivity of nanoparticles, *D*, as there is an ongoing debate on it for particles smaller than 5 nm. Worth noting here is already it is recognized that the force field of such small nanoparticles has to be accounted for to describe their interaction with gases even though the force field of gases seems to have been ignored until now. This could be useful also in description of compressible flows, where knowing the local mean free path can help determine whether one breaks the continuum/local equilibrium assumption within a grid cell in the vicinity of shocks and slip lines. Also, the lower mean free path determined here means that the drag force currently calculated for nanoparticle transport in neutral, magnetic, or electric fields is an underestimate. This can be quite important also in the development of instruments for characterization of nanoparticles from measurements of their transport properties (i.e., mobility).

## SUPPLEMENTARY MATERIAL

See the supplementary material for information regarding the impact of MD on aerosol nanoparticle formation, simulation details on all employed force fields and validation against experimental measurements, theoretical models for the calculation of the mean free path of air and nitrogen, Maxwell–Boltzmann velocity distribution, methodology followed in selecting the appropriate MD cell size, description of hazard plot analysis to track spurious collisions, statistics of multi-body and spurious collisions, methodology to calculate the mean free path, and the impact of collision distance and force field on the mean free path.

## ACKNOWLEDGMENTS

We cheerfully recognize the stimulating suggestions for future research and applications by one of the journal referees. This research was funded in part by Swiss National Science Foundation (Grant Nos. 175754 and 182668) and ETH Zurich.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Dimitrios Georgios Tsalikis:** Formal analysis (equal); Methodology (equal); Software (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). **Vlasis G. Mavrantzas:** Conceptualization (equal); Supervision (supporting); Writing – original draft (equal); Writing – review & editing (equal). **Sotiris E. Pratsinis:** Conceptualization (equal); Funding acquisition (lead); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this paper (if not present herein) are available on reasonable request from the authors.

## REFERENCES

*Smoke, Dust and Haze: Fundamentals of Aerosol Dynamics*

*Lectures on Gas Theory*

*Principles of the Kinetic Theory of Gases*

*The Mathematical Theory of Non Uniform Gases*

*Molecular Gas Dynamics and the Direct Simulation of Gas Flows*

*Nonequilibrium Gas Dynamics and Molecular Simulation*

*Molecular Theory of Gases and Liquids*

_{2}nanoparticles by molecular dynamics

*An Introduction to the Kinetic Theory of Gases*

*Kinetic Theory of Gases with an Introduction to Statistical Mechanics*

*The Mechanics of Aerosols*

_{2}, O

_{2}, H

_{2}, CO and CO

_{2}

_{2}, O

_{2}, O

_{3}, OH, H

_{2}O, HO

_{2}, and H

_{2}O

_{2}

*ab initio*ground-state potential energy surface of N

_{4}

_{2}O

_{2}

_{2}collisions

*ab initio*potential energy surface

_{2}–N

_{2}collisions at hyperthermal temperatures

_{2}collisions

*Computer Simulation of Liquids*

*Large-Scale Atomic/Molecular Massively Parallel Simulator*

*Methods for Statistical Analysis of Reliability and Life Data*

*Diffusion: Mass Transfer in Fluid Systems*

*Aerosol Technology: Properties, Behavior, and Measurement of Airborne Particles*