key: cord-0756420-07jfiz1k
authors: White, Nathan; Seelig, John-David; Loyalka, Sudarshan K.
title: Computation of drag and diffusion coefficient for coronavirus: I
date: 2021-05-07
journal: J Aerosol Sci
DOI: 10.1016/j.jaerosci.2021.105806
sha: 3bd70914745bb06e6ea392d7b70909ac568eb8c1
doc_id: 756420
cord_uid: 07jfiz1k

Monte Carlo simulations and integral equation techniques allow for the flexible and efficient computation of drag and diffusion coefficients for virus mimetic particles. We highlight a Monte Carlo method that is useful for computing the drag on biomimetic particles in the free-molecular regime and a numerical technique to compute a (boundary integral solution) of the Stokes equation in the hydrodynamic limit. The free-molecular and the continuum results allow the construction of an approximation for the drag applicable over the full range of Knudsen numbers. Finally, we outline how this work and will be useful in modeling viral transport in air and fluids and in viral morphology measurements and viral separations via electrospray-differential mobility analyzers (ES-DMA).

The coronavirus is believed to be mainly transported in air as embedded in drops and droplets (see an earlier study on pathogen transport in air by Vejerano and Marr (2018) and a recent study by Lieber, et al. (2021) . 1, 2 The evaporation and the drying of these drops and droplets can lead to viral (or clumps of them) transport in the later stages or through filters, masks, and human respiratory tract as partially wet or dry.

While convective processes are primarily responsible for transporting viruses and other bioaerosols in bulk, diffusion plays a crucial role in their deposition and capture on surfaces. The morphology of the coronavirus and other bioaerosols (such as pollen) are quite complex; thus, analytical solutions or approximations for the drag coefficient or diffusion coefficients of such particles are of interest. Additionally, greater understanding of the impact complex shape/morphology has on diffusion also applies to related areas of study involving general (bio) particle transport.

An area of considerable current interest, is the use of electrospray differential mobility analyzers, referred to as ES-DMA. [3] [4] [5] [6] [7] [8] [9] In this technique, a fine spray from a liquid-virus solution is introduced into an ionizing chamber and into a flowing gaseous system (air or other gasses at atmospheric pressure). As the droplets change size, the bioparticles accumulate charge and begin to approximate a Boltzmann distribution; and sensors can then measure the bioparticle mobility (DMA) and concentration (with a Condensation Nuclei Counter). The technique is useful for virus sizes from 20 to 800 nm, and the mobility measurements can be quite precise. Since the mobility depends on the viral structure, the technique allows for both structural determination and separation of the viruses.

The inverse problem of computing viral structure from mobility measurements requires knowledge of mobility dependence on the structure itself. The diffusion coefficient, D is related to the particle mobility B through the well-known Einstein relation. The mobility itself can be obtained from the calculation or measurements of drag force D F on a particle, and this has been an area of considerable interest in aerosol science and engineering. [10] [11] [12] [13] [14] [15] [16] [17] [18] While the early work was focused on spherical particles and simple nonspherical geometries, there has been work also on non-spherical particles and particle fractal aggregates through use of Monte-Carlo techniques. [19] [20] [21] Realistic viral morphology is non-spherical; even approximately spherical capsids have myriad protein protrusions. Indeed, in coronaviruses the surface spike glycoproteins give the viruses their distinctive, eponymous shape. To our knowledge there has not been a thorough investigation of the drag force these spike proteins induce on viral motion in fluids.

Since viral particle diameters span a range of sizes compared to the mean free path of air (~0.06µm), it is useful to have available results that can span the entire size range (or alternatively, Knudsen number). We have explored a range of computations in this regard in both in the free molecular and the continuum limits and explored approximations over the entire Knudsen number range. For computational efficiency, we have approximated spikes on a virus surface by small spheres, but we will be exploring spheroidal, conical, and other shapes in the future. 7 

Following earlier works, we consider a particle situated in an infinite expanse of a gas flowing with a macroscopic velocity u (we will choose its direction as the x-direction). Then the distribution far away from it is described by a Maxwellian: 10

With nondimensional velocity scaling factors: (2) In which 0 n is the number density, m is the mass of a gas molecule, B k is the Boltzmann constant, and T is the temperature. The vector u is the molecular-velocity while c refers to the particle-velocity, respectively, and the starred variants stand for nondimensionalizations of those quantities. Also,

is an absolute Maxwellian.

We will also assume that any molecules reflected from a surface follow Maxwell's diffuse-specular reflection model:

, 1 :

Where     0 , f rc is a Maxwellian (with a density normalized by conservation of mass), and R is the reflection operator. Heren is the normal unit vector to the surface pointing in the gas, and α is the coefficient of diffuse specular reflection (α is one for purely diffuse reflection and zero for purely specular reflection).

We assume that each gas molecule infinitesimally contributes to the net drag, via the relation:

In the free molecular regime intermolecular collisions are exceedingly rare, and collisions with convex surfaces are non-reentrant. In this case, one may compute the incipient distribution directly, through equation (1) . Particularly, for a spherical particle of radius R one gets:

Note that in equation (6) , the scalar pressure is:

The mobility then is,   (9) Note that these results depend on both the molecular density and the molecular mass, among other parameters. We will express the nondimensionalized drag as: 2 DD F F pR u  åå (10) Using equation (6) , the drag (for a sphere of radius R ) in equation (10) becomes:

And has the values 13.16 for purely diffuse scattering and 9.45 for scattering that is purely specular. We note that Williams gave an alternate derivation of this result for the case of diffuse reflection, via a calculation of the collisional cross-sections. Williams also considered simultaneous heat transfer and further gas-surface interaction models. 18 For particles smaller than the molecular mean free path, intermolecular collisions are effectively disregarded, but for non-reentrant surfaces (meaning that molecules reflected from surface point to surface point; either directly, through concave geometry, or indirectly via sky-shine physics), especially ones with non-convex geometries (as viral surfaces can be), it is not possible to obtain helpful analytical expressions.

The reentrant case involves difficult integral equations that account for reflections, and most authors have either resorted to numerical techniques in the face of intractability or used numerical methods such as Monte Carlo algorithms for the flexibility they offer in exploring a wide range of shapes and surface distributions such as found on viral surfaces. Geometric cells as generally used in DSMC, are not needed as intermolecular collisions are neglected. 19 We also note that both Chan and Dahneke and Mackowski have previously used the TPMC technique to calculate drag on straight chain of uniform particles and agglomerates, respectively. We have found the procedure of Mackowski particularly convenient.

We have constructed a Mathematica program using this procedure to explore dependence of drag on a spherically shaped virus, with varying distributions of smaller spheres (to approximate spikeglycoproteins) on its surface. In constructing this program, we were able to repurpose Wolfram Language software we constructed previously to model condensational growth on aerosol chains and agglomerates. 23, 24 Since Mackowski described the procedure in detail, we only note essential points. Following Mackowski, we enclose the subject particle(s) in an artificial rectangular (or square) box, and distinguish between surface areas that are parallel, A or perpendicular A  , to the direction of motion, and sample the starting molecules from these surfaces in proportion to their particle fluxes.

We sample the molecular speeds (note that Mackowski instead used a mean speed approximation) to support generality and flexibility in computing different gas/surface molecular scattering kernels. We sampled the speed, c, by simply using a root finding algorithm on an equation of the form: 

is a uniform random variate on the unit interval, and nZ . This is (theoretically) equivalent to sampling from a particular scaling of the Probability Distribution Function (PDF) of a gamma distribution, or more specifically, the PDF of Nakagami's one-sided distribution, which is itself a special case of the gamma distribution. For an arbitrary mean and standard deviation this distribution is:

We experimented with writing our own sampling routines, as well as using the provided Wolfram Language random variate routines. Our experience was that the provided routines offered superior speed compared to our own, and for our test cases the difference between results obtained from using equation (12) or (13) was not significant (see also the Appendix).

We verified the program against known results for a single sphere first and have obtained results in agreement with those from equation (10), for five different values of α. Next, we considered a number, N, of equal sized spheres, in a straight chain. We report these results in Table 1 .

Subtracting the single sphere values from the two sphere values we deduced results for the Basic Chain Unit (BCU). We report these results in Table 2 along with a benchmark comparison from Chan and Dahneke.

In terms of the BCU, the drag for N spheres in a straight chain is approximately:

And the results of Table 1 are in accord with this.

J o u r n a l P r e -p r o o f (10), for a the case of a linear chain of spheres (i.e., a group of 2spheres embedded in a 3 dimensional Euclidean space), similar to a straight line of croquet or billiard balls touching one another, but not resting on any surface. The results are as expected; note that the drag increases slightly with the accommodation coefficient for a given geometry. In this work we approximated coronavirus bioparticles by modeling spike glycoproteins as spherical protrusions from a primary body meant to approximate the viral capsid. We first took each of these protrusions to be 1/16 th the radius of the primary sphere and considered sixteen different distributions of protrusions as shown in Figure 1 .

The ratio between the radius of the main body and the radius of the protrusions was selected based on some preliminary findings of a coronavirus structure that relate the length of peplomers to capsid body. 25 The choice of these configurations was incremental and meant to extend complexity from a bare central single sphere (in geometry one) to the fully covered, "bumpy" sphere (in geometry sixteen).

The geometries that explored here are the result of trial and error in creating an artistic rendition, of the novel coronavirus. The result of that process was a stochastic geometry that randomizes itself during each rendering; we show one rendition in Figure 3 . However, even if the randomized geometry resembles illustrations of coronaviruses that have appeared in online and in popular media, this induces too much complexity for a first study.

Rather than performing hundreds or thousands of simulations and reporting statistical distributions that are themselves statistical results on Monte Carlo simulations, we opted to simulate a bare sphere, then iteratively add protrusions to incrementally examine the effect of adding these protrusions on the drag experienced by the particle system.

We report our corresponding TPMC results for drag in Table 3 . We additionally computed the nondimensionalized drag obtained by the Volume Equivalent Sphere (VES) approximation. The radius of the volume equivalent sphere is and the corresponding drag, (16) Where , DR F å is the nondimensionalized drag as given by equation (5), with a radius R. Note that the volume equivalent drag is: (17) and these results are straightforward to compute from the data in Table 3 . We see that the VES approximation is quite good (we note however that for the straight chain of uniform spheres, the VES approximation is better when the number of spheres is small (say, four or less) and worsens when n is larger).

To explore the usefulness of the volume equivalent sphere approximation further, we have obtained added results by scaling the radii of the geometry by a factor of ¾ and performing the same computations for the geometric layouts defined by cases two through sixteen in Figure 1 . We report these results in Table 4 .

We present an image of the geometries used for the Monte Carlo simulations and various calculations throughout the paper. Note that case number one is unshown as it consisted of a single, isolated sphere. Geometries two and three are identical save for the orientation of their single protrusions. In each of the geometries pictured, under the nondimensionalizations used in this paper, the larger, central body-particle has a radius of one. The exception to this would be the geometries used in the construction of Table 3 . The Table 3 Results of the drag Monte Carlo simulation across five different coefficients of diffuse/specular reflection, for sixteen different geometries meant to mimic the novel coronavirus SARS-CoV-2. Please note that the simulation index, N, is arbitrary, and the sphere number, S, is the number of protrusions from the main body. Within each geometry set, the top lines are the results of our Monte Carlo Simulations, the second lines are the drag experienced by a volume equivalent sphere, and the third lines are the ratios of the first two, respectively. The drag for the volume equivalent sphere cases was computed with the aid of equations Drag computed via equations (5) - (11) . The results show that for each simulation geometry, an increase in the quantity of the capsid surface corresponds to an increase in the drag experienced by the particle. The uncertainty in the drag, which as the result of the Monte Carlo simulation, is the less precise than values obtained via evaluation of analytic formulae, drives the uncertainty of the calculations in general. Table 3 , the primary difference being that we scaled down the spheres used to model the particles by a factor of ¾, to ease some comparison on the relationship between corona size and drag. Please note that we did not duplicate experiment one replicated, as it consisted of a single bare sphere. Otherwise, the descriptions from Table 3 are still applicable, as well as the pictures in panel graphic Figure 1 . We discuss here first the results relating to Figure 1 and Table 3 . We reference each computation by our so-called simulation number, S, which is a completely arbitrary index, and the number of small spheres on the larger sphere as the index N. However, not all geometries are alike and for visual purposes readers may reference Figure 1 . Exact coordinate chart data showing the specifications for each geometry are available online via the publisher.

We note from Table 3 that the addition of spheres to the surface increases the drag. The drag also increases as diffuse reflection begins to dominate over specular reflection. Comparing simulation five to simulation twelve (both have eight protrusions, with a different spatial distribution), we see trivial difference between the two drag results.

However, extending the number of spheres and the complexity we see differences increase such that disorder and spread in sphere location increases drag, as seen specifically in comparing eleven with sixtytwo small spheres and fourteen with forty-four small spheres. Despite having less spheres, simulation fourteen had higher drag and deviates further from the VES approximation.

Finally, simulation sixteen with 144 spheres has significant deviation from VES approximation greater than 10%. This simulation acts as the extreme example that a spread-out layout of smaller spheres and highlights the advantage of performing the computations.

Comparing the geometries with the larger geometries (small-sphere radius of 1/16) and the smaller geometries (small-sphere radius of 3/64) is insightful. We see that in the smaller cases there is much more uniformity with respect to the number and distribution of the small spheres, compared to the larger case where there is a much more noticeable increase in drag force as the number of protrusions increases. We noticed the same increase in drag with the accommodation coefficient, α, for the scaled-down simulation.

We have computed drag on a range of virus configurations in the free-molecular limit by constructing and using a versatile Test Particle Monte Carlo (TPMC) program. Our results show that increasing protrusions result in higher drag. Our results also show that diffuse molecular-reflection increases drag, but also has a higher sensitivity to the number (and placement) of protrusions. Finally, we find that the volume equivalent sphere approximation supplies a good estimation of the drag for the types of geometries we studied, under the assumptions we made. The volume equivalent sphere approximation however gives an increasingly worse estimate of the drag as the number of protrusions increases.

We should note that in our program, we used about four million particles and parallel processing (4 cores). The program can be slow for complicated geometry computations (partly due to incorporation of visualization and use of Mathematica). We have not explored its optimization as it is in the developments stage. Its use on a hundred or thousand core processor (for example) does not require any change.

Further, we have reported above drag results only along the x-axis (in line with our assumption of gas flow on this axis only), and the geometries we considered are spherically symmetric. However, we also computed the full mobility tensor and computed the averaged drag force via the harmonic mean:

Where ii b are the eigenvalues of the mobility tensor; we have reported typical values in Table 5 .

J o u r n a l P r e -p r o o f 16 

We further note that the ion mobility, K, is

where ze is the charge on the ion, and B is the mobility defined earlier,

During the initial review of the paper, we were informed Larriba-Andaluz and Hogan Jr (2013, 2013, 2014) had computed K for ions composed of aggregates of spheres, in an executable form of their program IMoS (Ion Mobility Software), Coots, et al. (2020) , which is available for download and public use. [25] [26] [27] [28] [29] They had considered two approaches for calculating K (essentially the mobility, B, using TPMC (following Mackowski, with provisions for both cuboid and spherical enclosures), and the Mason-Schamp first approximation of the two-temperature theory: 26 

And hence,

Where, n is the gas (molecular) density,  is the reduced mass ( the two results agree only, when α is zero, i.e., the case of specular reflection. The differences are likely due to different assumptions regarding molecule-molecule (ion) interactions (diffuse scattering vs elastic collisions), and have led to several options for elastic, inelastic scattering, etc. in IMoS). 26 The issue is not the first order Chapman-Enskog approximation in this case (it does not differ much from the exact result for the case of rigid spheres. 32 (A Mathematica program for evaluation of collision integrals for elastic scattering for Lennard Jones potentials is given in the book by Ivchenko, 10 and can be changed to suit other potentials and/or inelastic scattering).

Note that, in terms of the collisional integral, one can express the nondimensionalized drag force as:

We report results from the IMoS program for the non-dimensional drag in Table 6 ; the results are comparable with the results we obtained from our program. The table shows results from the IMoS program for the nondimensionalized drag for the biomimetic geometries presented in Figure 1 . The results show that for each simulation geometry an increase in the total surface area of the particle corresponds to an increase in drag experienced by the particle. 

Applications of both TPMC and the DSMC to computations for arbitrary Knudsen number would be most useful, but the Monte Carlo methods have limitations for infinite domain problems because of 1/ r decay (of moments such as velocity) at large distances and the convergence for Knudsen numbers in the continuum can be extremely slow or deceptive. If one can truncate the problem geometry without affecting the accuracy of the simulation, and if particle scattering events are isotropic, then one can use the so-called Walk on Spheres (WoS) method(s) to perform a particularly efficient (grid free) form of Monte Carlo by recasting integral equations as stochastic processes. They are particularly suitable for graphics problems; see the work of Sawhney and Crane (2020) for a recent review. 33 Due to the restrictions of WoS methods, stemming from the subtlety of the connection between the Itô calculus and PDEs, especially their inability to handle infinite domains, and Neumann-or Robin-type boundary conditions, an alternate method was desired.

For more background information about the techniques discussed here, we recommend (in addition to other works cited in this section) the reference texts by Wyld and Powell (2020), Leal (2007), Pozrikidis (1992) , respectively. [34] [35] [36] In the continuum regime, the drag and torque can be computed through solutions of the Stokes Equation, and for a rigid body(s) of arbitrary shape by solving the boundary integral equation corresponding to it: 36

Where the Green's function (tensor) is: (27) Here u is the fluid translational velocity, and Ω is the angular velocity far from the fixed rigid body. Details of the derivation of equation (26) from the Stokes equation are well explained by Leal (we refer the reader to chapter eight of that text). 35 For the computations of interest here we will assume, 0  Ω , i.e. that the rotational component vanishes. Note that in the literature, authors may label the equation (27) as the Oseen tensor (especially when it is derived from Oseen's equations; fundamental solutions may be further called Oseenlets) or the Stokes tensor (with fundamental solutions sometimes called Stokeslets). To be specific, the difference between the two is that technically within an Oseen analysis there is consideration of convective (pure fluid) or advective (fluid as a carrier) acceleration, while in a Stokes analysis these effects are unconsidered. Here,   s T r is the stress tensor, and       n   s s s φT r r e r is the vector of stress components we need to compute the drag. Note that in the equation above s r is a point on the body, and ˆn e is a unit vector normal to the surface at s r directed from the fluid (gas) to the body. As usual, μ is the viscosity. We obtain the drag then by integrating   s r φ over the surface(s) of the body. We use the subscript 'cont.' to show that this is the continuum result, and to distinguish it from the free molecular result, which we will show with a 'fm.' subscript.

Equation (26) is a system of three singular integral equations; Youngren and Acrivos (1975) solved it numerically, by using a boundary element method for certain spheroids. 36 We have found it far simpler to solve it by using a singularity subtraction technique (for details of the technique we refer the reader to previous work on condensation rate on aggregates of spheres, spheroids and ellipsoids in the continuum regime, and on the spinning rotor gauge in the slip regime). [38] [39] [40] We write equation (26) as:

Where we have defined the integrated Green's tensor:

We evaluate the tensor   s r g numerically, and use Gaussian quadrature for the second integral, noting that integrand is zero for ss   rr , and is unevaluated at that point. For Gaussian quadrature, of order Q, we use QQ  points corresponding to the polar and azimuthal integrations over (0, )  and (0,2 )  , respectively. We use a collocation method (over all surfaces in the geometry) to solve the resulting set of rank-3 equations (there are three equations corresponding to each surface point, one for each spatial direction, and the total number of equations is thus

where N is the number of bodies. We have not assumed any spatial symmetry, and there is considerable flexibility in choosing the integration and quadrature schemes; for example, one could use separate quadrature meshes in the polar and azimuthal integrations).

In the cases (spheres, spheroids including needles, discs, linear chains of spheres), that we used for verification, we found very satisfactory results even with low order quadrature. Additionally, we have successfully used the program for ellipsoids also, but are still testing and improving it. We will be reporting separately the full details of the technique with analytical benchmarking for both translation and rotation of arbitrary bodies (see the monograph by Williams and Loyalka (1991) for a detailed description of cases and results). 12 Note that for a single sphere of radius R, one has the well-known Stokes result: 

We report numerical results that we obtained for the case of a linear chain of equal sized spheres (assuming fluid flow 1 u  ue , i.e. flow solely along the x-axis) in Table 7 . Our numerical results agree with equations (31) and (32), almost precisely, even with a low order quadrature. We have shown added results also in Table 8 .

We have reported results for drag, for selected configurations of the virus mimetic particle of Table 3 and corresponding to , 1 u  ue in Table 8 . These results are helpful in tandem with those for the respective free molecular cases, in the approximations discussed in the next section. We have chosen here configurations with fewer spheres as the computations are time consuming, despite parallelization. We plan to report more extensive results with spikes shaped as thin needles, in future papers.

We have also calculated drag tensor through variations of velocity u , and obtained the orientation averaged drag, but are not reporting those results here. We plan to discuss these, together with alternate methods of calculating (approximating) the average, separately. , for a sample of the biomimetic geometries considered in this study. We selected geometries for inclusion in this table if the number of bodies protruding from the main bulk was eight or less. Note that it appears in the continuum regime, increasing the bumpiness of the particle system increases the drag, but not dramatically so. The symbols G, and Q keep their meanings as the index of the geometry and the number of quadrature points, respectively. Like the results in the last section, one can examine a VES for computations in the continuum range. We have shown the results of this approximation for all sixteen geometries in Table 7 . Note that for the geometries 1-5, 8 and 12 the results of Table 7 are very close to those of Table 6 , and it seems that the VES approximation is quite good for all the sixteen geometries considered here. Thus, further numerical computations for other geometries via the numerical technique do not appear necessary.

J o u r n a l P r e -p r o o f 21 Figure 1 , corresponding to all sixteen geometries used in Table 3 . This table is the product of an indirect computation; first, we computed the drag on a bare sphere of varying radius, in the vicinity of one (to get good resolution around R VES ), the result of which we interpolated to generate the results seen here. Note that in the table G is the geometry index, N is the number of body protrusions, V is the total volume of the geometry, R is the VES radius, and F is the nondimensionalized drag. We note that the original computations used eight quadrature points per sphere. We completed the same calculations with a Q of sixteen and found the results identical to at least four decimal places. Admittedly, we could have obtained the results for the bare sphere from equation (31) 

Since computations over the entire Knudsen number range for particles of arbitrary shape are challenging. Note, Monte Carlo simulations as well as other direct numerical methods have extremely poor convergence in the continuum limit. It is practical to consider an approximate expression for the drag for all Knudsen numbers as derived from a combination of the free molecular and the continuum results. One such expression is based on the harmonic approximation, known as the Sherman's approximation in Rarefied Gas Dynamics: 14-16, 42

Note that the approximation (33) appears similar to a harmonic mean, but is in fact different, as it is meant to be an interpolation, not a statistic computed over a set of values. The subscripts 'fm.' and 'cont.' once again denote the free molecular and the continuum values of the drag, respectively. Defining the socalled Cunningham slip (or correction) factor as the ratio:

We can express the drag force as:

Note that using the free molecular and the continuum drag expression for a single sphere, as in the Sherman approximation, we have (we have displayed below the Kn dependence of the Cunningham Slip Factor):

The Knudsen number is once again: 37) and the mean free path is

We note that different authors distinct definitions for the mean free path, differing by multiplicative factors (such as /2  ), and one needs to take care to use consistent definitions when comparing the work of different authors.

It has been found that the Sherman approximation, (equations (33) - (36) , α equal to of one) for the drag and Cunningham's Slip Factor may differ from (numeric) results obtained from the BGK model or the Boltzmann Equation (again with an accommodation coefficient, α, of one), and Millikan's data by at most ten percent in the transition regime. 14- 16, 42 Since we have reported both the free molecular and the continuum results in Table 3 and Table 8 , respectively, it is straightforward to use the equation (33) for the entire size range by using the appropriate values for a particular configuration of the particle for which the values have been given here; i.e., with the free-molecular drag: J o u r n a l P r e -p r o o f 

That is Kn a will depend on all the factors (shape, orientation) describing the body. We can now estimate the drag on the body using equation (50) .

The ASA approximation is easy to implement, and we have found that both the Sherman and the ASA approximation yield close results. We show results for the Cunningham slip factor in Sherman's approximation and the ASA approximation, in Figure 2 , for a needle like spheroid (we take in Dahneke's nomenclature the parameters of the needle as a = 20, b=1 and , 0   , and f the accommodation coefficient as =0.9). Zhang, et al. (2012) further verified that the ASA approximation gives results that are comparable to the DSMC results for particle aggregates over a range of transition-regime Knudsen numbers. 19 We believe both the Sherman and the ASA approximations are especially useful, but their foundation is heuristic, not theoretical. Neither of them yields correct slip coefficients when examined in the near continuum regime, and their derivatives in the near free molecular regimes do not have logarithmic terms as they should. However, the error induced by these limitations is typically small despite them. Efforts to improve upon these approximations would remain of interest (see Buckley and Loyalka (1989) for an exploration). 42 We have further explored the effect of varying the radius of the small protrusions, while keeping a constant radius of the larger particle main body (we use this radius as the basic length unit, so varying it at the same time is unnecessary). The results of this comparison are shown in Table 11 for continuum computations, respectively. The results show that for small protrusions from a large body, the drag is insensitive to these perturbations, but the effect becomes larger for geometries with more protrusions. We will explore these effects) more in a later paper where we will consider nonspherical protrusions (needle like, among others) of varying sizes and examine the usefulness of the VES approximation further.

