Run-and-Tumble particle in one-dimensional confining potential: Steady state, relaxation and first passage properties

# Run-and-Tumble particle in one-dimensional confining potential: Steady state, relaxation and first passage properties

Abhishek Dhar International Centre for Theoretical Sciences, TIFR, Bangalore 560089, India    Anupam Kundu International Centre for Theoretical Sciences, TIFR, Bangalore 560089, India    Satya N. Majumdar LPTMS, CNRS, Univ. Paris-Sud, Université Paris-Saclay, 91405 Orsay, France    Sanjib Sabhapandit Raman Research Institute, Bangalore 560080, India    Grégory Schehr LPTMS, CNRS, Univ. Paris-Sud, Université Paris-Saclay, 91405 Orsay, France
September 14, 2019
###### Abstract

We study the dynamics of a one-dimensional run and tumble particle subjected to confining potentials of the type , with . The noise that drives the particle dynamics is telegraphic and alternates between values. We show that the stationary probability density has a rich behavior in the -plane. For , the distribution has a finite support in and there is a critical line that separates an active-like phase for where diverges at , from a passive-like phase for where vanishes at . For , the stationary density collapses to a delta function at the origin, . In the marginal case , we show that, for , the stationary density is a symmetric exponential, while for , it again is a delta function . For the special cases and , we obtain exactly the full time-dependent distribution , that allows us to study how the system relaxes to its stationary state. In addition, in these two cases, we also study analytically the full distribution of the first-passage time to the origin. Numerical simulations are in complete agreement with our analytical predictions.

## I Introduction

There has been a surge of interest in understanding the statics and the dynamics of a class of non-equilibrium systems, commonly called “active matter”. Such systems appear in a variety of contexts ranging from self-catalytic colloids Fodor17 , living cells, biological processes, in the growth of microbial colonies Magistris15 , in the physics of self-motile systems Magistris15 , active gels Sriram10 , motility induced phase transition Cates15 to flocking of birds Toner05 . While the collective properties of such systems are of great interest Fodor17 ; bechinger_active_2016 ; cates_motility-induced_2015 , they exhibit a rich static and dynamical behaviour even at a single-particle level. This is due to the fact that the single-particle stochastic dynamics is typically non-Brownian: the effective noise in the Langevin equation of an active particle is typically “coloured”, i.e., the autocorrelation of the noise decays exponentially with a finite correlation time . The limit would correspond to the Brownian motion, while the limit corresponds to purely ballistic motion.

The most common example of such an active particle dynamics is the well known “persistent Brownian motion”, renamed recently as the “run and tumble particle” (RTP). In one-dimension, the overdamped stochastic dynamics of an RTP, starting at , is governed by the Langevin equation

 dxdt=v0σ(t), (1)

where the noise is a dichotomous or telegraphic noise which takes values and flips its sign with a constant rate , as shown in Fig 1. The auto-correlation function of this telegraphic noise is given by

 ⟨η(t1)η(t2)⟩=v20e−2γ|t1−t2|. (2)

In the limit , keeping the ratio fixed, the noise becomes a white noise with correlator and consequently represents the position of a Brownian particle. The probability distribution of the particle’s position is known exactly Othmer ; Weiss ; Martens . Very recently, the first-passage properties were also computed exactly Malakar2018 .

What happens if this RTP is subjected to an external potential ? The corresponding Langevin equation now reads

 dxdt=f(x)+v0σ(t), (3)

where is the external force. Once again, in the Brownian limit, when , keeping the ratio fixed, the probability distribution converges at long times to the equilibrium Boltzmann-Gibbs distribution . However, when the correlation time of the noise is finite, i.e. the noise has a finite memory, the stationary state is nontrivial and is typically non-Boltzmann. While an explicit formula is known for this , originally derived in the quantum optics literature Horsthemke84 ; Klyatskin78a ; Klyatskin78b ; Lefever80 ; Hanggi95 and re-derived more recently in the active-matter literature (see for instance the Appendix in Cates_Nature ), the structure of the stationary state has not been analysed in detail for different types of potentials. One of the goals of this paper is to study this stationary state for a class of confining potentials of the type with and . Indeed, keeping fixed, and varying and , we uncover an interesting phase diagram in the plane. For , we show that while the stationary solution has a finite support for any , there is a critical line

 α=αc(p)=(v0p)2−p(γp(p−1))p−1, (4)

