Steady state, relaxation and first-passage properties of a run-and-tumble particle in one-dimension

It is useful to define the Fourier transforms $\tilde{P}_\pm(k)=\int_{-\infty}^\infty P_\pm(x, t) {\rm e}^{i k x}~{{\rm d}x}$. Fourier transforming equation (3b) with respect to x, we obtain (in matrix form):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & \frac{\rm d}{{\rm d}t} \left(\begin{array}{@{}l@{}} \tilde{P}_{+}(k,t) \nonumber \\ \tilde{P}_{-} (k,t) \end{array}\right) = \mathbb{A}_k \left(\begin{array}{@{}l@{}} \tilde{P}_{+} (k,t) \nonumber \\ \tilde{P}_{-} (k,t) \end{array}\right), \label{Ppm_Eqn} \nonumber \end{align} \tag{ 6 }$$

where

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \mathbb{A}_k= \left(\begin{array}{@{}cc@{}} -1-{\rm i}k-\mathcal{D}k^2 & \gamma \nonumber \\ \gamma & -1+{\rm i}k-\mathcal{D}k^2 \end{array}\right). \nonumber \end{align*}$$

Diagonalizing the matrix $\mathbb{A}_k$ for each k and solving the resulting linear equations gives

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left(\begin{array}{@{}l@{}} \tilde{P}_{+}(k,t) \nonumber \\ \tilde{P}_{-} (k,t) \end{array}\right) &= \mathbb{W}_k \left(\begin{array}{@{}l@{}} \tilde{P}_{+} (k,0) \nonumber \\ \tilde{P}_{-} (k,0) \end{array}\right), \label{Ppm_Eqn-2} \nonumber \end{align} \tag{ 7 }$$

where

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \mathbb{W}_k &= \left(\begin{array}{@{}cc@{}} \frac{{\rm e}^{\alpha_+ t}}{\mathcal{N}_+^2} + \frac{{\rm e}^{\alpha_-t}}{\mathcal{N}_+^2} \left({\rm i}k+\sqrt{1-k^2} \right)^2 & \frac{\gamma \left({\rm e}^{\alpha_+ t} - {\rm e}^{\alpha_- t} \right)}{2 \sqrt{1 - k^2}} \nonumber \\ \frac{\gamma \left({\rm e}^{\alpha_+ t} - {\rm e}^{\alpha_- t} \right)}{2 \sqrt{1 - k^2}} & \frac{{\rm e}^{\alpha_+ t}}{\mathcal{N}_-^2} + \frac{{\rm e}^{\alpha_-t}}{\mathcal{N}_-^2} \left({\rm i}k-\sqrt{1-k^2} \right)^2 \end{array}\right), \nonumber \end{align*}$$

with $\mathcal{N}_\pm =\sqrt{2}~(1 -k^2 \pm {\rm i}k \sqrt{1-k^2}){}^{1/2}$ and $\alpha_{\pm}=-(1 +\mathcal{D}~k^2) \pm \sqrt{1-k^2}$.

Consider the natural initial condition in which the particle starts at x  =  0, with equal probability to be either in the  +  or the  −  state. The Fourier transform of the initial probability is then $\tilde{P}(k, 0)= (4 \pi){}^{-1}$. Using this in (7) and simplifying, we find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:Pk} \qquad \tilde{P_{\pm}} (k,t) = \frac{{\rm e}^{\alpha_+ t}}{2} \left(\frac{1}{\mathcal{N}_\pm^2}+\frac{1}{2\sqrt{1-k^2}} \right) + \frac{{\rm e}^{\alpha_- t}}{2} \left(\frac{\left({\rm i}k+\sqrt{1-k^2} \right)^2}{\mathcal{N}_\pm^2}-\frac{1}{2\sqrt{1-k^2}} \right). \nonumber \end{align} \tag{ 8 }$$

From (8), the Fourier transform of the total probability $\tilde{P}(k, t) = \tilde{P}_+(k, t)+\tilde{P}_-(k, t)$ is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{totalP_infiniteline} \tilde{P} (k,t) = {{\rm e}^{-\left(1 +\mathcal{D}k^2\right) t}} \left[ {\rm cosh} \left(t\sqrt{1-k^2} \right) + \frac{1 }{\sqrt{1-k^2}}\, {\rm sinh} \left(t\sqrt{1-k^2} \right) \right]. \nonumber \end{align} \tag{ 9 }$$

We can alternatively derive this result as follows: The displacement of an RTP that starts at x  =  0, can be written formally by integrating the Langevin equation (1) to give $\newcommand{\e}{{\rm e}} x(t) = \int_0^t \sigma(t) {\rm d}t + \sqrt{2\mathcal{D}} \int_0^t \xi(t) {\rm d}t\equiv A(t) +B(t)$. Since the random processes $A(t)$ and $B(t)$ are independent of each other, $\tilde{P}(k, t) = \left\langle {\rm e}^{{\rm i}kA(t)} \right\rangle \left\langle {\rm e}^{{\rm i}kB(t)} \right\rangle$. It is easy to see that the expression in (9) is actually in this product form once one identifies $\left\langle {\rm e}^{{\rm i}kB(t)} \right\rangle = {\rm e}^{-\mathcal{D} k^2 t}$ for the Brownian motion $B(t)$. The process $A(t)$ is the motion of an RTP with $\mathcal{D}=0$, whose dynamics is described by the telegrapher's equation, for which $\left\langle {\rm e}^{{\rm i}kA(t)} \right\rangle$ can be computed explicitly (see e.g. [31, 32]). This immediately leads to the expression in (9).

Using the product structure of $\tilde{P}(k, t)$, we invert the Fourier transform in (9) to derive the probability $P(x, t)$ in the convolution form $P(x, t)= \int_{-\infty}^\infty g(x-y, t)h(y, t)$ where $\newcommand{\e}{{\rm e}} g(x, t)= \exp(-x^2/4\mathcal{D}t)/\sqrt{4 \pi \mathcal{D} t}$ is the inverse Fourier transform of $\left\langle {\rm e}^{{\rm i}kB(t)} \right\rangle = {\rm e}^{-\mathcal{D} k^2 t}$ and $h(x, t)$ is the inverse Fourier transform of $\left\langle {\rm e}^{{\rm i}kA(t)} \right\rangle$. Using the explicit expression of $h(x, t)$ from [31, 32], we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P(x,t)&=\frac{{\rm cosh}(x/2\mathcal{D})}{\sqrt{4 \pi \mathcal{D} t}}{\rm e}^{-t-(x^2+t^2)/4\mathcal{D}t} \nonumber \\ & + \frac{{\rm e}^{-t}}{2} \int_{-\infty}^\infty \!\!{\rm d}y \, \frac{{\rm e}^{-(x-y)^2/4\mathcal{D}t}}{\sqrt{4 \pi \mathcal{D} t}}\left [I_0\left (\sqrt{t^2-y^2} \right) +\frac{t}{\sqrt{t^2-y^2}}I_1\left(\sqrt{t^2-y^2}\right) \right ]\Theta(t-|y|), \label{p-inf} \nonumber \end{align} \tag{ 10 }$$

