Asymmetric Lévy Flights Are More Efficient in Random Search

Abstract

We study the first-arrival (first-hitting) dynamics and efficiency of a one-dimensional random search model performing asymmetric Lévy flights by leveraging the Fokker–Planck equation with a $\delta$-sink and an asymmetric space-fractional derivative operator with stable index $\alpha$ and asymmetry (skewness) parameter $\beta$. We find exact analytical results for the probability density of first-arrival times and the search efficiency, and we analyse their behaviour within the limits of short and long times. We find that when the starting point of the searcher is to the right of the target, random search by Brownian motion is more efficient than Lévy flights with $\beta \le 0$ (with a rightward bias) for short initial distances, while for $\beta >0$ (with a leftward bias) Lévy flights with $\alpha \to 1$ are more efficient. When increasing the initial distance of the searcher to the target, Lévy flight search (except for $\alpha =1$ with $\beta =0$) is more efficient than the Brownian search. Moreover, the asymmetry in jumps leads to essentially higher efficiency of the Lévy search compared to symmetric Lévy flights at both short and long distances, and the effect is more pronounced for stable indices $\alpha$ close to unity.

1. Introduction

Search strategies, in general, and random search processes, in particular, have recently attracted interest among scientists and practitioners who seek to understand food search strategies of animals [1,2] and diffusion control of molecular processes in biological cells [3], or aim at improving the performance of computer search algorithms [4,5] and optimise search processes for military tasks (hunting for submarines, locating enemy vessels and mines) [6,7]. In the concept of random search, the searcher has no prior knowledge about the location of its target and performs a random walk until it encounters the target—this process is called saltatory motion or blind search in the ecology literature [8]. In foraging theory, this saltatory search is distinguished from cruise search, in which the searcher looks out for targets (i.e., resource patches, mates, etc.) during its movements [9,10]. An essential part of all these processes and applications are various forms of random walk models, continuous or discrete in time.

The first studies on random search considered Brownian motion of the searcher as a default strategy. A prominent early example, motivated by his contact with the malaria researcher Ross, is Pearson’s idea of a drunken man performing a “random walk” on the street grid of a city as a model for the spreading of malaria by mosquitoes in previously non-infested areas [11]. Later, Shlesinger and Klafter proposed that Lévy motions with a scale-free, power-law distribution of jump lengths, represent a much more efficient strategy than the standard random walk or (in the continuum limit) Brownian motion, in the search for sufficiently sparse targets. Indeed, this higher efficiency is due to the combination of the thorough local search with occasional long excursions and hence the exploration of previously unvisited areas, thereby reducing the tendency towards “oversampling” of Brownian motion in one and two spatial dimensions [12]. The works on Lévy flight search of albatrosses [13], and on the general optimality of Lévy search [1], have attracted broad attention of researchers in the field of animal motility patterns and optimal search strategies. Viswanathan and colleagues [1,14], in particular, proposed superdiffusive Lévy flights as an optimal search strategy when searching for sparse targets in their “Lévy flight hypothesis”.

For the random search processes based on the Lévy flight model, the distribution of displacement lengths $\left|x\right|$ has the asymptotic power-law form $\lambda \left(x\right)\simeq {\left|x\right|}^{-1-\alpha }$ with $0<\alpha <2$ for which the second moment of the jump lengths diverges, i.e., $\langle x^2\left(t\right)\rangle \to \infty$ [15]. The resulting motion is spatially scale-free, and fractional, qth-order moments $\langle |x{|}^q\left(t\right)\rangle$ with $0<q<\alpha$ exist. Lévy flights have found various applications, ranging from the famous flight of the albatross [13], spreading of spider-monkeys [16], grazing patterns of bacteria [17], over economical data [18,19] to molecular collisions [20], plasmas [21,22,23] and lattice gas automata [24,25]. In recent years, the search processes that optimise random blind searches for sparse targets in different settings have been discussed [9,26,27]. In particular, combinations of Brownian search and Lévy flights were studied [28].

In the last decade, other random search processes have been also proposed as alternatives to Lévy flights and have attracted considerable attention in the scientific community. Such models include intermittent dynamics that switch between local Brownian search events and ballistic relocation phases [29,30], for which the relocation time probability density may also be of the power-law form [27,31]. The intermittent strategy minimises the search time [32], and, for instance, is observed for the dinoflagellate Oxyrrhis marina when it preys on a microzooplankton [33], as well as for Nahua mushroom gatherers [34] following the so-called area-restricted search strategy [35]. In computer science and data mining (information retrieval) [36] one particularly popular search method is based on a random walk model, in which a combination of local and long-range searches is implemented, while in complex network theory, the calculation of the mean first-passage time of a searcher to find the target on a given complex networks is of high importance [37,38] for the description of a large number of real-world phenomena. Lévy flights have been promoted as a preferred strategy when there is insufficient prior knowledge on the search space. In particular, for sufficiently sparse targets, several analyses claim that the optimal value for the power-law exponent is $\alpha =1$ [13,27,39,40], so that the jump lengths are distributed according to the Cauchy law. The central advantage of the Lévy flight’s strategy in comparison with the other strategies is its robustness: while other models work best when their parameters are optimised for specific environmental conditions such as the target density, Lévy flights remain close to optimal even when these conditions are altered [27]. Random search processes based on Lévy flights in the presence of an external drift (underwater current, atmospheric wind, a preference of the walker owing to prior experience or a general bias in an abstract search space), have also attracted attention in recent years [28].

