Generalized diffusion and random search processes

Searching problems are ubiquitous in the natural world. They occur on many fields. In our daily life, searching for a lost key is a representative search problem [1]. In order to produce a chemical reaction, transcription factor proteins search for a specific binding spot on a DNA chain in biology [2, 3]. In addition, the behaviors of animals foraging are a branch of the search problem [4, 5]. Mathematically, computer algorithms searching for minima in a complex search space [6] are classified as a search problem. Brownian motion with Gaussian distributed jump length was considered as an optimal search strategy until Shlesinger and Klafter suggested that Lévy flight is more efficient when searching for sparse targets in 1986 [7]. The typical characteristic of one- and two-dimensional Brownian search lies in oversampling, which frequently arrives at the previously visited point. The Brownian search strategy is sufficient when we have a non-destructive search scenario where the target does not get depleted from being visited by the searcher, but in the case of destructive foraging, where the target becomes undetectable after visiting it once, this strategy becomes inefficient [5, 7–9]. This oversampling behavior can be avoided by Lévy flight search and in cases of destructive foraging it becomes the optimal search strategy since the scale free of jump length [10, 11], where the distribution of jump length satisfies power-law distribution, i.e. λ(x) ∼ 1/|x|^1+β with 0 < β < 2 [12–14].

We focus on the search problems in the discipline of movement ecology [4, 15], where the movement patterns of animals are saltatory, occupied with long excursions. One of the classical examples is albatross birds searching for food [16, 17]. Within this movement ecology the Lévy flight hypothesis is formulated: it states that the scale-free search strategies called Lévy flight search minimise the search time [5, 8, 16] and is the optimal search strategy under a broad set of conditions, such as target resource distribution (sparse/abundant resources), target revisitability (destructive/non-destructive targets) and prey targeting [18]. For albatross birds searching for sparse food, it is a good choice to search with power-law distributed jump length. It is discovered that there are a large number of animals searching with such scale-free strategies. For instance, deer and goats [19, 20], bees [21] as well as marine predators [22, 23]. Though there are some questions [17, 24, 25] about the Lévy flight search strategies of albatross birds, in [26] it implies that the searching strategies of individual albatross birds are indeed Lévy flights.

In order to verify whether the searcher can locate the target or not, search reliability is defined. However, it cannot be used to judge a search strategy since it does not provide any information about the efficiency of a search strategy. Correspondingly, search efficiency is defined, which is used to choose the optimal search strategies and is a crucial quantity in search problem [27]. Based on the space-fractional Fokker–Planck equation as well as search reliability and search efficiency, [28] compares the Brownian search with the Lévy flight search, it found that the search efficiency depends on the initial searcher–target distance, Brownian search is a good choice when the initial searcher–target distance is smaller than some critical value, otherwise Lévy flight optimizes the search process. [28] also discusses the problem of searching for multiple targets. Moreover, the search efficiencies are compared in [29] in the presence of an external drift. It shows that Brownian searcher is always the most efficient for short initial separation or when the target is in the downhill case for larger initial separation, whereas Lévy flight turns out to be efficient search strategy when the target is in the uphill case. In addition, due to the complex biological reality, the search reliability and search efficiency of combined Lévy–Brownian motion are discussed in [30].

Beginning with the Montroll–Weiss equation [31, 32], which is derived in the frame of continuous time random walk (CTRW) model, in [33] a generalized diffusion equation in modified or Riemann–Liouville form is obtained, which generalizes the time fractional diffusion equation [14, 34]. In this paper we aim to study the search process based on the generalized diffusion equation for different memory kernels. The aim of introducing such search problems is manifold. Apart from the application in animal foraging, the considered models could be of interest for analysis of the propagation of contamination particles in the environment. The CTRW model [14, 36, 37] was applied, for example, to infiltration of anomalously diffusing particles from one material to another through a biased interface [38], and recently in the long-term prediction of ¹³⁷Cs concentration in the lakes after the Fukushima Daiichi Nuclear Power Plant accident [39], where the contamination particle ${\left(\,\right.}^{137}\mathrm{C}\mathrm{s}\left.\right)$ follows a long-tailed waiting time between jumps. Additionally to the ground and underground water pollution, such CTRW models could apply to model propagation of sandstorm particles, including pathogens or microbes carried by the sandstorm [40–42]. The sandstorm particles fly over a long distance before they arrive at the particular area and stay around for some time while the next flight takes place. The long flight of sandstorm particles may be interrupted by long waiting times. By studying the first passage and/or arrival time distribution of the contamination particle to the particular point of interest, one can map the area affected by the contamination particles at some given time.

This paper is organized as follows. In section 2, we introduce the Brownian search model with trapping in terms of the generalized diffusion equation, calculate the search reliability and search efficiency for different memory kernels. In section 3, we turn to discuss the Lévy flight search with trapping in the presence of one target. Finally, we conclude the paper with some discussions in section 4. In the appendices, we give the properties of Mittag–Leffler and Fox H-functions, as well as general results for Lévy flight search with trapping in the presence of multiple targets and combined Lévy flight search in the presence of a single target.

We consider one dimensional search with trapping, which is given by the following generalized diffusion equation with a δ-sink of a strength ℘_fa(t)

$$\begin{equation}\frac{\partial }{\partial t}P(x,t)=\mathcal{D}\frac{\partial }{\partial t}{\int }_{0}^{t}\eta (t-{t}^{\prime })\frac{{\partial }^{2}}{\partial {x}^{2}}P(x,{t}^{\prime })\,\mathrm{d}{t}^{\prime }-{\wp }_{\text{fa}}(t)\delta (x-X),\end{equation} \tag{ 1 }$$

where the initial condition meets P(x, t = 0) = δ(x − x₀). Here, $\mathcal{D}$ is a diffusion coefficient and η(t) represents the memory kernel, satisfying ${\mathrm{lim}}_{s\to 0}\frac{1}{\hat{\eta }(s)}=0$ and $\frac{1}{\hat{\eta }(s)}$ is a complete Bernstein function [33, 43–45]. In this equation, δ-sink means that the random searcher positioned at x = x₀ at the beginning, will be removed when it arrives for the first time at x = X, i.e. P(x = X, t) = 0. Thus, ℘_fa(t) represents the first arrival time density (FATD) [46, 47].

Here we note that in the absence of a sink, equation (1) is reduced to the generalized diffusion equation [33, 43, 48]

$$\begin{equation}\frac{\partial }{\partial t}P(x,t)=\mathcal{D}\frac{\partial }{\partial t}{\int }_{0}^{t}\eta (t-{t}^{\prime })\frac{{\partial }^{2}}{\partial {x}^{2}}P(x,{t}^{\prime })\,\mathrm{d}{t}^{\prime },\end{equation} \tag{ 2 }$$

describing the anomalous diffusive process. This equation can be obtained from the CTRW model [37] with a memory-dependent waiting time probability density function (PDF) ψ(t), which in Laplace space is given by $\hat{\psi }(s)=1/(1+{\tau }_{\eta }{\hat{\eta }}^{-1}(s))$, and jump length PDF λ(x) with finite variance, which in Fourier space is $\tilde{\lambda }(k)\sim 1-\frac{{\sigma }^{2}}{2}{k}^{2}$. The generalized diffusion coefficient in equation (2) reads $\mathcal{D}={\sigma }^{2}/[2{\tau }_{\eta }]$.

Integrating equation (1) with respect to x from −∞ to ∞, we get the expression of FATD, which reads

$$\begin{equation}{\wp }_{\text{fa}}(t)=-\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{-\infty }^{\infty }P(x,t)\,\mathrm{d}x=-\frac{\mathrm{d}}{\mathrm{d}t}\mathcal{S}(t).\end{equation} \tag{ 3 }$$

It is the negative time derivative of the survival probability $\mathcal{S}(t)={\int }_{-\infty }^{\infty }P(x,t)\,\mathrm{d}x$. Two important quantities characterising the diffusion process respectively are search reliability and search efficiency. The search reliability quantifies the extent to which the search process can locate the trap, which is accounted as the cumulative arrival probability [47]

$$\begin{equation}\mathcal{P}={\int }_{0}^{\infty }{\wp }_{\text{fa}}(t)\,\mathrm{d}t.\end{equation} \tag{ 4 }$$

If the search reliability equals one, $\mathcal{P}=1$, the searcher will find the target with the probability 1. With the aid of Laplace transform ⁷ the search reliability can be simplified as

$$\begin{equation}\mathcal{P}={\hat{\wp }}_{\text{fa}}(s=0),\end{equation} \tag{ 5 }$$

where ${\hat{\wp }}_{\text{fa}}(s)$ is the Laplace image of ℘_fa(t). The search reliability may be relevant for choosing the search strategy. A large value of search reliability corresponds to a high success probability to which the searcher can locate the target. It can be smaller than unity for example for Brownian search with drift, which pushes the particle in the opposite direction from the target location [47].

The search efficiency can be defined in different ways either as a ratio of the visited number of targets and the number of steps or as a ratio of the visited number of targets and the distance travelled. The first definition can be applied when a searcher moves with jumps, while the second definition applies for example to Brownian motion search and finite-velocity Lévy walks search [28–30]. We focus on searching for a single or a finite number of targets with saltatory motion in this paper. In our continuous time model, we argued that the search efficiency can be described as a ratio of the visited number of targets and the time of the process with a proper averaging. Naturally, it relates to the mean first arrival time, $\langle t\rangle ={\int }_{0}^{\infty }t{\wp }_{\text{fa}}(t)\,\mathrm{d}t=-{\left.\frac{\partial {\hat{\wp }}_{\text{fa}}(s)}{\partial s}\right\vert }_{s=0}$, in two ways, as $\mathcal{E}=1/\left\langle t\right\rangle$, or as the averaged inverse search time [47]

$$\begin{equation}\mathcal{E}=\left\langle \frac{1}{t}\right\rangle ={\int }_{0}^{\infty }\frac{{\wp }_{\text{fa}}(t)}{t}\mathrm{d}t.\end{equation} \tag{ 6 }$$

This definition of efficiency (6) was used for Brownian motion with drift [47], symmetric Lévy flight search with drift [47], for combined Lévy–Brownian and Lévy–Lévy search [28–30] of single and multiple targets, as well as asymmetric Lévy flight search [49], and in all cases such an approach gives reasonable results. Moreover, ⟨t⟩ diverges in some situations since if we consider, for example, Brownian search, then some of the particles will go in the opposite direction of the target and will never reach it, which by averaging will lead to infinite mean time to reach the target [50]. Also, ⟨t⟩ is diverging when an external bias initially pushes the Lévy flight searcher towards the target [29]. In this case, due to the presence of leapovers a Lévy flight searcher may miss the target [29]. Furthermore, the efficiency of a given search strategy could be also analyzed in terms of calculation of the most probable value of the first arrival time, which can be obtained by finding the first derivative with respect to time t of the FATD. Thus, finding the most probable value of the first arrival time may require numerical calculations. In this paper, we calculate the most probable value of the first arrival time for some cases of the memory kernel analytically and for some numerically. However, detailed analysis and comparison between different definitions of the efficiency require further extensive numerical study and we leave this problem for some future work.

Here, we utilize (6) as a definition of search efficiency. A search strategy is optimal when the search efficiency of the corresponding diffusion process is maximal. Using the properties of Laplace transform

$$\begin{equation}\mathcal{L}\left\{\frac{f(t)}{t}\right\}(s)={\int }_{0}^{\infty }{\mathrm{e}}^{-st}\frac{f(t)}{t}\,\mathrm{d}t={\int }_{s}^{\infty }\hat{f}(u)\mathrm{d}u,\end{equation} \tag{ 7 }$$

