key: cord-0329470-jfyng37g
authors: Khadem, S. Mohsen J.; Klapp, Sabine H.L.; Klages, Rainer
title: Search efficiency of discrete fractional Brownian motion in a random distribution of targets
date: 2020-07-15
journal: nan
DOI: 10.1103/physrevresearch.3.023169
sha: 834972069593af70868414af4e48ec876f22a05c
doc_id: 329470
cord_uid: jfyng37g

Efficiency of search for randomly distributed targets is a prominent problem in many branches of the sciences. For the stochastic process of L'evy walks, a specific range of optimal efficiencies was suggested under variation of search intrinsic and extrinsic environmental parameters. In this article, we study fractional Brownian motion as a search process, which under parameter variation generates all three basic types of diffusion, from sub- to normal to superdiffusion. In contrast to L'evy walks, fractional Brownian motion defines a Gaussian stochastic process with power law memory yielding anti-persistent, respectively persistent motion. Computer simulations of search by time-discrete fractional Brownian motion in a uniformly random distribution of targets show that maximising search efficiencies sensitively depends on the definition of efficiency, the variation of both intrinsic and extrinsic parameters, the perception of targets, the type of targets, whether to detect only one or many of them, and the choice of boundary conditions. In our simulations we find that different search scenarios favour different modes of motion for optimising search success, defying a universality across all search situations. Some of our numerical results are explained by a simple analytical model. Having demonstrated that search by fractional Brownian motion is a truly complex process, we propose an over-arching conceptual framework based on classifying different search scenarios. This approach incorporates search optimisation by L'evy walks as a special case.

Finding randomly located objects is a challenge for every human being, be it the search for mushrooms [1] , for lost keys [2] , or for food [3] . Within the context of modern society, attempts to solve this problem fuelled the development of operations research, which aims at optimising tasks such as search games [4] , locating submarines [5] , or human rescue missions [6] . For biological organisms, the successful location of food sources is crucial for their survival, as is addressed within movement ecology [7] [8] [9] . On a theoretical level, biological foraging processes are typically modelled in terms of stochastic dynamics [2, 8, 10 ]. An important paradigm for this modeling was put forward by Karl Pearson, who suggested at the beginning of the last century that organisms may migrate according to simple random walks [11] characterised by Gaussian position distributions in a suitable scaling limit. This paradigm was challenged two decades ago by the experimental observation [12] and a corresponding theory [13] that wandering albatrosses searching for food performed flights according to non-Gaussian step length distributions [14] . In this case, the mean square displacement (MSD) of an ensemble of moving agents may not grow linearly in time like for Gaussian spreading generated by random walks or Brownian motion. Instead, it may grow nonlinearly, x 2 (t) ∼ t α with α = 1, where x(t) is the position of an agent in space at time t. This phenomenon is known as anomalous diffusion, where α > 1 denotes superdiffusion and α < 1 subdiffusion while α = 1 refers to normal diffusion [15] [16] [17] [18] [19] .

Motivated by these developments, much recent research was devoted to explore the relevance of more nontrivial diffusion processes for modeling foraging. Inspired by Refs. [12, 13] , the focus was on superdiffusive Lévy walks (LWs) determined by power-law step length distributions [10] . Along similar lines normal diffusive intermittent motion [2] and correlated random walks [8, 20] have been analysed. Nevertheless, an over-arching framework for stochastic search in movement ecology is still missing. Especially for the wide variety of anomalous stochastic processes beyond LWs [19] , efficiency of search has not very much been investigated. This applies particularly to fractional Brownian motion (fBm), a paradigmatic stochastic process that, in contrast to Lévy dynamics, exhibits Gaussian distributions and power-law memory by generating the whole spectrum of anomalous diffusion [21] . The goal of our article is to explore in a very systematic way, based on extensive computer simulations and simple analytical arguments, the complexity of search exhibited by fBm. We hope that our work will set the scene for further studies to understand biological foraging on the basis of stochastic theory.

Already at the beginning of the 90's it was suggested that superdiffusive motion obtained from Lévy flights and walks may optimise the search for targets by increasing the search efficiency compared to Brownian motion [22] .

