The spreading of viruses by airborne aerosols: lessons from a first-passage-time problem for tracers in turbulent flows

We study the spreading of viruses, such as SARS-CoV-2, by airborne aerosols, via a new firstpassage-time problem for Lagrangian tracers that are advected by a turbulent flow: By direct numerical simulations of the three-dimensional (3D) incompressible, Navier-Stokes equation, we obtain the time tR at which a tracer, initially at the origin of a sphere of radius R, crosses the surface of the sphere for the first time. We obtain the probability distribution function P(R, tR) and show that it displays two qualitatively different behaviors: (a) for R LI, P(R, tR) has a power-law tail ∼ t−α R , with the exponent α = 4 and LI the integral scale of the turbulent flow; (b) for LI . R, the tail of P(R, tR) decays exponentially. We develop models that allow us to obtain these asymptotic behaviors analytically. We show how to use P(R, tR) to develop social-distancing guidelines for the mitigation of the spreading of airborne aerosols with viruses such as SARS-CoV-2.


INTRODUCTION
By 1 June 2020 (14:31 GMT) the COVID-19 coronavirus pandemic had affected 213 countries and territories and 2 international conveyances; the numbers of cases and deaths were, respectively, 6, 300, 444 and 374, 527 [1]. Social distancing has played an important role in mitigation strategies that have been used in several countries to arrest the spread of COVID- 19 [2]. To optimise social-distancing guidelines we must ask: How far, and how fast, do small respiratory droplets or virusbearing aerosols spread in turbulent flows? Given the ongoing COVID-19 pandemic, it is important and extremely urgent to have at least a semi-quantitative answer to this question. SARS-CoV-2, the virus that causes COVID-19, spreads, principally, in two different ways: (1) First: Respiratory droplets, ejected by the sneeze or cough of a patient, fall on nearby surfaces or persons; in this case, approximate estimates of the distance, over which droplets are likely to travel, are available [3][4][5].
(2) Second: Transmission of this virus can occur because of airborne aerosols, such as, (a) a cloud of fine droplets, with diameters smaller than 5 micrometers, emitted by an infected person while speaking loudly [6] or (b) the SARS-CoV-2 RNA on fine, suspended particulate matter [7]. These aerosols may remain suspended in the air for a long time. Indeed, they have been reported in two hospitals in Wuhan [8]; and there is growing evidence that the SARS-CoV-2 virus could also spread via airborne aerosols [6-10], typically indoors [11]. Other diseases can also spread because of airborne aerosols; examples include measles [12], chickenpox [13], tuberculosis [14], and avian flu [15].
The typical sedimentation speed for such aerosols is comparable to their thermal speed. Therefore, at the simplest level, it is natural to model these aerosol parti-cles as neutrally-buoyant Lagrangian tracers, which are advected by the flow, but are passive, in the sense that they do not affect the flow velocity. We can then study the spread of viruses, such as SARS-CoV-2, via the airborne-aerosol route, by considering the advection of such tracers by turbulent fluid flows. There have been extensive studies of such tracers in the fluid-dynamics literature [16,17]; and models for such tracers have been used, inter alia, to model the dispersion of pheromones by lepidoptera [18].
We would like to determine the time that an aerosol particle (one of the red particles in the schematic diagram of Fig. 1) takes to travel a distance R from its source (the man at the centre of Fig. 1). In a turbulent flow, this time is random; furthermore, a tracer particle can go past the distance R, turn back, and reach R again. It is important, therefore, to calculate the time it takes for an aerosol particle to reach the distance R for the first time and to calculate the probability distribution function (PDF) of the first-passage time of a tracer in a turbulent flow. We carry out this calculation below.
Specifically, we consider Lagrangian tracer particles that emanate from a point source in a turbulent fluid. If t R is the time at which a tracer, initially at the origin of a sphere of radius R, crosses the surface of the sphere for the first time, what is the probability distribution function (PDF) P(R, t R )? The answer to this question is of central importance in both fundamental nonequilibrium statistical mechanics [19][20][21][22][23] and in understanding the dispersal of tracers by a turbulent flow, a problem whose significance cannot be overemphasized, for it is of relevance to the advection of (a) airborne viruses, as we have noted above, and (b) pollutants in the atmosphere. First-passage-time problems have been studied extensively [20][21][22][23] and they have found applications in a variety of areas in physics and astronomy, chemistry [24], biology [25], and finance [26]. In the fluid-turbulence context, different groups have studied zero crossings of velocity fluctuations [27] or various statistical measures of two-particle dispersion, including exit-time statistics for such dispersion in two-and three-dimensional (2D and 3D) turbulent flows [28,29]. In contrast to these earlier studies (e.g., Refs. [28][29][30]), the first-passage-time problem we pose considers one tracer in a turbulent flow that is statistically homogeneous and isotropic. To the best of our knowledge, this first-passage problem has not been studied hitherto. For such a particle we show, via extensive direct numerical simulations (DNSs), that P(R, t R ) displays a crossover between two qualitatively different behaviors: (a) for R L I , P(R, t R ) ∼ t −α R , with L I the integral scale of the turbulent flow and the exponent α = 4; (b) for L I R, P(R, t R ) has an exponentially decaying tail (Fig. 2). We develop models that allow us to obtain these two asymptotic behaviors analytically. Most important, we show how to use P(R, t R ) to obtain estimates of social-distancing guidelines for the mitigation of the spreading of airborne aerosols with viruses such as SARS-CoV-2.