From the mathematical point of view, the problem of calculating the first-hitting or first-arrival time of a point or an interval in the class of stable processes was the subject of several studies. We mention studies of the first-hitting time density of symmetric stable processes [41,42,43,44,45,46,47,48,49] and spectrally positive stable process [50,51] in which the authors obtained series expressions for the first-hitting density. Moreover, for asymmetric stable processes, with the help of the Lamperti–Kiu representation of self-similar Markov processes, the Mellin transform of the first-hitting time and a series expression for the first-hitting time density was obtained [52]. With the use of space-fractional Fokker–Planck equations, the first-hitting problems of symmetric Lévy flights and combined Lévy-Brownian motions were studied [9,28,53,54,55,56,57]. This method is also employed in the present paper to study the first-hitting dynamics of asymmetric Lévy flights. Another method based on the Langevin equation, known as the simulated annealing algorithm, which is a semi-local search consisting of occasional long jumps based on the Cauchy distribution and Lévy stable laws, was studied in [58,59,60], respectively. Moreover, the first-passage problem of symmetric and asymmetric Lévy flights was studied by leveraging the Skorokhod theorem, and the analytical results were confirmed by extensive numerical simulations based on the numerical solutions of space-fractional Fokker–Planck equation and stochastic Langevin equations [61,62,63,64]. We finally note studies of the related problem of barrier crossing of Lévy flights [65,66,67,68,69,70,71].

By providing analytical and numerical investigations of random search based on asymmetric Lévy flights, we here aim at understanding this search strategy by analysis of the first-arrival time, the reliability and the efficiency. We seek to clarify under which conditions this type of asymmetric motion optimises random search. This will be achieved by investigation of both the (deterministic) Fokker–Planck equation with a $\delta$-sink and the associated (stochastic) Langevin equation. The structure of the paper is as follows. In Section 2, we define the process of asymmetric Lévy flight search and compute the first-arrival density in terms of Fox H-functions. We provide a detailed analysis of the asymptotic behaviour of the first-arrival density in the short and long time limits. In Section 3, the analytical results are confirmed by numerical use of the Langevin equation for asymmetric Lévy flights. We also obtain the search efficiency and compare it with the corresponding efficiencies for Brownian and Lévy search in Section 4. The Summary is provided in Section 5. Details of the definition of $\alpha$-stable processes, the definition and properties of the Fox H-function, as well as technical details of derivations are presented in the Appendix A, Appendix B, Appendix C, Appendix D, Appendix E, Appendix F and Appendix G.

2. Formulation of the Problem and Solution

We consider a space-fractional Fokker–Planck equation of the form:

$$\frac{\partial }{\partial t}f(x,t)=K_{\alpha }{}_x{\mathrm{D}}_{\alpha }^{\beta }f(x,t)-{\wp }_{\mathrm{fa}}\left(t\right)\delta \left(x\right),f(x,0)=\delta (x-x_0),$$

(1)

where $x_0>0$ is the initial point, $K_{\alpha }$ is the generalised diffusion coefficient and ${}_x{\mathrm{D}}_{\alpha }^{\beta }$ for $0<\alpha \le 2$ is an asymmetric space-fractional derivative operator of order $\alpha$ and asymmetry parameter $\beta$, which is defined in the Fourier space as follows [72,73,74,75],

$$\mathcal{F}\left[{}_x{\mathrm{D}}_{\alpha }^{\beta }f\left(x\right)\right]\left(k\right)=-{\psi }_{\alpha }^{\beta }\left(k\right)\mathcal{F}\left[f\left(x\right)\right]\left(k\right),$$

(2)

where:

$${\psi }_{\alpha }^{\beta }\left(k\right)={\left|k\right|}^{\alpha }\left(1-\mathrm{i}\beta \mathrm{sgn}ktan\left(\alpha \pi /2\right)\right),0<\alpha \le 2,-1\le \beta \le 1.$$

(3)

In this paper, we exclude the special case $\alpha =1$ and $\beta \ne 0$ that requires special attention. For convenience, in the analytical point of view we use the characteristic function in the form:

$${\psi }_{\alpha }^{\rho }\left(k\right)={\left(\mathrm{i}k\right)}^{\alpha }{\mathrm{e}}^{-\mathrm{i}\pi \alpha \rho \mathrm{sgn}k},1-\min (1,1/\alpha )<\rho <\min (1,1/\alpha ),$$

(4)

where the following relation between the parameters is established (see details in Appendix A) [63,72]:

$$\rho =\frac{1}{2}+\frac{1}{\alpha \pi }arctan\left(\beta tan(\alpha \pi /2)\right).$$

(5)

Equation (1) is a generalisation of the space-fractional diffusion equation for symmetric Lévy flights [15] in the presence of a $\delta$-sink of strength ${\wp }_{\mathrm{fa}}\left(t\right)$ [53], and describes a target search by a walker exploring its accessible one-dimensional space by performing left-right asymmetrical random jumps represented by the asymmetric space-fractional derivative operator (2). The scale parameter $K_{\alpha }$ (along with the stable index $\alpha$) physically sets the size of the Lévy flight jumps. The skewness parameter $\beta$ may be related to an effective drift or counter-gradient effect [76,77].

The consequence of the point sink at $x=0$ is that the random walker is removed once the target is hit (see Figure 1). Thus, the time dependent weight ${\wp }_{\mathrm{fa}}\left(t\right)$ is the first-arrival probability density. By integration of Equation (1) over the position co-ordinate x, we obtain that ${\wp }_{\mathrm{fa}}\left(t\right)$ is the negative time derivative of the survival probability ${\int }_{-\infty }^{\infty }\mathrm{d}xf(x,t)$ [53], i.e.,

$${\wp }_{\mathrm{fa}}\left(t\right)=-\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{-\infty }^{\infty }\mathrm{d}xf(x,t).$$

(6)

This means that ${\wp }_{\mathrm{fa}}\left(t\right)$ is the probability density function (PDF) of first-arrival: once a random walker arrives at the sink it is annihilated. In what follows, we solve Equation (1) with the help of the Fourier–Laplace transform and various properties of the Mittag–Leffler and Fox H-functions (see Appendix B).

2.1. Diffusion Equation for Asymmetric Lévy Flights

Application of the Fourier transform to Equation (1) yields:

$$\frac{\partial }{\partial t}f(k,t)=-{\mathcal{D}}_{\alpha }{\psi }_{\alpha }^{\rho }\left(k\right)f(k,t)-{\wp }_{\mathrm{fa}}\left(t\right),$$

(7)

where [63,72]:

$${\mathcal{D}}_{\alpha }=\frac{K_{\alpha }}{cos\left(\alpha \pi \right(\rho -1/2\left)\right)}.$$

(8)

Let us first consider Equation (7) without the $\delta$-sink, for the initial condition $f(x,0)=\delta (x-x_0)$. In the Fourier–Laplace space, the solution is given by:

$$f(k,s)=\frac{exp\left(\mathrm{i}k{x}_0\right)}{s+{\mathcal{D}}_{\alpha }{\psi }_{\alpha }^{\rho }\left(k\right)},$$

(9)

from which it follows that:

$$f(k,t)=exp\left(-{\mathcal{D}}_{\alpha }{\psi }_{\alpha }^{\rho }\left(k\right)t+\mathrm{i}k{x}_0\right).$$

(10)

Note that solution (10) can be easily obtained from Equation (7) without applying the Laplace transform. We add here Equation (9) in order to be consistent with the formula, including the sink, see Equation (14) below. Applying the inverse Fourier transform, we find (see details in Appendix C) [73]:

$$f(x,t)=\frac{1}{{\left({\mathcal{D}}_{\alpha }t\right)}^{1/\alpha }}{L}_{\alpha }^{\rho }\left(\frac{x-x_0}{{\left({\mathcal{D}}_{\alpha }t\right)}^{1/\alpha }}\right),0<\alpha \le 2,$$

(11)

where $L_{\alpha }^{\rho }\left(z\right)$ is the Lévy $\alpha$-stable PDF defined as (A22) in terms of the H-function and $z>0$. In Figure 2, we plot the PDF of Lévy flights with different stable index $\alpha$ and asymmetry parameter $\beta$. From Equation (11), by help of the Mellin transform of the H-function (A17) and the symmetry property $L_{\alpha }^{\rho }\left(-z\right)=L_{\alpha }^{1-\rho }\left(z\right)$ [72], we obtain the qth-order moment:

$${\langle \left|x\right|}^q\left(t\right)\rangle =\frac{sin\left(\pi \rho q\right)+sin\left(\pi (1-\rho )q\right)}{\pi q}\Gamma (1+q)\Gamma (1-q/\alpha ){\left({\mathcal{D}}_{\alpha }t\right)}^{q/\alpha },-1<q<\alpha .$$

(12)

One can check that the moments in (12) are positive (see details in Appendix D). Note that for the case $0<\alpha <1$ with $\beta =-1,1$, corresponding to $\rho =0,1$ respectively, the qth-order moment reads:

$${\langle \left|x\right|}^q\left(t\right)\rangle =\frac{\Gamma (1-q/\alpha )}{\Gamma (1-q)}{\left({\mathcal{D}}_{\alpha }t\right)}^{q/\alpha },-\infty <q<\alpha ,$$

(13)

