First-passage functionals of Brownian motion in logarithmic potentials and heterogeneous diffusion

We study the statistics of random functionals $\mathcal{Z}=\int_{0}^{\mathcal{T}}[x(t)]^{\gamma-2}dt$, where $x(t)$ is the trajectory of a one-dimensional Brownian motion with diffusion constant $D$ under the effect of a logarithmic potential $V(x)=V_0\ln(x)$. The trajectory starts from a point $x_0$ inside an interval entirely contained in the positive real axis, and the motion is evolved up to the first-exit time $\mathcal{T}$ from the interval. We compute explicitly the PDF of $\mathcal{Z}$ for $\gamma=0$, and its Laplace transform for $\gamma\neq0$, which can be inverted for particular combinations of $\gamma$ and $V_0$. Then we consider the dynamics in $(0,\infty)$ up to the first-passage time to the origin, and obtain the exact distribution for $\gamma>0$ and $V_0>-D$. By using a mapping between Brownian motion in logarithmic potentials and heterogeneous diffusion, we extend this result to functionals measured over trajectories generated by $\dot{x}(t)=\sqrt{2D}[x(t)]^{\theta}\eta(t)$, where $\theta<1$ and $\eta(t)$ is a Gaussian white noise. We also emphasize how the different interpretations that can be given to the Langevin equation affect the results. Our findings are illustrated by numerical simulations, with good agreement between data and theory.


