A new equation for the mean free path of air

Abstract A recent rigorous methodology for determining the mean free path, λ, of air at ambient conditions [Tsalikis, D. G., Mavrantzas, V. G., & Pratsinis, S. E. (2023). Physics of Fluids, 35, 097131] is extended to pressures, P, from 0.5 to 5 atm and temperatures, T, from 100 to 3000 K that are of environmental and industrial relevance, allowing to derive a new simple equation for λ as a function of T and P. This entails molecular dynamics (MD) simulations of air accounting for the actual shape and force field of nitrogen and oxygen and analysis of the computed microcanonical ensemble of free paths to derive λ at any P and T. Simulations are rigorously validated by comparing MD-predicted air densities, diffusivities and viscosities at various temperatures and pressures against experimental measurements, theoretical expressions and ab initio simulations. At all T and P, the MD-computed λ are systematically smaller (at least 40%) than those from classic kinetic theory and its variants. The new equation for λ is: λ(T,P)=0.033916×T1.23P−1 and can be also expressed in terms of the standard Jennings’ expression, λ=π/8μu1ρP [Jennings, S. G. (1998). Journal of Aerosol Science, 19(2), 159–166], relating λ to air viscosity, μ, and density, ρ, using just a new value for its numerical factor u, u = 0.81475 ± 0.00288, which is 63% larger than the current u = 0.4987445.


Introduction
The mean free path of air, k, the average length traveled by air molecules between their successive collisions, is a fundamental quantity in aerosol science.By definition, k is intimately connected with collisions between gas molecules (Maxwell 1860); however, through kinetic theory it also enters well-established mathematical expressions defining important macroscopic properties of air (Jeans 1982).The mean free path of air further dictates if particle transport takes place in the free molecular, transient, or continuum regimes (Fuchs 1964).Distinctly different theories of particle dynamics hold in the three regimes, therefore, an accurate k is important for correctly assessing the relevant transport regimes in aerosol processes.Even more important, aerosol dynamics in each regime are described by equations that make direct use of k.This is the case, e.g., with the drag force or the particle diffusivity in neutral, magnetic and/or electric fields (Larriba-Andaluz and Carbone 2021).Clearly, an accurate knowledge of k is a prerequisite for choosing the appropriate theoretical framework for studying aerosol dynamics in the atmosphere (Seinfeld and Pandis 2016), for process design and/or scale-up of aerosol synthesis of materials (Ulrich 1971) and aerosol instrumentation.
Unfortunately, measuring directly the mean free path of air is not feasible.On the other hand, k can be indirectly estimated or implicitly extracted using kinetic theory equations that relate it to transport properties that are rather straightforward to measure (e.g., viscosity) (Jeans 1982).However, kinetic theory (see Section S1 in the online supplementary information [SI]) expressions for the viscosity do not describe well its temperature dependence (Boyd and Schwartzentruber 2017), which necessitated improvements of the theory (Bird 1994;Boyd and Schwartzentruber 2017) known as the Variable Hard Sphere (VHS) model (Bird 1983), the Variable Soft Sphere (VSS) model (Koura andMatsumoto 1991, 1992), and the Generalized Soft Sphere (GSS) model (Hassan and Hash 1993).Overall, the indirect way of computing k entails significant uncertainties because of the approximations invoked in the corresponding theories, the most important being the assumption that gas molecules are spheres colliding elastically with a potential energy function that can be described by a hard sphere or a Lennard-Jones (LJ) model with fitted (optimized) parameters.It is thus not surprising that k values estimated from these theories (Bird 1983;Jennings 1988;Koura andMatsumoto 1991, 1992) are not consistent but vary from about 54 to about 67 nm at ambient conditions of temperature, T ¼ 300 K and pressure, P ¼ 1 atm (Table S5 in Tsalikis, Mavrantzas and Pratsinis 2023).
To overcome such uncertainties in the calculation of k from approximations in the corresponding theoretical treatments, one can resort to detailed molecular simulations at the atomistic level where the true atomic structure is maintained and accurate force fields are available (Bender et al. 2015;Bouanich 1992;Varga, Paukku and Truhlar 2017;Wang, Hou and Heinz 2021).In contrast to theoretical treatments that assume variable (i.e., weakly dependent) collision diameters with temperature to capture the correct temperature scaling of air viscosity (Boyd and Schwartzentruber 2017), atomistic simulations utilize universal parameters (e.g., atomic diameters) independent of temperature and pressure.
Recently, such an approach was introduced in order to compute k at room temperature and pressure by combining explicit-atom molecular dynamics (MD) simulations with a detailed collision algorithm capable of identifying molecular collisions (binary and multiple), spotting spurious events and filtering them out, and computing k by analyzing the ensemble of free paths accumulated in the course of such simulations (Tsalikis, Mavrantzas and Pratsinis 2023).It was found that air molecules experience a wide variety of collision patterns, including spurious (i.e., successive collisions for the same pair of molecules until they split apart moving in different directions) and multibody ones (i.e., collisions involving more than two molecules), both attributed mainly to the two-pronged nature (attractive and repulsive) of the LJ potential and, to a lesser extent, to the diatomic shape of the molecules.Spurious collisions between two molecules were counted as one collision event in the calculation of collision densities and k.The mean free path of air was determined as k ¼ 38.5 ± 1 nm at T ¼ 300 K and P ¼ 1 atm, almost 43% smaller than the 67.3 nm k value that is widely accepted today at these conditions (Jennings 1988).This result was independent of the employed force field, as long as this could reproduce the macroscopic values of air density, diffusivity and viscosity.Such a considerably smaller k implies that the Navier-Stokes equations (i.e., the continuum regime) can be safely employed for even smaller particles than had been thought so far.
Here, the above simulations are extended to a wide range of temperatures and pressures, namely from 0.5 to 5 atm and from 100 to 3000 K, to cover the entire regime of relevance to atmospheric and industrial aerosol processes.As in our previous study, the employed force field (Zambrano, Walther and Jaffe 2014) was extensively re-validated by comparing predictions for several macroscopic observables of air to experimental measurements, theoretical expressions and ab initio simulation results over the above range of T and P (see section S2 in the SI).The simulation results for k are further analyzed as a function of T and P, and the respective scalings are compared with those from kinetic theory (Jeans 1982) and its variants (Bird 1983;Koura andMatsumoto 1991, 1992), but above all against Jennings' widely-accepted analytical formula providing k in terms of air density and viscosity (Jennings 1988).A new equation is thus proposed for k as a function of T and P which is further mapped onto the Jennings formula.