where $\hat{f}(u)$ is the Laplace image of f(t). The search efficiency (6) can be written in the following form [47]

$$\begin{equation}\mathcal{E}={\int }_{0}^{\infty }{\hat{\wp }}_{\text{fa}}(s)\,\mathrm{d}s.\end{equation} \tag{ 8 }$$

Thus, both the search reliability and the search efficiency are related to the Laplace image of FATD, i.e. ${\hat{\wp }}_{\text{fa}}(s)$. In the following, we aim to derive the expression of ${\hat{\wp }}_{\text{fa}}(s)$.

By Laplace and Fourier transformation ⁸ of equation (1), we find

$$\begin{equation}s\tilde{\hat{P}}(k,s)-{\mathrm{e}}^{{\imath}k{x}_{0}}=-\mathcal{D}{k}^{2}s\hat{\eta }(s)\tilde{\hat{P}}(k,s)-{\hat{\wp }}_{\text{fa}}(s){\mathrm{e}}^{{\imath}kX},\end{equation} \tag{ 9 }$$

from there we have

$$\begin{equation}\tilde{\hat{P}}(k,s)=\frac{1}{s\hat{\eta }(s)}\frac{{\mathrm{e}}^{{\imath}k{x}_{0}}-{\hat{\wp }}_{\text{fa}}(s){\mathrm{e}}^{{\imath}kX}}{\frac{1}{\hat{\eta }(s)}+\mathcal{D}{k}^{2}}.\end{equation} \tag{ 10 }$$

The inverse Fourier transform ⁹ of equation (10) yields

$$\begin{equation}\hat{P}(x,s)=\frac{{s}^{-1}}{2\sqrt{\mathcal{D}\hat{\eta }(s)}}\left[{\mathrm{e}}^{-\frac{\vert x-{x}_{0}\vert }{\sqrt{\mathcal{D}\hat{\eta }(s)}}}-{\hat{\wp }}_{\text{fa}}(s)\,{\mathrm{e}}^{-\frac{\vert x-X\vert }{\sqrt{\mathcal{D}\hat{\eta }(s)}}}\right].\end{equation} \tag{ 11 }$$

From the condition $\hat{P}(x=X,s)=0$ for the FATD we obtain

$$\begin{equation}{\hat{\wp }}_{\text{fa}}(s)={\mathrm{e}}^{-\frac{\vert X-{x}_{0}\vert }{\sqrt{\mathcal{D}\hat{\eta }(s)}}}.\end{equation} \tag{ 12 }$$

By substitution of FATD in equation (11), for the $\hat{P}(x,s)$ we find

$$\begin{equation}\hat{P}(x,s)=\frac{{s}^{-1}}{2\sqrt{\mathcal{D}\hat{\eta }(s)}}\left[{\mathrm{e}}^{-\frac{\vert x-{x}_{0}\vert }{\sqrt{\mathcal{D}\hat{\eta }(s)}}}-{\mathrm{e}}^{-\frac{\vert x-X\vert +\vert X-{x}_{0}\vert }{\sqrt{\mathcal{D}\hat{\eta }(s)}}}\right].\end{equation} \tag{ 13 }$$

This result can be obtained if one uses the subordination approach, such as in [51].

In addition, combined with (5) we find that the search reliability equals $\mathcal{P}={\hat{\wp }}_{\text{fa}}(s=0)=1$ since the memory kernel satisfies ${\mathrm{lim}}_{s\to 0}\frac{1}{\hat{\eta }(s)}=0$, which means the searcher will always arrive at the target eventually. We calculate the search efficiency based on (12), which is expressed by

$$\begin{equation}\mathcal{E}={\int }_{0}^{\infty }{\mathrm{e}}^{-\frac{\vert X-{x}_{0}\vert }{\sqrt{\mathcal{D}\hat{\eta }(s)}}}\,\mathrm{d}s.\end{equation} \tag{ 14 }$$

Next, we will calculate the search efficiency for some representative memory kernels η(t).

Here we note that for the case without sink, we find the solution of equation (2) in the form (x₀ = 0)

$$\begin{equation}\hat{P}(x,s)=\frac{{s}^{-1}}{2\sqrt{\mathcal{D}\hat{\eta }(s)}}{\mathrm{e}}^{-\frac{\vert x\vert }{\sqrt{\mathcal{D}\hat{\eta }(s)}}}.\end{equation} \tag{ 15 }$$

From here, one can find the corresponding mean squared displacement (MSD) [33], i.e.

$$\begin{equation}\left\langle {x}^{2}(t)\right\rangle ={\int }_{-\infty }^{\infty }{x}^{2}P(x,t)\,\mathrm{d}x=2\mathcal{D}\,{\mathcal{L}}^{-1}\left[{s}^{-1}\hat{\eta }(s)\right].\end{equation} \tag{ 16 }$$

2.1. Brownian search

The standard case of Brownian search is obtained for η(t) = 1, i.e. $\hat{\eta }(s)=1/s$. In this situation, the waiting time PDF behaves as Poisson distribution $\psi (t)={\mathcal{L}}^{-1}\left[\frac{1}{1+\tau \,s}\right]=\frac{1}{\tau }{\mathrm{e}}^{-t/\tau }$ and MSD in the absence of a sink has a linear dependence on time $\left\langle {x}^{2}(t)\right\rangle =2\mathcal{D}t$. The FATD then becomes

$$\begin{equation}{\hat{\wp }}_{\text{fa}}(s)={\mathrm{e}}^{-\frac{{s}^{1/2}}{\sqrt{\mathcal{D}}}\vert X-{x}_{0}\vert },\end{equation} \tag{ 17 }$$

i.e.

$$\begin{equation}{\wp }_{\text{fa}}(t)=\frac{\vert X-{x}_{0}\vert }{\sqrt{4\pi \mathcal{D}{t}^{3}}}{\mathrm{e}}^{-\frac{{(X-{x}_{0})}^{2}}{4\mathcal{D}t}},\end{equation} \tag{ 18 }$$

which is the Lévy–Smirnov distribution with power-law decay t^−3/2 in the long time limit [46], see figure 1. This result coincides with the first passage time density [50] since the particle does not perform the jump-like motion and the first passage time will be the same as the first arrival/hitting time. The reliability equals one, $\mathcal{P}=1$, while the efficiency is [47]

$$\begin{equation}\mathcal{E}={\int }_{0}^{\infty }{\mathrm{e}}^{-\frac{{s}^{1/2}}{\sqrt{\mathcal{D}}}\vert X-{x}_{0}\vert }\mathrm{d}s=\frac{2\mathcal{D}}{{(X-{x}_{0})}^{2}}.\end{equation} \tag{ 19 }$$

We note that the mean time of the Brownian particle to hit the target, $\langle t\rangle ={\int }_{0}^{\infty }t\,{\wp }_{\text{fa}}(t)\,\mathrm{d}t$, diverges [50] since some particles can move in the wrong direction and never hit the target. By taking the first derivative with respect to time t of the FATD, we figure out the most probable value of the first arrival time t_*, which behaves as ${t}_{\ast }=\frac{{(X-{x}_{0})}^{2}}{6\mathcal{D}}$. Thus, if we define the efficiency as the reciprocal of the most probable value of first arrival time t_*, we find ${\mathcal{E}}_{\ast }=\frac{1}{{t}_{\ast }}=\frac{6\mathcal{D}}{{(X-{x}_{0})}^{2}}$. Moreover, for a fixed initial searcher–target distance, increasing diffusivity of the searcher can improve the search efficiency.

Zoom In Zoom Out Reset image size

2.2. Power-law memory kernel: anomalous diffusive search

We consider the power-law memory kernel $\eta (t)=\frac{{t}^{\alpha -1}}{{\Gamma}(\alpha )}$, 0 < α < 1, i.e. $\hat{\eta }(s)={s}^{-\alpha }$. In this case the waiting time PDF behaves as the two parameter Mittag–Leffler distribution $\psi (t)={\mathcal{L}}^{-1}\left[\frac{1}{1+{\tau }^{\alpha }{s}^{\alpha }}\right]=\frac{1}{\tau }{\left(\frac{t}{\tau }\right)}^{\alpha -1}{E}_{\alpha ,\alpha }\left(-{\left[\frac{t}{\tau }\right]}^{\alpha }\right)$ [52], while the MSD in the absence of a sink has a power law dependence on time, $\left\langle {x}^{2}(t)\right\rangle =2\mathcal{D}\frac{{t}^{\alpha }}{{\Gamma}(1+\alpha )}$. This gives an anomalous diffusive process [14], i.e. subdiffusion, since 0 < α < 1, which has been observed in many physical and biological systems, such as in charge carrier motion in amorphous semiconductors [53], anomalous transport in biological cells [54], subdiffusion in artificially crowded systems [55], etc.

For the FATD in the Laplace space we yield

$$\begin{equation}{\hat{\wp }}_{\text{fa}}(s)={\mathrm{e}}^{-\frac{{s}^{\alpha /2}}{\sqrt{\mathcal{D}}}\vert X-{x}_{0}\vert }={H}_{0,1}^{1,0}\left[\begin{matrix}\hfill \left.\frac{{s}^{\alpha /2}}{\sqrt{\mathcal{D}}}\vert X-{x}_{0}\vert \right\vert \begin{matrix}\hfill -\hfill \\ \hfill (0,1)\hfill \end{matrix}\hfill \end{matrix}\right],\end{equation} \tag{ 20 }$$

where ${H}_{p,q}^{m,n}(z)$ is the Fox H-function (A.5). The above equals sign holds by using relation (A.16). Taking the inverse Laplace transform (A.11) for the FATD gives

$$\begin{equation}{\wp }_{\text{fa}}(t)=\frac{1}{t}{H}_{1,1}^{1,0}\left[\begin{matrix}\hfill \left.\frac{\vert X-{x}_{0}\vert }{\sqrt{\mathcal{D}{t}^{\alpha }}}\right\vert \begin{matrix}\hfill (0,\alpha /2)\hfill \\ \hfill (0,1)\hfill \end{matrix}\hfill \end{matrix}\right],\end{equation} \tag{ 21 }$$

with 0 < α < 1. This is the case of anomalous diffusive search with waiting time PDF, which scales as ψ(t) ∼ t^−1−α . Graphical representation of the FATD (21) is given in figure 2, showing its long time behavior of form ℘_fa(t) ∼ t^−α/2−1. We note that the FATD has the same behavior as the first passage time density (see equation (54) in [35]), since the particle does not perform jump-like motion and thus the FATD and the first passage time density coincides. The reliability equals one, while the efficiency is

$$\begin{equation}\mathcal{E}={\int }_{0}^{\infty }{\mathrm{e}}^{-\frac{{s}^{\alpha /2}}{\sqrt{\mathcal{D}}}\vert X-{x}_{0}\vert }\mathrm{d}s=\frac{{\mathcal{D}}^{1/\alpha }}{\vert X-{x}_{0}{\vert }^{2/\alpha }}{\Gamma}\left(\frac{2+\alpha }{\alpha }\right).\end{equation} \tag{ 22 }$$

We note that the case α = 1/2 can be related to the Brownian backbone search on a two-dimensional comb [56], where the efficiency behaves as $\mathcal{E}\sim \frac{1}{\vert X-{x}_{0}{\vert }^{4}}$, while for α = 1/4 to the Brownian backbone search on a three-dimensional comb with $\mathcal{E}\sim \frac{1}{\vert X-{x}_{0}{\vert }^{8}}$ [57]. The result of one dimension Brownian search is recovered when α = 1. Since the MSD of the searcher in the absence of a sink has a power law dependence on time, we call it a subdiffusive search. We also note that from the definition $\langle t\rangle =-{\left.\frac{\partial {\hat{\wp }}_{\text{fa}}(s)}{\partial s}\right\vert }_{s=0}$ and equation (20), one finds that the mean time the particle to hit the target diverges, since some of the particles move in the wrong direction and never hit the target.