where I_n is the nth-order modified Bessel function of the first kind and Θ is the Heaviside step function. Note that in the limit $x, t \to \infty$, $P(x, t)$ reduces to a simple Gaussian with diffusion constant $(\mathcal{D}+1/2)$. This can be easily seen from (9) where $\newcommand{\e}{{\rm e}} \tilde{P}(k, t) \to \exp [- (\mathcal{D}+1/2)k^2 t ]$ as $t \to \infty$ and $k\to 0$.

It is instructive to examine the spatial moments of the probability distribution. All odd moments are zero by symmetry. Formally, the even moments of the distribution are given by

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \langle x^{2n}(t)\rangle = (-1)^n\frac{\partial^{2n} \tilde{P}(k,t)}{\partial k^{2n}}\Big|_{k=0}. \nonumber \end{align*}$$

For the second moment, we find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{xt} \langle x^{2}(t)\rangle = (2\mathcal{D} + 1) \, t \; - \frac{(1-{\rm e}^{-2 t})}{2}\, \to (2D + v^2/\gamma) \, t \; - \frac{v^2}{2\gamma^2} (1-{\rm e}^{-2 \gamma t}) \nonumber \end{align} \tag{ 11 }$$

where in the last simplification we have put in all the dimensional parameters. The above result has two non-trivial limiting cases. For $\gamma \neq 0$ and $t\to\infty$, (11) reduces to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{xta} \langle x^{2}(t)\rangle\simeq (2\mathcal{D}+1)t \to (2D+v^2/\gamma)t. \nonumber \end{align} \tag{ 12 }$$

In the $t\to\infty$ limit, the finite switching rate γ leads to an enhancement of the microscopic diffusion coefficient in a manner that is reminiscent of hydrodynamic dispersion [42–45]. On the other hand, in the limit $\gamma\to 0$, we find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \langle x^{2}(t)\rangle\to 2\mathcal{D}t + v^2 t^2. \nonumber \end{align} \tag{ 13 }$$

Thus the mean-square displacement crosses over from growing linearly with t to quadratically with t as $\gamma\to 0$.

We can also compute higher-order derivatives of $\tilde{P}(k, t)$ from which higher moments of the displacement can be deduced. The fourth moment is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \langle x^4(t)\rangle = 3t^2\left(2D + \frac{v^2}{\gamma}\right)^2 - \frac{3v^2 \, t}{\gamma^3}\left[2D\gamma \left(1 - {\rm e}^{-2\gamma t}\right) - v^2\left(2 + {\rm e}^{-2\gamma t}\right)\right] +\frac{9 v^2}{2\gamma^4}\left(1 - {\rm e}^{-2\gamma t}\right). \nonumber \end{align} \tag{ 14 }$$

The important feature of this last result is that as $t\to\infty$, $\langle x^4\rangle/ 3 \langle x^2\rangle^2\to 1$, which is just the relation between the fourth and second moments for a Gaussian distribution. The behavior of the higher moments also conforms to those of the Gaussian distribution as $t\to\infty$.

In figure 1, we plot the temporal evolution of this occupation probability $P(x, t)$ for different values of the dimensionless diffusion coefficient $\mathcal{D}$. For $\mathcal{D}$ greater than a critical value $\mathcal{D}_c$, the probability distribution is unimodal for all times. However, for $\mathcal{D} < \mathcal{D}_c~(\approx 0.175)$, the occupation probability evolves from a unimodal distribution at short times, to a multimodal distribution at intermediate times and finally back to a unimodal distribution at long times. This non-trivial behavior of $P(x, t)$ for small $\mathcal{D}$ arises from the competition between the stochastic flipping of particle states (at rate γ) and translational diffusion. For a small $\mathcal{D}$, since $P_{\pm}(x, 0)=\delta(x)/2$, the RTPs in the  +  (−) states move to the right (left) in an almost ballistic manner. This splits the initial unimodal distribution into a bimodal distribution with two symmetric peaks (see figure 1). While these two peaks are moving ballistically in opposite directions, they are also broadening because of the true diffusion term. As a result, at intermediate times when the tails of these two separated peaks meet at the centre, there again starts accumulation of particles (see figure 1(a)). This may lead to a central peak before the two ballistically moving side peaks disappear, which depends on the relative strengths of γ and $\mathcal{D}$. Once developed, the central peak starts continuously broadening and on time scales much longer than the stochastic flipping rate $\gamma^{-1}$, the RTPs remix, leading to an effective diffusion constant as discussed above. As a result, the multimodal distribution at intermediate times merges into a unimodal distribution, which as $t\to \infty$, converges to a Gaussian. On the other hand, for large $\mathcal{D}$, the split peaks of the two RTP states overlap to such an extent that the full distribution always remains unimodal. This behavior suggests that there exists a critical $\mathcal{D}$ where the effects of translational diffusion and stochastic flipping balance each other.

Zoom In Zoom Out Reset image size

To understand this transition, we plot $P(x=0, t)$ for various $\mathcal{D}$ values in figure 2(a), and notice that the occupation probability at x  =  0 is higher for smaller $\mathcal{D}$ at short times as compared to that for larger $\mathcal{D}$. This then crosses over to a lower value at intermediate times and finally becomes larger at long times. Furthermore, to investigate the nature of the occupation probability at x  =  0, we plot the second derivative $\partial_x^2 P(x, t) \vert_{x=0}$. We find that for small $\mathcal{D}$, $P(x=0, t)$ has a maximum at short times, crosses over to a minimum at intermediate times, and finally crosses over to a maximum again at long times. However, for the critical value of $\mathcal{D}_{c}\approx 0.175$, these two crossover times merge, resulting in a unimodal distribution. For $\mathcal{D} \geqslant 0.175$, we find that $\partial_x^2 P(x, t) \vert_{x=0}$ is always negative and hence $P(x, t)$ is always unimodal.