We examine the effect of varying the size of the protrusions (shown in the column headings), while keeping the main-body radius constant for three different geometries in the free-molecular regime. Within each case there are three lines, the first line shows the computed nondimensionalized drag, the second line shows the VES drag for the system, and the third line is the ratio of the first two. The results show a slight increase in the error of the VES approximation with the accommodation coefficient, and a variance with respect to the size and number of protrusions that is not simple to quantify. We present a comparison of Sherman's approximation and Dahneke's ASA for the Cunningham Slip Factor for a specific spheroid. Note that for vanishing Knudsen number, (the continuum limit) and for diverging Knudsen number, (the free molecular limit) the difference between the two approximations vanishes. It is only in the transition regime that there is a small, but noticeable difference.

J o u r n a l P r e -p r o o f 30

For the free-molecular computations, we have used here the classic diffuse-specular reflection model of Maxwell (1879) to explore the range of the drag values, but other models based on based on accommodation coefficients are also available. 47 Williams (1973) , had considered various models for simultaneous drag and heat transfer on a sphere, and Chernyak and Sograbi (2019) also reported results for a sphere for models of interest. 18, 48 Then the relationship between the inbound and outbound particle distributions are related by:

In equation (47), the symbol A represents a suitable gas-surface scattering operator. Under these conditions, equation (5) can be evaluated either analytically or numerically. Thus, for example, for the model of Cercignani and Lampis (1971) , Chernyak and Sograbi (assuming α n and α t close to 1), have reported: 39, 48, 49 