I. INTRODUCTION
Consider the stochastic trajectory of a one-dimensional particle described by the Langevin equation where µ(x) = −V ′ (x) represents a deterministic force derived from an external time-independent potential V (x), ).Here F (x) = 1/x and Z thus corresponds to the area under the graph of 1/x(t), where x(t) is the stochastic trajectory displayed in the inset.The trajectory starts from x0 = 1 and the motion is evolved up to the first-exit time from Ω.The diffusion constant is equal to one and V0 = 0.5.η(t) is a Gaussian white noise with zero mean and autocorrelation η(t)η(t ′ ) = δ(t − t ′ ) and D(x) is the spacedependent diffusion coefficient.Suppose that the motion generated by (1) starts from a point x 0 inside a given interval Ω, and the first passage outside Ω occurs after a random time T , which we call the first-passage time.Define where F (x) is, in principle, an arbitrary function that makes the integral convergent.Such a random variable is known as first-passage functional.Quantities of this kind have been extensively studied in the case of free Brownian motion, i.e., for µ(x) = 0 and D(x) = D, for the simple reason that many problems may be formulated in terms of first-passage Brownian functionals [1].
In this paper we wish to consider a subclass of firstpassage functionals, where the integral in ( 2) is evaluated for F (x) = x γ−2 , with γ ∈ R, over the trajectory of a Brownian particle with constant diffusion coefficient D in a logarithmic potential V (x) = V 0 ln(x/κ), where κ is a length scale that we can conveniently set to one.This choice is motivated by the fact that many interesting problems can be mapped to the study of functionals of this kind.For example, for γ = 2 one has F (x) = 1 and thus Z simply corresponds to the firstpassage time T , which is a stochastic quantity relevant for a plethora of applications [15,16].For γ = 3, Z is equivalent to the first-passage area A, i.e., the area swept by the trajectory x(t) in the xt plane till the first-passage time.This quantity has attracted a lot of interest and was studied for instance in the case of Brownian motion [3,5,17], Brownian motion with drift [3,4,6], Brownian motion with stochastic resetting and jump-diffusion processes [2,9,18], Orstein-Uhlenbeck process with and without resetting [8,12,14], Lévy processes [7], with applications in queueing theory and combinatorics [3], percolation [19], animal movements [20], snow melt [21] and DNA breathing dynamics [22], to cite a few examples.Other nontrivial and interesting cases are γ = 3  2 , which is related to the oscillation period in the underdamped one-dimensional Sinai model [23], and γ = 1  2 , which is associated with the lifetime of a comet in the solar system [1,24].Remarkably, in the case of free Brownian motion diffusing in Ω = (0, ∞), it is possible to obtain the distribution of Z for any γ > 0 [17].It is natural to try to extend this result to more general situations, for example by adding the presence of an external driving force.The specific case of a logarithmic potential is interesting for several reasons: first, it has been extensively studied in the literature [25][26][27][28][29][30] and recognized as a model naturally appearing in different contexts, such as stochastic thermodynamics [31][32][33][34], vortex dynamics [35,36], long-range interacting systems [37][38][39], ion condensation on a long polyelectrolyte [40], sleep-wake transitions [41], DNA denaturation [22] and diffusion of cold atoms in optical lattices [42][43][44][45][46][47]; in particular, in the latter two cases the first-passage area [22] and the area under an excursion [48,49], namely a trajectory that begins and ends at the origin without crossing it at intermediate times, are of particular interest.Second, there exists a discrete counterpart, known as the Gillis random walk [50,51], which can be solved exactly and has been considered in some recent work [52][53][54][55].This model is a critical case for the study of recurrence in stochastic processes [51,[56][57][58], with unique first-passage properties that are also recovered in the continuous system.Third, it has been shown that certain models of heterogeneous diffusion can be mapped to the dynamics of Brownian motion (with constant diffusion coefficient) in a logarithmic potential [59][60][61].Hence, obtaining the distribution of Z in the latter case allows us to derive also the solution of the problem in the case of a spatially-varying diffusion coefficient.We remark that heterogeneous diffusion has attracted a lot of interest in the statistical physics community [59,[62][63][64][65][66][67][68] and not only, as situations where D is nonconstant are ubiquitous: examples include contexts related to biology [69][70][71][72], finance [73], solute transport in heterogeneous media [74] and Richardson diffusion in turbulence [75].
The outline of the paper is the following: in the next section we use the method of [1] to write a backward evolution equation for the Laplace-transformed probability density function of Z, when evaluated along a trajectory generated by Eq. (1).Then in section III we summarize the main results for the particular case of a logarithmic potential and a constant diffusion coefficient.As a corollary, we also obtain the distribution of Z when x(t) is generated by ẋ(t) = √ 2Dx θ η(t), with θ < 1, which is a model for heterogeneous diffusion.In sections IV, V and VI we derive the results by providing detailed calculations.Finally, in section VII we draw our conclusions.

II. BACKWARD EQUATION FOR THE PROBABILITY DENSITY FUNCTION
Let us call p(z, x 0 ) the probability density function (PDF) of Z, knowing that the trajectory started from x 0 ∈ Ω.The idea is to derive a backward evolution equation for the Laplace transform of the PDF, which corresponds to the expected value of e −wZ , where w is the Laplace variable: Here the expected value is taken over all realizations that start from x 0 and leave Ω for the first time at T .To do this, one can rewrite Eq. ( 2) as [1] and note that the second integral at the right-hand side (rhs) corresponds to the definition of Z, but for a trajectory that starts from a random position x(dt) = x 0 + dx(0).Hence by using p(w, x 0 ) = E e −wZ , we have where the average at the rhs is taken over all possible x(dt), viz., over all possible dx(0).According to Eq. ( 1), for any t the displacement dx(t) = x(t+dt)−x(t) is given by where dW (t) is the increment of a Wiener process of variance dt and, more importantly, x * is a point between x(t) and x(t + dt).The choice of the point depends on the interpretation given to the Langevin equation (1), and different choices lead to different solutions [76][77][78][79][80].In other words, if we set x * = αx(t + dt) + (1 − α)x(t), with 0 ≤ α ≤ 1, the value of α determines the "rule" to integrate (1), and the choice is often motivated by physical reasons.The interpretations considered most significant in the physics literature are those of Itô (α = 0), 2 ) and Hänggi-Klimontovich (α = 1) [81][82][83][84].In our case, it is useful to make the nonanticipating choice α = 0 (Itô).Nevertheless, we are not bound to consider exclusively the Itô interpretation, as any other interpretation can be recovered by inserting an additional drift term dependent on α.In other words, (7) is equivalent to where Now, from Itô formula [77] we can write d p(w, x 0 ) = p(w, x 0 + dx) − p(w, x 0 ) as thus by inserting this in Eq. ( 6), taking the average over dW and discarding terms that are o(dt), we obtain which is the backward evolution equation for p(w, x 0 ), to be accompanied by the appropriate boundary conditions and the normalization condition p(0, x 0 ) = 1.If D(x) is always bigger than zero in Ω, we can define where the lower bound of integration can be any point of Ω, and then (11) can be written as where We can therefore set to obtain a simpler equation for ψ(w, x 0 ): Note that the normalization condition p(0, x 0 ) = 1 imposes ψ(0, x 0 ) = N (x 0 ).

III. SUMMARY OF THE MAIN RESULTS
For V (x) = V 0 ln x, with −∞ < V 0 < ∞, and D(x) = D, Eq. ( 16) simplifies to where we have introduced the parameter We start by considering the dynamics in an interval Ω = (a, b), with 0 < a < b.Then we generalize to intervals of the kind Ω = (0, a) and Ω = (b, ∞).Finally, we will consider the problem in Ω = (0, ∞).In the last scenario, we will also provide the solution when the dynamics is generated by a Langevin equation of the kind with θ < 1, and underline how it depends on different interpretations.To simplify the following formulas, it is convenient to define for γ = 0 the exponent and use the notation q to indicate the scaled variable where w will be the Laplace variable.

A. Finite intervals left-bounded by a positive number
Consider Ω = (a, b), with 0 < a < x 0 < b.For γ = 0, the Laplace transform p(w, x 0 ) is given by where H ν (x, ŷ) is defined as Here I ν (z) and K ν (z) are the modified Bessel functions of the first and second kind, respectively [85].Similarly, for γ = 0 we define and have which can be inverted, yielding the PDF Furthermore, if we consider the set of trajectories that leave Ω = (a, b) from b, which has probability then the distribution of Z measured only on those trajectories has the normalized density Analogously, the probability E a (b) of leaving from a can be obtained from Eq. ( 27), and the corresponding normalized density from Eq. ( 28), by exchanging a and b.

B. Finite intervals left-bounded by the origin, or infinite intervals left-bounded by a positive number
Now take Ω = (0, r) or Ω = (r, ∞), with r > 0.
There is a correspondence between the solutions in the two cases.More precisely: 1.When β > 0 and Ω = (0, r), the functional Z is well-defined only for γ > 0. In this case, the Laplace transform p(w, x 0 ) is given by Nevertheless, if we examine only the set of trajectories that leave from r, which has probability E R = (x 0 /r) β , then Z is well-defined for any γ, and the corresponding normalized conditional PDF is whereas for γ = 0 the conditional PDF is given explicitly by Similarly, if we take now Ω = (r, ∞) and exchange β → −β and γ → −γ, the set of trajectories that leave the interval in a finite time has probability E = (x 0 /r) β and the corresponding normalized conditional PDF is given again by Eqs.(30) and (31).
2. When β ≤ 0 and Ω = (0, r), a trajectory leaves the interval from r with probability one, thus Z is welldefined for any γ.The Laplace transform p(w, x 0 ) is given by and for γ = 0 the PDF is given explicitly by Similarly, if we take now Ω = (r, ∞) and exchange β → −β and γ → −γ, the set of trajectories that leave the interval in a finite time has probability one and the corresponding PDF is given again by Eqs.(32) and (33).

C. Positive real axis
For the positive real axis Ω = (0, ∞), the functional Z is well-defined only when both β and γ are positive.However, in this case the PDF can be computed explicitly.By defining the PDF can be written as where Γ(ν) is the Euler Gamma function.The result for free Brownian motion [17] is recovered by setting V 0 = 0, i.e., by putting β = 1, which yields the exponent ν = 1/γ.As a corollary, consider now Z evaluated on trajectories generated by the Langevin equation with θ < 1, that we may interpret with any 0 ≤ α ≤ 1. where The same applies to 1 2 < α ≤ 1, if we add the condition θ < 1 2α .By way of illustration, the first-passage time density is recovered from Eq. ( 37) by setting γ = 2: with ν α = (1 − 2αθ)/(2 − 2θ), which agrees perfectly with recent results [86].
In the following, we go into the details of the derivation and present plots in which we compare our findings with numerical simulations.

IV. FINITE INTERVALS LEFT-BOUNDED BY A POSITIVE NUMBER
In this section we deal with intervals of the kind Ω = (a, b).This case can be treated for any value of γ, but we must distinguish γ = 0 and γ = 0.
Let us begin with γ = 0. Eq. ( 17) can be brought back to the modified Bessel equation: To see this, we make the ansatz ψ(w, x 0 ) = x ρ 0 ϕ(λx σ 0 ), where λ depends on w and arrive at the following equation for ϕ(z): where z = λx σ 0 .Then, by choosing we obtain the modified Bessel equation (40), with ν = β/γ, which admits the general solution where c 1 and c 2 are coefficients that depend on a, b and w, and we recall see Eq. ( 21).The function ψ(w, x 0 ) is thus of the form ψ(w, x 0 ) = √ x 0 ϕ(2x 0 ) and by recalling Eq. ( 15), we must have p(w, x 0 ) = x To determine the correct boundary conditions, we just note that when the starting point of the trajectory is close to one of the boundaries, the first-passage time tends to zero, and so does the integral in Eq. ( 2).Hence p(w, x 0 ) = E(e −wZ ) must be equal to one for x 0 equal to a or b: Then c 1 and c 2 can be determined from the simple linear system where and The solution can be finally written as One can verify that p(w, x 0 ) satisfies the normalization condition p(0, x 0 ) = 1.Now we consider the case γ = 0, which corresponds to Equation ( 16) reads and now we seek solutions of the form ψ(w, x 0 ) = ϕ(λ ln x 0 ).The corresponding equation for ϕ(z) is just a second order linear ordinary differential equation with constant coefficients.The characteristic roots are and the solution can be thus written as where By using ψ(w, x 0 ) = ϕ(λ ln x 0 ), we see that the value of λ is arbitrary, so we can set it to one.Finally, recalling p(w, x 0 ) = x where c 1 and c 2 have to be determined once again in such a way that the boundary conditions p(w, a) = 1 and p(w, b) = 1 are satisfied.Similarly to the previous case, we have to solve an equation of the kind but this time with Once the coefficients have been determined, we find that the solution can be written as which has the same structure as Eq. ( 22), with H ν (x, ŷ) replaced by A. Conditioning on leaving the interval from a given boundary The structure of ( 50) and ( 61) allows us to easily solve the problem with the additional condition that x(t) leaves the interval from a chosen boundary.This request is relevant, for instance, in extreme value theory [53,[87][88][89][90][91][92], where the probability of leaving the interval from a given boundary is related to the statistics of the maximum or the minimum of the process.
Let us take s > 0 and q > 0, and assume min{s, q} < x 0 < max{s, q}.Define p s (z, q, x 0 ) as the PDF of the functional Z under the condition that x(T ) = s, which corresponds to requiring the process to leave the interval from s.Note that, according to this definition, where E s (q) is the splitting probability, namely, the probability of leaving the interval from s.The Laplace transform must thus satisfy the following boundary conditions: a) p s (w, q, s) = 1: just as the case considered previously, if the dynamics starts close to the boundary s, the first-passage time and thus also Z tend to zero, hence E(e −wZ ) → 1; b) p s (w, q, q) = 0: if the dynamics starts instead close to the other boundary q, the probability of leaving from s tends to zero, and so does the integral in Eq. ( 63), from which it follows that also p s (z, q, x 0 ) and its Laplace transform must vanish.
From ( 50) and ( 61), we see that the solutions we found previously are written as the sum of two terms.Taking for example (50), it is easy to verify that and the same is true for Eq. ( 61).Therefore The splitting probability E s (q) can be computed by using the results of Appendix A. We find and consequently, the complementary probability of leaving from the other boundary is E q (s) = 1 − E s (q).Before continuing, let us illustrate the results of this section more explicitly.When γ = 0, the Laplace transform of Eq. ( 66) can be inverted for some particular values of the exponent ν.For example, when ν = ± 1 2 , one can write the modified Bessel functions in terms of elementary functions [85] which makes the inversion particularly easy.Indeed, the function H ν (x, y) becomes thus where we have defined The poles of p s (w, q, x 0 ) are w n = −(nπ/L) 2 , so the inversion yields From p(z, x 0 ) = p s (z, q, x 0 ) + p q (z, s, x 0 ) it is then possible to obtain the full distribution.Recall that ν = ± 1 2 implies β = ± γ 2 , hence the validity of this result is limited to those cases.When γ = 0 instead, p s (z, q, x 0 ) can be inverted for any value of β.The poles are now w n = −D[nπ/ ln(q/s)] 2 − Dβ 2 /4, so we get p s (z, q, x 0 ) = 2πD ln 2 (q/s) and from p(z, x 0 ) = p s (z, q, x 0 ) + p q (z, s, x 0 ), with s = a and q = b, we obtain Eq. ( 26).
In Fig. 2 we show examples of p(z, x 0 ) for γ = 3, γ = −3 and γ = 0. Note that in the first two cases, the condition ν = ± 1 2 requires us to choose β = ± 3 2 .The datasets are obtained by measuring Z over trajectories in Ω = ( 1 2 , 5 2 ) starting from x 0 = 1 and evolved up to the first-exit time.The details on the numerical simulations are given in Appendix C. We find that the comparison between the theoretical densities and the numerical results is good for any combination of the parameters.Interestingly, while the cases with ν = 1 2 display a unimodal distribution, for ν = − 1 2 the distribution can become bimodal, as shown in panel (a).