Zoom In Zoom Out Reset image size

In summary, we find that, in contrast to the Gaussian form for a Brownian particle, the probability distribution for an RTP can be multimodal depending on the value of the dimensionless diffusion coefficient $\mathcal{D}$. This diversity in the probability distribution also occurs in other systems in which the motion of a diffusing particle is influenced by an interplay with a convection field that changes sign [48].

We now treat an RTP in the interval $\newcommand{\e}{{\rm e}} x \in [-\ell, \ell]$. In this case, the probability distribution will reach a steady-state in the long time limit. For $v=0$, the particle performs pure Brownian motion and reaches a spatially uniform steady-state at long times. On the other hand if $\mathcal{D}=0$, the particle is subjected to only the dichotomous noise $\sigma(t)$. Here, the particle reaches a different steady-state in which, for any finite flipping rate, there is an accumulation of particles at the boundaries. However, for very large flipping rates, one regains a spatially uniform distribution with an diffusion effective coefficient $v^2/\gamma$. In the case where both $v$ and $\mathcal{D}$ are nonzero, we anticipate a steady-state which is intermediate to these two extreme cases. We first solve for the probability distribution in the steady-state and then we turn to the more complicated time-dependent solution.

In the steady-state, the dimensionless Fokker–Planck equation (3b) reduce to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:Pss} \begin{array}{@{}l@{}} \mathcal{D} \, \partial_x^2 P_{+} - \partial_x P_{+} - P_{+} + P_{-} =0 , \nonumber \\ \mathcal{D} \, \partial_x^2 P_{-} + \partial_x P_{-} + P_{+} -P_{-} =0 \, \end{array} \nonumber \end{align} \tag{ 15 }$$

which we have to solve subject to the boundary conditions (5). The details of this calculation are given in appendix A. The final result for the probability distribution is:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P(x)={\left[\frac{\tanh\left(\frac{\sqrt{2\mathcal{D}+1} }{\mathcal{D}} \ell\right)}{\sqrt{2\mathcal{D}+1}} +2\ell\right]}^{-1}\left[\frac{\cosh\left(\frac{\sqrt{2\mathcal{D}+1}}{\mathcal{D}} x\right)}{2 \mathcal{D} \cosh\left(\frac{\sqrt{2\mathcal{D}+1}}{\mathcal{D}} \ell \right)} +1\right]. \label{expression_P} \nonumber \end{align} \tag{ 16 }$$

In figure 3, we compare (16) for the steady-state probability distribution $P(x)$ with results of simulation of the Langevin equation (1), and find nice agreement. We observe that probabilities are higher near the boundaries than at the center of the interval, in contrast to the uniform density which one would observe if there was no activity i.e. $v=0$. Such accumulation of active particles near the boundaries of a confined domain is quite generic and has been observed in experimental systems such as motile rods [49] and bacterial suspensions [50].

Zoom In Zoom Out Reset image size

In the limit $\mathcal{D} \to 0$, the peaks near the boundaries become progressively sharper, and eventually become delta-function peaks. The full distribution is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P(x)|_{\mathcal{D} \to 0} =\frac{2 + \delta(x-\ell)+\delta(x+\ell) }{2(1 + 2\ell)}. \nonumber \end{align} \tag{ 17 }$$

We observe that the probability is uniform everywhere except for the delta function peaks at the boundaries. This $\mathcal{D} \to 0$ case has recently been considered in a similar context [40] and our method reproduces their results.

Let us now turn to the full time-dependent solution. We are interested in how the distribution $P(x, t)$ approaches the steady-state in the $t \to \infty$ limit. To this end, we have to solve the coupled time-dependent Fokker–Planck equation (3b) within the interval $\newcommand{\e}{{\rm e}} x \in [-\ell, \ell]$, subject to the boundary conditions (5). For the time-dependent solution, we expand $P(x, t)$ in terms of the complete set of basis functions as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left(\begin{array}{@{}l@{}} P_+(x,t) \nonumber \\ P_-(x,t)\end{array}\right) = \sum_{n} a_n {\rm e}^{\lambda_n t}~ \left(\begin{array}{@{}l@{}} \phi^+_n(x) \nonumber \\ \phi^-_n(x) \end{array}\right), \label{seriessol} \nonumber \end{align} \tag{ 18 }$$

where $\lambda_n$ are the eigenvalues and the eigenfunctions $[\phi^+_n(x), \phi^-_n(x)]$ satisfy

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:lambda} \begin{array}{@{}l@{}} \mathcal{D}~\partial_{xx} {\phi}^+_{n} - \partial_x {\phi}^+_{n} - {\phi}^+_{n} + {\phi}^-_{n} = \lambda_n {\phi}^+_{n}, \nonumber \\ \mathcal{D}~\partial_{xx} \phi^{-}_n + \partial_x {\phi}^{-}_n + {\phi}^{+}_n - {\phi}^{-}_n = \lambda_n {\phi}^-_{n}, \end{array} \nonumber \end{align} \tag{ 19 }$$

subject to the boundary conditions (5). The coefficients a_n are given in terms of the left eigenvectors $\langle \chi_n|=[\chi_n^+(x), \chi_n^-(x)]$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle a_n=\langle \chi_n|P(x,t=0)\rangle=\int_{-\ell}^\ell {{\rm d}x}~ \left[~\chi_n^+(x) P_+(x,t=0)+ \chi_n^-(x) P_-(x,t=0)~\right]. \nonumber \end{align} \tag{ 20 }$$

The left eigenvectors can be obtained as solutions of the adjoint Fokker–Planck operator. It can be shown that this has the same form as that in equation (19), with the sign of the $\partial_x$ term changed, and with Neumann boundary conditions for both $\chi_n^+(x)$ and $\chi_n^-(x)$.

