First passage properties of asymmetric L\'evy flights

L\'evy Flights are paradigmatic generalised random walk processes, in which the independent stationary increments---the"jump lengths"---are drawn from an $\alpha$-stable jump length distribution with long-tailed, power-law asymptote. As a result, the variance of L\'evy Flights diverges and the trajectory is characterised by occasional extremely long jumps. Such long jumps significantly decrease the probability to revisit previous points of visitation, rendering L\'evy Flights efficient search processes in one and two dimensions. To further quantify their precise property as random search strategies we here study the first-passage time properties of L\'evy Flights in one-dimensional semi-infinite and bounded domains for symmetric and asymmetric jump length distributions. To obtain the full probability density function of first-passage times for these cases we employ two complementary methods. One approach is based on the space-fractional diffusion equation for the probability density function, from which the survival probability is obtained for different values of the stable index $\alpha$ and the skewness (asymmetry) parameter $\beta$. The other approach is based on the stochastic Langevin equation with $\alpha$-stable driving noise. Both methods have their advantages and disadvantages for explicit calculations and numerical evaluation, and the complementary approach involving both methods will be profitable for concrete applications. We also make use of the Skorokhod theorem for processes with independent increments and demonstrate that the numerical results are in good agreement with the analytical expressions for the probability density function of the first-passage times.


