aa r X i v : . [ c ond - m a t . s t a t - m ec h ] J u l Vicious L´evy flights
Igor Goncharenko and Ajay Gopinathan
School of Natural Sciences, University of California, Merced, California, 95343, USA
We study the statistics of encounters of L´evy flights by introducing the concept of vicious L´evyflights - distinct groups of walkers performing independent L´evy flights with the process terminatingupon the first encounter between walkers of different groups. We show that the probability thatthe process survives up to time t decays as t − α at late times. We compute α up to the secondorder in ε -expansion, where ε = σ − d , σ is the L´evy exponent and d is the spatial dimension. For d = σ , we find the exponent of the logarithmic decay exactly. Theoretical values of the exponentsare confirmed by numerical simulations. PACS numbers: 64.60.ae, 64.60.F-, 05.40.Jc, 64.60.Ht
Diffusive processes with long range jumps play an important role in many physical, chemical and biological phenom-ena. A L´evy flight is an example of such a process where the probability distribution of the length of an individualstep, r , is governed by the power-law r − d − σ , where d is the dimension of the space and σ is the L´evy exponent.Smaller values of σ therefore produce longer range jumps while for σ ≥
2, the mean jump length is finite and simplediffusive behavior is recovered. L´evy flights have been used to describe a wide range of processes including epidemicspreading, transcription factor proteins binding to DNA, kinetic Ising models with long-range interactions, foraginganimals and light propagation in disordered optical materials [1–6]. While individual L´evy flights have been studiedin great detail, the same is not true if we consider several distinct groups of L´evy flights. One could, for example, beinterested in the statistics of encounters between members of different groups. This question is relevant for processeswhere the outcome depends on the occurrence of such encounters. Examples include sharks and other marine animalssearching for prey [7], chemical reactions in turbulent environments [8], electron-hole recombination in disorderedmedia [9] and even male spider-monkeys encountering their mates or other aggressive males in the forest [10].In this Letter we compute the survival probability, i.e. the probability that no two members of different groupsof L´evy flights have met up to time t . For the case of simple diffusion with exactly one particle in each group, thiscorresponds to the classic problem of Gaussian vicious walks [11], i.e. walks that are prohibited from being on thesame site at the same time, but remain independent otherwise. Here we generalize this concept to groups of L´evyflights under the same constraints. We term them vicious L´evy flights (VLF). We consider p sets of particles with n i particles in each set, i = 1 . . . p , that are driven by L´evy noise on the d dimensional regular lattice. A pairwiseinterset short-range (delta-function) interaction is introduced to guarantee that trajectories which continue beyondan intersection are discarded, i.e. have zero statistical weight. This terminates the process at the first encounterbetween members of different groups. Particles belonging to the same set do not interact. We note that L´evy flightsare allowed to jump over each other, unlike ordinary random walks which can only jump to neighboring sites and cannot intersect with the vicious constraint. In d = 1 this means that the ordering is preserved for vicious walks butnot for VLF. For simplicity we assume that L´evy exponents for all flights are the same. Generalization to the case ofdifferent L´evy exponents will be done elsewhere. At time t = 0 all particles start in the vicinity of the origin. We areinterested in the survival probability of this system at late times.We start with a field theoretic formulation of the problem. Methods to formulate field theories for such stochasticsystems are well established [12, 13]. Specifically for Gaussian vicious walks, such a formulation exists [14] and theform of the action is known. We can adapt the action to our case by replacing the Laplacian ∇ with the operator ∇ σ that generates long-range jumps. This gives S ( φ i , φ † i ) = Z dtd d x p X i =1 [ φ † i ∂ t φ i + φ † i ∇ σ φ i ] + X ≤ i 2. VLF exhibit different phases (see inset Figure 1) depending on the values of d and σ . In the mean field phase for d > d c (Region I), the survival probability of VLF is non-zero at infinite timebecause the walks become non-recurrent and particles can avoid each other for all time. For σ ≥ FIG. 1: a. ln( S ) vs ln( t ) for two identical VLF in d = 1. σ values from top tp bottom are 1 . , . , . , . , . σ − d plane. reproduce Gaussian vicious walks. For d < σ < ε -expansion ( ε = σ − d ) around meanfield theory in two-loop approximation.We now turn to the renormalization group analysis. The propagator given by (1) is Γ (1 , j ( s, k ) = ( s + k σ ) − . Theparticular form of the vertex in (1) leads to the fact that there are no diagrams which dress the propagator. Thisimplies that the bare propagator is the exact propagator for the theory. The proper vertex is defined by factoringout external legs from the ordinary 4-point Green’s function of (1) G (2 , ij ( s l , k l ; λ ). Here λ = { λ ij } is the collection ofcoupling constants and ( s l , k l ) for l = 1 . . . (2 , ij ( s l , k l ; λ ) = G (2 , ij ( s l , k l ; λ ) Q m =1 Γ (1 , ( s m , k m ) . (2)The renormalized coupling, λ Rij , is the value of the proper vertex at ( s l = µ, k l = 0) , for all l , where µ is arenormalization group flow scaling parameter. It is possible to sum all diagrams in the series with the result λ Rij = λ ij (1 + λ ij I ) − , (3)where I = (2 π ) − d R d d q (2 µ + 2 q σ ) − . The value of the integral I is given by I = Aµ − ε/σ ε − , where A = ε Γ( d/σ )Γ( ε/σ )2 d π d/ σ Γ( d/ 2) = 2 − σ π − σ/ Γ( σ/ 2) + O ( ε ) , (4)is the geometric factor at the leading order in ε = σ − d . By introducing the redefined coupling constant g Rij = λ Rij µ − ε/σ , we obtain the following renormalization group flow equations (see Appendix A for details): µ ∂g Rij ∂µ = ( − ε + Ag Rij ) g Rij /σ. (5)The fixed point of this flow is g ∗ = εA − . We note that this value of the fixed point is exact to all orderssince all diagrams were taking into account in (3). The stability of the fixed point follows from the fact that − ∂β ij /∂g Rij | g Rij = ε/A = − ε/σ < β ij = µ∂g Rij /∂µ is renormalization group beta function.We now consider the survival probability which is defined as the correlation function [14] S ( t ; λ ) = Z p Y i =1 n i Y α i =1 d d x i,α i h φ i ( t, x i,α i )( φ † i (0 , n i i , (6)with the measure R D φ † D φ exp[ −S ]. The Feynman diagram of (6) at zero order is a vertex with 2 N external legs.Similar to the case of the (2 , s l , k l ; λ ) = FIG. 2: a. 1-loop diagram corresponding to the integral I b., c. 2-loop intgrals corresponding to I and I respectively. S ( s l , k l ; λ ) / Q Nm =1 Γ (1 , ( s m , k m ). The finite renormalized truncated correlation function Γ R ( s l , k l ; λ R , µ ), where λ R = { λ Rij } is a collection of renormalized coupling constants, is related to the bare truncated correlation function byΓ R ( s l , k l ; λ R , µ ) = Z ( λ, µ )Γ( s l , k l ; λ ), where Z ( λ, µ ) is the scaling function. From this one obtains the renormalizationgroup equation for Γ R ( s l , k l ; λ R , µ ) using the chain rule( µ ∂∂µ + β ij ∂∂g Rij + γ )Γ R ( s l , k l ; λ R , µ ) = 0 , (7)where γ = µ∂ ln Z/∂µ. At the fixed point (7) reduces to ( ∂/∂ ln( µ ) + γ ∗ )Γ R ( s, µ ) = 0 whose solution isΓ R ∼ exp( Z µ γ ∗ d (ln( µ ′ ))) , (8)where γ ∗ = µ∂ ln Z/∂µ | g Rij = g ∗ . Since γ ∗ is constant at the fixed point we have Γ R ∼ µ − γ ∗ . The fact that thedimensions of field and the action are [ φ † ] = [ φ ] = k d/ and [ S ] = 1 implies that [ S ] = 1. Thus it follows that thesurvival probability can only be a function of the dimensionless product µt . From this one infers that the asymptoticbehavior of the survival probability is S ( t ) ∼ t − γ ∗ which gives α = γ ∗ . In order to find Z one uses a normalizationcondition on Γ R that fixes the value of Z . This can be chosen as Γ R ( s l , k l ; λ R , µ ) = 1 when s l = µ and k l = 0 for all l . This implies that Z = Γ( µ, λ ) − . Γ( µ, λ ) can be expressed as a series, up to two-loop order, of the integralscorresponding to the diagrams shown in Fig.2 with appropriate combinatorial factors originating from the numberof distinct ways in which propagators can be assigned to the same diagram. I was evaluated before. The integralcorresponding to the diagram 2c is I = Z d d kd d q (2 µ + 2 k σ )(3 µ + q σ + k σ + | k + q | σ ) (9)This can be evaluated using Mellin-Barnes representation [18] which replaces the sum in the denominator of | k + q | σ and q σ by the product of these terms raised to some power. The result for I up to the leading order reads I = 2 − σ π − σ Γ( σ/ µ − ε/σ (cid:18) ε + 2( − log(3 / / C ) σε (cid:19) , (10)where C = [ ψ (0) ( σ/ 2) + log(4 π )] / ψ (0) ( x ) is standard digamma function. The details of the calculations aresummarized in Appendix B. Knowing I , I and the appropriate combinatorial factors allows us to evaluate Γ( µ, λ )and therefore Z ( µ, λ ) as a series. Differentiating ln Z with respect to µ and substituting λ with λ R (inverting (3))and then taking the value at the fixed point gives us the survival probability exponent α = γ ∗ , with the final result(see Appendix C): α = X ≤ i 2. For σ < σ ≥ Now we describe the details of the numerical simulation that we used to confirm our results in d = 1. At t = 0 westart with N = P pi =1 n i particles belonging to p distinct sets placed equidistantly on the lattice. At each time step wegenerate N random variables, x j , drawn from the uniform distribution on the interval (0 , l j = x − /σj with equal probability to the left or to the right. This procedure generates an independent L´evyflight trajectory for each particle. The process stops whenever particles from different sets land on the same site. Weperform ∼ iterations for each set of parameter values. The survival probability S ( t ) is defined by the numberof processes that survived beyond time t divided by the total number of iterations. Figure 1 shows the plot of thesurvival probability as a function of time for N = 2 for different values of σ . It is clear that at late times ln S ( t ) islinear in ln t verifying our predicted power-law decay S ( t ) ∼ t − α . The critical exponent α is evaluated from the slopeof the best fit line of the late time data.We first consider systems with exactly one particle in each set. Figure 3 shows the value of exponent α for variousvalues of σ with the total number of particles N = 2 , d = 1. Values of σ ≥ N = 2 case is equivalent to the first returnprobability of a single L´evy flight to the origin after time t which scales as t − d/σ [19] and matches our results. For σ < 1, or d > d c , we expect mean field behavior where survival probability at late times approaches a non-zero valueimplying that there is a finite probability that L´evy flights with σ < α for σ = 0 . 9. For 1 < σ < σ close to one, the 1-loop result is in good agreement with the simulation. For larger values of σ , the 2-loopcorrections perform better (see Fig.3 inset). It is to be noted that the discrepancy between theory and simulationbecomes large for higher values of N because the combinatorial factors in (11) become large and we therefore need tokeep higher order terms in ε for the same degree of accuracy. It is interesting that the 1-loop approximation worksreasonably well over the entire range of 1 < σ ≤ σ = 2. We notice that in all cases the value of the survival probability exponent α increaseswith σ starting from α ∼ σ < σ − α FIG. 4: α vs σ for predator and prey problem in d = 1. Symbols represent simulation data for 4 predators and 1 prey (circles)and 3 predators and 2 prey (squares). Lines are 2-loop approximation from formula (11). decays faster for smaller values of σ [16].We now consider a system that consists of 2 sets of identical VLF with different numbers of independent particles ineach set. We shall call one set predators and the other set prey. Figure 4 compares the values of the 2-loop exponentsto the simulation results for various values of σ for two different cases: 4 predators - 1 prey and 3 predators - 2prey. Similar to the previous case we have mean field and Gaussian behaviors for σ < σ ≥ σ and total number of predators plus prey, the number of potentially lethal encounters is maximizedwhen the difference between the number of predators and prey is the smallest implying that the survival probabilitywill decay faster as seen in Fig.4.Our results suggest that it is interesting to solve the problem in the general case of particles with different diffusionconstants and L´evy exponents. The predictive power of the ε -expansion for VLF, that we have demonstrated, shouldbe useful in many applications of practical importance. Examples include the optimization of the predator-preysearch [21] or trapping probabilities [22]. Generalization to the case of intelligent predators, i.e. interacting with aprey by means of the long-range potential, may lead to different critical behavior [23–25]. Simple diffusion processesin power-law small world networks are effectively L´evy flights [26] with the exponent σ controlling the distributionof long-range links. Our work could be used to understand what network structure, or what σ , would optimize thesearch and how much more efficient several independent searchers will be.The authors would like to acknowledge UC Merced start-up funds and a James S. McDonnell Foundation Awardfor Studying Complex Systems. APPENDIX A Here we derive the 1-loop integral I = I ( µ ) = (2 π ) − d Z d d k (2 µ + 2 k σ ) − . (14)We will use dimensional regularization. First we notice that there is no angle dependence under the integral thus onecan integrate out d − X − λ = Γ( λ ) − Z + ∞ dαα λ − exp( − αX ) (15)to handle 1d momenta integral: I = S d + ∞ Z k d − dkµ + k σ = S d + ∞ Z ∞ Z dαdkk d − exp( − αµ − αk σ ) = S d σ Γ( d/σ ) + ∞ Z dαα − d/σ exp( − αµ ) , (16)where S d = 2 π d/ / Γ( d/ 2) is the area of the d-dimensional unit sphere. After taking the integral over α one has I = S d σ Γ( d/σ )Γ( ε/σ ) µ − ε/σ = Aµ − ε/σ ε − + O ( ε ) , (17)where A has been defined by the formula (4).Now we show details of deriving renormalization group flow equations (5). We start with equation (3) and express λ ij in terms of λ Rij . The result reads λ ij = λ Rij − λ Rij I . (18)Multiplying left and right hand side of the last equation on µ − ε/σ and redefining the coupling constant g Rij = µ − ε/σ λ Rij we infer that µ − ε/σ λ ij = g Rij / (1 − g Rij Aε − ) . (19)Differentiating left and right hand side of (19) with µ ∂∂µ we obtain( − ε/σ ) µ − ε/σ λ ij = − β ij g − Rij / ( g − Rij − Aε − ) (20)Now we substitute (19) into (20) and find beta function up to second order in small ε and g Rij expansion: β ij = ( − εg Rij + Ag Rij ) /σ + O ( εg ) (21) APPENDIX B Here we derive the 2-loop integral I = I ( µ ) = (2 π ) − d Z d d kd d q [(2 µ + 2 k σ )(3 µ + k σ + q σ + | k + q | σ )] − . (22)The term | k + q | σ leads to the appearance of angle integration. Nevertheless it is possible to avoid angle integration.The key idea is to use Mellin-Barnes representation [18]:1( X + Y ) λ = Z + i ∞− i ∞ dz πi Y z X λ + z Γ( λ + z )Γ( − z )Γ( λ ) (23)Applying MB formula twice we split the sum of two terms containing q integration into the factor of these termsraised to some power: I = Z d d kd d q π ) d + i ∞ Z − i ∞ dz πi Γ(1 + z )Γ( − z ) µ + k σ (3 µ + k σ + q σ ) z | k + q | σ (1+ z ) = Z d d kd d q π ) d + i ∞ Z − i ∞ dz dz (2 πi ) Γ(1 + z )Γ( − z + z )Γ( − z ) µ + k σ (3 µ + k σ ) z | k + q | σ (1+ z ) q σ ( − z + z ) , (24)Now integral over q becomes standard: I q = Z d d q ( q ) a (( k + q ) ) a = π d/ k d − a + a ) Γ( a + a − d/ d/ − a )Γ( d/ − a )Γ( a )Γ( a )Γ( d − a − a ) , (25)where a = σ ( − z + z ) / a = σ (1 + z ) / 2. Thus we will be left with integral over k of the form: I k = Z d d kk − ε − σz (3 µ + k σ )2 µ + 2 k σ (26)The function under the integral does not depend on the angle and therefore I k can be cast into one dimensionalintegral over momenta: I k = S d σ + ∞ Z dkk − ε/σ − z (3 µ + k ) z µ + k (27)We will compute this integral using alpha representation. I k = S d σ + ∞ Z dkdα dα α − z − k − ε/σ − z Γ( − z ) exp( − µα − α k − α µ − α k ) (28)After momenta integration we obtain I k = S d Γ(1 − ε/σ − z )2 σ Γ( − z ) + ∞ Z dα dα ( α + α ) ε/σ + z − α − z − exp( − µα − α µ ) (29)First we will take care the integral over α . We do substitution ˜ α = α + α I k = S d Γ(1 − ε/σ − z )2 σ Γ( − z ) + ∞ Z dα α − z − e − µα + ∞ Z α d ˜ α ˜ α ε/σ + z − e − ˜ α µ = S d Γ(1 − ε/σ − z ) µ − z − ε/σ σ Γ( − z ) + ∞ Z dα α − z − e − µα Γ(2 ε/σ + z , α µ ) , (30)where Γ( λ, x ) is incomplete gamma function. The value of the last integral can be found in [27]. The final result or I k reads I k = S d σ Γ(1 − ε/σ − z )Γ(1 − z ) Γ(2 ε/z ) µ − ε/σ − ε/σ F (1 , ε/σ, − z , / 3) (31)Inserting (25) and (31) into (24) we infer I = S d π d/ σ (2 π ) d Γ( σ/ − ε/ − ε/σ − µ − ε/σ − ε/σ Γ(2 ε/σ ) + i ∞ Z − i ∞ dz dz (2 πi ) F (1 , ε/σ, − z , / 3) Γ(1 − ε/σ − z ) − z Γ(1 + z )Γ( σ (1 + z ) / 2) Γ( − z + z )Γ( σ ( − z + z ) / 2) Γ( ε/ σz / σ/ − ε − σz / 2) Γ( − ε/ − σz / − σ ( − z + z ) / − ε/ σ/ . (32)First we sum over all poles of Γ( − ε/ − σz / 2) and then over pole at z = 0. The result reads I = S d π d/ σ (2 π ) d µ − ε/σ − ε/σ Γ(2 ε/σ ) F (1 , ε/σ, , / 3) Γ(1 − ε/σ )Γ( σ/ − ε ) Γ( ε/ + ∞ X n =0 ( − n n ! Γ(1 − ε/σ + 2 n/σ )Γ( σ (1 − ε/σ + 2 n/σ ) / 2) Γ( ε/σ − n/σ )Γ( σ ( ε/σ − n/σ ) / 2) Γ( − ε + nσ/ . (33)We will look the final result in the form I = µ − ε/σ e − Bε ( c − ε − + c − ε − ) (34)To obtain the divergent part of I it is convenient to use MATHEMATICA . The result for coefficients c − and c − aregiven by the formula (10). APPENDIX C