where we use Euler’s reflection formula. For the symmetric case ($\rho =1/2$), we recover the result for symmetric Lévy flights [78,79,80].

2.2. Presence of the $\delta$-Sink and Calculation of the First-Arrival Density

Let us now consider the model, including the $\delta$-sink. The solution of Equation (1), for the initial condition $f(x,0)=\delta (x-x_0)$, in Fourier–Laplace space is given by:

$$f(k,s)=\frac{exp\left(\mathrm{i}k{x}_0\right)-{\wp }_{\mathrm{fa}}\left(s\right)}{s+{\mathcal{D}}_{\alpha }{\psi }_{\alpha }^{\rho }\left(k\right)}.$$

(14)

By integration of Equation (14) over k we find:

$${\int }_{-\infty }^{\infty }\mathrm{d}kf(k,s)=f(x=0,s)=0.$$

(15)

Thus, the first-arrival density in the Laplace domain reads:

$${\wp }_{\mathrm{fa}}\left(s\right)=\frac{I_2\left(s\right)}{I_1\left(s\right)},$$

(16)

where for $I_1\left(s\right)$ and $I_2\left(s\right)$, we obtain (see details in Appendix E):

$$I_1\left(s\right)={\int }_{-\infty }^{\infty }\mathrm{d}k\frac{1}{s+{\mathcal{D}}_{\alpha }{\psi }_{\alpha }^{\rho }\left(k\right)}=\frac{2\pi }{\alpha }\frac{s^{1/\alpha -1}}{{\mathcal{D}}_{\alpha }^{1/\alpha }}\frac{sin\left(\pi \rho \right)}{sin(\pi /\alpha )},\alpha >1$$

(17)

and,

$$\begin{array}{c}\hfill I_2\left(s\right)={\int }_{-\infty }^{\infty }\mathrm{d}k\frac{exp\left(\mathrm{i}k{x}_0\right)}{s+{\mathcal{D}}_{\alpha }{\psi }_{\alpha }^{\rho }\left(k\right)}=\frac{2\pi }{\alpha x_0}\frac{1}{s}{H}_{2,3}^{2,1}\left[\frac{x_0}{{\mathcal{D}}_{\alpha }^{1/\alpha }}{s}^{1/\alpha }\left|\begin{array}{c}(1,1/\alpha ),(1,1-\rho )\hfill \\ (1,1/\alpha ),(1,1),(1,1-\rho )\hfill \end{array}\right.\right].\end{array}$$

(18)

We mention that for $\alpha <1$, $I_2\left(s\right)$ is finite while $I_1\left(s\right)$ diverges for any s, i.e., ${\wp }_{\mathrm{fa}}\left(s\right)=0$. For $\alpha =1$, with $\beta =0$ ($\rho =1/2$), $I_2\left(s\right)$ is finite for $s\ne 0$, whereas $I_1\left(s\right)$ still diverges for $s>0$, thus we have ${\wp }_{\mathrm{fa}}\left(s\right)=0$ again. In the case $\alpha >1$ both $I_1\left(s\right)$ and $I_2\left(s\right)$ converge for finite s, thus, ${\wp }_{\mathrm{fa}}\left(s\right)$ is non-zero [28,55]. Equation (16) then assumes the form:

$$\begin{array}{c}\hfill {\wp }_{\mathrm{fa}}\left(s\right)=\frac{{\mathcal{D}}_{\alpha }^{1/\alpha }}{x_0{s}^{1/\alpha }}\frac{sin(\pi /\alpha )}{sin\left(\pi \rho \right)}{H}_{2,3}^{2,1}\left[\frac{x_0}{{\mathcal{D}}_{\alpha }^{1/\alpha }}{s}^{1/\alpha }\left|\begin{array}{c}(1,1/\alpha ),(1,1-\rho )\hfill \\ (1,1/\alpha ),(1,1),(1,1-\rho )\hfill \end{array}\right.\right],1<\alpha \le 2.\end{array}$$

(19)

By help of the inverse Laplace transform (see Equation (A16)) the first-arrival density yields as:

$${\wp }_{\mathrm{fa}}\left(t\right)=\frac{t^{-1+1/\alpha }}{x_0{\mathcal{D}}_{\alpha }^{-1/\alpha }}\frac{sin(\pi /\alpha )}{sin\left(\pi \rho \right)}{H}_{3,3}^{2,1}\left[\frac{x_0}{{\left({\mathcal{D}}_{\alpha }t\right)}^{1/\alpha }}\left|\begin{array}{c}(1,1/\alpha ),(1,1-\rho ),(1/\alpha ,1/\alpha )\hfill \\ (1,1/\alpha ),(1,1),(1,1-\rho )\hfill \end{array}\right.\right],\alpha >1.$$

(20)

Note that for $\alpha =2$, i.e., $\rho =1/2$, ${\mathcal{D}}_2\equiv K_2$, by using the reduction property of the H-function (see relations (A13) and (A14)), we obtain the Lévy–Smirnov distribution:

$$\begin{array}{cc}\hfill {\wp }_{\mathrm{fa}}\left(t\right)& =\frac{K_2^{1/2}{t}^{-1/2}}{x_0}{H}_{1,1}^{1,0}\left[\frac{x_0}{{\left(K_{2}t\right)}^{1/2}}\left|\begin{array}{c}(1/2,1/2)\hfill \\ (1,1)\hfill \end{array}\right.\right]=\frac{2x_0}{\sqrt{K_2{t}^3}}{H}_{1,1}^{1,0}\left[\frac{x_0^2}{K_{2}t}\left|\begin{array}{c}(-1/2,1)\hfill \\ (-1,2)\hfill \end{array}\right.\right]\hfill \\ \hfill & =\frac{x_0}{\sqrt{4\pi K_2{t}^3}}{H}_{0,1}^{1,0}\left[\frac{x_0^2}{4K_{2}t}\left|\begin{array}{c}\\ \stackrel{\_}{(0,1)}\hfill \end{array}\right.\right]=\frac{x_0}{\sqrt{4\pi K_2{t}^3}}exp\left(-\frac{x_0^2}{4K_{2}t}\right).\hfill \end{array}$$

(21)

In the latter equality, we used Equation (A10) and the duplication rule $2^{2z}\Gamma \left(z\right)\Gamma (z+1/2)=2\sqrt{\pi }\Gamma \left(2z\right)$ of the $\Gamma$-function. By help of Equations (A11) and (A12) the first-arrival density (20) can be expressed in the following form:

$${\wp }_{\mathrm{fa}}\left(t\right)=\frac{t^{-1+1/\alpha }}{x_0{\mathcal{D}}_{\alpha }^{-1/\alpha }}\frac{\alpha sin(\pi /\alpha )}{sin\left(\pi \rho \right)}{H}_{3,3}^{1,2}\left[\frac{{\mathcal{D}}_{\alpha }t}{x_0^{\alpha }}\left|\begin{array}{c}(0,1),(0,\alpha ),(0,\alpha -\alpha \rho )\hfill \\ (0,1),(0,\alpha -\alpha \rho ),(1-1/\alpha ,1)\hfill \end{array}\right.\right].$$

(22)

Inserting the series expansion of the H-function (A8) we obtain:

$$\begin{array}{cc}\hfill {\wp }_{\mathrm{fa}}\left(t\right)& =\frac{t^{-1+1/\alpha }}{x_0{\mathcal{D}}_{\alpha }^{-1/\alpha }}\frac{\alpha sin(\pi /\alpha )}{sin\left(\pi \rho \right)}\sum _{k=1}^{\infty }\frac{\Gamma \left(1+k\right)\Gamma \left(1+\alpha k\right)}{\Gamma \left(1+\alpha k-\alpha \rho k\right)\Gamma \left(1/\alpha +k\right)\Gamma \left(\alpha \rho k-\alpha k\right)}\frac{{(-1)}^k}{k!}{\left(\frac{{\mathcal{D}}_{\alpha }t}{x_0^{\alpha }}\right)}^k\hfill \\ \hfill & =\frac{t^{-1+1/\alpha }}{x_0{\mathcal{D}}_{\alpha }^{-1/\alpha }}\frac{\alpha sin(\pi /\alpha )}{\pi sin\left(\pi \rho \right)}\sum _{k=1}^{\infty }{(-1)}^{k-1}sin\left(\pi \alpha (1-\rho )k\right)\frac{\Gamma \left(1+\alpha k\right)}{\Gamma \left(1/\alpha +k\right)}{\left(\frac{{\mathcal{D}}_{\alpha }t}{x_0^{\alpha }}\right)}^k.\hfill \end{array}$$

(23)

This form can be shown to be equivalent to the corresponding expression in Theorem 3.14 of [52], in a dimensionless form. Therefore, from the equation above, the short time limit of the first-arrival density has the asymptotic behaviour

$${\wp }_{\mathrm{fa}}\left(t\right)\sim \frac{\alpha sin(\pi /\alpha )sin\left(\pi \alpha \right(1-\rho \left)\right)}{\pi sin\left(\pi \rho \right)}\frac{\Gamma (1+\alpha )}{\Gamma (1/\alpha +1)}\frac{{\mathcal{D}}_{\alpha }^{1+1/\alpha }}{x_0^{1+\alpha }}{t}^{1/\alpha }.$$

(24)

We also provide another method to get the above expression in Appendix F.

In order to obtain the long time limit, we use the Laplace transform of the first-arrival density (19) at small s, and then apply the Tauberian theorem [81] to get (see Appendix for details):

$${\wp }_{\mathrm{fa}}\left(t\right)\sim \frac{sin\left(\pi \alpha \right(1-\rho \left)\right)}{sin\left(\pi \rho \right)}\frac{\Gamma (2-\alpha )}{{\Gamma }^2(1/\alpha )\Gamma (1-1/\alpha )}\frac{{\mathcal{D}}_{\alpha }^{-1+1/\alpha }}{x_0^{1-\alpha }}{t}^{-2+1/\alpha }.$$