such that for the stationary state distribution vanishes at the edges, while for the distribution diverges at the edges. Thus, while for the particle behaves more like an “active” one, it behaves more as a “passive” particle for . We call this transition from the active to the passive phase across a “shape transition” of the stationary profile. The behaviour of around depends on and is also interesting. We find that for , the function near the origin is convex (i.e., there is a local minimum at ) while for it is concave (i.e., there is a local maximum at the origin), irrespective of whether or . A weak fingerprint of activeness is thus manifest, for , even in the passive phase, leading to bimodal distributions (see for instance Fig. 7). For , we show that the RTP collapses to the origin at long times, leading to trivial stationary state . The case turns out to be the marginal case. In this case there is a nontrivial stationary state only for , while the stationary distribution is again for . These results lead to the phase diagram shown in Fig 2.

In the presence of the external confining potential , it is also equally interesting to study the dynamics of the system away from the stationary state, i.e., how the system relaxes to the stationary state at long times. For a simple hard-box potential, where the RTP is confined inside a finite box, e.g. when for while outside the box, this relaxation to the stationary state was recently studied in Ref. Malakar2018 . This case corresponds to the limit of the potential . It is thus interesting to study the relaxational dynamics for a finite value of . In this paper, we provide a detailed study of the relaxational dynamics for the case , i.e. the harmonic potential with . Interestingly, our results show that, if the particle starts exactly at the origin, i.e., , the inverse relaxation time characterizing the exponential decay of to its stationary distribution undergoes a transition at : for (passive phase) , while for (active phase) we show that (see Fig. 10 ). For a generic starting point , the inverse relaxation time is always given by . Thus if one wants to detect the signature of active to passive transition in the relaxational dynamics, one needs to start at the special initial position .

Another related interesting observable is the first-passage probability which denotes the probability density that the RTP, starting at , crosses the origin for the first time at time . For arbitrary confining , there exist extensive results in the literature on the mean first-passage time Lindenberg ; Weiss_Lind . However, there is hardly any result on the full distribution of the first-passage time. Only very recently exact results were obtained for a free RTP, i.e. for , as well as for the box potential Malakar2018 . In this paper, we present exact results for for two cases and .

The rest of the paper is organised as follows. In Section II, we provide a derivation of the steady state probability density for arbitrary confining potential and then discuss the form of the stationary distribution for potentials of the form , with . We show that there is a shape transition in the -plane from an active profile to a passive profile across a critical line . The resulting phase diagram is shown in Fig. 2. In Section III, we discuss the relaxation to the steady state for the harmonic potential . First-passage properties for this harmonic case are discussed in Section IV. We conclude in Section V with a summary and open questions. Some details have been relegated to two Appendices. In Appendix A we provide some details of the calculations of the first-passage properties for the harmonic case . In Appendix B we briefly discuss the stationary-state, relaxation to the stationary state as well as the first-passage properties for the marginal case where .

We consider the Langevin equation (3) for a particle moving in an arbitrary confining potential with the force and subjected to the telegraphic noise . For sufficiently confining potential, e.g. with , the system reaches a stationary state at long times. This steady state distribution is actually known for generic confining potential . It was initially derived in the quantum optics literature Horsthemke84 ; Klyatskin78a ; Klyatskin78b ; Lefever80 , later to study the role of colored noise in dynamical systems Hanggi95 and more recently it has been re-derived in the context of active particles systems Cates_Nature . Even though an explicit formula existed for the stationary state, the physical consequences of this formula has not been analysed in detail in the literature, to the best of our knowledge. In this section, we re-derive this formula and discuss its implications for confining potentials of the type . One of the interesting outcomes of our analysis is a novel shape transition in the plane, for a fixed (see Fig. 2).

