Self-propelled particle in a nonconvex external potential: Persistent limit in one dimension

TL;DR

This paper develops a method to predict the density profile of a one-dimensional active particle in arbitrary nonconvex potentials in the persistent limit, highlighting the influence of inflection points and nonlocal effects.

Contribution

It introduces a novel approach to analyze active particles in nonconvex potentials under high persistence, extending beyond equilibrium mapping techniques.

Findings

Predicts density profiles in nonconvex potentials in the persistent limit.

Highlights the importance of inflection points in the potential.

Shows nonlocal dependence of density on potential features.

Abstract

Equilibrium mapping techniques for nonaligning self-propelled particles have made it possible to predict the density profile of an active ideal gas in a wide variety of external potentials, however they fail when the self-propulsion is very persistent and the potential is nonconvex, which is precisely when the most uniquely active phenomena occur. Here we show how to predict the density profile of a 1D active Ornstein-Uhlenbeck particle in an arbitrary external potential in the persistent limit and discuss the consequences of the potential's nonconvexity on the structure of the solution, including the central role of the potential's inflection points and the nonlocal dependence of the density profile on the potential.

Equations

$$\dot{x}$$

$$\dot{f}$$

$$f(t)=f(0)+{\int }_0^{t}d{t}^{\prime }\frac{\sqrt{D}}{\tau }\xi (t^{\prime })e^{(t^{\prime }-t)/\tau }$$

$${\partial }_{t}p={\partial }_f\left(\frac{1}{\tau }fp+\frac{D}{{\tau }^2}{\partial }_{f}p\right)$$

$$p(f)=\sqrt{\frac{\tau }{2\pi D}}\exp \left[-\frac{\tau f^2}{2D}\right]$$

$$\rho (x)=\frac{df}{dx}p(f)=\frac{U^{\prime \prime }(x)\sqrt{\tau }}{\sqrt{2\pi D}}\exp \left[-\frac{\tau U^{\prime }(x{)}^2}{2D}\right]$$

$$\displaystyle\rho(x)=\dfrac{1}{Z}\,\Big{[}{\color[rgb]{0,0,1}1}+{\color[rgb]{0,190,0}\tau U^{\prime\prime}(x)}\Big{]}\,\exp\left[{\color[rgb]{0,0,1}-\dfrac{U(x)}{D}}{\color[rgb]{0,190,0}-\dfrac{\tau U^{\prime}(x)^{2}}{2D}}\right]$$

$${\partial }_t{p}_i={\partial }_f\left(\frac{1}{\tau }f{p}_i+\frac{D}{{\tau }^2}{\partial }_f{p}_i\right)+\sum _k{ϵ}_{ij}{J}_k\delta (f-f_k)$$

$$J_k=-\frac{1}{\tau }f{p}_{i(k)}(f_k)-\frac{D}{{\tau }^2}{\partial }_f{p}_{i(k)}(f_k)$$

$$\frac{1}{\tau }f{p}_i+\frac{D}{{\tau }^2}{\partial }_f{p}_i=-\sum _k{ϵ}_{ik}{J}_k\Theta (f-f_k)+a_i$$

$$p_i(f)=\left(b_i-\sum _{f_k<f}{ϵ}_{ik}\frac{{\tau }^2{J}_k}{D}{\int }_{f_k}^f\frac{du}{\bar{p}(u)}\right)\bar{p}(f)$$