V. FINITE INTERVALS LEFT-BOUNDED BY THE ORIGIN AND INFINITE INTERVALS LEFT-BOUNDED BY A POSITIVE NUMBER
We now want to generalize the treatment to intervals of the type Ω = (0, b) or Ω = (a, ∞).These cases may be interpreted as the limit for a → 0 or b → ∞ of the results of Sec.IV, and both introduce difficulties not present before.For example, studying the problem in (0, b) allows the trajectory to hit the origin, so there may be values of γ for which the integral defining Z does not converge, see Eq. ( 2).For (a, ∞), since the motion occurs in an infinite domain under the action of an external potential, there may be realizations for which the first-passage time T is not finite.For these reasons, it will be also necessary to rediscuss the boundary conditions that the solutions must satisfy.
A. Finite intervals of the kind Ω = (0, b) When we consider diffusion in a logarithmic potential V (x) = V 0 ln(x), the first-passage problem to the origin must be treated with special attention.Indeed, the nature of this point depends on the relative magnitude of the potential V 0 with respect to the diffusion constant D, which in our discussion is measured by the parameter β.A detailed analysis following Feller's classification scheme can be found in Ref. [28], according to which the origin is an exit boundary for β ≥ 2, a regular boundary for 0 < β < 2 and an entrance boundary for β ≤ 0. Exit and regular boundaries are both accessible; entrance boundaries are inaccessible [93], meaning that they can not be reached in finite time from the interior of the state space (in our case, from any point inside Ω).We can see this by evaluating the splitting probability E a (b) in the limit a → 0, see Eq. ( 67): Hence when β ≤ 0 the probability of hitting the origin vanishes, and a trajectory will leave Ω from b with probability one.From a physical point of view, this can be motivated by noting that for β ≤ 0 there is a strong repulsive potential, with V 0 ≥ D, pushing the particle away from the origin.As a consequence, the boundary condition p(w, 0) = 1 is no longer correct: even if the process starts very close to the origin, it is immediately pushed inside Ω and the motion goes on up to the first-passage to b.Thus, contrarily to what we had previously, T does not vanish in this case.The fact that the trajectory may or may not leave the origin affects the value of γ we can choose so that the integral in Eq. ( 2) will be convergent.Therefore in the following we differentiate the treatment depending on the sign of γ.
1. Functionals with γ > 0 When γ > 0, the limit a → 0 corresponds to the limit â → 0. From the results in Appendix A, we have thus if we evaluate p b (z, a, x 0 ) as a → 0, see Eq. ( 66), we obtain which is the contribution of trajectories hitting b, whereas the contribution of trajectories hitting the origin behaves for small a as whose limit depends on the sign of β.We can use [85] to see that for β > 0, since we have |ν| = ν = β/γ, both terms give a nonvanishing contribution in the limit, and the solution (50) converges thus to the function which is normalized, since p(0, x 0 ) = 1, and furthermore satisfies the boundary conditions p(w, 0) = 1 and p(w, b) = 1.For β ≤ 0 instead, the rhs of (80) vanishes, which is consistent with the fact that the origin is an entrance boundary.Thus in the limit a → 0 only the contribution of p b (w, a, x 0 ) survives and the solution converges to We remark that, even if this expression originates only from p b (w, a, x 0 ), it actually represent the full distribution, as can be seen from the fact that it satisfies the normalization condition p(0, x 0 ) = 1.At the boundaries, we have p(w, b) = 1, as expected, whereas for x 0 = 0 we obtain which is always smaller than one for b > 0, as one can verify by using the series expansion of the modified Bessel function.Hence, consistently with the fact that the origin is an entrance boundary when β ≤ 0, the functional Z is strictly positive even when measured on trajectories that start very close to x = 0.