Introduction
Normal Brownian motion described by Fick's second law, the diffusion equation, is characterised by the linear time dependence x 2 (t) t of the mean squared displacement (MSD) [1]. Deviations from this Fickean time dependence typically occur in the power-law form x 2 (t) t κ of the MSD [2,3]. Depending on the value of the anomalous diffusion exponent κ we distinguish between subdiffusion for 0 < κ < 1, normal, Fickean diffusion for κ = 1, superdiffusion for 1 < κ < 2, and ballistic, wavelike motion for κ = 2. The range κ > 2 is sometimes referred to as hyperdiffusion. The theoretical description of anomalous diffusion phenomena of physical particles (passive or active) often requires a radical departure from the classical formalism for Brownian motion. Namely, effects of energetic or spatial disorder, collective dynamics, or nonequilibrium conditions need to be addressed with more complex approaches [2,4]. For instance, fractional Brownian motion [5] is a process in which the Langevin equation is driven by Gaussian yet power-law correlated noise (fractional Gaussian noise) effecting both sub-and superdiffusion. The generalised Langevin equation [6] includes a memory integral with a kernel, that balances the input fractional Gaussian noise and effects a thermalised process. Processes with explicitly time or position dependent diffusion coefficients such as scaled Brownian motion [7] or heterogeneous diffusion processes [8], respectively, also lead to sub-and superdiffusion. Diffusion on fractals [9] due to the fact that the particle in the highly ramified environment often has to back-track its motion, has similar characteristics as subdiffusive fractional Brownian motion. We finally mention the continuous time random walk model, in which the standard Pearson walk was generalised to include continuous waiting times [10]. When the distribution of waiting times becomes scale-free, with diverging characteristic waiting time, the continuous time random walk process is subdiffusive [11,12]. Conversely, when the continuous time random walk has a finite characteristic waiting time but is equipped with a scale-free distribution of jump lengths with power-law asymptote λ(x) ∼ |x| −1−α (0 < α < 2) the resulting process is a "Lévy flight" (LF). As in this case the variance of the process diverges, diffusion can be characterised in terms of rescaled fractional order moments |x(t)| η 2/η t 2/α with 0 < η < α [3,13,14,15,16,17]. Mathematically, asymptotic power-law forms of the jump length distribution can be explained by the generalised central limit theorem [2,18,19,20], which gives rise to the much higher likelihood for extremely long jumps [21,22,23] in comparison to conventional Pearson random walks.
Lévy stable laws play a crucial role in the statistical description of scale invariant stochastic processes [21,24] not only in physical contexts but also in biological, chemical, geophysical, sociological, economical or financial systems, among others. In physics Lévy statistics were demonstrated to explain deviations of complex systems from the Gaussian paradigm, inter alia, for the power-law blinking of nano-scale light emitters [25], diffusive transport of light [26], photons in hot atomic vapours [27], tracer particles in a rotating flow [28], passive scalars in vortices in shear [29], anomalous diffusion in disordered media [2], weakly chaotic and Hamiltonian systems [30,31], in the divergence of kinetic energy fluctuations of a single ion in an optical lattice [32], fluctuations in the transition energy of a single molecule embedded in a solid [33], in the interaction of two level systems with single molecules [34], the distribution of random single molecule line shapes in low temperature glasses [35], the diffusion of a collection of ultra-cold atoms and single ions in optical lattices [36], but also in reaction diffusion systems [37,38]. Lévy statistics was observed experimentally in tokamak and stellarator fusion devices [39,40,41]. It was also shown that the phenomenon of L-H transitions observed in the stellarator is accompanied by the crossover from Lévy to Gaussian fluctuation statistics [42]. Numerous examples for Lévy statistics exist in the dynamics of plasmas, including the anomalous transport in magnetic confinement [43,44,45,46,47], the dynamics of a charged particle in a plasma [48,49], anomalous transport of ions and electrons in solar winds [50], nonlocal transport in plasma turbulence [51,52,53,54,55], and heat transport in magnetised plasmas [56,57].
In the biological sciences, many organisms from bacteria to humans use Lévy stable relocation statistics in their search for resources [58], tracer motion in living cells [59,60,61,62], or the superdiffusive motion of bacteria within a swarm [63]. On long, fast-folding polymers the search process of a binding molecule is based on Lévy motion [64,65,66,67]. In geoscience paleoclimate time series show signatures of Lévy noise [68], earthquake statistics exhibit distinct power-law behaviour [69], as well as tracer plumes in heterogeneous aquifers [70,71,72], and the transport of ensembles of particles on the Earth surface [73]. The mechanisms of the worldwide spread of infectious diseases [74], pollen dispersal by bees [75] and human mobility patterns and social interactions revealed by tracing mobile phones or following banknotes [76,77,78,79,80,81] also reveal Lévy statistics. Evidence of Lévy stable laws was also unveiled in the human cognition for the retrieval dynamics of memory [82] and in human mental search [83,84,85], as well as search in multiplex networks [86]. The optimal search patterns of robots were shown to be based on Lévy stable laws [87]. In finance and economical contexts [88,89,90,91] Lévy statistics govern the distribution of trades. A particular area in which Lévy relocation statistics has been widely explored is movement ecology [92]. Search patterns of foraging animals [93] that follow Lévy statistics include marine predators [94,95], albatross birds [96], Agrotis segetum moths [97], fruit flies (Drosophila melanogaster) [98], bumblebees [99], jellyfish [100], goats [101,102], immune function of human T cells [103] and human hunter-gatherers [104]. We hasten to note that in the context of animal movement there exist some debates on the predominance of Lévy search patterns, especially the disqualification of Lévy statistics for albatrosses [105,106] in [107] became a strong argument against the LF hypothesis. However, there is strong evidence that for individual albatross birds LFs are indeed a real search pattern [108]. Moreover, [109] reported that for spider monkeys the foraging pattern is deterministic, mussel movements are rather multimodal [110,111] and black bean aphids individually move in a predominantly diffusive manner [112].
The efficiency of search processes is typically benchmarked by the time it takes the searcher to reach a certain region in space. One relevant measure is therefore the first-passage time, quantifying the time it takes from the original position to first cross a point located at a given distance away. For instance, the first-passage time in a financial time series could be defined by a preset increase or decrease in the price of a given stock.
Once this threshold value is reached, the stock is sold or bought. Similarly, we could talk about the instant of time when foraging animals first randomly locate a resourcerich area away from their original location. Such first-passage times in a stochastic search process will vary from realisation to realisation, and can be quantified by the first-passage time density ℘(t). While the mean first-passage time t = ∞ 0 t℘(t)dt can capture some aspects of this dynamics, ‡ the full information encoded in ℘(t) provides significant additional insight [114,115,116]. Here we study the first-passage properties for a general class of α-stable Lévy laws. We go beyond previous approaches [113,117,118,119,120,121,122,123,124,125,126,127,128] focusing on symmetric and one-sided α-stable relocation distributions and consider α-stable laws with arbitrary asymmetry in semi-infinite and bounded domains. Our approach is based on the convenient formulation of LFs in terms of the space-fractional diffusion equation. We derive these integro-differential equations for LFs based on general asymmetric α-stable distributions of relocation lengths in finite domains, and thus go beyond studies of the exit time and escape probability in bounded domain for symmetric LFs [129,130,131]. An important aspect of this study is that we complement our results with numerical analyses of the (stochastic) Langevin equation for LFs and show how both approaches complement each other.
The paper is organised as follows. In section 2 we define Lévy stable laws and the associated fractional diffusion equation. In section 3 we set up our numerical model for the fractional diffusion equation and the associated Langevin equation. Moreover, a comparison between the numerical method and α-stable distributions for symmetric and asymmetric density functions is presented. Section 4 then presents the numerical results for the survival probability and first-passage time density for both symmetric and asymmetric probability density functions. Our findings are compared with results derived from the Skorokhod theorem for symmetric, one-sided, and extremal two-sided stable distributions. We draw our conclusion in section 5. In the Appendices, we present details of several derivations. is any real number, and the phase factor ω is defined as In the real space-time domain the PDF of the α-stable distribution can be expressed via elementary functions for the following three cases: (a) Gaussian distribution, α = 2, β irrelevant: (b) Cauchy distribution, α = 1, β = 0: In physics, the Cauchy distribution is also often called a Lorentz distribution.
(c) Lévy-Smirnov distribution, α = 1/2, β = 1: Schneider reported the representation of general Lévy stable densities in terms of Fox H-functions [133,134,135]. Somewhat simpler representations for rational indices α and β are given in [136,137]. More information on Lévy stable densities and their parametrisation are provided in Appendix A.
Physically, the parameter µ accounts for a constant drift in the system. In this paper we consider the first-passage process in the absence of a drift, thus in what follows we set µ = 0. Let us first consider the case α = 1 and −1 ≤ β ≤ 1. The corresponding space-fractional diffusion equation for the PDF α,β (x, t) then reads where D α x is the space-fractional derivative operator, that we compose of −∞ D α x and x D α ∞ , the left and right hand side space-fractional operators, respectively. We use the Caputo form of the operators defined by (n − 1 < α < n) [138] and L α,β and R α,β are the left and right weight coefficients, defined as [52,53] The Fourier transforms of the operators (6a) and (6b) have the forms [138] where ÷ defines the Fourier transform pairs. For the case α = 1 and β = 0 we have L 1,0 = R 1,0 = 1/π [139], and instead of equation (5) we find where H is the Hilbert transform in terms of the principle value integral − . In Fourier space the Hilbert transform has the simple form In what follows we do not consider the particular case α = 1, β = 0 since it cannot be described in terms of a space-fractional operator. For all other choices of the parameters by substitution of relations (8a), (8b), and (11) into equation (4) we recover the characteristic function (1) of the α-stable process after Fourier transform.