Where α n and α t are, respectively, the normal energy and tangential momentum accommodation coefficients. For values of α n and α t not close to one, the integration of equation (16) requires numerical techniques. We computed these integrals and found that the above result is correct even for values of α n and α t close to 0.1 (within 2% or so of values obtained from numerical integrations). We will be reporting these numerical results/comparisons separately.

We should note that all gas-surface models are phenomenological, and most useful values for such coefficients come from experiments with careful control of surface conditions, and the inverse problem (for extracting coefficients and fitting the theoretical results to the data) is not ill-posed. But for most applications, one deals with natural (or engineered) surfaces and particulate surface control is often unavailable.

One can make the problem less ambiguous by assuming relationships between the coefficients. This normally requires (unavailable, or at least inconvenient) experimental data, but datasets do exist in the literature. For example, one may obtain the tangential momentum accommodation coefficient and the coefficient of diffuse specular refection, α, from Millikan's data.

Additional measurements made by Dahneke with an aerosol beam, 44-46 and measurements with levitated and rotating steel sphere in gases such as He, Ar, Kr, Ne, N 2 , CH 4 (and some mixtures of these) also provide experimental data. 41, [50] [51] [52] [53] All these experiments suggest that assuming a value for the accommodation coefficient, α, in its upper range (0.8-0.95) is proper for commonly encountered gasses. For example, the rotating sphere measurements showed α values in the range of 0.88 to 0.97 for He, 0.85 to 0.93 for Ar, 0.80 to 0.95 for Kr, and finally 0.83 to 0.89 for nitrogen, all at STP.