Functionals with γ < 0
For γ < 0, the limit a → 0 is equivalent to the limit â → ∞, which yields (see Appendix A) then as a → 0 (â → ∞) whereas which vanishes when â → ∞, due to the exponential divergence of I ν (2â).Hence, differently from the case γ > 0, the contribution of the walks that hit the origin vanishes independently of the value of β, and the only relevant term is We note that this function satisfies p b (w, b) = 1, but vanishes for any ν in the limit x 0 → 0, due to the behavior of K ν (z) at infinity [85] K Moreover, by computing p b (0, x 0 ) we obtain which corresponds to the splitting probability E R = 1 − E L of leaving Ω from b, see Eq. ( 77) for the expression of E L .We can interpret this result as follows: when β ≤ 0, as we mentioned earlier, the origin is an entrance boundary, hence it can not be reached from the interior of Ω and all the trajectories starting from x 0 > 0 leave the interval from b.Therefore, the functional Z can be measured over each trajectory and Eq. ( 88) describes the full distribution, that is, p(w, x 0 ) = p b (w, x 0 ).On the other hand, when β > 0, a trajectory can leave the interval from any of the two boundaries, but those that leave Ω from the origin yield a diverging Z for γ < 0. In other words, Z is not well-defined if we allow the particle to hit the origin.As we can deduce from the fact that it is normalized to E R , Eq. ( 88) in this case describes the distribution of Z measured on the set of trajectories that leave the interval from b, namely, it is the conditional distribution.We can then define p(w, x 0 ) = p b (w, x 0 )/E R , so that p(w, x 0 ) always denotes the normalized PDF.In summary, we have The fact that p(w, x 0 ) tends to zero as x 0 approaches the origin means that Z is diverging, so contributions from paths passing near x = 0 can be expected to cause a heavy-tailed decay of the distribution.Indeed, by expanding in powers of w, we find with logarithmic corrections appearing for ν = 0 and ν = ±1.By using Tauberian arguments [94], we can conclude that the PDF is characterized by a power-law decay as p(z, x 0 ) ∼ z −1−|ν| .This can be shown even more explicitly if we consider ν = ± 1 2 , for which the inversion can be carried out easily.Indeed, in both cases we have and the inversion yields which indeed decays as p(z, x 0 ) ∼ z −3/2 .