(25)

This expression is in agreement with the results in [53] for the symmetric case ($\beta =0$). Note that the short and long time limit behaviours are also in agreement with intuition, namely, the first-arrival density initially increases from zero as a power of time t, and asymptotically decreases with a negative power of t (we note that this derivation assumed $\alpha >1$). Note also that, with the variation of $\rho$ (see Equation (5)), Equations (24) and (25) are non-negative.

The fractional-order moments of the first-arrival density with $1<\alpha <2$ read:

$$\begin{array}{cc}\hfill \langle t^q\rangle ={\int }_0^{\infty }\mathrm{d}t{t}^q{\wp }_{\mathrm{fa}}\left(t\right)& =\frac{sin(\pi /\alpha )sin\left(\pi \right(1-\rho \left)\right(1+\alpha q\left)\right)}{sin\left(\pi \rho \right)sin\left(\frac{\pi }{\alpha }(1+\alpha q)\right)}\frac{\Gamma (1-\alpha q)}{\Gamma (1-q)}\frac{x_0^{\alpha q}}{{\mathcal{D}}_{\alpha }^q},\hfill \\ \hfill & -1-\frac{1}{\alpha }<q<1-\frac{1}{\alpha },\beta \ne 1,\hfill \end{array}$$

(26)

where we use the Mellin transform of the H-function (see Equation (A17)). This form can be shown to be equivalent to the corresponding result in [52] (see Theorem 3.9). Due to the definition of $\rho$ (Equation (5)), and from variations of q, it is obvious that the moments are positive. For the case $\beta =1$ with $1<\alpha <2$, we have $\rho =1-1/\alpha$, and thus arrive at:

$$\langle t^q\rangle =\frac{\Gamma (1-\alpha q)}{\Gamma (1-q)}\frac{x_0^{\alpha q}}{{\mathcal{D}}_{\alpha }^q},-\infty <q<\frac{1}{\alpha },\beta =1,$$

(27)

while for the limit of Brownian motion ($\alpha =2$), we get:

$$\langle t^q\rangle =\frac{\Gamma (1-2q)}{\Gamma (1-q)}\frac{x_0^{2q}}{K_2^q},-\infty <q<\frac{1}{2}.$$

(28)

In Figure 3, we plot the fractional-order moments of symmetric ($\beta =0$) and the extremal two-sided ($\beta =1$) Lévy flights as functions of the parameter q for different values of the stable index $\alpha$. As can be seen, the PDFs are normalised, $\langle t^{q=0}\rangle =1$, but have diverging means for all $\alpha$.

3. Langevin Equation Approach

A complementary view to the calculation of time dependent probabilities is the numerical study of individual search paths. For such a simulation of asymmetric Lévy flights, we use the (stochastic) Langevin equation for the particle position x [67],

$$\frac{\mathrm{d}x\left(t\right)}{\mathrm{d}t}=K_{\alpha }^{1/\alpha }\zeta \left(t\right),$$

(29)

driven by random noise $\zeta \left(t\right)$ with Lévy stable distribution, which is characterised by the stable index $\alpha$ and the asymmetry parameter $\beta$. For simulation purposes we pass to the discretised version of the Langevin equation,

$$x_{n+1}-x_n={(K_{\alpha }\Delta t)}^{1/\alpha }{\zeta }_t,$$

(30)

taken at times $t_n=n\Delta t$, $n=0,1,2,\cdots$, for more details, see [67]. The noise ${\zeta }_t$ is computed according to the algorithm described in [82]. For each fixed value of $\alpha$, $\beta$ and $K_{\alpha }$, the random walker starts its motion at $x_0$. When it hits the vicinity of the origin $x=0$, the corresponding event time $t=n\Delta t$ is recorded and a new particle is released. This procedure is repeated $2\times {10}^6$ times. The time steps to determine the long-time and short-time behaviours of the first-arrival density were $\Delta t=0.001$ and $0.0001$, respectively. For the very short-time behaviour when $t\to 0$ we chose $\Delta t=0.00001$. While we study analytically a point-like target, in the simulations, we have to endow the target with a finite size d. This size should be sufficiently large to be detected by the simulated walk—and small enough to be considered as a point-like target and thus warrant proper statistics. This issue was discussed in detail in [9,54]. We find that the appropriate size of the target should be $d\sim 2{({\mathcal{D}}_{\alpha }\Delta t)}^{1/\alpha }$, as a function of the time step $\Delta t$, stable index $\alpha$, asymmetry parameter $\beta$ and the generalised diffusion coefficient $K_{\alpha }$ (see Equations (5) and (8)). Note that in contrast to the symmetric case $\beta =0$ studied in [9,54], here, the coefficient ${\mathcal{D}}_{\alpha }$ depends on both $\alpha$ and $\beta$, such that we should choose a smaller time step for simulations with $\alpha$ closer to unity.