We start by deriving the Fokker-Planck equation associated to the Langevin equation (3). Since the telegraphic noise in (3) takes only two values , it is natural to define and denoting the probability densities for the particle to be at position at time with velocity and , respectively. Clearly, the main object of interest is the full probability density . It is quite easy to see that these two densities and evolve according to

 ∂P+∂t = −∂∂x[(f(x)+v0)P+]−γP++γP−, (5) ∂P−∂t = −∂∂x[(f(x)−v0)P−]+γP+−γP−. (6)

where the first term in both the equation appear from the drift term in the Langevin Eq. (3), while the remaining terms appear due to the stochastic change in the direction of the velocity with rate . For sufficiently confining potentials, we expect the system to reach a stationary state which is obtained by setting the time derivative to be zero on the left hand side in Eqs. (5) and (6). We denote the stationary distributions (dropping the -dependence) which then satisfy a pair of coupled first-order differential equations

 ddx[(f(x)+v0)P+]+γP+−γP− = 0, (7) ddx[(f(x)−v0)P−]−γP++γP− = 0. (8)

To solve these differential equations, we need to specify the boundary conditions for and . For that, we first need to know where the boundaries are. In fact, it turns out that in the steady state, for sufficiently confining potentials, these active particles accumulate over a finite region of space around the minimum of the confining potential (assuming that it has a single minimum, such as in a harmonic potential). This can be easily seen by examining the Langevin equation (3). We first note that the right hand side of Eq. (3) has two fixed points such that

 f(x−)=v0,f(x+)=−v0,withx+>x−. (9)

For example, for a harmonic potential , and hence . Consider a particle with position . Then the drift is always towards the center of the trap, no matter what the sign of is, as long as (see Fig. 3). Thus, even if the particle starts with initial position , it will eventually arrive at a position . Similarly, for any particle with position , the drift is always towards the center, irrespective of the sign of . Hence, eventually, the position of the particle will be larger than , even if it starts to the left of (see Fig. 3). Thus it is clear that, in the steady state, the probability distribution will have a finite support where it is non-zero, while it vanishes outside this box (see Fig. 3). Therefore, the steady state equations in (7) and (8) are valid for and we need to specify the boundary conditions at . To find these boundary conditions, we note that if the particle is at with a negative velocity, it always moves to the left of . Therefore, in the stationary state, there can not be any particle at with a negative velocity, indicating that

 P−(x+)=0. (10)

Note that remains unspecified at . A similar argument at shows that

 P+(x−)=0, (11)

while remains unspecified at . With this pair of boundary conditions (10) and (11), the differential equations in Eqs. (7) and (8) admit a unique solution that we determine below.

To proceed, we add and subtract these two equations (7) and (8). Defining

 P(x)=P+(x)+P−(x), (12) Q(x)=P+(x)−P−(x), (13)

they satisfy, respectively, the equations

 ddx[f(x)P+v0Q]=0, (14) ddx[f(x)Q+v0P]+2γQ=0. (15)

Eq. (14) implies

 [f(x)P+v0Q]=C, (16)

where is a constant. To determine this constant, we now invoke the boundary conditions (10) and (11). Clearly, Eq. (16) holds for all . Setting we get

 [f(x−)+v0]P+(x−)+[f(x−)−v0]P−(x−)=C. (17)

Using from Eq. (9), and the boundary condition from Eq. (10), we get . Hence from Eq. (16), we get

 Q(x)=−1v0f(x)P(x). (18)

We now substitute this equation in Eq. (15) and obtain

 ddx[(v20−f2(x))P(x)]−2γf(x)P(x)=0. (19)

Solving this equation, we get

 P(x)=Av20−f2(x)exp(2γ∫x0dyf(y)v20−f2(y)), (20)

valid for . The overall constant in (20) is determined from the normalisation . This provides the general steady state solution for arbitrary , provided is normalisable. Below, we analyse this formula for the stationary distribution (20) for confining potentials of the type , i.e., with . We will see below that the situation is quite different for and , including the marginal case . We first start with the case, which corresponds to a harmonic potential .

Harmonic potential (). To analyse the implications of this solution in Eq. (20), let us first focus on the specific example of a harmonic potential and for using the standard notations of the harmonic oscillator we set where is the standard spring constant. In this case, using , and , one gets from Eq. (20)

 (21)