$$\rho (x)=\{\begin{array}{ll}0& \text{if }{U}^{\prime \prime }(x)<0\\ U^{\prime \prime }(x)p_{i(x)}[U^{\prime }(x)]& \text{if }{U}^{\prime \prime }(x)>0\end{array}$$

$$\frac{p_1(f)}{\bar{p}(f)}$$

$$\frac{p_2(f)}{\bar{p}(f)}$$

$$-J_C+J_B$$

$$J_C-J_B$$

$$b_1$$

$$b_1-J_C{\int }_{U_C^{\prime }}^{U_B^{\prime }}\frac{{\tau }^{2}du}{D\bar{p}(u)}$$

$$b_2$$

$$b_2+J_C{\int }_{U_C^{\prime }}^{U_B^{\prime }}\frac{{\tau }^{2}du}{D\bar{p}(u)}$$

$$\frac{p_1(f)}{\bar{p}(f)}$$

$$\frac{p_2(f)}{\bar{p}(f)}$$

$$\frac{\rho (x)}{\bar{\rho }(x)}$$

$$a_i$$

$$\sum _k{ϵ}_{ik}{J}_k$$

$${\delta }_{1,i}+\sum _k{ϵ}_{ik}{J}_k{\int }_{f_k}^{f_i^R}\frac{\tau du}{D\bar{p}(u)}$$

$$\begin{gathered} \sum _{i,k}{\int }_{f_k}^{f_i^R}\frac{\tau du}{D\bar{p}(u)}\sum _k{ϵ}_{ik}{J}_k= \\ \sum _i{\int }_0^{f_i^R}\frac{\tau du}{D\bar{p}(u)}\sum _k{ϵ}_{ik}{J}_k \\ -\sum _k{\int }_0^{f_k}\frac{\tau du}{D\bar{p}(u)}\sum _i{ϵ}_{ik}{J}_k \end{gathered}$$

$$\begin{gathered} \rho (x)=\bar{\rho }(x)\times \Theta [U^{\prime \prime }(x)] \\ \times \left({\delta }_{i(x),1}-\sum _{f_k<f}{ϵ}_{i(x),k}\frac{{\tau }^2{J}_k}{D}{\int }_{f_k}^f\frac{du}{\bar{p}(u)}\right). \end{gathered}$$

$$\sum _k{ϵ}_{ik}{J}_k$$

$${\delta }_{1,i}+\sum _k{ϵ}_{ik}{J}_k{\int }_{f_k}^{{\max }_k|{ϵ}_{ik}|f_k}\frac{\tau du}{D\bar{p}(u)}$$

Full text

I Introduction

Active matter is a class of driven nonequilibrium systems in which the driving forces are controlled locally rather than globally. The development of a unified theoretical framework for active matter has transformed the way we understand a wide range of biological systems (e.g., the cell’s cytoskeleton, bacterial colonies, animal flocks) and allowed the emergence of new types of biomimetic materials (e.g., self-propelled colloids, vibrated grains, motor-filament suspensions) Marchetti et al. (2013). On the other hand, the nonequilibrium nature of active systems makes it extremely challenging to provide exact answers to even the most basic questions.

Active systems often exhibit strong, counter-intuitive responses to confinement including rectification, curvature-dependent attraction towards walls, spontaneous motion and deformation of passive objects immersed in or containing an active fluid, lack of an equation of state for pressure, dynamic geometric frustration, and more Galajda et al. (2007); Di Leonardo et al. (2010); Sokolov et al. (2010); Kaiser et al. (2014); Guidobaldi et al. (2014); Mallory et al. (2014); Fily et al. (2014, 2015, 2016, 2017); Mallory et al. (2015); Nikola et al. (2016); Yan and Brady (2015); Takatori and Brady (2016); Solon et al. (2015); Fily et al. (2018); Keber et al. (2014); Sknepnek and Henkes (2015). Yet, the density profile of a single active particle in an external potential in unknown in the general case, even in one dimension. In some cases the problem can be mapped onto an equilibrium problem and exact or good approximate results exist Tailleur and Cates (2009); Nash et al. (2010); Maggi et al. (2015). Still, some of active systems’ most striking behaviors occur precisely in the regime where this technique fails Fily et al. (2017). In this paper, we show how to derive exact analytical results for the density profile in an external potential in the persistent regime where current equilibrium mapping techniques least apply.

A number of recent advances on mapping active systems onto equilibrium ones have been driven by the Gaussian colored noise model of self-propulsion, also known as the active Ornstein-Uhlenbeck particle (AOUP) model. It allows the mapping and subsequent approximations to be done in a well controlled manner in a wide variety of systems Farage et al. (2015); Maggi et al. (2015); Marini Bettolo Marconi and Maggi (2015); Wittmann et al. (2017a, b). Its biggest caveat is that the derivation of steady-state quantities such as the density profile relies on the positive-definiteness of the operator $(1+\tau\mu\mathcal{H})$ where $\tau$ is the persistence time, i.e., the correlation time of the self-propulsion force, $\mu$ is the mobility, and $\mathcal{H}$ is the Hessian of the potential. If the violations of this condition are mild, or if they only occur in regions the particles hardly visit, it is possible to mend the theory and make predictions based on a modified equilibrium mapping Fily et al. (2017); Wittmann et al. (2017a). Conversely, strongly negative eigenvalues of $(1+\tau\mu\mathcal{H})$ in regions the particles do explore can result in unusual properties, e.g., an ideal active gas whose density profile depends nonlocally on the external potential Fily et al. (2017). Interestingly, Ref. Fily et al. (2017) derives the exact density profile in a situation in which $(1+\tau\mu\mathcal{H})$ has infinitely negative eigenvalues in locations that matter, about as far from the normal range of validity of the equilibrium mapping approach as it gets. It does so by restricting itself to the infinite persistence limit and using a quasistatic approximation wherein the self-propulsion force and the force exerted by the wall always cancel, effectively slaving the spatial dynamics to that of the self-propulsion force Fily et al. (2014, 2015, 2016).

In this paper, we show that a similar quasistatic approximation can be used to predict the steady-state density profile of a self-propelled particle in an arbitrary external potential in one dimension in the persistent limit. The predicted density profile in turn exhibits some of the same features previously predicted under nonconvex hard-wall confinement Fily et al. (2015, 2017). We first describe the AOUP model and the meaning of the persistent ($\tau\rightarrow\infty$) and quasi-thermal ($\tau\rightarrow 0$) limits (section II). We then discuss the quasistatic approximation and its generalization from hard-wall confinement Fily et al. (2015, 2017) to arbitrary external potentials (section III). In convex potentials, we show that the quasistatic approach predicts a density profile obtained by stripping the equilibrium mapping prediction of every term that survives in the quasi-thermal limit (section IV). In nonconvex potentials, we explicitly derive the density profile in a double well potential (equation (17)), then describe an algorithm to derive the density profile in an arbitrary confining potential (section V). Finally in section VI we discuss three key features of the density profile in nonconvex external potentials: the emptiness of concave regions, the crucial importance of the potential’s inflection points, and the nonlocal dependence of the profile on the potential.

II Model

We consider an active Ornstein-Uhlenbeck particles (AOUP) with position $x$ and self-propulsion force $f$ in an external potential $U$ in one dimension. The equations of motion are

[TABLE]

Dots and primes denote time and space derivatives, respectively. The first equation describes an overdamped particle with mobility $\mu$ and self-propulsion force $f$ in an external potential $U$. The second equation describes an Ornstein Uhlenbeck process, i.e., a random walk in a harmonic potential. $\tau$ is the persistence time. It is the typical time it takes $f$ to change significantly. $D$ is the active diffusion constant. The long-time ($t\gg\tau$) behavior of a free particle ($U=0$) is diffusive with diffusion constant $D$. $\xi$ is a white Gaussian noise with zero mean and variance $\langle\xi(t)\xi(t^{\prime})\rangle=2\delta(t-t^{\prime})$ where $\delta$ is the Dirac delta function.

Three additional quantities, derived from the ones above, are particularly useful to our analysis: the AOUP’s root mean squared self-propulsion force $f_{0}=\sqrt{D/(\tau\mu^{2})}$, its root mean squared velocity in the absence of external force $v_{0}=\mu f_{0}$, also called active velocity, and its persistence length $\ell=v_{0}\tau$.

For convenience we work in a unit system in which $\mu=1$ and omit $\mu$ in the rest of the paper. This effectively makes velocities and forces interchangeable.

Single active particle vs. ideal active gas. The statistical properties of a collection of noninteracting active particles, or ideal active gas, follow straightforwardly from those of a single particle. Typically, the dynamics is better understood in terms of a single particle whereas the statistics, including the density profile, are better understood in terms of a collection of particles. In this paper we use the two points of view interchangeably, choosing whichever we believe makes a specific point easier to understand. It should be kept in mind that, when referring to multiple particles, we always mean noninteracting particles. Furthermore, all normalization factors are meant for a single particle.

Thermal limit, persistent limit, and relevant variables. Only two of $D$, $\tau$, $f_{0}$ are needed to characterize the active noise. Which two are most relevant depends on the regime one is interested in.

When the persistence time $\tau$ is negligible, the active noise $f(t)$ is equivalent to a thermal noise with temperature $T=D/(\mu k)$ where $k$ is the Boltzmann constant. In this case $D$ is the only relevant variable and $f_{0}=\sqrt{D/\tau}$ is infinite. If $\tau$ is small but not negligible, $D$ and $\tau$ are the right variables for a perturbative expansion: $D$ defines the thermal solution to expand about while $\tau$ is the expansion parameter.

The main purpose of this paper is to derive the density profile in the opposite limit, or persistent limit, when $\tau$ is large compared to every other time scale in the system. We choose to derive this density profile in terms of $D$ and $\tau$ in order to 1) remain consistent with equations (1)-(2), which are our starting point, and 2) simplify the comparison with the thermal limit. However, it should be noted that the final result (equations (6) and (22)-(24)) only depends on $D/\tau$, i.e., on $f_{0}$. In particular, the persistent limit should be taken at constant $f_{0}$, with both $\tau$ and $D=\tau f_{0}^{2}$ going to infinity.