Figure 4 shows the long-time asymptotic behaviour of the first-arrival PDF along with the short-time behaviour and the behaviour around $t=0$, respectively. The coloured solid curves in the left and centre panels are obtained from the exact analytical solution (20), and the symbols show the results based on numerical solution of the Langevin Equation (30). The short black lines in the right panels represent the long-time asymptotic behaviour $t^{-2+1/\alpha }$ (see Equation (25)). From the short-time behaviour of the first-arrival density in Figure 4 (left and centre panels), it can be seen that, initially (at $t\le 0.1$), particles moving with negative skewness ($\beta <0$) have a higher chance to find the target as compared with particles with non-negative $\beta$. In physical terms, the negative skewness corresponds to shorter jumps to the right occurring with higher frequency and longer jumps occurring with lower frequency to the left. However, it is seen from the centre panels that at times longer than $0.1$ the random walker with positive $\beta$ will detect the target with higher probability. This is intuitively clear since such a walker performs more short jumps to the negative direction (towards the target) for positive $\beta$, and when $\alpha$ tends to $1^{+}$ the number of short jumps increases. Thus, the area on the left of the initial position $x_0$ is densely covered, whereas the random searcher with $\beta <0$ performs long jumps to the left direction and may more easily overshoot the target.

The short-time behaviour of the first-arrival density of asymmetric Lévy flights is depicted in Figure 5 for different values of the stable index $\alpha$ and diffusion coefficient $K_{\alpha }=10$. The results obtained from simulations (coloured symbols) are in good agreement with the theoretical result (20). One can see that for the case $\beta =1$ (left panel), the searcher with smaller stable index $\alpha$ will find the target faster since, in this case, there are more short jumps in the direction of the target. Respectively, in case of the negative skewness $\beta =-0.5$ (right panel) the searcher with larger $\alpha$ will find the target faster. These observations—valid for close initial distance between the searcher and the target ($x_0=1$)—are confirmed by calculation of the search efficiency in the next Section. We note that, based on the simulations whose results are shown in Figure 4 and Figure 5 for different time scales, the first-arrival probability densities are unimodal.

4. Calculation of the Search Efficiency

In order to study the search efficiency, we make use of the definition proposed in [9], which is defined by:

$$\mathcal{E}=〈\frac{1}{t}〉={\int }_0^{\infty }\mathrm{d}t\frac{{\wp }_{\mathrm{fa}}\left(t\right)}{t}={\int }_0^{\infty }\mathrm{d}s{\wp }_{\mathrm{fa}}\left(s\right).$$

(31)

A given search strategy is optimal when the efficiency $\mathcal{E}$ of the corresponding search process reaches its maximum. The constructive meaning of such a definition of efficiency was tested in [9] for Brownian motion with drift towards and away from the target. The expediency of the definition was addressed in the first arrival problems for symmetric Lévy flights with drift in [9], for combined Lévy–Brownian and Lévy–Lévy search in [55,83] of single and multiple targets. In all the cases, such a definition gives reasonable results, which allow explanation at an intuitive level. Another strong motivation is its analytical simplicity, which allows for getting exact analytical formulas in cases where the most probable value of the first arrival time can be calculated only numerically.

The exact result for the search efficiency can be obtained by substituting Equation (19) or (20) into the above expression and then using equation (A17). We find that:

$$\mathcal{E}=\frac{sin\left[\pi (1-\rho )(\alpha -1)\right]}{sin\left(\pi \rho \right)}\Gamma (\alpha +1)\frac{{\mathcal{D}}_{\alpha }}{x_0^{\alpha }}.$$

(32)

Consequently, in general, the efficiency depends on the positive parameter $\rho$ and, hence, on the skewness parameter $\beta$. In particular, for the symmetric case $\rho =1/2$ or $\beta =0$ it turns out that [9]:

$$\mathcal{E}=sin\left(\frac{\pi }{2}(\alpha -1)\right)\Gamma (\alpha +1)\frac{{\mathcal{D}}_{\alpha }}{x_0^{\alpha }},$$

(33)

and in the limit $\alpha =2$ for a Brownian walker, the efficiency is given by:

$$\mathcal{E}=\frac{2{\mathcal{D}}_2}{x_0^2}.$$

(34)

In the Cauchy limit $\alpha =1$ ($\beta =0$), the efficiency goes to zero. For the case $\alpha \to 1^{+}$ with $\beta \ne 0$, substituting Equations (5) and (8) into Equation (32) yields:

$$\begin{array}{c}\hfill \mathcal{E}=\frac{sin\left[\left(\frac{\pi }{2}-\frac{1}{\alpha }arctan\left(\beta tan\left(\frac{\alpha \pi }{2}\right)\right)\right)\left(\alpha -1\right)\right]}{cos\left[\frac{1}{\alpha }arctan\left(\beta tan\left(\frac{\alpha \pi }{2}\right)\right)\right]cos\left[arctan\left(\beta tan\left(\frac{\alpha \pi }{2}\right)\right)\right]}\frac{\Gamma (\alpha +1)K_{\alpha }}{x_0^{\alpha }}.\end{array}$$