Molecular dynamics (MD)
The simulations were conducted using LAMMPS (Thompson et al. 2022) with the same atomistic model (Zambrano, Walther and Jaffe 2014) as in our previous study (Tsalikis, Mavrantzas and Pratsinis 2023) that incorporates both intramolecular (i.e., bond stretching) and intramolecular (i.e., attractive and repulsive) LJ interactions: with r and e denoting the atomic diameter and depth of the energy well, respectively.The MD runs were performed with rather large cubic cells containing 5 � 10 4 air molecules (79 mol % N 2 and 21 mol % O 2 ) (Tsalikis, Mavrantzas and Pratsinis 2023).For example, at T ¼ 300 K and P ¼ 1 atm, the simulation cell had a length of 126.7 nm, corresponding to an air concentration of about 0.041 M. At least five different cell sizes (38 − 126.7 nm) were explored to ensure the absence of system-size effects.To ensure the reproducibility of the MD results, simulations were repeated three times, utilizing initial configurations with entirely different atomic positions and velocities.In all cases, fully consistent results were obtained.To avoid perturbations of the atomistic trajectories mainly by the barostat (Figure S1 in Tsalikis, Mavrantzas and Pratsinis 2023) and, to a lesser extent, by the thermostat, all simulations were carried out in the microcanonical (NVE) ensemble using a timestep of 1 fs for the integration of Newton's equations of motion.Initial configurations for the NVE simulations were generated after long NPT runs at the T and P conditions of interest.The NVE MD runs were conducted at density and energy conditions corresponding to temperatures from 100 to 3000 K and pressures from 0.5 to 5 atm.The pressure and temperature in each simulation were controlled extremely well as the corresponding statistical fluctuations were always within less than 1% of the mean.

Collisions
A collision is detected based on a threshold distance r m , utilizing a distance criterion (Tsalikis, Mavrantzas and Pratsinis 2023).This is the interatomic separation below which the corresponding interatomic force switches from negative to positive, i.e., from attractive to repulsive.Since a 12-6 LJ potential is used to describe the interatomic force (Allen and Tildesley 2017), r m ¼ 2 1/6 r ij where r ij denotes the characteristic LJ diameter for the colliding pair of atoms i and j.
Table S1 shows the r m for collisions between all combinations of oxygen and nitrogen atoms.