MODELS, METHODS, AND RESULTS
The 3D incompressible, Navier-Stokes equation is and Here, u(x, t) is the Eulerian velocity at position x at time t, p(x, t) is the pressure field, and ν is the kinematic viscosity of the fluid; the constant density is chosen to be unity. Our direct numerical simulation (DNS) uses the pseudo-spectral method [31], with the 2/3 rule for dealiasing, in a triply periodic cubical domain with N 3 collocation points; we employ the second-order, exponential, Adams-Bashforth scheme for time stepping [32]. We obtain a nonequilibrium, statistically stationary turbulent state via a forcing term f , which imposes a constant rate of energy injection [33,34], in wave-number shells k = 1 and k = 2 in Fourier space; this turbulent state is statistically homogeneous and isotropic.
To obtain the statistical properties of Lagrangian tracers, which are advected by this turbulent flow, we seed the flow with N p independent, identical, tracer particles. If the Lagrangian displacement of a tracer, which was at position r 0 at time t 0 , is r(t | r 0 , t 0 ), then its temporal evolution is given by where v is its Lagrangian velocity. In Eq.
(2), we need the Eulerian flow velocity at off-grid points; we obtain this by tri-linear interpolation; and we use the first-order Euler method for time marching (see, e.g., Ref. [32]). We give important parameters for our DNS runs in Table I. These include the time step dt, the en- is the integral length scale and T eddy = L I /u rms is the integral-scale eddy-turnover time; η = (ν 3 / ) 1/4 and τ η = (ν/ ) 1/2 are, respectively, the Kolmogorov dissipation length and time scale; and k max is the maximum wave number that we use in our DNS.
Clearly, t R is the first time at which |r| becomes equal to R. Instead of computing the PDF (or histogram) of t R numerically, we calculate the complementary cumulative probability distribution function (CPDF) Q(t R ), by using the rank-order method [35], to circumvent binning errors. In Fig. 2, we present log-log and semi-log plots of ; note that, in this power-law scaling regime, the complementary CPDFs for different values of R/L I collapse onto a universal scaling form, if we plot Q( t R /T eddy R/LI ). In contrast, Fig. 2 (b) shows that, for L I R, the tail of Q(t R /T eddy ) decays exponentially. For the first-passagetime PDF, these results imply that (3) We now develop models that allow us to understand these two asymptotic behaviors analytically.   I. Parameters (see text for definitions) for our DNS runs: N 3 is the total number of collocation points; ν is the kinematic viscosity; dt is the time step; Re λ is the Taylor-microscale Reynolds number; is the energy-dissipation rate; η is the Kolmogorov dissipation length scale; kmax is the maximum wave number that we use; LI is the integral length scale; T eddy is the integral-scale eddy-turnover time; τη = (ν/ ) 1/2 is the Kolmogorov dissipation time scale; and Np is the number of tracer particles.
For the power-law behavior of P(R, t R ), in the range R L I , we construct the following, natural, ballistic model: Tracer particles emanate from the origin with (a) a velocity whose magnitude v is a random variable with a PDF p(v); and (b) when it starts out from the origin, the tracer's velocity vector points in a random direction. Tracers move ballistically, for short times. Therefore, for R L I , the first-passage time t R = R/v; and the firstpassage PDF is In statistically homogeneous and isotropic and incompressible-fluid turbulence, each component of the Eulerian velocity has a PDF that is very close to Gaussian [36], so p(v) has the Maxwellian [37] form where C d depends on the spatial dimension d, and C d = 4π for d = 3, and σ = v 2 = u rms . We substitute Eq. (5) in Eq. (4); then, by integrating over v, we obtain Therefore, in the limit of small R and large t R , the firstpassage-time probability is this power-law exponent is the same as the one we have obtained from our DNSs above (Table I and Fig. 2). We can obtain the tail P(R, t R /T eddy ) ∼ exp(−(t R /T eddy )) for L I R as follows. At times that are larger than the typical auto-correlation time of velocities in the Lagrangian description, we follow Taylor [38] and assume that the motion of a tracer particle is diffusive. Therefore, we consider a Brownian particle in three dimensions (3D). To calculate the first-passage-time PDF, we must first obtain the survival probability S(t, R|0), i.e., the probability that the particle has not reached the surface of the sphere of radius R up to time t, if it has started from the origin of this sphere. We start with the forward Fokker-Planck equation [22,39] for the PDF of finding the particle at a distance r from the origin at time t: where K is the diffusion constant; this PDF satisfies the initial condition, P (r, 0) = δ(r)/(4πr 2 ) and the absorbing boundary condition P (R, t) = 0, for all t at r = R. We obtain the following solution: where, in the last step, we have used Eq. (9). The firstpassage-time probability is At large times, the first term (n = 1) is the dominant one; therefore, the exponential form that we have obtained from our DNS (Fig. 2 (b)); the 1/R 2 pre-factor cannot be extracted reliably from our DNS data, because this requires much longer runs than are possible with our computational resources. We now show that both the small-and large-R/L I behaviors of P(R, t R )(R, t R ) in Eq. (3) can be obtained from one stochastic model for the motion of a particle. The simplest such model uses a particle that obeys the following Ornstein-Uhlenbeck (OU) model: Here, γ and Γ are positive constants; x i and v i are the Cartesian components of the position and velocity of the particle; in three dimensions, i = 1, 2, and 3; ζ i (t) is a zero-mean Gaussian white noise with ζ i = 0 and ζ i (t)ζ j (t ) = δ ij δ(t − t ); this noise is such that the fluctuations-dissipation theorem (FDT) holds. Note that there is no FDT for turbulence. However, for the one-particle statistics that we consider, the simple OU model is adequate. We use N p = 50, 000 particles; for each particle, the initial-position components x i (t = 0) are distributed randomly and uniformly on the interval [0, 2π]; and the velocity components v i (t = 0) are chosen from a Gaussian distribution. For each particle, we obtain, numerically, the time t R at which it reaches a distance R from the origin for the first time. We then obtain the first-passage-time complementary CPDF Q(t R ), which we plot in Fig. 3, for R L and L R, where L = Γ γ 3 , the natural length scale for Eq. (13), plays the role of L I in our DNSs above (Table I and Fig. 2). We find these are the OU-model analogs of our DNS results Eq. (3). We have carried out two OU-model simulations: (a) we have designed the first, with γ = 0.01, to explore the form of P(R, t R ) in the ballistic regime R L; (b) the second, with γ = 30, allows us to uncover the form of P(R, t R ) in the diffusive regime L R. (From a numerical perspective, it is expensive to obtain the precise form of P(R, t R ) in both ballistic and diffusive regimes, with one value of γ.) We now explore in detail the forms of P(R, t R ) in these two regimes. In Fig. 3(a), we present log-log plots of the complementary CPDFs of the scaled first-passage time t R /R, for R L and γ = 0.01. The complementary CPDFs of t R /R, for R/L = 0.0002, R/L = 0.00035, and R/L = 0.0005, collapse onto one curve; i.e., in this regime, t R scales as R, which is a clear manifestation of ballistic motion. In Fig. 3(b), we present semi-log plots of the complementary CPDFs of the scaled first-passage time t R /R 2 , for L R and γ = 30. The complementary CPDFs of t R /R 2 , for R/L = 10, R/L = 14, R/L = 18, and R/L = 20, collapse onto one curve; from this we conclude that, in this regime, t R scales as R 2 , which is a clear signature of diffusive motion.