Larriba-Andaluz and Carbone (2021) have recently reviewed values inferred from experiments including with differential mobility analyzers and other means, noting a general range of 0.8-0.9, with the values increasing as the particle size increases. 54 There remain interesting possibilities for analyzing the existing experimental data of Millikan and others (See Rader (1990) for a compilation for several gases), further with the use of the Cercignani-Lampis model in that the two accommodation coefficients in there can be determined by using the free molecular and the slip limits (these will provide two equations). 49, 55 We have noted the free molecular limit result above inequation (26) . Accurate results for the velocity slip coefficient for the Cercignani-Lampis model and general intermolecular interaction laws can be obtained for example, from the variational and similar general expressions given by Loyalka (1971) , and further discussed in Ivchenko, et al. (2007) . 10, 43 Since in the near term at least, there will remain uncertainties with the gas-surface interaction models and parameters, the practical solution is to have theoretical results available for a range of values of these parameters and assess the sensitivity of particle mobility and overall aerosol dispersion/transport and ES-DMA analysis to these variations.

It is also pertinent to note here that in interpretations of experiments, analysts often ignore wall effects. While small in absolute size, the nature of the confined (chamber) geometries may lead to these effects becoming significant. Thus, extensions of the present computations to particle diffusion in confined geometries will be of interest. 

