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

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 rm, and the impact parameter b of the collision (or rm in the absence of molecular interactions).3 

At a higher level of analysis,3 one solves the non-linear integrodifferential Boltzmann equation4 with the Chapman–Enskog method5,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 extensions7–11 of the hard-sphere model, such as the variable hard-sphere (VHS) model7 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 algorithms17 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 σ i j = ( σ i + σ 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.

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 Leroi33 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 fields44–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 nitrogen41,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 field32 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 thermostat56,57 and the Parrinello–Rahman barostat58 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 method59 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 × 104 air molecules (consisting of 79 mol. % N2 and 21 mol. % O2 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 O2 and N2 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.

The HS model includes only repulsive interactions between oxygen and nitrogen rigid spheres, and their collisions are purely elastic, consistent with classic kinetic theory,3,
u HS ( r ) = { for r σ i j 0 for r > σ i j ,
(1)
where r denotes the separation distance. To avoid the discontinuity at r = σ i j, the HS potential is approximated with a “softer” version for r σ i j,60,
u soft ( r ) = A [ 1 + cos ( π r σ i j ) ] for r σ i j ,
(2)
where A = 5000 kcal/mol, an energy prefactor.60 The van der Waals diameters are15  σ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 rules54 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.

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 rm is equal to the van der Waals sum of radii of the two molecules, σij. For the LJ and FA models, it is equal54 to 21/6σij, where i and j refer to entire nitrogen or oxygen molecules (for the LJ model) or their atoms (for the FA model). The rm 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 collisions15) 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.

Let h ( λ ) d λ be the probability that a molecule which has traveled length λ from its last collision experiences a new collision after having traveled length λ + dλ. The cumulative probability (cumulative hazard) will be H ( λ ) = 0 λ h ( λ ) d λ . Then, the expectation of the cumulative hazard at the length of the mth ordered collision (out of n collisions) is62,
H ¯ m = l = 0 l = m 1 n l .
(3)
For a set of uncorrelated collisions described by an exponential distribution (such as suggested by kinetic theory for molecular collisions) with a certain rate (say equal to 1/ξ), the cumulative hazard is a straight line62 (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.

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, tobs), 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 λ ( t obs ) = ( b 1 + b 2 t obs ) / ( 1 + a t obs ) with a, b1, and b2 numerical constants, and then taking the limit tobs → ∞ to get λ = lim t obs λ ( 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-sphere7 (VHS) and variable soft-sphere8,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

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 expression29 as shown in Fig. S2. The employed atomistic force field was validated by comparing the MD-predicted air density (1.180 ± 0.01 kg/m3) to the experimentally measured one for dry air (1.177 kg/m3).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 O2 and N2 molecules in the Fickian (linear) regime.64 For the viscosity, η, the corresponding Green–Kubo equation54 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 cm2/s and 18.1 ± 0.2 μPas, respectively, both in excellent agreement with experimental65,66 or theoretical studies64,67 giving D = 0.203 cm2/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/m3, 0.185 ± 0.005 cm2/s, and 18.0 ± 0.5 μPas from the HS model, and 1.18 ± 0.01 kg/m3, 0.2 ± 0.005 cm2/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.

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 rm (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 θm 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 χ + 2 θ m = π; 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.3 

FIG. 1.

Collision trajectories by the three models in relative coordinates: (a) HS model, (b) LJ model, and (c) FA model. In each case, the reference (stationary) molecule is white and the colliding one blue. The bold line shows the trajectory followed by the moving molecule before and after collision. The thin lines parallel to the trajectory of the moving molecule help define the impact parameter b. Multimedia available online .

FIG. 1.

Collision trajectories by the three models in relative coordinates: (a) HS model, (b) LJ model, and (c) FA model. In each case, the reference (stationary) molecule is white and the colliding one blue. The bold line shows the trajectory followed by the moving molecule before and after collision. The thin lines parallel to the trajectory of the moving molecule help define the impact parameter b. Multimedia available online .

Close modal

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 B1 hits atom A2 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.

FIG. 2.

Spurious collisions (successive ones between the same molecules) in relative coordinates during a single otherwise collision by (a) the LJ model with the pair of colliding molecules (depicted with white and blue beads) experiencing seven spurious collisions within 17.1 ps during a collision event that constitutes an “orbiting” collision.15 (b) The FA model with the pair of colliding molecules experiencing three successive (spurious) collisions within 4 ps. For clarity, the rotational motion of the reference molecule (white) has been suppressed here. In both models, the line (bold black before spurious collisions begin and after they end, and red in-between) shows the trajectory of the moving molecule. The spurious collisions are not counted as independent events but they are part of a single collision. Multimedia available online .

FIG. 2.

Spurious collisions (successive ones between the same molecules) in relative coordinates during a single otherwise collision by (a) the LJ model with the pair of colliding molecules (depicted with white and blue beads) experiencing seven spurious collisions within 17.1 ps during a collision event that constitutes an “orbiting” collision.15 (b) The FA model with the pair of colliding molecules experiencing three successive (spurious) collisions within 4 ps. For clarity, the rotational motion of the reference molecule (white) has been suppressed here. In both models, the line (bold black before spurious collisions begin and after they end, and red in-between) shows the trajectory of the moving molecule. The spurious collisions are not counted as independent events but they are part of a single collision. Multimedia available online .

Close modal

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 rm, 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.

FIG. 3.

(a) Evolution of separation distances between pairs of N2 and O2 molecules involved in a three-body collision: Distances r12, r13, and r23 indicate inter-atomic separations between the 12, 13, and 23 pairs of molecules, respectively, for the constituent atoms in the colliding molecules. Indexes 1 and 2 refer to the red and green N2 molecules, while index 3 refers to the blue O2 molecule. (b) The evolution of separation distances between colliding pairs of atoms in molecules 1, 2, 3, and 4 in a four-body collision: Distances r12, r13, r14, r23, r24, and r34 indicate inter-atomic separations between the 12, 13, 14, 23, 24, and 34 pairs of molecules, respectively, for the constituent atoms in the colliding molecules. Indexes 1, 2, 3, and 4 refer to the red, green, blue, and yellow N2 molecules. Black and orange horizontal dashed lines mark the threshold distances for collision between nitrogen–nitrogen and nitrogen–oxygen pairs of atoms (equal to 3.726 and 3.541 Å, respectively, in Table S2 of the supplementary material). Multimedia available online .

FIG. 3.

(a) Evolution of separation distances between pairs of N2 and O2 molecules involved in a three-body collision: Distances r12, r13, and r23 indicate inter-atomic separations between the 12, 13, and 23 pairs of molecules, respectively, for the constituent atoms in the colliding molecules. Indexes 1 and 2 refer to the red and green N2 molecules, while index 3 refers to the blue O2 molecule. (b) The evolution of separation distances between colliding pairs of atoms in molecules 1, 2, 3, and 4 in a four-body collision: Distances r12, r13, r14, r23, r24, and r34 indicate inter-atomic separations between the 12, 13, 14, 23, 24, and 34 pairs of molecules, respectively, for the constituent atoms in the colliding molecules. Indexes 1, 2, 3, and 4 refer to the red, green, blue, and yellow N2 molecules. Black and orange horizontal dashed lines mark the threshold distances for collision between nitrogen–nitrogen and nitrogen–oxygen pairs of atoms (equal to 3.726 and 3.541 Å, respectively, in Table S2 of the supplementary material). Multimedia available online .

Close modal

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, rm, between colliding molecules.

Figure 3(a) shows that at t ≈ 2.5 ps, the inter-atomic separations r12, r13, and r23 of one O2 and two N2 molecules involved in a collision lie in pairs simultaneously below the threshold distance rm for collision as indicated by the dashed orange and black lines (rm = 3.726 Å for two nitrogen atoms, and rm = 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 t ≈ 1.7 ps, three out of the six pairs of four molecules (all N2) that are involved in the event have inter-atomic separations below the shortest threshold distance for collision (rm = 3.726 Å for two nitrogen atoms), a direct manifestation of the 4-body character of this event.

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

TABLE I.

Collision densities by the hard sphere (HS), Lennard–Jones (LJ), and fully atomistic (FA) models in comparison with the kinetic theory of gases that assumes purely elastic spheres.

Collision densities (1033 m−3s−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 (1033 m−3s−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 expressions1,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, tobs, over which molecular paths are counted was progressively increased until convergence is attained (see Fig. S7 of the supplementary material). Then, at each tobs 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 tobs = 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.

FIG. 4.

Probability distributions p(λ′) of free paths λ′ from the HS (squares), LJ (circles), and FA (triangles) models at an observation time of tobs = 6 ns. For the HS model, also the best fit to the data is shown with an exponential function (solid line) following kinetic theory.

FIG. 4.

Probability distributions p(λ′) of free paths λ′ from the HS (squares), LJ (circles), and FA (triangles) models at an observation time of tobs = 6 ns. For the HS model, also the best fit to the data is shown with an exponential function (solid line) following kinetic theory.

Close modal

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.

FIG. 5.

Mean free path as a function of observation time (symbols) from the three models obtained by directly averaging over the corresponding density distributions of free paths, λ′. The errors in the calculations in the three different models are smaller than the symbol size. The solid lines are the best fits to the MD data with a hyperbolic function. The broken line for the HS model corresponds to extracting the q by first regressing the distribution with an exponential function (as suggested by kinetic theory) and then computing the average of this function.

FIG. 5.

Mean free path as a function of observation time (symbols) from the three models obtained by directly averaging over the corresponding density distributions of free paths, λ′. The errors in the calculations in the three different models are smaller than the symbol size. The solid lines are the best fits to the MD data with a hyperbolic function. The broken line for the HS model corresponds to extracting the q by first regressing the distribution with an exponential function (as suggested by kinetic theory) and then computing the average of this function.

Close modal

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 literature1,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 model12 [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 λD = 63.7 nm, respectively. Similarly, the VSS model12,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.

TABLE II.

Comparison of the λ of air at 1 atm and 300 K from the HS, LJ, and FA models by (a) Direct averaging over the total population of MD-obtained free paths followed by asymptotic analysis; (b) analysis of MD hazard rates; and (c) the MD-predicted collision densities [in Table I, using Eq. (S9) of the supplementary material] of kinetic theory.

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 21/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.

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 literature35–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 Siepmann39), these estimates agree with literature values and are fully consistent with those from the Zambrano et al. model32 employed here. On the other hand, the Wang et al. potential35 underestimates the density by ∼16% and overestimates the diffusivity and viscosity by ∼19% and ∼8%, respectively. Finally, the Vácha et al. potential38 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. potential38 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 Matsumoto41 and another from Kosyanchuk and Yakunchikov43) 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. model32 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 ToolBox63), diffusivity72 and viscosity65,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. model32 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' equation70 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.

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 O2 and N2, 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 O2 and N2 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).

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.

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.

The authors have no conflicts to disclose.

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

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

1.
S. K.
Friendlander
,
Smoke, Dust and Haze: Fundamentals of Aerosol Dynamics
(
John Wiley & Sons
,
New York
,
1977
).
2.
J. C.
Maxwell
, “
V. Illustrations of the dynamical theory of gases.—Part I. On the motions and collisions of perfectly elastic spheres
,”
London, Edinburgh Dublin Philos. Mag. J. Sci.
19
,
19
(
1860
).
3.
D. A.
McQuarrie
,
Statistical Mechanics
(
Harper and Row Publishers
,
New York
,
1976
).
4.
L.
Boltzmann
,
Lectures on Gas Theory
(
Dover Publications
,
New York
,
1964
).
5.
H.
Grad
,
Principles of the Kinetic Theory of Gases
(
Springer
,
Berlin
,
1958
).
6.
S.
Chapman
and
T.
Cowling
,
The Mathematical Theory of Non Uniform Gases
(
Cambridge University Press
,
Cambridge, UK
,
1970
).
7.
G.
Bird
, “
Definition of mean free path for real gases
,”
Phys. Fluids
26
,
3222
(
1983
).
8.
K.
Koura
and
H.
Matsumoto
, “
Variable soft sphere molecular model for inverse‐power‐law or Lennard‐Jones potential
,”
Phys. Fluids A
3
,
2459
(
1991
).
9.
K.
Koura
and
H.
Matsumoto
, “
Variable soft sphere molecular model for air species
,”
Phys. Fluids A
4
,
1083
(
1992
).
10.
H.
Hassan
and
D. B.
Hash
, “
A generalized hard‐sphere model for Monte Carlo simulation
,”
Phys. Fluids A
5
,
738
(
1993
).
11.
B. L.
Haas
,
D. B.
Hash
,
G. A.
Bird
,
F. E.
Lumpkin
III
, and
H. A.
Hassan
, “
Rates of thermal relaxation in direct simulation Monte Carlo methods
,”
Phys. Fluids
6
,
2191
(
1994
).
12.
G. A.
Bird
,
Molecular Gas Dynamics and the Direct Simulation of Gas Flows
(
Clarendon Press
,
Oxford UK
,
1994
).
13.
I. D.
Boyd
and
T. E.
Schwartzentruber
,
Nonequilibrium Gas Dynamics and Molecular Simulation
(
Cambridge University Press
,
2017
).
14.
E. A.
Mason
and
W. E.
Rice
, “
The intermolecular potentials for some simple nonpolar molecules
,”
J. Chem. Phys.
22
,
843
(
1954
).
15.
J. O.
Hirschfelder
,
C. F.
Curtiss
, and
R. B.
Bird
,
Molecular Theory of Gases and Liquids
(
Wiley
,
New York
,
1964
).
16.
J.
Huang
,
S.
Rauscher
,
G.
Nawrocki
,
T.
Ran
,
M.
Feig
,
B. L.
De Groot
,
H.
Grubmüller
, and
A. D.
MacKerell
, Jr.
, “
CHARMM36m: An improved force field for folded and intrinsically disordered proteins
,”
Nat. Methods
14
,
71
(
2017
).
17.
P.
Eastman
,
J.
Swails
,
J. D.
Chodera
,
R. T.
McGibbon
,
Y.
Zhao
,
K. A.
Beauchamp
,
L.-P.
Wang
,
A. C.
Simmonett
,
M. P.
Harrigan
, and
C. D.
Stern
, “
OpenMM 7: Rapid development of high performance algorithms for molecular dynamics
,”
PLoS Comput. Biol.
13
,
e1005659
(
2017
).
18.
S.
Chmiela
,
H. E.
Sauceda
,
K.-R.
Müller
, and
A.
Tkatchenko
, “
Towards exact molecular dynamics simulations with machine-learned force fields
,”
Nat. Commun.
9
,
3887
(
2018
).
19.
J. C.
Phillips
,
D. J.
Hardy
,
J. D.
Maia
,
J. E.
Stone
,
J. V.
Ribeiro
,
R. C.
Bernardi
,
R.
Buch
,
G.
Fiorin
,
J.
Hénin
, and
W.
Jiang
, “
Scalable molecular dynamics on CPU and GPU architectures with NAMD
,”
J. Chem. Phys.
153
,
044130
(
2020
).
20.
E.
Blaisten-Barojas
and
M.
Zachariah
, “
Molecular-dynamics study of cluster growth by cluster-cluster collisions
,”
Phys. Rev. B
45
,
4403
(
1992
).
21.
A.
Violi
, “
Modeling of soot particle inception in aromatic and aliphatic premixed flames
,”
Combust. Flame
139
,
279
(
2004
).
22.
B.
Buesser
,
A.
Grohn
, and
S. E.
Pratsinis
, “
Sintering rate and mechanism of TiO2 nanoparticles by molecular dynamics
,”
J. Phys. Chem. C
115
,
11030
(
2011
).
23.
B.
Buesser
and
S. E.
Pratsinis
, “
Morphology and crystallinity of coalescing nanosilver by molecular dynamics
,”
J. Phys. Chem. C
119
,
10116
(
2015
).
24.
E.
Goudeli
and
S. E.
Pratsinis
, “
Surface composition and crystallinity of coalescing silver-gold nanoparticles
,”
ACS Nano
11
,
11653
(
2017
).
25.
G. A.
Sotiriou
,
G. D.
Etterlin
,
A.
Spyrogianni
,
F.
Krumeich
,
J.-C.
Leroux
, and
S. E.
Pratsinis
, “
Plasmonic biocompatible silver-gold alloyed nanoparticles
,”
Chem. Commun.
50
,
13559
(
2014
).
26.
N.
Eom
,
M. E.
Messing
,
J.
Johansson
, and
K.
Deppert
, “
General trends in core–shell preferences for bimetallic nanoparticles
,”
ACS Nano
15
,
8883
(
2021
).
27.
A.
Sharma
,
K. M.
Mukut
,
S. P.
Roy
, and
E.
Goudeli
, “
The coalescence of incipient soot cluster
,”
Carbon
180
,
215
(
2021
).
28.
V. G.
Mavrantzas
and
S. E.
Pratsinis
, “
The impact of molecular simulations in gas-phase manufacture of nanomaterials
,”
Curr. Opin. Chem. Eng.
23
,
174
(
2019
).
29.
J.
Jeans
,
An Introduction to the Kinetic Theory of Gases
(
Cambridge University Press
,
Cambridge
1982
).
30.
E. H.
Kennard
,
Kinetic Theory of Gases with an Introduction to Statistical Mechanics
(
McGraw-Hill
,
New York
,
1947
).
31.
N. A.
Fuchs
,
The Mechanics of Aerosols
(
Dover Publications Inc
.,
Mineola, NY
,
1964
).
32.
H.
Zambrano
,
J. H.
Walther
, and
R.
Jaffe
, “
Molecular dynamics simulations of water on a hydrophilic silica surface at high air pressures
,”
J. Mol. Liq.
198
,
107
(
2014
).
33.
J.
Laufer
and
G.
Leroi
, “
Calculation of the zero wave vector lattice frequencies of alpha‐and beta‐oxygen
,”
J. Chem. Phys.
55
,
993
(
1971
).
34.
C.
Murthy
,
K.
Singer
,
M.
Klein
, and
I.
McDonald
, “
Pairwise additive effective potentials for nitrogen
,”
Mol. Phys.
41
,
1387
(
1980
).
35.
S.
Wang
,
K.
Hou
, and
H.
Heinz
, “
Accurate and compatible force fields for molecular oxygen, nitrogen, and hydrogen to simulate gases, electrolytes, and heterogeneous interfaces
,”
J. Chem. Theory Comput.
17
,
5198
(
2021
).
36.
J.-P.
Bouanich
, “
Site-site Lennard-Jones potential parameters for N2, O2, H2, CO and CO2
,”
J. Quant. Spectrosc. Radiat. Transfer
47
,
243
(
1992
).
37.
J.
Vrabec
,
J.
Stoll
, and
H.
Hasse
, “
A set of molecular models for symmetric quadrupolar fluids
,”
J. Phys. Chem. B
105
,
12126
(
2001
).
38.
R.
Vácha
,
P.
Slavíček
,
M.
Mucha
,
B. J.
Finlayson-Pitts
, and
P.
Jungwirth
, “
Adsorption of atmospherically relevant gases at the air/water interface: Free energy profiles of aqueous solvation of N2, O2, O3, OH, H2O, HO2, and H2O2
,”
J. Phys. Chem. A
108
,
11573
(
2004
).
39.
L.
Zhang
and
J. I.
Siepmann
, “
Direct calculation of Henry's law constants from Gibbs ensemble Monte Carlo simulations: Nitrogen, oxygen, carbon dioxide and methane in ethanol
,”
Theor. Chem. Acc.
115
,
391
(
2006
).
40.
K.
Koura
, “
Monte Carlo direct simulation of rotational relaxation of diatomic molecules using classical trajectory calculations: Nitrogen shock wave
,”
Phys. Fluids
9
,
3543
(
1997
).
41.
T.
Tokumasu
and
Y.
Matsumoto
, “
Dynamic molecular collision (DMC) model for rarefied gas flow simulations by the DSMC method
,”
Phys. Fluids
11
,
1907
(
1999
).
42.
P.
Valentini
,
C.
Zhang
, and
T. E.
Schwartzentruber
, “
Molecular dynamics simulation of rotational relaxation in nitrogen: Implications for rotational collision number models
,”
Phys. Fluids
24
,
106101
(
2012
).
43.
V.
Kosyanchuk
and
A.
Yakunchikov
, “
A detailed multiscale study of rotational–translational relaxation process of diatomic molecules
,”
Phys. Fluids
33
,
022003
(
2021
).
44.
Y.
Paukku
,
K. R.
Yang
,
Z.
Varga
, and
D. G.
Truhlar
, “
Global ab initio ground-state potential energy surface of N4
,”
J. Chem. Phys.
139
,
044309
(
2013
).
45.
Z.
Varga
,
R.
Meana-Pañeda
,
G.
Song
,
Y.
Paukku
, and
D. G.
Truhlar
, “
Potential energy surface of triplet N2O2
,”
J. Chem. Phys.
144
,
024310
(
2016
).
46.
Z.
Varga
,
Y.
Paukku
, and
D. G.
Truhlar
, “
Potential energy surfaces for O+ O2 collisions
,”
J. Chem. Phys.
147
,
154312
(
2017
).
47.
A.
Yakunchikov
,
V.
Kosyanchuk
,
A.
Kroupnov
,
M.
Pogosbekian
,
I.
Bryukhanov
, and
A.
Iuldasheva
, “
Potential energy surface of interaction of two diatomic molecules for air flows simulation at intermediate temperatures
,”
Chem. Phys.
536
,
110850
(
2020
).
48.
P.
Valentini
,
T. E.
Schwartzentruber
,
J. D.
Bender
,
I.
Nompelis
, and
G. V.
Candler
, “
Direct molecular simulation of nitrogen dissociation based on an ab initio potential energy surface
,”
Phys. Fluids
27
,
086102
(
2015
).
49.
D. A.
Andrienko
and
I. D.
Boyd
, “
Vibrational energy transfer and dissociation in O2–N2 collisions at hyperthermal temperatures
,”
J. Chem. Phys.
148
,
084309
(
2018
).
50.
M. S.
Grover
,
T. E.
Schwartzentruber
,
Z.
Varga
, and
D. G.
Truhlar
, “
Vibrational energy transfer and collision-induced dissociation in O + O2 collisions
,”
J. Thermophys. Heat Transfer
33
,
797
(
2019
).
51.
P.
Valentini
,
A. M.
Verhoff
,
M. S.
Grover
, and
N. J.
Bisek
, “
First-principles predictions for shear viscosity of air components at high temperature
,”
Phys. Chem. Chem. Phys.
25
,
9131
(
2023
).
52.
A.
Van der Avoird
,
P.
Wormer
, and
A.
Jansen
, “
An improved intermolecular potential for nitrogen
,”
J. Chem. Phys.
84
,
1629
(
1986
).
53.
J.-P.
Ryckaert
,
G.
Ciccotti
, and
H. J.
Berendsen
, “
Numerical integration of the cartesian equations of motion of a system with constraints: Molecular dynamics of n-alkanes
,”
J. Comput. Phys.
23
,
327
(
1977
).
54.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
(
Oxford University Press
,
Oxford
,
2017
).
55.
A. P.
Thompson
,
H. M.
Aktulga
,
R.
Berger
,
D. S.
Bolintineanu
,
W. M.
Brown
,
P. S.
Crozier
,
P. J.
in't Veld
,
A.
Kohlmeyer
,
S. G.
Moore
, and
T. D.
Nguyen
, “
LAMMPS-A flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales
,”
Comput. Phys. Commun.
271
,
108171
(
2022
).
56.
S.
Nosé
, “
A molecular dynamics method for simulations in the canonical ensemble
,”
Mol. Phys.
52
,
255
(
1984
).
57.
W. G.
Hoover
, “
Canonical dynamics: Equilibrium phase-space distributions
,”
Phys. Rev. A
31
,
1695
(
1985
).
58.
M.
Parrinello
and
A.
Rahman
, “
Polymorphic transitions in single crystals: A new molecular dynamics method
,”
J. Appl. Phys.
52
,
7182
(
1981
).
59.
L.
Verlet
, “
Computer “experiments” on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules
,”
Phys. Rev.
159
,
98
103
(
1967
).
60.
LAMMPS Users Manual
.
Large-Scale Atomic/Molecular Massively Parallel Simulator
(
Sandia National Labs
,
2003
).
61.
N. R.
Mann
,
R. E.
Schafer
, and
N. D.
Singpurwalla
,
Methods for Statistical Analysis of Reliability and Life Data
(
John Wiley and Sons, Inc
.,
New York
,
1974
).
62.
E.
Helfand
, “
Brownian dynamics study of transitions in a polymer chain of bistable oscillators
,”
J. Chem. Phys.
69
,
1010
(
1978
).
63.
Engineering ToolBox
, “
Air–density, specific weight and thermal expansion coefficient vs. temperature and pressure
” (
2003
).
64.
E. L.
Cussler
,
Diffusion: Mass Transfer in Fluid Systems
(
Cambridge University Press
,
New York
,
2007
).
65.
J.
Kestin
and
W.
Leidenfrost
, “
An absolute determination of the viscosity of eleven gases over a range of pressures
,”
Physica
25
,
1033
(
1959
).
66.
J.
Kestin
and
J.
Whitelaw
, “
The viscosity of dry and humid air
,”
Int. J. Heat Mass Transfer
7
,
1245
(
1964
).
67.
E. W.
Lemmon
and
R. T.
Jacobsen
, “
Viscosity and thermal conductivity equations for nitrogen, oxygen, argon, and air
,”
Int. J. Thermophys.
25
,
21
(
2004
).
68.
W. C.
Hinds
,
Aerosol Technology: Properties, Behavior, and Measurement of Airborne Particles
(
John Wiley & Sons
,
New York
,
1982
).
69.
J.
Xie
,
M. K.
Borg
,
L.
Gibelli
,
O.
Henrich
,
D. A.
Lockerby
, and
J. M.
Reese
, “
Effective mean free path and viscosity of confined gases
,”
Phys. Fluids
31
,
072002
(
2019
).
70.
S. G.
Jennings
, “
The mean free path in air
,”
J. Aerosol. Sci.
19
,
159
(
1988
).
71.
N.
Dongari
,
Y.
Zhang
, and
J. M.
Reese
, “
Molecular free path distribution in rarefied gases
,”
J. Phys. D
44
,
125502
(
2011
).
72.
E. B.
Winn
, “
The temperature dependence of the self-diffusion coefficients of argon, neon, nitrogen, oxygen, carbon dioxide, and methane
,”
Phys. Rev.
80
,
1024
(
1950
).
73.
J.
Kestin
,
S.
Ro
, and
W.
Wakeham
, “
Viscosity of the binary gaseous mixture helium‐nitrogen
,”
J. Chem. Phys.
56
,
4036
(
1972
).
74.
G.
Bird
, “
Direct simulation of gas flows at the molecular level
,”
Commun. Appl. Numer. Methods
4
,
165
(
1988
).
75.
S. E.
Pratsinis
, “
Aerosol-based technologies in nanoscale manufacturing: From functional materials to devices through core chemical engineering
,”
AIChE J.
56
,
3028
(
2010
).

Supplementary Material