Mean free path
The methodology employed to compute the mean free path entails the following steps (Tsalikis, Mavrantzas and Pratsinis 2023): first, from the accumulated MD trajectories, all (literally all) free paths between successive collisions are computed by tracking particle dynamics in every integration step.Subsequently, spurious collisions are detected and filtered using Hazard plot analysis (see Sections S3 and S4 in the SI).Then, for a given simulation or observation time, t obs , the corresponding probability density of free paths and their mean value are computed, to obtain the mean free path at infinite observation time.For example, the temporal evolution of k with t obs is regressed with a hyperbolic function of the form: with a, b 1 , and b 2 being numerical constants, and the limit t obs ! 1 is taken to get:

Validation: Density, diffusivity, and viscosity as a function of temperature and pressure
Before calculating the mean free path, a careful validation of the force field used was performed by comparing MD-predicted macroscopic thermodynamic and transport properties of air against experimental measurements and theoretical predictions.Figure 1 shows that the MD-extracted air density (symbols) is in excellent agreement with experimental measurements (continuous line (Engineering ToolBox 2003)) for dry air at T ¼ 100 -3000 K and P ¼ 1 atm.The MD-extracted density values in Figure 1, originate from simulations in the isothermal-isobaric ensemble (see Methods), conducted to equilibrate air energy and density at the T and P conditions of interest, prior to the subsequent production runs in the microcanonical ensemble (NVE).Also shown in Figure 1, are density predictions by the thermodynamic formulation of Lemmon et al. (2000) according to which, the excess Helmholtz free energy of air is given by:  Lemmon et al. (2000).The statistical uncertainty of the MD-predicted density is smaller than the size of the symbols.
where d is the reduced density d ¼ q=q j and s the reduced temperature s ¼ T j =T: Symbols T j and q j correspond to the maxcondentherm (maximum coexistence gas and liquid) temperature (i.e., T j ¼ 123:6312 K) and density (i.e., q j ¼ 10:4477 mol=dm 3 ), respectively.All parameters (N k , i k , j k and l k ) in Equation ( 3) are reported in Table S2 of the SI.Then the density can be estimated through derivation of Equation (3) as: where R is the ideal gas constant.For all temperatures examined, the MD-extracted densities are found in excellent agreement with the theoretical predictions by the Lemmon et al. model.Similarly, the diffusion coefficient of air, D, as was estimated from the MD simulations with the Einstein equation was compared to the widely-used Chapman-Enskog (Chapman and Cowling 1970) equation: and the well-trusted empirical correlation by Fuller, Schettler and Giddings (1966): where T is the temperature in Kelvin, P the pressure in atm, M 1 and M 2 the molecular weights of the gas species (N 2 and O 2, respectively), r 12 the LJ diameter computed as the arithmetic mean, r 12 ¼ 0.5(r 1 þ r 2 ), of the corresponding r 1 and r 2 diameters for oxygenoxygen and nitrogen-nitrogen species, and X the collision integral which depends on species composition and temperature (Cussler 2007).In Equation ( 6), v i,1 and v i,2 are the diffusion volumes of N 2 and O 2 (Table S3).The MD-computed D values agree well with those from Equation ( 6), as shown in Figure S1 in the SI, particularly between 200 and 1000 K, where the corresponding discrepancy is less than 8%.In this temperature region, the D from Equation (5) are higher than those from Equation ( 6) and the MD simulations by 17% and 26%, respectively.The temperature scaling of D from the MD simulations is described well by a power law of the form D � T 1.70 ± 0.02 , which is in reasonable agreement with the scalings suggested by Equations ( 5) and ( 6), i.e., D � T 1.67 and D � T 1.75 , respectively.Note that temperature scaling from Equation ( 5) is stronger than T 1.5 (i.e., D � T 1.67 ) due to the temperature dependence of X (Cussler 2007).
A third validation involved the comparison of the MD-predicted viscosity of air, l, at various temperatures (under constant pressure, P ¼ 1 atm), which was obtained with the help of the corresponding Green-Kubo expression (McQuarrie 1976) through the integral of the time autocorrelation function of the off-diagonal components of the atomic stress tensor, with data from several sets of experimental measurements (Hellemans, Kestin and Ro 1973;Matthews et al. 1976) and recent simulation results utilizing a complex ab initio force field (Valentini et al. 2023).For each temperature, averaging over at least forty independent trajectories was performed to calculate the viscosity.The MD results are also compared with the predictions from the theoretical models of Lemmon and Jacobsen (2004) and Kadoya, Matsunaga and Nagashima (1985).The air viscosity in the Lemmon and Jacobsen (2004) model is computed as the sum of two terms, l 0 and l r , corresponding, respectively, to the dilute gas viscosity and the residual fluid viscosity (i.e., l ¼ l 0 þ l r ).The dilute gas viscosity is (Lemmon and Jacobsen 2004): where M is the molar mass, r the diameter of the molecule and X the collision integral: where T � ¼ T=ðe=k B Þ, e the Lennard-Jones energy parameter and k B the Boltzmann constant.The residual fluid viscosity is given by: with s ¼ T c =T and d ¼ q=q c where q is the density.The model parameters r, e, M, T c , and q c are given in Table S4.The collision integral parameters (i.e., b i 's) in Equation ( 8) are given in Table S5.The coefficients (i.e., N i , c i ) and exponents (i.e., t i , l i ) utilized in Equation ( 9) are reported in Table S6.As far as the model of Kadoya, Matsunaga and Nagashima (1985) is concerned, similarly to the work of Lemmon and Jacobsen, the air viscosity is the sum of the dilute l 0 and excess viscosity Dl : and with q r ¼ q=q 0 : The parameters A 1 , A 0.5 and A i as well as B i and q 0 are reported in Tables S7 and S8 of the SI.Then, the total viscosity according to Kadoya, Matsunaga and Nagashima (1985) is: with T r ¼ T=T 0 : Parameters H and T 0 are reported in Table S9 of the SI.
The temperature dependence of the air viscosity from 200 to 4000 K at 1 atm is presented in Figure 2 together with the corresponding predictions from the aforementioned models (Kadoya, Matsunaga and Nagashima 1985;Lemmon and Jacobsen 2004), experimental measurements (Hellemans, Kestin and Ro 1973;Matthews et al. 1976), and simulations using ab-initio force fields (Valentini et al. 2023).
The agreement between MD-extracted results, experimental measurements and various theoretical expressions appears to be quite satisfactory.For temperatures above 4000 K, in particular, ab initio simulations have shown that the role of electronic excitation has to be taken into account.For example, Valentini et al. (2023) reported that at T ¼ 10,000 K, ab initio simulations of the ground state underestimate the viscosity by 15% compared to theoretical data that include electronically excited states.
Overall, the consistency between the MD-extracted values for the density, diffusion coefficient and dynamic viscosity of air and experimental measurements or theoretical expressions is very good and fully validates the employed atomistic force field in the simulations.

Multi-particle collisions: Impact of temperature and pressure
Table 1 reports the occurrence probabilities of 2-and multi-body (i.e., 3-, 4-, and 5-body) collisions for temperatures in the regime 100 -3000 K at 1 atm.Contributions from spurious collisions have been removed from the analysis.At ambient conditions most collisions are bimolecular (�99.16%),but there also exist a small (but non-negligible) fraction of 3-(�0.84%)and 4-body ones (0.005%) which arise mostly from the attractive part of the interatomic potential.With decreasing temperature, the fraction of the multi-body collisions increases.In particular, at the cryogenic temperature of T ¼ 100 K, the fraction of 3-and 4-body collisions increases to 10.48% and 0.45% respectively, whereas 5body collisions (0.02%) are also observed.On the other hand, all multi-body collisions decrease with increasing temperature due to the associated density reduction and the enhancement of the thermal motion of molecules (higher molecular velocities).Indeed, with decreasing density, the mean distance between air molecules increases, thus, the probability of observing triplets or quadruplets of air molecules at small separations (e.g., less than 12 Å, namely the cutoff length of the LJ potential (Tsalikis, Mavrantzas and Pratsinis 2023)) for possible collision decreases.On the other hand, as the temperature increases, so do the molecular velocities (according to the Maxwell-Boltzmann distribution), implying higher kinetic energies and thus a higher probability of escaping the attractive well of the LJ potential exerted on them by  neighboring molecules as they pass close to each other, thus helping them to avoid collisions.
Results for the pressure dependence of the occurrence probabilities for two-and many-body collisions at constant temperature T ¼ 300 K are reported in Table 2.
As already discussed, bimolecular collisions dominate collision statistics.However, with increasing pressure, the probability of observing multi-body collisions increases considerably.For example, at the highest pressure here (P ¼ 5 atm), the total fraction of multi-body collisions reaches almost 5.1%.The fraction of five-body collisions (see Figure 3 below and Figure S2 in the SI, and Supplementary Movie FA-5-body) also increases significantly, while many more six-body collisions (see Figure S3 in the SI, and Supplementary Movie FA-6-body) are observed.
A representative example of such a collision (at T ¼ 500 K and P ¼ 5 atm) between five N 2 molecules is shown in Figure 3 (see also Supplementary Movie 1).In panels a-c of Figure 3, the five N 2 molecules approach each other while in panel (d) they start colliding in pairs (the red collides with the green, and the orange with the yellow), in triplets (the yellow with the blue and the orange, and the green with the blue and the red), even in quadruples.After the five-body collision (panel d in Figure 3), the five N 2 molecules are splitting apart, thus moving away (Figures 3e-g) from each other.At even lower temperatures (T ¼ 200 K) multi-body collisions become evidently more pronounced with increasing pressure (see Table S10 in the SI).Following the above reasoning, the pressure dependence of the many-body collisions can be understood in terms of the density increase accompanying the increase in pressure which reduces intermolecular separation and, therefore, increases the likehood of finding more than two molecules simultaneously at close enough separations to experience collision.

Spurious collisions: Impact of temperature and pressure
A representative example of a spurious collision (five successive collisions between the same molecules) at T ¼ 500 K and P ¼ 0.5 atm is presented in Figure S4 (see also Supplementary Movie 3).The population of spurious collisions (n spurious ) is obtained and compared to the total number of collisions, namely the sum of spurious and regular collisions, n total ¼ n spurious þ n regular .The successive collisions involved in a spurious event are explicitly accounted for the calculation of the spurious collision population, n spurious .Figure 4 shows the fraction of spurious collisions, u spurious ¼ n spurious /n total , as a function of T and P. The u spurious is mainly temperature-dependent, in agreement with Hirschfelder et al. (Hirschfelder, Curtiss and Bird 1964) who showed that orbiting collisions are possible when the relative kinetic energy between two colliding molecules is below a limiting value.Indeed, at the level of interatomic interactions discussed here, for a two-body collision to occur, two atoms must be at a separation distance r smaller than 2 1=6 r to experience a net repulsive force.On the other hand, at higher separations in the range 2 1=6 r < r � r c (r c ¼ 12 Å is the cutoff of the LJ interatomic potential), the force exerted on them is attractive, which increases their relative velocity.The molecules grapple against the attractive force, and if their relative energy (that also depends on the impact parameter, b (Hirschfelder, Curtiss and Bird 1964)) is sufficiently large to overcome it, they can escape and move away.On the other hand, if their relative kinetic energy is not that large, attraction prevails and the two molecules will collide.In fact, two molecules will remain together until their relative energy increases to such a level that they can escape their attraction or until a third molecule intervenes.Overall, spurious collisions have little effect on the transport properties at high densities, while at high temperatures, their population is dramatically reduced, thus their impact on transport properties is also reduced (Hirschfelder, Curtiss and Bird 1964).

New equation for the mean free path as a function of temperature and pressure
The evolution of the probability density distribution of free paths for observation times t obs from 1 to 30 ns at T ¼ 3000 K and P ¼ 1 atm is shown in Figure S5.With increasing t obs , the distribution of sampled paths shifts clearly to larger values, but with a decreasing rate.With increasing observation time, longer and longer free paths are sampled but with a significantly lower probability of occurrence (practically the relative population of free paths decreases exponentially with increasing path length), thus eventually reaching convergence (Tsalikis, Mavrantzas and Pratsinis 2023).Indeed, in Figure S5, the distributions of free paths for t obs ¼ 20 ns and t obs ¼ 30 ns practically overlap.Figure S6, on the other hand, shows the k at P ¼ 1 atm, for T ¼ 200 -3000 K, together with the corresponding hyperbolic function lines (Equation ( 2)).
At lower temperatures and high pressures, convergence occurs within a few nanoseconds of t obs , but at the reverse conditions, longer paths must be sampled, thus, convergence occurs after several tens of nanoseconds of t obs (see Figure S7).Then, by regressing the MD extracted k for different t obs with a hyperbolic function (Equation ( 2)), one can obtain the mean free path for practically infinite t obs .Results from the simulations at the various T and P performed here can then be fitted to a simple expression of the form kðT, PÞ ¼ bT c P d with b, c, and d numerical constants.During the regression, the constants b and c were allowed to vary so as to minimize the statistical uncertainties between the MD data and the regression equation at 200-3000 K which is of environmental and industrial relevance, whereas d (the exponent for the power-law pressure scaling) was set to d ¼ −1.This led to: with k in nm, T in K, and P in atm.This equation describes well the pressure and temperature dependence of the MD-computed k as shown in Figure 5 below.A three-dimensional graphical representation of Equation ( 13) is shown in Figure 6, including the raw MD-extracted k values (symbols).The corresponding determination coefficient is R 2 ¼ 99.97%, suggesting a rather accurate description of the simulation data, with the mean and maximum deviations between MD-measured values and Equation (13) being equal to 1.8 and 4.3%, respectively, in the region 200 K < T � 3000 K for all pressures examined.It is also nice that the MDderived k shown in Figure 6 (symbols) lie almost everywhere on the (T, P) surface of Equation ( 13), which indicates that the new equation provides a reliable description of the raw data independent of the specific (T, P) pair in the range of phase state points covered.
Figure 7 shows the (a) temperature (at constant P ¼ 1 atm) and (b) pressure (at constant T ¼ 300 K) dependence of the MD-extracted k. Figure 7 shows also the corresponding k from Equation ( 13), the kinetic theory of gases (Jeans 1982) and its variants (Bird 1983;Koura andMatsumoto 1991, 1992), and from the well-established Jennings' equation (Jennings 1988).Jennings made use of a more straightforward expression for the mean free path of air by minimizing the dependence on physical quantities (Jennings 1988): where u a numerical factor equal to 0.4987445 derived assuming gas molecules as rigid elastic spheres (Pekeris and Alterman 1957).
The MD-extracted k values in Figure 7 are systematically lower than those predicted by all other methods (the kinetic theory with its VHS and VSS variants, and by Jennings).On the other hand, the temperature dependence of k from Equation (13) (i.e., k � T 1.23 ) agrees with the Jennings one (i.e., k � T 1.23 ± 0.01 ).In contrast, the kinetic theory scaling (k � T 1 ), which suggests a weaker temperature dependence, is not supported by the MD data and is clearly connected with the inability of the theory to follow the temperature dependence of the viscosity (Bird 1994).The variants of the kinetic theory (VHS and VSS) had resolved this issue; thus, they predict larger scalings, closer to the Jennings' and present MD ones.Concerning the pressure dependence of k, Figure 7b shows that all approaches point out to the same scaling, namely, k � 1/P.Also, it becomes clear that state points at the same density do not have the same k, in sharp contrast to the kinetic theory predictions.13) with the MD predictions at the simulated T and P conditions (symbols).The coefficient of determination of the fitting of the MD data with Equation ( 13) is R 2 ¼ 99.97%, implying excellent matching.

New equation for the mean free path as a function of viscosity and pressure
The remarkable agreement of the temperature and pressure scalings of the MD-predicted mean free paths with the trend of the Jennings expression (Jennings 1988) shown in Figure 7 implies that the latter expression may be used to describe the present results for k with just a different numerical factor.Figure S8 shows the ratio of all MD-extracted k to those obtained from Jennings' equation (Equation ( 14)) for all T and P conditions examined here.This ratio is rather a constant and equal to k MD /k Jennings ¼ 0.583 ± 0.021, suggesting that the MD-obtained k can be described with Equation ( 14) but with a different numerical factor u. To compute this u, the air densities obtained from the MD simulations were used, but for the corresponding viscosities, due to the remarkable agreement between MD and theoretical models (Kadoya, Matsunaga and Nagashima 1985;Lemmon and Jacobsen 2004) (Figure 2), the viscosities of Lemmon and Jacobsen (2004) were adopted.Subsequently, the following equation is proposed for k over P ¼ 0.5 -5 atm and T ¼ 100 -3000 K: A direct comparison of the MD-derived k (symbols) with the above updated Jennings expression for k is shown in Figure 8 where k is plotted as a function of the air viscosity l and of the combined variable q P to enable a graphical representation of the two sets of data as a 2-D surface in 3-D space, i.e., k ¼ k(l, q P).Clearly, the updated expression describes accurately the MD-extracted mean free paths all over the entire (l, q P) space.This is also reflected in the corresponding coefficient of determination, R 2 ¼ 99.92%, indicating perfect matching between MD and the updated Jennings expression, Equation ( 15).
The precise new value of the numerical constant is u ¼ 0.81475 ± 0.00288, almost 63% larger than that proposed by Jennings (u ¼ 0.4987445) based on the earlier calculations for the viscosity coefficient of air at different state points by Pekeris and Alterman (1957).(Bird 1983) and VSS (Koura andMatsumoto 1991, 1992), are shown by the yellow and green lines, respectively.The red line shows the prediction from the new equation for k (Equation ( 13)).In all cases, the corresponding scaling exponents for the dependence of k on T and P are also shown.

Conclusions
A systematic analysis of molecular collisions in dry air over a broad range of temperatures and pressures was performed using a recently proposed methodology that accounts for the force field and diatomic structure of oxygen and nitrogen molecules (Tsalikis, Mavrantzas and Pratsinis 2023).By systematically analyzing atomistic trajectories, several collision patterns were identified, including grazing, head-on, and orbiting ones as had been envisioned many years ago by Hirschfelder, Curtiss and Bird (1964).The analysis also revealed spurious events (namely, multiple successive collisions between the same molecules before they split apart) and multi-body collisions (i.e., involving 3, 4, 5 and even 6 molecules).The classic kinetic theory does not account for the latter since its basis, the Boltzmann equation (Boltzmann 1964;Grad 1958), considers only binary interactions even though there are such models for chemical reactions (Boyd and Schwartzentruber 2017) or granular particles (Duan and Feng 2019) that consider non-binary collisions.These complex collision patterns are attributed primarily to the attractive part of the LJ potential and, to a lesser extent, to the non-spherical shape of air molecules (i.e., their rotating motion), both neglected by the classic kinetic theory and its main variants.Multi-body collisions are strongly affected by P and T. For example, at T ¼ 300 K, with increasing P from 0.5 to 5 atm, the fraction of multi-body collisions increases from �0.5 to �10%; and the reverse is observed with increasing temperature.Spurious collisions are practically pressure-independent because pressure alters the average distance between molecules as it affects the system volume.On the other hand, such collisions strongly depend on temperature, that affects the relative kinetic energy between colliding molecules (as implied from the Maxwell-Boltzmann molecular velocities).As spurious collisions occur when that energy drops below a limiting value, temperature strongly affects their frequency.
By detecting and filtering spurious collisions, probability density free path distributions were also computed, which served to extract the mean free path as a function of temperature and pressure.Due to the presence of attractive interactions between molecules, the computed mean free path values were systematically smaller than those predicted by the classic kinetic theory and the well-accepted correlation by Jennings (1988) as well as from the variants (Bird 1983;Koura andMatsumoto 1991, 1992) of the kinetic theory.
The temperature scaling (at constant pressure) of the MD-computed k values (�T 1.23 ) is almost identical to that of Jennings and agrees reasonably well with those suggested by all kinetic theory variants.On the other hand, the temperature scaling from the classic kinetic theory is incorrect, reflecting the well-known inability of the latter to correctly describe the scaling of the viscosity with temperature.The kinetic theory variants had resolved this issue (Bird 1994).On the other hand, all models agree with the MD calculations concerning the pressure scaling (k � P −1 ).
A new equation for the mean free path of air is also proposed as a function of temperature and pressure, namely kðT, PÞ ¼ 0:033916 � T 1:23 P −1 , which describes the MD-obtained k remarkably well.Due to the similarity of the pressure and temperature scalings between MD and the expression by Jennings involving the density and viscosity of air, the MD-computed k were also regressed with the Jennings expression (Jennings 1988), allowing its numerical factor u to vary.The new value is u ¼ 0.81475 ± 0.00288 and, when inserted in the Jennings expression, describes the MD-derived k very well.The new value for u is about 63% higher than that reported by Jennings, u ¼ 0.4987445 (Pekeris and Alterman 1957).
The new equation allows to reliably compute the k of air across a range of conditions of industrial and atmospheric relevance, spanning the regime from cryogenic temperatures (e.g., 100-250 K of interest for ice nucleation (Haag et al. 2003;Mahrt et al. 2020)) all the way up to high subsonic air flow temperatures (e.g., 3000 K).For even higher temperatures, it is expected that electronic excitation becomes important, thus classical force fields such as the one employed here to compute k become inapplicable (Valentini et al. 2023).
Our work highlights the critical importance of the force field in describing interactions between gas molecules, thus also between gas molecules and nanoparticles as well as between nanoparticles (particularly less than 5 nm, in the free molecular regime).Although for nanoparticle interactions, accounting for the true interatomic potential has recently received recognition, the detailed potential for interactions involving gases had been neglected until now (i.e., Larriba-Andaluz and Carbone 2021).The interactions of such nanoparticles with the 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) (Pratsinis 2010).
Having an equation that provides the correct k as a function of temperature and pressure is not only critical to precisely distinguishing between the various dynamic regimes (free molecular, transition, and continuum) for particle transport but also to describing important phenomena in each regime such as particle diffusion, condensation, coagulation and agglomeration, to name a few.The present results are relevant also to gas diffusion in porous media (Freeman, Moridis and Blasingame 2011;Lu and Zhao 2004) and to drag force and mobility of aerosol nanoparticles (Larriba-Andaluz and Carbone 2021).