where is the beta function and the exponent

 ϕ=γμ−1=γ2α−1, (22)

where we used . Clearly, this solution is symmetric around and has a finite support over . We have verified this prediction in Eq. (21) numerically for two different values of (see Fig. 5). Interestingly, the shape of the solution in Eq. (21) depends on the exponent . If , the solution vanishes at the two ends of the support, while for , the solution diverges at the two ends. This shape transition occurs at which corresponds to a critical value . Indeed, the solution is concave-shaped for and convex-shaped for . To appreciate this shape transition, it is instructive to compare this solution for the active particle to that of a passive particle described by the Langevin equation (3) with the same force , but driven by a delta-correlated white noise. In this passive case, Eq. (3) just describes an Ornstein-Uhlenbeck (OU) process, whose steady state has the Boltzmann distribution

 POU(x)=√μ2πDe−μ2Dx2, (23)

valid over the full space. It is easy to check that the active solution in Eq. (21) actually approaches the Boltzmann distribution in (23) in the limit , but keeping the ratio fixed. Indeed, in the limit, the telegraphic noise approaches the delta-correlated noise and naturally the two stationary solutions also coincide. However, for finite , the active noise has an exponentially decaying memory and this leads to highly non-Boltzmann distributions in Eq. (21). While for , the stationary distribution still retains a single peak structure (as in the Boltzmann case), it undergoes a “shape transition” at . For , the stationary distribution has a double-peaked structure (with diverging peaks at ), with a minimum at the center of the trap at . Thus, the single-peak structure is a signature of a passive phase, while the double-peak structure signifies an active phase. By reducing , one crosses over from this passive phase () to an active phase (). Alternatively, this shape transition can also be viewed as a function of (for fixed ). There is a critical value at which the exponent in Eq. (22) changes sign. The active phase (where ) corresponds to and the passive phase (where ) corresponds to .

Finally, we note that Eq. (21) gives the expression for the total probability distribution . It is instructive to compute ans separately also, as they have rather different behaviours because the boundary conditions they satisfy are quite different [see Eqs. (10) and (11)]. Indeed from Eq. (18), for , we have

 P+(x)−P−(x)=Q(x)=μxv0P(x). (24)

This, along with , allows to express and in terms of in Eq. (21). This gives

 P±(x)=A2(1±μxv0)γμ(1∓μxv0)γμ−1, (25)

where is the normalization constant appearing in Eq. (21), fixed by . Note that the these distributions satisfy the boundary conditions in Eqs. (10) and (11) respectively. In fig. 4 we verify these expressions of numerically and observe a nice agreement.

Anharmonic confining potentials with . We now analyse the general stationary solution in Eq. (20) for the case , with . In this case . Since the general solution in Eq. (20) is symmetric around , i.e., , it is sufficient to consider for . Substituting with in Eq. (20), we get

 P(x)=Bpb2p−x2p−2exp[−2γαp∫x0yp−1b2p−y2p−2dy], (26)

where the parameter is given by

 bp=v0αp (27)

and the overall normalisation constant is such that . The solution is actually symmetric around with a finite support . We next analyse Eq. (26) to see how the solution behaves near the upper support end at . Clearly, by symmetry, it has the same behavior at the lower support . To proceed, we set where is small. The integral inside the exponential in Eq. (26) then reads

 I(ϵ)=∫b1p−1p−ϵ0yp−1b2p−y2p−2dy. (28)

We next make a change of variable in the integral in Eq. (28), and expand the integrand to leading order for small . Performing the integral gives, to leading order for small ,

 I(ϵ)≈−b2−pp−1p2(p−1)lnϵ. (29)

Substituting this result in Eq. (26), we then find that the stationary solution behaves, near the upper edge, as

 P(x=b1p−1p−ϵ)∼ϵϕ, (30)

with a non-universal (i.e., parameter dependent) exponent

 ϕ=b2−pp−1pαp(p−1)γ−1=1p(p−1)(v0p)2−pp−1γα1p−1−1. (31)