Zoom In Zoom Out Reset image size

The left panel of figure 3 shows the efficiency as a function of the initial searcher–target distance |X − x₀| for fixed exponent α. For close initial searcher–target distance |X − x₀|, subdiffusion search with small exponent α is the best search strategy. However, with increasing |X − x₀|, subdiffusive search with larger exponents α becomes more efficient, thus for large initial searcher–target distance the Brownian search is the most efficient search strategy. This behavior is consistent with our intuition since for short initial separation, Brownian search may start the search in the opposite direction of the target, while the subdiffusive search due to the waiting time in the traps reduces the probability that searcher moves away from the target in the wrong direction in a short time of the search. For longer initial searcher–target distances, the Brownian search is more efficient since with subdiffusive search the particle moves slower due to the waiting time in the traps.

Zoom In Zoom Out Reset image size

In addition to initial searcher-trap distance |X – x₀|, the search efficiency has a close relation with α. As [29] does, we also introduce the relative efficiency to discuss the dependence of search efficiency on α for fixed |X – x₀|. The relative efficiency is defined as

$$\begin{equation}{\mathcal{E}}_{\text{rel}}=\frac{\mathcal{E}(\alpha )}{\mathcal{E}({\alpha }_{\text{opt}})},\end{equation} \tag{ 23 }$$

where $\mathcal{E}({\alpha }_{\text{opt}})$ stands for the maximum efficiency for this initial separation and the corresponding value of power-law memory kernel exponent is α_opt. Obviously, ${\mathcal{E}}_{\text{rel}}=1$ if α = α_opt. As shown in the right panel of figure 3, for close initial separation the functional form of ${\mathcal{E}}_{\text{rel}}$ is non-monotonic. It drops with increasing α until it reaches the minimum, and then grows with α. However for the longer initial distance between the searcher and target |X – x₀|, the functional form of ${\mathcal{E}}_{\text{rel}}$ is monotonic. The search efficiency attains the maximum when α = 1, which implies the Brownian search is the most efficient among subdiffusive search for longer initial distance between the searcher and the target.

The efficiency as the reciprocal of the most probable value of first arrival time for the case of power-law memory kernel cannot be obtained analytically and requires numerical calculations. In figure 4, we give a comparison between the calculated efficiency (22) and the reciprocal of the most probable value of first arrival time.

Zoom In Zoom Out Reset image size

2.3. Exponential memory kernel

We consider the exponential memory kernel η(t) = e^−rt , i.e. $\hat{\eta }(s)=\frac{1}{s+r}$. Without the sink term, we recover the Brownian motion with exponential resetting to the initial position of the searcher (see for example [58]), which can be rewritten in the form [59]

$$\begin{equation}\frac{\partial }{\partial t}{P}_{r}(x,t)=\mathcal{D}\frac{{\partial }^{2}}{\partial {x}^{2}}{P}_{r}(x,t)-r\,{P}_{r}(x,t)+r\,\delta (x-{x}_{0}),\end{equation} \tag{ 24 }$$

or in the renewal form [60]

$$\begin{equation}{P}_{r}(x,t)={\mathrm{e}}^{-rt}P(x,t)+{\int }_{0}^{t}r\,{\mathrm{e}}^{-r{t}^{\prime }}P(x,{t}^{\prime })\,\mathrm{d}{t}^{\prime }.\end{equation} \tag{ 25 }$$

Here, P(x, t) is the PDF for r = 0. We note that due to the resetting mechanism, the condition ${\mathrm{lim}}_{s\to 0}\frac{1}{\hat{\eta }(s)}=0$ is no longer valid in this case.

Taking $\hat{\eta }(s)=\frac{1}{s+r}$ into FATD, we find

$$\begin{equation}{\hat{\wp }}_{\text{fa}}(s)={\mathrm{e}}^{-\frac{{(s+r)}^{1/2}}{\sqrt{\mathcal{D}}}\vert X-{x}_{0}\vert },\end{equation} \tag{ 26 }$$

from which, by inverse Laplace transform, we find the exponentially truncated Lévy–Smirnov distribution

$$\begin{equation}{\wp }_{\text{fa}}(t)=\frac{\vert X-{x}_{0}\vert }{\sqrt{4\pi \mathcal{D}{t}^{3}}}{\mathrm{e}}^{-\frac{{(X-{x}_{0})}^{2}}{4\mathcal{D}t}-rt},\end{equation} \tag{ 27 }$$

shown in figure 5. The reliability then becomes $\mathcal{P}={\mathrm{e}}^{-\sqrt{\frac{r}{\mathcal{D}}}\vert X-{x}_{0}\vert }$, which equals unity only for r = 0. Therefore, the resetting mechanism which brings the particle back to the initial position decreases the cumulative arrival probability for the particles to reach the target. For the efficiency, we find

$$\begin{equation}\mathcal{E}=\frac{2\sqrt{\mathcal{D}}(\sqrt{\mathcal{D}}+\sqrt{r}\vert X-{x}_{0}\vert )}{{(X-{x}_{0})}^{2}}{\mathrm{e}}^{-\sqrt{\frac{r}{\mathcal{D}}}\vert X-{x}_{0}\vert }.\end{equation} \tag{ 28 }$$

This result reduces to the efficiency (19) of Brownian search if r = 0. For a fixed initial distance between the searcher and the target, the search efficiency decreases with increasing reset rate r. Significantly different from the Brownian search strategy, when the memory kernel behaves as an exponential distribution, we conclude that the mean time of the particle to hit the target becomes finite, i.e. $\langle t\rangle =\frac{\vert X-{x}_{0}\vert }{\sqrt{4\mathcal{D}r}}{\mathrm{e}}^{-\sqrt{\frac{r}{\mathcal{D}}}\vert X-{x}_{0}\vert }$. This is due to the resetting mechanism which prevents the particles from moving in the wrong direction and never finding the target. Here, we note that in terms of the most probable value of first arrival time, for efficiency we find ${\mathcal{E}}_{\ast }=\frac{4r\sqrt{\mathcal{D}}}{2\sqrt{4r{(X-{x}_{0})}^{2}+9\mathcal{D}}-3\sqrt{\mathcal{D}}}$.

Zoom In Zoom Out Reset image size

2.4. Truncated power-law memory kernel

For the truncated (or tempered) power-law memory kernel $\eta (t)={\mathrm{e}}^{-rt}\frac{{t}^{\alpha -1}}{{\Gamma}(\alpha )}$, 0 < α < 1, r > 0, i.e. $\hat{\eta }(s)={(s+r)}^{-\alpha }$, in the absence of sink, we recover the subdiffusive search process with exponential resetting to the initial position. The corresponding diffusion equation for the PDF P_r (x, t) can be written in renewal form (25) through the PDF P(x, t) in the absence of resetting mechanism (r = 0) [61], as well as in the form of the tempered fractional Fokker–Planck equation [62, 63]

$$\begin{align}\hfill & {\int }_{0}^{t}{\mathrm{e}}^{-r(t-{t}^{\prime })}\frac{{(t-{t}^{\prime })}^{-\alpha }}{{\Gamma}(1-\alpha )}\frac{\partial }{\partial {t}^{\prime }}{P}_{r}(x,{t}^{\prime })\,\mathrm{d}{t}^{\prime }\hfill \\ \hfill & \hspace{25.0pt}=\mathcal{D}\frac{{\partial }^{2}}{\partial {x}^{2}}{P}_{r}(x,t)-{\int }_{0}^{t}{\mathrm{e}}^{-r(t-{t}^{\prime })}\frac{{(t-{t}^{\prime })}^{-\alpha }}{{\Gamma}(1-\alpha )}\left[r\,{P}_{r}(x,t)-r\,\delta (x-{x}_{0})\right]\mathrm{d}{t}^{\prime }.\hfill \end{align} \tag{ 29 }$$

For the FATD in Laplace space, we have

$$\begin{equation}{\hat{\wp }}_{\text{fa}}(s)={\mathrm{e}}^{-\frac{{(s+r)}^{\alpha /2}}{\sqrt{\mathcal{D}}}\vert X-{x}_{0}\vert }={H}_{0,1}^{1,0}\left[\begin{matrix}\hfill \left.\frac{{(s+r)}^{\alpha /2}}{\sqrt{\mathcal{D}}}\vert X-{x}_{0}\vert \right\vert \begin{matrix}\hfill -\hfill \\ \hfill (0,1)\hfill \end{matrix}\hfill \end{matrix}\right],\end{equation} \tag{ 30 }$$

which by inverse Laplace transform becomes

$$\begin{equation}{\wp }_{\text{fa}}(t)=\frac{{\mathrm{e}}^{-rt}}{t}{H}_{1,1}^{1,0}\left[\begin{matrix}\hfill \left.\frac{\vert X-{x}_{0}\vert }{\sqrt{\mathcal{D}{t}^{\alpha }}}\right\vert \begin{matrix}\hfill (0,\alpha /2)\hfill \\ \hfill (0,1)\hfill \end{matrix}\hfill \end{matrix}\right].\end{equation} \tag{ 31 }$$

The FATD in the long time limit is of exponentially truncated power-law form ℘_fa(t) ∼ t^−α/2−1e^−rt , which can be seen in figure 6.

Zoom In Zoom Out Reset image size

The reliability equals $\mathcal{P}={\mathrm{e}}^{-\frac{{r}^{\alpha /2}}{\sqrt{\mathcal{D}}}\vert X-{x}_{0}\vert }$, which equals unity again only for r = 0 since the resetting decreases the cumulative arrival probability for the particles to reach the target. For the efficiency we find

$$\begin{align}\hfill \mathcal{E}& ={\int }_{0}^{\infty }{\mathrm{e}}^{-\frac{{(s+r)}^{\alpha /2}}{\sqrt{\mathcal{D}}}\vert X-{x}_{0}\vert }\,\mathrm{d}s=-\frac{2}{\alpha }\frac{{\mathcal{D}}^{1/\alpha }}{\vert X-{x}_{0}{\vert }^{2/\alpha }}{\left.{\Gamma}\left(\frac{2}{\alpha },\frac{{(s+r)}^{\alpha /2}\vert X-{x}_{0}\vert }{\sqrt{\mathcal{D}}}\right)\right\vert }_{s=0}^{s=\infty }\hfill \\ \hfill & =\frac{2}{\alpha }\frac{{\mathcal{D}}^{1/\alpha }}{\vert X-{x}_{0}{\vert }^{2/\alpha }}{\Gamma}\left(\frac{2}{\alpha },\frac{{r}^{\alpha /2}\vert X-{x}_{0}\vert }{\sqrt{\mathcal{D}}}\right),\hfill \end{align} \tag{ 32 }$$

where ${\Gamma}(a,z)={\int }_{z}^{\infty }{y}^{a-1}{\mathrm{e}}^{-y}\,\mathrm{d}y$ is the upper incomplete gamma function. This result reduces to the efficiency (22) of subdiffusive random search when r = 0. From figure 7 we see that the subdiffusion search with truncated power-law memory kernel with smaller exponents α and r is more efficient. Furthermore, the mean time for the particle to hit the target becomes finite as well, i.e. $\langle t\rangle =\frac{\alpha \vert X-{x}_{0}\vert }{\sqrt{4\mathcal{D}{r}^{2-\alpha }}}{\mathrm{e}}^{-\sqrt{\frac{{r}^{\alpha }}{\mathcal{D}}}\vert X-{x}_{0}\vert }$, since the resetting prevents the particle from moving far away from the target in the wrong direction.