Numerical scheme
To determine the first-passage properties of α-stable processes we will employ different numerical schemes based on the space-fractional diffusion equation and the Langevin equation for LFs. We here detail their implementation.

Diffusion description
There are several numerical methods to solve space-fractional diffusion equations, such as the finite difference [140,141] and finite element [142,143,144] methods as well as the spectral method [145,146]. In this paper we use the finite difference method that uses differential quotients to replace the derivatives in the differential equations. The domain is partitioned in space and time, and approximations of the solution are computed. Due to causality we use forward differences in time on the left hand side of equation (4), where f j i = f (x i , t j ), x i = (i − I/2)∆x, and t j = j∆t, where ∆x and ∆t are step sizes in position and time, respectively. The i and j are non-negative integers, i = 0, 1, 2, . . . , I, and further x 0 = −L, x I = L, and ∆x = 2L/I. Similarly, j = 0, 1, 2, . . . , J − 1, t 0 = 0, t J = t, and ∆t = t/J. Absorbing boundary conditions for the determination of the first-passage event imply f j 0 = f j I = 0 for all j. The integrals on the right hand side of equation (4) are discretised as follows. For 0 < α < 1, for the left side derivative, and for the right side derivative. Thus the idea is to approximate only the derivative by the differences. The integral kernel is then calculated explicitly. For the estimation of the error in this scheme we refer to Appendix C. For the case 1 < α < 2 we use the central difference approximation, namely, for the left side, and for the right side. For the special case α = 1, β = 0 by using the discrete Hilbert transform [147] we approximate the derivative in space, namely For further details of the numerical scheme we refer the reader to Appendix B. To improve the stability we use the Crank-Nicolson method. By substitution of equations (12) to (15) into equation (4) we obtain in which the coefficients A and B have matrix form with dimension (I + 1) × (I + 1) and j = 0, 1, 2, . . . , J − 1. In the numerical scheme the initial condition f (x, 0) = δ(x) is approximated as For the setup used in our numerical simulations, see section 4 below, the initial point is the distance of x(0) from the right interval boundary (see figure 3). For this case we implement the initial condition as In the next step the time evolution of the PDF is obtained by using the absorbing boundary conditions f j 0 = f j I = 0 for all j and applying to the matrix coefficients A and B.

Langevin dynamics
The fractional diffusion equation (4) can be related to the Langevin equation [14,123,148] d dt where ζ(t) is white Lévy noise characterised by the same α and β parameters as the space-fractional operator (5) and with unit scale parameter. The Langevin equation (18) provides a microscopic (trajectory-wise) representation of the space-fractional diffusion equation (4). Therefore, from an ensemble of trajectories generated from equation (18) it is possible to estimate the time dependent PDF whose evolution is described by equation (4). In numerical simulations Lévy flights can be described by the discretised form of the Langevin equation where ζ stands for the α-stable random variable with a unit scale parameter [18,149] and the same index of stability α and skewness β parameters as in equation (18). Relation (19) is exactly the Euler-Maruyama approximation [150,151,152] to the general αstable Lévy process. From the trajectories x(t), see equations (18) and (19), it is also possible to estimate the first-passage time τ as From the ensemble of first-passage times it is also possible to estimate the survival probability S(t), the complementary cumulative density of first-passage times. More precisely, the initial condition is S(0) = 1, and at every recorded first-passage event at time τ i , S(t) is decreased by the amount 1/N , where N is the number of records of first-passage times. If a given estimation of the first-passage time is recorded k times the survival probability is decreased by k/N . For a finite set of first-passage times there exists a small fraction of very large first-passage times. Therefore, this estimation becomes poorer with increasing time t. In the next section we present a comparison between the numerical solution of equation (4) and the α-stable probability laws with characteristic function (1).