Lévy flights can be modelled as a Markov process, where the instantaneous jumps over certain distances are sampled randomly from a power law distribution [17, 23] . In the special case of LWs, any jump length is coupled with the time to perform the jump by assuming a constant speed [24] . In order to explain the experimental albatross data of Ref. [12] , Ref. [13] proposed a simple two-dimensional stochastic search model. It consists of a Lévy walker searching for targets, which are disks uniformly random distributed in the plane. For a sparse field of replenishing, immobile targets a suitably defined search efficiency yielded a maximum for power law jumps, while Brownian motion was optimal when the density of the target distribution increased [13, 25, 26] . This result became famous as the Lévy flight foraging hypothesis (LFFH) [10, 13] . The LFFH initiated a long debate, particularly when applied to finding many targets under biologically realistic search conditions [2, 8-10, 14, 24, 27-41] . A special case of the LFFH is single target search in simplified theoretical settings, which defines mathematically solvable first passage (FP) and first arrival (FA) problems [8, [42] [43] [44] [45] . 'Passage' corresponds to the situation of a biological cruise forager, who can perceive a target while moving. Thus, a target is found whenever a cruise searcher 'passes' it. A saltatory forager, on the other hand, does not scan for a target while moving but has to land, respectively to 'arrive' on it, or within a suitable neighborhood of it, to perceive it after performing a jump [2, 33, 44, 46] . Examples of cruise foragers are large fish, such as tuna and sharks. Saltatory foraging was observed among smaller fish, ground foraging birds, lizards and insects [47] .

The LFFH created awareness that apart from models related to classical random walks and Brownian motion, which yield normal diffusion with Gaussian distributions, more advanced stochastic processes are available, and needed, in order to understand data of foraging organisms (see, e.g., Refs. [8, 9, 24, 27, [34] [35] [36] and respective discussions). Most notably, it motivated the experimental biological foraging community to look for power laws in data. Over the past two decades many analyses of foraging data indeed suggested the existence of dynamics governed by power laws [10, [28] [29] [30] [31] 36] . At the same time, however, evidence accumulated that in many cases questionable data analyses were performed by checking for Lévy dynamics, partially due to a lack of full appreciation of the relevant theoretical background [2, 8, 9, 14, 24, 27, 34, 37, 40, 46, 48] . Another fundamental problem was missing knowledge whether the observed movements are intrinsic or extrinsic to a forager [7, 8] , being induced, say, by the food source distribution or other environmental conditions [9, 12, 27, 29, 46, 48, 49] , which by themselves could be governed by power laws. The complex interaction between forager and environment during search defines a topic of more recent research [20, 26, 27, 38, 41, [50] [51] [52] [53] [54] [55] . A question related to this is to which extent the optimality of a specific search strategy as predicted by a given model is robust with respect to varying parameters of the process and the given environment. Further, it remains unclear whether this optimality depends on other details of the search scenario. Along these lines the original theoretical results leading to the LFFH [13, 25, 26] were critically investigated by limiting their range of application [2, 8, 27, 33, 38, 46, 48, [56] [57] [58] . In very recent work they were eventually largely refuted [37] ; however, see Refs. [39, 40] for an associated reply and comment, and also Ref. [38] . We will come back to these important points at the end of our work in the light of our new results.

Apart from using simple random walks, classical Brownian motion or Lévy dynamics for modelling foraging, other theoretical studies considered intermittent motion [2, 59] , correlated random walks [20, 60, 61] , multiscale random walks [50, 57, 58] , more non-trivial onedimensional motion [62] and generalisations of the original [13] LW search model [51, 53] . Yet beyond these models, there exist numerous other types of, in particular, anomalous stochastic processes [15-17, 19, 63] which, to our knowledge, have not been assessed for optimising search. One of the most famous and important classes is fractional Brownian motion (fBm), which was originally studied by Kolmogorov [64] and Yaglom [65] and has become more widely known through the work of Mandelbrot and Van Ness [21] . FBm defines a Gaussian stochastic process generating anomalous diffusion [43, [66] [67] [68] [69] . However, while for non-Gaussian Lévy dynamics anomalous diffusion originates from sampling power law distributions, in the case of fBm it is generated by a power law correlation function decay in time yielding non-Markovian dynamics. This enables fBm to produce the whole spectrum of anomalous diffusion under parameter variation, from sub-to normal, to superdiffusion. LWs, on the other hand, are purely superdiffusive. Hence these two stochastic processes represent fundamentally different classes of anomalous dynamics. FBm has been widely used to describe the experimentally observed anomalous diffusion of a tracer particle in a visco-elastic system [70] , in artificially crowded environments [71] [72] [73] , in complex intracellular media [74, 75] and in living cells [76] ; for a wide range of further applications of fBm see Ref. [77] . On the other hand, so far fBm has been rarely used to understand the search by biological organisms. Motivated by insect motion and other problems in ecology, the efficiency of search of fBm and LWs has been compared in computer simulations [78] . In Ref. [32] , non-trivial autocorrelation function decay was observed in experimental data of bumblebee flights, with the conclusion that correlation decay is important to characterise biological search. While both exponential and power law decay of correlations has been reported for biological cell migration [79] , extracting correlation functions from experimental data has not been very prominent in the foraging literature. By studying search in fBm we thus further advocate the use of correlation functions, or memory, for understanding biological foraging. We may also remark that very recently, subd-iffusive search related to fBm has even been studied on the human proteomic network, as a possible explanation that via protein expression the COVID-19 virus only attacks a certain subset of organs [80] . FP problems of fBm on both finite and infinite domains have been analysed rigorously mathematically [81] , as well as in relation to physical and biological applications [43, 77, [82] [83] [84] [85] . However, apart from Ref. [78] , more complex search scenarios have not been investigated for fBm.

In this article we compute the efficiency of a searcher moving by fBm in a two-dimensional array of uniformly random distributed targets. The target distribution is determined by the number density and the radius of the circular targets as two extrinsic parameters. In addition, there are two parameters that are intrinsic to fBm, the Hurst exponent determining the type of diffusion generated by fBm and the mean jump length of the process characterising the strength of the diffusive spreading. Under these conditions, we consider both cruise and saltaltory foragers (viz. FP and FA problems) for finding targets. We compare these two different settings by using two generically different types of efficiencies. The first one employs an ensemble average of moving fBm particles for finding one target only, while the second one is given by the time-ensemble average of searchers for finding many targets along their whole trajectories [13, 25, 26, 33, 37, 44, 46, 51, 53, 56-58, 61, 62] . These efficiencies are numerically computed under variation of the two intrinsic parameters of fBm.

In the former case of single target search we define the setting so that boundaries play no major role. In the latter case of long trajectories hitting many targets, however, we inevitably run into the problem of boundary conditions. They are more than a technicality, as in reality all media in which search may take place (such as a lake, a forest or living cells) are surrounded by boundaries that the searching element can (or will) not exceed. In our simulations we model bounded regions by a box with reflecting boundaries. For foraging we found this setup more realistic than using periodic boundaries, apart from the problem that the long time memory in superdiffusive fBm may lead to questionable effects in combination with periodic boundaries. Interestingly, recent experimental and theoretical studies have demonstrated that so-called active particles, which model self-propelled biological motion [86, 87] , spend most of their time close to reflecting boundaries exhibiting 'stickiness' to the walls [53, 54, [87] [88] [89] [90] [91] . A similar phenomenon has recently been reported for superdiffusive fBm [92, 93] that by definition displays strong persistence, similar to self-driven active Brownian particles. Therefore we explore the interplay between non-Markovian persistence in fBm and reflecting boundary conditions in comparison to using periodic ones. For reflecting boundaries we find very intricate phenomena determining the success or failure of search.

The main lesson to be learned from our research is that for fBm we do not observe any universal search optimisation in terms of maximising a search efficiency. On the contrary, exploring a range of different search scenarios by both varying parameters and search settings we identify different mechanisms determining search success. Boundary effects play a crucial role for search, which to our knowledge has not been sufficiently appreciated in previous work. Our studies demonstrate that search is a very flexible, complex process that sensitively depends on the interplay between its different ingredients. We believe that these findings further open up the field of biologically inspired search research [5] . In particular, our work suggests to shift the focus from finding simple universalities to developing a much broader picture of search. Within this general framework the LFFH takes its place as a special case.

Our article is organised as follows: In Sec. II we briefly review the concept of fBm and define our basic search setting by introducing all relevant model parameters. In Sec. III we study FA and FP problems, viz. saltatory and cruise foraging, of finding only the first target in a field of resources under variation of both intrinsic parameters. Numerical findings are explained by a simple analytical argument. Section IV reports results for multi-target search along a trajectory, both in saltatory (arrival at targets) and cruise (passage through targets) mode. This is done for finding either replenishing or non-replenishing targets, which in this setting becomes a non-trivial variation. Within this context, the impact of boundary conditions turns out to be crucial, yielding a wealth of different search mechanisms. In our concluding Sec. V, we give a coherent overview of all different search scenarios that we have investigated in our work and the quantities they depend on by summarising our main results.

In this section, we first review the stochastic process of fBm by explaining its main characteristics. We then define our specific search problem by identifying all relevant model parameters.

FBm in d dimensions can be generated by the stochastic equation of motion [21, 43, 64, 65, 67] x H (t) = B H (t) (1) with position x H ∈ R d at time t, where H ∈ (0, 1) is the Hurst exponent. Here, x H (t) holds for a Gaussian stochastic process with zero mean,

and position autocorrelation function

where K H yields a generalized diffusion coefficient [68] .

Here we represent d-dimensional fBm by d independent one-dimensional fBms, i.e., one for each component. Similar to recent work on two-dimensional Lévy walks [94] , other definitions of two-dimensional fBm may be possible. The power law decay of the correlation function makes the process non-Markovian in terms of a very slow decay of memory with time. As we will argue below, this memory can be physically understood in terms of nontrivial position correlation decay in time within a heat bath in which a particle moving according to fBm is immersed. Equation (3) results via the Taylor-Green-Kubo formula in the MSD

with α = 2H as the exponent of anomalous diffusion. Depending on the value of the Hurst exponent, fBm thus leads to subdiffusion, H < 1/2, or to superdiffusion, H > 1/2. This corresponds to anti-persistent, respectively persistent motion of fBm particles, as can be seen from calculating the velocity autocorrelation function of the process [67] . That is, for H < 1/2 it decays to zero from negative values for long times reflecting anticorrelations. Topologically this shows up as trajectories that display a lot of turns, see the example for H = 0.25 in Fig. 1 . In contrast, for H > 1/2 this function decays to zero from positive values yielding positive correlations. Accordingly particles move more in one direction displaying trajectories that are more elongated, see the example for H = 0.75 in Fig. 1 . These two types of correlation function decay have been experimentally observed for bumblebee flights [32] . In the limiting case of H = 1 fBm generates ballistic motion with α = 2. For H = 1/2 the correlation function decay Eq. (3) boils down to a delta function, and one recovers the Markovian Wiener process with normal diffusion, α = 1; see again Fig. 1 for a third example. Hence there are three generic cases of fBm dynamics, anti-persistent subdiffusion, normal diffusion and persistent superdiffusion depending on the value of H as displayed in Fig. 1 . This discussion suggests that fBm is a generalisation of Brownian motion by correlating the position autocorrelation function. To see this more clearly one rewrites Eq. (1) in the form of a stochastic differential equation [18, 68] ,ẋ

where f H (t) denotes d-dimensional fractional, i.e., power law correlated Gaussian noise (fGn) [21, 64, 65, 95] ; see Refs. [43, [67] [68] [69] for details. Equation (5) can be understood as an overdamped Langevin equation [43, 67] . Within this framework the left hand side can be interpreted as a constant friction term without memory and the right hand side as a correlated random force. The latter models collisions of a tracer with heat bath particles that perform dynamics with power law memory decay. We remark that related equations have recently been used as models of active Brownian particles [96, 97] . These particles are self-propelled due to the fact that by definition, as in Eq. (5), the fluctuation-dissipation relation is broken [98] . That is, here we only have memory in the noise but not in the friction [67] . This type of dynamics has also been applied to model biological cell migration [79] . FBm can be generated numerically by different methods either (theoretically) exactly [99] [100] [101] or approximately [102] [103] [104] . In this work we obtain fBm from discrete time fGn using Hosking's method [99] . According to Eq. (5), the increments of fBm are then computed by the integral B H (t) = t 0 f H (t ′ )dt ′ . In practice, we discretize the motion of an fBm particle to a series of successive jumps at unit time steps, t 0 = 1. Then the integral breaks down to the discrete sum

This means that such a timediscrete fBm particle performs ballistic jumps during unit times t 0 according to the increments generated from fGn at time t. The direction of the subsequent ballistic step is then determined by the next increment of the fGn. The underlying Hurst parameter controls whether these consecutive steps are positively of negatively correlated. Importantly, this time-discrete set-up allows us to study both FP and FA search. In the latter case, a searcher may cross a target without finding it. This is in contrast to rigorous time-continuous fBm, where due to the strict self-similarity of the trajectories FP is identical to FA.

Evidently, this self-similarity is absent for time-discrete trajectories generated at short times t ∼ t 0 ; see [78] for a related discussion. We note, however, that in a suitable scaling limit the time-discrete fBm as defined above converges to the mathematically exact fBm.

In analogy to Ref. [13] , we consider the basic setting where a searcher moves in a plane, d = 2, to find uniformly random distributed targets, see again Fig. 1 . The motion of the searcher is obtained from the aforementioned fBm process. This means that according to Eq. (3), the searcher has a memory of its past positions. The memory could arise from extrinsic environmental conditions, for instance when a biological cell diffuses in a visco-elastic medium where the diffusing element is part of a larger evolving system [70] . Or, it could be an intrinsic property of a searcher, like internal memory during a foraging process [7, 8, 79] . Within our fBm framework we consider the memory to be intrinsic. In terms of movement, one may think of memory as generating persistence or anti-persistence between the different steps in the process. We say the search process exhibits non-renewal if the memory lasts during our whole measurement time.

If the memory has a certain duration, i.e., the memory kernel in an fBM process possesses a cut-off, we say the process has a renewal. Accordingly, below we may speak of resetting the process to some initial condition if we truncate it, or of non-resetting. The cut-off time can be a constant value or the process can be renewed after visiting a target. The latter situation was considered in Ref. [13] by truncating a LW upon hitting a target. This aspect will become important when we define two different types of search efficiencies further in this article.

Within this setting we consider the two different search modes of cruise and saltatory foraging [2, 9, 33, 44-46, 56, 57] briefly mentioned in Sec. I A. In detail they are defined as follows:

1. A cruise forager perceives a target while moving.

Accordingly, whenever a target is passed it is found. The special problem that a cruise forager detects only the first target yields mathematically a FP problem [2, 25, 26, 37, 45, 46, 50, 62, 85 ].

A saltatory forager cannot find a target while moving. It has to land on the target after performing a jump, or sufficiently close to it, to find it. After a jump the searcher typically changes direction. In the case of first target search this can be formulated as a FA (or hitting) problem [33, 38, 44, 45, [56] [57] [58] .

However, our setting is not yet complete, as we have to specify how a searcher 'perceives' a target. This can be modelled in two ways [2, 13, 37, 46] :

A. We may assign a radius of perception R p to a searcher. This means that within R p it perceives a target with certainty. In this case even point targets are found. Here, the perception radius is a parameter intrinsic to the searcher.

B. The searcher is blind and only finds a target when hitting it. In this case the searcher is reduced to a point without any perception. It thus can only find targets that have an extension.

Combining case 1 with cases A and B yields the following two rules for cruise search [27, 44, 46] :

with a tube of radius R p around the searcher's trajectory it is found.

If the searcher's trajectory crosses a target it is found.

For saltatory search we get accordingly [44, 46] : 2.A. A perceptive saltatory forager: If after a jump of duration t 0 the target intersects with a circle of perception radius R p around the searcher's position it is found.

If after a jump of duration t 0 the searcher's position is within a target it is found. Figure 1 depicts the special case of 1.B. where a blind cruise forager searches for disks of radius R e . This suggests that R e defines an extrinsic parameter specifying the environment. But for circular targets one can scale away either R p or R e by combining both parameters. As we have already encountered in the case of memory, speaking here of an extrinsic parameter is consequently ambiguous. These interpretations thus depend on the specific situation at hand. Consequently, for disks case 1.A. is mathematically equivalent to 1.B. and 2.A. to 2.B. In the search scenario shown in Fig. 1 we therefore only need to distinguish between cruise and saltatory foraging, cases 1. and 2., which is what we study in the following.

Our second extrinsic parameter quantifies the density of targets. We consider the case where N point targets are uniformly random distributed in a two-dimensional box of side length L. This results in a target number density of n = N/L 2 . The mean (free flight) distance l f between a point target and its nearest neighbor is thus calculated to

see Fig. 1 for the pictorial meaning. If not said otherwise, we choose reflecting boundaries for the simulation box. That is, when a particle crosses a boundary during a displacement drawn from fGn, it is specularly reflected back to the bulk. But this does not affect its next step governed by the memory in the fGn. This means we do not 'truncate' the trajectory by hitting a wall but keep its memory. An illustration of this process can be found in the inset of Fig. 6 (a). An important intrinsic length scale determining such a search process is the mean jump length of a searcher. Note that this only comes into play because of our choice of time-discrete fBm thus yielding a slight variation of mathematically rigorous fBm. If we assume that the trajectory of a searcher consists of a series of consecutive jumps with time steps t 0 , cf. Fig. 1 , using the MSD suggests an average displacement of a searcher during t 0 as

which is the standard deviation of fGn, see Eq. (5). By Eq. (4) one can express l d in terms of the generalized diffusion coefficient K H and the time step t 0 as

As mentioned before, in our simulations we set t 0 = 1. This is convenient, since that way in Eq. (8) we decouple l d from α = 2H. We thus vary l d by varying K H , which in turn depends on the fGn strength f H (t). While l d relates to the strength of diffusion, the second intrinsic parameter α, which we introduced before, determines the type of diffusion. For our studies, we keep the two extrinsic parameters R e and l f fixed, which defines a specific search environment. We then explore the impact of varying the two intrinsic parameters α and l d on two generic types of foraging efficiencies, which we define later in the text. But first we select the basic environmental regime viz. the properties of the targets that we focus on in this work in relation to all three length scales R e , l f and l d introduced above. The situation of scarce targets is perhaps the most interesting one as it poses the challenge to efficiently locate a target after many time steps. This regime of target densities translates geometrically to the condition that the effective target radius is much smaller than the mean distance between two point targets, R e ≪ l f . This means that after finding a target, there exists no other target for a searcher within radius R e . Furthermore, for effectively modelling low target densities, the mean jump length should be much smaller than the mean target distance, l d ≪ l f . We remark that it is typically assumed for LWs that l d ≪ R e [13, 26, 37, 53] . For all our simulations, the box size, effective radius and mean target distance were set to L = 10000, R e = 1 and l f = 40. Instead we varied the jump length 0.04 ≤ l d ≤ 24 and 0 < α < 2.

In this section, we study the problem of finding only the first target. For this purpose, we consider an ensemble of searchers and numerically compute two suitably defined search efficiencies. We do so under variation of the exponent α = 2H of anomalous diffusion for different mean jump lengths l d , in the case of both saltatory and cruise foragers. We explain our numerical results heuristically and by a simple analytical approximation.

For finding a single target, a search process consists of a starting point represented by a respective initial condition of a searcher and an end point, which is when a target has been found. For a saltatory searcher, such a search process is characterised by the mean FA time [33, 38, 44, 45, [56] [57] [58] . For a cruise searcher the corresponding quantity is the FP time [2, 8, 25, 26, 37, 45, 46, 50, 62, 85] . The situation of single target search applies, for instance, to non-recurrent chemical processes, where as soon as a reactant finds a target, the chemical reaction takes place and the search ends [2, 10, 42, 43] . To simulate this search scenario, we let a searcher start from a randomly chosen initial position. Note that the choice of initial condition is non-trivial [20, 33, 44, 45, [56] [57] [58] 62] , as we will discuss in Sec. IV for replenishing targets. When the searcher finds the first target, the process starts again from another random initial position. This can be considered as a resetting procedure, however, as eventually we average over all initial conditions, in effect this yields an ensemble average of searchers with respect to random initial conditions. In order to avoid here, as far as possible, the impact of boundaries on the results, the initial position is chosen from a small box in the center of the main simulation box. The size of the box is such that for a smaller jump length l d the probability that the searcher will find a target before reaching the boundaries is very large. For larger l d , however, boundary effects do come into play, as we discuss below. They will be investigated in full detail in Sec. IV.

Based on the mean FA (FP) time written as t A/P , where the angular brackets denote an ensemble average over both random initial conditions of the searcher and different target positions, we define the corresponding efficiency η A/P as [8, 10, 13, 25, 37, 46, 62] 

For the ensemble average we considered 10 5 simulation runs. After 500 runs, we regenerated the (uniform) distribution of point targets thereby averaging over this distribution. Since we can follow trajectories numerically for only finite times, we set the maximum time for each search process to T = 10 5 t 0 . If until then no target was found, we stop the search and use T for the corresponding trajectory to calculate the average in Eq. (9) . In that sense, modulo small statistical errors we obtain an upper bound for the efficiency η A/P . The results of these simulations for η A/P under variation of α for different jump lengths l d are displayed in Fig. 2 . Panel (a) depicts the efficiencies calculated using FA times, (b) is for FP times. Each curve in the figure is normalised with respect to the maximum value of the efficiency obtained at the corresponding jump length, i.e.,

The reason for the normalisation is that the unnormalised efficiencies vary over several orders of magnitude with l d , see Fig. 3 (b) that we discuss afterwards. Figure 2 shows that for very small jump lengths l d compared to the effective radius R e = 1 and the mean distance l f = 40, which are both held constant, the exponent α of anomalous diffusion that optimises both efficienciesη A andη P lies in the superdiffusive regime close to the ballistic limit α = 2. The physical explanation is that a searcher needs to compensate for small values of l d by performing quasiballistic motion in order to move at all through in space. If l d increases for both efficiencies, persistent motion close to α = 2 is not optimal anymore for finding targets. This reflects undersampling, i.e., instead of efficiently exploring a given area, a persistent searcher with large α, which performs in addition jumps of large length l d , moves immediately to another area by starting there its search again from scratch. But since it encounters a new random target distribution in this area, the searcher loses time. This has two important consequences for bothη A andη P . First, for larger l d both efficiencies develop a minimum at larger α values. Second, accordingly a maximum emerges in bothη A andη P for smaller α. In this region there appears to be an optimal interplay between a larger jump length l d and less persistence in the motion of a searcher leading to more frequent turns, which compensates for too large jumps by yielding an optimal scanning of a given area. While this is observed for both FA and FP, a crucial difference between both foraging modes is that for large l d , FA exhibits leapovers [33, 44, 45, 56] , where a searcher misses a target by jumping across the target without landing on it. We will discuss this phenomenon in more detail in relation to Fig. 3 (b) . Another observation is that, as l d increases by approaching R e ≃ 1, the optimal exponent α forη P seems to accumulate around α = 1.5, while forη A it exhibits a shift towards the normal diffusive regime around α = 1. This indicates that for FP, there is a special regime of parameters l d and α that optimises search related to the scanning procedure discussed above. This regime avoids both undersampling for large α, as well as oversampling [2, 20, 33, 37, 46, 53, 56, 59 ] of a small region by too many nearby trajectories for small α. In contrast, forη A there is a less pronounced accumulation of the maximum for larger l d with the peaks shifting and flattening out. That is,η A becomes to some extent independent of α around α = 1 for larger l d . This suggests that oversampling for small α is less of a problem for FA processes. which intuitively looks plausible, as this search process is characterised by random points in space and not by the full trajectories. We call this flattening out ofη A in α the 'paradise' regime, as for l d ≥ l f /2 the jumps of the searcher become large enough that it can always find a target irrespective of its persistence if α is not too large or too small. The factor of 1/2 is due to choosing random initial positions between two targets. Note that, to some extent, a paradise regime also exists for FP search at the largest l d = 24, however, it is not as pronounced as for FA processes.

For large l d and large α > 1.8, boundary effects come into play. However, here we argue that they do not play a role for the results presented in Fig. 2 and their interpretation, as follows: In Fig. 3 (a) we plot the fraction of search processes that failed for finding a target during the given maximum time T . Accordingly, this quantity yields the survival probability P , i.e., the probability of finding no target until time T , which in a way is an inverse measure of search efficiency. Shown are FP results for P as a function of α for different l d . But we remark that P is essentially the same for FA, as it is significantly different from zero only for small α or small l d , where there are no leapovers. The most important result of Fig. 3 (a) is the existence of two maxima at small and large α, which correspond to the respective minima of η A/P in Fig. 2 . The family of curves displaying small maxima around α = 2 for large l d values is indeed due to boundary effects, which will be explicitly investigated in Sec. IV. In this case, if a searcher hits the reflecting wall with a large velocity component perpendicular to it, due to the strong persistence in the fBm motion for large α, it will bounce back and forth off the wall for a long time. Consequently it is essentially stuck in one area by the wall, which means that no target may be found anymore. This effect is worse for smaller l d values, while for larger l d the searcher can still explore larger areas, which explains why the small peak around α = 2 decreases for larger l d values. Note that in any case, the fraction of searchers affected by this boundary effect is very small (typically much below 5% for l d > 0.4, cf. Fig. 3 (a) ). Therefore, this does not explain the drop-off in the efficiencies in Fig. 2 at larger α values, which conversely becomes stronger for larger l d values. For small α values, in most cases of l d the subdiffusive search times are so long that they go beyond the numerically accessible regime. But as explained above, for α < 1 all efficiencies safely yield at least an upper bound for the exact efficiency values. Figure 3 (b) supplements the analysis of normalised efficienciesη A/P presented in Fig. 2 by showing the unnormalised counterparts η A/P , here as functions of l d for three different α. As we mentioned before, both quantities vary over orders of magnitude. Most notably, η P increases monotonically close to a power law in l d (for exponents see Fig. 5 (a)) while η A saturates for larger l d . The saturation is a consequence of the leapover phenomenon discussed above, which exists for FA but not for FP. Indeed, the saturation of η A sets in around l d ≃ 2 = 2R e , which matches exactly the condition where a searcher can jump over the full diameter of a target without finding it. We also note that η A seems to decrease slightly for larger l d , which might be due to further undersampling related to the leapovers. Furthermore, both efficiencies decrease with larger α. This is in line with Fig. 2 , where it was explained by undersampling.

We now relate these results to previous works on LWs in which similar search scenarios have been studied [10, 13, 37] . As explained before, LWs and fBm define in principle two fundamentally different stochastic processes. However, in both cases there are important parameters that govern the type of diffusive spreading. For fBm, this is the exponent 0 < α < 2 in the correlation function decay Eq. (3), which determines the MSD Eq. (4). For LWs, in turn, the crucial quantity defining this process is the distribution P (r) from which the jump lengths r are sampled randomly at each time step, which is assumed to follow a power law, P (r) ∝ |r| −µ with 1 ≤ µ ≤ 3. Via continuous time random walk theory, the MSD can be calculated for this process, which depends on µ in a more complicated way [24] . Crucially, for µ = 3, a LW reproduces normal diffusion, while µ = 1 yields the limit of ballistic motion. As far as diffusion is concerned, one may compare the LW parameter regime of 1 ≤ µ ≤ 3 with the respective parameter regime of 0 ≤ α ≤ 2 for fBm. A basic difference is that LWs do not generate subdiffusion, hence there is no matching for 0 ≤ α < 1 in fBm to a corresponding parameter regime for LWs. However, one may identify qualitatively α = 2 with µ = 1 and α = 1 with µ = 3, since in the case of the former parameter values both processes yield purely ballistic motion, while in the latter case they reproduce normal (Brownian) diffusion.

It has now been claimed in the literature that for LWs, a universal exponent around µ = 2 yields an op-timal efficiency for finding sparsely distributed targets [10, 13, 25, 26, 28, 29, 61, 62] . Intuitively, a LW for µ = 2 suitably combines properties of the two extreme cases of ballistic and Brownian motion generating trajectories that explore a given area by avoiding oversampling (too frequent turns), as well as undersampling (too straight trajectories). However, this points to our previous argument about the accumulation of optimal efficiencies around α = 1.5 for fBm with FA for which we have given an analogous microscopic explanation. One may thus conclude that, under certain conditions, a search process that is intermediate between two extreme cases may indeed optimise search efficiencies. We emphasize, however, that we see no universality of such exponents for optimising search, as their values depend very much on the precise conditions of the search problem at hand. For example, our FA efficiency depicted in Fig. 2 does not exhibit any clear accumulation of maximal efficiencies around a particular α, and even for FP this effect is quite washed out. This is fully in line with the conclusions in Ref. [1] that the optimal exponent µ for LWs largely depends on scaling, as has been substantiated in the very recent Refs. [37, 40] . For a LW search in more complex settings, it has been found that while for similar qualitative reasons an optimal exponent typically appears to exist, its value depends on the variation of intrinsic and extrinsic parameters, as well as the precise topography at hand [51, 53] .

For calculating search efficiencies, Eq. (9) does not always yield a sensible definition, as in certain situations (like single target search on the line by Brownian motion) the mean FA and/or FP times may diverge [42, 56] . Hence, in this case one would only obtain the trivial result η A/P = 0. It was thus proposed to redefine efficiencies by using the different average [33, [56] [57] [58] 

Note that this definition implies a completely different weighting of search times compared to Eq. (9): While long times now only mildly suppress the value ofη A/P yielding small inverse values, the large inverse values for short search times contribute to the efficiency more profoundly. This relates to the specific properties of the tails of FA and FP time distributions which, however, are numerically difficult to obtain because of their extreme statistics. Figure 4 shows simulation results for the renormalised efficiencyη A/P of both FA and FP searches, where we set the time for failed processes to infinity. The data is, as before, for 10 5 runs with a length of T = 10 5 t 0 time steps. By comparing these curves with the ones in Fig. 2 we see, first, that for the efficiency obtained from Eq. (11), regardless of the jump length, the most efficient way of finding a target is to perform ballistic motion.

We also observe that by decreasing α the efficiency is decreasing, which for α < 1 reproduces roughly the same trend as in Fig. 2 . However, the new efficiency wipes out any non-trivial dependencies due to long search times. This means that the decay for large α is completely gone and correspondingly there is no local maximum anymore. Note that numerically there is a sharp drop ofη A/P for α → 0, which is very difficult to resolve computationally due to difficulties with properly capturing very long search processes, see Fig. 3 (a) . For practical purposes, one may thus need to make a sensible choice as to which definition to use for calculating efficiencies depending on the situation at hand. We also remark that this example shows very clearly how profoundly optimality depends on the definition of efficiency that one employs. 

In view of the non-Markovian nature of fBm and the complexity of two-dimensional multi-target search, providing a theory for respective FA and FP problems is a non-trivial task [10, 13, 26, 37-40, 85, 105] . Here, we outline an extremely simple, handwaving argument that analytically reproduces some of the parameter dependencies for the FP efficiencies as seen before. We emphasize that this argument has to be handled with much care, as we discuss in the following subsection, where we relate it to more substantial, general theories published in previous literature [85, 105] .

We start by simplifying our two-dimensional FP problem to a one-dimensional setting as follows: We interpret the mean free distance l f between two targets as the length of a line between two absorbing boundaries. Let the searcher start the diffusion process right in the middle of the line. Then it needs to travel a distance l f /2 towards the left or the right to hit a target. For estimating the mean FP time, one may now simply use the MSD of the fBm process given by Eq. (4) . Replacing x H (t) = l f /2 and solving for time t we obtain

The generalized diffusion constant K α can be in turn approximated to, see Eq. (8) with d = 1,

where in our case t 0 = 1. We can now calculate the FP efficiency η P in terms of all relevant parameters by substituting the FP time in Eq. (9) by the one of Eq. (12) supplemented by Eq. (13). This yields

Our simple argument may be understood as a kind of mean field approximation, in the sense that we use a single-target picture to approximate multi-target search in a low density limit. It should also apply to FA search as long as leapovers do not dominate, see Fig. 3 

We test the validity of our approach by comparing the dependencies of η P on l d and α as predicted by Eq. (14) with numerical data. Figure 5 (a) shows the unnormalized efficiency η P as a function of the mean jump length l d for different exponents α of anomalous diffusion (cp. with Fig. 3 (b) ). While due to the simplicity of the theoretical argument we may not expect a full quantitative matching between data and approximation, at least the power law dependence of η P with l d for different α is reproduced surprisingly well. Note that for large α we have restricted our fit region to smaller l d , as we may not expect our theory to capture the effect of undersampling. For smaller α, we need to go to larger l d due to the otherwise high percentage of failed searches, see Fig. 3 (a) . We did not include results for α < 0.91, as here the efficiencies become so small that they are difficult to compute numerically, cf. again Fig. 3 (a) and our respective diuscussion. Figure 5 (b) displays results for the normalised efficiencyη P as a function of α for different l d (selected range of data from Fig. 2 (b) ). The fits have been adjusted to an intermediate region of α values, as for small α the data is not reliable due to many failed searches, see Fig. 3 (a), which yields insufficient statistics to check for fine details. For larger α we may not expect Eq. (14) to reproduce the non-monotonic dependencies ofη P due to undersampling, see in Fig. 2 (b) . Again, in this intermediate regime our handwaving theory predicting an exponential dependence on 2/α works surprisingly well. 

Here, we first show that our simple approximation Eq. (14) forms a special case of a much more advanced theory for calculating mean FP times of fBm [85] . This theory, on the other hand, is recovered within the framework of a more general theory [105] which, if applied to our particular setting, indicates important limits of validity of our approach.

The first layer of embedding is provided by the theory of Guérin et al. [85] , who derived analytical results for the mean FP time t P of compact (as explained below) non-Markovian random walkers for finding a single target in a d-dimensional volume V confined by reflecting boundaries. By applying a generalised form of a renewal equation, for fBm in one dimension the mean FP time was calculated to

where x 0 is the initial condition of the searcher with a target at position x = 0. The constant β H is a nontrivial quantity that captures the non-Markovianity of the process and can only be computed numerically [85] . Comparing our trivial formula Eq. (12) with Eq. (15), one can see that the former yields the same dependence of the mean FP time on K H as the latter. Furthermore, replacing x 0 = l f /2 in Eq. (15) according to our simplified assumption, as well as using V = l f for our onedimensional setting, for the FP time in Eq. (12) we obtain exactly the same scaling with l f as in Eq. (15) . To what extent the scaling of the FP time with α for fixed K α and l f is reproduced is not so clear, as in detail this depends on β H . The relation between Eqs. (12) and (15) also explains why our analytical approximation cannot be applied to the efficiency definition Eq. (11) . We thus conclude that our approximation Eq. (12) corresponds to a simplified version of Eq. (15) if strictly constrained to one dimension. A more general theory for calculating mean FP times of non-Markovian scale-invariant, aging diffusion processes, which includes the one of Ref. [85] as a special case, was developed by Levernier et al. [105] . A crucial feature of this theory is to distinguish between compact and non-compact stochastic processes, as these two cases lead to different formulas for the corresponding mean FP times. Compactness (or recurrence) means that a process comes back to its starting point with probability one while non-compactness (or transience) implies the opposite [42, 43, 106, 107] . Consequently, this mathematical property intimately relates to what before we called oversampling, respectively undersampling. For non-ageing dynamics [105] , a simple criterion for compactness is derived from the walk dimension d w of a stochastic process, defined by x 2 (t) ∼ t 2/dw (t → ∞) [43, 105, 107] . If d w is greater than the dimension d of the embedding space of the walk, d w > d, a process is called compact. Conversely, non-compactness means d w < d; the case d w = d is called marginal [43, 106, 107] . For fBm with d w = 1/H, cf. Eq.(4), we thus have that in dimension d = 1 it is compact for all 0 < H < 1. Hence, strictly speaking Eqs. (12) and (15) only hold for compact one-dimensional fBm.

In higher dimensions, the theoretical results for mean FP times of both Ref. [85, 105] are in many respects much more complicated. For the particular case of fBm in d = 2 dimensions, the above definition yields that fBm is compact for H < 1/2 while for H > 1/2 it is noncompact. This is reflected in different formulas for t P in both regimes [105] . By using the notation of Ref. [105] (see the results in Tables I and II therein) , let R denote the confining domain size, r the distance the walker starts from a target and a the target size. Working in the limits R ≫ 1 and a ≪ r, it was shown that for compact fBm in d dimensions

while in the non-compact case

Let us now assume that R = l f and r = l f /2, as in our approximation, and in our setting we have a = 2R e . For one-dimensional (compact) fBm the scaling of t P ∼ l 1/H f is then recovered from Eq. (16) , in agreement with the one-dimensional fBm results Eqs. (12) and (15) . The very same scaling holds for d = 2 but only in the compact case of H < 1/2 while for H > 1/2 Eq. (17) leads to t P ∼ (l f /R e ) d R 1/H e . For two-dimensional fBm the theory of Ref. [105] thus predicts different laws for t P in the compact and non-compact regimes. In contrast, our numerical results Fig. 5 (b) display a smooth transition around H = 1/2 viz. α = 1 under variation of α, well fitted by the single formula Eq. (12) that according to the above theory should only hold for H < 1/2 in one dimension.

These seemingly different facts can be reconciled as follows: References [85, 105] consider the case of a walker starting at a given initial condition that searches for a single target at another given position. In contrast, we reported results for finding one out of many randomly distributed targets by averaging over random initial conditions of the searcher. In our approximation, this is roughly captured by assuming R = l f and r = l f /2 (see above). But these two particular assumptions cancel out the dimensionality dependence d in the compact case Eq. (16) . Hence, the respective one-dimensional result carries over to two dimensions, which in turn is what we assumed for our approximation. We also note that both Eqs. (16) and (17) are proportionalities (up to a prefactor independent of geometric parameters) yielding no information about the generalised diffusion coefficient K H and the associated jump length l d that we varied. Finally, for our choice of fixed parameters we have a = 2R e = 2 while r = l f /2 = 20, which does not quite match to the condition a ≪ r underlying the derivation of both Eqs. (16) , (17) , nor do we necessarily work at low densities of targets. In that respect the theory developed in Ref. [105] describes a different asymptotic parameter regime compared to our setting. We thus argue that it cannot be applied directly to understand the transition in the efficiency around α = 1 displayed in Fig. 5 (b) . To explain in further detail why for our data the transition in α looks smooth and is reproduced by the formula that according to Eqs. (16) , (17) should only hold for the compact regime thus remains an interesting open question.

In view of the above theoretical results it becomes clear, however, that in more general situations our simple approximation Eq. (12) can only be of rather limited validity. Typically, the dimensionality dependence of t P will not cancel out, and one may thus not trivially extrapolate a result for one to higher dimensions; see Refs. [37, 39, 40] for an analogous discussion in the case of LWs. Exactly for that reason we expect Eq. (12) to be wrong under variation of l f . To study the dependence of t P on the density of targets hence yields another interesting problem for further study. While thus there remain challenging open questions, we see no contradiction between neither our numerical results nor our handwaving approximation and the theories developed in Refs. [85, 105] .

We now study the problem of a searcher that, without resetting, finds many targets during one run. First we introduce the basic framework of this search problem by defining a suitably adapted search efficiency. This type of search very much depends on details of the environment, in particular properties of the target and boundary conditions. We first explore the search for targets in the bulk, here both replenishing and non-replenishing resources, before we investigate the impact of boundary conditions on finding replenishing targets. The latter case establishes a cross-link to the very recent field of active particles.

In numerous realistic situations, such as recurrent chemical reactions or animals looking for food [2, 8, 10, 43] many targets need to be found by a single searcher. Compared to the search problem of Sec. III, one may call this scenario non-resetting, as here a searcher consumes targets along a single path generated by its continuous movements in time. In Sec. III, efficiencies were calculated by the inverse of FA and FP times, which in turn were defined as ensemble averages over many searchers starting at different initial conditions, see Eq. (9) . In order to adequately describe multi-target search along a trajectory, one replaces this average by a combined timeensemble average as follows: First one counts the number N of targets that have been visited along a trajectory during a given time T [8, 10, 13, 37, 46, 51, 61, 62] . The fraction T /N then yields the time average for finding N targets along the trajectory of a searcher. This is easily seen by < t >= 1/N N −1 k=0 (T k+1 − T k ), where T 0 = 0 is the initial time and T k the time to find the k-th target. For large N we can neglect (T − T N )/(N + 1) and have < t >≃ T /N . If not said otherwise, for the following results we fix T = 10 6 . Since for this finite search time the outcome may still depend on the specific path of the searcher, in addition we average over an ensemble of searchers starting at random initial positions. And we typically choose 10 4 simulation runs. In analogy to Eq. (9) , an adequate definition of efficiency is then obtained by [8, 37, 51, 61, 62] 

As before, the angular brackets denote the average over an ensemble of searchers starting at random initial positions, here with respect to the number N A/P of targets found during the total time T while arriving (A) at or passing (P) through targets. It is trivially clear that the above equation gives nothing else than the inverse of the combined time-ensemble average for the average time T A/P = T / N A/P along a path. In analogy to Sec. III, in what follows we numerically investigate this search situation for both saltatory and cruise searchers, viz. arrival at and passage through targets, by varying the exponent α of anomalous diffusion. We do so for the same fixed extrinsic parameters L, R e , l f as before, cf. Sec. II, by choosing different jump lengths l d . The key quantity to compute is again the efficiency, in this case defined by Eq. (18) . However, in contrast to first target search, for consecutively finding many targets one needs to distinguish between two different types of resources [8, 10, 13, 37, 51, 60 ].

1. Non-destructive, or replenishing targets: After a target is visited by a searcher, it remains intact and can be revisited. Typically, such a target is modelled as a replenishing resource, i.e., it reappears either when a certain delay time has passed after its consumption [26] , or after the searcher has passed a certain cut-off distance away from it [27, 37] . This prevents the searcher from artifically consuming the same target over and over again. For our subsequent studies, we choose a cut-off distance of d c = 2R e = 2.

Destructive, or non-replenishing targets [53, 61] : After being visited by a searcher the target disappears forever.

These two extrinsic environmental conditions for the targets define completely different search scenarios [10, 13, 37, 60] . An additional crucial complication that we will explore in detail is the interplay between these two different target types and the boundary conditions.

We start by investigating the situation of replenishing target search in the bulk, i.e., without elaborating on the impact of boundaries. Figure 6 shows simulation results for the efficiencies η * A/P , Eq. (18), of both arrival (a) at and passage (b) through replenishing targets under variation of the exponent α of anomalous diffusion for different jump lengths l d . We see that for the smallest l d = 0.04, in both cases ballistic motion outperforms any other type of motion, in analogy to Fig. 2 . The physical explanation is the same as for Fig. 2 : A large α needs to compensate for a small l d for the searcher to move anywhere. Furthermore, if a target is found, for small jump lengths the cut-off distance d c translates into long delay times before the visited target reappears. Hence, the return time to revisit the same target is very long, which explains why efficiencies are close to zero in the subdiffusive regime of α < 1, where returns dominate the search due to antipersistence in fBm. This particular type of return dynamics, which we call the revisiting target mechanism, will subsequently become very important. Note that not any subdiffusive dynamics is characterised by returns. A counterexample is a continuous time random walk, where subdiffusion originates from power law waiting times at a given position [15] . This dynamic should thus yield very different search efficiencies compared to fBm.

For the slightly larger jump length l d = 0.4, the region of maximal values of the two efficiency curves widens a bit by including slightly smaller α values. One may speculate that, as in Fig. 2 , the maximal efficiency now starts to shift to smaller α values due to undersampling around α = 2. Overall, both efficiencies are getting larger for all α values. That this happens in the region of α < 1 can again be explained by the revisiting target mechanism described above: Notably, for larger l d it now switches from slowing down search to enhancing it, because the searcher leaves the target more quickly after finding it, leading it to replenish quickly. However, due to the anti-persistence of fBm, for α < 1 the searcher more frequently returns to the same target. Hence, by increasing l d subdiffusive dynamics yields a new important search strategy for exploiting replenishing targets. That subdiffusion can enhance search success has also been reported for a very different search setting in Ref. [82] . Very recently, by analysing experimental data it has been found that subdiffusion adequately describes the area-restricted search of avian predators [41] .

For the next largest jump length l d = 1, instead of the maximum around α = 2 further shifting to smaller α values as in Fig. 2 , there is a dramatic increase of both efficiencies towards α = 2 again, which is in sharp con- It shows a corner of the wall confining the simulation box. The black line with the arrows represents the trajectory of an fBm particle that has hit the wall at least once already. The blue dotted lines depict how it would continue without the presence of the wall. The red dotted line indicates how the particle is specularly reflected at the wall. For a particle with persistence, this leads to the phenomenon of 'stickiness' at the wall. The grey disks are scatterers, as in Fig. 1. trast to the first target search. We will argue in the next Sec. IV C that this is due to an interplay between persistence in fBm and the reflecting boundaries that we have chosen. But as this is a highly non-trivial problem by itself, we now primarily focus on α < 1.5 where searchers typically do not hit the boundaries and bulk dynamics dominates the search. Here, both efficiency curves fur-ther increase by increasing l d compared to the two previous smaller jump lengths by flattening out up to α < 1.5. This increase may again be explained by smaller delay times for a visited target to reappear when l d is getting larger, which enables normal and sub-diffusive search to more strongly contribute to search success due to the revisiting target mechanism. Note that at the smallest α values the search times become very large and have to be numerically truncated, cf. Fig. 3 (a) and our respective discussion. Hence, beyond the numerically supported general trend to small efficiencies, we cannot resolve whether the shown fluctuations are due to statistical errors or hint at more subtle local non-monotonicities.

If we further increase the jump length to l d > 1, we obtain two different families of efficiencies for arrival compared to passage search for all α values. This is because of the onset of leapovers for arrival search as discussed in Sec. III, see particularly Fig. 3 (b) . That is, while η * A generally starts to slightly decrease for all α by increasing l d , in line with our results for first target search in Fig. 3 (b), η * P generally keeps slightly increasing with l d until a quite perfect plateau region has been reached for α < 1.5 up to the largest l d considered here. It is intuitively clear that for passage search, a larger l d should generally increase the search efficiency unless there are other effects mitigating this mechanism.

Remarkably, we observe two rather spectacular transitions at the largest α values, from maximal efficiencies for l d < 1 to minimal values for l d > 1, and then the reverse between l d = 2.8 and l d = 4. This happens for both arrival and passage search and thus cannot be attributed to the onset of leapovers as discussed above. In the following section, we will argue that these two transitions are again subtle boundary effects. Correspondingly, for 1 < l d < 2.8 now non-ballistic search with α < 1.5 yields the largest efficiencies for both passage and arrival search. This can be understood again by the revisting target mechanism introduced above. We argue that, surprisingly, subdiffusive search maximises our search efficiencies within this l d parameter regime, cf. again Ref. [82] for related results. Figure 7 demonstrates that this mechanism provides indeed the correct explanation. Similar to Fig. 6 , it displays results for both efficiencies η * A and η * P as functions of α, here for the two particular jump lengths l d = 1, 2. But in this case, we have non-replenishing targets that are destroyed after a visit, hence the revisiting target mechanism cannot contribute to search success anymore. Note that the y-axis for the values of the efficiencies is scaled by a factor of 10 −6 , in sharp contrast to Fig. 6 where the scale is 10 −3 . Overall, both search efficiencies are diminuished dramatically in the case of nonreplenishing targets. In the subdiffusive regime of α < 1, the efficiencies are indeed close to zero, which confirms the importance of the revisiting target mechanism. We furthermore observe, at least with respect to two values of l d , that both efficiencies decrease for all α by increasing l d > 1. While for arrival this is in line with Figs. 6 (a) and 3 (b) due to the onset of leapovers, for passage this is exactly the opposite to Fig. 6 (b) . One may speculate that this reflects some interplay between the periodic boundaries and FP search. However, based on only two curves more detailed conclusions cannot be drawn and remain open for further research. The existence of local maxima in 1 < α < 1.5 viz. the suppression of efficiencies for α → 2 reflects again boundary effects, as will be explained in the next section. 

We now explain the two transitions between maxima and minima around α = 2 in both η * A and η * P shown in Fig. 6 to which we referred in passing above. In the strongly superdiffusive, quasi-ballistic regime of α > 1.5, a searcher quickly approaches the boundaries of the system. But since we have chosen reflective walls, the searcher is thrown back into the bulk after hitting a boundary. However, since for α > 1.5 fBm displays strong persistence in the motion, the searcher will immediately move back to the wall by getting reflected again, and so on. This microscopic mechanism is illustrated in the inset of Fig. 6 (a) . It leads to the phenomenon of stickiness to the wall, see the blue trajectory in Fig. 8 , recently investigated for fBm in Refs. [92, 93] . This is more widely known for active Brownian particles that by definition of activity exhibit persistence in their motion [53, 54, [87] [88] [89] [90] [91] . In more detail, when hitting the wall one can decompose the velocity of a searcher into a component parallel and one perpendicular to it. Since the probability that a searcher hits the wall strictly perpendicular to it is zero (with respect to Lebesgue measure in angular space), there will always be a component of the velocity parallel to the wall. Typically, an fBm searcher will thus for a long time move along the wall by displaying zig-zag quasi one-dimensional quasi-ballistic creeps. These are constrained to a boundary layer of an approximate width l d , which defines a crucial boundary length scale. But on top of that, two further length scales come into play, governing the search process in this boundary layer. The second one is the effective (perception) radius R e = 1, which defines the relevant parameter for first finding a target. If a searcher is moving in a boundary layer of width l d in which there are randomly distributed (point) targets, the efficiency for finding targets will be best if approximately l d ≤ R e , where R e can alternatively be interpreted as determining the average extension of a circular target, cf. Sec. II. We call this maximisation of efficiency the pac-man effect, in analogy to an old computer game where a searcher subsequently eats targets by moving along channels in a maze, as this channeling helps to locate targets by increasing the search efficiency; see again the blue colored regions at the boundaries in Fig. 8 . However, once l d > R e , the l d boundary layer is getting too wide compared with R e , and the searcher starts to miss targets. This explains the breakdown of the maximum at α = 2 from l d = 1 to l d = 1.24 by a breakdown of the pac-man effect. Note that for l d ≤ R e the pac-man effect equally applies to arrival and passage search, which elucidates why here the maxima are very similar for both efficiencies (for l d = 1 the maxima look even identical). But this mechanism does not clarify why the maxima when l d ≤ 1 become minima for 1 < l d < 4. However, there is a third length scale that plays a crucial role at the boundary, which is the cut-off dis-tance d c = 2 specific to replenishing targets. If we now consider the motion approximately perpendicular to the wall, a searcher must leave a target region of approximately 2d c = 4 before the target can replenish. But this defines yet another boundary layer of respective width that a searcher should leave as quickly as possible to benefit from the revisiting target mechanism, here induced by anti-persistence due to the reflecting boundaries. Indeed, a target will not replenish at all for a searcher bouncing multiple times perpendicularly to the wall with jump length l d < 4. We call this the bouncing fly effect, as this is similar to a fly hitting a window many times, sometimes at almost the same spot, by trying to escape. While this effect is always present when l d < 4, for l d < 1, it seems to be dominated by the above pac-man effect. However, as pac-man breaks down for l d > 1, our explanation of the minimum at α = 2 for 1 < l d < 4 is that here a searcher moves in a boundary layer that is deficient for revisiting target search induced by anti-persistence due to the reflecting boundaries.

This effect minimising efficiencies breaks down again when l d > 4, as then for the first time a searcher can jump over a distance larger than the cut-off replenishing target region. This now reactivates the revisiting target mechanism as a beneficient search strategy since during a jump, a target replenishes and is available again to be found. This explains why for l d = 4, we have again maximal efficiencies at α = 2. In the transition from smaller to larger α values at l d = 4, there even appears to be a slight minimum around α ≃ 1.7 in η * P while, conversely, the corresponding efficiency for FP search in Fig. 2 (b) is rather maximal in the same parameter region. This demonstrates again the sensitivity of optimal search on the variation of both internal and external parameters of the whole process.

Similarly, the minima in Fig. 7 for the l d = 1 curves around α = 2 can be explained. Since in this case the targets are non-replenishing, the boundary layer of width l d will become depleted of targets, which is detrimental to repeated pac-man search success by generating very small efficiencies. Minima for l d = 2 around α = 2 existed before already both in Fig. 6 and in Fig. 7 , and there is no other mechanism in place that could yield any larger efficiency here.

Finally, to confirm the impact of the boundary conditions on the efficiencies, Fig. 9 shows again η * A/P for reflecting boundaries, cf. η * A/P at l d = 1, 2 in Fig. 6 , in comparison to the ones for periodic boundaries. One can see that the periodic boundary conditions eliminate the maxima and minima at α = 2 in all cases. This unambiguously demonstrates that all these extrema are indeed due to boundary effects, as argued above. We furthermore remark that while for passage search ballistic motion with α = 2 now yields optimal efficiencies for replenishing targets in Fig. 9 (b) (if we neglect the fluctuations at small α values), this is not the case for arrival processes with l d = 2. This might be due to leapovers that become stronger for larger α values. 

In this work, we studied the efficiency of search generated by fBm in a random field of targets. In more general terms we explored the sensitivity of search succes on the specific setting and the parameters defining the search process, the environment and the interaction of both with each other. That way our essentially computational study suggests a conceptual framework that should apply to any theoretical description of a respective search problem, irrespective of whether one considers fBm, LWs or other types of motion. Figure 10 identifies important ingredients on which such a generic search setting depends. They can be broadly classified as intrinsic to a searcher by characterising its dynamics, or extrinsic to it by defining the search environment [7, 8] . The results also depend on the quantity by which search success is assessed in terms of statistical analysis [33] . For fBm we have investigated all of the conditions marked by (black) stars, as we will briefly summarise below by going through this figure. We emphasize that the picture put forward in Fig. 10 provides only a first sketch, which invites to be amended in future research.

The given properties of the searcher (light blue box to the left and associated tree structure in Fig. 10 ) define a key aspect of search. For search dynamics one needs to distinguish between the motion of the searcher and its associated modes of target detection. In our case, the motion was stochastic and governed by the two intrinsic parameters defining fBm, the jump length l d Eq. (8) and the exponent α of the mean square displacement Eq. (4). The latter determines in turn the memory of the process via the position autocorrelation function decay Eq. (3). In case of intermittent search, one has a separation of the search dynamics into local exploitation and longrange exploration [20, 50] . Concerning target detection, one distinguishes between cruise searchers that perceive a target while moving and saltatory foragers that only find targets after landing within the perception radius R e next to it [46] . The former relate mathematically to (first) passage problems of finding a target while the latter are (first) arrival, or (first) hitting problems [44] . Yet another distinction is whether one looks at the problem of finding only one target, which we modelled by a resetting procedure of the searcher to a random initial position, or many of them along a single trajectory, which we denoted as non-resetting after visiting a target.

These properties in turn determine how a search process is assessed by evaluating observables (turquoise box at the bottom in Fig. 10 ). Here we defined two generically different efficiencies depending on the statistical averaging applied, see Eqs. (9) [13] , respectively (18) [61] . The averaging was even obtained in a third way leading to yet another type of efficiency, see Eq. (11) [33] . We found that for non-Markovian processes with memory these three definitions yield very different values. Recent results showed that for LWs the dimensionality of the search process plays a crucial role for determining the values of efficiencies [37] , which however we did not explore in this work.

The search environment (green-grey box at the top) we modelled by a disordered, i.e., uniformly random distribution of homogeneous targets in the plane. Computer simulations are constrained to a finite area or volume determined by the length of the simulation box L as an extrinsic parameter. This quantity becomes a non-trivial parameter when an fBm searcher interacts with the walls of the system depending on the boundary conditions [92] . Here we considered primarily reflecting boundaries but compared some of our results to the situation of periodic ones. In some works, the impact of drift on (Lévy flight) search has been investigated [33, 56] . As targets (red-brown box to the right) we chose disks of radius R e , see Sec. II B, with a density measured by the mean target distance l f Eq. (6), which yields another two extrinsic parameters. As explained before, R e can be reinterpreted as the radius of perception of a searcher for finding point targets, hence it is ambiguous to categorise it as an intrinsic or extrinsic parameter. The disks were immobile, and we had many of them. One also needs to choose the type of resource [13] : A target may be nondestructive in the sense that it replenishes after having been found, or it is destroyed upon finding it and does not replenish. In the former case, one needs to introduce another extrinsic parameter, which in our case was the cut-off distance d c the searcher has to be away from the target after having found it, in order for it to replenish [37] , see Sec. IV.

We performed computer simulations to study the dependence of our three different search efficiencies on the above two intrinsic parameters by keeping the three extrinsic ones fixed. We also tested the impact of the boundary conditions. Our main results are summarised as follows:

1. For FA and FP search in a replenishing field of targets, we observe the existence of maxima in both search efficiencies for intermediate α values, 0.5 < α < 1.7, see Fig. 2 . The maxima are especially sharp in this regime if R e ≃ l d ≪ l f , and for FP with R e < l d < l f they accumulate around α ≃ 1.5. The latter seems to reflect an optimal sampling of the target space in-between over-and undersampling.

Similar results have been reported for FP search by LWs leading to the LFFH, i.e., that search strategies right between ballistic dynamics and normal diffusion are optimal to find sparse, replenishing targets in bounded domains [13] . The specific search scenario yielding the LFFH is marked by red stars in Fig. 10 . It is thus recovered as a special case in our over-arching conceptual framework.

2. For FA search, we encounter a 'paradise regime' when the jump length starts to exceed the mean target distance, R e ≪ l d ≃ l f . In this case, l d is large enough for the searcher to always find a target. This is essentially independent of the persistence in fBm if α is not too large or too small. In Fig. 2 (a) , this is represented by the efficiency curve flattening out, which implies that targets are found with maximal efficiency over a wide range of α values, 0.5 < α < 1.5.

3. In contrast to FP problems, FA processes exhibit leapovers when l d ≥ R e , that is, a searcher can jump over a target without finding it [33] . This mechanism diminuishes search efficiencies, as is clearly seen in Fig. 4 (b) for FP compared to FA, and also in Fig. 7 by comparing the efficiencies for arrival (saltatory) search to the ones for passage (cruise) search under variation of l d .

Optimising search by maximising efficiencies depends very much on the definition of efficiency that one chooses. This is demonstrated in Figs. 2, 6 and 7, which display results for the three different efficiency definitions Eqs. (9), (11) and (18), respectively. One can see that these three different efficiencies yield totally different results. 5 . We put forward a very simple analytical argument, which may be considered as a boiled-down version of the theory in Ref. [85] , that analytically reproduced the functional forms of the FP efficiencies under variation of l d and α, respectively, see Fig. 5 .

6. Subdiffusion can optimise search efficiencies for multi-target search along a trajectory of both arrival and passage processes in an area with reflecting boundaries when R e < l d ≪ l f , see Fig. 6 . This is due to a revisiting target mechanism, which in turn is generated by the anti-persistence in fBm for α < 1 [82] . The mechanism also holds for arrival processes under periodic boundary conditions in this regime of l d as shown in Fig. 9 (a).

7. Reflecting boundaries can generate very intricate memory effects in search governed by fBm [92] . This is represented by multiple transitions between maxima and minima in the efficiencies of multitarget search along a path for both arrival and passage processes, cf. again Fig. 6 for α > 1.5. We explained these variations microscopically by what we called pac-man and bouncing fly effects. There is a cross-link between these effects and the wellknown stickiness of active particles to walls due to self-driven persistent motion [87] .

We conclude that a seemingly simple problem of stochastic search in a random distribution of targets delivered highly non-trivial numerical results. Search turned out to be an extremely sensitive process, exhibiting all signatures of a complex system, where the whole is more than the sum of its single parts: Testing different search scenarios by differently combining search process, environment and their interaction with each other yielded entirely different results for respective search efficiencies, as is reflected in a high sensitivity of all results on variation of model parameters. A somewhat related sensitivity was already observed in foraging experiments [32, 41] . Our simulation results are thus in sharp contrast to claims of robust, universal optimal search strategies suggested by the LFFH and reported verifications of it by experiments conducted in the wild [10] . They are, however, fully in line with recent work limiting the range of validity of the LFFH [37, 40] by re-evaluating the theoretical model underpinning it. While here we did not investigate LWs, for the LW search that led to the LFFH the general framework is the same as the one summarised in Fig. 10 . The specific search situation that applies to the LFFH does have its place in this figure, but it only defines a small subset in it, as marked by the (red) crosses. We believe that this gives proper credit to the LFFH by adequately embedding it into a more general framework. For a similar search setting, we have found that maximisations of search efficiencies for fBm are not robust under variation of the search conditions and are thus not universally valid. In Refs. [37, 40] it was shown that the same holds to quite some extent for LWs.

To our knowledge, this is the first work where it has been explored how search by fBm depends on the variation of a large number of important conditions defining a specific search scenario. Our main results summarised above, in combination with Fig. 10 , reveal the need to investigate further search scenarios on a case by case basis. The most straightforward extension would be to study search under variation of the three external parameters of our model, which deliberately we kept fixed, i.e., the target size R e , the density of targets governed by l f , and the cut-off distance for the replenishing of targets d c . While we believe that we have identified generic boundary effects in our work, checking for their robustness under parameter variation would be interesting. More complicated problems would be search for moving targets, as studied to some extent already in Ref. [4] . Biologically, one might be interested in search restricted to the home range of a forager [108] . And investigating swarms of searchers is important for robotic applications [109] . We furthermore remark that stochastic dynamics with persistence became very prominent as models for active Brownian particles reproducing self-propelled biological motion [96, 97] . Within the latter context, it might be interesting to study biologically relevant search problems in more detail [52, 54] . Finally, a more comprehensive analytical description of fBm search [85, 105] going beyond our simple handwaving argument would be highly desirable. This theory should explain the intricate dependence of search efficiencies on the variation of model parameters and other settings as observed in our simulations.

Cover times of random searches

Intermittent search strategies

Evidence of Lévy walk foraging patterns in human hunter-gatherers

The theory of search games and rendezvous, International Series in Operations Research and Managment Science

Search research

Theory of Optimal Search

A movement ecology paradigm for unifying organismal movement research

Stochastic Foundations in Movement Ecology

Search for food of birds, fish and insects

The Physics of Foraging

Mathematical contributions to the theory of evolution -a mathematical theory of random migration

Lévy flight search patterns of wandering albatrosses

Optimizing the success of random searches

Revisiting Lévy flight search patterns of wandering albatrosses, bumblebees and deer

The random walk's guide to anomalous diffusion: A fractional dynamics approach

The Langevin Equation

Anomalous transport: Foundations and Applications

Models of anomalous diffusion in crowded environments

Anomalous diffusion models and their properties: nonstationarity, non-ergodicity, and ageing at the centenary of single particle tracking

Foraging success under uncertainty: search tradeoffs and optimal space use

Fractional Brownian motions, fractional noises and applications

Lévy flights and related topics

Strange kinetics

Lévy walks

Average time spent by Lévy flights and walks on an interval with absorbing boundaries

Dynamical robustness of Lévy search strategies

How many animals really do the Lévy walk?

Scaling laws of marine predator search behaviour

Environmental context explains Lévy and Brownian movement patterns of marine predators

Lévy walks evolve through interaction between movement and environmental complexity

Generalized Lévy walks and the role of chemokines in migration of effector cd8+ t cells

Spatiotemporal dynamics of bumblebees foraging under predation risk

Lévy flights do not always optimize random blind search for sparse targets

Understanding movements of organisms: it's time to abandon the Lévy foraging hypothesis

Liberating Lévy walk research from the shackles of optimal foraging

Current status and future directions of Lévy walk research

Inverse square Lévy walks are not optimal search strategies for d ≥ 2

Intermittent inversesquare Lévy walks are optimal for finding targets of all sizes

Comment on "inverse square Lévy walks are not optimal search strategies for d ≥ 2

Reply to "comment on 'inverse square Lévy walks are not optimal search strategies for d ≥ 2

Ergodicity breaking and lack of a typical waiting time in area-restricted search of avian predators

A guide to first-passage processes

First-Passage Phenomena and Their Applications

First passage and first hitting times of Lévy flights and Lévy walks

Universality classes of hitting probabilities of jump processes

Efficient or inaccurate? Analytical and numerical modelling of random search strategies

Search strategies of foraging animals

Assessing Lévy walks as models of animal foraging

Scale-free foraging by primates emerges from their interaction with a complex environment

First-passage times in multiscale random walks: The impact of movement scales on search efficiency

Optimal Lévy-flight foraging in a finite landscape

Active Brownian particles and run-andtumble particles separate inside a maze

The topography of the environmkent alters the optimal search strategy for active particles

Active Brownian particles moving in a random Lorentz gas

Exact spatiotemporal dynamics of confined lattice random walks in arbitrary dimensions: A century after Smoluchowski and Pólya

Space-fractional Fokker-Planck equation and optimization of random search processes in the presence of an external bias

Search reliability and search efficiency of combined Lévy-Brownian motion: Long relocations mingled with thorough local exploration

Comparison of pure and combined search strategies for single and multiple targets

Lévy strategies in intermittent search processes are advantageous

Animal search strategies: a quantitative random-walk analysis

Fractal reorientation clocks: Linking animal behavior to statistical patterns of search

Stochastic optimal foraging: Tuning intensive and extensive dynamics in random searches

First Steps in Random Walks: From Tools to Applications

The wiener spiral and some other interesting curves in hilbert space

Correlation theory of processes with random stationary nth increments

Selfsimilar processes

Fractional Brownian motion and motion governed by the fractional Langevin equation in confined geometries

Single particle tracking in systems showing anomalous diffusion: the role of weak ergodicity breaking

Fractional Brownian motions: memory, diffusion velocity, and correlation functions

Viscoelastic subdiffusion: Generalized langevin equation approach

Single-particle tracking data reveal anticorrelated fractional Brownian motion in crowded fluids

What information is contained in the fluorescence correlation spectroscopy curves, and where

Elucidating the origin of anomalous diffusion in crowded fluids

In vivo anomalous diffusion and weak ergodicity breaking of lipid granules

Probing the interior of living cells with fluorescence correlation spectroscopy

Anomalous subdiffusion is a measure for cytoplasmic crowding in living cells

First passage in an interval for fractional Brownian motion

Scale-free animal movement patterns: Lévy walks outperform fractional Brownian motions and fractional lévy motions in random search scenarios

Cellular automaton models for timecorrelated random walks: Derivation and analysis

Analyzing the impact of sars cov-2 on the human proteome

Maximum of a fractional Brownian motion: Probabilities of small values

Sampling the cell with anomalous diffusion -The discovery of slowness

First passage behaviour of fractional Brownian motion in twodimensional wedge domains

First passage times for a tracer particle in single file diffusion and fractional Brownian motion

Mean first-passage times of non-markovian random walkers in confinement

Active Brownian particles

Active particles in complex and crowded environments

Active Brownian particles near straight or curved walls: Pressure and boundary layers

Boundaries can steer active Janus spheres

Run-and-tumble dynamics of self-propelled particles in confinement

How to capture active particles

Fractional Brownian motion in a finite interval: correlations effect depletion or accretion zones of particles near boundaries

Reflected fractional Brownian motion in one and higher dimensions

Superdiffusive dispersals impart the geometry of underlying random walks

Processes with long-range correlations: Theory and applications

How far from equilibrium is active matter?

Langevin dynamics driven by a telegraphic active noise

Normal and anomalous fluctuation relations for Gaussian stochastic dynamics

Modeling persistence in hydrological time series using fractional differencing

Stochastic simulation with a view towards stochastic processes

Tests for Hurst effect

Self-similarity through high-variability: statistical analysis of ethernet lan traffic at the source level

Simulation of fractional Brownian motion with conditionalized random midpoint displacement

On spectral simulation of fractional Brownian motion

Universal first-passage statistics in aging media

Geometrycontrolled kinetics

From first-passage times of random walks in confinement to geometry-controlled kinetics

From micro-to-macro: How the movement statistics of individual walkers affect the formation of segregated territories in the territorial random walk model

Random walks in swarm robotics: An experiment with kilobots

Acknowledgements: This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -Projektnummer 163436311 -SFB 910. RK thanks the SFB for a Mercator Visiting Professorship, which initiated this project and during which most parts of it have been performed. He is also grateful for funding from the Office of Naval Research Global and acknowledges an External Fellowship from the London Mathematical Laboratory. All three authors thank two anonymous referees for their very helpful comments that led to a profoundly improved version of their manuscript.