The details of calculating the eigenstates $\phi_n^\pm(x)$ are given in appendix B. Here we compare the time-evolution of the probability density obtained from a numerical solution of the Fokker–Planck equation (3b) using the boundary conditions (5), with the spectral expansion given by (18). The ground state eigenvalue $\lambda_0=0$, and the corresponding eigenstate (the steady-state) is known exactly and given by equation (16). In figure 4, we show the time-evolution of $P(x, t)$ obtained from a numerical solution of equation (3b). The unimodal to bimodal crossover discussed in the unbounded system can be seen in the figure. We also observe that our numerical solution converges to the exact steady-state. All the other eigenvalues $\lambda_n, ~ n>0$ have to be found numerically from the zeros of the determinant of the matrix ${M}$ in (B.6). The long-time relaxation of the system to the steady-state would be determined by the eigenvalue with the largest non-zero real part, and the corresponding eigenfunction. As an illustrative example, we choose $\mathcal{D}=0.5$ and $\newcommand{\e}{{\rm e}} \ell=1$, in which case the dominant eigenvalue is given by $\lambda_1 \approx -1.556\, 84$, while the corresponding eigenfunction is given by (B.5), with

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle &(\beta^{(1)}_{++},\beta^{(1)}_{-+},\beta^{(1)}_{+-},\beta^{(1)}_{--}) = (1.664\,09,-1.664\,09,-0.998\,284 i, 0.998\,284 i), \nonumber \\ & (\alpha^{(1)}_{++},\alpha^{(1)}_{-+},\alpha^{(1)}_{+-},\alpha^{(1)}_{--})~ \nonumber \\ &=(-0.277\,351, -3.605\,54, -0.058\,5543 - 0.998\,284 i, -0.058\,5543 + 0.998\,284 i), \nonumber \\ &(C^{(1)}_{++},C^{(1)}_{-+},C^{(1)}_{+-},C^{(1)}_{--}) \nonumber \\ &= (0.169\,039, 0.046\,8831,-0.196\,156 + 0.184\,987 i,-0.196\,156 - 0.184\,987 i). \nonumber \end{align*}$$

Zoom In Zoom Out Reset image size

Thus we know the functions $\phi_1^+(x)$, $\phi_1^-(x)$ explicitly, and in terms of these, we expect at long times

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P(x,t)=P_+(x,t)+P_-(x,t)=P_{SS}+ a_1 {\rm e}^{\lambda_1 t} [\phi_1^+(x)+\phi_1^-(x)]+\ldots.\label{spectralR} \nonumber \end{align} \tag{ 21 }$$

The parameter a₁ depends on initial conditions and can be obtained from the corresponding left eigenvector $\langle \chi_1 |$ and we find $a_1=\langle \chi|P(t=0)\rangle=0.7789...$. In figure (5(a)), we plot the evolution of $P(x, t)$ obtained from a direct numerical solution of equation (3b) and (5), starting from an initial condition $\delta(x-1/2)$. The plot in figure (5(b)) shows P(x,t)  −  P_SS and compares this with the prediction from the first term in the spectral representation equation (21). It is seen that agreement is very good. In general, we find that eigenvalues can have imaginary parts (in which case they come in complex conjugate pairs) and so one can see oscillatory relaxation.

Zoom In Zoom Out Reset image size

We are interested in the probability for an RTP, which starts from a point x on the semi-infinite line with velocity $\pm 1$, to arrive at the origin for the first time at time t. This quantity is directly related to the survival probability of the RTP in the same geometry in the presence of an absorbing boundary at x  =  0. Let S₊ (x,t) [S₋(x,t)] denote the probability that the RTP, starting initially at $x\geqslant 0$ with a positive [negative] velocity, survives being absorbed at the origin x  =  0 until time t, i.e. it does not cross the origin up to time t. Given the Langevin equation (1), it is convenient to write the backward Fokker–Planck equations for the evolution of $S_{\pm}(x, t)$, where the initial position x is treated as a variable [53]. Consider the evolution over the time window $[0, t+{\rm d}t]$ and break it into two sub-intervals $[0, {\rm d}t]$ and $[{\rm d}t, t+{\rm d}t]$. It follows from equation (1) that in a small time ${\rm d}t$ following t  =  0 (i.e. during the first interval $[0, {\rm d}t]$), the position of the particle evolves to a new position $\newcommand{\e}{{\rm e}} x'= x+ v\sigma(0)\, {\rm d}t + \sqrt{2 {D}}\eta(0)\, {\rm d}t$, where $\sigma(0)$ and $\newcommand{\e}{{\rm e}} \eta(0)$ are the initial noises. For the subsequent evolution in the time interval $[{\rm d}t, t+{\rm d}t]$, the new starting position is then $x'$. Thus the survival probability satisfies the evolution equations

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{nbfp} \begin{array}{@{}l@{}} S_{+}(x, t+{\rm d}t) = (1-\gamma {\rm d}t) \langle S_+(x+ v\, {\rm d}t+ \sqrt{2D}\,\eta(0)\, {\rm d}t, t)\rangle + \gamma {\rm d}t~S_-(x,t), \nonumber \\ S_{-}(x,t+{\rm d}t) = (1-\gamma {\rm d}t) \langle S_-(x- v\, {\rm d}t+ \sqrt{2D}\,\eta(0)\, {\rm d}t, t)\rangle + \gamma {\rm d}t~S_+(x,t), \end{array} \nonumber \end{align} \tag{ 22 }$$

where the $\langle\, \rangle$ denotes the average over $\newcommand{\e}{{\rm e}} \eta(0)$. Expanding in Taylor series for small ${\rm d}t$, using the properties of $\newcommand{\e}{{\rm e}} \eta(0)$ and taking ${\rm d}t\to 0$ limit, one directly arrives at a pair of backward equations which read, in dimensionless units,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}l@{}} \label{eqnQ} \partial_t S_{+}(x,t)= - S_{+}+ S_{-} + \partial_x S_{+} + \mathcal{D}\, \partial_x^2 S_{+}, \nonumber \\ \partial_t S_{-}(x,t)= S_{+} - S_{-} - \partial_x S_{-} + \mathcal{D}\, \partial_x^2 S_{-}. \end{array} \nonumber \end{align} \tag{ 23 }$$

These equations are valid for $x\geqslant 0$, with the initial conditions $S_{\pm}(x, 0)= 1$ for all x  >  0. In addition, we need to specify the boundary conditions. As the starting point $x\to \infty$, it is clear that $S_{\pm}(x\to \infty, t)=1$, since the particle will surely not cross the origin in a finite time. In contrast, the boundary condition at x  =  0 is subtle: it depends on whether ${{\mathcal D}}=0$ or ${{\mathcal D}}>0$. Consider first the case ${{\mathcal D}}=0$, i.e. in absence of normal diffusion. In this case, if the particle starts at x  =  0 with a negative velocity, it will surely cross the origin in a finite time. Hence

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S_{-}(x=0,t)=0 \quad \quad {\rm when} \quad {{\mathcal D}}=0 \, . \label{cald0} \nonumber \end{align} \tag{ 24a }$$