Zoom In Zoom Out Reset image size

In figure 8 as in the case of power-law memory kernel, we give a comparison between the efficiency calculated with (32) and the efficiency given with the reciprocal of the most probable value of first arrival time for the truncated power-law memory kernel. Here again as it is with the power-law memory kernel, the first arrival time cannot be obtained analytically and the time is calculated using numerical methods.

Zoom In Zoom Out Reset image size

We may also consider the first arrival process of Lévy flight search in the presence of the target located in position X. The particle will be removed from the location X once it arrives for the first time. The corresponding Fokker–Planck equation contains a δ-sink term with strength ℘_fa(t),

$$\begin{equation}\frac{\partial }{\partial t}P(x,t)={\mathcal{D}}_{\beta }\frac{\partial }{\partial t}{\int }_{0}^{t}\eta (t-{t}^{\prime })\frac{{\partial }^{\beta }}{\partial \vert x{\vert }^{\beta }}P(x,{t}^{\prime })\,\mathrm{d}{t}^{\prime }-{\wp }_{\text{fa}}(t)\delta (x-X),\end{equation} \tag{ 33 }$$

where $\frac{{\partial }^{\beta }}{\partial \vert x{\vert }^{\beta }}$ is the Riesz fractional derivative of order 0 < β ⩽ 2. It is defined as a pseudo-differential operator with the Fourier symbol |k|^β , i.e. [64]

$$\begin{equation}\mathcal{F}\left\{\frac{{\partial }^{\beta }}{\partial \vert x{\vert }^{\beta }}f(x)\right\}(k)=-\vert k{\vert }^{\beta }\tilde{f}(k).\end{equation} \tag{ 34 }$$

In the following discussion, we only consider 1 < β ⩽ 2 since the search reliability equals zero for one-dimension Lévy flight search without a bias, which is explained in [29, 46]. In addition, ${\mathcal{D}}_{\beta }$ represents the generalized diffusion coefficient [14]. For β = 2, it corresponds to a Brownian search with trapping. Integrating both sides of (33) over x from −∞ to ∞, the same formula for the FATD is obtained as the one given by equation (3).

Equation (33) without a sink gives the generalized, nonlocal in space and time, diffusion equation

$$\begin{equation}\frac{\partial }{\partial t}P(x,t)={\mathcal{D}}_{\beta }\frac{\partial }{\partial t}{\int }_{0}^{t}\eta (t-{t}^{\prime })\frac{{\partial }^{\beta }}{\partial \vert x{\vert }^{\beta }}P(x,{t}^{\prime })\,\mathrm{d}{t}^{\prime }.\end{equation} \tag{ 35 }$$

This equation can be obtained from CTRW with the same waiting time PDF as we consider before, and jump length PDF with infinite variance, i.e. $\tilde{\lambda }(k)\sim 1-\frac{{\sigma }^{\beta }}{2}\vert k{\vert }^{\beta }$, 0 < β < 2, giving the generalized diffusion coefficient ${\mathcal{D}}_{\beta }={\sigma }^{\beta }/[2{\tau }_{\eta }]$. Therefore, equation (35) combines Lévy flights (represented by the Riesz space fractional derivative) with subdiffusive dynamics (represented by the memory kernel).

By Fourier–Laplace transformation we have

$$\begin{equation}s\tilde{\hat{P}}(k,s)-{\mathrm{e}}^{{\imath}k{x}_{0}}=-{\mathcal{D}}_{\beta }\vert k{\vert }^{\beta }s\hat{\eta }(s)\tilde{\hat{P}}(k,s)-{\hat{\wp }}_{\text{fa}}(s){\mathrm{e}}^{{\imath}kX},\end{equation} \tag{ 36 }$$

from which we find

$$\begin{align}\hfill \tilde{\hat{P}}(k,s)& =\frac{{\mathrm{e}}^{{\imath}k{x}_{0}}-{\mathrm{e}}^{{\imath}kX}{\hat{\wp }}_{\text{fa}}(s)}{s}\frac{1}{1+{\mathcal{D}}_{\beta }\vert k{\vert }^{\beta }\hat{\eta }(s)}\hfill \\ \hfill & =\frac{{\mathrm{e}}^{{\imath}k{x}_{0}}-{\mathrm{e}}^{{\imath}kX}{\hat{\wp }}_{\text{fa}}(s)}{s}{H}_{1,1}^{1,1}\left[{\mathcal{D}}_{\beta }\hat{\eta }(s)\vert k{\vert }^{\beta }\left\vert \begin{array}{l}(0,1)\hfill \\ (0,1)\hfill \end{array}\right.\right],\hfill \end{align} \tag{ 37 }$$

where we use relation (A.17). From the condition P(x = X, s) = 0, by inverse Fourier transform of equation (37), for the FATD we find

$$\begin{equation}{\hat{\wp }}_{\text{fa}}(s)=\frac{{\int }_{0}^{\infty }\mathrm{cos}\left(k\vert X-{x}_{0}\vert \right){H}_{1,1}^{1,1}\left[\left.{\mathcal{D}}_{\beta }\hat{\eta }(s){k}^{\beta }\right.\left\vert \begin{array}{l}(0,1)\hfill \\ (0,1)\hfill \end{array}\right.\right]\mathrm{d}k}{{\int }_{0}^{\infty }{H}_{1,1}^{1,1}\left[\left.{\mathcal{D}}_{\beta }\hat{\eta }(s){k}^{\beta }\right.\left\vert \begin{array}{l}(0,1)\hfill \\ (0,1)\hfill \end{array}\right.\right]\mathrm{d}k}.\end{equation} \tag{ 38 }$$

It can be simplified as