The case γ = 0
We now analyze the particular case γ = 0.For any β, the limit a → 0 of p b (w, a, x 0 ) yields with whereas p a (w, b, x 0 ) vanishes in the limit.So we are in the same situation as the case γ < 0: for β > 0, there is a positive probability E L = 1 − E R that a trajectory hits the origin, yielding a diverging Z, hence the distribution must be measured only on the walks that leave Ω from b.
For β ≤ 0 instead, a trajectory leaves from b with probability one, hence p b (w, x 0 ) corresponds to the full distribution.In both cases, we can set p(w, x 0 ) = p b (w, x 0 )/E R and write which satisfies p(0, x 0 ) = 1 and p(w, b) = 1 and vanishes for x 0 → 0. The inverse transform of Eq. ( 97) is We see that for z → 0 the PDF goes to zero as Hence, for β = 0 there is an exponential cut-off ensuring the convergence of all moments, while for β = 0 we observe a pure power-law decay p(z, x 0 ) ∼ z −3/2 as z → ∞.
The PDF is displayed in Fig. 3 and compared to the results obtained from numerical simulations, showing good agreement.The chosen values of β cover all the cases: for β = 0.5 the theoretical result of (98) corresponds to the full distribution, while for β = 0 and β = −0.5 it is the PDF of Z measured on walks that leave the interval from b, normalized dividing by E R .Note that p(z, x 0 ) only depends on the sign of β, hence the PDFs of the data with β = 0.5 and β = −0.5 are described by the same theoretical curve, as we observe in the figure .B. Infinite intervals of the kind Ω = (a, ∞) The case of infinite intervals marks a difference with the previous one in that the particle is not guaranteed for trajectories in Ω = (0, 1).The starting point is x0 = 0.7 and the diffusion coefficient is set to one.The solid black curves are the theoretical predictions given by Eq. ( 98), the symbols represent the numerical distributions obtained from simulations.We have considered β = −0.5 (magenta squares), β = 0 (blue circles) and β = 0.5 (turquoise asterisks), with good agreement in all cases.Note that data corresponding to the same |β| overlap, confirming the symmetry in β of the PDF.In the inset, a plot in logarithmic scale of the case β = 0, which shows that the simulations also capture the heavy tail of the distribution.The number of simulations is 10 8 for β = 0.5 and 10 7 in the other cases.The trajectories are evolved up to the first-passage time with small time step ∆t = 10 −5 .
to leave Ω in a finite time.Indeed, if we consider E a (b) given by (67), and take the limit b → ∞, we get meaning that for negative values of β, there is a nonzero probability 1 − E to observe an infinite first-passage time.
Note that if we take also the limit a → 0, then E converges to 0 for β < 0, i.e., the set of trajectories that do not leave Ω = (0, ∞) has probability one, as is known [35].To avoid having to deal with generalized functionals of the form in the following we restrict ourselves to the case where actually T < ∞.Note that it would not be appropriate to speak about first-passage functionals if the first-passage time is not finite.
1. Functionals with γ > 0 When γ has positive sign, the limit b → ∞ corresponds to the limit b → ∞.Since then as b → ∞ the function p a (w, b, x 0 ), see Eq. ( 66), converges to (104) which is the same result obtained for the problem in (0, b) in the case γ < 0, see Eq. ( 88).The limiting function satisfies p a (w, a) = 1 and vanishes for x 0 → ∞, i.e., p a (w, ∞) = 0.The latter condition may be explained by the fact that if the motion starts very far from a one observes a very large first-passage time, and thus we should expect larger and larger values of Z, with p(w, x 0 ) = e −wZ consequently vanishing.Regarding the normalization, we get p a (0, x 0 ) = E, therefore we conclude that Eq. ( 104) is the PDF describing the full distribution of Z for β ≥ 0, while for β < 0 it describes the distribution of Z measured only on the trajectories that actually leave Ω = (a, ∞) at some finite T .Hence we define again the normalized distribution as Note that this is equivalent to (91), hence the same considerations follow.
2. Functionals with γ < 0 Here the limit b → ∞ is equivalent to b → 0. By using we see that as b → ∞ ( b → 0), the Laplace transform p a (w, b, x 0 ) goes to which satisfies p a (w, a) = 1, while p a (0, x 0 ) = E, see the expression for E in Eq. ( 101).Therefore, when E = 1, i.e., for β ≥ 0, we have p(w, x 0 ) = p a (w, x 0 ), whereas in the opposite case we set p(w, x 0 ) = p a (w, x 0 )/E, that is Remarkably, for x 0 → ∞, this function converges to which implies the convergence of all the moments even in the large-x 0 limit.