However, note that if the particle starts with a positive velocity, it can survive up to finite t, hence the boundary condition S₊ (0,t) is unspecified. We will see below that just one boundary condition in equation (24a) is sufficient to make the solution of equation (23) unique. Under normal diffusion, it is well known that if a particle crosses the origin at some time, it recrosses it immediately infinitely often [52]. Hence, if the particle starts at the origin, no matter whether the initial velocity is positive or negative, it will surely cross zero within a short time ${\rm d}t$, provided ${{\mathcal D}}>0$. This follows from the fact that in the ${\rm d}t \to 0$ limit, the Brownian noise dominates over the drift term irrespective of its sign. Hence, in this case, we have the two boundary conditions

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S_{\pm}(0, t) =0 \quad \quad {\rm when} \quad {{\mathcal D}}>0 \, . \label{caldnonzero} \nonumber \end{align} \tag{ 24b }$$

We will see later that indeed for ${{\mathcal D}}>0$, we will need both boundary conditions in equation (24b) to fix the solutions of equation (23) uniquely.

It is convenient to first take a Laplace transform, with respect to time t, of the pair of equation (23). Using the initial conditions $S_{\pm}(x, 0)=1$, it is easy to see that the Laplace transforms satisfy

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{lbfp} \begin{array}{@{}l@{}} -1 + s\, \tilde{S}_{+}(x,s) = {{\mathcal D}}\, \partial_x^2 \tilde{S}_{+} + \partial_x \tilde{S}_{+} - \tilde{S}_{+}+ \tilde{S}_{-} , \nonumber \\ -1 + s\, \tilde{S}_{-}(x,s) = {{\mathcal D}}\, \partial_x^2 \tilde{S}_{-} - \partial_x \tilde{S}_{-} - \tilde{S}_{-}+ \tilde{S}_{+} , \end{array} \nonumber \end{align} \tag{ 25 }$$

where $\tilde{S}_{\pm}(x, s)= \int_0^{\infty} {\rm d}t\, {\rm e}^{-st}\, S_{\pm}(x, t)$ is the Laplace transform.

These equations can be made homogeneous by the shift: $\tilde{S}_{\pm}(x, s)= 1/s + U_{\pm}(x, s)$, where $U_{\pm}$ satisfy

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{upm1} \begin{array}{@{}l@{}} \left[{{\mathcal D}}\partial_x^2 + \partial_x - (1+s)\right] U_{+}(x,s) = - U_{-}(x,s), \nonumber \\ \left[{{\mathcal D}}\partial_x^2 - \partial_x - (1+s)\right] U_{-}(x,s) = - U_{+}(x,s). \end{array} \nonumber \end{align} \tag{ 26 }$$

Furthermore, by differentiating twice, one can write closed equations for U₊ and U₋

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{upm2} \begin{array}{@{}l@{}} \left[{{\mathcal D}}\partial_x^2 - \partial_x - (1+s)\right]\left[{{\mathcal D}}\partial_x^2 + \partial_x - (1+s)\right]U_+(x,s)= U_+(x,s) , \nonumber \\ \left[{{\mathcal D}}\partial_x^2 + \partial_x - (1+s)\right]\left[{{\mathcal D}}\partial_x^2 - \partial_x - (1+s)\right]U_-(x,s)= U_-(x,s). \end{array} \nonumber \end{align} \tag{ 27 }$$

Below, we first solve the simpler case ${{\mathcal D}}=0$, followed by the more complex ${{\mathcal D}}>0$ case.

4.1.1. The case ${{\mathcal D}}=0$.

This particular case has been considered earlier with space-dependent transition rates, where only the mean first-passage time was computed [56]. Here we are interested in the full first-passage time distribution, which we can obtain using the above backward Fokker–Plank equation approach. For ${{\mathcal D}}=0$, equation (27) are ordinary second-order differential equations with constant coefficients. Hence, we can try solutions of the form: $U_{\pm}(x, s)\sim {\rm e}^{-\lambda x}$. Substituting this in either of equation (27), we find that λ satisfies the quadratic equation, $(\lambda+1+s)(\lambda-1-s)=1$, which gives two roots: $\lambda(s)= \pm \sqrt{s^2+2s}$. Obviously, the negative root is not admissible, since the solution must remain finite as $x\to \infty$. Retaining only the positive root, the general solutions of equation (27) can be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle U_{+}(x,s) = B\, {\rm e}^{-\lambda(s)\, x};\quad U_{-}(x,s) = A\, {\rm e}^{-\lambda(s)\, x}\quad \, {\rm where}\quad \lambda(s)=\sqrt{s^2+2s}. \label{AB} \nonumber \end{align} \tag{ 28 }$$

The two unknown constants B and A are however related, as they must also satisfy the pair of first-order equation (26) (upon setting ${{\mathcal D}}=0$). This gives $B= A/[1+s+\lambda(s)]$. Hence, finally, using $\tilde{S}_{\pm}(x, s)=1/s+ U_{\pm}(x, s)$, we get

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{spm1} \begin{array}{@{}l@{}} \tilde{S}_{+}(x,s) = \frac{1}{s}+ \frac{A}{1+s+\lambda(s)}\, {\rm e}^{-\lambda(s)\, x}, \nonumber \\ \tilde{S}_{-}(x,s) = \frac{1}{s}+ A\, {\rm e}^{-\lambda(s)\, x}, \end{array} \nonumber \end{align} \tag{ 29 }$$

where $\lambda(s)=\sqrt{s^2+2s}$. It remains to fix the only unknown constant A. This is done by using the boundary condition $\tilde{S}_{-}(0, s)= 0$ which fixes A  =  −1/s. Hence, we obtain the final solutions

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{spm2} \begin{array}{@{}l@{}} \tilde{S}_{+}(x,s) = \frac{1}{s}\left[1- \frac{1}{1+s+\lambda(s)}\, {\rm e}^{-\lambda(s)\, x}\right] , \nonumber \\ \tilde{S}_{-}(x,s) = \frac{1}{s}\left[1- {\rm e}^{-\lambda(s)\, x}\right], \end{array} \nonumber \end{align} \tag{ 30 }$$

with $\lambda(s)=\sqrt{s^2+2s}$.