We have discussed measurements of gas-surface interaction parameters above. As we have noted in the introduction, one can measure particle mobility via ES-DMA, and such measurements with well-defined virus like particles would be of considerable value.

There have been measurements on non-spherical particles with improved Millikan cells also, discussed by Cheng, et al. (1988) , and it would again be of interest to measure drag and hence the mobility of virus like particles in such cells. Levitation experiments, with a rotor gauge or similar device, to measure frictional torque on such particles should also be of interest. 56 One could also use electrodynamic balances as described by Davis and Schweiger 57 .

We should also note that in the continuum regime, one can additionally measure drag on particles via settling experiments (in high viscosity fluids), and for common cases, computations can be both verified and confirmed.

J o u r n a l P r e -p r o o f 33

Our aim in this paper was to elucidate the role of virus structure on virus drag, diffusion, and mobility. The methods and the techniques considered also potentially apply to interpretation of virus structures, as well as virus separations through use of ES-DMS techniques and devices. We have considered model virus mimetic geometries (composed of spheres) and obtained free-molecular and continuum results. In the free molecular limit, we have used the Maxwell diffuse-specular reflection model and a Monte Carlo technique, and our results are in good agreement with those obtained from the IMoS program. In the continuum limit, we have solved the boundary integral equation associated with the Stokes equation by using a singularity subtraction and collocation technique and obtained satisfactory results. In both the free molecular and the continuum limits we have found that a volume equivalent sphere approximation gives satisfactory results for the geometries considered in this paper. We have described construction of approximate expressions, using the free molecular and continuum results, which are applicable over the entire size range of the virus size as compared to molecular mean free path (the Knudsen number). We have also discussed gas-surface interaction models and values of the accommodation coefficients reported in the literature. We have further briefly noted the experiments that would be useful in supplying useful data and verification of the theoretical results.