III Quasistatic dynamics

A key property of the AOUP model is that the self-propulsion force is entirely decoupled from the environment. As a result, one can integrate the equation of motion for $f$ independently of the position $x$:

[TABLE]

The other key ingredient of our theory is an approximation: in the persistent limit ($\tau\rightarrow\infty$), the self-propulsion force $f$ varies quasistatically, i.e., we can solve equation (1) as if $f$ were constant and discard the transient part of the solution on the basis that its duration is negligible compared to the time scale over which $f$ varies.

Setting $f$ constant in equation (1) yields a first order ordinary differential equation for $x$ with no explicit time dependence. A particle initially at $x_{0}$ then moves in the direction indicated by the sign of $f-U^{\prime}(x_{0})$ until it reaches the next stable fixed point $x^{*}$, namely $x^{*}=\min\{x>x_{0}{\ |\ }U^{\prime}(x)=f\text{ and }U^{\prime\prime}(x)>0\}$ if $f>U^{\prime}(x_{0})$ or $x^{*}=\max\{x<x_{0}{\ |\ }U^{\prime}(x)=f\text{ and }U^{\prime\prime}(x)>0\}$ if $f<U^{\prime}(x_{0})$. In steady-state ($\dot{x}=0$) the particle simply remains at that fixed point. Graphically, the particle moves along the $U^{\prime}(x)$ curve shown in the bottom left panel of figure 1 towards the horizontal line $U^{\prime}=f$. If it is below the line, it moves to the right. If it is above the line, it moves to the left. If it reaches an intersection point, it stops.