The first-passage time distribution is simply related to the survival probability via

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_+(x,t)=-\partial_tS_+(x,t),\qquad f_-(x,t)=-\partial_tS_-(x,t), \nonumber \end{align} \tag{ 31 }$$

or in Laplace variables

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{fpm1} \begin{array}{@{}l@{}} \tilde{f}_+(x,s)= 1-s \tilde{S}_+(x,s)= \frac{1}{(s+1+\lambda(s))} \,\,{\rm e}^{-\lambda(s) x}, \nonumber \\ \tilde{f}_-(x,s)= 1-s \tilde{S}_-(x,s)={\rm e}^{-\lambda(s) x}, \end{array} \nonumber \end{align} \tag{ 32 }$$

where we recall $\lambda(s)=\sqrt{s^2+2s}$.

It turns out that the Laplace transforms in equation (32) can be exactly inverted. Before doing so, it is useful to extract the long-time asymptotics directly from the Laplace transforms in equation (32), by considering the $s\to 0$ limit. A scaling limit then naturally emerges where $s\to 0$, $x\to \infty$ but with the product $x\sqrt{s}$ fixed. This corresponds, in the time domain, to the scaling limit $t\to \infty$, $x\to \infty$, but keeping the ratio $x/\sqrt{t}$ fixed. In this limit, $\lambda(s)=\sqrt{s^2+2s}\to \sqrt{2s}$ as $s\to 0$. Then, using the Laplace inversion

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{L}^{-1}\left[ {\rm e}^{-a \sqrt{s}}\right]= \frac{a}{\sqrt{4\pi t^3}}\,\, {\rm e}^{-a^2/{4t}} \label{li.1} \nonumber \end{align} \tag{ 33 }$$

we find that both $f_{\pm}(x, t)$ converge, in the scaling limit, to the Holtsmark distribution

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_{\pm}(x,t) \to \frac{x}{\sqrt{4\pi\, D_0\, t^3}}\, {\rm e}^{-x^2/{4\,D_0\,t}}\quad {\rm where}\quad D_0=\frac{1}{2}. \label{fpm.2} \nonumber \end{align} \tag{ 34 }$$

Now we recall that for a Brownian particle evolving via $\newcommand{\e}{{\rm e}} {{\rm d}x}/{\rm d}t= \sqrt{2D_0}\, \eta(t)$, the first-passage probability $f(x, t)$ is given precisely [51] by the formula in equation (34). Hence, for our RTP that evolves via the telegraphic noise in equation (1) with ${{\mathcal D}}=0$, its first-passage probability in the scaling limit is equivalent to that for a normally diffusing particle with diffusion constant D₀  =  1/2. Indeed, this result is also consistent with our findings in equation (12), where we showed that, for ${{\mathcal D}}=0$, $\langle x^2(t)\rangle \to t$, which also corresponds to an effective normal diffusion at late times, with diffusion constant D₀  =  1/2.

To find the behavior of $f_{\pm}(x, t)$ for finite t, we need to invert the Laplace transforms in equation (32) exactly. Fortunately, this can be done using the following Laplace inversions

$$\begin{align*} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \mathcal{L}^{-1} \left(\frac{{\rm e}^{-x \lambda(s)}}{\lambda(s)[s+1+\lambda(s)]}\right) = &~{\rm e}^{-t}\frac{\sqrt{t-x}}{\sqrt{t+x}}~ I_1 \big(\sqrt{t^2-x^2}\big), \nonumber \\ \mathcal{L}^{-1}\left(\frac{{\rm e}^{-x\lambda(s)}}{\lambda(s)}\right) =&~ {\rm e}^{- t} \,I_0 \big(\sqrt{t^2-x^2}\big)\, , \nonumber \end{align*}$$

where I_0,1(t) are modified Bessel functions and $\lambda(s)=\sqrt{s^2+2s}$. Taking the derivative with respect to x, we obtain

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{fptsol} \begin{array}{@{}l@{}} f_+(x,t)=\frac{{\rm e}^{-t}}{t+x}\left[ x ~I_0\big(\sqrt{t^2-x^2} \big) + \frac{\sqrt{t-x}}{\sqrt{t+x}}~I_1\big(\sqrt{t^2-x^2} \big)\right]\theta(t-x), \nonumber \\ f_-(x,t)={\rm e}^{- t} \frac{x}{\sqrt{t^2-x^2}} \,I_1 \big(\sqrt{t^2-x^2} \big)~ \theta(t-x) + {\rm e}^{- t} \delta(t-x). \end{array} \nonumber \end{align} \tag{ 35 }$$

These results match with those obtained by Orsingher [57] using a different approach.

In figure 6 we verify these results for the first-passage time distributions with simulations. For any given x, the large t limit of (35) can be taken by using the asymptotic behavior $I_1(z)\to e^z/\sqrt{2\pi z}$ as $z\to \infty$. This yields, as $t\to\infty$, for any x,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_-(x,t)\simeq \frac{1}{\sqrt{2\pi}}\frac{x}{t^{3/2}}\, \quad {\rm{and}}\quad f_+(x,t)\simeq \frac{1}{\sqrt{2\pi}}\frac{x+1}{t^{3/2}}. \label{asym-f} \nonumber \end{align} \tag{ 36 }$$

Zoom In Zoom Out Reset image size

While the tail of f₋(x,t) behaves exactly as in the case of Brownian diffusion with a diffusion coefficient D₀  =  1/2 with a starting position x, the tail of f₊ (x,t) is equivalent to that in a Brownian diffusion with a starting position x  +  1. The extra length $1 \, (=v/\gamma)$ is the average distance the RTP with a positive velocity moves before taking the first turn. In figure 6 we compare the asymptotic results of (36) with numerical simulations and find very good agreement. From (36), the large time behavior of the survival probabilities are given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S_-(x,t) \sim \frac{x}{\sqrt{t}}\quad {\rm{and}}\quad S_+(x,t)\sim \frac{x+1}{\sqrt{t}}. \nonumber \end{align} \tag{ 37 }$$

In comparison with a particle with the negative starting velocity, a particle with the positive starting velocity has a higher probability of survival.

4.1.2. The case ${{{\mathcal D}}\ne 0}$.