Comparison with α-stable distributions
We now show that the difference scheme for the space-fractional diffusion equation provides excellent agreement with the theoretical results for the shapes of α-stable probability densities. To this end we use a MATLAB code to obtain the inverse Fourier transform of the characteristic function [153]. This programme employs Zolotarev's so-called M-form for the parametrisation of α-stable distributions with parameters α, β M , µ M , and K M α , while in the main text we use the A-form with parameters α, β A , µ A and K A α [154]. Thus, along with the code we use the corresponding change of the distribution's parameters, see Appendix A and, in particular, equation (A.8) for details. We use absorbing boundary conditions and a finite domain with half length L = 100 in one dimension, the initial condition is a Dirac δ-function located at x(0) = 0. The probability density function for β = 0 and different sets of the index of stability α at t = 1 is displayed in figure 1 on the left. The tails display the correct power-law scaling. In the right panel of figure 1 we show the PDF for α = 0.5 and β = 0 at different times. The insets focus on the central part of the PDFs. In all cases and over the entire plotted range the agreement between the numerical solutions and the theoretical densities is excellent.
3.3.2. Asymmetric α-stable distributions. Asymmetric stable distributions with nonzero values of the skewness β may occur in various situations, for instance, when in a random walk the left and right diffusion coefficients are different. In figure 2 (top) we plot the PDF with α = 1/2 at t = 1 for different values of the skewness parameter β.
On the left side in the main panel the negative side of the tails is shown, on the right side we display the positive side of the tails. Figure 2 (bottom) analogously shows the PDF for α = 1.3 and different β at time t = 1. Note that for β = 1 and 0 < α < 1 the PDF is completely one-sided on the positive axis and does not possess a left tail.
By comparison we see that again the numerical scheme for solving the spacefractional diffusion equation produces solutions that are in very good agreement with the numerical results for the α-stable laws. In the following we study the first-passage properties of random walk processes with α-stable jump length distribution obtained from two numerical methods: the space-fractional diffusion equation and the Langevin approach. We also compare the numerical method for solving the space-fractional diffusion equation with results following from Skorokhod's theorem.

Survival probability and first-passage properties
The survival probability and the first-passage time are observable statistical quantities characterising the stochastic motion in bounded domains with absorbing boundary conditions. In the following we investigate the properties of the survival probability and the first-passage time density in a finite interval for symmetric and asymmetric α-stable laws underlying the space-fractional diffusion equation. To this end we use the setup shown in figure 3, in which absorbing boundaries are located at −L and L, and the initial point of the initial δ-distribution is located the distance d away from the right boundary.
The probability that at time t the random walker is still present in the interval [−L, L] is defined as the survival probability [155] and the first-passage time PDF is given by the negative time derivative, Therefore, in Laplace domain with initial condition S(0) = 1 the relation between the survival probability and the first-passage time reads We now first consider the survival probability for symmetric α-stable laws in a semiinfinite and finite interval and demonstrate how the asymptotic properties change with the length of the interval. Furthermore, we compare the results obtained from the numerical difference scheme in section 3.1 with the Langevin equation approach, before embarking for the study of asymmetric α-stable laws.