3.
The case γ = 0 The limit b → ∞ of (66) yields which satisfies p a (w, a) = 1 and vanishes in the limit x 0 → ∞, for the same reason of the case γ > 0. We have once again p a (0, x 0 ) = E, hence the normalized distribution is p(w, x 0 ) = p a (w, x 0 )/E, which can be inverted, yielding This result is equivalent to what we obtained in the interval (0, b), therefore the considerations made in that case still hold.

VI. POSITIVE REAL AXIS
It is straightforward now to obtain the solution of the problem in the positive real axis Ω = (0, ∞).It should be clear that Z is well-defined only for β > 0 and γ > 0. Indeed, as we discussed previously, the condition β > 0 is necessary for the first-passage time to be finite, while γ > 0 ensures that Z does not diverge when a trajectory hits the origin.Hence we restrict to β > 0 and γ > 0, viz., ν > 0.
The solution can be obtained by simply considering the limit b → ∞ of (29), which yields This is basically equivalent to the result obtained for free Brownian motion, see [17], which can be recovered by setting V 0 = 0, viz., β = 1, yielding the exponent ν = 1/γ.By introducing a logarithmic potential, one thus obtains a generalized exponent ν = β/γ.The Laplace transform can be inverted exactly by using [95] valid for R(y) > 0, yielding We note that the PDF can be written in scaling form as p(z, x 0 ) = P (ζ)/Z D , where T 0 [x(t)] −1 dt for trajectories in the positive real axis starting from x0 = 1 and evolved up to the first-passage time to the origin.The diffusion coefficient is set to one.The solid black curves are the theoretical predictions given by Eq. ( 115), the symbols represent the numerical distributions obtained from simulations.For both values of γ we considered different values of β: (a) in the case γ = 3 we choose β = 1.5 (red squares), β = 4 (blue circles) and β = 7 (green triangles), (b) in the case γ = 1 we take β = 0.8 (red squares), β = 1.5 (blue circles) and β = 4 (green triangles).The number of simulations depends on β: for β = 0.8 and β = 1.5 it is 3 • 10 5 , for β = 4 and β = 7 it is 3 • 10 6 .All the trajectories are evolved with small time step ∆t = 10 −5 .
For ζ → 0 the function P (ζ) vanishes displaying an essential singularity, while for ζ → ∞ we get a power-law decay P(ζ) ∼ ζ −1−ν .Therefore, the m-th moment of the distribution is finite only for m < ν, in which case is equal to Note that ν > m means β > mγ, hence for fixed γ we can tune β so that all moments up to the m-th are finite.On the other hand, for fixed β the m-th moment is finite only if γ < β/m, therefore we can for instance have a finite T but a diverging A (first-passage area).
In Fig. 4 we present the distribution of the scaled variable Z/Z D given by Eq. ( 115), and compare it with numerical data.We show the cases γ = 3 and γ = 1, each with three different values of β, chosen so that for both functionals a case with 0 < ν < 1 (infinite mean and variance), one with 1 < ν < 2 (finite mean, infinite variance) and one with ν > 2 (finite mean and variance) is displayed.The agreement between data and theory is evident in all cases.

A. Heterogeneous diffusion
We now extend the results of this section to the case where Z is measured over stochastic trajectories generated by with θ < 1.As discussed in Sec.II, different interpretations can be assigned to this Langevin equation, and we will see how the results are affected by the interpretation.
One possible approach to this problem would be to write down Eq. ( 16) for D(x) = Dx 2θ and µ α (x) = αD ′ (x), and then solve the resulting equation, namely, accompanied by the appropriate boundary conditions.The solution may be then used to obtain the Laplace transform of the PDF.However, as it has already been pointed out [59][60][61], there exists a close relation between Brownian motion in logarithmic potentials and heterogeneous diffusion which we may exploit to obtain the solution in a much more straightforward way.Let us call x(t) a trajectory generated by with x(0) = x 0 > 0 and evolving in Ω = (0, ∞) till the first-passage time to the origin T .Let θ < 1 and define the following transformation on the trajectory: By applying Itô formula on y(t), we see that the transformed trajectory evolves according to which is interpreted in Itô scheme.We recall β = 1 + V 0 /D.This is the Langevin equation of a system with a  117) and evolved up to the first-passage time to the origin.The starting point is y0 = 1 and the space-dependent diffusion coefficient is D(x) = x 2/3 .The solid black curves are the theoretical predictions given by Eq. ( 128), the symbols represent the numerical distributions obtained from simulations.For both values of γ we considered three different interpretations, corresponding to Itô (α = 0, green triangles), Stratonovich (α = 1 2 , blue circles) and Hänggi-Klimontovich (α = 1, red squares).The insets, which are presented in log-log scale in (a) and semilog scale in (b), display the heavy-tails of the distributions.All datasets are obtained by evaluating Z over 10 6 trajectories, with small time step ∆t = 10 −4 .The integration scheme used to integrate Eq. ( 117) is the Itô scheme, and the desired interpretation is obtained by adding the appropriate drift µα(x) = αD ′ (x).
space-dependent diffusion coefficient D(x) = Dx 2θ and a drift term that may be written as Thus, by setting the coefficient in front of D ′ (x) equal to α, we obtain exactly the Itô form of ẏ = 2D(y)η(t) in the α-interpretation.It is immediate to see that α = 1 2 , which corresponds to Stratonovich interpretation, is recovered by setting β = 1, i.e., by identifying x(t) as free Brownian motion, with a free choice of θ < 1.This mapping between Brownian motion and heterogeneous diffusion with Stratonovich interpretation is well-known, see for instance Ref. [62].All other interpretations can be obtained by observing that the parameters α, β and θ are related by This also means that for fixed α we can tune θ by changing the value of β in the original model.One must recall, however, that since θ < 1, one is limited to take β < 2α when α > 1 2 and β > 2α when α < 1 2 .It is clear that the first-passage time to the origin of the original trajectory x(t), starting from x 0 , is the same as the transformed trajectory y(t), with the initial condition y 0 = [(1 − θ)x 0 ] 1/(1−θ) .Hence by using Eq.(120) we can write with If the functional at the rhs has a proper distribution, with PDF p(z, x 0 ), then the lhs has a proper distribution too, with PDF We recall that this is true for x(t) if both β and γ ′ are positive.The first condition is always met when α < 1 2 , because θ < 1 implies β > 2α > 0; when α > 1 2 instead, this must be added to the previous condition β < 2α, obtaining 0 < β < 2α, which implies that we are limited to θ < 1 2α .The condition hence the lower bound for γ is the exponent appearing in the expression of the diffusion coefficient.By using Eq.(114), we obtain where Hence, the interpretation given to the Langevin equation strongly affects the distribution by changing the powerlaw decay exponent of the PDF.
The results are displayed in Fig. 5 for the cases γ = 3 and γ = 3  2 , with a diffusion coefficient D(x) = Dx 2/3 .For each γ, we choose three possible interpretations: α = 0 (Itô), α = 1 2 (Stratonovich) and α = 1 (Hänggi-Klimontovich), corresponding to the exponents Note that for θ > 0 we have ν I > ν S > ν HK , whereas the opposite happens for θ < 0. For the chosen values of θ and γ, we obtain in every case an heavy-tailed distribution: for the first-passage area we have ν I = 3 7 , ν S = 2 7 , and ν HK = 1 7 , while for γ = 1.5 we get ν I = 6 5 , ν S = 4 5 and ν HK = 2 5 .The agreement between theory and numerical results is generally good.The data can replicate all the features of the PDF, including the tails, see the insets in both panels.We remark that the numerical results have been obtained by measuring Z over trajectories generated by Eq. ( 117), hence they are independent of the method we discussed here.