$$\begin{align}\hfill {\hat{\wp }}_{\text{fa}}(s)& =\frac{\mathrm{sin}(\pi /\beta ){[{\mathcal{D}}_{\beta }\hat{\eta }(s)]}^{1/\beta }}{\vert X-{x}_{0}\vert }{H}_{2,3}^{2,1}\left[\left.\frac{\vert X-{x}_{0}\vert }{{[{\mathcal{D}}_{\beta }\hat{\eta }(s)]}^{1/\beta }}\right.\left\vert \begin{array}{l}(1,1/\beta ),(1,1/2)\hfill \\ (1,1),(1,1/\beta ),(1,1/2)\hfill \end{array}\right.\right]\hfill \\ \hfill & =\mathrm{sin}(\pi /\beta ){H}_{2,3}^{2,1}\left[\left.\frac{\vert X-{x}_{0}\vert }{{[{\mathcal{D}}_{\beta }\hat{\eta }(s)]}^{1/\beta }}\right.\left\vert \begin{array}{l}(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \\ (0,1),(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \end{array}\right.\right],\hfill \end{align} \tag{ 39 }$$

where we used equation (A.7), or

$$\begin{align}\hfill {\hat{\wp }}_{\text{fa}}(s)& =\frac{\mathrm{sin}(\pi /\beta ){[{\mathcal{D}}_{\beta }\hat{\eta }(s)]}^{1/\beta }}{\sqrt{\pi }\vert X-{x}_{0}\vert }{H}_{1,3}^{2,1}\left[\left.\frac{\vert X-{x}_{0}\vert }{2{[{\mathcal{D}}_{\beta }\hat{\eta }(s)]}^{1/\beta }}\right.\left\vert \begin{array}{l}(1,1/\beta )\hfill \\ (1/2,1/2),(1,1/\beta ),(1,1/2)\hfill \end{array}\!\right.\!\right]\hfill \\ \hfill & =\frac{\mathrm{sin}(\pi /\beta )}{2\sqrt{\pi }}{H}_{1,3}^{2,1}\left[\left.\frac{\vert X-{x}_{0}\vert }{2{[{\mathcal{D}}_{\beta }\hat{\eta }(s)]}^{1/\beta }}\right.\left\vert \begin{array}{l}(1-1/\beta ,1/\beta )\hfill \\ (0,1/2),(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \end{array}\right.\right],\hfill \end{align} \tag{ 40 }$$

where we used equation (A.8). For β = 2 this result turns to ${\hat{\wp }}_{\text{fa}}(s)={\mathrm{e}}^{-\frac{\vert X-{x}_{0}\vert }{\sqrt{\mathcal{D}\hat{\eta }(s)}}}$, as it should be. Notably, Lévy flight search diminishes the oversampling due to the scale-free of jump length. However, the leapovers of Lévy flight caused by the extremely long jumps may make the searcher overshoot the target. Since the length of leapovers is wider than the jump length distribution, rendering the first passage of Lévy flight, which behaves as t^−3/2 different from the first arrival process.

We note that in the absence of a sink the solution of equation (35) in Laplace space is given by (x₀ = 0)

$$\begin{align}\hfill \hat{P}(x,s)& ={\mathcal{F}}^{-1}\left[\frac{1}{s}{H}_{1,1}^{1,1}\left[{\mathcal{D}}_{\beta }\hat{\eta }(s)\vert k{\vert }^{\beta }\left\vert \begin{array}{c}\hfill (0,1)\hfill \\ \hfill (0,1)\hfill \end{array}\right.\right]\right]\hfill \\ \hfill & =\frac{{s}^{-1}}{\vert x\vert }{H}_{2,3}^{2,1}\left[\left.\frac{\vert k{\vert }^{\beta }}{{\mathcal{D}}_{\beta }\hat{\eta }(s)}\right\vert \begin{array}{l}(1,1),(1,\beta /2)\hfill \\ (1,\beta ),(1,1),(1,\beta /2)\hfill \end{array}\right].\hfill \end{align} \tag{ 41 }$$

From here one can calculate the qth moment for 0 < q < β ⩽ 2, which reads

$$\begin{align}\hfill \left\langle \vert x(t){\vert }^{q}\right\rangle & ={\int }_{-\infty }^{\infty }\vert x{\vert }^{q}P(x,t)\,\mathrm{d}x\hfill \\ \hfill & =\frac{2{\mathcal{D}}_{\beta }^{q/\beta }}{\beta }\frac{{\Gamma}(1+q){\Gamma}(1+q/\beta ){\Gamma}(-q/\beta )}{{\Gamma}(-q/2){\Gamma}(1+q/2)}{\mathcal{L}}^{-1}\left[{s}^{-1}{\left(\hat{\eta }(s)\right)}^{q/\beta }\right].\hfill \end{align} \tag{ 42 }$$

Therefore, the second moment, i.e. MSD (q = 2) diverges, and thus one defines fractional moments ⟨|x(t)|^q ⟩ to characterize the transport, from which as MSD one uses ${\langle \vert x(t){\vert }^{q}\rangle }^{2/q}$. For β = 2 and q = 2 we recover the MSD given by equation (16).

3.1. Lévy flight search

The typical case of Lévy flight search is derived when η(t) = 1, i.e. $\hat{\eta }(s)=1/s$. In the absence of a sink, from (42) for the qth moment one finds $\left\langle \vert x(t){\vert }^{q}\right\rangle \sim {t}^{q/\beta }$ and thus ${\left\langle \vert x(t){\vert }^{q}\right\rangle }^{2/q}\sim {t}^{2/\beta }$. This is the case of superdiffusion (2/β > 1), which is characteristic for Lévy flights [14]. The FATD reads

$$\begin{equation}{\hat{\wp }}_{\text{fa}}(s)=\mathrm{sin}(\pi /\beta ){H}_{2,3}^{2,1}\left[\left.\frac{{s}^{1/\beta }\vert X-{x}_{0}\vert }{{\mathcal{D}}_{\beta }^{1/\beta }}\right.\left\vert \begin{array}{l}(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \\ (0,1),(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \end{array}\right.\right].\end{equation} \tag{ 43 }$$

The inverse Laplace transformation yields

$$\begin{equation}{\wp }_{\text{fa}}(t)=\mathrm{sin}(\pi /\beta )\frac{1}{t}{H}_{3,3}^{2,1}\left[\left.\frac{\vert X-{x}_{0}\vert }{{[{\mathcal{D}}_{\beta }t]}^{1/\beta }}\right.\left\vert \begin{array}{l}(1-1/\beta ,1/\beta ),(1/2,1/2),(0,1/\beta )\hfill \\ (0,1),(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \end{array}\right.\right].\end{equation} \tag{ 44 }$$

Using the asymptotic expression of H-function in [65], for the long-time limit we have

$$\begin{equation}{\hat{\wp }}_{\text{fa}}(s)\simeq 1-\mathrm{sin}\left(\frac{\pi }{\beta }\right)\mathrm{sin}\left(\frac{\pi \beta }{2}\right)\frac{\beta {\Gamma}[2-\beta ]\vert X-{x}_{0}{\vert }^{\beta -1}}{\pi (\beta -1)}{\left(\frac{s}{{\mathcal{D}}_{\beta }}\right)}^{\frac{\beta -1}{\beta }},\end{equation} \tag{ 45 }$$

from where it follows that ℘_fa(t) ∼ t^−2+1/β in the long time limit [46]; see figure 9. The searcher will always find the target after a sufficiently long time since the search reliability is one. The search efficiency is

$$\begin{align}\hfill \mathcal{E}& =\beta \,\mathrm{sin}(\pi /\beta ){\int }_{0}^{\infty }{H}_{2,3}^{2,1}\left[\left.\frac{s\vert X-{x}_{0}{\vert }^{\beta }}{{\mathcal{D}}_{\beta }}\right.\left\vert \begin{array}{l}(1-1/\beta ,1),(1/2,\beta /2)\hfill \\ (0,\beta ),(1-1/\beta ,1),(1/2,\beta /2)\hfill \end{array}\right.\right]\mathrm{d}s\hfill \\ \hfill & =\frac{{\mathcal{D}}_{\beta }\beta }{\vert X-{x}_{0}{\vert }^{\beta }}\mathrm{cos}\left(\pi \left[1-\frac{\beta }{2}\right]\right){\Gamma}[\beta ].\hfill \end{align} \tag{ 46 }$$

This type of Lévy flight search was first considered and in details analyzed in [66], and then additionally investigated in [29], as well. As shown in figure 10, the efficiency of Lévy flight search is symmetric with respect to target position X. For short initial searcher–target distances, Brownian search is a more efficient process and the functional form of relative efficiency ${\mathcal{E}}_{\text{rel}}$ is completely monotonous, see the left panel of figure 11. However, for large initial searcher–target distances, Lévy flight search with smaller power-law exponent β is dominant. This behavior is consistent with our intuition since the jump length of Lévy flight increases with decreasing power-law exponent β, leading to the searcher arrives at the target faster for large initial distances. In addition, the power-law exponent β of optimal Lévy flight search strategy shifts to left with increasing initial searcher–target distance |X – x₀|.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2. Power-law memory kernel. Lévy flight search with trapping

Similarly, taking $\eta (t)={t}^{-\alpha }/{\Gamma}[1-\alpha ]\ (\hat{\eta }(s)={s}^{-\alpha })$ for 0 < α < 1, we have the situation of Lévy flight search (due to the Riesz space fractional derivative in the equation) with trapping (due to the power-law memory kernel, which gives long tailed waiting time PDF between the jumps). In the absence of a sink, for the qth moment we have $\left\langle \vert x(t){\vert }^{q}\right\rangle \sim {t}^{\alpha q/\beta }$, i.e. ${\left\langle \vert x(t){\vert }^{q}\right\rangle }^{2/q}\sim {t}^{2\alpha /\beta }$. This is the case of typical competition between long jumps, characterized by parameter β and long waiting times, characterized by parameter α, see [14]. For $\frac{2\alpha }{\beta }\lesseqqgtr 1$ one observes either subdiffusion, normal diffusion or superdiffusion. Thus, from equation (39), we get the expression of FATD

$$\begin{equation}{\hat{\wp }}_{\text{fa}}(s)=\mathrm{sin}(\pi /\beta ){H}_{2,3}^{2,1}\left[\left.\frac{{s}^{\alpha /\beta }\vert X-{x}_{0}\vert }{{\mathcal{D}}_{\beta }^{1/\beta }}\right.\left\vert \begin{array}{l}(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \\ (0,1),(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \end{array}\right.\right].\end{equation} \tag{ 47 }$$

By inverse Laplace transformation, one obtains

$$\begin{equation}{\wp }_{\text{fa}}(t)=\mathrm{sin}(\pi /\beta )\frac{1}{t}{H}_{3,3}^{2,1}\left[\left.\frac{\vert X-{x}_{0}\vert }{{\mathcal{D}}_{\beta }^{1/\beta }{t}^{\alpha /\beta }}\right.\left\vert \begin{array}{l}(1-1/\beta ,1/\beta ),(1/2,1/2),(0,\alpha /\beta )\hfill \\ (0,1),(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \end{array}\right.\right].\end{equation} \tag{ 48 }$$

For s → 0 we have

$$\begin{equation}{\hat{\wp }}_{\text{fa}}(s)\simeq 1-\mathrm{sin}\left(\frac{\pi }{\beta }\right)\mathrm{sin}\left(\frac{\pi \beta }{2}\right)\frac{\beta {\Gamma}[2-\beta ]\vert X-{x}_{0}{\vert }^{\beta -1}}{\pi (\beta -1)}{\left(\frac{{s}^{\alpha }}{{\mathcal{D}}_{\beta }}\right)}^{\frac{\beta -1}{\beta }},\end{equation} \tag{ 49 }$$

which for the long time limit yields ℘_fa ∼ t^{−1−α+α/β} . The reliability equals one since for s → 0 from (49) we have $\mathcal{P}=1$. Moreover, for the efficiency one obtains

$$\begin{align}\hfill \mathcal{E}& =\frac{\beta }{\alpha }\mathrm{sin}(\pi /\beta )\hfill \\ \hfill & \quad \times {\int }_{0}^{\infty }{H}_{2,3}^{2,1}\left[\left.\frac{s\vert X-{x}_{0}{\vert }^{\beta /\alpha }}{{\mathcal{D}}_{\beta }^{1/\alpha }}\right.\left\vert \begin{array}{l}(1-1/\beta ,1/\alpha ),(1/2,\beta /(2\alpha ))\hfill \\ (0,\beta /\alpha ),(1-1/\beta ,1/\alpha ),(1/2,\beta /(2\alpha ))\hfill \end{array}\right.\right]\mathrm{d}s\hfill \\ \hfill & =\frac{\beta }{\alpha }\mathrm{sin}\left(\frac{\pi }{\beta }\right)\frac{{\mathcal{D}}_{\beta }^{1/\alpha }}{\vert X-{x}_{0}{\vert }^{\beta /\alpha }}\frac{\mathrm{cos}\left(\frac{\pi \beta }{2\alpha }\right){\Gamma}[\beta /\alpha ]}{\mathrm{sin}\left(\frac{\pi }{\beta }-\frac{\pi }{\alpha }\right)}.\hfill \end{align} \tag{ 50 }$$

For α = 1, we recover the efficiency of Lévy flights search (46), discussed in the previous subsection 3.1, while for β = 2 and 0 < α < 1 we recover the efficiency (22) for subdiffusive search, considered in subsection 2.2. For a fixed power-law exponent α, the power-law exponent β of maximum efficiency decreases with increasing the initial searcher–target distance, as shown in figure 12. For example for α = 1, it shifts from β ≈ 1.5 for x₀ = 10 to β ≈ 1.25 for x₀ = 100. From figure 13, we observe that for close initial searcher–target distance |X − x₀|, the optimal power-law exponent β_opt is 2, and it drops for growing initial searcher–target distance since its long jumps.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.3. Exponential memory kernel

When the memory kernel behaves as an exponential function η(t) = e^−rt , i.e. $\hat{\eta }(s)=\frac{1}{s+r}$, in the absence of the sink the corresponding equation describes Lévy flight with resetting [67],

$$\begin{equation}\frac{\partial }{\partial t}{P}_{r}(x,t)={\mathcal{D}}_{\beta }\frac{{\partial }^{\beta }}{\partial \vert x{\vert }^{\beta }}{P}_{r}(x,t)-r\,{P}_{r}(x,t)+r\,\delta (x-{x}_{0}).\end{equation} \tag{ 51 }$$

In the presence of the sink, the FATD in the Laplace space becomes

$$\begin{equation}{\hat{\wp }}_{\text{fa}}(s)=\mathrm{sin}(\pi /\beta ){H}_{2,3}^{2,1}\left[\left.\frac{{(s+r)}^{1/\beta }\vert X-{x}_{0}\vert }{{\mathcal{D}}_{\beta }^{1/\beta }}\right.\left\vert \begin{array}{l}(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \\ (0,1),(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \end{array}\right.\right],\end{equation} \tag{ 52 }$$

which by inverse Laplace transform reads

$$\begin{equation}{\wp }_{\text{fa}}(t)=\mathrm{sin}(\pi /\beta )\frac{{\mathrm{e}}^{-rt}}{t}{H}_{3,3}^{2,1}\left[\left.\frac{\vert X-{x}_{0}\vert }{{[{\mathcal{D}}_{\beta }t]}^{1/\beta }}\right.\left\vert \begin{array}{l}(1-1/\beta ,1/\beta ),(1/2,1/2),(0,1/\beta )\hfill \\ (0,1),(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \end{array}\right.\right].\end{equation} \tag{ 53 }$$

The reliability then becomes

$$\begin{equation}\mathcal{P}=\mathrm{sin}(\pi /\beta ){H}_{2,3}^{2,1}\left[\left.\frac{{r}^{1/\beta }\vert X-{x}_{0}\vert }{{\mathcal{D}}_{\beta }^{1/\beta }}\right.\left\vert \begin{array}{l}(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \\ (0,1),(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \end{array}\right.\right].\end{equation} \tag{ 54 }$$

The search reliability has a close relation with the power-law exponent β and the resetting rate r, as well as the initial searcher-trap distance |X − x₀|. For a fixed initial separation and fixed r, we view search reliability as a function of power-law exponent β. As shown in figure 14, we find the search reliability improves for growing power-law exponent β. In addition, for fixed power-law exponent β, the search reliability grows with decreasing the resetting rate r. According to the definition of the search efficiency, we get

$$\begin{align}\hfill \mathcal{E}& ={\int }_{0}^{\infty }\frac{{\wp }_{\text{fa}}(t)}{t}\mathrm{d}t=\beta \,\mathrm{sin}(\pi /\beta )\frac{{\mathcal{D}}_{\beta }^{2}}{\vert X-{x}_{0}{\vert }^{2\beta }}\hfill \\ \hfill & \quad \times {\int }_{0}^{\infty }{\mathrm{e}}^{-rt}{H}_{3,3}^{1,2}\left[\left.\frac{{\mathcal{D}}_{\beta }t}{\vert X-{x}_{0}{\vert }^{\beta }}\right.\left\vert \begin{array}{l}(1-2\beta ,\beta ),(-2+1/\beta ,1),(1/2-\beta ,\beta /2)\hfill \\ (-2+1/\beta ,1),(1/2-\beta ,\beta /2),(-1,1)\hfill \end{array}\right.\right]\mathrm{d}t\hfill \\ \hfill & =\mathrm{sin}(\pi /\beta )r{H}_{4,3}^{1,3}\left[\left.\frac{{\mathcal{D}}_{\beta }^{1/\beta }{r}^{-1/\beta }}{\vert X-{x}_{0}\vert }\right.\left\vert \begin{array}{l}(2,1/\beta ),(1,1),(1/\beta ,1/\beta ),(1/2,1/2)\hfill \\ (1/\beta ,1/\beta ),(1/2,1/2),(1,1/\beta )\hfill \end{array}\right.\right]\hfill \\ \hfill & =\mathrm{sin}(\pi /\beta )r{H}_{3,4}^{3,1}\left[\left.\frac{\vert X-{x}_{0}\vert {r}^{1/\beta }}{{\mathcal{D}}_{\beta }^{1/\beta }}\right.\left\vert \begin{array}{l}(1-1/\beta ,1/\beta ),(1/2,1/2),(0,1/\beta )\hfill \\ (-1,1/\beta ),(0,1),(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \end{array}\right.\!\right].\hfill \end{align} \tag{ 55 }$$

Zoom In Zoom Out Reset image size

3.4. Truncated power-law memory kernel

The case with a truncated power-law memory kernel $\eta (t)={\mathrm{e}}^{-rt}\frac{{t}^{\alpha -1}}{{\Gamma}(\alpha )}$, 0 < α < 1, i.e. $\hat{\eta }(s)={(s+r)}^{-\alpha }$, in the absence of a sink, corresponds to the Lévy flight motion with trapping in the presence of a resetting mechanism, governed by the Fokker–Planck equation

$$\begin{align}\hfill & {\int }_{0}^{t}{\mathrm{e}}^{-r(t-{t}^{\prime })}\frac{{(t-{t}^{\prime })}^{-\alpha }}{{\Gamma}(1-\alpha )}\frac{\partial }{\partial {t}^{\prime }}{P}_{r}(x,{t}^{\prime })\,\mathrm{d}{t}^{\prime }\hfill \\ \hfill & \hspace{25.0pt}={\mathcal{D}}_{\beta }\frac{{\partial }^{\beta }}{\partial \vert x{\vert }^{\beta }}{P}_{r}(x,t)-{\int }_{0}^{t}{\mathrm{e}}^{-r(t-{t}^{\prime })}\frac{{(t-{t}^{\prime })}^{-\alpha }}{{\Gamma}(1-\alpha )}\left[r\,{P}_{r}(x,t)-r\,\delta (x-{x}_{0})\right]\mathrm{d}{t}^{\prime }.\hfill \end{align} \tag{ 56 }$$

Substituting $\hat{\eta }(s)={(s+r)}^{-\alpha }$ into (39), the FATD is

$$\begin{equation}{\hat{\wp }}_{\text{fa}}(s)=\mathrm{sin}(\pi /\beta ){H}_{2,3}^{2,1}\left[\left.\frac{{(s+r)}^{\alpha /\beta }\vert X-{x}_{0}\vert }{{\mathcal{D}}_{\beta }^{1/\beta }}\right.\left\vert \begin{array}{l}(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \\ (0,1),(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \end{array}\right.\right],\end{equation} \tag{ 57 }$$

where 0 < α < 1, r > 0. By inverse Laplace transformation we yield

$$\begin{equation}{\wp }_{\text{fa}}(t)=\mathrm{sin}(\pi /\beta )\frac{{\mathrm{e}}^{-rt}}{t}{H}_{3,3}^{2,1}\left[\left.\frac{\vert X-{x}_{0}\vert }{{\mathcal{D}}_{\beta }^{1/\beta }{t}^{\alpha /\beta }}\right.\left\vert \begin{array}{l}(1-1/\beta ,1/\beta ),(1/2,1/2),(0,\alpha /\beta )\hfill \\ (0,1),(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \end{array}\right.\right].\end{equation} \tag{ 58 }$$

The reliability equals

$$\begin{equation}\mathcal{P}=\mathrm{sin}(\pi /\beta ){H}_{2,3}^{2,1}\left[\left.\frac{{r}^{\alpha /\beta }\vert X-{x}_{0}\vert }{{\mathcal{D}}^{1/\beta }}\right.\left\vert \begin{array}{l}(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \\ (0,1),(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \end{array}\right.\right],\end{equation} \tag{ 59 }$$

while efficiency reads

$$\begin{equation}\mathcal{E}=\mathrm{sin}(\pi /\beta )r{H}_{3,4}^{3,1}\left[\frac{\vert X-{x}_{0}\vert {r}^{\alpha /\beta }}{{\mathcal{D}}_{\beta }^{1/\beta }}\left\vert \begin{array}{l}(1-1/\beta ,1/\beta ),(1/2,1/2),(0,\alpha /\beta )\hfill \\ (-1,\alpha /\beta ),(0,1),(1-1/\beta ,1/\beta ),(1/2,1/2)\hfill \end{array}\right.\,\right],\end{equation} \tag{ 60 }$$

which are derived on account of its definition.

We discuss the search problem from the perspective of the generalized diffusion equation. First, we consider one-dimensional Brownian search with trapping, i.e. the MSD of jump length is finite, combined with some representative cases of the memory kernel, which enter the generalized waiting time PDF, we derive the concrete form of the FATD, search reliability and efficiency. Next, we consider a one-dimensional Lévy search with trapping, i.e. the jump length of a one-dimensional searcher is power-law distribution with 1 < β ⩽ 2. We also derive the general form of the FATD, search reliability and efficiency. Finally, we generalize the above method to analyze the FATD when there are multiple targets as well as combined search strategies.

For the case of one-dimensional Brownian search with trapping, the search reliability equals one, which means the searcher will always arrive at the position of the target. For a close initial distance between the searcher and the target, the efficiency of the subdiffusion search process with small memory kernel exponent is maximal, which states the fact that the subdiffusion search process with small memory kernel exponent is the best search strategy. By the same method we find that the Brownian search outperforms for longer initial search-trap distance. In addition, for long initial separation we also find that the efficiency of subdiffusion search with large memory kernel exponent is bigger than the one with small memory kernel exponent. However, this situation changes when tempering of the memory kernel takes action, i.e. the efficiency of the subdiffusion search process with small memory kernel exponent is bigger than the one with large memory kernel exponent.

For the case of one-dimensional Lévy search with trapping, the search reliability equals one. For close initial distance between the searcher and the target, Brownian search outperforms than Lévy flight search. Combined with the above conclusion, subdiffusion search process with small memory kernel exponent is the best strategy for close initial separation. For growing initial separation, the search efficiency is maximal for Lévy flight search with smaller power-law exponent β due to its long jump. In addition, there are two optimal power-law exponents when the memory kernel of the searcher is power-law with smaller exponent α. For Lévy flight search with exponential memory kernel, the search reliability is less than one. Moreover, the search reliability improves for decreasing resetting rate r and increasing power-law exponent β, respectively. By the same method, we derive the general form of the FATD both for Lévy search in presence of multiple targets and combined Lévy searches in the presence of one target, see appendices B and C.

The authors acknowledge support from the bilateral Macedonian-Chinese research project 20-6333, and 6-13, funded under the intergovernmental Macedonian-Chinese agreement. P T, T S and L K acknowledge financial support by the German Science Foundation (DFG, Grant No. ME 1535/12-1) T S was supported by the Alexander von Humboldt Foundation. T Z was supported by the Outstanding Graduates Student Project of Gansu Education Department (Grant No. 2022CXZX-051). T Z and W D was supported by the National Natural Science Foundation of China (Grant No. 12071195).

The three parameter Mittag-Leffler function is defined by [68]

$$\begin{equation}{E}_{\rho ,\mu }^{\gamma }(t)=\sum\limits _{k=0}^{\infty }\frac{{(\gamma )}_{k}}{{\Gamma}(\rho k+\mu )}\frac{{t}^{k}}{k!},\end{equation} \tag{ A.1 }$$

where ${(\gamma )}_{k}=\frac{{\Gamma}(\gamma +k)}{{\Gamma}(\gamma )}$ is the Pochhammer symbol. Its Laplace transform reads

$$\begin{equation}\mathcal{L}\left\{{t}^{\mu -1}{E}_{\rho ,\mu }^{\gamma }(v{t}^{\rho })\right\}=\frac{{s}^{\rho \gamma -\mu }}{{({s}^{\rho }-v)}^{\gamma }}.\end{equation} \tag{ A.2 }$$

For 0 < α < 2, by using the formula [69],

$$\begin{equation}{E}_{\alpha ,\beta }^{\gamma }(-z)=\frac{{z}^{-\gamma }}{{\Gamma}(\gamma )}\sum\limits _{n=0}^{\infty }\frac{{\Gamma}(\gamma +n)}{{\Gamma}(\beta -\alpha (\gamma +n))}\frac{{(-z)}^{-n}}{n!},\quad z > 1,\end{equation} \tag{ A.3 }$$

one finds the asymptotic behavior of the three parameter M–L function. For large z (z ≫ 1) one finds

$$\begin{equation}{E}_{\alpha ,\beta }^{\gamma }(-z)\sim \frac{{z}^{-\gamma }}{{\Gamma}(\beta -\alpha \gamma )}-\gamma \frac{{z}^{-(\gamma +1)}}{{\Gamma}(\beta -\alpha (\gamma +1))},\quad z\gg 1.\end{equation} \tag{ A.4 }$$

The Fox H-function is defined as the inverse Mellin transform for a set of gamma functions [70]

$$\begin{equation}{H}_{p,q}^{m,n}\left[z\left\vert \begin{array}{l}({a}_{p},{A}_{p})\hfill \\ ({b}_{q},{B}_{q})\hfill \end{array}\right.\right]={H}_{p,q}^{m,n}\left[z\left\vert \begin{array}{l}({a}_{1},{A}_{1}),\dots ,({a}_{p},{A}_{p})\hfill \\ ({b}_{1},{B}_{1}),\dots ,{b}_{q},{B}_{q}\left.\right)\hfill \end{array}\right.\right]=\frac{1}{2\pi {\imath}}{\int }_{{\Omega}}\theta (s){z}^{-s}\,\mathrm{d}s,\end{equation} \tag{ A.5 }$$

where

$$\begin{equation}\theta (s)=\frac{{\prod }_{j=1}^{m}{\Gamma}({b}_{j}+{B}_{j}s){\prod }_{j=1}^{n}{\Gamma}(1-{a}_{j}-{A}_{j}s)}{{\prod }_{j=m+1}^{q}{\Gamma}(1-{b}_{j}-{B}_{j}s){\prod }_{j=n+1}^{p}{\Gamma}({a}_{j}+{A}_{j}s)},\end{equation} \tag{ A.6 }$$

with 0 ⩽ n ⩽ p, 1 ⩽ m ⩽ q, a_i , b_j ∈ C, A_i , B_j ∈ R⁺, i = 1, ..., p, and j = 1, ..., q. The contour Ω starting at c − i∞ and ending at c + i∞ separates the poles of the function Γ(b_j + B_j s), j = 1, ..., m from those of the function Γ(1 − a_i − A_i s), i = 1, ..., n.

The Mellin-cosine transform of Fox H-function is given by [70]

$$\begin{align}\hfill & {\int }_{0}^{\infty }{\kappa }^{\rho -1}\mathrm{cos}(\kappa x){H}_{p,q}^{m,n}\left[a{\kappa }^{\delta }\left\vert \begin{array}{l}({a}_{p},{A}_{p})\hfill \\ ({b}_{q},{B}_{q})\hfill \end{array}\right.\right]d\kappa \hfill \\ \hfill & \hspace{25.0pt}=\frac{\pi }{{x}^{\rho }}{H}_{q+1,p+2}^{n+1,m}\left[\frac{{x}^{\delta }}{a}\left\vert \begin{array}{l}(1-{b}_{q},{B}_{q}),\left(\frac{1+\rho }{2},\frac{\delta }{2}\right)\hfill \\ (\rho ,\delta ),(1-{a}_{p},{A}_{p}),\left(\frac{1+\rho }{2},\frac{\delta }{2}\right)\hfill \end{array}\right.\right].\hfill \end{align} \tag{ A.7 }$$

$$\begin{align}\hfill & {\int }_{0}^{\infty }{\kappa }^{\rho -1}\mathrm{cos}(\kappa x){H}_{p,q}^{m,n}\left[a{\kappa }^{\delta }\left\vert \begin{array}{l}({a}_{p},{A}_{p})\hfill \\ ({b}_{q},{B}_{q})\hfill \end{array}\right.\right]d\kappa \hfill \\ \hfill & \hspace{25.0pt}=\frac{{2}^{\rho -1}\sqrt{\pi }}{{x}^{\rho }}{H}_{p+2,q}^{m,n+1}\left[a{\left(\frac{2}{x}\right)}^{\delta }\left\vert \begin{array}{l}\left(\frac{2-\rho }{2},\frac{\delta }{2}\right),({a}_{p},{A}_{p}),\left(\frac{1-\rho }{2},\frac{\delta }{2}\right)\hfill \\ ({b}_{q},{B}_{q}),\hfill \end{array}\right.\right]\hfill \\ \hfill & \hspace{25.0pt}=\frac{{2}^{\rho -1}\sqrt{\pi }}{{x}^{\rho }}{H}_{q,p+2}^{n+1,m}\left[\frac{1}{a}{\left(\frac{x}{2}\right)}^{\delta }\left\vert \begin{array}{l}(1-{b}_{q},{B}_{q})\hfill \\ \left(\frac{\rho }{2},\frac{\delta }{2}\right),(1-{a}_{p},{A}_{p}),\left(\frac{1+\rho }{2},\frac{\delta }{2}\right)\hfill \end{array}\right.\right].\hfill \end{align} \tag{ A.8 }$$

The Mellin transform of the Fox H-function yields

$$\begin{equation}{\int }_{0}^{\infty }{x}^{\xi -1}{H}_{p,q}^{m,n}\left[ax\left\vert \begin{array}{l}({a}_{p},{A}_{p})\hfill \\ ({b}_{q},{B}_{q})\hfill \end{array}\right.\right]\mathrm{d}x={a}^{-\xi }\theta (\xi ),\end{equation} \tag{ A.9 }$$

where θ(ξ) is defined in equation (A.6). The Laplace transform of the Fox H-function is

$$\begin{equation}\mathcal{L}\left\{{t}^{\delta }{H}_{p,q}^{m,n}\left[a{t}^{\sigma }\left\vert \begin{array}{l}({a}_{p},{A}_{p})\hfill \\ ({b}_{q},{B}_{q})\hfill \end{array}\right.\right]\right\}={s}^{-1-\delta }{H}_{p+1,q}^{m,n+1}\left[a{s}^{-\sigma }\left\vert \begin{array}{l}(-\delta ,\sigma ),({a}_{p},{A}_{p})\hfill \\ ({b}_{q},{B}_{q})\hfill \end{array}\right.\right],\end{equation} \tag{ A.10 }$$

where $\mathfrak{R}\left\{\delta +\sigma {\mathrm{min}}_{1\leqslant j\leqslant m}\left(\frac{{b}_{j}}{{B}_{j}}\right)\right\} > -1$, δ ∈ C, $s\in C(\mathfrak{R}(s) > 0)$, and a, σ, θ are positive. The parameter θ defined by $\theta ={\sum }_{j=1}^{n}\,{A}_{j}-{\sum }_{j=n+1}^{p}\,{A}_{j}+{\sum }_{j=1}^{m}\,{B}_{j}-{\sum }_{j=m+1}^{q}\,{B}_{j}$. The inverse Laplace transform of the H-function is as follows

$$\begin{equation}{\mathcal{L}}^{-1}\left\{{s}^{-\rho }{H}_{p,q}^{m,n}\left[a{s}^{\sigma }\left\vert \begin{array}{l}({a}_{p},{A}_{p})\hfill \\ ({b}_{q},{B}_{q})\hfill \end{array}\right.\right]\right\}={t}^{\rho -1}{H}_{p+1,q}^{m,n}\left[a{t}^{-\sigma }\left\vert \begin{array}{l}({a}_{p},{A}_{p}),(\rho ,\sigma ),\hfill \\ ({b}_{q},{B}_{q})\hfill \end{array}\right.\right],\end{equation} \tag{ A.11 }$$

where $\mathfrak{R}\left\{\rho +\sigma {\mathrm{max}}_{1\leqslant i\leqslant n}\left(\frac{1}{{A}_{i}}-\frac{{a}_{i}}{{A}_{i}}\right)\right\} > 0$, ρ, a, s ∈ C ($\mathfrak{R}\left\{s\right\} > 0$), $\vert \mathrm{arg}\,a\vert < \frac{1}{2}\pi (\theta -\sigma )$ and σ > 0 [70].

Fox H-function has the following properties [70]

$$\begin{equation}{H}_{p,q}^{m,n}\left[{z}^{\delta }\left\vert \begin{array}{l}({a}_{p},{A}_{p})\hfill \\ ({b}_{q},{B}_{q})\hfill \end{array}\right.\right]=\frac{1}{\delta }\cdot {H}_{p,q}^{m,n}\left[z\left\vert \begin{array}{l}({a}_{p},{A}_{p}/\delta )\hfill \\ ({b}_{q},{B}_{q}/\delta )\hfill \end{array}\right.\right],\end{equation} \tag{ A.12 }$$

$$\begin{equation}{z}^{\sigma }{H}_{p,q}^{m,n}\left[z\left\vert \begin{array}{l}({a}_{p},{A}_{p})\hfill \\ ({b}_{q},{B}_{q})\hfill \end{array}\right.\right]={H}_{p,q}^{m,n}\left[z\left\vert \begin{array}{l}({a}_{p}+\sigma {A}_{p},{A}_{p})\hfill \\ ({b}_{q}+\sigma {B}_{q},{B}_{q})\hfill \end{array}\right.\right],\end{equation} \tag{ A.13 }$$

$$\begin{align}\hfill & {H}_{p,q}^{m,n}\left[z\left\vert \begin{array}{l}({a}_{1},{A}_{1}),({a}_{2},{A}_{2}),\dots ,({a}_{p},{A}_{p})\hfill \\ ({b}_{1},{B}_{1}),\dots ,({b}_{q-1},{B}_{q-1}),({a}_{1},{A}_{1})\hfill \end{array}\right.\right]={H}_{p-1,q-1}^{m,n-1}\left[z\left\vert \begin{array}{l}({a}_{2},{A}_{2}),\dots ,({a}_{p},{A}_{p})\hfill \\ ({b}_{1},{B}_{1}),\dots ,({b}_{q-1},{B}_{q-1})\hfill \end{array}\right.\right],\hfill \end{align} \tag{ A.14 }$$

$$\begin{equation}{H}_{p,q}^{m,n}\left[z\left\vert \begin{array}{l}({a}_{p},{A}_{p})\hfill \\ ({b}_{q},{B}_{q})\hfill \end{array}\right.\right]={H}_{q,p}^{n,m}\left[\frac{1}{z}\left\vert \begin{array}{l}(1-{b}_{q},{B}_{q})\hfill \\ (1-{a}_{p},{A}_{p})\hfill \end{array}\right.\right].\end{equation} \tag{ A.15 }$$

The Fox H-function as a special cases transforms to other known functions:

$$\begin{equation}{\mathrm{e}}^{-z}={H}_{0,1}^{1,0}\left[z\left\vert \begin{array}{ll}\hfill & -\hfill \\ \hfill & (0,1)\hfill \end{array}\right.\right],\end{equation} \tag{ A.16 }$$

$$\begin{equation}\frac{1}{{(1\pm z)}^{a}}=\frac{1}{{\Gamma}(a)}{H}_{1,1}^{1,1}\left[\pm z\left\vert \begin{array}{l}(1-a,1)\hfill \\ (0,1)\hfill \end{array}\right.\right].\end{equation} \tag{ A.17 }$$

We may also consider the problem of Lévy search in the presence of multiple targets,

$$\begin{equation}\frac{\partial }{\partial t}P(x,t)={\mathcal{D}}_{\beta }\frac{\partial }{\partial t}{\int }_{0}^{t}\eta (t-{t}^{\prime })\frac{{\partial }^{\beta }}{\partial \vert x{\vert }^{\beta }}P(x,{t}^{\prime })\,\mathrm{d}{t}^{\prime }-\sum\limits _{j=1}^{N}\,{\wp }_{\text{fa},j}(t)\delta (x-{x}_{j}).\end{equation} \tag{ B.1 }$$

By Fourier–Laplace transformation we obtain

$$\begin{equation}s\tilde{\hat{P}}(k,s)-{\mathrm{e}}^{{\imath}k{x}_{0}}=-{\mathcal{D}}_{\beta }\vert k{\vert }^{\beta }s\hat{\eta }(s)\tilde{\hat{P}}(k,s)-\sum\limits _{j=1}^{N}\,{\hat{\wp }}_{\text{fa},j}(s){\mathrm{e}}^{{\imath}k{x}_{j}},\end{equation} \tag{ B.2 }$$

from where we find

$$\begin{align}\hfill \tilde{\hat{P}}(k,s)& =\frac{{\mathrm{e}}^{{\imath}k{x}_{0}}-{\sum }_{j=1}^{N}{\hat{\wp }}_{\text{fa},j}(s){\mathrm{e}}^{{\imath}k{x}_{j}}}{s+{\mathcal{D}}_{\beta }\vert k{\vert }^{\beta }s\hat{\eta }(s)}=\frac{{\mathrm{e}}^{{\imath}k{x}_{0}}-{\sum }_{j=1}^{N}{\hat{\wp }}_{\text{fa},j}(s){\mathrm{e}}^{{\imath}k{x}_{j}}}{s}\frac{1}{1+{\mathcal{D}}_{\beta }\vert k{\vert }^{\beta }\hat{\eta }(s)}\hfill \\ \hfill & =\frac{{\mathrm{e}}^{{\imath}k{x}_{0}}-{\sum }_{j=1}^{N}{\hat{\wp }}_{\text{fa},j}(s){\mathrm{e}}^{{\imath}k{x}_{j}}}{s}{H}_{1,1}^{1,1}\left[\left.{\mathcal{D}}_{\beta }\hat{\eta }(s)\vert k{\vert }^{\beta }\right.\left\vert \begin{array}{l}(0,1)\hfill \\ (0,1)\hfill \end{array}\right.\right].\hfill \end{align} \tag{ B.3 }$$

From the condition $\hat{P}(x={x}_{j},s)=0$, by taking the inverse Fourier transform to equation (B.3) we get for j = 1, 2, ..., N

$$\begin{align*}\hfill & {\int }_{-\infty }^{\infty }{\mathrm{e}}^{{\imath}k({x}_{0}-{x}_{j})}{H}_{1,1}^{1,1}\left[\left.{\mathcal{D}}_{\beta }\hat{\eta }(s)\vert k{\vert }^{\beta }\right.\left\vert \begin{array}{l}(0,1)\hfill \\ (0,1)\hfill \end{array}\right.\right]\mathrm{d}k=\sum\limits _{i=1}^{N}{\hat{\wp }}_{\text{fa},i}(s){\int }_{-\infty }^{\infty }{\mathrm{e}}^{{\imath}k({x}_{i}-{x}_{j})}{H}_{1,1}^{1,1}\left[\left.{\mathcal{D}}_{\beta }\hat{\eta }(s)\vert k{\vert }^{\beta }\right.\left\vert \begin{array}{l}(0,1)\hfill \\ (0,1)\hfill \end{array}\right.\right]\mathrm{d}k.\hfill \end{align*}$$

It can be simplified as

$$\begin{align}\hfill & {\int }_{0}^{\infty }\mathrm{cos}(k\vert {x}_{0}-{x}_{j}\vert ){H}_{1,1}^{1,1}\left[\left.{\mathcal{D}}_{\beta }\hat{\eta }(s){k}^{\beta }\right.\left\vert \begin{array}{l}(0,1)\hfill \\ (0,1)\hfill \end{array}\right.\right]\mathrm{d}k\hfill \\ \hfill & \hspace{25.0pt}=\sum\limits _{i=1}^{N}{\hat{\wp }}_{\text{fa},i}(s){\int }_{0}^{\infty }\mathrm{cos}(k\vert {x}_{i}-{x}_{j}\vert ){H}_{1,1}^{1,1}\left[\left.{\mathcal{D}}_{\beta }\hat{\eta }(s){k}^{\beta }\right.\left\vert \begin{array}{l}(0,1)\hfill \\ (0,1)\hfill \end{array}\right.\right]\mathrm{d}k.\hfill \end{align} \tag{ B.4 }$$

With the help of the Fox-H function, we have

$$\begin{align}\hfill & \frac{\sqrt{\pi }}{\vert {x}_{0}-{x}_{j}\vert }\frac{1}{\beta }{H}_{1,3}^{2,1}\left[\frac{\vert {x}_{0}-{x}_{j}\vert }{2{[{\mathcal{D}}_{\beta }\hat{\eta }(s)]}^{\frac{1}{\beta }}}\left\vert \begin{array}{l}(1,1/\beta )\hfill \\ (1/2,1/2),(1,1/\beta ),(1,1/2)\hfill \end{array}\right.\right]\hfill \\ \hfill & \hspace{25.0pt}={\hat{\wp }}_{\text{fa},j}(s)\frac{1}{\beta {[{D}_{\beta }\hat{\eta }(s)]}^{1/\beta }}\frac{\pi }{\mathrm{sin}\frac{\pi }{\beta }}+\sum\limits _{i\ne j}{\hat{\wp }}_{\text{fa},i}(s)\frac{\sqrt{\pi }}{\vert {x}_{i}-{x}_{j}\vert }\frac{1}{\beta }{H}_{1,3}^{2,1}\hfill \\ \hfill & \hspace{25.0pt}\quad \times \left[\frac{\vert {x}_{i}-{x}_{j}\vert }{2{[{\mathcal{D}}_{\beta }\hat{\eta }(s)]}^{\frac{1}{\beta }}}\left\vert \begin{array}{l}(1,1/\beta )\hfill \\ (1/2,1/2),(1,1/\beta ),(1,1/2)\hfill \end{array}\right.\right].\hfill \end{align} \tag{ B.5 }$$

Further, we assume that W(x, t) is the solution of (B.1) without the sink term and the initial position is assumed to be 0, as above we find it reads

$$\begin{equation}W(x,s)=\frac{2}{s}\frac{1}{\beta }\frac{\sqrt{\pi }}{\vert x\vert }{H}_{1,3}^{2,1}\left[\frac{\vert x\vert }{2{[{\mathcal{D}}_{\beta }\hat{\eta }(s)]}^{\frac{1}{\beta }}}\left\vert \begin{array}{l}(1,1/\beta )\hfill \\ (1/2,1/2),(1,1/\beta ),(1,1/2)\hfill \end{array}\right.\right].\end{equation} \tag{ B.6 }$$

Let us use the simplified notation W(x_i − x_j , s) = W_ij . Hence (B.5) can be simplified as

$$\begin{equation}{W}_{0j}={\hat{\wp }}_{\text{fa},j}(s){W}_{0}+\sum\limits _{i\ne j}\,{\hat{\wp }}_{\text{fa},i}(s){W}_{ij},\end{equation} \tag{ B.7 }$$

where ${W}_{0}=W(0,s)=\frac{2}{s\beta {[{D}_{\beta }\hat{\eta }(s)]}^{1/\beta }}\frac{\pi }{\mathrm{sin}(\pi /\beta )}$.

The FATD is the sum of fluxes to both targets

$$\begin{equation}{\hat{\wp }}_{\text{fa}}(s)=\sum\limits _{i=1}^{N}\,{\hat{\wp }}_{\text{fa},i}(s).\end{equation} \tag{ B.8 }$$

We can derive all the splitting FATDs ${\hat{\wp }}_{\text{fa},i}(s)$ from (B.7) since it has a unique solution for a system of N equations, so the FATD can be obtained. Since the calculation for N targets is too complicated, we only focus here on two targets, respectively, placed at x₁ and x₂. Then the FATD becomes

$$\begin{equation}{\hat{\wp }}_{\text{fa}}(s)=\frac{{W}_{01}+{W}_{02}}{{W}_{21}+{W}_{0}}=\frac{\frac{{W}_{01}}{{W}_{0}}+\frac{{W}_{02}}{{W}_{0}}}{\frac{{W}_{21}}{{W}_{0}}+1}.\end{equation} \tag{ B.9 }$$

For W_ij the limit of small s reads

$$\begin{equation}{W}_{ij}\simeq \frac{2}{s\beta {[{\mathcal{D}}_{\beta }\hat{\eta }(s)]}^{1/\beta }}\frac{\pi }{\mathrm{sin}\frac{\pi }{\beta }}+\frac{\sqrt{\pi }}{s}\frac{{\Gamma}[1/2-\beta /2]}{{\Gamma}[\beta /2]}\frac{\vert {x}_{i}-{x}_{j}{\vert }^{\beta -1}}{{2}^{\beta -1}{\mathcal{D}}_{\beta }\hat{\eta }(s)}.\end{equation} \tag{ B.10 }$$

For the long-time radio of $\frac{{W}_{ij}}{{W}_{0}}$, we have

$$\begin{equation}\frac{{W}_{ij}}{{W}_{0}}\simeq 1-{\gamma }_{\beta }\vert {x}_{i}-{x}_{j}{\vert }^{\beta -1}{[{\mathcal{D}}_{\beta }\hat{\eta }(s)]}^{1/\beta -1},\end{equation} \tag{ B.11 }$$

where ${\gamma }_{\beta }=\frac{\beta {\Gamma}[2-\beta ]}{\pi (\beta -1)}\mathrm{sin}(\frac{\pi }{\beta })\mathrm{sin}(\frac{\pi \beta }{2})$. Correspondingly,

$$\begin{equation}{\hat{\wp }}_{\text{fa}}(s)\simeq 1-\frac{{\gamma }_{\beta }{[{\mathcal{D}}_{\beta }\hat{\eta }(s)]}^{1/\beta -1}}{2}(\vert {x}_{1}-{x}_{0}{\vert }^{\beta -1}+\vert {x}_{2}-{x}_{0}{\vert }^{\beta -1}-\vert {x}_{2}-{x}_{1}{\vert }^{\beta -1}).\end{equation} \tag{ B.12 }$$

From the expression, one can see that the search reliability $P={\hat{\wp }}_{\text{fa}}(s=0)=1$ since ${\mathrm{lim}}_{s\to 0}[1/\hat{\eta }(s)]=0$.

For the special case of η(t) = 1, i.e. $\hat{\eta }(s)=1/s$, we recover the results of Brownian search and Lévy flight search for two targets considered in [28].

We may also consider the problem of combined Lévy search in the presence of one target,

$$\begin{equation}\frac{\partial }{\partial t}P(x,t)=\frac{\partial }{\partial t}{\int }_{0}^{t}\eta (t-{t}^{\prime })\sum\limits _{j=1}^{n}\,{\mathcal{D}}_{{\beta }_{j}}\frac{{\partial }^{{\beta }_{j}}}{\partial \vert x{\vert }^{{\beta }_{j}}}P(x,{t}^{\prime })\,\mathrm{d}{t}^{\prime }-{\wp }_{\text{fa}}(t)\delta (x-X).\end{equation} \tag{ C.1 }$$

By Fourier–Laplace transform we have

$$\begin{equation}s\tilde{\hat{P}}(k,s)-{\mathrm{e}}^{{\imath}k{x}_{0}}=-\sum\limits _{j=1}^{n}\,{\mathcal{D}}_{{\beta }_{j}}\vert k{\vert }^{{\beta }_{j}}s\hat{\eta }(s)\tilde{\hat{P}}(k,s)-{\hat{\wp }}_{\text{fa}}(s){\mathrm{e}}^{{\imath}kX},\end{equation} \tag{ C.2 }$$

from where it follows

$$\begin{equation}\tilde{\hat{P}}(k,s)=\frac{{\mathrm{e}}^{{\imath}k{x}_{0}}-{\hat{\wp }}_{\text{fa}}(s){\mathrm{e}}^{{\imath}kX}}{s+{\sum }_{j=1}^{n}\,{\mathcal{D}}_{{\beta }_{j}}\vert k{\vert }^{{\beta }_{j}}s\hat{\eta }(s)}=\frac{1}{s}\times \frac{{\mathrm{e}}^{{\imath}k{x}_{0}}-{\hat{\wp }}_{\text{fa}}(s){\mathrm{e}}^{{\imath}kX}}{1+{\sum }_{j=1}^{n}\,{\mathcal{D}}_{{\beta }_{j}}\hat{\eta }(s)\vert k{\vert }^{{\beta }_{j}}}.\end{equation} \tag{ C.3 }$$

Following the similar way, the FATD takes on the form

$$\begin{equation}{\hat{\wp }}_{\text{fa}}(s)=\frac{{\int }_{0}^{\infty }\frac{\mathrm{cos}\left(k\vert X-{x}_{0}\vert \right)}{\frac{1}{\hat{\eta }(s)}+{\sum }_{j=1}^{n}{\mathcal{D}}_{{\beta }_{j}}{k}^{{\beta }_{j}}}\mathrm{d}k}{{\int }_{0}^{\infty }\frac{1}{\frac{1}{\hat{\eta }(s)}+{\sum }_{j=1}^{n}{\mathcal{D}}_{{\beta }_{j}}{k}^{{\beta }_{j}}}\mathrm{d}k}.\end{equation} \tag{ C.4 }$$

The search reliability is simplified as

$$\begin{equation*}P=\frac{{\int }_{0}^{\infty }\frac{\mathrm{cos}\left(k\vert X-{x}_{0}\vert \right)}{{\sum }_{j=1}^{n}{\mathcal{D}}_{{\beta }_{j}}{k}^{{\beta }_{j}}}\mathrm{d}k}{{\int }_{0}^{\infty }\frac{1}{{\sum }_{j=1}^{n}{\mathcal{D}}_{{\beta }_{j}}{k}^{{\beta }_{j}}}\mathrm{d}k},\end{equation*}$$

since we consider ${\mathrm{lim}}_{s\to 0}[1/\hat{\eta }(s)]=0$. We recover the results of search reliability in [28] if we consider a combination of two Lévy processes, which states the search reliability is 1 if both 1 < β_j ⩽ 2. In addition, the search reliability and efficiency for combined Lévy–Brownian search were considered in [30] by the method of numerical integration when the memory kernel behaves as η(t) = 1.