Activated escape of a self-propelled particle from a metastable state

TL;DR

This paper derives exact escape rates for active particles from metastable states, revealing how self-propulsion and potential shape influence escape behavior and induce phase transitions.

Contribution

It provides the first exact expressions for escape rates of active particles in arbitrary dimensions, including explicit solutions for RTPs and ABPs, and uncovers novel escape phenomena.

Findings

Escape rate depends on full potential shape, not just barrier height.

Particles can prefer higher barriers under certain conditions.

Escape routes can switch discontinuously, causing dynamical phase transitions.

Abstract

We study the noise-driven escape of active Brownian particles (ABPs) and run-and-tumble particles (RTPs) from confining potentials. In the small noise limit, we provide an exact expression for the escape rate in term of a variational problem in any dimension. For RTPs in one dimension, we obtain an explicit solution, including the first sub-leading correction. In two dimensions we solve the escape from a quadratic well for both RTPs and ABPs. In contrast to the equilibrium problem we find that the escape rate depends explicitly on the full shape of the potential barrier, and not only on its height. This leads to a host of unusual behaviors. For example, when a particle is trapped between two barriers it may preferentially escape over the higher one. Moreover, as the self-propulsion speed is varied, the escape route may discontinuously switch from one barrier to the other, leading to a dynamical phase transition.

Tables

Table 1. Table 1: Time step sizes and numbers of samples of the simulations in Fig. (1) of the main text. The tumbling rate is α = 1 𝛼 1 \alpha=1 .

$v$	$\delta t$	$D=0.08$	$D=0.0675$	$D=0.058$	$D=0.051$
0.3	0.1	1117	237
0.4	0.1	19329	2000	98
0.45	0.1				396
0.5	0.005	73319	5000	2100	2000
0.6	0.005	70099	20000	20000	20000
0.7	0.002	62639	40000	40000	40000
0.8	0.002	100000	200000	40000	40000
0.9	0.001	100000	200000	200000	200000
1	0.001	100000	200000	200000	200000
1.1	0.001	100000	200000	200000	200000

Equations

$$\dot{\mathbf{x}}=-\mu \nabla V+v\mathbf{u}(\theta )+\sqrt{2D}\mathbf{\xi }(t).$$

$$P({\mathbf{x}}_2,t|{\mathbf{x}}_1,0)={\int }_{{\mathbf{x}}_1}^{{\mathbf{x}}_2}\mathcal{D}\left[\mathbf{x}(t),\theta (t)\right]e^{-\frac{1}{D}\mathcal{A}[\mathbf{x},\theta ]}\mathcal{P}[\theta (t)].$$

$$\mathcal{A}[\mathbf{x},\theta ]=\frac{1}{4}{\int }_0^t{\Vert \dot{\mathbf{x}}+\mu \nabla V(\mathbf{x})-v\mathbf{u}(\theta )\Vert }^2\mathrm{d}{t}^{\prime }.$$

$$\int \mathcal{D}[\theta (t)]e^{-\frac{1}{D}\mathcal{A}[\mathbf{x},\theta ]}\mathcal{P}[\theta (t)]\underset{D\to 0}{\asymp }{e}^{-\frac{1}{D}\mathcal{A}[\mathbf{x},\tilde{\theta }]}$$

$$\mathcal{A}[\mathbf{x},\tilde{\theta }]=\underset{\theta }{\inf }\left\{\frac{1}{4}{\int }_0^t{\Vert \dot{\mathbf{x}}+\mu \nabla V(\mathbf{x})-v\mathbf{u}(\theta )\Vert }^2\mathrm{d}{t}^{\prime }\right\}.$$

$$\mathbf{u}(\tilde{\theta })=\frac{\dot{\mathbf{x}}+\mu \nabla V(\mathbf{x})}{\Vert \dot{\mathbf{x}}+\mu \nabla V(\mathbf{x})\Vert }.$$

$$\mathcal{A}[\mathbf{x}]=\frac{1}{4}{\int }_{-\infty }^{\infty }{\left(\Vert \dot{\mathbf{x}}+\mu \nabla V(\mathbf{x})\Vert -v\right)}^2\mathrm{d}{t}^{\prime },$$