We will be reporting results for more complicated virus structures in the future as shown in Figure 3 below, in aiding understandings of virus and other biological particle transport, both in gases and in liquids. We would be particularly attentive to the applications of this work, and extensions of it, to virus structure determination and separations via the ES-DMA technique.

Shown is a coronavirus mimetic geometry with usage of stochastically placed cones tipped with spheres as computational phantoms for the spike glycoproteins, and small, hemispherical protrusions randomly placed to model the membrane glycoproteins embedded in the capsid surface. Unlike the geometries represented in this study, the location Unlike the geometries represented in this study, we placed the surface features non-deterministically via sampling random matrices, in a procedure we intend to further explore in a later paper. 

We include here a discussion (a bit involved) of the molecular speed sampling from the gamma distribution. There is a generalization of the gamma distribution, fittingly called the "generalized gamma distribution," or (at the US Department of Energy) the Amoroso, (see LLNL (arXiv:1005.3274v2) for a review, or Crooks (Crooks, Gavin E. "Field guide to continuous probability distributions." Berkeley Institute for Theoretical Sciences, Berkeley (2019).) for a uniquely thorough summary of parametric probability distributions) that has historical roots in many problems. The Amoroso has over fifty named special cases in different application areas: particular cases include, but are not limited to, the normal, log-normal, power-law, Weibull, Maxwellian, Nakagami, exponential, chi, chi-square and gamma distributions. The Amoroso has four formal parameters and is defined as: Referring to inbound and outbound particle distributions, respectively. Originally, we sampled these two distributions first by looking at the CDFs:

To be clear regarding notation, in the expression above, the exponent '3' corresponds to the case for F  , taking the cube of the variable  å , and the exponent '4' similarly applies to the F  case. As is traditional, to sample the probability distribution, we take the CDF and invert it with any proper root-finding technique (e.g., Newton-Raphson, etc.):

This is the standard method and is correct. This procedure is slow however, in that root finding algorithms are expensive; the CDF function is costly to compute in addition to constant need to resample added PRNGs on the unit interval. This is when we noticed that the CDFs that we are sampling are distributed according to the square root of the sum of a collection of random squared normal distributions, i.e., a chi (not chi-square) distribution:

    

Which are of course the so-called regularized gamma distributions. At this point, it is important to note a theorem from mathematical statistics, which is the family of gamma distributions is closed under certain arithmetic operations. Specifically, we omit proof of it here. However, examine the random variable scaling transformation:

After some tedious manipulation one finds that f  is identical a Nakagami distribution with both formal parameters equal to 2, and that f  is a Nakagami distribution with both formal parameters equal to 5/2.

In the figures above, the histograms each hold 10k randomly generated variables, the LHS column in each chart was produced over 1k times faster than the direct CDF sampling method on the RHS. The precise details of the sampling technique are restricted. However, we do note the algorithm is hardware dependent, and that if one is using a scaler processor the function will consist of a modified "Ziggurat" algorithm, conversely, if one is operating a vector processor the function will use a "Box Muller" technique. Tables   Table 1: 

We report the nondimensionalized drag, defined via equation (10), for a the case of a linear chain of spheres (i.e., a group of 2spheres embedded in a 3 dimensional Euclidean space), similar to a straight line of croquet or billiard balls touching one another, but not resting on any surface. The results are as expected; note that the drag increases slightly with the accommodation coefficient for a given geometry. benchmarks for the Stokes (continuum) problem for the geometric case of linear chains of touching spheres. In the table, S is the number of spheres in that linear chain, and Q corresponds to 2 Q quadrature points taken per sphere. To ease comparison with other results presented in this paper, we took the radius of the spheres to be one (under our nondimensionalization). For two spheres, the exact result from equation (32) , for a sample of the biomimetic geometries considered in this study. We selected geometries for inclusion in this table if the number of bodies protruding from the main bulk was eight or less. Note that it appears in the continuum regime, increasing the bumpiness of the particle system increases the drag, but not dramatically so. The symbols G, and Q keep their meanings as the index of the geometry and the number of quadrature points, respectively. Figure 1 , corresponding to all sixteen geometries used in Table 3 . This table is the product of an indirect computation; first, we computed the drag on a bare sphere of varying radius, in the vicinity of one (to get good resolution around R VES ), the result of which we interpolated to generate the results seen here. Note that in the table G is the geometry index, N is the number of body protrusions, V is the total volume of the geometry, R is the VES radius, and F is the nondimensionalized drag. We note that the original computations used eight quadrature points per sphere. We completed the same calculations with a Q of sixteen and found the results identical to at least four decimal places. Admittedly, we could have obtained the results for the bare sphere from equation (31) Table 10 : Protrusions of varying size in the free molecular regime We examine the effect of varying the size of the protrusions (shown in the column headings), while keeping the main-body radius constant for three different geometries in the free-molecular regime. Within each case there are three lines, the first line shows the computed nondimensionalized drag, the second line shows the VES drag for the system, and the third line is the ratio of the first two. The results show a slight increase in the error of the VES approximation with the accommodation coefficient, and a variance with respect to the size and number of protrusions that is not simple to quantify. 

We present an image of the geometries used for the Monte Carlo simulations and various calculations throughout the paper. Note that case number one is unshown as it consisted of a single, isolated sphere. Geometries two and three are identical save for the orientation of their single protrusions. In each of the geometries pictured, under the nondimensionalizations used in this paper, the larger, central body-particle has a radius of one. The exception to this would be the geometries used in the construction of Table 3 . The Table 3 geometries have main spheres with a radius of ¾. Primarily, and exceptions are noted where appropriate, the small spheres surrounding the main bodies all have radii of 1/16.

Physico-Chemical Characteristics of Evaporating Respiratory Fluid Droplets

Insights into the Evaporation Characteristics of Saliva Droplets and Aerosols: Levitation Experiments and Numerical Modeling

Parallel Differential Mobility Analysis for Electrostatic Characterization and Manipulation of Nanoparticles and Viruses

Physical Analysis of Virus Particles Using Electrospray Differential Mobility Analysis

Charge Reduced Electrospray Size Spectrometry of Mega-and Gigadalton Complexes: Whole Viruses and Virus Fragments

Physical Characterization of Icosahedral Virus Ultra Structure Stability, and Integrity Using Electrospray Differential Mobility Analysis

Virus Size Analysis by Gas-Phase Mobility Measurements: Resolution Limits

Exceeding a Resolving Power of 50 for Virus Size Determination by Differential Mobility Analysis