The quasistatic approximation assumes permanent steady-state, i.e., $U^{\prime}(x)=f$ and $U^{\prime\prime}(x)>0$ are satisfied at all times. When the self-propulsion force changes, the particle jumps instantaneously from the fixed point corresponding to its old self-propulsion force to the one corresponding to its new self-propulsion force. If the potential is convex, the location is entirely determined by the current self-propulsion force through $x=(U^{\prime})^{-1}(f)$. Conversely, if the potential is locally concave 111We do not consider potentials that are concave everywhere. They have no stable fixed point, causing particles to run away to $x=\pm\infty$., $U^{\prime}$ is not monotonic and there are multiple locations with the same external force $-U^{\prime}$, thus multiple possible fixed points for some values of $f$. Which of those fixed points the particle jumps to then depends on its previous location.

The bottom left panel of figure 1 illustrates this situation in the case of a double well potential. $B$ and $C$ are the inflexion points of the potential. Each marks the boundary between the central concave region and one of the two convex regions. A particle with $U^{\prime}_{C}<f<U^{\prime}_{B}$ has two possible stable fixed points, one between $A$ and $B$ and one between $C$ and $D$. A particle previously left of $B$ jumps to the stable point between $A$ and $B$. A particle previously right of $C$ jumps to the stable point between $C$ and $D$. In other words, the particle stays in the convex region it was previously in. When reaching the end of a convex region, the particle jumps to a new location with the same $U^{\prime}$ in a different convex region. A particle initially at $B$ whose self-propulsion force increases just past $U^{\prime}_{B}$ jumps to $D$. A particle at $C$ whose self-propulsion force decrease just below $U^{\prime}_{C}$ jumps to $A$.