Clearly, for , the exponent , indicating that in Eq. (30) is integrable near the edge. For , setting , we recover as in Eq. (22). For any , clearly there is a shape transition in the plane (see Fig. 2) accross the line

 α=αc(p)=(v0p)2−p(γp(p−1))p−1 (32)

such that for and for . Thus, as in the case, there is a shape transition (from converging edges to diverging edges) for generic as increases through . The existence of this critical line in the plane, separating an active and a passive phase, is indeed the main new result of this section.

Indeed the shape-transition for , where the exponent in Eq. (31) changes sign, can be accessed by varying either or . Indeed from Eq. (31), the exponent only depends on the ratio . Thus, small (respectively large) corresponds to large (respectively small) for the same value of . While the phase diagram in Fig. 2 is presented in the -plane, for , we can equivalently probe this transition by varying , with fixed. In fact, in the numerical simulations for (Fig. 5) and (Fig. 6), we indeed kept fixed, while varying . In both cases, one clearly sees the transition and an excellent agreement between theoretical predictions and simulations.

The behaviour near is also interesting. Performing an analysis around , in the same was as was done near the right edge at , one finds for and small

 P(ϵ)≃Bpb2p[1−2γαp2b2p |ϵ|p+1b2p|ϵ|2(p−1)+O(|ϵ|3p−2)].

We see that, as , the dominant sub-leading term is either the second (if ) or the third term (if ). In the former case, the sign of the sub-leading term is negative, while in the latter case it is positive. Consequently, the distribution , near , is convex for whereas for it is concave, irrespective of whether is greater or less than . The behaviour for is demonstrated in Fig. 7 where the distribution is convex near for both and .

The marginal case where . In this case, setting in Eq. (20), we find that a non-localised stationary solution exists only for where it is given by

 P(x)=γαv20−α2exp(−2γαv20−α2|x|), −∞≤x≤∞. (33)

Note that, unlike in the case, the stationary solution is no-longer supported on a finite interval but extends over the full space. From this solution in Eq. (33), we also see that, when approaches from below, the stationary distribution approaches to a delta function centred at the origin . Indeed, for any , the system collapses to a delta function at , at long times. This can be seen by analysing the Langevin equation (3) which reads, in this case,

 dxdt=−αsgn(x)+v0σ(t). (34)

Imagine that the particle starts at some position with positive velocity . If this velocity stays at , and , the particle gets driven towards the origin with . Thus, in a finite time, it will arrive to the origin and subsequently it can not escape the origin and eventually collapses to the origin, giving rise to a delta function at . If the initial velocity is negative , it arrives at the origin even faster and, consequently, the same collapse phenomenon to the origin occurs.

The case . This case is very similar to the marginal case with , as discussed above. Indeed, in the large time limit, the stationary solution is given by a delta-function at the origin , for all and arbitrary . To see this, we again consider the Langevin equation that now reads

 dxdt=−α psgn(x)|x|p−1+v0σ(t). (35)

We now consider a similar argument as in the marginal case. We first assume that the particle starts at with . In this case, irrespective of the sign of , the particle feels a force toward the origin, which gets stronger and stronger as . Thus, in this case, the particle eventually collapses to . If the initial position and if the initial noise is positive the particle initially feels a force away from the origin. However, when the noise changes sign, it reverses the direction. Since the walk, at long times, is expected to be recurrent, at some time, it will definitely arrive at , and with probability one it will cross to the other side. Once , then again the particle gets strongly attracted to the origin, leading to an eventual collapse to the origin. A similar argument can be made for . Thus, in all cases, we would expect that, for , no matter what the starting position is, the stationary solution is just . We have verified this numerically for different values of and different values of . In Fig. 8 we show the stationary for and for different values of (setting , , in which case since ).

Thus all the cases , and can be summarised succinctly in the phase diagram in the plane as shown in Fig. 2.

## Iii Relaxation to the steady state

After having discussed the stationary state in the previous section, it is natural to ask next how the time dependent probability distribution relaxes to this stationary state at long times. In principle, this relaxation dynamics can be studied for arbitrary confining potentials . However, for simplicity we focus here on the simple harmonic case with , i.e., for which Eqs. (5) and (6) can be solved explicitly. These equations are valid on the whole real line and satisfy the following boundary and initial conditions

 P±(|x|→∞,t)=0,   and   P±(x,t=0)=12δ(x). (36)

