Turbulent transport of droplets in a jet/puff is a problem of paramount importance in science and engineering that nowadays has become even more important given the global COVID-19 emergency caused by the SARS-CoV-2 infection; for a recent review see [1] [2] [3] . The relationship between COVID-19 and turbulent transport of droplets stems from the fact that the dominant route of SARS-CoV-2 spread is via small virus-containing respiratory droplets that the infected person exhales when coughing, sneezing or talking [4] . Different factors make the fate of liquid droplets hard to predict. The exhalation is far from a homogeneous fluid. It rather consists of a two-phase mixture of liquid droplets dispersed into a fluid phase which is usually at a higher temperature and humidity than the ambient air. Evaporation thus occurs making the droplets lighter (and less inertial) than they were at the expulsion stage. The exhaled air is also turbulent because of the large velocities reached during violent expiratory events. The Reynolds number is about 10 4 [5, 6] for cough and even larger (of about a factor 4) for sneeze [6] . Turbulence also characterizes the space-time evolution of the humidity field which, in turns, affects the droplet evaporation. In this Letter, we aim at elucidating the role of turbulence on the trajectories and on the evaporation of the respiratory droplets in a cough. By means of fullyresolved direct numerical simulations (DNS) we show that turbulence strongly affects droplet evaporation provided they are sufficiently small: with respect to coarsegrained descriptions (including the mean-field approach where turbulence is neglected altogether), fully-resolved turbulent fluctuations cause a delay of the evaporation process which, in turn, causes droplets to remain heavier for longer. This delay increases the inertia of the liquid droplets which fail to follow the initial accelerated phase [5] of the exhaled air. The net result is that droplets spuriously remain airborne for longer times if turbulence is not properly accounted for. Because virus-containing airborne particles are responsible for the SARS-CoV-2 spread, turbulence has a direct impact on the 'social distancing' issue [7] .
Air exhaled from the mouth is ruled by the incompressible Navier-Stokes equations
with ν being the air kinematic viscosity and ρ a the air density. The list of all relevant parameters used in this study are reported in [8] .
Instead of simulating the evolution of the absolute humidity field (the exhaled air is saturated, or close to saturation [9] ) it is more convenient to model directly the supersaturation field (i.e. s = RH − 1, RH being the relative humidity). Indeed, the supersaturation dictates the evaporation/condensation process, as it appears in the evolution equation for droplet radius [10] . The supersaturation field is ruled by the advectiondiffusion equation [11] : for the velocity field in Ref. [12] which also holds for the supersaturation field.
D v being the water vapor diffusivity. Equation (2) assumes that the saturated vapor pressure is constant, an assumption that holds as long as the ambient is not much colder than the exhaled air, which is at about 30 o C according to [9] . We assume an inlet air velocity profile representative of cough (supposed to be saturated, i.e. s = 0) as reported in [5] , exiting from a mouth opening of area of 4.5 cm 2 . The duration of the expulsion is 0.4 s and the peak velocity is 13 m/s. The resulting Reynolds number (based on the peak velocity and on the mouth average radius) is about 9000. The flow field is thus fully turbulent as one can easily realize by looking at Fig. 1 . Before discussing how the liquid part of the two-phase mixture is modelled, let us first validate the puff dynamics of the exhaled air. Kovasznay et al. [12] predicted the powerlaw behavior t −3/4 for the self-similar decay of the puff velocity. Simple scaling arguments applied to a passive scalar advected by the turbulent puff lead to the same t −3/4 power law. This scaling is thus expected to hold for the decay of the supersaturation field. The reliability of our puff dynamics is demonstrated in Fig. 2 which clearly shows the expected scaling laws for more than two decades with high accuracy.
We are now ready to introduce the model for the liquid part of the two-phase mixture. It is described as an ensemble of N inertial particles ruled by the well-known set of equations [13] 
(4) where N is the number of exhaled droplets (here N∼ 5000 according to [14] ), X i is the position of the i-th droplet and U i its velocity, and, finally, g is the gravitational acceleration. Each droplet is affected by a Brownian contribution via the white-noise process η i . Here, ρ D i is the density of the i-th droplet. Because the volume fraction of the liquid phase for cough is always smaller than 10 −5 [6, 15] , the back-reaction of the droplets to the flow is irrelevant. Droplets are assumed to be made by salty water (water and NaCl) and a solid insoluble part (mucus) [16] [17] . Finally, τ i is the Stokes time of the i-th droplet and R i is its radius. Droplet radii evolve according to the ruling equation [10] 
No feedback of this equation to Eq. (2) is considered here because of the very small values of the liquid volume fraction we have already stated above. In Eq. (5), C R is the droplet condensational growth rate, s crh = −0.55 (CRH = 0.45, the so-called crystallization RH or efflorescence RH) for the NaCl [18] . Fig. 3 of [19] and Ref. [20] show the weak dependence of CRH on temperature. r N i is the radius of the (dry) solid part of the i-th droplet when the salt is totally crystallized (i.e. below CRH). The dependence of r N i on physical/chemical/geometrical properties of the exhaled droplets is reported in [21] together with the expressions of parameters A and B. On the basis of the parameters assumed here, the ratio r N i /R i (0) is 0.16 [22] which agrees with the estimations discussed in [23] . As far as the initial size of each exhaled droplet is concerned, we assume here for the sake of example, the one from the seminal paper by Duguid [14] , a paper still considered as a reference report on the subject. Exhaled droplets enter the ambient considered at rest with a relative humidity RH = 60% (i.e. s a = −0.4), larger than the crystallization RH. States of local equilibrium are possible from Eq. (6) owing to the solute effect [10] . Because the supersaturation field evolves as a passive scalar in a turbulent field, it exhibits the well known "plateaux-and-cliffs" structures [24] [25] [26] [27] [28] . Accordingly, the scalar field is turbulent with very strong fluctuations occurring in small regions (called cliffs or fronts) separating larger areas where the scalar is well mixed (also called plateaux). Because small droplets and supersaturation are transported by the same velocity field (this does not occur for the large droplets which are affected by inertial effects), correlations occur between droplet trajectories and supersaturation values [11] . This phenomenon causes droplets of sufficiently small size to stay long in the large well-mixed regions where they can locally equilibrate with the (local) value of the supersaturation. The droplet evaporation process is thus expected to behave in time by alternating phases of equilibrium with phases of rapid evaporation: a sort of stop-and-go process. The same type of structures are expected also for the decay of droplet radii. This phenomenon can be clearly detected in Fig. 3 where the temporal behavior of the supersaturation field along the Lagrangian trajectory of a small droplet is reported (group of lines denoted by St < 1, St being the Stokes number) together with the time evolution of the corresponding droplet radius (inset). The time history with the fully resolved DNS (blue, continuous line) clearly shows the effect of the plateaux-and-cliffs structures on the evaporation process which is however absent for the larger droplet affected by inertia (group of lines denoted by St > 1). The fact that the radius closely follows the temporal behavior of the supersaturation field (see inset of Fig. 3) is the signature of a quasi-adiabatic picture for the evaporation process (i.e. the process of radius adjustement due to evaporation is much faster than the corresponding variation of the supersaturation field). It is worth noting that if one considers the smaller droplet evolving in coarse grained fields (long dashed line in red, where both velocity and supersaturation have been coarse grained in space as discussed in [29] ), the effect of the plateaux-and-cliffs structures on the evaporation process reduces and vanishes when the turbulent fields are replaced by their mean field components (green dashed line). Having shown that sufficiently small droplets correlate with the supersaturation field, let us now discuss the consequences on droplet motion. For smaller droplets remaining for a sufficiently long time in regions where the supersaturation field is locally constant, with a value larger (smaller) than the mean, the evaporation takes place more slowly (rapidly) than what would be for the same droplet experiencing smoother fluctuations as in the filtered DNS or in the mean-field approach. The two effects, i.e. reduction vs increase in evaporation time, are however not symmetric due to the fact that the supersaturation field is decaying in time as we have already shown in Fig. 2 . The net result caused by turbulent fluctuations on the fate of small droplets is thus to increase their evaporation time.
Let us now quantify the delay caused by turbulence in the evaporation process by comparing, for an observation time of 60 s, the time it takes for each airborne droplet to shrink to their final equilibrium radius. Let us denote those typical evaporation times as τ evap . All droplets which sedimented within the observation time of 60 s were not included in this analysis. The sole airborne droplets were selected here, thus automatically satisfying the requirement of having a sufficiently small radius. The results are presented in Fig. 4 where the probability density functions of τ evap are reported both for the fully resolved case and for the evolution with the sole mean fields (of both the carrying flow and the supersaturation field) and with the filtered DNS. The corresponding mean evaporation times are reported in Tab. I. The role of turbulence clearly emerges, both causing broader probability density functions, the fingerprint of fluctuations, and to delay the evaporation process. Finally, we show that the observed delay in the evaporation affects droplet motion. This is depicted in Fig. 5 where we report, for the same airborne droplets considered to evaluate the probability density functions of τ evap , the streamwise coordinate, x(t), of the center of mass of the cloud of airborne droplets as a function of time. Shown in this figure are the fully resolved DNS, the filtered DNS, and the mean-field approach. In the two cases where turbulent fluctuations are either coarse grained or entirely neglected, droplets travel further than in the fully resolved DNS. This is the fingerprint of the reduced inertia of the droplets evolving in the filtered fields.
In the initial stage of their evolution, these droplets are indeed spuriously lighter than the droplets evolving in the fully-resolved DNS. Being lighter, they are carried more efficiently by the underlying rapidly accelerating flow thus reaching longer distances before touching the floor.
In order to ascertain whether the observed delay of trajectories of small droplets is a genuine effect caused by the interplay between turbulence and inertia, a subset of idealized simulations have been performed where monodisperse droplets of R i (0) = 5 µm have been considered, with and without inertia (i.e. simply switching on/off inertia in the ruling equations (3) and (4)). This size is close to the peak of the droplet size distribution of [14] we have used in the previous analysis, and corresponds to droplet neither too large to be insensitive to turbulence, nor too small to make the mass loss due to evaporation negligible. The results are shown in Fig. 6 . Both in the presence and in the absence of droplet inertia we found the turbulence-induced broadening of the probability density functions of the evaporation time. The inset of Fig. 6 shows this fact for the simulations without inertia. Filtering the turbulence fluctuations (long-dashed black curve in the inset) reduces the broadening as observed for the polydisperse case with inertia. It is now worth remarking that the observed difference between the mean evaporation time measured from the DNS and the one measured from the filtered DNS does not produce any relevant effect on the droplet motion when inertia is switched off in the droplet ruling equations. The similarity in the main frame between the continuous gray curve and the black long-dashed curve confirms this fact. Switching-on the inertia, the effect of the delayed evaporation in the DNS case becomes apparent (see in the main frame the differences between the continuous blue curve and the red long-dashed curve). Fig. 6 confirms that turbulence is the root cause of the broadening of evaporation times, whereas inertia causes differences in the trajectories.
In conclusion, turbulence increases the droplet evaporation time thus reducing the flight time of airborne droplets. Because this implies shorter distances travelled before reaching the floor, turbulence enters forcefully into the current debate on the 'social distancing' issue.