The bottom right panel of figure 1 shows an alternate representation of the problem that will prove most useful to derive the density profile. It is motivated by the observation that the dynamics of $f$ is much simpler to study than that of $x$. Most notably, it does not depend on the potential. In convex potentials, the one-to-one mapping $f=U^{\prime}(x)$ means one can solve in $f$-space without worrying about $x$, then port the solution to $x$-space at the very end. A similar approach is possible in nonconvex potentials. However, for the mapping to remain one-to-one we must keep track of which convex region the particle is in as well as its self-propulsion force. In other words, we must solve for $f$ (or its probability distribution) in the space represented in the bottom right panel of figure 1. Each vertical line corresponds to a convex region. The horizontal arrows correspond to the jumps the particle experiences when it reaches the end of a region. The concave region is not represented because it contains no stable fixed point.

Finally, the dotted lines connecting the three panels of figure 1 highlight the way nonconvexity implies overlaps between the convex regions and jumps from the end of one region to the interior of another.

IV Density profile in a convex potential

To obtain the density profile $\rho(x,t)$, we first solve for the density profile $p(f,t)$ in the space of the self-propulsion force, i.e., the probability distribution of $f$, then use the one-to-one mapping $f=U^{\prime}(x)$, whose bijectivity is guaranteed by the convexity of the potential, to get $\rho(x,t)$. Applying standard Itō calculus to equation (2) yields the usual Ornstein-Uhlenbeck Fokker-Planck equation:

[TABLE]

The normalized steady-state solution is Gaussian:

[TABLE]

A simple change of variable then yields the steady-state density profile:

[TABLE]

In the rest of the paper we refer to those convex potential solutions as $\bar{p}(f)$ and $\bar{\rho}(x)$.

*Comparison with UCNA. * Before moving on to nonconvex potentials, it is instructive to compare this result with the equilibrium mapping prediction from the unified colored noise approximation (UCNA) Maggi et al. (2015) 222One may also use the Fox approximation, which yields the same steady-state profile Wittmann et al. (2017a).:

[TABLE]

where $Z$ is a normalization constant. This expression is exact in both the thermal limit ($\tau\rightarrow 0$) and the persistent (or quasistatic) limit ($\tau\rightarrow\infty$). In the termal limit, it reduces to the usual Boltzmann weight: $\rho(x)=\exp[-U(x)/kT]/Z$ where $k$ is the Boltzmann constant and $T=D/(k\mu)$ is the effective temperature. In the persistent limit, equation (7) reduces to equation (6), confirming the exactness of our quasistatic approach when $\tau\rightarrow\infty$. Interestingly, every term in (6) comes either from the thermal limit (blue terms) or the quasistatic limit (green terms). This decomposition makes sense if one recalls that the UCNA approximation is designed to be exact in both limits, and merely hopes that the intermediate $\tau$ regime can be inferred by extrapolating between them. What it does hightlight, though, is the need for a persistent limit result in nonconvex potentials, even if the gap between the thermal and persistent regimes cannot be bridged yet.

V Density profile in a nonconvex potential

V.1 Dynamics

In a nonconvex potential, the mapping $f=U^{\prime}(x)$ is not one-to-one anymore. We cannot simply derive $\rho(x)$ from $p(f)$ by change of variable. On the other hand, the mapping remains one-to-one within each convex region, and concave regions are empty because they host no stable fixed point. The solution, therefore, is to solve for the distribution $p_{i}(f)$ in each convex region $i$. The dynamics of $p_{i}(f)$ has two contributions: changes in the self-propulsion force of the particles within the convex region, and particles jumping from one region to another as illustrated in figure 1. The first contribution alone leads to equation (2). The second contribution takes the form of a point sink at each edge of each region (i.e., at every extremum of $U^{\prime}$), where particles jump off, and a point source where the jump lands. Let $k$ index the jumps. Let $f_{k}$ be the value of $U^{\prime}$ at extremum $k$. Let $J_{k}$ be the number of particles jumping from extremum $k$ per unit time. Let $\epsilon_{ik}=-1$ if jump $k$ starts in region $i$, $+1$ if jump $k$ lands in region $i$, and [math] otherwise. Then:

[TABLE]

Note that the probability density of $f$ regardless of the region, $p(f)=\sum_{i}p_{i}(f)$, still obeys equation (2) because each jump appears twice with opposite signs, once as a sink in the region where the jump originates, and once as a source in the region where it lands (same $f_{k}$, same $J_{k}$, opposite $\epsilon_{ik}$). In other words, $p(f)$ is equal to the convex potential solution $\bar{p}(f)$ from equation (5).

Since every particle crossing the edge of a convex region jumps to a different convex region, the jump rate $J_{k}$ is equal to the probability flux through the edge from the interior of the region:

[TABLE]

where $i(k)$ is the region where the jump originates. Substituting equation (9) into equation (8) yields a closed set of coupled partial differential equations for the distributions $p_{i}$. This is a convenient form to, e.g., integrate the dynamics numerically.

For the purpose of deriving the steady-state density profile, however, it is more convenient to disregard equation (9) and treat the $J_{k}$’s as unknowns for a little longer.

V.2 Steady-state

To get the steady-state solution we set $\partial_{t}p_{i}=0$ in equation (8) and integrate once with respect to $f$:

[TABLE]

where $\Theta$ is the Heaviside function and $a_{i}$ is an integration constant. Next we rewrite the left-hand side as $\dfrac{D\bar{p}}{\tau^{2}}\partial_{f}\left(\dfrac{p_{i}}{\bar{p}}\right)$ where $\bar{p}(f)$ is the convex potential solution from equation (5), multiply by $\tau^{2}/(D\bar{p})$, and integrate with respect to $f$ to get

[TABLE]

where $a_{i}$ has been set to [math] to ensure $\lim_{f\rightarrow\infty}p_{i}<\infty$ and $b_{i}$ is a new integration constant. Assuming the potential $U$ and the self-propulsion parameters $\tau$ and $D$ are known, the only unknowns left in equation (11) are the $b_{i}$’s and the $J_{k}$’s, which we must now solve for. Since equation (11) is linear in the unknowns, we need as many linearly independent equations as there are unknowns.

A first set of equations comes from imposing either that $p_{i}(f)$ can’t diverge at $+\infty$ or that the total jump rate of any region has to be zero for the region’s population to remain constant. Both yield $\sum_{k}\epsilon_{ik}J_{k}=0$.

Another set of equations comes from the boundary conditions at each end of each region. The most common boundary condition is $p_{i}=0$ at the boundary between a convex region and a convex region, which follows from imposing finite jump rates. However, a confining potential also has two convex regions that extend to $f=\pm\infty$, one on each side. Such a region has a range of $f$ in which it is the only available region. In that case the boundary condition is $p_{i}(f)=\bar{p}(f)$, which can be written at any $f$ in that range. Some types of potentials (e.g., nonconfining or noncontinuous) may require other types of boundary conditions, however we do not discuss them in this paper.

Those two sets of equations are typically sufficient to solve for the $b_{i}$’s and the $J_{k}$’s, thus for the self-propulsion force distribution in each convex region through equation (11). The mapping to position-space $f=U^{\prime}(x)$ is one-to-one in each convex region, therefore $\rho(x)=\frac{df}{dx}p_{i}(f)=U^{\prime\prime}(x)\,p_{i}[U^{\prime}(x)]$ for $x$ in region $i$. Finally, the full density profile is

[TABLE]

where $i(x)$ is the index of the convex region $x$ belongs to (when $U^{\prime\prime}(x)>0$) and $p_{i(x)}$ is the density in $f$-space in that region.

V.3 Double well potential

In the case of a double well potential like the one shown in figure 1, there are only two convex regions and two jumps. Region $1$ extends from $x=-\infty$ to $x=x_{B}$, or, in $f$-space, from $f=-\infty$ to $f=U^{\prime}_{B}$. Region $2$ extends from $x=x_{C}$ to $x=\infty$ and from $f=U^{\prime}_{C}$ to $f=\infty$. One jump goes from $B$ to $D$ with current $J_{B}$ and force $f=U^{\prime}_{B}=U^{\prime}_{D}$. The other goes from $C$ to $A$ with current $J_{C}$ and force $f=U^{\prime}_{C}=U^{\prime}_{A}$. Equation (11) reads