As explained in the beginning of this subsection, the survival probabilities $S_{\pm}(x, t)$ satisfy the boundary conditions in equation (24b), i.e. in the Laplace domain, $\tilde{S}_{\pm}(x=0, s)=0$. In this case, the shifted functions $U_{\pm}(x, s)$ each satisfy a fourth-order ordinary differential equation with constant coefficients, as seen from equation (27). Trying again a solution of the form: $U_{\pm}(x, s)\sim {\rm e}^{-\lambda x}$, we find that λ has now 4 possible values that are the roots of the fourth-order polynomial

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[{{\mathcal D}}\lambda^2-\lambda-(1+s)\right]\left[{{\mathcal D}}\lambda^2+\lambda-(1+s)\right]=1\, . \label{poly4} \nonumber \end{align} \tag{ 38 }$$

There are 4 solutions given by $\pm \lambda_1(s)$ and $\pm \lambda_2(s)$ where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{l12s} \begin{array}{@{}l@{}} \lambda_1(s) = \left[\frac{1+{{\mathcal D}}(1+s)-\sqrt{1+4{{\mathcal D}}(1+s)+4{{\mathcal D}}^2}}{2{{\mathcal D}}^2}\right]^{1/2}, \nonumber \\ \lambda_2(s) = \left[\frac{1+{{\mathcal D}}(1+s)+\sqrt{1+4{{\mathcal D}}(1+s)+4{{\mathcal D}}^2}}{2{{\mathcal D}}^2}\right]^{1/2}\, . \end{array} \nonumber \end{align} \tag{ 39 }$$

Evidently, $\lambda_1(s)<\lambda_2(s)$. Again discarding the negative roots $-\lambda_1(s)$ and $-\lambda_2(s)$ (since the solution cannot diverge as $x\to \infty$), the general solutions of equation (27) can be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}l@{}} U_+(x,s) = B_1\, {\rm e}^{-\lambda_1(s)\, x} + B_2\, {\rm e}^{-\lambda_2(s)\, x} , \nonumber \\ U_{-}(x,s)= A_1\, {\rm e}^{-\lambda_1(s)\, x} + A_2\, {\rm e}^{-\lambda_2(s)\, x} , \end{array} \nonumber \end{align} \tag{ 40 }$$

where $\lambda_{1, 2}(s)$ are given in equation (39). However, these solutions must also satisfy the individual second-order equation (26). This indicates that A₁, A₂ are related to B₁ and B₂. Indeed, by substituting these solutions in equation (26) gives the following relations

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{ba12} \begin{array}{@{}l@{}} B_1 = -\frac{A_1}{{{\mathcal D}}\lambda_1^2-\lambda_1-(1+s)}= -A_1 \left({{\mathcal D}} \lambda_1^2+\lambda_1-(1+s)\right) , \nonumber \\ B_2 = -\frac{A_2}{{{\mathcal D}}\lambda_2^2-\lambda_2-(1+s)}= -A_2 \left({{\mathcal D}} \lambda_2^2+\lambda_2-(1+s)\right). \end{array} \nonumber \end{align} \tag{ 41 }$$

Note that we have used equation (38) to obtain the last two relations.

Hence, the solutions for the survival probabilities are given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}l@{}} \tilde{S}_+(x,s)= \frac{1}{s} +B_1\, {\rm e}^{-\lambda_1\, x}+ B_2\, {\rm e}^{-\lambda_2\, x}, \nonumber \\ \tilde{S}_-(x,s)= \frac{1}{s}+ A_1\, {\rm e}^{-\lambda_1\, x} + A_2\, {\rm e}^{-\lambda_2\, x}\, \end{array}\label{survpm1} \nonumber \end{align} \tag{ 42 }$$

where B₁, B₂ are related to A₁ and A₂ via equation (41). We are still left with two unknown constants A₁ and A₂. To fix them, we use the two boundary conditions: $\tilde{S}_{\pm}(x=0, s)=0$. This gives two linear equations for A₁ and A₂ whose solution is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{a12} \begin{array}{@{}l@{}} A_1 = -\frac{1}{s} \frac{{{\mathcal D}}\lambda_2^2 + \lambda_2 -s}{[{{\mathcal D}} (\lambda_2^2-\lambda_1^2)+ \lambda_2-\lambda_1]}, \nonumber \\ A_2 = \frac{1}{s} \frac{{{\mathcal D}}\lambda_1^2 + \lambda_1 -s}{[{{\mathcal D}} (\lambda_2^2-\lambda_1^2)+ \lambda_2-\lambda_1]}\, . \end{array} \nonumber \end{align} \tag{ 43 }$$

This then uniquely determines the solutions for the survival probabilities $\tilde{S}_{\pm}(x, s)$. The corresponding first-passage probabilities are given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{fpml2} \begin{array}{@{}l@{}} \tilde{f}_{+}(x,s)= 1-s \tilde{S}_{+}(x,s)= -s\left[ B_1\, {\rm e}^{-\lambda_1\, x} + B_2\, {\rm e}^{-\lambda_2\, x}\right] , \nonumber \\ \tilde{f}_{-}(x,s)= 1-s \tilde{S}_{-}(x,s)= -s\left[ A_1\, {\rm e}^{-\lambda_1\, x} + A_2\, {\rm e}^{-\lambda_2\, x}\right] , \end{array} \nonumber \end{align} \tag{ 44 }$$

where the constants A₁, A₂, B₁ and B₂ are determined explicitly above and $\lambda_{1, 2}(s)$ are given in equation (39).

The first nontrivial check is the limit ${{\mathcal D}}\to 0$. In this limit, it is easy to verify, from equation (39), that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \lambda_1(s) \to \lambda(s)=\sqrt{s^2+2s},\qquad {\rm and} \qquad \lambda_2(s)\to \frac{1}{{{\mathcal D}}}\to \infty\, . \label{lD0} \nonumber \end{align} \tag{ 45 }$$

In addition, one finds that as ${{\mathcal D}}\to 0$,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle A_1\to - \frac{1}{s}, \qquad {\rm and} \qquad A_2\to 0\, , \label{a1a2D0} \nonumber \end{align} \tag{ 46 }$$

and consequently

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle B_1\to -\frac{1}{s (1+s+ \sqrt{s^2+2s})}, \qquad {\rm and} \qquad B_2 \to 0\, . \label{b1b2D0} \nonumber \end{align} \tag{ 47 }$$

We therefore recover the ${{\mathcal D}}=0$ results in equation (32).