(35)

Finally, taking the limit $\alpha \to 1^{+}$, one gets (see details in Appendix G):

$$\begin{array}{c}\hfill \mathcal{E}=\left\{\begin{array}{cc}\infty \hfill & \beta >0\hfill \\ \frac{2\left|\beta \right|}{\pi }\frac{K_{\alpha }}{x_0}\hfill & \beta <0\hfill \end{array}\right..\end{array}$$

(36)

In Figure 6, Figure 7 and Figure 8, we use Equation (32) to compare the search efficiency for different values of the stable index $\alpha$, skewness $\beta$ and initial distance $x_0$. In Figure 6, we show the dependence of the search efficiency on the stable index $\alpha$ for different values of $\beta$. We observe that when the target is located on the left of the random searcher at the initial distance $x_0$, it is more probable to be detected by Brownian motion than by Lévy flights with non-positive asymmetry parameter ($\beta \le 0$). This is intuitively clear: in contrast to the Brownian walker, such a Lévy searcher performs long jumps to the left, which may easily lead him far away from the target located at $x=0$. Vice versa, Lévy flights with $\beta >0$ are always more efficient than Brownian search since for such searchers long jumps to the left are forbidden, while they perform more shorter jumps to the left than the Brownian walker. Such an intuitive picture stems from the shapes of the $\alpha$-stable PDFs depicted in the left panel of Figure 2.

In Figure 7, we illustrate the dependency of the search efficiency on the asymmetry parameter $\beta$ for different stable index $\alpha$ and for two different initial target distances. It can be seen that for the close target ($x_0=0.5$) Brownian motion is more efficient than the Lévy search (except for Lévy flights with $\beta >0$ and $\alpha \to 1$). By increasing the initial distance, Lévy flights are more efficient as compared to the Brownian motion (except for $\alpha =1$ with $\beta =0$). Moreover, the asymmetry in the jumps leads to a higher efficiency, and this effect is more pronounced for stable indices close to unity.

Finally, we study the effect of the initial distance $x_0$ on the search efficiency for different values of the asymmetry $\beta$ in Figure 8. It was observed that for smaller values of $\alpha$, the asymmetry of the jumps led to higher search efficiency, and by increasing $\alpha \to 2$ at a fixed value of $x_0$, all Lévy flights have almost the same search efficiency.

Another important quantity is the search reliability, defined as the cumulative arrival probability [9],

$$\mathcal{P}={\int }_0^{\infty }\mathrm{d}t{\wp }_{\mathrm{fa}}\left(t\right).$$

(37)

From relation (20), we find that $\mathcal{P}=1$, which means that in the case of asymmetric Lévy flights, as in the case of symmetric Lévy flights, the searcher will always find the target, and this statement is valid for $\alpha >1$ with all $\beta$ ($-1\le \beta \le 1$).

5. Conclusions

We investigated the dynamics and performance of asymmetric Lévy flights as a random search strategy by solving the deterministic Fokker–Planck equation with asymmetric space-fractional derivative in the presence of a sink, which mimics a point target. We found exact results for the resulting first-arrival or first-hitting density, its fractional order moments and search efficiency in terms of the stable index $\alpha$ and the skewness parameter $\beta$ of the asymmetric Lévy flight. We observed that the first-arrival density is identically zero for processes with $\alpha \le 1$. For $\alpha >1$ the first-arrival density increases at short times as $t^{1/\alpha }$ and decreases at long times as $t^{-2+1/\alpha }$. Thus, the exponents of the short- and long-time asymptotics are independent of the skewness $\beta$.

The efficiency of the search depends on both parameters, the stable index and skewness, as well as on the initial distance $x_0$ of the searcher to the target and the diffusion coefficient $K_{\alpha }$. We demonstrated that for a short initial distance $x_0$, the Brownian strategy is more efficient than Lévy flights, except when $\beta >0$ with stable index $\alpha \to 1$. For long initial distances, Lévy flights become more efficient than the Brownian motion. These observations extend the results obtained in [9] to the asymmetric case. We also found that the asymmetry in jump yields higher efficiency, especially for stable indices close to unity. In the limiting case $\alpha \to 1^{+}$, the search efficiency is infinite for positive skewness, and the situation drastically changes for a negative skewness. It is observed that the search efficiency in this limit case has a linear behaviour with respect to $\left|\beta \right|$.

It will be interesting to analyse the case of asymmetric Lévy flight search on multiple targets, as was done for symmetric Lévy flights [83]. Another interesting search strategy could be the case of combined search [55]: asymmetric Lévy flights with local Brownian search. It is also of interest to study the random Lévy search in two and three dimensions, as well as applications of such non-local asymmetric search processes in finding the global extremum of potential functions [59]. Different definitions of the search efficiency are also worth being analysed and compared.