CONCLUSIONS AND DISCUSSION
We have defined and studied a new first-passage-time problem for Lagrangian tracers that are advected by a 3D turbulent flow that is statistically steady, homogeneous and isotropic. Our work shows that the first-passage-time PDF P(R, t R ) has tails that cross over from a power-law form to an exponentially decaying form as we move from the regime R L I to L I R (Eq. (3)). We develop ballistic-transport and diffusive models, for which we can obtain these limiting asymptotic behaviors of P(R, t R ) analytically. We also demonstrate that an OU model, with Gaussian white noise, which mimics the effects of turbulence, suffices to obtain the crossover between these limiting forms. Of course, such a simple stochastic model cannot be used for more complicated multifractal properties of turbulent flows [29,36,40].
Earlier studies have concentrated on two-particle relative dispersion by using doubling-time statistics, in 2D fluid turbulence; in particular, they have shown that the PDF of this doubling time has an exponential tail [28]. Studies of velocity zero crossings [27], in a turbulent boundary layer, have shown that PDFs of the zerocrossing times have exponential tails.
The single-particle first-passage-time statistics that we study have not been explored so far. Furthermore, P(R, t R ) can be used to to develop social-distancing guidelines for the mitigation of the spreading of airborne aerosols with viruses such as SARS-CoV-2 as we show below.
Given a pseudospectral DNS, of the type we have carried out, we can obtain the integral scale L I and u rms from the energy spectrum E(k), as we have noted above. A recent study of COVID-19 in 320 municipalities in China suggests that a very large fraction of COVID-19 infections occur because of indoor transmission of the SARS-CoV-2 virus [11]. Therefore, it is important to study such transmission in rooms and offices; a comprehensive DNS study of the Navier-Stokes equation, with the correct boundary conditions enforced at every wall and surface in the room and accurate forcing functions (dictated, e.g., by fans and vents), is a considerable challenge. Furthermore, it is not possible to carry out such a DNS for every room with a different arrangement of the furniture in it. Hence, it is important to come up with semi-quantitative criteria that help us to understand, and mitigate, the indoor transmission of such virues. Turbulence models have been used to study the flow of air in rooms and offices [41,42]; from these models and related experiments we obtain the estimate u rms 0.05 m/s in a typical office. We must also estimate L I , for it is an important crossover length scale in our analysis of P(R, t R ). In our DNS L I is 0.1L, where L is the linear size of our simulation domain. In a typical office or a train, with fixed forcing, via fans or vents, we use L I to be approximately a few meters; of course, L I must depend on the degree of crowding on a train or the number of cubicles in a large office room. Now consider one infected person who is at a distance R from another person. The probability of virus-laden aerosol particles not reaching the second person, up until time t is related to P(R, t R ) as follows: which we calculate, by using Eq. (11), and depict in  Table II we give the values of W (R, t) for different values of R and t. These figures and Table II lead to three clear observations: 1. If the separation R L I , i.e., we have to consider the ballistic regime (see the Supplementary Material), then W (R, t) goes very rapidly to 0 (i.e., the aerosol particle reaches the second person), even if t is very small.
2. The smaller the separation R, between two persons, the shorter the time t in which W (R, t) becomes very small, i.e., the aerosol particles reach from one person to the other.
3. Our calculation leads to quantitative predictions; e.g., if the separation L I R, i.e., we have to con-sider the diffusive regime, then W (R, t) goes to 0 in tens of seconds, if R = 2 m, and in hundreds of seconds, R = 10 m, for the representative parameters that we use to obtain Table II. A recent study [47] has suggested that the SARS-CoV-2 virus remains viable in aerosols for nearly 3 hours. Therefore, if the concentration of virus-laden aerosols is high in a poorly ventilated room, then we must employ more stringent social-distancing norms than are in place now, even if people spend only tens of minutes together in such a room.
The methods that we have developed can be applied, mutatis mutandis, (a) in sophisticated models for virus particles or droplets, e.g., those that use inertial particles [43,44] or multi-phase flows [45,46] and (b) in turbulent flows that are not statistically homogeneous and isotropic. We will examine these in future work. At the moment, it is important to use our minimal model to obtain semi-quantitative for social-distancing guidelines, as we have done above.
D.M. and A.B. thank John Wettlaufer for useful discussions. A.K.V. and R.P. thank Jaya Kumar Alageshan for discussions and for help with   , t), the probability of virus-laden aerosol particles not reaching a person (at a distance R from an infected person), up until time t (the first-passage time) in the diffusive region, for the representative values K = 0.1m 2 /s, urms = 0.05 m/s, and LI = 2m.