[TABLE]

The total jump rate equations are $-J_{C}+J_{B}=0$ (region $1$) and $J_{C}-J_{B}=0$ (region $2$).

If $f<U^{\prime}_{A}$, the only available stable fixed point is in region $1$. As a result, every particle with $f$ in that range must be in region $1$: $\forall f<U^{\prime}_{A},\ p_{1}(f)=\bar{p}(f)$. In particular, $p_{1}(U^{\prime}_{A})=\bar{p}(U^{\prime}_{A})$, which we use as the “left end” boundary condition for region $1$. Conversely, a particle with $f<U^{\prime}_{A}$ cannot be in region $2$: $\forall f<U^{\prime}_{A},\ p_{2}(f)=0$. In particular $p_{2}(U^{\prime}_{C})=0$, which is our left end boundary condition for region $2$. Similarly, $\forall f>U^{\prime}_{D},\ p_{1}(f)=0$ and $p_{2}(f)=\bar{p}(f)$, which yields the right end boundary conditions $p_{1}(U^{\prime}_{B})=0$ and $p_{2}(U^{\prime}_{D})=\bar{p}(U^{\prime}_{D})$.

After dividing the boundary condition equations by $p(f)$, the complete system (jump rates and boundary conditions) reads

[TABLE]

There are four unknowns and six equations but only four linearly independent equations. The (unique) solution is $J_{B}=J_{C}=\left(\int_{U^{\prime}_{C}}^{U^{\prime}_{B}}\frac{\tau^{2}du}{D\bar{p}(u)}\right)^{-1}$, $b_{1}=1$, $b_{2}=0$. The corresponding density in $f$-space is

[TABLE]

where $\bar{p}$ is given by equation (5). The density profile in position space follows by applying the change of variable $f=U^{\prime}(x)$ in each region separately:

[TABLE]

where $\bar{\rho}$ is given by equation (6).

V.4 Arbitrary confining potential

Consider a potential whose derivative is continuous and goes from $-\infty$ to $+\infty$. The potential is confining in the sense that, for any probability $P$, no matter how close to $1$, there is a finite region of space that contains the particles with probability $P$. More importantly for the density profile, there is a stable fixed point at any self-propulsion force $f$.

The derivation of the steady-state density profile follows the same steps as in the double well case above. There is now an arbitrary number $N$ of convex regions. The $(i,f)$ representation in the top left panel of figure 2 has $N$ branches. Let region $1$ be the (only) one that extends to $f=-\infty$ and region $N$ the (only) one that extends to $f=+\infty$. The left boundary condition for region $1$ yields $b_{1}=1$ the same way it did for the double well’s region $1$. The left boundary condition for every other region yields $a_{i}=0$ the same way it did for the double well’s region $2$. Similarly, the right boundary conditions are $p_{i}=0$ for $i<N$ and $p_{N}=p$. They are all evaluated at the region’s rightmost jump $f_{i}^{R}\equiv\max_{k}|\epsilon_{ik}|f_{k}$, which coincides with either the end of the region ($i<N$) or the point where region $N$ no longer overlaps with any other region. After including the total jump rate equations we get:

[TABLE]

$\forall i\in[1,N]$ where $\delta$ is the Kronecker delta.

Equation (18) is already a solution for the $a_{i}$’s. Equations (19) and (20) each provide $N$ linear equations for the $J_{k}$’s. Since $\sum_{i}\epsilon_{ik}J_{k}=0$ (the source and sink corresponding to the same jump cancel), summing equation (19) over $i$ yields zero and (19) provides at most $N-1$ linearly independent equations. Similarly, summing equation (20) over $i$ yields

[TABLE]

The first inner sum on the right-hand side is (19). The second inner sum is zero because $\sum_{k}\epsilon_{ik}J_{k}=0$. Thus equations (19)-(20) provide at most $2N-2$ linearly independent equations. This is also the number of jumps, thus of jump rates: one from each of regions $1$ and $N$, and two from each of the other regions. Numerically constructing the system (19)-(20) for some specific potentials suggests that eliminating one equation from (19) and one equation from (20) does typically yield a nonsingular system that can be inverted to obtain the jump rates, then the density profile by equation (12). In other words, this is a systematic route to deriving the density profile in an arbitrary external potential.