First-passage processes for symmetric α-stable laws
Here we study the properties of α-stable processes in domains restricted by one or two boundaries. In figure 4 we show the survival probability for different α and two Figure 3. Schematic of the setup used in our approach: in the interval of length 2L the initial probability density function is given by a δ-distribution located at where d is the distance from the right boundary. At both interval boundaries we implement absorbing boundary conditions, that is, when the particle hits the boundaries or attempts to move beyond them, it is absorbed.   figure 4 clearly show an exponential decay (in analogy to the escape dynamics of Lévy flights from a confining potential [156,157]). For the short interval the exponential behaviour sets in almost immediately on the linear time axis in the plot, while for the longer interval a crossover behaviour is visible, as we will see below. Interestingly, we see in figure 4 that the trend of decay reverses with respect to the stable index α. To understand this behaviour of the survival probability we use the following approximation of the survival probability of symmetric Lévy flights in finite intervals, where the mean first-passage time is given by [158,159] (24) and (26) agree very well with the numerical results for L = 1: an increase of the interval length L leads to a decrease of τ x(0) −1 for larger α (right panel of figure 5). For a semi-infinite domain with absorbing boundary condition it is well known that the first-passage time density for any symmetric jump length distribution in a Markovian setting has the universal Sparre Andersen asymptotic ℘(t) t −3/2 (and thus S(t) t −1/2 [155,160,161]. This is exactly our setting here for the symmetric case with β = 0, and the Sparre Andersen universality was consistently corroborated for symmetric LFs in a number of works, inter alia, [117,123,124,125].
In figure 6 we study what happens at intermediate times in the case of a finite interval, before the terminal exponential shoulder cuts off the survival probability, as shown in figure 4. On the log-log scale of figure 6, we indeed recognise a transient power-law scaling with the universal Sparre Andersen exponent 1/2 for the survival probability. The onset of the hard exponential cutoff is shifted to longer times with increasing interval size L, in which, on average, it takes the particles longer to explore the full extent of the domain. This, of course, is fully consistent with results for normal diffusion as well as continuous time random walk subdiffusion subordinated to regular random walks [114,115,116,155,162], compare also the discussion of the area filling dynamics of LFs [163]. Moreover, we see that the results from numerical solution of the fractional diffusion equation and simulations of the Langevin equation almost perfectly agree with each other. The lines without symbols in the top left of figure 6 correspond to cases when the numerical approach based on the fractional diffusion equation did not converge well. We note that one has to increase the value of L with decreasing α in order to meet the Sparre-Andersen scaling for a semi-infinite interval. This is intuitively clear, as smaller α enhances the likelihood of longer jumps and thus effects a higher probability of interaction with the interval boundaries at fixed L.
In figure 7 for the interval size L = 100 we show the survival probability in the left panel along with the the first-passage time density in the right panel, for different α and β = 0. Consistently, the transient Sparre Andersen scaling is passed on from the powerlaw exponent 1/2 for the survival probability to the exponent 3/2 of the first-passage density.
In the theory of a general class of Lévy processes, that is, homogeneous random processes with independent increments, there exists a theorem, that provides an analytical expression for the PDF of first passage times in a semi-infinite interval, often referred to as the Skorokhod theorem [132,164]. Based on this theorem the asymptotic expression for symmetric α-stable laws, the first-passage time PDF is (Appendix D) [125] which specifies an exact expression for the prefactor of the Sparre Andersen power-law. For Brownian motion (α = 2), the PDF for the first-passage time has the well known Lévy-Smirnov form [165] that also emerges from the Skorokhod theorem in the limit α = 2. Equation (28) is exact for all times [165,166], and apart from the Sparre Andersen law ℘(t) t −3/2 it includes the hard short time exponential cutoff reflecting the fact that it takes the particle a typical time ∝ d 2 /K 2 to reach the absorbing boundary.

First-passage processes for asymmetric α-stable laws
The case of asymmetric α-stable laws is mathematically more involved and also has been less well studied numerically. We now present results for the survival probability and first-passage time PDF for different skewness parameters β, addressing first the special cases of completely one-sided and extremal two-sided laws.
4.2.1. One-sided α-stable laws. One-sided α-stable laws with α ∈ (0, 1), β = 1 satisfy the non-negativity of their increments. Physically, such laws are appropriate for jump processes that always move in the same direction, for instance, as a generalisation of shot noise. Applying the Skorokhod theorem to this case one can show that the firstpassage time PDF in the permitted interval α ∈ (0, 1) has the exact analytical solution (Appendix D) [125] and where M α is the Wright M -function (also called Mainardi function) [138,167] with series representation (E.9) and asymptotic exponential decay (E.11). At sufficiently long times the first-passage time PDF reads where we used the coefficients From equation (31) we see that for smaller α the first-passage time density decays faster which is intuitively clear since Lévy flights with smaller α feature longer jumps in the positive direction. The value of ℘(t) at t = 0, demonstrates that, in contrast to the Gaussian case, the probability density does not vanish at t = 0, indicating the possibility of immediate escape.
Using equation (22) the survival probability for 0 < α < 1, β = 1 can be expressed exactly in terms of the Wright function W a,b (see equation (E.7)) or its series expansion, In the limit α = 1/2 this expression can be reduced to the simple form and the corresponding first-passage time density is the half-sided Gaussian [124,168]

4.2.2.
Extremal two-sided α-stable probability laws. Stable probability laws with stability index 1 < α < 2 and skewness β = 1 or β = −1 are called extremal twosided skewed α-stable laws. When β = 1 the PDF of an α-stable random variable has a positive power-law tail x −1−α , and a negative tail which is lighter than Gaussian [169], see figure 2. For β = −1 the properties of the tails are exchanged. In Appendix D (see equation (D.54)) by applying the Skorokhod theorem it is shown that for β = −1 the PDF of the first-passage time for extremal two-sided stable probability laws has the exact form in terms of the Wright M -function M 1/α . In the limit α = 2 we recover the PDF (28) for a Gaussian process. Moreover by using equation (22) the survival probability can be transformed to Equation (37) has the following limiting behaviours: at short times t → 0 corresponding to the asymptotic of large argument in the Wright function, by using the asymptotic expression (E.11) we find where At long times, t → ∞, consistent with the result in [124].
For the extremal two-sided α-stable probability laws with index 1 < α < 2 and skewness β = 1, by using the Skorokhod theorem the PDF of first-passage times has the following series representation (see Appendix D, equation (D.72)) For α = 2 we again consistently recover result (28). The asymptotic behaviour of equation (42) at long times is and in the limit t → 0, we find the finite value By using the relation between the survival probability and the first-passage time PDF in Laplace space (equation (23)) and applying the inverse Laplace transform we obtain a series representation for the survival probability for extremal two-sided α-stable probability laws (1 < α < 2, β = 1) in the form The first-passage time PDFs for extremal two-sided α-stable probability laws are displayed in figure 9.

4.2.
3. α-stable probability laws in general asymmetric form. We finally study the firstpassage behaviour for asymmetric Lévy stable laws of arbitrary skewness β, again based on the comparison between the numerical solution of the space-fractional diffusion equation and simulations of the Langevin approach for different stable index α. The results are shown in log-log scale in figures 10 and 11 in the range 1 < α < 2. Figure  10 depicts the cases β = 1 and β = −1, while figure 11 shows the cases β = 0.5 and β = −0.5. As can be seen from both figures, a positive value of the skewness leads to a slower decay (shallower slope) than for the Gaussian case, and opposite for negative values of β. Indeed, this behaviour is not immediately intuitive, as a positive skewness means that the stable law has a longer tail on the positive axis. However, what matters for the short and intermediate first-passage time scales are values of the stable law around the most likely value, and for positive skewness this is located on the negative axis (compare bottom panels in figure 2). Thus, Lévy flights with a positive skewness experience an effective drift to the left, in our setting away from the absorbing boundary. From the Skorokhod theorem for α ∈ (0, 1) and β ∈ (−1, 1), α = 1 and β = 0, as well as α ∈ (1, 2] and β ∈ [−1, 1] we obtained the power-law decay (see Appendix D.6) where the positive parameter ρ is defined as Following [170] (page 218) we write this in a form with a general β. This is the direct generalisation of the classical Sparre Andersen universality for asymmetric stable jump length distributions. It is easy to check that this result reduces to the known cases for vanishing skewness. We note that the inapplicability of the Sparre Andersen law for asymmetric jump length distributions was already pointed out by Spitzer ([171], page 227). The analytical predictions from relations (46) and (47) are in excellent agreement with the data shown in figures 10 and 11.

Discussion and unsolved problems
LFs are Markovian stochastic processes that are commonly used across disciplines as models for jump processes that exhibit the distinct propensity for very long jumps. The scale-free nature of the underlying, Lévy stable PDF of jump lengths with its heavy power-law tail has been shown to effect a more efficient random search strategy than the more conventional Brownian search processes. We here combined numerical inversion methods of the solution of the space-fractional diffusion equation and simulation of the Langevin equation fuelled by α-stable white noise to quantify the first-passage dynamics of LFs with a general asymmetric jump length PDF. In particular, we demonstrated that in all cases both approaches are in excellent agreement. As both approaches have advantages and disadvantages, it is very useful to have available two equally potent methods. In addition, we verified the crossover to an exponential behaviour of the (0,1) 1 Andersen theorem are no longer fulfilled, we derived the analytical behaviour from the Skorokhod theorem for specific values of the skewness. In the general case the direct extension of this analytical law was shown to be fully consistent with the numerical and simulations results. The results obtained here will be of use in applications, as these typically are involved with search processes and thus measure first-passage times. Concurrently, these results also further complete the mathematical theory of Lévy stable processes.
The first-passage time properties of general α-stable probability laws are summarised in table 1. For the known cases we include the references to the relevant equations of the exact PDF as well as the asymptotic prefactor. Some special cases are included, as well. Those cases with previously known results refer to the relevant references.
It is possible to extend the studied setup to higher dimensions [175,176,177,178]. In this case, the scalar noise term ζ(t) in the Langevin equation (18) for x(t) is replaced § Note that the result (41) differs from that of [172] by a factor which appears due to the use of two different forms of the characteristic function for the α-stable process. by multivariate Lévy noise in the higher-dimensional Langevin equation for x(t). Here, multivariate α-stable variables are characterised by a spectral measure defined on the unit circle [18]. For the numerical scheme of the multi-dimensional space-fractional diffusion equation we refer to [179,180,181,182,183,184,185,186]. We note that to the best knowledge of the authors no multi-dimensional generalisation of Skorokhod's theorem exists. Thus, the extension of the analytical and numerical approaches presented here to higher dimensions represents a challenging problem requiring further studies.
Generally the formulation of non-local and/or correlated stochastic processes is not always an easy task and, in some cases, still not fully understood. Apart from LFs, we may allude to the debate on the formulation and solution of boundary value problems for fractional Brownian motion, a process fuelled with Gaussian yet longrange correlated noise [187,188,189]. For LFs, in addition to the results obtained here it will be interesting to generalise the results obtained for symmetric α-stable laws in the presence of an external drift in [113]. Similarly, it will be of interest to investigate the PDF of first arrival times, related to the probability of hitting a small target in an otherwise unbounded environment, for the general case of asymmetric Lévy stable laws.

Appendix A. Parametrisation of characteristic function for α-stable processes
There are several forms for the parametrisation of α-stable laws appearing in literature, basically because of historical reasons. Each form might be useful in a particular situation. For example, one of them is preferable for analytical calculations, whereas the other ones can be more convenient for numerical purposes or for fitting of data. Following the exposition of the various forms of stable laws in [154,190], we here present four parameterisation forms for the characteristic functions.
In the main text we use the A-form of the characteristic function, In this paper we exclude the case α = 1 and β = 0. The B-form is helpful from an analytical point of view, it is given bŷ where (for α = 1) and K(α) = α − 1 + sgn(1 − α). The parameters have the same domains of variation as in the A-form, The M-form is used in numerical simulations and readŝ The domain of variation of the parameters in the A-and M-forms is the same and connected by the following relations Finally, the Z-form is represented bŷ where the parameters α and ρ vary within the bounds 0 < α ≤ 2, 1 − min(1, 1/α) ≤ ρ ≤ min(1, 1/α), t > 0, (A. 10) and the relation between the parameters in the A-and Z-forms is as follows, In the A-, B-, and M-forms β A = 1 corresponds to β M = 1 and β B = 1, while in the case of the Z-form the value β A = 1 corresponds to the value ρ = 1 if α < 1 and to the value ρ = 1 − 1/α if α > 1.

Appendix B. Some details of the numerical scheme
With the use of equations (5) and (12) we can write equation (4) on a discrete space-time grid as f j+1 Here we consider discretisation schemes for the three cases 0 < α < 1, 1 < α ≤ 2, and α = 1 separately.
Appendix B.1. 0 < α < 1 By substitution of equations (12) to (13b) into equation (B.1) and using the relation where the sign + is taken for x > b > a and − for x < a < b, we obtain Then the matrices A and B in equation (16) are with the definition (B.12) and changing the notation as above, we obtain Then the elements of the matrix A and B in equation (16) are and By substituting equations (15) and (12) into equation (9) we obtain Defining λ n = 1 2n + 1 , Ω L = 2L α,β K α ∆t ∆x , Ω R = 2R α,β K α ∆t ∆x , (B.17) changing notation as above, we obtain Then the elements of the matrices A and B in equation (16) are and (B.20)

Appendix C. Error estimation of the difference scheme
We here provide some details on the error estimate of our difference scheme. For the case 0 < α < 1 (equation (13)) for the left side derivative, and for the right side derivative. This efficient way to approximate the Caputo derivative of order α (0 < α < 1) is called L1 scheme [191,192,193] and its error estimate is O(∆x 2−α ) (see figure 2 top left panel) [192,194,195]. For the case 1 < α < 2 the suitable method to discretise the Caputo derivative is the L2 scheme [191,193,196] x for the left side, and for the right side. The truncation error for this scheme is O(∆x) (see figure 2 top right panel) [196,197]. For the special case α = 2 the central difference scheme has truncation error O(∆x 2 ) (see figure C1 top right panel). For comparison, in [194] an error estimate of order O(∆x 3−α ) is presented for 1 < α < 2. In [198] a computational algorithm for approximating the Caputo derivative was developed, and the convergence order is O(∆x 2 ) for 0 < α ≤ 2. Another difference method of order two was derived in [197] for 1 < α ≤ 2.
For the special case α = 1, β = 0, the truncation error is O(∆x 2 ) (see figure C1 top left panel). To evaluate the truncation error we used the relation is the exact solution and f j i is the approximated solution of function f (x, t). For e(t) 2 this is similar and we use The results of the error analysis for α,β (x, t), survival probability and the first-passage time PDF are shown in figure C1.

Appendix D. A short review of the Skorokhod theorem
The Skorokhod theorem establishes a general formula for the Laplace transform ℘(λ) of the first-passage time PDF in the semi-infinite domain for a broad class of homogeneous random processes with independent increments and thus has a pivotal role in the theory of first-passage processes [132,164]. For the process starting at x = 0 with a boundary at x = d, Here the auxiliary measure p + (λ, x) is defined via its Fourier transform as where the function f (x, t) is the PDF of the process, that is α,β (x, t) in our case. The boundary condition reads p + (λ, x) = 0 at x = 0. Below, for didactic purposes we first calculate ℘(t) for Brownian motion and then proceed to symmetric (0 < α ≤ 2 and β = 0), one-sided (0 < α < 1 and β = ±1), extremal two-sided (1 < α < 2 and β = ±1), and finally to the general case (0 < α < 2, α = 1, and −1 < β < 1).

Appendix D.1. First passage time PDF for Brownian motion
For Brownian motion the PDF reads

Substitution into equation (D.2) yields
Now, we use the relation to get to After an inverse Fourier transform according to equation (D.2) we arrive at For the first integral we have With the residue theorem, For the second integral, we write Again, with the residue theorem, Therefore, by substitution of equations (D.10) and (D.12) into equation (D.8) we obtain (x > 0) Using the boundary condition we get and thus with the help of equation (D.1), Finally, by inverse Laplace transform we get This is the famous Lévy-Smirnov law representing a well-known result of Brownian motion [165]. This derivation is, of course, overly complicated for the Gaussian case, but the same procedure can be applied to the general case of an asymmetric Lévy stable law, as we now show.
Appendix D.2. First passage time PDF for symmetric α-stable processes Due to the symmetry of the PDF the function ln q + (k, λ), equation (D.2), can be written as To find B(λ, k) at small λ the self-similar property of the α-stable process comes in useful, where α,0 (y, 1) = α,0 (−y, 1) =  Then the inverse Fourier transform of the above equation renders and, after applying the boundary condition, Recalling equation (D.1), Now, with the help of the Tauberian theorem [165] (Chapter XIII, section 5) we find that the small-λ asymptotic of the Laplace transform corresponds to the long-time asymptotic of the PDF ( [199], chapter 3) Therefore, the long-time asymptotic of the first-passage time PDF for symmetric αstable process has the form [125] Appendix D.3. First passage time PDF for one-sided α-stable processes, 0 < α < 1 and β = 1 Due to the monotonic growth of the process in this case there exists a simple relation between the cumulative probabilities of the first-passage time and the α-stable process itself, see, e.g., [168]. However, for a didactic purpose in this Appendix we obtain the result by the use of Skorokhod's method. Since for one-sided α-stable processes with 0 < α < 1 and β = 1 the PDF α,1 (x, t) vanishes for x < 0, we get After plugging this into equation (D.2), Therefore, ∂p + (λ, x)/∂x follows from result (D.32) by inverse Fourier transform, Recalling relation (E.5) we obtain d dx where E α (z) = ∞ n=0 z n /Γ(1 + αn) is the Mittag-Leffler function, see [138,167] and Appendix E. With the boundary condition p + (λ, x ≤ 0) = 0 we get Thus, for the Laplace transform of the first-passage time PDF ℘(λ) = 1 − p + (λ, d) we have The Laplace inversion is then immediately accomplished in terms of the Wright function (see equation (E.12)) for 0 < α < 1 [125], (D.39) By using relation (E.14) we make sure that ℘(t) is normalised, This can be also shown by taking the integral form (E.10) of the M -function. By changing the order of integration we get Appendix D.4. First passage time PDF for extremal two-sided α-stable process, 1 < α < 2 and β = −1 To apply the Skorokhod theorem we need to calculate the following integral To this end we use the Laplace transform of α-stable processes with 1 < α ≤ 2 and β B = −1, which is derived in [154] (page 169, equation (2.10.9)) in dimensionless B-form with K B α = 1 and t = 1, 43) In dimensional variables this equation reads and after plugging this expression into equation (D.2), (D.48) To calculate expression (D.48) we first find where we employ the Laplace transform (E.4) of the Mittag-Leffler function. By taking the indefinite integral over λ we obtain and then from equation (D.50), by inverse Fourier transform, With the boundary condition p + (λ, x = 0) = 0 we get Thus, for the Laplace transform ℘(λ) we obtain which is of a stretched exponential form. Recalling now the Laplace transform pair (E.13) for the M -function, we finally arrive at the first-passage time PDF for the extremal α-stable process with 1 < α ≤ 2 and β = −1, Let us show the normalisation of this function. By using the integral form (E.10) of the M -function and changing the order of integration we have (D.55) By change of variable u = [K α /d α | cos (απ/2)|σ] −1/α t −1/α , performing the inner integral and using the Hankel formula (E.3) for the Gamma function we obtain the necessary normalisation condition. Now, if we employ the series expansion (E.9) for the Mfunction, from equation (D.54) we arrive at a series which corresponds to that in equation (2.25) of [172]. Note that in our case the additional factor K α /| cos(απ/2)| appears due to a different starting form for the characteristic function of the α-stable process.
Appendix D.5. First passage time PDF for extremal two-sided α-stable processes, 1 < α < 2, β = 1 Similar to above, at first we obtain the Laplace transform of the α-stable PDF with 1 < α ≤ 2 and β = 1. We write 56) and then use the property α,β (x, t) = α,−β (−x, t) to get The second integral on the right side can be written as To take the first integral on the right hand side of equation (D.57) we use the characteristic function in the A-form. To take the second integral (D.58) we employ the Laplace transform of the PDF with 1 < α < 2 and β = −1 given by equation (D.45) with s = ik. Thus, relation (D.57) in the A-form has the shape With the help of equation (D.59) we write (again the index A is omitted in what follows) By substituting this expression into equation (D.2) we get The derivative with respect to λ reads where for the second term we employ the Laplace transform (E.4) of the Mittag-Leffler function. By taking the indefinite integral over λ we obtain Then, dp + (λ, x)/dx follows from equation (D.63) by inverse Fourier transform, and with the boundary condition p + (λ, x = 0) = 0 we get Thus, for the Laplace transform ℘(λ) = 1 − p + (λ, d) we obtain By applying the inverse Laplace transform with respect to λ and using the series representation (E.  Now, the generalised four-parametric Mittag-Leffler function has the series representation (page 129, equation (6.1.1) [200]) E α 1 ,β 1 ;α 2 ,β 2 (z) = ∞ k=0 z k Γ(α 1 k + β 1 )Γ(Γ(α 2 k + β 2 ) , z ∈ C, (D. 74) for α 1 , α 2 ∈ R (α 2 1 + α 2 2 = 0) and β 1 , β 2 ∈ C. It is an entire function, and if α 1 + α 2 > 0 it has the Mellin-Barnes integral form (page 132, equation (6.1.14) of [200]) where L = L −∞ is a contour running in a horizontal strip, from −∞ + iφ 1 to −∞ + iφ 2 , with −∞ < φ 1 < 0 < φ 2 < ∞. This contour separates the poles of the Gamma functions Γ(s) and Γ(1 − s). The function E α 1 ,β 1 ;α 2 ,β 2 (z) with α 1 + α 2 > 0 converges for all z = 0. For real values of the parameters α 1 , α 2 ∈ R and complex values of β 1 , β 2 ∈ C the four-parametric Mittag-Leffler function E α 1 ,β 1 ;α 2 ,β 2 can be represented in terms of the generalised Wright function and the Fox H-function. If α 1 > 0, α 2 < 0 and the contour of integration in expression (D.75) is chosen as L = L −∞ for α 1 + α 2 > 0, by identification with the corresponding Mellin-Barnes integral definiont of the H-function one can obtain the following representation of E α 1 ,β 1 ;α 2 ,β 2 in terms of the H-function We now employ theorem 4 from [201], which says that the Laplace transform with respect to x of an α-stable law in the Z-form of the characteristic function has the form In the next step we again change the order of integration, and thus (D.91) Going back to equation (D.2) by inverse Fourier transform we find (x > 0) .