$$\{\begin{array}{ll}⟨\tau ⟩& \underset{D\to 0}{\asymp }{e}^{\frac{\varphi }{D}},\\ \varphi & =\underset{\{{\mathbf{x}}_1\in {\mathbf{C}}_1,{\mathbf{x}}_2\in {\mathbf{C}}_2\}}{\inf }\underset{\mathbf{x}(t)}{\inf }\mathcal{A}[\mathbf{x}(t)].\end{array}$$

$$\mathcal{A}[x]=\frac{1}{4}{\int }_{-\infty }^{\infty }{\left[\dot{x}+\mu V^{\prime }(x)-v\right]}^2\mathrm{d}{t}^{\prime }.$$

$$\dot{x}=\partial_{x}\big{[}\mu\,(V(x)-V(C_{1}))-v\,(x-C_{1})\big{]}\equiv\partial_{x}\varphi(x)\;,$$

$$\varphi =\mu \left[V(C_2)-V(C_1)\right]-v\left(C_2-C_1\right).$$

$$A$$

$$\mathfrak{F}{\int }_{-\infty }^{C_1}\frac{\mu V^{\prime }(y)}{v^2-{\left(\mu V^{\prime }(y)\right)}^2}\mathrm{d}y$$

$$\displaystyle\underset{x\to C_{1}}{\lim}\Big{\{}\int_{-\infty}^{x}\frac{\mu V^{\prime}(y)}{v^{2}-\left(\mu V^{\prime}(y)\right)^{2}}{\rm d}y$$

$$\varphi (v)=\mu \left[V(C_2(v))-V(C_1(v))\right]-v[C_2(v)-C_1(v)].$$

$$V(x,y)={\lambda }_m\frac{x^2}{2}+{\lambda }_M\frac{y^2}{2},$$

$$\varphi(\mathbf{x}_{\text{f}})\equiv\underset{\mathbf{x(t)}}{\inf}\big{\{}\mathcal{A}[\mathbf{x}(t)]\>\big{|}\>\mathbf{x}(-\infty)\in\mathbf{C}_{1},\mathbf{x}(\infty)=\mathbf{x}_{\text{f}}\big{\}}.$$

$$\varphi =\mu V_0{\left(1-\sqrt{\frac{v^2}{2{\mu }^2{\lambda }_m{V}_0}}\right)}^2.$$

$$\mathcal{P}[\mathbf{\xi }(t)]=\mathcal{N}{e}^{-\frac{1}{2}{\int }_0^T{\Vert \mathbf{\xi }(t)\Vert }^2\mathrm{d}t},$$

$$\frac{1}{\sqrt{2D}}\left[\dot{\mathbf{x}}(t)-\left(-\mu \nabla V(\mathbf{x}(t))+v\mathbf{u}(\theta (t))\right)\right]=\mathbf{\xi }(t),$$

$$\mathcal{P}[\mathbf{x}(t)|\theta (t)]={\mathcal{N}}^{\prime }{e}^{-\frac{1}{4D}{\int }_0^T{\Vert \dot{\mathbf{x}}(t)-\left(-\mu \nabla V(\mathbf{x}(t))+v\mathbf{u}(\theta (t))\right)\Vert }^2\mathrm{d}t},$$

$$\{\begin{array}{ll}{\partial }_t{P}_{+}=& -{\partial }_x\left[\left(v-\mu {\partial }_{x}V\right)P_{+}\right]-\frac{\alpha }{2}\left(P_{+}-P_{-}\right)+D{\partial }_x^2{P}_{+},\\ {\partial }_t{P}_{-}=& -{\partial }_x\left[\left(-v-\mu {\partial }_{x}V\right)P_{-}\right]-\frac{\alpha }{2}\left(P_{-}-P_{+}\right)+D{\partial }_x^2{P}_{-},\end{array}$$

$$\{\begin{array}{ll}{\partial }_t\rho =& -{\partial }_x\left[vm-\rho \mu {\partial }_{x}V-D{\partial }_x\rho \right],\\ {\partial }_{t}m=& -{\partial }_x\left[v\rho -v\mu {\partial }_{x}V\right]-\alpha m+D{\partial }_x^{2}m.\end{array}$$

$$j(x)=vm-\rho \mu {\partial }_{x}V-D{\partial }_x\rho ,$$

$$vm-\rho \mu {\partial }_{x}V-D{\partial }_x\rho =J.$$

$$-{\partial }_x\left[v^2\rho -\left(\rho \mu {\partial }_{x}V+D{\partial }_x\rho \right)\mu {\partial }_{x}V\right]-\alpha \left(\rho \mu {\partial }_{x}V+D{\partial }_x\rho \right)+D{\partial }_x^2\left(\rho \mu {\partial }_{x}V+D{\partial }_x\rho \right)=\left(\alpha -\mu {\partial }_x^{2}V\right)J.$$

$$\{\begin{array}{l}\rho (x)\underset{x\to -\infty }{⟶}{\rho }_0,\\ \rho (C_3)=0.\end{array}$$

$$-{\partial }_x\left[\left(v^2-{\left(\mu {\partial }_{x}V\right)}^2\right)\rho \right]-\alpha \rho \mu {\partial }_{x}V=0,$$

$$\rho (x)\underset{x\to -\infty }{⟶}{\rho }_0.$$

$$-{\partial }_x\left[\left(v^2-{\left(\mu {\partial }_{x}V\right)}^2\right)\rho \right]-\alpha \rho \mu {\partial }_{x}V=\left(\alpha -\mu {\partial }_x^{2}V\right)J,$$

$$\rho (x)=\frac{{\rho }_0{v}^2}{v^2-{\left(\mu {\partial }_{x}V\right)}^2}{e}^{-\alpha {\int }_{-\infty }^x\frac{\mu {\partial }_{y}V}{v^2-{\left(\mu {\partial }_{y}V\right)}^2}\mathrm{d}y}.$$

Full text

In equilibrium, most of our intuition regarding such events is based on Kramers seminal work kramers1940brownian on Brownian particles (see hanggi1990reaction for a review). When the thermal energy is much lower than the potential barriers, there is a time-scale separation between rapid equilibration within metastable states and rare noise-induced transitions between them, a simple physical picture which is at the root of the modern view on metastability gaveau_theory_1998 ; bovier_metastability_2015 . In this limit, the mean escape time over a potential barrier of height $\Delta V$ is given by $\left\langle\tau\right\rangle{\sim}\exp(\frac{\Delta V}{k_{\text{B}}T})$. At the exponential level, the crossing time over a potential barrier only depends on its height.

To develop a corresponding intuition for activated processes in active matter, we follow Kramers and consider the dynamics of an active particle confined in a metastable well described by a potential $V$:

[TABLE]

Here, ${\bf x}$ is the position of the particle, $v$ its self-propulsion speed, and $\mu$ its mobility. The orientation of the particle ${\bf u}(\theta)$ evolves stochastically with a persistence time $1/\alpha$. Here, $\theta$ is a generalized angle parametrizing the $d-1$ dimensional unit sphere. Finally, ${\boldsymbol{\xi}}(t)$ is a Gaussian white noise which may stem from either thermal fluctuations, in which case $D=\mu\,k_{\text{B}}T$, or from fluctuations of the activity. As we show below, the escape of such an active particle from a metastable state is very different from the equilibrium case, leading to a host of interesting phenomena. For example, direct simulations of Eq. (1) show that active particles confined between two barriers may preferentially escape over the higher one, depending on the self-propulsion $v$ (See Fig. 1).

In what follows, we provide a complete solution of the Kramers problem for active particles described by Eq. (1), in any dimension, using a path-integral formalism. In contrast to existing works on first-passage times angelani2017confined ; dhar2018run ; caprini2019active , we focus on cases in which the potential is strictly confining at $D=0$ and the barrier can only be crossed using fluctuations. We refer to such case as confining potentials. We give an explicit expression for the mean escape time in terms of a variational problem for run-and-tumble particles (RTPs) berg_chemotaxis_1972 ; schnitzer_theory_1993 and active Brownian particles (ABPs) schimansky-geier_structure_1995 , the latter being studied only in $d\geq 2$ dimensions. In one dimension, RTPs had previously been studied in the limits $\alpha\to 0$ and $\alpha\to\infty$ geiseler2016kramers ; Here, we provide the full solution of the activation time for RTPs for all $\alpha$, including its sub-exponential prefactor. In cases with multiple competing reaction paths, our results provide the selection principle for the most likely escape route. In particular, we explain the dynamical phase transition observed in Fig. 1.

For confining potentials, it is natural to divide the barrier into separate regions depending on whether the force $|\nabla V|$ is larger or smaller than the propulsion force $f_{p}=v/\mu$. Consider, for instance, the escape in one dimension from a metastable well, see Fig. 2. We can identify four different regions separated by three points $\{C_{1},C_{2},C_{3}\}$ satisfying $|V^{\prime}(C_{i})|=f_{p}$. In regions (i) and (iii), when $x\leq C_{1}$ or $C_{2}\leq x\leq C_{3}$, the particles feel a force $-V^{\prime}$ smaller in magnitude than $f_{p}$. In the $D\to 0$ limit the contribution of the noise $\xi(t)$ to the dynamics can be neglected. In region (ii), where $C_{1}\leq x\leq C_{2}$, the particles cannot climb the potential without the noise $\xi(t)$. Crossing this region is therefore a rare event which controls the escape from the metastable state. In region (iv), where $x>C_{3}$, the particles would need the noise to come back to region (i), were they to reverse direction. This is a rare event and the particle has thus effectively crossed the barrier once it has reached $C_{3}$. The generalization of these points to lines or surfaces in higher dimensions (denoted ${\bf C}_{i}$) is straightforward and an example is displayed in Fig. 2 111Note that in case of saddles the surfaces ${\bf C}_{2}$ and ${\bf C}_{3}$ may merge into a single surface with two distinct faces.. Note that the problem is activated only if region (ii) exists. Otherwise, the problem, as considered for example in 1d in angelani2014first , is a first-passage problem with no instanton physics.

The activated process only corresponds to moving across region (ii) so that the crossing probability is given, to leading order, by histories connecting points on ${\bf C}_{1}$ and ${\bf C_{2}}$. To obtain the escape time we then write the transition probability $P(\mathbf{x_{2}},t|\mathbf{x_{1}},0)$ to be at $\mathbf{x_{2}}\in{\bf C}_{2}$ at time $t$ starting at $\mathbf{x_{1}}\in{\bf C}_{1}$ as a path integral in its Onsager–Machlup form onsager_fluctuations_1953

[TABLE]

${\cal P}[\theta(t)]$ is the probability of a history of the angle $\theta$. For example, ABPs in 2d with rotational diffusivity $\alpha$ lead to ${\cal P}[\theta(t)]\propto e^{-\int_{0}^{t}\dot{\theta}^{2}/(4\alpha){\rm d}t^{\prime}}$. In Eq. (2), the action $\mathcal{A}[\mathbf{x},\theta]$ is given by

[TABLE]

We first integrate expression (2) over the paths $\theta(t)$ to obtain an effective action for the probability of a path $\mathbf{x}(t)$. In the limit $D\rightarrow 0$, we use a saddle-point approximation in (2) to get:

[TABLE]

where $\asymp$ stands for logarithmic equivalence and $\tilde{\theta}(t)$ is the path satisfying the variational problem

[TABLE]

Note that ${\cal P}[\theta(t)]$ is a subdominant contribution and any cost to the action arising from it can be ignored to leading order 222The results might change in cases where, say for ABPs, $\alpha$ is proportional to $D$.. Clearly, the optimum requires $\mathbf{u}(\theta)$ to be in the same direction as $\mathbf{\dot{x}}+\mu\nabla V(\mathbf{x})$ so that

[TABLE]

Using Eqs. (6) and (3), we find that the transition probability between ${\bf x}_{1}$ and ${\bf x}_{2}$ is dominated by paths which minimize the action

[TABLE]

where we have sent the limits of the integral to $\pm\infty$, using the fact that extremal trajectories start and end at stationary points (see for instance Tailleur2008JPA ; baek2015singularities ). Finally, the escape time is given by

[TABLE]

The inner minimization corresponds to optimizing the action over different paths; it is realized by an instanton $\mathbf{x}(t)$ which connects $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$. The outer minimization corresponds to optimizing over all possible initial and final positions of the instanton. Eq. (8) provides a full solution to the escape problem for both ABPs and RTPs as a variational problem. It generalizes the Kramers law and we discuss the physics of the quasi-potential barrier $\phi$ below. Note that when $v=0$ the minimizers of the action are $\mathbf{\dot{x}}=\mu\nabla V(\mathbf{x})$ and we recover the usual Kramers law with $\phi=\mu\,\Delta V$, where $\Delta V$ is the minimal potential difference across the barrier. We now turn to apply our results to a general one-dimensional potential barrier and to an elliptic well in two dimensions.

RTPs in one dimension: Here, ${\bf u}(\theta)$ is replaced by a binary variable $u=\pm 1$ which flips with rate $\alpha/2$. As in Fig. 2, the barrier is located on the right of the metastable well. $\mathbf{x}_{1,2}$ are then given by $C_{1,2}$. Clearly, the minimal action is obtained by particles with $u=1$: particles which reverse their motion in the middle of the instanton are exponentially less likely to cross the barrier. The action then reduces to

[TABLE]

It is thus equivalent to an equilibrium problem in an effective titled potential $\varphi(x)/\mu$; the instanton solution obeys

[TABLE]

which gives, for the quasi-potential barrier introduced in Eq. (8),

[TABLE]

Our predictions (8) and (11) are verified in Fig. 3 using direct simulation of Eq. (1) with a single barrier.

Using asymptotic techniques supp ; bouchet2016generalisation , we also obtain the leading sub-exponential amplitude of the transition time (8). For simplicity we consider a boundary condition in which the potential is flat on the left of the barrier and the density of particles in that region is $\rho_{0}$; other boundary conditions are discussed in supp . The mean time between particles crossing the barrier is then given by $\langle\tau\rangle\underset{D\rightarrow 0}{\sim}Ae^{\frac{\phi}{D}}$, where

[TABLE]

Here, ${\cal T}_{\rm inst}=\mathfrak{F}\int_{C_{1}}^{C_{2}}\frac{{\rm d}y}{\partial_{y}\varphi}$ is the duration of the instanton, $k_{i}=\mu V^{\prime\prime}(C_{i})$, $\Gamma(x)$ is the Euler Gamma function, and $\mathfrak{F}$ denotes the finite part of the integral, defined by removing the logarithmic divergences occurring at $C_{1}$ and $C_{2}$, e.g.

[TABLE]

The term $e^{-\frac{\alpha}{2}{\cal T}_{\rm inst}}$ has a simple interpretation: it is the probability that the particle does not flip along the instanton. Note that the $v=0$ limit is singular: all histories of $u(t)$ are then equally likely, a degeneracy which otherwise does not exist.

Equations (11) and (29) provide an explicit solution to the Kramers problem in one dimension. Note that the effect of the activity cannot be cast into a simple description with an effective temperature. Both $\phi$ and the prefactor indeed depend on the full functional form of the potential $V$.

Dynamical phase transition: We now show how the analysis of the quasi-potential accounts for the non-trivial choice of escape routes when the particle is trapped between two potential barriers. In the small $D$ limit, the escape time is controlled by the quasi-potential (11) of each barrier, which we can study separately. For the right barrier, the explicit dependence of $\phi$ on $v$ reads

[TABLE]

When $v=0$, we recover the standard Kramers result $\phi(0)=\mu\left(V(C_{2})-V(C_{1})\right)$. Using $\mu V^{\prime}(C_{1})=\mu V^{\prime}(C_{2})=v$, one has $\phi^{\prime}(v)=-(C_{2}-C_{1})$, which implies that $\phi$ is a decreasing function of $v$. When $v>v_{\text{cr}}\equiv\underset{x}{\max}\left\{\mu|V^{\prime}(x)|\right\}$, the particle can cross the barrier without thermal activation so that $\phi(v_{\text{cr}})=0$. $\phi(v)$ thus decreases from the equilibrium $v=0$ value to zero. The initial decrease of the escape time is given by $\phi^{\prime}(0)=-[C_{2}(0)-C_{1}(0)]\equiv-\ell$ which is nothing but the distance between the maxima and the minima of the potential $V$, i.e. the width of the barrier. The same construction holds for the second barrier.

Next, consider the two potential barriers $V_{\text{R},\text{L}}(x)$ of equal height described in Fig. 4. The right barrier is wider, $\ell_{\text{R}}>\ell_{\text{L}}$, but has a larger maximal slope than the left barrier so that $v_{\text{cr}}^{\text{R}}>v_{\text{cr}}^{\text{L}}$. To leading order, the escape rates over the two barriers for $v=0$, $\phi_{\text{L}}(0)$ and $\phi_{\text{R}}(0)$, are equal. Following the above discussion, $\phi_{\text{R}}(v)$ decreases faster than $\phi_{\text{L}}(v)$ near $v=0$ because the right barrier is wider than the left one: for small $v$, the particle is more likely to escape over the right barrier. $\phi_{\text{R}}(v)$ however vanishes at a value $v_{\text{cr}}^{\text{R}}$ larger than $v_{\text{cr}}^{\text{L}}$ due to the existence of a steeper portion in the right barrier. For large $v$, the escape is thus more likely through the left barrier. Hence, there exists a critical self-propulsion speed at which the most likely escape route changes discontinuously. The physics presented in Fig. 1 can be understood from the above discussion, the sole difference being that the escape rates are different at $v=0$ due to the different barrier heights. In the $D\to 0$ limit, the sigmoid function presented in Fig 1 hence converges to a discontinuous step function. In fact, it is straightforward to see that one could also observe not one but two successive dynamical phase transitions if the larger and steeper barrier were also higher. Interestingly, the dependence of the escape time on $v$ can be used to sort active particles depending on their velocities (See Supplementary movie).

Escape from two-dimensional elliptic potentials: We now consider the escape of active particles from a two-dimensional potential well of the form

[TABLE]

with $\lambda_{M}>\lambda_{m}$ (for an analysis of the steady-state distribution for the case $\lambda_{m}=\lambda_{M}$, see malakar2019exact ). We assume that particles escape when they reach a given height $V(x,y)=V_{0}$. This level line $\mathbf{C}$ replaces $\mathbf{C}_{2}$ of the general discussion, see Fig. 5.

The most-probable escape routes can be computed by solving the Euler-Lagrange equations for the action given in Eq. (7), as detailed in the SI. Following the previous argument we introduce

[TABLE]

which yields, at the exponential level, the probability to reach any point $\mathbf{x}_{\text{f}}$ on the boundary. This log-probability, which we compute in supp , is plotted, in Figure 5, as a function of the angular parametrization of $\mathbf{x}_{\text{f}}$, and compared with numerics. Interestingly, the quasi-potential is not constant over the boundary: the particles have a much larger probability to escape in the direction of the major axis of the ellipse. This is the most striking difference with the equilibrium problem: For passive Brownian particles, the quasi-potential is $\varphi(\mathbf{x})=\mu V(\mathbf{x})$, so that particles have an equal probability (at the exponential level) to escape through any point along the boundary ${\bf C}$. Activity thus breaks the equilibrium quasi-potential symmetry.

Furthermore, one can compute explicitly the full expression of $\phi$ given by the minimum of the function $\varphi(\mathbf{x}_{\text{f}})$ along $\mathbf{C}$:

[TABLE]

The escape time from the elliptical well is then given by $\langle\tau\rangle\asymp\exp(\phi/D)$. It solely depends on the potential height, the particle speed, and the semi-axis corresponding to the most likely exit direction. As expected, we recover the standard equilibrium result $\phi=\mu V_{0}$ when $v=0$.

By providing a full solution to the Kramers problem for both ABPs and RTPs in any dimensions, we have highlighted how the physics of these non-equilibrium systems is very different from that of the equilibrium problem. In particular, the activation barrier, encoded in the quasi-potential, is not solely defined by the height of the potential well. Instead, it corresponds to the region where the self-propelling force fails to overcome the confining one, leading to activation paths and times that depend in a non-trivial way on both the self-propelling speed and the full shape of the potential, and to a wealth of unusual features. Our results also highlight why an effective equilibrium approach is inappropriate. Beyond the case addressed here of an external potential, escape problems play an important role in a host of collective phenomena, from nucleation to glassy physics. It will thus be very interesting to see how the phenomena uncovered in this manuscript play a role in these more complicated systems.

Acknowledgments: YK & EW acknowledge support from I-CORE Program of the Planning and Budgeting Committee of the Israel Science Foundation and an Israel Science Foundation grant. JT is funded by ANR Bactterns. JT & YK acknowledge support a joint CNRS-MOST grant. VL is supported by the ERC Starting Grant No. 68075 MALIG, the ANR-18-CE30-0028-01 Grant LABS and the ANR-15-CE40-0020-03 Grant LSD.

Appendix A Path integral formulation

In this section, we give a simple derivation of the path integral formulation Eq. (2) in the main text. Let $\boldsymbol{\xi}(t)$ be the standard Gaussian white noise in $d$ dimensions with correlation function $\left\langle\xi_{i}(t^{\prime})\xi_{j}(t)\right\rangle=\delta_{ij}\delta(t-t^{\prime})$. The probability of a given realization of $\boldsymbol{\xi}(t)$ on the interval $[0,T]$, denoted $\mathcal{P}[\boldsymbol{\xi}(t)]$, can be formally written as

[TABLE]

where $\mathcal{N}$ is the normalization factor. From Eq. (1) in the main text, it appears that the probability distribution of any path $\mathcal{P}[\mathbf{x}(t)]$ can be expressed from the distribution of the noise $\mathcal{P}[\boldsymbol{\xi}(t)]$ through the change of variable

[TABLE]

where $\theta(t)$ is a given realization of the angle history. Combining (18-19), the probability of a given path $\mathbf{x}(t)$ conditioned on a realization $\theta(t)$ becomes

[TABLE]

where $\mathcal{N}^{\prime}$ is the new normalization factor. Note that in the Itô convention of stochastic calculus, the Jacobian of the change of variable is a constant independent of the field $\mathbf{x}(t)$, hence the new normalization factor $\mathcal{N}^{\prime}$. As the dynamics of $\theta(t)$ is decoupled from that of $\mathbf{x}(t)$, the joint probability of $[\mathbf{x}(t),\theta(t)]$ is given by $\mathcal{P}[\mathbf{x}(t),\theta(t)]=\mathcal{P}[\left.\mathbf{x}(t)\right|\theta(t)]\mathcal{P}[\theta(t)]$. Integrating the joint probability $\mathcal{P}[\mathbf{x}(t),\theta(t)]$ over all possible angle histories and all possible paths $\mathbf{x}(t)$ joining $\mathbf{x}_{1}$ to $\mathbf{x}_{2}$ gives the result presented in Eq. (2) and Eq. (3) in the main text.

Appendix B Run-and-Tumble particles in one-dimension

The most straightforward method to obtain the leading order behavior of the escape rate is presented in the main text. However, in order to find the sub-leading correction it is useful to employ a different approach involving asymptotic matching of solutions. In what follows this approach, whose final results is Eq. (12) of the main text, is detailed. In addition to the result of the main text we also provide here the prefactor for the mean escape time from a metastable well in Eq. (37) of Sec. B.4.

B.1 Description of the problem and main equations

We study RTPs particles in one-dimension. The particles experience a driving force $v/\mu$ which reverses its direction with rate $\alpha/2$. In addition, they are subject to an external potential $V$. Denoting by $P_{+}(x,t)$ and $P_{-}(x,t)$ the probability density of particles moving to the right and left respectively, the Fokker–Plank equation for $P_{+},P_{-}$ is

[TABLE]

Here, as in the main text, $\mu$ is the mobility and $D$ is the diffusion coefficient. We are interested in the limit $D\rightarrow 0$, which can be interpreted physically as the asymptotic regime $D\ll v\ell$ where $\ell$ is the barrier length. Let $\rho=P_{+}+P_{-}$ be the total density of active swimmers in the medium, and $m=P_{+}-P_{-}$. From Eq. (20)

[TABLE]

The first equation describes the mass conservation with a flux

[TABLE]

which is constant $j(x)=J$ in the steady-state. This gives

[TABLE]

Using this relation in the second equation of (21) we have

[TABLE]

We now solve this equation using standard asymptotic matching techniques with the boundary conditions

[TABLE]

In this configuration, the transition of particles across the barrier is a Poisson process with rate $J$. The mean waiting time between two particles crossing the barrier is given by $\left\langle\tau\right\rangle=\frac{1}{J}$. As stated above, we will also consider the situation where the particles start from a metastable state (instead of the boundary conditions described by Eq. (24)) and provide an explicit expression of the mean escape time in that case.

B.2 Methods

To proceed we solve the problem in the three regions (i), (ii), and (iii) defined in the main text and then match the solutions. We first note the following about the different regions.

**region (i): **The flux $J$ is so small compared to the other terms in Eq. (23) that the solution is given by the steady-state with $D=0$ and $J=0$. The corrections are of order $D$. Namely, we solve

[TABLE]

together with the boundary condition

[TABLE] 2. 2.

region (ii): Here we use the WKB-like Ansatz $\rho(x)=C_{D}(x)e^{-\frac{\varphi(x)}{D}}$ in Eq. (23). The expression for $\varphi(x)$ is identical to that obtained using the methods of the main text. Note that also here to leading order $J=0$. 3. 3.

region (iii): As in region (i) the contribution of diffusion terms $\propto D$ to the dynamics can be neglected. However, since the density of particles is now very low, the current $J$ is no longer negligible and one has to be accounted for it. Therefore, here we solve

[TABLE]

with the absorbing boundary condition $\rho(C_{3})=0$.

The solutions found separately in regions (i), (ii) and (iii) have to match together at the two points $C_{1}$ and $C_{2}$. To do this we have to calculate the structure of the solution near the the two points $C_{1}$ and $C_{2}$. These are given, as we detail below, by boundary layers of size $\sqrt{D}$ which can be matched to the solutions in the different regions.

B.3 Solutions

We next carry out the calculation outlined above in detail.

B.3.1 Region (i)

The explicit solution of Eq. (25) is

[TABLE]

In order to match this solution we have to understand how it behaves near $C_{1}$. To do this we make the change of variable $x\leftarrow x-C_{1}$. The force can be expanded according to

[TABLE]

The equivalent of the integral in the exponential of Eq. (27) is

[TABLE]

where $\gamma_{1}$ is a finite constant that depends explicitly on the potential through the relation

[TABLE]

where the last equality defines, as in the main text, the finite part of the diverging integral. Note that there is some arbitrariness in the definition of the finite part. Any function of the form $\frac{1}{2k_{1}}\log\left(\frac{\left|x\right|}{L}\right)$, where $L$ is some arbitrary length scale, could be removed from the integral to define the finite part. The above choice $L=\frac{v}{k_{1}}$ has been used in order to make the final expression for the mean escape time more compact.

Restoring the original coordinate $x$, we therefore find

[TABLE]

The solution has two different behaviors depending on the value of the second derivative $k_{1}=\mu\left.\partial_{y}^{2}V\right|_{x=C_{1}}$. The density diverges at the critical point $x=C_{1}$ if $k_{1}>\frac{\alpha}{2}$, and vanishes if $k_{1}<\frac{\alpha}{2}$. Since $k_{1}>0$ the diverging solution remains integrable at $C_{1}$. We comment that it is straightforward to see that

[TABLE]

This implies that only right moving particles reach $C_{1}$.

B.3.2 Region (ii)

In region (ii), we use the WKB-like Ansatz

[TABLE]

where the large deviation pre-factor function can be expanded in powers of $D$ as

[TABLE]

To leading order it is easy to check that, as expected, this reproduces the Eq. (10) of the main text for $\varphi(x)$. Using this solution with the expansion of the pre-factor we find to next order

[TABLE]

whose solution is

[TABLE]

with $x_{0}$ an arbitrary point between $C_{1}$ and $C_{2}$.

Again to match this solution we have to consider its behavior close to the two critical points $C_{1}$ and $C_{2}$. To this end, we make the change of variable $x\leftarrow x-C_{1}$ and study the behavior of $C_{D}^{0}(x)$ close to $C_{1}$. Close to $x=0$, we use the expansion

[TABLE]

which shows that the integral in (30) can be expanded around $x=0$ as

[TABLE]

where $\gamma_{2}$ is a finite constant, and we used Eq. (10) of the main text. Using the same notations as in section B.3.1, we have

[TABLE]

Coming back to the original variable $x$, this gives

[TABLE]

The same line of arguments, gives the equivalent of the pre-factor close to $C_{2}$ as

[TABLE]

where $k_{2}=\mu\partial_{x}^{2}V(C_{2})$ is the (negative) second derivative of the potential, and $\mathfrak{F}\int_{x_{0}}^{C_{2}}\frac{{\rm d}y}{\partial_{y}\varphi}$ is the finite part defined as

[TABLE]

In what follows to match this solution with the other regions we note that the above results imply that close to $C_{1}$

[TABLE]

and close to $C_{2}$

[TABLE]

These specify the boundary layers at the edges of region (ii). Their typical extension is $\sqrt{\frac{D}{k_{1}}}$ and $\sqrt{\frac{D}{\left|k_{2}\right|}}$ at $C_{1}$ and $C_{2}$ respectively.

B.3.3 Region (iii)

Eq. (26) can be solved to give

[TABLE]

Note that this expression is well defined, because $e^{\alpha\int_{x}^{y}\frac{\mu\partial_{z}V}{v^{2}-\left(\mu\partial_{z}V\right)^{2}}{\rm d}z}$ is integrable close to $C_{3}$. Using this we find that near $C_{2}$ the solution can be written as

[TABLE]

where again, the notation $\mathfrak{F}$ means

[TABLE]

B.3.4 Matching at the boundary layers

We now have to match all the solutions (28,31,32,33) at the two critical points $C_{1}$ and $C_{2}$. To do this we need to solve the Fokker-Planck equation in the boundary layers around $C_{1}$ and $C_{2}$. To this end we define the variables $x-C_{i}=\sqrt{\frac{D}{|k_{i}|}}y_{i}$. Using this in Eq. (23) we obtain to zeroth order in $D$

[TABLE]

where ${\rm sgn}(k)=\pm 1$ denotes the sign of $k_{i}$. The solutions of this equation for large positive or negative values of $y_{i}$ have to be matched with the solutions in the different region. The solution for $i=1$ satisfies

[TABLE]

and for $i=2$

[TABLE]

where $A_{D_{1}}$ and $A_{D_{2}}$ are two undetermined constant. By matching the asymptotic behavior (36) of the boundary layer solution with the behavior of the solutions (28,31,32,33) in the different regions close to $C_{1}$ and $C_{2}$ one finds after a lengthy calculations Eq. (12) of the main text.

B.4 Mean escape time from a metastable well

We now generalize our result to the mean escape time from a metastable well. We introduce the critical point $C_{0}$ on the left of $C_{1}$ such that $\mu\partial_{x}V(C_{0})=-v$. The metastable well is represented in Fig. 6.

According to expression (27), the zero-fluctuations solution in region (i) writes

[TABLE]

where $x_{b}$ is some arbitrary point between $C_{0}$ and $C_{1}$, and $N$ is a constant, given by the normalization constrain $\int_{C_{0}}^{C_{1}}\rho(x){\rm d}x=1.$ We find

[TABLE]

The mean escape time is then simply given by Eq. (12) replacing the term

[TABLE]

by

[TABLE]

which can be equivalently written as

[TABLE]

We obtain the formula

[TABLE]

Appendix C Escape from a two-dimensional elliptic potential

This section presents the computation of the quasi-potential for the active escape problem out of the two-dimensional elliptic barrier described in the main text. The potential can be written as

[TABLE]

where $A=\begin{pmatrix}\mu\lambda_{m}&0\\ 0&\mu\lambda_{M}\end{pmatrix}$ is a symmetric matrix of the second derivatives of the potential. We consider without loss of generality that $0<\lambda_{m}<\lambda_{M}$.

Using the results of the main text, the fluctuation paths between specified initial and final positions are minimizers of the action

[TABLE]

To compute the fluctuation paths, we solve the Euler-Lagrange equation. As will become clear, it is useful to consider the momentum

[TABLE]

Interestingly for a quadratic potentials we find from Eq. (38)

[TABLE]

Using this relation, the Euler-Lagrange equations then translate into an equation for the momentum $\mathbf{p}$

[TABLE]

whose solution is $\mathbf{p}(t)=e^{At}\mathbf{p}_{0}$. In the present problem, $\mathbf{p}_{0}$ should be understood as the momentum at the final position of the trajectory $\mathbf{x}(t=0)=\mathbf{x}_{\text{f}}$. The explicit solution of $\mathbf{p}(t)$ together with Eq. (39) gives the first order equation for ${\bf x}$

[TABLE]

To solve Eq. (40), we first take the norm of both sides of the equality to get

[TABLE]

with the implicit assumption that the instanton path satisfies the condition $\left\|\dot{\mathbf{x}}+A\mathbf{x}\right\|>v$. Using (41) in Eq. (40), we have

[TABLE]

The boundary conditions for this equation are

[TABLE]

Eq. (42) together with the constraints (43) can only be satisfied if $\mathbf{x}_{1}$ is an eigenvector of $A$. To see this, we expand of the right-hand side of Eq. (42) in the limit $t\rightarrow-\infty$. Because the two eigenvalues of $A$ satisfy $\lambda_{m}<\lambda_{M}$, we have $\left\|e^{At}\mathbf{p}_{0}\right\|\underset{t\rightarrow-\infty}{\sim}\ p_{0}^{x}e^{\lambda_{m}t}$ where $(p_{0}^{x},p_{0}^{y})$ are the two components of $\mathbf{p}_{0}$. We further have $e^{At}\mathbf{p}_{0}=p_{0}^{x}e^{\lambda_{m}t}\mathbf{e}_{x}+p_{0}^{y}e^{\lambda_{M}t}\mathbf{e}_{y}$. For $p_{0}^{x}\neq 0$, the first term in Eq. (42) thus gives

[TABLE]

which proves, using $\dot{\mathbf{x}}(t)\underset{t\rightarrow-\infty}{\longrightarrow}0$, that $\mathbf{x_{1}}=\left(\pm\frac{v}{\mu\lambda_{m}},0\right)$. When $p_{0}^{x}=0$ we find $\mathbf{x_{1}}={\bf x}_{1}^{y}=(0,\pm\frac{v}{\mu\lambda_{M}})$. This has a simple geometric interpretation. Generically $\mathbf{x}_{1}$, sitting on the $x$-axis, is a local extremum of $V(\mathbf{x})$ on the curve ${\bf C}_{1}$ (see Fig. (4) of the main text). With the exception of fluctuation paths which end on the $y$-axis all the paths start at one of the two local maxima located at $\mathbf{x_{1}}=\left(\pm\frac{v}{\mu\lambda_{m}},0\right)$. This clearly minimizes the cost of the path. Fluctuation paths which end on the $y$-axis start at ${\bf x}_{1}^{y}$.

We now turn to the full computation of the quasi-potential $\varphi(\mathbf{x}_{0})$, where $\mathbf{x}_{0}$ is the final position of the fluctuation path. The explicit expression of $\mathbf{x}_{0}$ can be computed from the general solution of Eq. (42)

[TABLE]

Using Eq. (41) in Eq. (38) and carrying out the integration in time we obtain the large deviation rate function as a function of $\mathbf{p}_{0}$

[TABLE]

Eqs. (44) and (45) can both be solved numerically to compute the quasi-potential $\varphi$ displayed in Fig. (4) of the main text.

Besides, it is straightforward use Eq. (44) to perform a small-fluctuations expansion around ${\bf p}_{0}$ in order to show that the action is minimal for paths moving only along the $x$-direction. Using this one can then easily compute the full expression for $\phi=\min\left\{\varphi(\mathbf{x})|V(\mathbf{x})=V_{0}\right\}$. We find

[TABLE]

As expected, we recover the standard equilibrium result $\phi=\mu V_{0}$ when $v=0$.

Appendix D Supplementary information for the figures

In this section, we provide the details about the potential $V(x)$ in each figure of the main text. We also describe the algorithm used for the numerics in Fig. 1, Fig. 3, and in the supplementary movie.

D.1 Algorithm

For all the simulations presented in the main text, we adapted the Heun algorithm to simulate, for each individual particle, the following over-damped stochastic differential equation (see Eq. (1) in the main text)

[TABLE]

Here, ${\bf x}$ is the position of the particle and $v$ is its self-propulsion speed. The orientation of the particle ${\bf u}(\theta)$ evolves stochastically with a persistence time $1/\alpha$. Note that compared to Eq. (1) of the main text, we have set $\mu=1$ everywhere. $\boldsymbol{\xi}(t)$ is a vector of Gaussian white noise, such that

[TABLE]

We discretized the time with time step $\delta t$, and updated the status of the particles according to

[TABLE]

where $W_{t,i}\sim\mathcal{N}(0,1)$ is a normal distributed random number.

The reorientation of the particle is independent from its position. We thus sample the next tumbling time of each particle from the exponential distribution $\alpha e^{-\alpha t}$. We split the time step where the tumbling happens into two smaller time steps: we first update the position of the particle until the time it tumbles, and then we uniformly randomly assign a new direction, and finish the remaining time.

We calculated the $C_{1}$, $C_{2}$, and $C_{3}$ numerically using false position method. The particles were considered to have escaped when their position reaches $C_{3}$.

D.2 Figure 1

The left potential barrier $V_{\text{L}}(x)$ in Fig. (1) is defined by

[TABLE]

where the coefficients are defined through $\Delta_{\text{L}}=2/(\ell_{\text{L}}Z_{\text{L}})$, $\ell_{\text{L}}=3$, $Z_{\text{L}}=1$ and $H_{\text{L}}=1.2$. The width of the barrier is $4\ell_{\text{L}}$, the height of the potential is $2H_{\text{L}}$ and the maximal slope is $H_{\text{L}}Z_{\text{L}}$.

The right potential barrier $V_{\text{R}}(x)$ in Fig. (1) is defined by

[TABLE]

where the coefficients are defined through $\Delta_{\text{R}}=2/(\ell_{\text{R}}Z_{\text{R}})$, $\ell_{\text{R}}=12$, $Z_{\text{R}}=2$ and $H_{\text{R}}=1$. The width of the barrier is $4\ell_{\text{R}}$, the height of the potential is $2H_{\text{R}}$ and the maximal slope is $H_{\text{R}}Z_{\text{R}}$.

The parameters for the Heun algorithm defined in section D.1 are listed in Table 1. We simulate each particles until it escapes from the barrier, that is, until it reaches $C_{3}$.

D.3 Figure 3

The potential in Fig. 3 is defined through

[TABLE]

where $A=1.5$, $B=1$, $C=2$. Those conditions correspond to a reflective boundary at $x=0$. We set $\alpha=1$, we use a time step $\delta t=0.001$, and $N_{\mathrm{samples}}=2\times 10^{6}$. The values of the particle’s velocity are given by $v=1,\ 1.5,\ 2,\ 2.5$ respectively. We simulated each particles until it reaches $C_{3}$.

D.4 Figure 4

The functional dependence of the two barriers is exactly the same as in Fig. 1 (see section D.2), with the parameters:

Left barrier: $\ell_{\text{L}}=3$, $Z_{\text{L}}=1$, $H_{\text{L}}=1$.

Right barrier: $\ell_{\text{R}}=8$, $Z_{\text{R}}=2$, $H_{\text{R}}=1$.

Contrary to the potential of Fig. (1), the left and right barriers have the same height here. The quasi-potential has thus the same value for the two barriers at $v=0$, and escape of passive particles is equally likely left and right, at least at the exponential level.

D.5 Movie

The potential $V(x)$ in the movie is built according to

[TABLE]

The two populations have $v=1.5$ (blue particles) and $v=3.3$ (red particles) respectively. Other parameters are $\alpha=1$, $D=0.06$, $\delta t=0.001$. $N_{\mathrm{samples}}=10^{5}$ particles are used to generate the histogram. The total time of the simulation is $t=20000$.

Active particles start in the metastable state around $x=0$, and escape left or right. When they escape, they are then trapped in the two deep wells located at $x=-3$ and $x=3$ respectively.