Note that this choice of initial condition ensures that the particle, initially at the origin at , always remains confined in the finite box . This can be easily seen by solving (3) for the extreme cases, or for , which yields the deterministic solution . So the distributions have a finite support for all time .

To solve the time-dependent equations Eqs. (5) and (6), with , it is convenient to first take their Laplace transforms with respect to time. We thus define

 ~P±(x,s)=∫∞0e−stP±(x,t)dt. (37)

Taking Laplace transforms of Eqs. (5) and (6), integrating by parts the left hand side using the initial conditions in Eq. (36), we get

 (μx−v0)∂~P+∂x+(μ−γ−s)~P++γ~P−=−12δ(x),(μx+v0)∂~P−∂x+(μ−γ−s)~P−+γ~P+=−12δ(x), (38)

for . It is convenient to first make a change of variable

 x=v0μ(1−2z), (39)

that transforms the domain to . We get,

 [z∂∂z+(1−¯γ−¯s)]~P+=−¯γ~P−−14δ(z−12),[(z−1)∂∂z+(1−¯γ−¯s)]~P−=−¯γ~P+−14δ(z−12), (40)

where and

 ¯γ=γ/μ,  and  ¯s=s/μ. (41)

Let us remark that the original normalisation condition

 ∫v0/μ−v0/μ[P+(x,t)+P−(x,t)]dx=1 (42)

translates, in terms of the coordinate defined in Eq. (39)

 ∫10[P+(z,t)+P−(z,t)]dz=μ2v0. (43)

Taking the Laplace transform with respect to yields

 ∫10[~P+(z,s)+~P−(z,s)]dz=μ2v0s=12v0¯s, (44)

with . We further choose, for convenience and without any loss of generality,

 v0=1.

We now solve these equations separately for and . Due to presence of the delta-function on the right hand side of Eq. (40), the solutions in these two disjoint regions are related via the matching condition at .

 ~P+(z=1/2+ϵ,¯s)−~P+(z=1/2−ϵ,¯s)=1/4,~P−(z=1/2+ϵ,¯s)−~P−(z=1/2−ϵ,¯s)=−1/4, (45)

where . From Eqs. (40) it is clear that the solution has the following symmetry property

 ~P±(z,¯s)=~P∓(1−z,¯s), (46)

and hence it is sufficient to solve Eqs. (40) in one region only, say . Using this symmetry (46), we can further replace by in the first line of Eq. (45), which then reads

 ~P−(z=1/2−ϵ,¯s)−~P+(z=1/2−ϵ,¯s)=1/4. (47)

We now focus on the case with the boundary condition as in Eq. (47). From the two first order differential equations in (40) one can eliminate (respectively ) to get a second order differential equation for (respectively )

 z(1−z)∂2~P±∂z2+[c±−(1+a+b)z]∂~P±∂z−ab~P±=0, (48) where, a=1−¯s,  b=1−2¯γ−¯s, (49) c+=2−¯γ−¯s,  and  c−=1−¯γ−¯s. (50)

Equation (48) is a hypergeometric differential equation which has two linearly independent solutions. The general solutions for are thus readily written as

 ~P+(z,¯s)= AF(1−¯s,1−2¯γ−¯s,2−¯γ−¯s,z) + B z¯γ+¯s−1 F(¯γ,−¯γ,¯γ+¯s,z), (51) ~P−(z,¯s)= C F(1−¯s,1−2¯γ−¯s,1−¯γ−¯s,z) + D z¯γ+¯s F(1+¯γ,1−¯γ,1+¯γ+¯s,z), (52)

where is a standard hypergeometric function and and are constants yet to be fixed. We first note that the constants and are related to and . This is because and satisfy coupled first-order differential equations. Indeed, by plugging these general solutions in the two first-order differential equations (40), for . It is then easy to see that

 C=−1−¯γ−¯s¯γ A   and  D=γ¯γ+¯s