Contact statistics in populations of noninteracting random walkers in two dimensions
FFebruary 25, 2021
Contact statistics in random walks: Noninteracting walkers
Mark Peter Rast
Department of Astrophysical and Planetary Sciences,Laboratory for Atmospheric and Space Physics, University of Colorado, Boulder, CO 80309,USA
We examine the contact statistics of noninteracting random walkers on a two-dimensional plane.We find that the waiting time between successive contacts for an individual (the interarrival time) isnon-exponentially distributed, with consequent non-Poison distributed contact counts. This mani-fests as over-dispersion of the contact-count probability mass function, caused by negative durationdependence of the waiting time interval when the walkers are in close proximity. While succes-sive contacts are independent, the probability of repeat contact between individuals decreases withtime after a previous contact. Further, contact duration depends sensitively on walker radius. Forwalker radii smaller than both the mean step length and the mean nearest neighbor separation,the probability density of the contact duration is strongly peaked at the ballistic crossing time forhead-on collisions. The contact statistics of individuals in biological populations, dilute componentsof chemical systems, or Lagrangian particles in turbulent flows depends critically on the behaviorat close proximity, whether it is attractive, repulsive, or neutral. This work clarifies those statisticsfor the neutral case, under the assumption that the underlying motions approximate a collection ofrandom walkers, and we suggest it may be useful in understanding contagion in human populations.
PACS numbers: 47.27.-i, 47.32.-yKeywords: contact statistics, random walks
I. INTRODUCTION
This is the first in a series of papers motivated by theSARS-CoV-2 pandemic that investigate contact statis-tics in populations with different underlying internal dy-namics. This first paper focuses on the simplest casein which each individual in the population undertakes arandom walk in two dimensions. Follow-on work looksat how the contact statistics change with more complexdynamics, first for particles embedded in an idealizedtwo-dimensional turbulent flow, and finally for humanmotions as captured by GPS positioning data.Modeling the contact interval between dynamical tra-jectories is important to a wide range of biological andphysical systems. Examples include contagion, preda-tion and reproduction in biological systems, moleculardynamics and reaction rates in chemical systems, andcloud droplet formation and particle aggregation in tur-bulent atmospheric and astrophysical systems. Some im-portant statistical properties to understand are the dis-tribution of waiting times between contacts (also calledthe interarrival time), the distribution of contact occur-rences over specified temporal intervals (count statistics),and the distribution of contact durations once they oc-cur. These address how often contact occurs, how manycontacts are achieved, and how long each contact is typ-ically sustained.While seemingly elementary, we find important sub-tleties in these measures even when the motions of indi-vidual members of a population follow simple random-walk dynamics. Multiple length-scales define a collectionof random walks: the step-size, the mean separation ornumber density of walkers, and the interaction distance(the distance between walkers at which contact occurs,equal to twice the walker radius). These can vary in their hierarchical relationships, and while in the simplerandom walk model presented here they are input pa-rameters to be varied, in more complex systems they areconstrained by the underlying dynamics of the motions.
II. RANDOM WALK MODEL
We look at contact statistics between random walk-ers on a periodic two-dimensional plane of unit length inwidth and height. Each successive step of each randomwalk occurs in a random direction (uniformly distributedin angle between zero and 2 π ) and has a random length(uniformly distributed between zero and a specified max-imum). The position of each walker is well resolved bynumerically advancing their positions in small sub-steps,and a small random correction is made to the size of thefinal sub-step to avoid having multiple random walk tra-jectories make directional changes simultaneously, whichwould otherwise occur due to the discrete nature of thenumerical steps. In other words, small equal-length nu-merical sub-steps are taken in the same direction for aspecified total number that is uniformly distributed be-tween zero and the maximum step-length. In more com-plex flows, the step length can be thought of as the au-tocorrelation length of the velocity along a trajectory (orthe Lagrangian integral scale) of the flow.We have carefully checked that none of the conclu-sions drawn in this paper depend on the details of howthe random walks were constructed or on the imposeddomain periodicity in which they occur. We cite some ofthose checks and instances of solution sensitivity in thediscussions below. The statistical properties reported inthis paper are the same as those which would be ob-tained in an infinite rather than periodic domain using a r X i v : . [ m a t h - ph ] F e b isotropically directed steps of constant step-size equal tothe mean of our step-size distribution, so long as thenumerical sub-steps employed are sufficiently small to re-solve contacts. The governing parameters of the multiple-walker random-walk simulations are thus the walker ra-dius a and mean step length δr , which are taken commonto all the walkers, and the number density of walkers n , or equivalently their mean nearest-neighbor separa-tion, which is equal to 0 . / √ n (see derivation followingA. Alexakis in [1]). It is the relative magnitudes of thesequatities that govern the solution statistics. We notethat the random walks are unchanged by contact betweenwalkers. No interaction between walkers is modeled. III. THE DISTANCE BETWEEN WALKERS
Underlying the contact statistics of random walks arethe statistics of pair separation. For example, the distri-bution of the waiting time interval between contacts (theinterarrival time distribution) for any individual can bethought of as the minimum statistic, over all other in-dividuals, of the first passage time for a specified pair-separation contact-distance. Working directly in two-dimensions proves problematic, for while a direct solu-tion of the diffusion equation readily yields the randomwalk first-passage time distribution to a point or spherein one or three dimensions (the Levy distribution) it failsto yield such a convenient solution in two-dimensions(see Appendix). A fundamental difficulty arises because,while the motion of each individual on a two dimen-sional plane can be described as a two-dimensional ran-dom walk, the distance between any two walkers onlyapproximates a one-dimensional constant diffusivity pro-cess in the short and long time limits. At intermediatetimes the pair separation probability distribution evolvesfrom approximately Gaussian to Rayleigh with a corre-sponding nonlinear change in its variance.In two dimensions, the distance r of a single unbiasedrandom walker from its origin is Rayleigh distributed, p ( r ) = rσ exp (cid:18) − r σ (cid:19) , (1)with the displacement variance σ r = (4 − π ) / σ . Inthe continuous time and space limit, the displacementdistribution evolves diffusively, with the diffusion coef-ficient D = δr / δt when steps of mean length δr aretaken isotropically in random directions. The displace-ment variance σ r scales linearly with elaspsed time sincethe start of the random walk. It is constant multiple ofthe scale parameter σ = 2 D t .The distance r s between two unbiased random walkers,as a function of their initial separation r at t = t , onthe other hand, is Rice distributed p ( r s | r ) = r s σ exp (cid:18) − r s + r σ (cid:19) I (cid:16) r r s σ (cid:17) , (2) FIG. 1: Standard deviation of the pair separation distributionas a function of time (in units of the number of steps taken)for random walkers of varying step length. Curves of solid col-ors indicate pair separation results for walks with step lengths( left to right ) δr = [0 . , . , . , . , . r = 0 . Black curve indicates stan-dard deviation of the separation for a pair with step length δr = 0 . r = 0 . Dotted curvesof the same color indicate the expected standard deviationsin the diffusion limit as discussed in main body of text. The orange and blue dashed curves plot the pair separation dis-tribution standard deviation for two pairs of random walkerswith step length δr = [0 . , . , ] and initial separation r = [0 . , . dotted vertical fiducial lines indicatethe time of first possible contact between the pair and theapproximate time of transition to the long-time Rayleigh dis-tribution (see Figure 2 and discussion in text). with I denoting the lowest order modified Bessel func-tion of the first kind [2]. The pair-separation variance(the variance of the Rice distribution) is σ s = 2 σ + r − πσ L / (cid:0) − r / σ (cid:1) , (3)where L / indicates the square of the L / Laguerrepolynomial. In the continuous time and space limit eachwalker individually diffuses from its initial position (asabove with D = δr / δt ), but the variance in their sepa-ration σ s is a nonlinear function (Equation 3) of the scaleparameter σ = 4 D ( t − t ), and thus time.The standard deviations of the separation betweenpairs of random walkers are plotted in Figure 1 as afunctions of time, with time measured in units of thenumber of steps taken (i.e., δt = 1). Curves of solid color show values from simulations of many initially wellseparated ( r = 0 .
7) random-walker pairs, with colors in-dicating the results for different mean random-walk steplengths decreasing from left to right . As expected, thedistance between two random walkers scales ballistically, σ s ∼ t , for times short compared to one step and diffu- FIG. 2: Snapshots of the pair separation distance probabil-ity density as a function of time for pairs of walkers whoseseparation variance is plotted with the dashed dark blue curvein Figure 1 at t − t = [0 . , . , . , r / (2 δr ) = 4 .
76] (narrowto wider pdfs respectively) in upper panel, and for t − t =[10 . , . , r /δr = 363] (narrow to wider pdfs respectively)in lower panel. In upper panel, black curves indicate best fitGaussian distributions, and in the lower panel black curvesindicate best fit Rice distributions. Best fit Rayleigh distri-bution are shown as red curves in the lower panel. sively, σ s ∼ t / , for times greater than this. These scal-ings are indicated with thin black dashed lines in Figure 1for reference, and the strict diffusive limit for each sim-ulation, σ s = σ = (cid:112) D ( t − t ), is plotted as a dotted curve of the same color. The offset between the theo-retical diffusive limit and numerically determined result(between the dotted lines and and solid curves at longtimes) reflects the discrete, initially ballistic, motion ofthe simulated walkers, which is resolved numerically fortimes less than one by taking small sub-steps. The sizeof the numerical sub-steps is constant across the simula-tions shown.Two of the numerical pair-separation experiments (thesecond and third from the top in Figure 1) were re-run with a smaller initial pair separation ( r = 0 . dashed dark or-ange and dashed dark blue linestyles. These two curvesillustrate how all the others would behave if extended to longer times. The curves deviate from diffusive scal-ing at t − t ≈ r / (2 δr ) (marked by the two lefthandshort dashed vertical fiducial lines). This is the ear-liest possible time that the distance between the tworandom walkers can equal zero. Diffusive scaling is re-established later, but with a reduced variance equal to σ s = (4 − π ) / σ (shown with dottted dark orange and dotted dark blue curves). By Equation 3, this occurs when r (cid:28) σ , and the approximate time at which this occurs t − t ≈ r /δr is marked by the two righthand short dashed vertical fiducial lines.The change in variance between t − t ≈ r / (2 δr ) and t − t ≈ r /δr reflects a change in the underlying sep-aration probability density function which in two dimen-sions depends critically on both the initial separation andthe elapsed time. Figure 2 shows the temporal evolutionof the normalized pair separation distribution for pairsof walkers whose separation variance is plotted with the dashed dark blue curve in Figure 1. For times shorterthan one step, the probability density function is non-Gaussian ( grey curve in top panel of Figure 2), reflectingthe ballistic separation of the walkers r s = r + 2 r δr t (cos θ − cos θ ) − δr t cos( θ − θ ) , (4)with θ and θ measured relative to the radius vectorconnecting the two walkers and uniformly distributedbetween 0 and 2 π and t measured in steps. This bal-listic phase will be important in understanding the con-tact duration distributions (Section V B). For somewhatlonger but still short-times, 1 . (cid:46) t − t ≤ r / (2 δr ) (illus-trated by the remaining, non- grey curves in upper panel),the pair separation is distributed as a truncated Gaus-sian, with the truncation occurring at a value of r = r ± δr ( t − t ), the maximum and minimum separationsthat the walkers can achieve over the elapsed time. Dur-ing this period the separation distribution variance scalesdiffusively, with value equal to twice that of the individ-ual walker displacement variance, σ s = σ = 4 D ( t − t ).For times longer than t − t = r / (2 δr ) ( lower panel),the pair separation is Rice distributed, with the meanshifting to larger separations with time and the variancegrowing more slowly than t because separation values re-flect across r = 0. At long times, t − t (cid:38) r /δr , theRice distribution asymptotes to Rayleigh, with the vari-ance again scaling with t − t σ s = (4 − π ) / σ . These behaviors make pair-separationdifficult to model as a one-dimensional random variable,but also provide an opportunity for simplified analysis inthe short and long time limits.It is important for the analysis of the multiple random-walker simulation results below to note that, for ini-tial separations smaller than the step size, the shift toRayleigh distributed walker separation occurs during theballistic phase of the motion (over times shorter thanabout one step). This is ilustrated by black curve inFigure 1). For that pair-separation example, the initialseparation was taken to be twice the typical simulationwalker radius (the walker interaction distance), the valueimmediately after contact between two walkers, and thestep size was take to be a factor of about 12.5 greaterthan this, a value typical of most the simulations un-dertaken for this paper. The pair separation becomesRayleigh as the variance transitions from ballistic to dif-fusive scaling, and the variance asymptotes directly tothe Rayleigh value. In most of the simulations we con-duct the walker interaction distance is taken to be muchsmaller than the step length, and thus, for simulationtimes longer than about one step after contact, the sepa-ration between walkers is Rayleigh distributed, with themean separation increasing as the square-root of timeand the variance increasing linearly with time. IV. WAITING TIME STATISTICS
We have conducted random walk simulations of a col-lection of discrete walkers on a periodic two-dimensionalplane, as described in Section II above, and have mea-sured the time interval between successive contacts foreach individual. Here we examine the distribution ofthese waiting (or interarrival) times, excluding intervalsof zero time, which indicate that the walker is in continu-ous contact with another walker. The probability densityof the duration of the continuous contact intervals, andits sensitivity to the random walk parameters, is consid-ered in Section V.
A. Contact interval
The normalized waiting-time probability densities p (∆ t ) for one series of such simulations are shown inFigure 3, with the waiting time ∆ t measured in units ofnumber of steps taken. The simulations in this series havedifferent walker number densities n , but the same meanstep-size and interaction distance. The waiting time dis-tributions are vertically shifted, each by a factor of ten forclarity, with uppermost curve plotting the unshifted dis-tribution for the simulation with the lowest walker num-ber density ( n = 4) and the lowermost plot showing thatfor the simulation with the highest ( n = 1600). We notethat for a given walker radius a (one-half the specified in-teraction distance), each of the distributions in Figure 3approximates a family of distributions sharing the samevalue of n δr ( n δr ≈ [0 . , . , . , . , ., . ] in Fig-ure 3, top to bottom ). To achieve the same waiting timedistribution with decreased random-walk step-size, thenumber density must increase so that the mean numberof collisions per step remains constant. The overlappingdistributions plotted with similar color schemes third andfifth from the top in Figure 3 resulted from pairs of walkswith the same product n δr , but for which n and δr werechanged by factors of 2 and 1/2 respectively. Other fac-tors were tested as well, but are not shown. The distri-butions are very similar for any given n δr , though there FIG. 3: Normalized probability densities p (∆ t ) of the timeinterval ∆ t (measured in units of mean step duration) be-tween two successive contacts for any individual randomwalker in a population of walkers. Distributions are off-set vertically by a factors of ten for clarity. Cases differin the ratio R of the mean-free-path to mean step-length,with results from simulations with values of approximately R ≈ [4980 , , , , . , .
45] shown top to bottom .The random walk simulations which produced the overlap-ping distributions in the upper two most plots were conductedin both 1 × × n and δr but differing values for these quan-tities independently. An analytic approximation (see text) forthe large R waiting time distribution is over plotted in black . are small differences that become more apparent in thecount statistics (Section IV B).Separately, sensitivity to domain periodicity was inves-tigated. The two uppermost pairs of waiting time distri-butions in Figure 3 each plot the distributions achievedby two random walk simulations sharing identical pa-rameter values but conducted in different sized domains,3 ×
3, shown in dark blue and magenta , and 1 ×
1, shownin light blue and purple (all other simulations for whichthe waiting time distribution is plotted in Figure 3 sharethe 1 × n = 4). This suggests a possible weak sensitiv-ity to domain periodicity that decreases with increasingdomain size and/or number density because edge effectsbecome less important relative to contacts in the bulk ofthe domain as the number of walkers simulated increases.This sensitivity to domain periodicity for small walkernumbers is again more apparent in the count statistics(Section IV B) than directly in the waiting time distribu-tions shown here.The waiting time distributions plotted in Figure 3change dramatically from top to bottom depending onthe ratio, R , of the mean-free-path (the mean distancethat would be traveled on ballistic trajectories betweencollisions) to the random-walk mean step-length. On atwo-dimensional plane, this ratio R = 1 / (4 n a δr ). Forhigh number densities (small R, short mean-free-pathcompared to the mean step-length), the walkers take fewsteps between collisions. In the limit where one step orless is taken before collision, the walkers move on ballistictrajectories between collisions and the interval betweencollisions is distributed as the mean-free-path collisionrate [e.g., 3], p (∆ t ) = 1 / ¯ τ exp( − ∆ t/ ¯ τ ), where the meanfree collision interval ¯ τ = 1 / (4 n a ¯ v ), with ¯ v = √ v theaverage relative ballistic velocity [e.g., 4] between walk-ers ( v = δr/δt ). This limiting exponential is over-plottedwith solid lines in Figure 3 for the three cases with highestwalker number densities. The waiting time distributionfor random walker collections with small R approach thisasymptote, though small deviations are apparent in thesimulation results. We note that, without good resolu-tion of the shortest waiting times and confidence in theshortest waiting time counts, the contact interval proba-bility densities obtained with larger R (smaller step-sizeto walker separation ratios) could also be misconstrued tobe exponential, with consequent implications for whetherPoisson count statistics are appropriate (Section IV B).In the opposite limit, low number densities (large R,long mean-free-path compared to the mean step-length),the walkers undertake many random walk steps betweencontacts. In this limit, the waiting time is the mini-mum statistic of the first-passage time for any individualwalker in the group to cross into a circle of radius 2 a sur-rounding a given individual. This minimum statistic canbe computed, most simply when that first-passage timedistribution is known analytically. Unfortunately, that isnot the case in two dimensions (see Appendix A 2). Al-ternatively, one can compute the first-passage probabilityof the walker pair separation in one-dimension, the first-passage probability to a pair separation of 2 a startingfrom a separation greater than that.The pair separation distance r s between two unbiasedrandom walkers, as a function of their initial separa-tion r at t = t , is Rice distributed (Equation 2 with σ = 4 D ( t − t )). As in the case of the fully two-dimensional formulation, this distribution does not read-ily yield an analytic expression for the first-passage timeprobability density (neither Equation A4 nor A5 can besolved analytically when free-particle probability densityis Rice distributed). However, as discussed in Section IIIand illustrated by Figure 2 above, Equation 2 has twolimiting forms. For 1 . (cid:46) t − t ≤ r / (2 δr ), the pair separation distribution evolves as a truncated Gaussian, p ( r s , t | r , t ) ≈ √ πσ exp (cid:18) − ( r s − r ) σ (cid:19) , (5)and for this the first-passage time distribution reducesto the one-dimensional solution derived in the appendix(Section A 1, Equation A9 with a → a and D → D ).While convenient, this solution based on the Gaussianseparation-distribution limit is not particularly relevant.It is very unlikely that the waiting-time interval in a col-lection of random walkers will be dominated by the min-imum time to contact for each initial separation.The second limiting form of the separation distributionis more useful. For t − t (cid:38) r /δr , the separationdistribution is approximately Rayleigh, p ( r s , t | r , t ) ≈ r s σ exp (cid:18) − ( r s − r ) σ (cid:19) . (6)This yields, via Equation A4 for example, the Laplacetransform of the first passage distribution˜ f ( s | r ) = K (cid:18)(cid:113) ( r s − r ) s D (cid:19) K (cid:18)(cid:113) ( r s − a ) s D (cid:19) , (7)with the Bromwich integral solution f ( t | r , t ) = 2 π (cid:90) ∞ u e − u ( t − t ) J ( Bu ) Y ( Au ) − J ( Au ) Y ( Bu )[ J ( Bu )] + [ Y ( Bu )] du , (8)where A = | r s − r | / √ D and B = | r s − a | / √ D , and K , J and Y denote the lowest order modified Besselfunction of the second kind and the lowest order Besselfunctions of the first and second kind respectively [2].Beyond this formal solution, the inverse Laplace trans-form of Equation 7 can be determined analytically in thelarge s (short t − t ) limit. In this limit and to low-est order, K ( x ) ≈ √ π/ x − / e − x [2], and the inversetransform of Equation 7 is f ( t | r , t ) ≈ ar ( r − a ) (cid:112) πD ( t − t ) exp (cid:18) − ( r − a ) D ( t − t ) (cid:19) , (9)when taken to be independent of r s , as appropriate. Un-der these approximations, the first contact time t − t be-tween any two walkers is L´evy distributed, and dependson their initial separation r and radii a . We note thecurious fact that, under the approximation that the sep-aration between pairs of random walkers r s is Rayleighdistributed, and in the limit of large s (short t ), the first-crossing time distribution in two dimensions is equivalentto that without approximation in the three-dimensions(Equation A12). FIG. 4: Probability density p ( r ) of the distance betweenpairs of random walkers on a periodic 1 × solid blue ) and that between two randomly chosenpoints on the same plane ( orange dashed ) [7]. Given this first-passage-time probability density(Equation 9), one can calculate the waiting time distri-bution. The waiting time for any individual walker in alarge R collection of walkers is the minimum first-passagetime to a of separation 2 a . It is the shortest time to reacha distance 2 a from any other member of the collectionof discrete walkers. For independent but non-identicallydistributed variates, the minimum oder-statistic is givenby F (1) ( x ) = 1 − (cid:81) Mi (1 − F i ( x )), where (cid:81) indicates theproduct and F i ( x ) are the cumulative distributions fromwhich M samples are drawn [e.g., 5, 6]. For the L´evy dis-tributed first-crossing time (Equation 9), the cumulativedistributions are complementary error functions [2], so F (1) (∆ t ) = 1 − N − (cid:89) i (cid:34) − ar i erfc (cid:32)(cid:114) ( r i − a ) D ∆ t (cid:33)(cid:35) , (10)where ∆ t = t − t is the waiting time interval for anindividual walker and r i is the separation between itand each of the N − F (1) (∆ t ) can be determined by sampling the probabil-ity distribution governing the distance between any tworandomly chosen points on the plane (Figure 4) [7], p ( r ) = πr − r + 2 r , r ≤ − (2 π +4) r + 8 r (cid:113) r − − r + 8 r arcsin (cid:18) r (cid:19) , < r ≤ √ , (11) N − r i . The average first-passage cumulative distributionis achieved by repeating this process many times, and thecorresponding average probability density is plotted as a black curve in Figures 3 and 5. FIG. 5: Probability density p (∆ t/ ¯ τ ) of the waiting time ∆ t scaled by the mean-free-collision interval ¯ τ = 1 / (4 n a ¯ v ), with¯ v = √ v . Same simulations and color scheme as FIG. 3, withthe analytic approximation for the waiting time distributionover plotted in black . Inset replicates plots for two pairs ofrandom walk simulations (upper pair offset vertically), illus-trating dependence on domain periodicity (lower pair) andthe small mismatch between solutions with the same value of R but differing number density n and mean step length δr (upper pair). B. Contact counts
The non-exponential waiting time distributions ob-served imply non-Poisson count statistics [e.g., 8]. Therapid decline in probability density at short waitingtimes yields a rapid decrease in the hazard function, h (∆ t ) = p (∆( t ) / (1 − P (∆ t )), where P (∆ t ) is the cu-mulative waiting time distribution, with increasing wait-ing time. This suggests negative duration dependenceand over-dispersion of the contact count probability massfunction [e.g., 9], p ( N t ), where N t is the number of con-tacts an individual walker has with any other over a timeinterval of t steps. Contacts between random walkers areindependent, but the probability of contact depends onthe on the time since the last occurrence. It is much morelikely for any two random walkers to contact each otheragain if they have recently been in contact because themean and variance of the separation between any pairof walkers increases with time (Section III). This in turnsuggests that, even in the case of small R (short mean-free-path compared to random-walk mean step-length),the waiting time probability density is not exponentialfor very short time intervals; that the waiting time is notexponentially distributed over waiting times for whichthe likelihood of recontact with a previously contactedwalker is comparable to the likelihood of contact with anew walker.This is indeed the case. Figure 5 re-plots Figure 3 afterscaling the waiting time by ¯ τ , where ¯ τ is the mean-free-collision time, the mean time between collisions if thewalkers moved ballistically (without change in direction)at the random walk step velocity δr/δt . For reference,the mean-free-collision time has a values ranging fromabout 12 .
45 random walk steps for the the lowest R casewe computed (that with n = 1600, the bottom most dis-tribution shown in Figure 3) to 4980 for the highest valueof R (that with n = 4, the upper most distribution shownin Figure 3). For the latter, the abscissal range in Fig-ure 5 is thus about twice that shown in Figure 3, while forthe former it is a small fraction. Note that, as ¯ τ → n → ∞ for a given finite interaction distance and step-size), the range of waiting times over which distributionis non-exponential goes to zero; the distribution becomesstrictly exponential for all waiting times. In that limit,all walker contacts occur on strictly ballistic trajectoriesand the negative duration dependence of the waiting timevanishes because and there is no chance of recontactinga previously contacted walker before another.Scaled by the mean-free-collision time, the waitingtime probability densities collapse to a nearly commondistribution (Figure 5). They are all non-exponential,and the differences between them are small, though wenote some differences: The lower pair of purple and ma-genta curves in the inset of Figure 5 illustrate the differ-ences induced by domain periodicity when a very smallnumber total of walkers are simulated (as discussed inSection IV A above). The upper vertically-offeset pink and red curves in that inset, on the other hand, illus-trate the waiting-time distribution differences which re-sult when the simulations share the same R (and thus¯ τ ) but differ in n and δr (holding a constant). Thedifferences are small but indicate that simulations arenot strictly similar. They are sensitive to the changesin the relationships between the interaction distant andthe step-size and between the interaction distant and themean spacing of the walkers in the domain. These sensi-tivities were quite apparent when constructed the com-mon R simulations we have been examining. Simulationswhich shared a common value for a while varying n or δ r and holding the n δr constant resulted in nearly identicalwaiting time distributions, as we have illustrated. How-ever, those sharing common n or δr while holding aδr or aδn constant showed very different waiting time distri-butions. At close range and over short waiting times, thefrequency of repeated walker contacts is very sensitive towalker size. This is further examined in the Section V.The mean-free-collision time scaling is also useful inunderstanding the non-Poisson count statistics observed.The upper panel of Figure 6 displays (for each of thesimulations whose waiting time distribution is shown inFigure 3) the number of contacts an individual walkerexperiences over a 2000 step interval. Unsurprisingly,higher counts occur in simulations with higher walkernumber densities (smaller R ). Perhaps less expected,significant non-Poisson statistics, over-dispersion of thecontact count probability mass function, are apparenteven in the highest number density case, a case for which FIG. 6: Probability mass function p ( N t ) of the number ofcontacts N t experienced by a walker in time interval t . Sim-ulation color coding is the same as used in Figures 3 and 5.In the upper panel, count distributions are plotted for 2000step intervals in each of the simulations. In the lower panel,count distributions are plotted for all simulations for four dif-ferent time intervals scaled by the mean-free-collision time:8, 4, 2, and 0.4 ¯ τ (offset top to bottom ). For reference thesecorrespond to 40000, 20000, 10000, and 2000 random-walksteps in the n = 4 simulation and 100, 50, 25, and 5 steps inthe n = 1600 run. Poison distributions based on the mean N t values are over plotted with solid-dotted curves in bothpanels. the waiting time distribution is exponential at all butthe very shortest waiting times. The short waiting timeexcess results in non-Poison count statistics in all thesimulations. In the lower panel of Figure 6, the countdistributions are plotted for four different time intervalseach scaled by the mean-free-collision time: 8 ¯ τ , 4 ¯ τ , 2 ¯ τ ,and 0.4 ¯ τ (offset top to bottom ). With the notable ex-ception of those simulations in which the domain peri-odicity is important because the total number of walk-ers is low (shown light blue and purple in Figure 6), thePoisson distributed cores of the count distributions forthe same mean-free-collision-frequency-weighted time in-tervals overlap. In the count-distribution wings, how-ever, this is not the case. The distribution wings shownon-Poisson count excesses that persist for count inter-vals spanning many mean-free-collision times and differin amplitude depending on R . Close range recontact be-tween walkers plays a proportionately larger role in the FIG. 7: Waiting-time distributions for simulations thatshare the same step size and number density but differin interaction distance 2 a . Main plot shows distributions( top to bottom , blue to brown ) for simulations with walkerradii a ≈ [0 . , . , . , . , . δr or equivalently about[0 . , . , . , . , . a ≈ [0 . , . δr orabout [0 . , .
5] times the mean nearest-neighbor distance, grey and yellow respectively. small R (high number density) simulations, and they con-sequently show somewhat greater over-dispersion of thecontact-count distribution than do simulations than withlarger R . V. CONTACT DURATIONA. Sensitivity to interaction distance
The non-exponential waiting-time and non-Poissoncount statistics described in the previous section resultfrom the negative duration dependency of the waitingtime interval. The interval between contacts is short-ened when the walkers are in closer proximity becausethe non-ballistic motions of the walkers brings them intorepeated contact. For times-short compared to the meanfree collision interval, the waiting time reflects the highlynonlinear first-passage probability to a random-walk pairseparation of 2 a . It depends in detail on the first passageintersection of two circular objects on a two-dimensionalplane. Unlike n and δr , the interaction distance 2 a doesnot enter as a simple scaling parameter. Simulations withdiffering values of walker radius a and the same value R do not shown the same waiting time distributions. More-over scaling the waiting time by the mean-free-collisiontime ¯ τ does not bring distributions with differing R intonear alignment, as it does for the constant a cases dis-cussed above (Figure 5). The waiting time distribution is very sensitive to thewalker radius, as illustrated by Figure 7. The simu-lations underlying the distributions plotted there werecomputed with identical walker step size and numberdensities but differing interaction distances. The orange-red curve plots the same waiting time distribution as thatof the same color in Figures 3, 5, and 6. The simula-tion responsible employed a walker radius equal to 4% ofthe mean step size and 0.7% the mean nearest-neighbordistance. The walker radii employed by the other sim-ulations for which distributions are shown range fromone-tenth ( blue curve) to seventy times ( yellow curve inthe inset) those values. In the latter case, the interactiondistance is equal to the mean nearest-neighbor separationof the walkers.It is apperent from Figure 7 that, as the walker ra-dius increases, the waiting time probability density in-creasingly deviates from exponential, but not in the self-similar fashion described for the constant a varying n δr cases discussed in Section IV. As the walker radii getlarger, the non-exponential behavior extends to longertimes, even when time is scaled by the mean-free-collisioninterval. It extends well beyond one mean-free-collisioninterval for even modest increase in the walker radii, withextended notable distortion of the waiting-time distribu-tion apparent when the radii exceed about 5 – 10% thestep size, or about 1 – 2 % the mean nearest neighborseparation. For even larger walkers, as the walker in-teraction distance approaches the mean nearest-neighborseparation, multi walker overlap becomes more likely, andthe shape of probability density function changes dramat-ically ( yellow curve Figure 7 inset). Even for very shortwaiting time intervals it no longer resembles the charac-teristic first-passage distribution derived in Section IV A.In the opposite limit, very small walker size, relegates thenon-exponential behavior to very short times. As a → n δr , the time over which the walkers are inclose enough proximity for non-ballistic motions to beimportant for recontact also goes to zero. The mean-free-collision interval then goes to infinity and the wait-ing time distribution becomes strictly exponential for allwaiting times. B. Contact duration
Similarly the contact duration t c is highly sensitive tothe interaction distance or walker radius. The contact-duration probability density is shown in the upper panelof Figure 8 for the same set of simulations as those forwhich the waiting-time distributions are plotted in Fig-ure 7. As expected, the mean contact duration increaseswith increasing walker size (interaction distance), but ad-ditionally the distributions are structured, showing a dis-crete peak. This can be understood as follows:In most of the simulations conducted for this paperthe mean random walk step size is much larger than theinteraction distance which defines contact. For the ex-amples shown in the main body of Figure 7, 2 a/δr ≈ [0 . , . , . , . , . a (first contact) to lessthan 2 a and back (last contact). It is given by the t c > t c = 2 aδr (cos θ − cos θ )1 − cos( θ − θ ) , (12)where θ and θ are the angles the walker velocity vec-tors make with the axis between them at first contactand time is measured in steps so that their velocitiesare equal to δr (since δt = 1). The crossing durationdistribution under this approximation is obtained inde-pendently and uniformly (between zero and 2 π ) sampling θ and θ and repeatedly evaluating Equation 12. Theprobability densities of resulting values are plotted in thelower panel of Figure 8. Each of them peaks at 2 a/δr ,the ballistic contact time for head-on collisions (markedwith vertical dashed fiducial lines in the figure), with thedistribution wing to lower duration times reflecting thechord length distribution for walkers moving in oppositedirections and that to long contact durations reflectingcontacts between walkers traveling in the same direction.The ballistic-crossing-duration distributions shown inthe lower panel of Figure 8 are self-similar, overlappingwhen t c is scaled by 2 a/δr . This is not the case for thecontact durations in the random walk simulations (up-per panel). While the simulation contact-duration dis-tributions share the same peak locations as the ballisticmodels, the shapes of the distribution wings are sensitiveto the walker radii (step size and number density heldconstant). As the interaction distance (walker radius)increases the probability of a directional change duringcontact also increases. Such non-ballistic changes in thedirection of motion flatten the distribution peak and atthe same time lessen the probability of very short or verylong duration contacts, shrinking its tails. Large val-ues of the interaction distance first exceed the step size( grey distribution) and then the mean nearest-neighborseparation ( yellow curve), with consequent distortion ofthe distributions. This ordering is not unique. Here itdepends on the specified parameters (the interaction dis-tance, step size, and number density) of the random walksimulations. In more complex and realistic settings it de-pend on the physical properties of the participants andthe underlying dynamics of their motions. VI. CONCLUSION
Here we examined the waiting-time (interarrival-time),count, and contact-duration statistics in a collection of
FIG. 8: Contact duration (in units of steps) distributions (up-per panel) for the same simulations and using the same colorscheme as those ilustrated by Figure 7. Duration of randomlyoriented ballistic intersections (lower panel) between objectsof the same size as the walkers in the simulations illustratedabove. Vertical fiducial lines indicate the overlap time forhead-on collisions, 2 a/δr . noninteracting random walkers on a two dimensionalplane. We uncovered non-exponential waiting time, non-Poission count, behavior associated with a negative du-ration dependency of the waiting time interval. Repeatcontacts between any two random walkers decreases withdistance between them, leading to enhancement of shortwaiting times and over-dispersion of contact counts; non-ballistic motions of walkers in close proximity shortenthe waiting time to repeat contact enhancing the overallnumber of contacts at short times. Thus contact statis-tics are highly sensitive to the random walk step length(correlation length of the motions), the interaction dis-tance, and the mean separation of the walkers (numberdensity), and their hierarchical relationships.The canonical random-walk induced by elastic colli-sions between ballistic trajectories is a special limitingcase. It displays strictly Poisson statistics because thedirectional change in the random walk is caused by thecollisional contact between particles themselves. In manynatural systems this may not be the case. The ran-dom walk characteristics may be determined separatelyfrom the contacts. Examples range from biological, inwhich the random walk may be determined by behavior,0to physical, in which the motion of a dilute componentis governed by collisions with the dominant component.In such cases the ratio of the mean-free-path length torandom-walk-step length, or step-duration to the mean-free-collision frequency, are critical to the contact statis-tics.One important application of this work may be to con-tagion in human populations. The motion of individualsin an urban setting for example, is, as in our model,largely determined independently from contact events.Individuals cross each others’ paths within an interactiondistance (contagion radius) of each other, often with littleinfluence on their motion, and that interaction distance isoften smaller than the correlation length of the motions.Additionally, the motions as a whole may well be approx-imated as a collection of random walks. We suggests thatunder these circumstances, contacts are unlikely to showexponential waiting-time distributions and correspond-ing Poisson counts. Moreover, the number density ofindividuals varies between populations, and this shouldinfluence the time scale over which non-exponential con-tact statistics are apparent. It would be interesting totest the limits of this hypothesis using high resolutioncell-phone location data. Acknowledgements:
Special thanks to Luke Rast andJacob Rast for stimulating discussions, incisive questions,and key pointers.
Appendix A: First-passage time distribution
For the reader’s convenience this appendix provides abrief review of the unbiased random-walk first-passagetime to a point, circle and sphere in one, two and threedimensions respectively. The discussion is particularlyintended to highlight the difficulties encountered in solv-ing the two-dimensional problem (Section A 2 below).There are two general approaches to the problem. Inone [e.g., 10, 11], the zero drift Smoluchowski equation,the Fokker-Planck or Kolmogorov forward equation gov-erning the evolution of the position probability density p ( x , t ), ∂p∂t = D ∇ p , (A1)is solved given an initial condition p ( x , t ) and subjectto the Dirichlet boundary condition p ( x , t ) = 0 on thefirst-passage boundary of interest x ∈ B . The Dirich-let boundary condition allows determination of the inte-grated zero-reflection flux of p across B as a function oftime. This is the first-passage time probability density, f ( t | x , t ) = − (cid:90) D ∇ p · ˆ n d B , (A2)where ˆ n is the outward directed normal to B . Thisapproach is valid in the Brownian limit, when a large number of unbiased steps of length δx are taken over ashort interval of time δt , so that the diffusion coefficient D ∝ δx /δt remains finite as both δx and δt approachzero.A second approach, works directly with the free parti-cle probability density itself [12–15]. For the unbiasedrandom walk in the Brownian limit, the free particleprobability density, p ( x , t | x , t ), is the solution to Equa-tion A1 without boundary B present, for initial condition p ( x , t ) = δ ( x − x ).For x outside of in the region bounded by the first-passage boundary of interest x / ∈ R at time t , the proba-bility density p ( x / ∈ R , t | x / ∈ R , t ) can be expressed asthe sum of two contributions: random walkers that havenot yet passed the boundary s ( x , t | x , t ) and walkersthat have crossed the boundary for a first time at someearlier time τ and returned to x . Thus, for x / ∈ R and x / ∈ R p ( x , t | x , t ) = s ( x , t | x , t ) + (cid:90) tt f ( τ | x , t ) p ( x , t | x ∈ B , τ ) dτ . (A3)This expression is true even if the particle path makesfurther crossings into and out of R on its way back to x / ∈ R , as accounted for by all paths contributing to p ( x , t | x ∈ B , τ ). Similarly, for x ∈ R and x / ∈ R at time t the particle path must have made a first crossing of B at an earlier time τ , so that p ( x , t | x , t ) = (cid:90) tt f ( τ | x , t ) p ( x , t | x ∈ B , τ ) dτ , (A4)where p ( x , t | x ∈ B , τ ) now captures the probability of ar-riving at x ∈ R after the first crossing of B at time τ .Equation A4 is a Volterra integral equation of the firstkind for the first-passage time distribution f ( t | x , t )given the know function p ( x , t | x ∈ B , τ ) and kernel p ( x , t | x ∈ B , τ ). Similarly, Equation A3 can be convertedinto an integrodifferential equation for the survival prob-ability [15], S ( t | x , t ) = (cid:82) s ( x , t | x , t ) d x , by integrat-ing over all x / ∈ R and noting that f ( τ | x , t ) = − dS/dτ , S ( t | x , t ) = (cid:90) p ( x , t | x , t ) d x + (cid:90) (cid:90) tt dS ( τ ) dτ p ( x , t | x ∈ B , τ ) dτ d x . (A5)The time derivative of − S yields the first passage timedistribution, since s ( x , t | x , t ) was defined as the prob-ability density of walkers having not yet crossed theboundary B .
1. First-passage in one and three dimensions
Consider a one-dimensional random walk with a firstpassage boundary point at x = a and initial condition1 p ( x, t ) = δ ( x − x ) with x > a . With a change of vari-ables x (cid:48) = x − a and t (cid:48) = t − t the initial conditionbecomes p ( x (cid:48) ,
0) = δ ( x (cid:48) − ( x − a )), with the first passageboundary now x (cid:48) = 0. Solving the Laplace transforma-tion of Equation A1 with this initial condition, D d ˜ pdx (cid:48) − s ˜ p = − δ ( x (cid:48) − ( x − a )) , (A6)subject to the boundary condition ˜ p (0 , s ) = 0 yields theGreen’s function solution [e.g., 11]˜ p ( x (cid:48) , s ) = 1 √ sD (cid:104) e ± √ sD ( x (cid:48) − ( x − a )) − e − √ sD ( x (cid:48) +( x − a )) (cid:105) , (A7)where the positive and negative exponents in the firstexponential term correspond to solutions for x (cid:48) ≤ x − a and x (cid:48) ≥ x − a respectively. Inverse Laplace transformyields, for all x (cid:48) , p ( x (cid:48) , t (cid:48) ) = 1 √ πDt (cid:48) (cid:20) e − ( x (cid:48)− ( x − a ))24 Dt (cid:48) − e − ( x (cid:48) +( x − a ))24 Dt (cid:48) (cid:21) , (A8)and consequently by Equation A2 a L´evy distributedfirst-passage time, f ( t | x , t ) = ( x − a ) (cid:112) πD ( t − t ) e − ( x − a )24 D ( t − t . (A9)The same result can be reached by the second routeoutlined above [15], using either Equation A4 or A5.For example, Equation A4 in one-dimension, with firstpassage boundary at x = a , initial condition p ( x, t ) = δ ( x − x ), and p ( x, t | x , t ) taken to be the solution tothe one-dimensional diffusion equation, is1 √ πDt (cid:48) e − ( x − x Dt (cid:48) = (cid:90) t (cid:48) f ( τ | x , (cid:112) πD ( t (cid:48) − τ ) e − ( x − a )24 D ( t (cid:48)− τ ) dτ , (A10)where t (cid:48) = t − t . Its Laplace transform (with x > a ) is˜ f ( s | x ,
0) = e − √ sD ( x − a ) , (A11)which when inverse transformed directly yields the L´evydistribution (Equation A9) as previously. It is impor-tant to note, that while in this example of an unbiasedrandom walk, p ( x, t | x , t ) is taken as the free-particlesolution to Equation A1 in the Brownian limit, and thusthe two solution methods yield the same first passagetime distribution, the second method is more broadlyapplicable and can be applied to situations in which thefree-particle probability density function is known but isnot necessarily the solution to Equation A1.In three-dimensions, the unbiased random walk first-passage time distribution to a sphere of radius r = a centered on the origin, from a point outside the sphere,can be similarly derived by solving Equation A1 assumingspherical symmetry, with the initial condition p ( r, t ) = δ ( r − r ) / (4 πr ) and boundary condition p ( a, t ) = 0. Thisis facilitated by noting that the spherically symmetricform of the diffusion equation can be converted to the one-dimensional Cartesian form, by a simple change ofvariables [e.g., 10, 11] u ( r, t ) = r p ( r, t ). The first-passagetime in three dimensions is again L´evy distributed, f ( t | r , t ) = ( r − a ) (cid:112) πD ( t − t ) ar e − ( r − a )24 D ( t − t , (A12)with the minor difference between Equations A9and A12, the presence of the radius ratio a/r in thelater, reflecting the spherical geometry.
2. First-passage in two dimensions
In two dimensions there is no such convenient changeof variables. Consider the unbiased random walk first-passage time distribution to a circle of radius r = a cen-tered on the origin from a point outside of the circle,initially at r = r . The free particle probability den-sity distribution, in the Brownian limit, is the solutionto Equation A1 in axisymmetric polar coordinates withinitial condition p ( r, t ) = δ ( r − r ) / (2 πr ). It is givenby the Rician distribution [e.g., 16–18] p ( r, t | r , t ) = r πr σ exp (cid:18) − r + r σ (cid:19) I (cid:16) r rσ (cid:17) , (A13)with σ = 2 D ( t − t ). Neither Equation A4 nor A5 canreadily be solved analytically for the first-passage timeusing this distribution.Similarly, solution of Equation A1 with boundary con-dition p ( a, t ) = 0 and initial condition p ( r, t ) = δ ( r − r ) / (2 πr ) poses difficulties. With t = 0, the Laplacetransform of Equation A1 in cylindrical coordinates, withthe initial condition above, can be written as a modifiedBessel equation of order zero [2], with solution for r < r ˜ p ( r, s | r ,
0) = 12 πr s K ( √ s/D r ) K ( √ s/D a ) K ( √ s/D a ) I ( √ s/D r ) − I ( √ s/D a ) K ( √ s/D r ) I ( √ s/D r ) K ( √ s/D r ) + K ( √ s/D r ) I ( √ s/D r ) . (A14)Evaluating the Bromwich integral for the inverse of theLaplace transform by contour integration on the complexplane, with a branch cut along the negative real axis, andevaluating Equation A2 using that result, yields the firstpassage time distribution f ( t | r ) = √ Dπ (cid:90) ∞ e − u t J ( r / √ D u ) J ( a/ √ D u ) + Y ( r / √ D u ) Y ( a/ √ D u )[ J ( a/ √ D u )] + [ Y ( a/ √ D u )] du , (A15)where J and Y are the lowest order Bessel functions ofthe first and second kind [2]. Unfortunately, this integralis not readily solvable analytically.2 [1] M. P. Rast and J.-F. Pinton, Phys. Rev. E , 046314(2009).[2] M. Abramowitz and I. A. Stegun, Handbook of mathemat-ical functions : With formulas, graphs, and mathematicaltables (US Department of Commerce, National Bureau ofStandards, Applied Mathematical Series, 55, 1972).[3] S. Chapman, T. Cowling, D. Burnett, and C. Cercig-nani,
The Mathematical Theory of Non-uniform Gases:An Account of the Kinetic Theory of Viscosity, ThermalConduction and Diffusion in Gases , Cambridge Mathe-matical Library (Cambridge University Press, 1990).[4] S. T. Paik, American Journal of Physics , 602 (2014),1309.1197.[5] G. Cao and M. West, Communications in Statistics -Theory and Methods , 755 (1997).[6] H. A. David and H. N. Nagaraja, Order statistics (JohnWiley, 2003).[7] J. Philip, TRITA MAT (2007), URL https://people.kth.se/~johanph/habc.pdf .[8] M. Taboga, Lectures on Probability Theory and Mathe-matical Statistics (CreateSpace Independent Publishing Platform, 2012).[9] R. Winkelmann, Journal of Business & Economic Statis-tics , 467 (1995).[10] S. Chandrasekhar, Reviews of Modern Physics , 1(1943).[11] S. Redner, A Guide to First-Passage Processes (Cam-bridge University Press, Cambridge, 2007).[12] E. Sch¨odinger, Physikalische Zeitschrift , 289 (1915).[13] A. J. Siegert, Physical Review , 617 (1951).[14] O. Lumpkin, Phys. Rev. A , 2758 (1995).[15] Z. Hu, L. Cheng, and B. J. Berne, J. Chem. Phys. ,034105 (2010).[16] S. O. Rice, Bell System Technical Journal , 46 (1945).[17] K. D. Cole, J. V. Beck, A. Haji-Sheikh, and B. Litk-ouhi, Heat Conduction Using Green’s Functions , Seriesin Computational Methods and Physical Processes inMechanics and Thermal Sciences (CRC Press, Taylor &Francis Group, 2011).[18] P. Agrawal, M. P. Rast, M. Goˇsi´c, L. R. Bellot Rubio, andM. Rempel, Astrophys. J.854