Virus-Like Particle Size and Molecular Weight/Mass Determination Applying Gas-Phase Electrophoresis

Aerosol Science: Theory and Practice: With Special Applications to the Nuclear Industry

Fundamentals of Aerosol Dynamics

The Mechanics of Aerosols

Flow of a Rarefied Gas Past an Axisymmetric Body. II. Case of a Sphere

Motion of a Sphere in a Gas: Numerical Solution of the Linearized Boltzmann Equation

Numerical Analysis of a Uniform Flow of a Rarefied Gas Past a Sphere on the Basis of the Boltzmann Equation for Hard-Sphere Molecules

On the Motion of Small Spheres in Gases III. Drag and Heat Transfer

Free-Molecular Drag on Straight Chains of Uniform Spheres

Monte Carlo Simulation of Hydrodynamic Drag and Thermophoresis of Fractal Aggregates of Spheres in the Free-Molecular Flow Regime

Determination of the Scalar Friction Factor for Nonspherical Particles and Aggregates across the Entire Knudsen Number Range by Direct Simulation Monte Carlo (DSMC)

Molecular Gas Dynamics and the Direct Simulation of Gas Flows

Computation of Fission Product Condensation on Chainlike Aerosols and Agglomerates--II: Role of Energy and Mass Dependance of Molecules

Computation of Fission Product Condensation on Chainlike Aerosols and Agglomerates

Coronavirus Envelope Protein: Current Knowledge

Free Molecular Collision Cross Section Calculation Methods for Nanoparticles and Complex Ions with Energy Accommodation

Ion Mobilities in Diatomic Gases: Measurement Versus Prediction with Non-Specular Scattering Models

Collision Cross Section Calculations for Polyatomic Ions Considering Rotating Diatomic/Linear Gas Molecules

A Parallelized Tool to Calculate the Electrical Mobility of Charged Aerosol Nanoparticles and Ions in the Gas Phase

Gaseous Ion Mobility and Diffusion in Electric Fields of Arbitrary Strength

The Mathematical Theory of Non-Uniform Gases

Chapman-Enskog Solution for Diffusion: Pidduck's Equation for Arbitrary Mass Ratio

Monte Carlo Geometry Processing: A Grid-Free Approach to PDE-Based Methods on Volumetric Domains

Advanced Transport Phenomena: Fluid Mechanics and Convective Transport Processes

Boundary Integral and Singularity Methods for Linearized Viscous Flow

Stokes Flow Past a Particle of Arbitrary Shape: Numerical Method of Solution

Condensation on Nonspherical Aerosol Particles: Numerical Solutions in the Continuum Regime

Condensation/Evaporation on Arbitrarily Shaped Aerosols

Theory of the spinning rotor gauge in the slip regime

Boundary Element Method Calculations of the Mobility of Nonspherical Particles 1. Linear Chains

Cunningham Correction Factor and Accomodation Coefficient: Interpretation of Millikan's Data

Slip Correction Factors for Nonspherical Bodies -III. The Form of the General Law

Slip Correction Factors for Nonspherical Bodies -I. Introduction and Continuum Flow

The Role of Molecule-Surface Interaction in Thermophoresis of an Aerosol Particle

The Spinning Rotor Gauge: Measurements of Viscosity, Velocity Slip Coefficients, and Tangential Momentum Accommodation Coefficients

The Spinning Rotor Gauge: Measurements of Viscosity, Velocity Slip Coefficients, and Tangential Momentum Accomodation Coefficients for N 2 and Ch 4

Viscosity and Velocity Slip Coefficients for Gas Mixtures: Measurements with a Spinning Rotor Gauge

Measurements of Viscosity, Velocity Slip Coefficients, and Tangential Momentum Accomodation Coefficients Using a Modified Spinning Rotor Gauge

The authors would like to thank the editor, and two anonymous reviewers for their helpful suggestions. 

Results of the drag Monte Carlo simulation across five different coefficients of diffuse/specular reflection, for sixteen different geometries meant to mimic the novel coronavirus SARS-CoV-2. Please note that the simulation index, N, is arbitrary, and the sphere number, S, is the number of protrusions from the main body. Within each geometry set, the top lines are the results of our Monte Carlo Simulations, the second lines are the drag experienced by a volume equivalent sphere, and the third lines are the ratios of the first two, respectively. The drag for the volume equivalent sphere cases was computed with the aid of equations Drag computed via equations (5) - (11) . The results show that for each simulation geometry, an increase in the quantity of the capsid surface corresponds to an increase in the drag experienced by the particle. The uncertainty in the drag, which as the result of the Monte Carlo simulation, is the less precise than values obtained via evaluation of analytic formulae, drives the uncertainty of the calculations in general. Table 3 , the primary difference being that we scaled down the spheres used to model the particles by a factor of ¾, to ease some comparison on the relationship between corona size and drag. Please note that we did not duplicate experiment one replicated, as it consisted of a single bare sphere. Otherwise, the descriptions from Table 3 are still applicable, as well as the pictures in panel graphic Figure 1 . 

Shown is a coronavirus mimetic geometry with usage of stochastically placed cones tipped with spheres as computational phantoms for the spike glycoproteins, and small, hemispherical protrusions randomly placed to model the membrane glycoproteins embedded in the capsid surface. Unlike the geometries represented in this study, the location Unlike the geometries represented in this study, we placed the surface features non-deterministically via sampling random matrices, in a procedure we intend to further explore in a later paper.J o u r n a l P r e -p r o o f

Declarations of interest: none.To the best of our knowledge, we have no conflicts of interest to report.