We now turn to the long-time asymptotic solutions of equation (44) for arbitrary ${{\mathcal D}}$. Hence we consider the $s\to 0$ limit, with finite ${{\mathcal D}}$. In this limit, it is easy to check that, to leading order for small s

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \lambda_1(s)\to \left[\frac{2}{1+2{{\mathcal D}}}\, s\right]^{1/2}, \qquad {\rm and}\qquad \lambda_2(s) \to \frac{\sqrt{1+2{{\mathcal D}}}}{{{\mathcal D}}}\, . \label{l1l2s0} \nonumber \end{align} \tag{ 48 }$$

Similarly, one can check that to leading order for small s

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle s\, A_1 \to -1, \qquad {\rm and}\qquad s\,A_2\to O(\sqrt{s}), \nonumber \end{align*}$$

and consequently

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle s\, B_1 \to -1 \qquad {\rm and}\qquad s\,B_2 \to O(\sqrt{s})\, . \nonumber \end{align*}$$

Substituting these results together in equation (44), we find that in the scaling limit ($s\to 0$, $x\to \infty$ with the product $x\sqrt{s}$ fixed)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{f}_{\pm}(x,s) \to \exp\left[- \sqrt{\frac{2s}{1+2{{\mathcal D}}}}\,\, x\right]\, . \label{fpmx.1} \nonumber \end{align} \tag{ 49 }$$

Upon inverting the Laplace transform using equation (33), we obtain our final results in the scaling limit

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_{\pm}(x,t) \to \frac{x}{\sqrt{4\pi\, D_1\, t^3}}\,\, {\rm e}^{-x^2/{4\,D_1\,t}}, \quad {\rm where} \quad D_1= {{\mathcal D}}+\frac{1}{2}\, . \label{fpmx.2} \nonumber \end{align} \tag{ 50 }$$

This result is precisely the same as the first-passage time density of an ordinary Brownian motion with diffusion constant $D_1= {{\mathcal D}} +1/2$. Note that for $D_1\to D_0=1/2$ as ${{\mathcal D}}\to 0$. Moreover, this effective diffusion constant $D_1= {{\mathcal D}}+1/2$ is consistent with our result $\langle x^2(t)\rangle \to (2{{\mathcal D}}+1)\, t$ in equation (12).

Unfortunately, unlike in the ${{\mathcal D}}=0$ case, for nonzero ${{\mathcal D}}$, we are not able to obtain the finite time result for $f_{\pm}(x, t)$ explicitly, due to the fact that the Laplace transforms are difficult to invert.

We now investigate an RTP on a finite interval and address two questions: (a) the probability for the particle, which starts at x, to eventually reach either of the boundaries, and (b) the mean time for the particle to exit the interval by either of the boundaries. Let E₊ (x) (E₋(x)) denote the exit probabilities, namely the probability for a particle that starts at x with velocity  +1 (−1) exits through the boundary at $\newcommand{\e}{{\rm e}} x=-\ell$. By comparison, the exit probability to $\newcommand{\e}{{\rm e}} x=-\ell$ for isotropic diffusion, $E(x)$, is simply $\newcommand{\e}{{\rm e}} \frac{1}{2}(1-\frac{x}{\ell})$; that is, the exit probability decreases linearly with the initial distance from the left edge.

It is easily seen that these hitting probabilities obey the backward equations [51]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}l@{}} \mathcal{D}\, \partial^2_{x}E_+ + \partial_x E_+ - (E_+-E_-) = 0 , \nonumber \\ \mathcal{D}\, \partial^2_x E_- - \partial_x E_- + (E_+-E_-) = 0 , \label{E+-} \end{array} \nonumber \end{align} \tag{ 51 }$$

subject to the appropriate boundary conditions, which are $\newcommand{\e}{{\rm e}} E_\pm(-\ell)=1$ and $\newcommand{\e}{{\rm e}} E_\pm(\ell)=0$. These boundary conditions fix the constants in $E_\pm$ and thus the problem is formally solved. The calculation is conceptually straightforward but tedious, and the details were performed by Mathematica. The basic steps and the final expressions for the exit probabilities are given by (C.5) in appendix C.

Figure 7 shows the exit probabilities $E_\pm(x)$ for representative values of the dimensionless diffusion coefficient $\mathcal{D}$. As one expects, for $\mathcal{D}\gg 1$, the exit probabilities are close to the isotropic random-walk form $\newcommand{\e}{{\rm e}} \frac{1}{2}(1-\frac{x}{\ell})$. However, for $\mathcal{D}\ll 1$, E₊ and E₋ become very distinct. Moreover, the exit probability E₊ decreases much more rapidly with x than $\newcommand{\e}{{\rm e}} (1-\frac{x}{\ell})/2$, while E₋ decreases much more slowly. Notice also that $\newcommand{\e}{{\rm e}} \mathcal{E}(x)\equiv \frac{1}{2}[E_+(x)+E_-(x)]\ne \frac{1}{2}(1-\frac{x}{\ell})$. That is, the exit probability, averaged over the two velocity states deviates significantly from the corresponding exit probability for unbiased diffusion.

Zoom In Zoom Out Reset image size

Let us now turn to the exit times. Let t₊ (x,t) (t₋(x,t)) be the mean first-passage time (to either boundary) for a particle that is at x and is also in the  +  (−) state. Again using the formalism given in [51], it is easily seen that these exit times obey the backward equations

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}l@{}} \label{t+-} \mathcal{D} \partial^2_{x}t_+ + \partial_x t_+ - (t_+-t_-) = -1 , \nonumber \\ \mathcal{D} \partial^2_x t_- - \partial_x t_- + (t_+-t_-) = -1 , \end{array} \nonumber \end{align} \tag{ 52 }$$

with boundary conditions $\newcommand{\e}{{\rm e}} t_\pm(\pm\ell)=0$, which corresponds to the particle being immediately absorbed if it starts at either end of the interval. The solution to equation (52) are obtained using the same approach as that given for equation (51) (see appendix D). While the resulting expressions for $t_\pm$ for the finite interval are too long to be displayed, the form of the first-passage times are easily visualized (figure 8). For small bias velocity the scaled diffusion constant $\mathcal{D}=D\gamma/v^2$ (see section 2) is large. In this limit the diffusive part of the dynamics for both types of the particles dominates. As a result, $t_+\approx t_-$, and both t₊ and t₋ become very close to the exit probability for isotropic diffusion $\newcommand{\e}{{\rm e}} t=(\ell^2-x^2)/2\mathcal{D}$. On the other hand, if $v$ is increased $\mathcal{D}$ is decreases and the active contribution to the motion dominates. In this limit the exit times t₊ and t₋ strongly deviate from each other.

Zoom In Zoom Out Reset image size