VI Discussion

Regardless of the specifics of the confining potential, the structure of the density profile derived in section V.4 contains valuable information. Let us first summarize the result:

[TABLE]

$U(x)$ is the external potential. The $f_{k}$’s are the values of the external force at the inflexion points of $U(x)$. $\bar{p}$ and $\bar{\rho}$ are given by equations (5)-(6). $i(x)$ is the label of the convex region $x$ is located in. $\epsilon_{ik}$ is $1$ if region $i$ ends at $k$, $-1$ if another convex region spills into region $i$ at $f_{k}$, and [math] otherwise. The jump rates $J_{k}$ obey

[TABLE]

Features inherited from convex potentials. The first term on the right-hand side of (22), $\bar{\rho}(x)$, is the density profile one would observe in a potential that would be convex everywhere and match $U$ around the point of interest $x$. From this term the density profile inherits most of the properties of convex-potential density profiles. Perhaps the most striking one is that the density does not depend on $U$ at all, only on $U^{\prime}$ and $U^{\prime\prime}$. In particular, particles do not accumulate where $U$ is lowest but where $\tau U^{\prime 2}-D\log U"^{2}$ is lowest, e.g., at walls Maggi et al. (2015); Fily et al. (2017).

Emptiness of concave regions. The second term on the right-hand side of (22), $\Theta[U^{\prime\prime}(x)]$, means that concave regions of $U$ are empty. This is because the particles only get stuck in convex regions. The time for which they remain stuck is controlled by their persistence time $\tau$, whereas the time they spend passing through concave regions is controlled by their self-propulsion speed $v_{0}$ and the size of the region. In the persistent limit $\tau\rightarrow\infty$, particles spend a vanishingly small fraction of their time in concave regions which end up empty.

Importance of the potential’s inflection points. The third and last term on the right-hand side of (22) captures the complex flow of particles between the convex regions of the potential. Strikingly, this term only depends on the potential $U$ through 1) the value of its derivative $U^{\prime}$ at its inflection points (the $f_{k}$’s), 2) which convex region spills into which convex region (the $\epsilon_{ik}$’s), and 3) the rate at which particles jump between regions (the $J_{k}$’s). The $J_{k}$’s in turn only depend on the $f_{k}$’s and the $\epsilon_{ik}$’s through equations (23)-(24), and the $\epsilon_{ik}$ are largely determined by the pair of $f_{k}$’s associated to each convex region. In other words, at the end of the daya it’s all about the inflection points, more specifically the value of the external force at each inflection point.

*Nonlocality of the nonconvex solution. * In a convex potential, the dependence of the density on the external potential is local: the density at position $x$ only depends on the potential through its value and the values of its first two derivatives at $x$ (see equation (6)). In a nonconvex potential, the density profile acquires a nonlocal dependence on the potential through the $f_{k}$’s, $J_{k}$’s, and $\epsilon_{ik}$’s. As we pointed out above, this dependence ends up being mostly on the $f_{k}$’s, i.e., the value of the external force at the inflection points of the potential, however far they may be from the point of interest.

VII Conclusion

Despite their considerable successes, the inability of the equilibrium mapping techniques initiated by Farage et al. (2015); Maggi et al. (2015) to provide insight into the density profile of very persistent self-propelled particles under nonconvex confinement is a major issue as this is precisely where uniquely active phenomena such as rectification occur. This paper describes an alternative approach which, although much less general, tackles the very regime in which previous methods fail most dramatically. Specifically, we derived the density profile of an active Ornstein-Uhlenbeck active particle in an arbitrary confining potential in one dimension in the persistent limit. We found it to exhibit very different properties from previous equilibrium mapping predictions, including a vanishing density in regions where the potential is concave, a nonlocal dependence of the density on the external potential, and the disproportionate influence of the external force at the potential’s inflection points. If this approach can be extended to finite persistence times and higher dimensions, it will no doubt offer significant novel theoretical insight into the surprising physics of confined active matter.