First-passage-time problem for tracers in turbulent ﬂows applied to virus spreading

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


I. 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 optimize social-distancing guidelines we must ask the following: How far, and how fast, do small respiratory droplets or virus-bearing aerosols spread in turbulent flows?Given the ongoing COVID-19 pandemic, it is important and extremely urgent to have at least a semiquantitative 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 μm, 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 particles 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 fluiddynamics 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 center 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.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; and (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.

II. 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 pseudospectral 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 trilinear 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 energy dissipation rate is the integral length scale and T eddy = L I /u rms is the integral-scale eddyturnover time; η = (ν 3 / ) 1/4 and τ η = (ν/ ) 1/2 are, respectively, the Kolmogorov dissipation length and timescale; 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 semilog plots of Q(t R ) vs t R /T eddy for several values of R. From Fig. 2(a) we conclude that, for , for large t R /T eddy , with α 4; 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/L I ).In contrast, Fig. 2(b) shows that, for L I R, the tail of Q(t R /T eddy ) decays exponentially.For the first-passage-time PDF, these results imply that (3) We now develop models that allow us to understand these two asymptotic behaviors analytically.
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 first-passage 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 autocorrelation 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, whence we get where, in the last step, we have used Eq. ( 9).The first-passagetime 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 prefactor 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 this noise is such that the fluctuation-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.000 35, 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 semilog 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.

III. CONCLUSIONS AND DISCUSSION
We have defined and studied a 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 zero-crossing times have exponential tails.
The single-particle first-passage-time statistics that we study still needs to be explored.Furthermore, P (R, t R ) can be used 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 furniture in it.Hence, it is important to come up with semiquantitative criteria that help us to understand and mitigate the indoor transmission of such viruses.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 Figs.

FIG. 1 .
FIG. 1.A schematic diagram illustrating how aerosols with viruses (red points) may be advected by a turbulent flow, from the person at the center (the source) to other persons at different distances from the center.

FIG. 2 . 3 (
FIG. 2. Plots of the complementary cumulative probability distribution functions (CPDFs) Q vs the scaled first-passage time t R (see text): (a) Log-log plots of Q( t R /T eddy R/L I ) for R/L I = 0.06 (blue), R/L I = 0.10 (purple), R/L I = 0.14 (green), and ( t R /T eddy R/L I ) −3 (black dashed line); the inset shows log-log plots of Q(t R /T eddy ) for the same values of R/L I .(b) Semilog plots of Q(t R /T eddy ) for R/L I = 1.02 (pink), R/L I = 1.43 (blue), R/L I = 1.84 (purple), and R/L I = 2.04 (green).

FIG. 3 .
FIG. 3. (a) Log-log plots of the complementary CPDFs Q( γ t R R/L ) of the scaled first-passage time γ t R R/L , for R L and γ = 0.01; the complementary CPDFs, for R/L = 0.0002 (green), R/L = 0.000 35 (blue), and R/L = 0.0005 (orange), collapse onto one curve; (b) semilog plots of the complementary CPDFs of the scaled first-passage time t R /R 2 , for L R and γ = 30.The complementary CPDF of t R /R 2 , for R/L = 10 (purple), R/L = 14 (green), R/L = 18 (blue), and R/L = 20 (orange), collapse onto one curve.Plots of the complementary CPDFs Q(γ t R ) vs γ t R are shown in the insets.

FIG. 4 .
FIG. 4. Surface plots of W (R, 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 (a) vs R and t, and (b) vs R/L I and t/T , in the diffusive region, for the representative values K = 0.1 m 2 /s, u rms = 0.05 m/s, and L I = 2 m.
4(a) and 4(b), in the diffusive regime; Fig. 4(a) is a surface plot of TABLE II.Table of values of W (R, 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.1 m 2 /s, u rms = 0.05 m/s, and L I = 2 m.

TABLE 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; k max is the maximum wave number that we use; L I is the integral length scale; T eddy is the integral-scale eddy-turnover time; τ η = (ν/ ) 1/2 is the Kolmogorov dissipation time scale; and N p is the number of tracer particles.