VII. CONCLUSIONS AND DISCUSSION
In this paper we have studied the statistical properties of random variables of the kind Z = T 0 [x(t)] γ−2 dt, where x(t) is a one-dimensional trajectory of Brownian motion with diffusion constant D evolving under the effect of a logarithmic potential V (x) = V 0 ln(x) that can be either attractive or repulsive.The trajectory starts from x 0 inside a given interval Ω and leaves it for the first time at some random instant T .We initially considered the problem for Ω = (a, b) entirely contained in the positive real axis, which can be treated for any γ and any value of V 0 .We then generalized to intervals of the kind Ω = (0, b) or Ω = (a, ∞).Both these generalizations introduce some limitations: in the former case, for γ < 0 the functional Z is defined in terms of a divergent integral when measured on trajectories hitting the origin.In the latter case, the presence of a repulsive potential may prevent the particle to leave Ω, which implies an infinite first-passage time.Interestingly, we have underlined that there is a correspondence between the solutions of the two cases, if we always restrict the study of Z on trajectories for which it is well-defined.Finally, we have computed exactly the density of Z when it is constructed on trajectories in Ω = (0, ∞), with γ > 0 and V 0 > −D.By using a close relation between Brownian motion in logarithmic potentials and heterogeneous diffusion, we have also obtained the distribution of Z measured on trajectories x(t) generated by ẋ = √ 2Dx θ η(t), with θ < 1.
This work extends some previously known results regarding first-passage functionals of Brownian motion [17].By introducing a potential, we were able to study how it affects the statistical properties of Z for a fixed value of γ.We emphasize that the logarithmic potential has unique properties, which stem from the fact that it grows as a slowly varying function for x → ∞, yielding a force that is proportional to 1/x.As already noted in the literature, this causes both the drift term and the diffusion term in the Fokker-Planck equation to scale as 1/x 2 [25,29].Therefore, the two effects (diffusion and drift) are comparable as long as the dynamics takes place away from the origin, and the system can be treated effectively as a perturbation of Brownian motion.Not surprisingly, the results we obtained in Sec.VI have the same functional form as those obtained for free Brownian motion [17].Nevertheless, the system is far from being trivial, as its behavior can be drastically modified by adjusting the parameters that govern the intensity of the drift and diffusion terms, namely the strength of the potential and the diffusion coefficient.This has consequences regarding for example the emergence of nonnormalizable steady states [25,26,29] or the recurrence properties, of which this system is a critical case study, as evident from the analysis of related discrete models [51,[56][57][58].For the problem considered in this paper, for instance, we found that in the case Ω = (0, ∞) the PDF has a power-law decay as p(z, x 0 ) ∼ z −1−ν , with ν = (D + V 0 )/(Dγ), which means that the distribution has infinite variance for V 0 < D(2γ − 1) and also infinite mean for V 0 < D(γ − 1).
Another interesting feature of Brownian motion in a logarithmic potential is that it is associated with heterogeneous diffusion, which is studied in many contexts.Our results can be easily extended to the case where the dynamics is generated by ẋ = 2D(x)η, with D(x) = Dx 2θ and θ < 1, as we have done for Ω = (0, ∞).In this context, a key role is played by the interpretation given to the Langevin equation, and we have seen how the value of the interpretation parameter α contributes, along with the exponent θ, to determine the statistics of Z for a given value of γ.
Finally, let us remark that the densities of Z over trajectories in Ω = (0, ∞) for logarithmic potentials and heterogeneous diffusion, given in Eq. (114) and Eq. ( 128) respectively, have the same structure, viz, in both cases we obtain the PDF of an Inverse-Gamma random variable [96].This fact is strictly connected to the property of selfsimilarity, which is shared by both models, as shown in Refs.[60,61].There the author proves that for any selfsimilar diffusion process the first-passage time to the origin has Inverse-Gamma statistics.We have found that the same statistics describes also Z = T 0 [x(t)] γ−2 dt, if it exists.Although this observation could be deduced from scaling arguments, at least with regard to the asymptotic behavior for large z [3,17], it was not trivial to determine how the entire distribution changes with γ.
As a future perspective, one can ask how the different functionals studied in this article are correlated.The correlation can be measured, for example, by computing the joint probability distribution between two observables Z 1 and Z 2 evaluated for two different γ.In particular, the case where one of them corresponds to the first-passage time may be particularly relevant, so that information on the correlation between spatial and temporal variables can be obtained directly.This type of joint distribution is indeed useful for comprehensively quantifying the properties of stochastic search processes, as recently observed and studied in [97].For the sake of completeness, we also consider p s (w, q, x 0 ) = x 0 s β/2 H ν (x 0 , q) H ν (ŝ, q) , (B9) and evaluate the coefficient of the linear term in the expansion in powers of w.In general, we can write p s (w, q, x 0 ) = E s (q) 1 + w Dγ 2 f 1 (x 0 , q; ν) f 0 (x 0 , q; ν) − f 1 (s, q; ν) f 0 (s, q; ν) + o(w) , (B10) where the splitting probability is given by E s (q) = x 0 s β/2 f 0 (x 0 , q; ν) f 0 (s, q; ν) .(B11) Then the conditional first moment is just Z s (q) = 1 Dγ 2 f 1 (s, q; ν) f 0 (s, q; ν) − f 1 (x 0 , q; ν) f 0 (x 0 , q; ν) , (B12) and one can verify that the previous general expression for Z can be obtained from Z = E s (q) Z s (q) +E q (s) Z q (s) .By skipping details, for ν = 0, ±1 we find Z s (q) = 1 Dγ 2 s γ 1 − ν 1 − (s/q) γ(ν−1) 1 − (s/q) γν − 1 − (x 0 /q) γ(ν−1) 1 − (x 0 /q) γν x 0 s γ + q γ 1 + ν 1 − (s/q) γ(ν+1) 1 − (s/q) γν − 1 − (x 0 /q) γ(ν+1) 1 − (x 0 /q) γν , (B13) the cases ν = ±1 are both covered by Z s (q) = 1 Dγ 2 s γ − x γ 0 2 + γq γ s γ ln(q/s) s γ − q γ − x γ 0 ln(q/x 0 ) and ν = 0 yields Z s (q) = 1 Dγ 2 s γ − x γ 0 + 2(q γ − s γ ) γ ln(s/q) − 2(q γ − x γ 0 ) γ ln(x 0 /q) .(B15) order 2 Runge-Kutta method [98,99]: where and ∆W n are all independent and identically distributed random variables drawn from a common distribution p(∆W ), such that E(∆W ) = 0 and E(∆W 2 ) = ∆t.For example, a popular choice is We recall that a discrete-time approximation x n is said to converge weakly to x(t) if for all polynomials q(z) [98]: In practice, weak convergence implies the convergence of all moments in the ∆t → 0 limit.The order of convergence m is defined by the order of the error in the moments with the step size: for sufficiently small ∆t [98].In the case of heterogeneous diffusion instead, we considered the Itô integration scheme and used the Euler-Maruyama method.By taking into account the possible interpretations of the Langevin equation, the method is implemented as Finally, to compute the functional defined by (2) in the main text, we approximate where T is the first-passage time outside a given interval Ω, and N is the random number of steps needed for the approximated trajectory to exit from Ω.

FIG. 1 .
FIG. 1. Example of first-passage functionalZ = T 0 F [x(t)]dt for Brownian motion in a potential V (x) = V0 ln(x) diffusing in Ω = ( 1 2 ,52 ).Here F (x) = 1/x and Z thus corresponds to the area under the graph of 1/x(t), where x(t) is the stochastic trajectory displayed in the inset.The trajectory starts from x0 = 1 and the motion is evolved up to the first-exit time from Ω.The diffusion constant is equal to one and V0 = 0.5.

FIG. 4 .
FIG. 4. PDF of the scaled variable Z/ZD, with Z representing the functionals (a) Z = T 0 x(t)dt and (b) Z =

FIG. 5 .
FIG. 5. PDF of (a) Z = T 0 y(t)dt and (b) Z = T 0 [y(t)] −1/2 dt for trajectories in the positive real axis generated by Eq. (117) and evolved up to the first-passage time to the origin.The starting point is y0 = 1 and the space-dependent diffusion coefficient is D(x) = x 2/3 .The solid black curves are the theoretical predictions given by Eq. (128), the symbols represent the numerical distributions obtained from simulations.For both values of γ we considered three different interpretations, corresponding to Itô (α = 0, green triangles), Stratonovich (α = 1 2 , blue circles) and Hänggi-Klimontovich (α = 1, red squares).The insets, which are presented in log-log scale in (a) and semilog scale in (b), display the heavy-tails of the distributions.All datasets are obtained by evaluating Z over 10 6 trajectories, with small time step ∆t = 10 −4 .The integration scheme used to integrate Eq. (117) is the Itô scheme, and the desired interpretation is obtained by adding the appropriate drift µα(x) = αD ′ (x).