Figure 1 .
Figure 1.Temperature dependence of the air density (red circles) at constant pressure P ¼ 1 atm.The line is the best power-law fit to the data reported in Engineering Toolbox 2003 and the blue crosses are the predictions using the theoretical formulation byLemmon et al. (2000).The statistical uncertainty of the MD-predicted density is smaller than the size of the symbols.

Figure 2 .
Figure 2.Temperature dependence (at P ¼ 1 atm) of the viscosity of air as extracted from the MD simulations (blue squares) and comparison against the theoretical expressions ofLemmon and Jacobsen (2004) (red line) andKadoya, Matsunaga and Nagashima (1985) (green line), experimental data(Hellemans, Kestin and Ro 1973;Matthews et al. 1976) (black symbols), and simulations(Valentini et al. 2023) using ab-initio force fields (cyan circles).The highlighted region shows the standard deviation in the MD-computed l, calculated by averaging over at least 40 independent trajectories at each temperature.

Figure 3 .
Figure3.Representative snapshots from the MD simulation at T ¼ 500 K and P ¼ 5 atm showing five N 2 molecules colliding simultaneously (five-body collision).These molecules are depicted by red, blue, green, yellow, and orange dumbbells.The molecules approach each other (a-c) and start colliding at t 3 ¼ 1.3 ps (d).The green nitrogen collides simultaneously with the red and blue nitrogens.At the same time, the blue and orange nitrogens collide with the yellow nitrogen.After the collision, the molecules drift apart (e-g).The center-of-mass trajectories of the colliding molecules are also shown with lines colored as the corresponding molecules.

Figure 4 .
Figure 4. Temperature dependence of the fraction of spurious collisions (ratio of spurious collisions to the total number of collisions) for pressures equal to 0.5, 1, 3 and 5 atm.

Figure 5 .
Figure 5. (a) Temperature and (b) pressure dependence of the mean free path, k, from the present MD simulations (symbols) and the corresponding fittings with Equation (13) (lines).

Figure 6 .
Figure 6.3-D representation of the temperature and pressure dependence of mean free path k from Equation (13) with the MD predictions at the simulated T and P conditions (symbols).The coefficient of determination of the fitting of the MD data with Equation (13) is R 2 ¼ 99.97%, implying excellent matching.

Figure 7 .
Figure 7. (a) Temperature dependence of k values at P ¼ 1 atm.(b) Pressure dependence of k at T ¼ 300 K.The black and blue lines are the k from Jennings (1988) and kinetic theory (KT) (Jeans 1982), respectively.The predictions of the two kinetic theory variants, VHS(Bird 1983) and VSS(Koura and Matsumoto 1991, 1992), are shown by the yellow and green lines, respectively.The red line shows the prediction from the new equation for k (Equation (13)).In all cases, the corresponding scaling exponents for the dependence of k on T and P are also shown.