The COVID-19 pandemic has made clear the fundamental role of airborne droplets and aerosols as potential virus carriers. The importance of studying the fluid dynamics of exhalations, starting from the formation of droplets in the respiratory tracts to their evolution and transport as a turbulent cloud, can now be recognized as the key step towards understanding SARS-CoV-2 transmission. Respiratory droplets are formed and emitted at high speed during a sneeze or cough [97] , and at lower speed while talking or breathing. The virus-laden droplets are then initially transported as part of the coherent gas puff of buoyant fluid ejected by the infected host [13] . The very large drops of O(mm) in size, which are visible to naked eye, are minimally affected by the puff.
They travel semi-ballistically with only minimal drag adjustment, but rapidly falling down due to gravitational pull. They can exit the puff either by overshooting or by falling out of the puff at the early stage of emission (Fig. 1 ). Smaller droplets ( O(100 µm)) that remain suspended within the puff are advected forward. As the suspended droplets steadily evaporate within the cloud, the virus takes the form of potentially inhalable droplets, or droplet residues when the evaporation of water is complete. Meanwhile, the velocity of the turbulent puff continues to decay both due to entrainment and drag. Once the puff slows down sufficiently, and its coherence is lost, the eventual spreading of the virus-laden droplet residue becomes dependent on the ambient air currents and turbulence.
The isolated respiratory droplet emission framework was introduced by Wells [115] in the 1930s and remains the framework used for guidelines by public health agencies, such as the WHO, CDC and others. However, it does not consider the role of the turbulent gas puff within which the droplets are embedded. Regardless of their size and their initial velocity, the ejected droplets are subject to both gravitational settling and evaporation [13] . Although droplets of all sizes undergo continuous settling, droplets with settling speed smaller than the fluctuating velocity of the surrounding puff can remain trapped longer within the puff (Fig. 2) . Furthermore, the water content of the droplets continuously decreases due to evaporation. When conditions are appropriate for near complete evaporation, the ejected droplets quickly become droplet residues of non-volatile biological material. The settling velocity of these droplet residues is sufficiently small that they can remain trapped as a cloud and get advected by ambient air currents and dispersed by ambient turbulence. Based on the above discussion, we introduce the following terminology that will be consistently used in this paper:
• Puff: Hot, moist air exhaled during breathing, talking, coughing or sneezing, which remains coherent and moves forward during early times after exhalation
• Cloud: The distribution of ejected droplets that remain suspended even after the puff has lost its coherence. The cloud is advected by the air currents and is dispersed by ambient turbulence
• Exited droplets: droplets that have either overshot the puff/cloud or settled down due to gravity
• Airborne (evaporating) droplets: droplets which have not completed evaporation and retained within the puff/cloud
• (Airborne) droplet residues: droplets that remain airborne within the puff/cloud and that have fully evaporated, which will also be termed aerosols.
Host-to-host transmission of virus-laden droplets and droplet residues occurs generally through direct and indirect routes [4, 12, 84] . The direct route of transmission involves the larger droplets that may ballistically reach the recipient's mucosa. This route is currently thought to involve either the airborne route or drops that have settled on surfaces. The settled drops remain infectious, to be later picked up by the recipient, and are generally thought to be localized to the vicinity or at close range of the original infectious emitter. With increased awareness and modified physical distancing norms it is possible to minimize the spreading of the virus by such direct route.
The indirect route of transmission is one that does not necessarily involve a direct or close interaction between the infectious individual and the recipient or for the two to be synchronously Figure 1 : Image reproduction showing the semi-ballistic largest drops, visible to the naked eye, and on the order of mm, which can overshoot the puff at its early stage of emission [14, 15] . The puff continues to propagate and entrain ambient air as it moves forward, carrying its payload of continuum of drops [13] , over distances up to 8 meters for violent exhalations such as sneezes [17] . present in the same contaminated space at the same time. Thus, the indirect route involves respiratory droplets and fully-evaporated droplet residues that are released to the surrounding by the infected individual, which remain airborne as the cloud carries them over longer distances. The settling speeds of the airborne droplets and droplet residues are so small, that they remain afloat for longer times, while being carried by the background turbulent airflow over distances that can span the entire room or even multiple rooms within the building (O(10 − 100) feet). A schematic of the two routes of transmission is shown in Fig. 2 and in this paper we will focus on the indirect airborne transmission.
Another factor of great importance is the possibility of updraft in the region of contamination, due to buoyancy of the virus-laden warm ejected air-mass. These slight updrafts can keep the virusladen droplets suspended in the air and enhance the inhalability of airborne droplets and droplet residues by recipients who are located farther way. The advection of airborne droplets and residues by the puff and subsequently as a cloud may represent transmission risk for times and distances much longer than otherwise previously estimated, and this is a cause of great concern [99, 100] . Note that if we ignore the motion of the puff of air carrying the droplets, as in the analysis of Wells, the airborne droplets and residues would be subjected to such high drag that they could not propagate more than a few cm away from the exhaler, even under conditions of fast ejections, such as in a sneeze. This illustrates the importance of incorporating the correct multiphase flow physics in the modeling of respiratory emissions [16] , which we shall discuss further here.
It has been recently reported that COVID-19 virus lives in droplets and aerosols for many hours in laboratory experiments [25] . At the receiving end, an increased concentration of virusladen airborne droplets and residues near the breathing zone increases the probability of them settling on the body or more importantly being inhaled. Depending on its material and sealing properties, the use of a mask by the infected host can help reduce the number of virus-laden droplets ejected into the air, and in a less effective way, the use of a mask or other protective devices by the receiving host may reduce the probability of inhalation of the virus-laden airborne droplets and residues.
The above description provides a clear sketch of the sequence of processes by which the virus is transferred host-to-host. This simplistic scenario, though pictorially evocative, is tremendously insufficient to provide science-based social distancing guidelines and recommendations. There is substantial variability (i) in the quantity and quality of contaminated droplets and aerosols generated by an infected person, (ii) in the manner by which the contaminated droplets and droplet residues remain afloat over longer distances and time, (iii) in the possibility of the contaminant being inhaled by a recipient and (iv) in the effectiveness of masks and other protection devices. Violent exhalations, such as sneezing and coughing, yield many more virus-laden droplets and aerosols than breathing and talking [13, 73] . All coughing and sneezing events are not alike -the formation of droplets by break up of mucus and saliva varies substantially between individuals. Significant variation in initial droplet size and velocity distribution has been reported in [3, 13, 44, 84] . The measured droplet size distribution, particularly for transient biological emissions such as respiratory exhalations, also depends on ambient temperature and humidity and on the methodology and instrumentation used to characterize the size distribution [4, 16, 73] . Furthermore, it is of importance to consider the volume of air, and the pathogen load, being inhaled during breathing by the receiving host. Thus, there is great variability in how much of the virus-laden aerosols reach from the infected host to the receiving host.
CDC guideline of social distancing of 2 meters (6 feet) is based on the disease transmission theory originally developed in 1930s and later improved by others [90, 114, 117] . The current recommendation of 6 feet as the safe distance can be improved in several ways: (i) by accurately accounting for the distance traveled by the puff and the droplets contained within it, while some continuously settling out of the puff, (ii) by accurately evaluating the evaporation of droplets and the subsequent advection and dispersal of droplet residues as a cloud [6] , (iii) by incorporating the effect of adverse flow conditions that prevail under confined indoor environment including elevators, aircraft cabins, and public transit, or favorable conditions of open space with good breeze or cross ventilation, and (iv) by correctly assessing the effectiveness of masks and other protective devices [20] . Thus, mechanistic, evidence-based understanding of exhalation and dispersal of expelled respiratory droplets, and their subsequent fate as droplet residues in varying scenarios and environments is important. We must therefore revisit the safety guidelines and update them to modern understanding. In particular, a multi-layered guideline that differentiates between crowed class rooms, auditoriums, buses, elevators and aircraft cabins from open outdoor cafes is desired. Only through reliable understanding of the underlying flow physics of virus transmission one can arrive at such nuanced guidance in designing case-specific social distancing guidelines.
The object of the paper is to present a coherent analytic and quantitative description of droplets generation, transport, conversion to droplet residues, and eventual inhalation. We will examine the available quantitative relationships that describe the above processes and adapt them to the present problem. The key outcomes that we desire are (i) A simple universal description of the initial droplet size spectrum generated by sneezing, coughing, talking and breathing activities. Such a description must recognize the current limitations of measurements of droplet size distribution under highly transient conditions of respiratory events. (ii) A first-order mathematical framework that describes the evolution of the cloud of respiratory droplets and their conversion to droplet residues, as a function of time, and (iii) A simple description of the inhalability of the aerosols along with corresponding evaluation of the effectiveness of different masks based on existing data reported to date. The physical picture and the quantitative results to be presented can then be used to study a statistical sample of different scenarios and derive case-specific guidelines. We anticipate the present paper to spawn future research in the context of host-to-host airborne transmission.
After presenting the mathematical framework in section 2, the three different stages of transmission, namely droplet generation, transport and inhalation will be independently analyzed in sections 3, 4 and 5. These sections will consider the evolution of the puff of exhaled air and the droplets contained within. Section 6 will put together the different models of the puff and droplet evolution described in the previous sections, underline their simplifications, and demonstrate their ability to make useful prediction. Finally, conclusions and future perspectives are offer in section 7. A schematic representation of size distribution at the infected source host, at an intermediate distance and at a receiving host located farther away is also shown.
We wish to describe the three main stages involved in the host-to-host transmission of the virus: droplet generation during exhalation, airborne transport, and inhalation by the receiving host. In the generation stage, virus-laden drops are generated throughout the respiratory tract by the exhalation air flow which carries them through the upper airway toward the mouth where they are ejected along with the turbulent puff of air from the lungs. The ejected puff of air can be characterized with the following four parameters: the volume Q pe , the momentum M pe , and the buoyancy B pe of the ejected puff, along with the angle θ e to the horizontal at which the puff is initially ejected. The initial momentum and buoyancy of the puff are given by M pe = ρ pe Q pe v pe and B pe = (ρ a − ρ pe )Q pe g, where v pe is the initial velocity of ejected puff, ρ pe and ρ a are the initial density of the puff and the ambient, respectively, and g is the gravitational acceleration. The ejected droplets are characterized by their total number N e , size distribution N e (d), droplet velocity distribution V de (d) and droplet temperature distribution T de (d), where d is the diameter of the droplet. To simplify the theoretical formulation, here we assume the velocity and temperature of the ejected droplets to depend only on the diameter and show no other variation. As we shall see in section 4, this assumption is not very restrictive, since the velocity and temperature of the droplets that remain within the puff very quickly adjust to those of the puff. Both the ejected puff of air and the detailed distribution of droplets depend on the nature of the exhalation event (i.e., breathing, talking, coughing or sneezing), and also on the individual. This is followed by the transport stage, where the initially ejected puff of air and droplets are transported away from the source. The volume of the puff of air increases due to entrainment of ambient air. The puff velocity decreases due to both entrainment as well as drag. Since the temperature and moisture content of the ejected puff of air is typically higher than the ambient, the puff is also subjected to a vertical buoyancy force, which alters its trajectory from a rectilinear motion. The exhaled puff is turbulent and both the turbulent velocity fluctuations within the puff and the mean forward velocity of the puff decay over time.
The time evolution of the puff during the transport stage can then be characterized by the following quantities: the volume Q p (t), the momentum M p (t), buoyancy B p (t) of the ejected puff, and ρ p (t) is the density of air within the puff which changes over time due to entrainment and evaporation. The trajectory of the puff is defined in term of the distance traveled s(t) and the angle to the horizontal θ(t) of its current trajectory. Following the work of Bourouiba et al. [13] we have chosen to describe the puff trajectory in terms of s(t) and θ(t). This information can be converted to horizontal and vertical positions of the centroid of the puff as a function time. If we ignore the effects of thermal diffusion and ambient stratification between the puff and the surrounding air, then the buoyancy of the puff remains a constant as B p (t) = B pe . Furthermore, as will be seen below, the buoyancy effects are quite weak in the early stages when the puff remains coherent and thus the puff to good approximation can be taken to travel along a straight line path, as long as other external flow effects are unimportant.
To characterize the time evolution of the virus-laden droplets during the transport stage we distinguish the droplets that remain within the puff, whose diameter is less than a cutoff (i.e., d < d exit ), from the droplets (i.e., d > d exit ) that escape out of the puff. As will be discussed subsequently in §4, the cutoff droplet size d exit decreases with time -at very early times it decreases as t −1/4 when the ejection behaves as a jet, then decays as t −1/2 when evolving as a puff and at later times the cutoff diameter decreases as t −3/8 . Thus, the total number of droplets that remain within the puff can be estimated as N (t) = d exit 0 N (d, t) dd. However, the size distribution of droplets at any later time, denoted as N (d, t), is not the same as that at ejection. Due to evaporation, size distribution shifts to smaller diameters over time. We introduce the mapping D(d e , t), which gives the current diameter of a droplet initially ejected as a droplet of diameter d e . Then, assuming well-mixed condition within the puff, the airborne droplet and residue concentration (number per volume) distribution can be expressed as
where the inverse mapping D −1 gives the original ejected diameter of a droplet whose current size is d. The prefactor 1/Q p (t) accounts for the decrease in concentration due to the enlargement of the puff over time. In this model the airborne droplets and residue that remain within the coherent puff are assumed to be in equilibrium with the turbulent flow within the puff. Under this assumption, the velocity V d (d, t) and temperature T d (d, t) of the droplets can be estimated with the equilibrium Eulerian approximation [32, 34] . When the puff's mean and fluctuating velocities fall below those of the ambient, the puff can be taken to loose its coherence. Thus, the puff remains coherent and travels farther in a confined relatively quiescent environment, such as an elevator, class room or aircraft cabin, than in an open outdoor environment with cross-wind or in a room with strong ventilation. We define a transition time t tr , below which the puff is taken to be coherent and the above described puff-based transport model applies. For t > t tr , we take the aerosol transport and dilution to be dominated by ambient turbulent dispersion. Accordingly, this late-time behavior of total number of airborne droplets and residue and their number density distribution are given by theory of turbulent dispersion. It should be noted that the value of transition time will depend on both the puff properties as well as the level of ambient turbulence (see section 4.4) .
We now consider the final inhalation stage. Depending on the location of the recipient host relative to that of the infected host, the recipient may be subjected to either the puff that still remains coherent, carrying a relatively high concentration of virus-laden droplets or residues, or to the more dilute dispersion of droplet residues, or aerosols. These factors determine the number and size distribution of virus-laden airborne droplets and residue the recipient host will be subjected to. The inhalation cycle of the recipient, along with the use of mask and other protective devices will then dictate the aerosols that reach sensitive areas of the respiratory tract where infection can occur. Following the above outlined mathematical framework we will now consider the three stages of generation, transport and inhalation.
Knowing the droplet sizes, velocities and ejection angles resulting from an exhalation is the key first step in the development of a predictive ability for droplet dispersion and evolution. Respiratory droplet size distributions have been the object of a large number of studies, as reviewed in [39] , and among them, those of Duguid [26] and Loudon & Roberts [69] have received particular scrutiny as a basis for studies of disease transmission by Nicas, Nazaroff & Hubbard [80] . There are substantial differences in the methodologies used for quantification of respiratory emission sprays. Few studies have used common instrumentation that have enough overlap to reconstruct the full distribution of sizes. For example, there are important gaps in reporting the total volume or duration of air sampling, in addition there are issues in reporting the effective evaporation rates used to back-compute the initial distribution and in the documentation of assumptions about optical or shape properties of the droplets being sampled. In addition, sensitivity analyses are often missing regarding the role of orientation or calibration of sensing instruments with respect to highly variable emissions from human subjects. Finally, regarding direct high-speed imaging methods [5, 97] , the tools for precise quantification of complex unsteady fragmentation and atomization processes are only now being developed [16, 110, 112] . There are far fewer studies on the velocities and angles of the droplets produced by atomizing flows.
The studies of Duguid and Loudon & Roberts have been performed by allowing the exhaled droplets to impact various sheets or slides, with different procedures being used for droplets smaller than 20 µm. The size of the stains on the sheets was observed and the original droplet size was inferred from the size of the stains. To account for the difference between the droplet and the stain sizes an arbitrary factor is applied and droplets smaller than 10 or 20 microns are processed differently than larger droplets. The whole process makes the determination of the number of droplets smaller than 10 microns less reliable. The data are replotted in Fig. 3 .
Many authors have attempted to fit the data with a log-normal probability distribution function. In that case, the number of droplets between diameter d and d + dd is N e (d) dd, and the frequency of ejected droplet size distribution is given by
where dd is a relatively small diameter increment or bin width, B is a normalization constant,μ is the expected value of ln d, also called the geometric mean andσ is the standard deviation of ln d, also called the geometric standard deviation (GSD).
On the other hand, there have been also numerous studies of the fragmentation of liquid masses in various physical configurations other than the exhalation of mucosalivary fluid [38, 97, 108, 111] . These configurations include spray formation on wave crests [106] , droplet impacts on solids and liquids [76] , wave impacts on vertical or finite walls/surfaces [60, 110, 113] , and jet atomization [29] . These studies reveal a number of qualitative similarities between the various processes, which can be best described as a cascade of events. Those events include a primary instability of sheared layers in high speed air flows [35] , and then the non linear growth of the perturbation into thin liquid sheets. The sheets themselves may be destabilized by two routes, one involving the formation of Taylor-Culick end rims [21, 103] , and their subsequent deformation into detaching droplets [110] . The other route to the formation of droplets is the formation of holes in the thin sheets [62, 83, 97] . The holes then expand and form free hanging ligaments which fragment into droplets through the Rayleigh-Plateau instability [29] .
Considering the apparent universality of the process, one may infer that a universal distribution of droplet sizes may exist. Indeed, the log-normal distribution has often been fitted to experimental [71] and numerical data on jet formation [47, 65] , for droplet impacts on solid surfaces [116] , and for wave impacts on solid walls [116] . The log-normal distribution is frequently suggested for exhalations [80, 115] . The fit of the numerical results of [65] is shown in Fig. 4 . However, this apparent universality of the log-normal distribution is questionable for several reasons. First, many other distributions, such as exponential, Poisson, Weibull-Rosin-Rammler, beta, or families of gamma or compound gamma distributions [61, 107] capture to some extent the complexity of atomization physics. Second, the geometrical standard deviation (GSD) of the log-normal fits to the many experimental measurements is relatively small (of the order of 1.2 [65] or 1.8 [71] ) while the wide range of scales in Fig. 3 seems to indicate a much larger GSD. Indeed Nicas, Nazarroff & Roberts [80] obtainσ 8−9. One explanation for the smaller GSD in jet atomization studies, both numerical and experimental, is that the numerical or optical resolution is limited at the small scales. Indeed, as grid resolution is increased, the observed GSD also increases [65] . Third, many authors observe multimodal or bimodal distributions, that can be obtained for example by the superposition of several physical processes. This would arise in a very simple manner if the Taylor-Culick rim route produced drops of a markedly different size than the holes-in-film route.
In order to elucidate this discrepancy, we take another look at the fit of the Duguid data in Fig. 5 . We replot the data that was provided in Table 3 of Duguid. Since the data are given as Fig. 5 , since if plotted in the variables x = ln d and y = ln[dN e (d)] the distribution (2) appears as a parabola. When one attempts to fit a parabola between 2 and 50 µm, one obtains a log-normal distribution withσ = 0.7 andμ = ln(12) (for diameters in microns). However the Frequency N e (d) (1/microns)
Duguid cough data Loudon and Roberts cough data Pareto B/d 2 fit Figure 3 : Frequency of droplet size distribution, replotted from Duguid [26] and Loudon & Roberts [69] . The Pareto distribution is also plotted.
data above 50 µm are completely outside this distribution. If instead the whole range from 1 to 1000 µm is fit to a log-normal distribution, one obtains a very wide log-normal or alternatively a Pareto distribution of power 2
In Figs. 3 and 5 we represent the Pareto distribution together with the Duguid and Loudon & Roberts data. It is especially clear from Fig. 3 that if one does not trust either data at d < 50 µm then both data sets are well described by the Pareto distribution. This however does not eliminate the possibility that more data with more statistical power could show deviations from Pareto, in particular as multimodal distributions. Nevertheless, the multimodal deviation from the Pareto distribution is difficult to characterize and will not be pursued in what follows for the sake of simplicity.
It is clear that the Pareto distribution cannot be valid at diameters that are either too large or too small. The equivalent diameter of the total mass of liquid being atomized is an obvious upper bound, but it is also very unlikely that droplets with d > h where h is the initial film thickness will be observed. It is reasonable to put this film thickness on the scale of 1 mm which corresponds to the upper bound on diameters in the data of Figs. 3 and 5. The lower bound on droplet diameter is much harder to determine. Exhalations are highly transient, or unsteady, processes involving complex multiscale geometry [97] , and thread breakup is a fractal multiscale process with satellite Fig. 4 . The fit is adequate only up to 50 µm. As a result, only a fraction of the reliable data fits the log-normal. The Pareto distribution is a reasonable capture of the data in the 10 to 1000 µm range. In the log-log coordinates the log-normal distribution appears as a parabola while the Pareto distribution is a straight line. droplets [28, 104] . Going down in scale, the fractal process repeats itself as long as continuum mechanics remains valid, to around 1 nm. This would not be relevant for viral disease propagation as a lot of the relevant viruses have sizes ranging from O(10 − 100 nm), with an estimated size for SARS-CoV-2 ranging from 60-120 nm, for example. If the smallest length scale is the thickness at which the thin liquid sheets will break, then experimental observations in water [83] suggest a scale of O(100) nm. Other fluids, including biological fluids or biologically contaminated fluids such as those investigated in [87, 88, 111] may yield different length scales.
Based on the above considerations, we take a histogram of droplet sizes that reads
where d 1 is set to O(100 nm) and d 2 to O(1 mm) for simplicity. The total volume of the droplets is
Since d 1 is four orders of magnitude smaller than d 2 , the total number of droplets is well approximated by
, the number of droplets with diameter smaller than x, is very well approximated by
so that f (10d 1 )/N e = 90% of the droplets are of size less that 10 d 1 1 µm. In other words, a numerical majority of the droplets are near the lower diameter bound. On the other hand, a majority of the volume of fluid is in the larger droplet diameters.
The distribution of velocities and ejection angles has been investigated in the atomization experiments of [23] which follow approximately the geometry of a high speed stream peeling by a gas layer. These experiments were qualitatively reproduced in the numerical simulations of [66] . To cite ref. [23] , "most of the ejection angles are in the range 0 • to 40 • , however it occurs occasionally that the drops are ejected with angles as high as 60 • ".
On the other hand, there are to our knowledge no experimental data on the velocity of droplets, as they are formed in an atomizing jet, that could be used directly to estimate the ejection speed of droplets in exhalation. There are however numerical studies [11, 53] in the limit of very large Reynolds and Weber number. The group velocity of waves formed on a liquid layer below a gas stream has been estimated by Dimotakis [24] as
where ρ d is droplet density. In [11, 53] it was shown that this was also the vertical velocity of the interface perturbation. It is thus likely that this velocity plays a role at the end of the first instability stage of atomization. After this stage, droplets are detached and immersed in a gas stream of initial ejection velocity v pe . Since the density ratio ρ p /ρ d is O(10 −3 ), we expect the initial velocity of the ejected droplets at the point of their formation to be small. As we show below, it is interesting to note that the large Reynolds number limit may apply at the initial injection stage to a wide range of droplets in the spectrum of sizes found above. Indeed the ejection Reynolds number of a droplet ejected at a velocity V de in a surrounding air flow of velocity v pe is
where ν a is the kinematic viscosity of the ejected puff of air (here taken to be the same as that of the ambient air). The largest Reynolds number is obtained for the upper bound of d = 1 mm. For example, if the droplet's initial velocity is set to V de ≈ 0, and the air flow velocity in some experiments [17] is as high as 30 m/s, we can estimate the largest ejection Reynolds number to be Re e ≈ 2000 and the Reynolds number will stay above unity for droplets down to micron size. But as the puff of air and the droplets move forward, the droplet Reynolds number rapidly decreases for the following reasons: (i) as will be seen in section 4.1 the puff velocity decreases due to entrainment and drag, (ii) as will be seen in section 4.2.1 the droplet diameter will decrease rapidly due to evaporation, (iii) as will be seen in section 4.2.2 the time scale τ V on which the droplet accelerates to the surrounding fluid velocity of the puff is quite small, and (iv) very large droplets quickly fall out of the puff and do not form part of airborne droplets. Thus, it can be established that droplets smaller than 100 µm quickly equilibrate with the puff within the first few cm after exhalation.
This section will consider the evolution of the puff of hot moist air with the droplets after their initial ejection. First in section 4.1 we will present a simple modified model for the evolution of the puff of exhaled air, evaluating the effects of drag and the inertia of the droplets within it. This will enable us, in section 4.2 to discuss the evolution of the droplet size spectrum, velocity and temperature distributions, with simple first order models. Additionally, section 4.3 will discuss the effect of non-volatiles on the droplet evolution and the formation of a fully evaporated droplet residue or aerosol particle. Late-time turbulent dispersion of the virus-laden droplet residues, when the puff of air within which they are contained stops being a coherent entity, is then addressed in section 4.4.
For the puff model we follow the approach of Bourouiba et al. [13] , but include the added effects of drag and the mass of the injected droplets. In addition, a perturbation approach is pursued to obtain a simple solution with all the added effects included. Fig. 6 shows the evolution of the puff along with quantities that define the puff [13] . We define t to be the time elapsed from exhalation and s(t) to be distance traveled by the puff since exhalation. For analytical considerations we define the virtual origin to be at a distance s e from the real source in the backward direction and t e to be the time it takes for the puff to travel from the virtual origin to the real source. We define t = t+t e to be time from the virtual origin and s = s + s e to the distance traveled from the virtual origintheir introduction simplifies the analysis.
From the theory of jets, plumes, puffs and thermals [105] the volume of the puff exhaled grows by entrainment. Bourouiba et al. [13] defined the puff to be spheroidal in shape with the transverse dimension to evolve in a self-similar manner as r (t ) = αs (t ), where α is related to entrainment coefficient. The volume of the puff is then Q p (t ) = ηr 3 (t ) = ηα 3 s 3 (t ) and the projected, or cross-sectional, area of the puff a(t ) = βr 2 (t ) = βα 2 s 2 (t ), where the constants η and β depend on the shape of the spheroid. For a spherical puff η = 4π/3 and β = π.
As defined earlier, the ejected puff at the real source (i.e., at t = t e ) is characterized by the volume Q pe = ηα 3 s 3 e , momentum M pe = ρ pe Q pe v pe , buoyancy B pe = Q pe (ρ a − ρ pe )g and ejection angle θ e . From the assumption of self-similar growth we obtain the virtual origin to be defined as
where the constant C depends on the drag coefficient of the puff and will be defined below. If we assume a spherical puff with an entrainment factor α = 0.1 [105] , the distance s e depends only on the ejected volume. Experimental measurements suggest Q pe to vary over the range 0.00025 to 0.0025 m 3 . Accordingly, s e can vary from 0.39 to 0.84 m. Similar estimates of t e can be obtained for a spherical puff: as Q pe varies from 0.00025 to 0.0025 m 3 and as the ejected velocity varies from 1 to 10 m/s the value of t e varies over the range 0.01 to 0.21 s. The horizontal and vertical momentum balances in dimensional terms are
In the above C D is the drag coefficient of the puff and M d is the momentum of droplets within the puff. While the puff velocity decreases rapidly over time, the velocity of the larger droplets will change slowly. Note that in the analysis to follow we take the velocity of those droplets that remain within the puff to be the same as the puff velocity. We use s e and t e as the length and time scales to define nondimensional quantities:s = s /s e andt = t /t e . With this definition the virtual origin becomest = 0 ands = 0 and the real source becomest = 1 ands = 1. In terms of non-dimensional quantities the governing momentum equations can be rewritten as
There are three nondimensional parameters: mass ratio of the initial ejected droplets to the initial air puff: r m = ρ d Q de /(ρ p Q pe ); the scaled drag coefficient: C = C D β/(2ηα); and the buoyancy parameter: A = B pe t 2 e /(ρ pe Q pe s e ). In the above equations, r m is defined in terms of the mass of Figure 6 : Evolution of a typical cloud of respiratory multiphase turbulent droplet-laden air following breathing, talking, coughing and sneezing activities. Image adapted from [13] .
the initial ejected droplets. This is an approximation since some of the droplets exit the puff over time. Even though the droplet mass decreases due to evaporation, the associated momentum is not lost from the system since it remains within the puff. In any case, soon it will be shown that the value of r m is small and the role of ejected droplets on the momentum balance is negligible. It should also be noted that under Boussinesq approximation the small difference in density between the puff and that the ambient is important only in the buoyancy term. For all other purposes the two will be taken to be the same and as a result the time variation of puff density is not of importance (i.e., ρ p = ρ pe = ρ a ).
The importance of inertia of the ejected droplets, drag on the puff and buoyancy effects can now be evaluated in terms of the magnitude of the nondimensional parameters. Typical experimental measurements of breathing, talking, coughing and sneezing indicate that the value of r m is smaller than 0.1 and often much smaller. Furthermore, as droplets fall out continuously [13] from the turbulent puff, this ratio changes over time. Here we will obtain an upper bound on the inertial effect of injected droplets by taking the value of r m to be 0.1.
The drag coefficient of a spherical puff of air is also typically small -again as an upper bound we take C D = 0.1, which yields C = 0.375 for a spherical puff. The value of the buoyancy parameter A depends on the density difference between the ejected puff of air and the ambient, which in turn depends on the temperature difference. For the entire range of ejected volumes and velocities, the value of A comes to be smaller than 0.01, for temperature differences of the order of ten to twenty degrees between the exhaled puff and the ambient.
Since all three parameters r m , C and A can be considered as small perturbations, the governing equations can be readily solved in their absence to obtain the following classical expressions for the nondimensional puff location and puff velocity:
With the inclusion of the drag term the governing equations become nonlinear. Nevertheless, they allow a simple exact solution which can be expressed as
Thus, as to be expected, the forward propagation of the puff slows down with increasing nondimensional drag parameter C. For small values of C the above can be expanded in Taylor series as
A comparison of the exact solution with the above asymptotic expansion shows its adequacy for small values of C.
For small non-zero values of r m , C and A, the governing equations can be solved using regular perturbation. The result can be expressed as
and the above expression is accurate to O(C 3 , r 2 m , A 2 ). Although the effect of buoyancy is to curve the trajectory of the puff, the leading order effect of buoyancy is to only alter the speed of rectilinear motion. Also, as expected, the effect of non-zero r m is to add to the total inertia and thereby slow down the motion of the puff. On the other hand, the effect of buoyancy is to slow down if the initial ejection is angled down (i.e., if θ e < 0) and to speed up if the ejection is angled up, provided the ejected puff is warmer than the ambient.
The time evolution of the puff as predicted by the above analytical expression is shown in Fig. 7 . Note that the point of ejection is given byt = 1,s = 1, and the initial non-dimensional velocityṽ(t = 1) = 1/4. The results for four different combinations of C and r m are shown. The buoyancy parameter has very little effect on the results and therefore is not shown. It should be noted that at late stages when the puff velocity slows down the effect of buoyancy can start to play a role as indicated in experiments and simulations. It can be seen that the effect of inertia of the ejected droplets, even with the upper bound of holding their mass constant at the initial value, has negligible effect. Only the drag on the puff has a significant effect in reducing the distance traveled by the puff. It can then be taken that the puff evolution to good accuracy can be represented by (16) . Over a time span of 10 nondimensional units the puff has traveled about 0.7s e and the velocity has dropped to about 15% of the initial velocity. By 100 nondimensional units the puff has traveled about 1.75s e and the velocity has dropped to about 2.5% of the initial velocity.
The ejected droplets are made of a complex fluid that is essentially a mixture of oral fluids, including secretions from both the major and minor salivary glands. In addition, it is added up with several constituents of non-salivary origin, such as gingival crevicular fluid, exhalted bronchial and nasal secretions, serum and blood derivatives from oral wounds, bacteria and bacterial products, viruses and fungi, desquamated epithelial cells, other cellular components, and food debris [54] . Therefore, it is not easy to determine precisely transport properties of the droplet fluid. Although surface tension is measured similar to that of water, viscosity can be one or two orders of magnitude larger [37] making drops less coalescence prone [92, 93] . In the present context, viscosity and surface tension might be of importance, because they can influence droplet size distribution specifically by controlling coalescence and breakage. These processes are important only during the ejection stage, and once droplets are in the range below 50 µm, coalescence and break up processes are impeded, allowing us to model drops as single, non-interacting small spheres. The ejected swarm of droplets is characterized by its initial size spectrum as given in (4). The time evolution of the spectrum of droplets that remain within the puff in terms of droplet size, velocity and temperature is the object of interest in this section. This evolution of the ejected droplets depend on the following four important parameters: the time scale τ V on which the droplet velocity relaxes to the puff fluid velocity (in the absence of other forcing), the time scale τ T on which the droplet temperature relaxes to the puff fluid temperature, the settling velocity W of the droplet within the puff fluid, and the Reynolds number Re based on settling velocity. These quantities are given by [9, 63, 64] 
where ρ ≈ 1000 is the droplet-to-air density ratio, C r ≈ 4.16 is the droplet-to-air specific heat ratio, g is acceleration due to gravity, ν p and κ p are the kinematic viscosity and thermal diffusivity of the puff. In the above, Φ = 1 + 0.15Re 0.687 and N u = 2 + 0.6Re 1/2 P r 1/3 are the finite Reynolds number drag and heat transfer correction factors; both of which simplify in the Stokes regime for drops smaller than about 50 µm. Here we take the Prandtl number of air to be P r = 0.72. In the Stokes limit, the velocity and thermal time scales, and the settling velocity of the droplet increase as d 2 , while Reynolds number scales as d 3 . The value of these four parameters for varying droplet sizes is presented in Fig. 8 , where it is clear that the effect of finite Re becomes important only for droplets larger than 50 µm. For smaller droplets τ V , τ T 1(s), W 1(m/s), and Re 1. 
The size of the droplets under investigation is sufficiently small, and the swarm is dilute to prevent their coalescence, and the only way in which droplets change their size is via evaporation. According to the analysis of Langmuir [58] , the rate of mass loss due to evaporation of a small sphere depends on the diffusion of the vapor layer away from the sphere surface, and under reasonable hypotheses [18, 58, 86, 95] , it can be expressed as :
where, m is the mass of a droplet of diameter d, D is the diffusion coefficient of the vapor, ρ p is the density of puff air and B m = (Y d − Y p )/(1 − Y s ) is the Spalding mass number, where Y d is mass fraction of water vapor at the droplet surface and Y p is mass fraction of water vapor in the surrounding puff. Under the assumption that N u and B m are nearly constant for small droplets, the above equation can be integrated to obtain the following law (mapping) for the evolution of the droplet:
where d e is the initial droplet diameter at ejection and k = 4DN u ln(1 + B m )/ρ has units of m 2 /s and thus represent an effective evaporative diffusivity. It is important to observe that (20) would predict a loss of mass per unit area tending to infinity as the diameter of the drop tends to zero. This implies that the droplet diameter goes to zero in a finite time and we establish the result
which for any time t yields a critical value of droplet diameter, and all droplets that were smaller, or equal, at exhalation (i.e., d e ≤ d e,evap ) would have fully evaporated by t. The only parameter is k . Assuming N u = 2 and D = 2.8 × 10 −5 m 2 /s, even for very small values of B m , we obtain the evaporation time for a 10 µm droplet to be less than a second. However, it appears that smaller than a certain critical size, the loss of mass due to evaporation slows down [18] . This could partly be due to the presence of non-volatiles and other particulate matter within the droplet, whose effects were ignored in the above analysis, and will be addressed in section 4.3. It seems that (20) can give reliable predictions for droplet diameter down to few µm with much slower evaporation rates for smaller size. Irrespective of whether water completely evaporates leaving only the non-volatile droplet residue, or the droplet evaporation slows down, the important consequence on the evolution of the droplet size distribution is that it is narrower and potentially centered around micron size.
We now consider the motion of the ejected droplets, while they rapidly evaporate. The equation of motion of the droplet is the Newton's law
where e e e z is the unit vector along the vertical direction, m p is the mass of puff displaced by the droplet, V V V d and v v v p are the vector velocity of the droplet and the surrounding puff. Provided the droplet time scale τ V is smaller than the time scale of surrounding flow, which is the case for droplets of diameter smaller than 50 µm, the above ODE can be perturbatively solved to obtain the following leading order solution [8, 32, 33] 
According to the above equation, the equilibrium Eulerian velocity of the droplet depends on the local fluid velocity plus the still fluid settling velocity W of the droplet plus the third term that arises due to the inertia of the droplet. Though at ejection the droplet speed is smaller than the surrounding gas velocity, as argued in section 3.2, the droplets quickly accelerate to approach the puff velocity. In fact, since the puff is decelerating (i.e. |dv v v p /dt| < 0), the droplet velocity will soon be larger than the local fluid velocity. As long as the droplet stays within the puff, the velocity and acceleration of the surrounding fluid can be approximated by those of the puff as |v v v p | = ds/dt and |dv v v p /dt| = d 2 s/dt 2 . This allows evaluation of the relative importance of the third term (inertial slip velocity) in terms of the puff motion, which is given in (16) as [63] 
This ratio takes its largest value at the initial time of injection and then decays as 1/t. Using the range of possible values of t e given earlier, this ratio is small for a wide range of initial droplet sizes. We thus confirm that for the most part droplet inertia can be ignored in its motion, and the droplet velocity can be taken to be simply the sum of local fluid velocity and the still fluid settling velocity of the droplet.
While the effect of buoyancy on the puff was shown to be small, the same cannot be said of the droplets. The vertical motion of a droplet with respect to the surrounding puff, due to its higher density, is dependent only on the fall velocity W , which scales as d 2 , which in turn decreases as given in (21) due to evaporation. The droplet's gravitational settling velocity can be integrated over time to obtain the distance over which it falls as a function of time. We now set this fall distance (left hand side) equal to the puff radius (right hand side) to obtain
where we have set the droplet diameter at exhalation to be d e,exit , indicating the fact that a droplet of initial diameter equal to d e,exit has fallen by a distance equal to the puff size at time t. Thus all larger droplets of size d e > d e,exit have fallen out of the puff by t and we have been referring to these as the exited droplets. The two critical initial droplet diameters, d e,evap and d e,exit are plotted in Fig. 9a as a function of t. The only other key parameter of importance is k whose value is varied from 10 −12 to 10 −6 m 2 /s. In evaluating d e,exit using (26), apart from the property values of water and air, we have used the nominal values of α = 0.1, s e = 0.5 m and t e = 0.05 s (as an example). The solid lines correspond to d e,exit , which decreases with increasing t and for each value of k there exits a minimum d e below which there is no solution to (26) since the droplet fully evaporates before falling out of the puff. The dotted lines correspond to d e,evap , which increases with t. The intersection of the two curves is marked by the solid square, which corresponds to the limiting time t lim (k ), beyond which the puff contains only fully-evaporated droplet residues containing the viruses. Correspondingly we can define a limiting droplet diameter d e,lim (k ). Given sufficient time, all initially ejected larger droplets (i.e., d e > d e,lim ) would have fallen out of the puff and all smaller droplets (i.e., d e ≤ d e,lim ) would have evaporated to become droplet residues.
At times smaller than the limiting time (i.e., for t < t lim ) we have the interesting situation of some droplets falling out of the puff (exited droplets), some still remaining as partially evaporated airborne droplets, and some fully-evaporated to become droplet residues. This scenario is depicted in Fig. 9a with an example of t = 3 s for k = 10 −10 m 2 /s plotted as a dash line.
There can be significant presence of non-volatile material such as mucus, bacteria and bacterial products, viruses and fungi, and food debris in the ejected droplets [54] . However, the fraction of ejected droplet volume Q de that is made up of these non-volatiles varies substantially from person to person. The presence of non-volatiles alters the analysis of the previous sections in two significant ways. First, each ejected droplet, as it evaporates, will reach a final size that is dictated by the amount of non-volatiles that were initially in it. The larger the droplet size at initial ejection, the larger will be its final size after evaporation, since it contains a larger amount of non-volatiles. If ψ is the volume fraction of non-volatiles in the initial droplet, the final diameter of the droplet residue after complete evaporation of volatile matter (i.e., water) will be d dr = d e ψ 1/3 . This size depends on the initial droplet size and composition. Note that even a small, for example 1%, non-volatile composition results in d dr being around 20% of the initial ejected droplet size. It has also been noted that the evaporation of water can be partial, depending on local conditions in the cloud or environment. We simply assume the fraction ψ to also account for any residual water retained within the droplet residue.
The second important effect of non-volatile is to reduce the rate of evaporation. As evaporation occurs at the droplet surface, a fraction of the surface will be occupied by the non-volatiles reducing the rate of evaporation. For small values of ψ, the effect of non-volatiles is quite small only at the beginning. The effect of non-volatiles will increase over time, since the volume fraction of nonvolatiles increases as the volatile matter evaporates. Because of this ever decreasing evaporation rate, it may take longer for a droplet to decrease from its ejection diameter of d e to its final droplet residue diameter of d dr , than what is predicted by (21) . It should be noted that intermittency of turbulence and heterogeneity of vapor concentration and droplet distribution within the puff will influence the evaporation rate [27, 30, 74] . Nevertheless, for simplicity, and for the purposes of the present first order mathematical framework, we use the d 2 -law given in (21) , but with a smaller value of effective k to account for the effect of non-volatiles and turbulence intermittency. This approximation is likely to be quite accurate in describing the early evolution of the droplet. Only at late stages as the droplet approaches its final diameter d dr , the d 2 -law will be in significant error.
Applying the analysis of the previous sections, taking into account the presence of non-volatiles, we separate the two different time regimes of t ≤ t lim and t ≥ t lim . In the case when t ≤ t lim , we have three types of droplets: (i) exited droplets whose initial size at injection is greater than d e,exit , (ii) droplets of size at ejection smaller than d e,evap that have completely evaporated to become droplet residues of size d dr and (iii) intermediate size airborne droplets that are within the puff and still undergoing evaporation. We assume an equation of the form (26) to approximately apply even in the presence of non-volatiles. With this balance between fall distance of a droplet and the puff radius we obtain the following expression d e,exit = 18ν a ρg αs e t t + t e t e 1/(4+C)
The corresponding limiting diameter of complete evaporation can be obtained from setting d = d e,evap ψ 1/3 and d e = d e,evap in (21) as
While the above two estimates are in terms of the droplet diameter at injection, their current diameter at t can be expressed as
Form the above expressions, we define t lim to be the time when d e,exit = d e,evap , which in terms of current droplet diameter becomes d exit = d evap . Beyond this limiting time (i.e., for t > t lim ) the droplets can be separated into only two types: (i) exited droplets whose initial size at injection greater than d e,exit = d e,evap , and (ii) droplets of size at ejection smaller that have become droplet residues. The variation of t lim and d e,lim as a function of k is presented in Fig. 9b . It is clear that as k varies over a wide range, t lim ranges from 0.01 s to 450 s, and correspondingly d e,lim varies from 415 to 7 µm.
We now put together all the above arguments to present a predictive model of the droplet concentration within the puff. The initial condition for the size distribution is set by the ejection process discussed in Section 3, and the simple Pareto distribution given in (4) provides an accurate description. Based on the analysis of the previous sections, we separate the two different time regimes of t ≤ t lim and t ≥ t lim .
In the case when t ≤ t lim the droplet/aerosol concentration (or the number per unit volume of the puff) can be expressed as
where we have recognized the fact that equation (21) is the mapping D between the current droplet size and its size at injection. Due to the turbulent nature of the puff, the distribution of airborne droplets and residues is taken to be uniform within the puff. Quantities such ass, d evap and d exit are as they have been defined above and the pre-factor 1/Q(t) accounts for the expansion of the puff volume. In the case of t ≥ t lim , the droplet number density spectrum becomes
and only droplet residues remain within the puff. Here, the size of the largest droplet residue within the puff is related to its initial unevaporated droplet size as d lim = d e,lim ψ 1/3 , and the plot of d e,lim as a function of k for a specific example case of puff and droplet ejection was shown in Fig. 9b .
In this subsection we will briefly consider droplet temperature, since it plays a role in determining saturation vapor pressure and the value of k . Following Pirhadi et al. [86] we write the thermal equation of the droplet as
where C pw is specific heat of water, k p is thermal conductivity of the puff air, L is the latent heat of vaporization, T d and T p are the temperatures of the droplet and the surrounding puff. The first term on the right accounts for convective heat transfer from the surrounding air and the second term accounts for heat needed for phase change during evaporation. It can be readily established that the major portion of heat required for droplet evaporation must come from the surrounding air through convective heat transfer. The equilibrium Eulerian approach [32] can again be used to obtain the asymptotic solution of the above thermal equation and the droplet temperature can be explicitly written as
where τ T is the thermal time scale of the droplet that was introduced earlier. The second term on the right is negative and thus contributes to the droplet temperature being lower than the surrounding puff. Simple calculation with typical values shows that the contribution of the third term is quite small and can be ignored. As a result, the temperature difference between the droplet and the surrounding is largely controlled by the evaporation rate dm/dt, which decreases over time.
Again, using the properties of water and air, and typical values for N u and B m , we can evaluate the temperature difference T p − T d to be typically a few degrees. Thus, the evaporating droplets need to be only a few degrees cooler than the surrounding puff for evaporation to continue.
When the puff equilibrates with the surrounding and its velocity falls below the ambient turbulent velocity fluctuation, the subsequent dynamics of the droplet cloud is governed by turbulent dispersion. This late-time evolution of the droplet cloud depends on many factors that characterize the surrounding air. This is where the difference between a small enclosed environment such as an elevator or an aircraft cabin or an open field matters, along with factors such as cross breeze and ventilation. A universal analysis of late-time evolution of the droplet residue cloud is thus not possible, due to problem-specific details. The purpose of this brief discussion is to establish a simple scaling relation to guide when the puff evolution model presented in the above sections gives way to advection and dispersion by ambient turbulence. It should again be emphasized that the temperature difference between the puff fluid containing the droplet residue cloud and the ambient air may induce buoyancy effects, which for model simplicity will be taken into account as part of turbulent dispersion.
We adopt the classical scaling analysis of Richardson [91] , according to which the radius of a droplet cloud, in the inertial range, will increase as the 3/2 power of time as given by
where c is a constant, is the dissipation rate, which will be taken to be a constant property of ambient turbulence, and t 0 is the time shift required to match the cloud size at the transition time between the above simple late time model and the puff model. In the above, the subscript lt stands for late-time behavior of the radius of the droplet-laden cloud. We now make a simple proposal that there exists a transition time t tr , below which the rate of expansion of the puff as given by the puff model is larger than dr lt /dt computed from the above expression. During this early time, ambient dispersion effects can be ignored in favor of the puff model. But for t > t tr droplet-laden cloud's ambient dispersion becomes the dominant effect. The constants t 0 and t tr can be obtained by satisfying the two conditions: (i) the size of the droplet-laden cloud given by (35) at t tr matches the puff radius at that time given by αs e ((t tr + t e )/t e ) 1/(4+C) , and (ii) the rate of expansion of the droplet-laden cloud by turbulent dispersion matches the rate of puff growth given by the puff model. This latter condition can be expressed as
From these two simple conditions we obtain the final expression for the transition time as
Given a puff, characterized by its initial ejection length and time scales s e and t e , and the ambient level of turbulence characterized by , the value of transition time can be estimated. If we take entrainment coefficient α = 0.1, the constant C = 0, and typical values of s e = 0.5 m and t e = 0.05 s, we can estimate t tr = 1.88 s for a dissipation rate of c = 10 −5 m 2 /s 3 . The transition time t tr increases (or decreases) slowly with decreasing (or increasing) dissipation rate. Thus, the early phase of droplet evaporation described by the puff model is valid for O(1) s, before being taken over by ambient turbulent dispersion. However, it must be stressed that the scaling relation of Richardson is likely an over-estimation of ambient dispersion, as there are experimental and computational evidences that suggest that the power-law exponent in (35) is lower than 3 [82] . But it must be remarked that even with corresponding changes to late-time turbulent dispersion, the impact on transition time can be estimated to be not very large. Also, it must be cautioned that according to classical turbulent dispersion theory, during this late-time dispersal, the concentration of virus-laden droplet residues within the cloud will not be uniform, but will tend to decay from the central region to the periphery. Nevertheless, for sake of simplicity here we assume (35) to apply and we take the droplet residue distribution to be uniform.
According to above simple hypothesis, the effect of late-time turbulent dispersion on number density spectrum is primarily due to the expansion of the could, while the total number of droplet residues within the cloud remains the same. Thus, the expressions (31) and (32) still apply. However the expression for the volume of the cloud must be appropriately modified as
The location of the center of the expanding cloud of droplets can be still taken to be given by the puff trajectory s(t), which has considerably slowed down during late-time dispersal. The strength of the above model is in its theoretical foundation and analytical simplicity. But, the validity of the approximations and simplifications must be verified in applications to specific scenario being considered. For example, considering variability in composition, turbulence intermittency, initial conditions of emissions and the state of the ambient, direct observations show that the transition between puff dominated and ambient flow dominated fate of respiratory droplets vary from O(1-100 s) [17] .
This section will mainly survey the existing literature on issues pertaining to what fraction of the droplets and aerosols at any location gets inhaled by the recipient host, and how this is modified by the use of masks. These effects modeled as inhalation and filtration efficiencies will then be incorporated into the puff-cloud model. The pulmonary ventilation (breathing) has a cyclic variation that varies markedly with age and metabolic activities. The intensity of breathing (minute ventilation) is expressed in L/min of inhaled and exhaled air. For the rest condition, the ventilation rate is about 5-8 L/min and increases to about 10-15 L/min for mild activities. During exercise, ventilation increases significantly depending on age and metabolic needs of the activity.
In the majority of earlier studies on airflow and particle transport and deposition in human airways, the transient nature of breathing was ignored for simplification and to reduce the computational cost. Haubermann et al. [45] performed experiments on a nasal cast and found that particle deposition for constant airflow is higher than those for cyclic breathing. Shi et al. [98] performed simulations on nanoparticle depositions in the nasal cavity under cyclic airflow and found that the effects of transient flow are important. Grgic et al. [40] and Horschler et al. [48] performed experimental and numerical studies, respectively, on flow and particle deposition in a human mouth-throat model, and the human nasal cavity. Particle deposition in a nasal cavity under cyclic breathing condition was investigated by Bahmanzadeh et al. [7] , Naseri et al. [78] , and Kiasadegh et al. [56] , where the unsteady Lagrangian particle tracking was used. They found there are differences in the predicted local deposition for unsteady and equivalent steady flow simulations.
In many of these studies, a sinusoidal variation for the volume of air inhaled is used. That is
Here Q max is the maximum flow rate, and T = 4 s is the period of breathing cycle for an adult during rest or mild activity. The period of breathing also changes with age and the level of activity.
Haghnegahdar et al. [42] investigated the transport, deposition, and the immune system response of the low-strain Influenza A Virus IAV laden droplets. They noted that the shape of the cyclic breathing is subject dependent and also changes with nose and mouth breathing. They provided an eight-term Fourier series for a more accurate descriptions of the breathing cycle. The hygroscopic growth of droplets was also included in their study. Analysis of aspiration of particles through the human nose was studied by Ogden and Birkett [81] and Armbruster and Breuer [2] . Accordingly, the aspiration efficiency η a is defined as the ratio of the concentration of inhaled particles to the ambient concentration. Using the results of earlier studies and also his works, Vincent [109] proposed a correlation for evaluating the inhalability of particles. That is, the aspiration efficiency η a of particles smaller than 100 µm is given as,
Here, d is an aerodynamic diameter of particles. While the above correlation provides the general trend that larger particles are more difficult to inhale, it has a number of limitations. It was developed for mouth-breathing with the head oriented towards the airflow direction with speeds in the range of 1 m/s to 4 m/s. The experimental investigation of aerosol inhalability was reported by Hsu and Swift [49] , Su and Vincent [101, 102] , Aitken et al. [1] , and Kennedy and Hinds [55] .
Dai et al. [22] performed in-vivo measurements of inhalability of large aerosol particles in calm air and fitted their data to several correlations. For calm air condition, they suggested, 
Computational modeling of inhalability of aerosol particles were reported by many researchers [19, 51, 52, 57, 75, 78] . Interpersonal exposure was studied by [36, 46] . The influence of thermal plume was studied by Salmanzadeh et al. [94] . Nasiri et al. [77] performed a series of computational modeling and analyzed the influence of the thermal plume on particle aspiration efficiency when the body temperature is higher or lower than the ambient. Their results are reproduced in Figure 10 .
Here the case that the body temperature T b = 26.6 • C and the ambient temperature T a = 21.3 • C (Upward Thermal Plume) and the case that T b = 32.2 • C and T a = 40.0 • C (Downward Thermal Plume) are compared with the isothermal case studied by Dai et al. [22] . It is seen that when the body is warmer than the surrounding, the aspiration ratio increases. When the ambient air is at a higher temperature than the body, the inhalability decreases compared to the isothermal case. In light of the results of the previous section it can be concluded that at a distance of O(1) m the ejected droplets have sufficiently reduced in size that these O(1) µm aerosols have near perfect inhalability. However, recall that this estimation must be adjusted to account for the possible slow down of evaporation due to variability in the initial conditions of emissions and the ambient, intermittency of cloud turbulence and drop concentration. Thus, the above conclusion represents a lower bound of timescale and distance of transition to fully evaporated droplet residues from the cloud model [17] . 
Using a respiratory face mask is a practical approach against exposure to airborne viruses and other pollutants. Among the available facepiece respirators, N95, and surgical masks are considered to be highly effective [41, 68] . N95 mask has a filtration efficiency of more than 95% in the absence of face leakage [79, 89] . Surgical masks are used extensively in the hospital and operating rooms [67] . Nevertheless, there have been concerns regarding their effective filtration of airborne bacteria and viruses [10, 59, 68] . There is often discomfort in wearing respiratory masks for extended durations that increases the risk of spread of infection. The breathing resistance of a mask is directly related to the pressure drop of the filtering material. The efficiency of respiratory masks varies with several factors, including the intensity and frequency of breathing as well as the particle size [119] . The filtration efficiencies of different masks under normal breathing conditions, as reported by Zhang et al. [118] and Feng et al. [31] , are shown in Figure 11 . It is seen that the filtration efficiencies of different masks vary significantly, with N95 having the best performance, which is followed by the surgical mask. It is also seen that all masks could capture large particles. The N95, Surgical, and Procedure masks remove aerosols larger than a couple of microns. Cotton and Gauze masks capture a major fraction of particles larger than 10 µm. The capture efficiency of all masks also shows an increasing trend as particle size becomes smaller than 30 nm due to the effect of the Brownian motion of nanoparticles. Figure 11 also shows that the filtration efficiencies of all respiratory masks drop for the particle sizes in the range of 80 nm to about 1 µm. This is because, in this size range, both the inertia impaction and the Brownian diffusion effect are small, and the mask capture efficiency reduces. Based on these results, and the earlier finding that most ejected droplets within the cloud have become sub-micron-sized aerosol particles by about O(1 − 10) m distance, it can be stated that only professional masks such as N95, Surgical, and Procedure masks provide reliable reduction in the inhaled particles. Hence, the importance for healthcare workers to have access to high-grade respirators upon entering a room or space with infectious patients [17] .
It should be emphasize that the concentration that a receiving host will inhale (φ inhaled ) depends on the local concentration in the breathing zone adjusted by the aspiration efficiency given by Equations (40) and (41) (or plotted in Figure 9 ). When the receiving host wears a mask, an additional important correction is needed by multiplying by a factor (1 − η f ), where η f is the filtration efficiency plotted in Figure 9 . That is,
where φ(d, t) is the droplet residue concentration at the breathing zone given in (31) or (32) . It is seen that the concentration of inhaled droplet larger than 10 microns significantly decreases when the mask is used. But the exposure to smaller droplets, particularly, in the size range of 100 nm to 1 µm varies with the kind of mask used.
The object of this section is to put together the different models of the puff and droplet evolution described in the previous sections, underline their simplifications, and demonstrate their ability to make useful prediction. Such results under varying scenarios can then be potentially used for science-based policy making, such as establishing multi-layered social distancing guidelines and other safety measures. In particular, we aim at modeling the evolution of the puff and the concentration of airborne droplets and residues that remain within the cloud so that the probability of potential transmission can be estimated.
As discussed in section 4.2, the virus-laden droplets exhaled by an infected host will undergo a number of transformations before reaching the next potential host. To prevent transmission, current safety measures impose a safety distance of two meters. Furthermore, cloth masks are widely used by public and their effectiveness has been shown to be questionable for droplets and aerosols of size about a micron. The adequacy of these common recommendations and practices can be evaluated by investigating the concentration of airborne droplets and residues at distances larger than one meter and the probability of them being around a micron in diameter, since such an outcome will substantially increase the chances of transmission. In the following we will examine two effects: the presence of small quantities of non-volatile matter in the ejected drops that remain as droplet residues after evaporation, and the adequacy of the log-normal or Pareto distribution to quantify the number of droplets in the lower diameter classes.
Let us consider the situation of speaking or coughing, whose initial puff volume and momentum are such that they yield s e 0.5 m and t e 0.05 s. Under this specific condition, as shown in Figure 7 the puff travels about 1 m in about 5 s 1 . For this simple example scenario, we will examine our ability to predict airborne droplet and residue concentration, as an important step towards estimating potential for airborne transmission in situations commonly encountered.
In most of the countries, current guidelines are based on the work by Xie et al. [117] , who revisited previous guidelines by [114] with improved evaporation and settling models. They identified the possibility that, due to evaporation, the droplets quickly become aerosolized before reaching a significant distance and thus may represent a minor danger for transmission due to their minimal virus loading. This scenario is shown in Figure 12 , where we present the evolution of the drop size Figure 12 : Evolution of the drop size distribution spectra according to the currently used evaporation models [114, 117] . spectrum while droplets are transported by the ejected puff. The initial droplet size distribution is taken to be that measured by Duguid [26] modeled with a log-normal distribution, which in the Monte-Carlo approach is randomly sampled with one million droplets divided into one thousand diameter classes. Each droplet is then followed while evaporating and falling. The evaporation model is taken to be (21) with the effective diffusion coefficient estimated as k 1 · 10 −8 m 2 /s. This value is computed under the assumption that drops are made of either pure water or a saline solution [117] and that air has about 98% humidity. Therefore, this is an environment unfavorable to evaporation and consequently drop size reduction happens relatively slowly. However, from the figure it is clear that, even in this extreme case, after few tens of centimeters, and within a second, all droplets have evaporated down to a size below 10 µm. This is in line with the predictions of Xie et al. [117] . Naturally, if the air is dryer, the effective evaporation coefficient will be larger (even as large as k 10 −5 m 2 /s) and the droplet size spectrum will evolve even faster, leaving virtually all droplets to be smaller than 1 µm in the puff. Recall that intermittency of turbulence with the puff can create clusters of droplets and concentration of vapor and thereby significantly alter the evaporation rate [27, 30, 74] . Hence, our estimate of evaporation time is a lower bound, as governed by the d 2 -law (21).
As discussed in section 4.3 there is current consensus that droplets ejected during sneezing or coughing contain, in addition to water, other biological and particulate non-volatile matter. Specifically, viruses themselves are just large protein chains of size almost 0.1 µm. Here we will examine the evolution of droplet size distribution in the presence of non-volatile matter. It will be clear in the following, that in this case, even a small amount of non-volatile matter plays an important role with the evaporation coefficient being a minor factor in deciding how fast the final state is reached. In Figure 13 , we show the final distribution of droplets under two scenarios where the initially ejected droplets contain 0.1% and 3.0% of non-volatile matter. In Figure 13a , the initial drop size distribution is modeled as a log-normal distribution (i.e. as in Fig. 12 ) whereas in Figure 13b , the initial drop size distribution is modeled according to the Pareto distribution with initial droplet size varying between 1 and 100 µm. This range is smaller than that suggested earlier in section 3. However drops that are larger than 100 µm fall out of the cloud and therefore are not important for airborne transmission and droplets initially smaller than 1 µm have much smaller viral load. Here "final droplet size distribution" indicates the number of droplets that remain within the puff after all the larger droplets have fallen out and all others have completed their evaporation to become droplet residues. This final number of droplet residues as a function of size does not vary with time or distance.
The size distribution is computed here as in Figure 12 , with a random sampling from the initial log-normal or Pareto distribution. As before, these computations used an evaporation coefficient of k = 10 −8 m 2 /s. However, there are two important differences: each droplet is allowed to fall vertically according to its time-dependent settling velocity, W , which decreases over time as the droplet evaporates. Integration of the fall velocity over time provides the distance traveled by the droplet relative to the puff. Droplets whose fall distance exceeds the size of the puff are removed from consideration. Second, each droplet that remains within the puff evaporates to its limiting droplet residue size that is dictated by the initial amount of non-volatile matter contained within the droplet. For ψ = 0.1% non-volatile matter, the final aerosol size cannot decrease below 10% of the initial droplet diameter, whereas for 3.0% of non-volatile matter, the final droplet size cannot decrease below 30% of the initial diameter. From Fig. 13 , it is clear that when eveporation is complete the drop size distribution rigidly shifts towards smaller diameters, with a cut-off upper diameter due to the settling of large drops. Essentially, it is clear that the initial number of viruses that were in droplets of size smaller than d e,exit still remain within the cloud almost unchanged, representing a more dangerous source of transmission than predicted by the conventional assumption of near-full evaporation. Again, it is important to note that the final droplet size distribution is established rapidly even with the somewhat lower effective evaporation diffusivity of k = 10 −8 m 2 /s, and when not accounting for the effect of localized moisture of the cloud in further reducing the rate. Figure  13 also illustrates the important difference in the drop size distribution. The Pareto distribution will predict a much larger number of drops in the micron and sub-micron range, possibly the most dangerous for both inhalation efficiency and filtration inefficiency.
In this section we will demonstrate the efficacy of the simple model presented in (31) and (32) for the prediction of droplet/aerosol concentration. In contrast to the Monte-Carlo approach of the previous subsection, where the evolution of each droplet was accurately integrated, here we will use the analytical prediction along with its simplifying assumptions. The cases considered are identical to those presented in Figure 13 for ψ = 0.1% and k = 10 −8 m 2 /s. The initial droplet size distributions considered are again log-normal and Pareto distributions. In this case, however, we underline that the quantity of importance in airborne transmission is not the total number of droplet residues, but rather their concentration in the proximity of a susceptible host. Accordingly, we plot in Figure 14 airborne droplet and residue concentration (per liter) of volume as a function of droplet size. Here the area under the curve between any two diameter yields the number of droplets within this size range per liter of volume within the cloud. At the early times of t = 0.025 and 0.2 s, we see that larger droplets above a certain size have fallen out of the cloud, while droplet residues smaller than d evap have fully evaporated and their distribution is a rigidly-shifted version of the original distribution. The distribution of intermediate size airborne droplets reflects the fact that they are still undergoing evaporation. Unlike in Figure 13 , the concentration continues to fall even after t lim 0.68 s when the number and size of droplets within the cloud have reached their limiting value. This is simply due to the fact that the volume of the puff continues to increase and this continuously dilutes the aerosol concentration. Most importantly, the results of the simple model presented in (31) and (32) are in excellent agreement with those obtained from Monte-Carlo simulation. The increasing size of the contaminated cloud with time can be predicted with (38) and the centroid is given by the scaling law (16) . 
The primary goal of this paper is to provide a unified theoretical framework that accounts for all the physical processes of importance, from the ejection of droplets by breathing, talking, coughing and sneezing to the inhalation of resulting aerosols by the receiving host. These processes include: (i) forward advection of the exhaled droplets with the puff of air initially ejected; (ii) growth of the puff by entrainment of ambient air and its deceleration due to drag; (iii) gravitational settling of some of the droplets out of the puff; (iv) modeling of droplets evaporation, assuming that the d 2 -law prevails; (v) presence of non-volatile compounds which form the droplet residue left behind after evaporation; (vi) late-time dispersal of the droplet residue-laden cloud due to ambient air turbulent dispersion.
Despite the complex nature of the physical processes involved, the theoretical framework results in a simple model for the airborne droplet and residue concentration within the cloud as a function of droplet diameter and time, which is summarized in equations (31), (32) and (38) . This framework can be used to calculate the concentration of virus-laden residues at the location of any receiving host as a function of time. As additional processes, the paper also considers (vii) efficiency of aspiration of the droplet residues by the receiving host; and (viii) effectiveness of different kinds of masks in filtering the residues of varying size.
It must be emphasized that the theoretical framework has been designed to be simple and therefore involves a number of simplifying assumptions. Hence, it must be considered as the starting point. By relaxing the approximations and by adding additional physical processes of relevance, more complex theoretical models can be developed. One of the primary advantages of such a simple theoretical framework is that varying scenarios can be considered quite easily: these different scenarios include varying initial puff volume, puff velocity, number of droplets ejected, their size distribution, non-volatile content, ambient temperature, humidity, and ambient turbulence.
The present theoretical framework can be, and perhaps must be, improved in several significant ways in order for it to become an important tool for reliable prediction of transmission. (i) Accurate quantification of the initially ejected droplets still remains a major challenge. Further high-quality experimental measurements and high-fidelity simulations are required, especially mimicking the actual processes of talking, coughing and sneezing, to fully understand the entire range of droplet sizes produced during the exhalation process. (ii) As demonstrated above, the rate at which an ejected droplet evaporates plays an important role in determining how fast they reach their fully-evaporated state. It is thus important to calculate more precisely the evaporation rate of non-volatile-containing realistic droplets resulting from human exhalation. The precise value of evaporation rate may not be important when droplets evaporate fast, since all droplets remaining within the puff would have completed their evaporation. But under slow evaporation conditions, accurate evaluation of evaporation is important. (iii) The assumption of uniform spatial distribution of droplets within the puff and later within the dispersing cloud is a serious approximation [96] . The intermittency of turbulence within the initial puff and later within the droplet cloud is important to understand and couple with the evaporation dynamics of the droplets. In addition to the role of intermittency, even the mean concentration of airborne droplets and residues may decay from the center to the outer periphery of the puff/cloud. Characterization of this inhomogeneous distribution will improve the predictive capability of the model. (iv) The presence of significant ambient mean flow and turbulence either from indoor ventilation or outdoor cross-flow will greatly influence the dispersion of the virus-laden droplets. But accounting for their effects can be challenging even in experimental and computational approaches. Detailed experiments and highly-resolved simulations of specific scenarios should be pursued. But it will not be possible to cover all possible scenarios with such an approach. A simpler approach where the above theoretical framework can be extended to include additional models such as random flight model (similar to those pursued in the calculation of atmospheric dispersion of pollutants [43] ) may be promising approaches.