Steady State of an Active Brownian Particle in Two-Dimensional Harmonic Trap

TL;DR

This paper derives an exact series solution for the steady-state distribution of a 2D active Brownian particle in a harmonic trap, revealing a re-entrant transition between active and passive regimes validated by simulations.

Contribution

It provides the first exact series solution for the steady-state distribution of an active Brownian particle in a harmonic trap, enabling detailed analysis of different regimes.

Findings

Identification of active and passive regimes

Prediction of a re-entrant active-passive transition

Validation through numerical simulations

Abstract

We find an exact series solution for the steady-state probability distribution of a harmonically trapped active Brownian particle in two dimensions, in the presence of translational diffusion. This series solution allows us to efficiently explore the behavior of the system in different parameter regimes. Identifying "active" and "passive" regimes, we predict a surprising re-entrant active-to-passive transition with increasing trap stiffness. Our numerical simulations validate this finding. We discuss various interesting limiting cases wherein closed form expressions for the distributions can be obtained.

Figures

Equations

$$\begin{array}{rl}\frac{d\mathbf{\rho }}{d\tau }& =u_0\hat{\mathbf{e}}(\theta )-\mu k\mathbf{\rho }+\sqrt{2D_t}{\mathbf{\xi }}_r(\tau ),\\ \frac{d\theta }{d\tau }& =\sqrt{2D_{\theta }}{\xi }_{\theta }(\tau ),\end{array}$$

$$\begin{array}{rl}\frac{d\mathbf{r}}{dt}& =\lambda \hat{\mathbf{e}}(\theta )-\beta \mathbf{r}+\sqrt{2}{\mathbf{\xi }}_r(t),\\ \frac{d\theta }{dt}& =\sqrt{2}{\xi }_{\theta }(t),\end{array}$$

$$\frac{\partial \mathcal{P}}{\partial t}=-\nabla \cdot [(\lambda \hat{\mathbf{e}}(\theta )-\beta \mathbf{r}-\nabla )\mathcal{P}]+\frac{{\partial }^2\mathcal{P}}{\partial {\theta }^2}.$$

$${\mathcal{L}}_0\mathcal{P}$$

$${\mathcal{L}}_0$$

$${\mathcal{L}}_1$$

$$\mathcal{P}(r,\chi )=\sum _{m=0}^{\infty }{\lambda }^m{\mathcal{P}}^{(m)}(r,\chi ).$$

$${\mathcal{P}}^{(m)}(r,\chi )=\sum _{n,l}{C}_{n,l}^{(m)}{\varphi }_{n,l}(r,\chi ).$$

$${\varphi }_{n,l}(r,\chi )$$

$$C_{n,l}^{(m)}$$

$$C_{n,0}^{(m)}$$

$$\begin{array}{rl}{\mathcal{P}}^{(0)}& =\frac{\beta }{2\pi }{e}^{-\frac{\beta r^2}{2}},{\mathcal{P}}^{(1)}=\frac{{\beta }^{2}r\cos \chi }{2\pi (\beta +1)}{e}^{-\frac{\beta r^2}{2}}\\ {\mathcal{P}}^{(2)}& =\frac{{\beta }^3{r}^2\cos (2\chi )+\beta (\beta +2)(\beta r^2-2)}{8\pi (\beta +1)(\beta +2)}{e}^{-\frac{\beta r^2}{2}.}\end{array}$$

$$P(r)=\sum _{n=0}^{\infty }{\lambda }^{2n}\sqrt{2\pi \beta }{C}_{n,0}^{(2n)}{e}^{-\frac{\beta r^2}{2}}{L}_n(\beta r^2/2),$$

$$P^{\prime \prime }{|}_{r=0}\equiv \sum _{n=0}^{\infty }(n+1){\lambda }^{2n}{C}_{n,0}^{(2n)}(\beta )=0.$$

$$\beta \gtrsim \frac{2+{\lambda }^2}{2{\lambda }^2}.$$

$$\beta \gtrsim {\lambda }^2.$$

$$\begin{array}{rl}\tilde{P}(\rho )& \simeq \frac{\mu k}{D_{\mathrm{eff}}}{e}^{-\frac{\mu k}{2D_{\mathrm{eff}}}{\rho }^2}[1-\frac{\mu k}{D_{\theta }}(\frac{3}{4}-\frac{5}{4}\frac{\mu k}{D_{\mathrm{eff}}}{\rho }^2\\ & +\frac{7}{32}\frac{{\mu }^2{k}^2}{D_{\mathrm{eff}}^2}{\rho }^4)+O\left(\frac{1}{D_{\theta }^2}\right)].\end{array}$$

$$\displaystyle\tilde{P}(\rho)\simeq\frac{1}{Z}\exp\bigg{(}-\frac{\mu k\rho^{2}}{2D_{t}}\bigg{)}~{}I_{0}\bigg{(}\frac{u_{0}\rho}{D_{t}}\bigg{)},$$

$$\begin{array}{rl}\frac{d\mathbf{r}}{dt}& =\lambda \hat{\mathbf{e}}(\theta )-\beta \mathbf{r}+\sqrt{2}{\xi }_r(t),\\ \frac{d\theta }{dt}& =\sqrt{2}{\xi }_{\theta }(t),\end{array}$$

$$\frac{\partial \mathcal{P}}{\partial t}=-\nabla \cdot (\lambda \hat{e}-\beta \mathbf{r}-\nabla )P+\frac{{\partial }^2\mathcal{P}}{\partial {\theta }^2}.$$

$$\frac{\partial \mathcal{P}}{\partial t}$$

$${\mathcal{L}}_0\mathcal{P}$$

$${\mathcal{L}}_0$$

$${\mathcal{L}}_1$$

$$\mathcal{P}(r,\chi )=\sum _{m=0}^{\infty }{\lambda }^m{\mathcal{P}}^{(m)}(r,\chi ).$$

$${\mathcal{L}}_0{\mathcal{P}}^{(0)}=0,{\mathcal{L}}_0{\mathcal{P}}^{(m)}={\mathcal{L}}_1{\mathcal{P}}^{(m-1)},m=1,2,\dots$$

$$\displaystyle\bigg{[}\frac{1}{r}\frac{\partial}{\partial r}\bigg{(}r\frac{\partial}{\partial r}\bigg{)}+\frac{1}{r^{2}}\frac{\partial^{2}}{\partial\chi^{2}}+\beta-\frac{\beta^{2}r^{2}}{4}+\frac{\partial^{2}}{\partial\chi^{2}}\bigg{]}\psi=\nu\psi.\quad~{}~{}~{}$$

$${\psi }_{n,l}(r,\chi )$$

$${\nu }_{n,l}$$

$${\varphi }_{n,l}=e^{-\beta r^2/4}{\psi }_{n,l}.$$

$${\tilde{\varphi }}_{n,l}=e^{\beta r^2/4}{\psi }_{n,l}^{\ast },$$

$${\int }_0^{\infty }drr{\int }_0^{2\pi }d\chi {\tilde{\varphi }}_{n,l}{\varphi }_{n^{\prime },l^{\prime }}={\delta }_{n,n^{\prime }}{\delta }_{l,l^{\prime }}.$$

Full text

active matter, Brownian motion, 2D harmonic trap

I Introduction

The study of active particles has seen an upsurge of interest in recent years for its relevance in describing many non-equilibrium processes Romanczuk et al. (2012); Bechinger et al. (2016). Active Brownian particles (ABPs) and run-and-tumble particles (RTPs) are minimal models for self-propelled active particles Cates and Tailleur (2013). The dynamics of such active particles breaks detailed balance at the microscopic level, and is characterized by properties which are remarkably different from equilibrium systems. Systems composed of interacting active particles are known to exhibit a plethora of exotic phenomenon including flocking Toner et al. (2005); Kumar et al. (2014), clustering Slowman et al. (2016, 2017), motility induced phase separation and segregation Cates and Tailleur (2015); Patch et al. (2017); Redner et al. (2013); Stenhammar et al. (2013), ratchet effects Reichhardt and Reichhardt (2017) etc.

While there exists a considerable body of work on the hydrodynamic description of active matter Marchetti et al. (2013), the number of exactly solvable models that illuminate the novel statistical physics of active particles are few. A number of recent studies, both experimental Li and Tang (2009); Maggi et al. (2014); Takatori et al. (2016); Walsh et al. (2017); Dauchot and Démery (2019) and theoretical Elgeti and Gompper (2015); Prymidis et al. (2015); Jahanshahi et al. (2017); Wagner et al. (2017); Pototsky and Stark (2012); Duzgun and Selinger (2018); Basu et al. (2018); Das et al. (2018), show that even a single active particle can exhibit rich and counter-intuitive physics such as non-Boltzmann distributions peaked away from potential minima Takatori et al. (2016); Dauchot and Démery (2019); Dhar et al. (2019) and clustering Tailleur and Cates (2009). For example, a recent experimental study of active Janus particles, in a two dimensional effectively harmonic trap, surprisingly observed that, in the dilute limit, the trap stiffness can be tuned to induce a crossover in the particle distribution from a Boltzmann-like distribution peaked at the trap center to a strongly active non-Boltzmann distribution with off-centered peaks Takatori et al. (2016). At first sight, the dynamics of a single active particle appears to be a small variation of passive Brownian motion. However, one finds that calculating even the steady state probability distribution is highly non-trivial. Some exact results obtained in Hänggi and Jung (1995); Solon et al. (2015); Malakar et al. (2018); Demaerel and Maes (2018); Dhar et al. (2019); Mallmin et al. (2019); Elgeti and Gompper (2015); Wagner et al. (2017); Pototsky and Stark (2012); Duzgun and Selinger (2018); Basu et al. (2018); Das et al. (2018) indicate the qualitatively rich physics that even single particle active systems can exhibit. Exactly solvable models of isolated active particles are thus important not only for understanding laboratory experiments of self-propelled particles in confined geometries, but they also provide a good starting point to study systems of weakly interacting active particles.

In this paper, we consider the stochastic dynamics of an active Brownian particle (ABP) confined to move in an isotropic potential. Our main result is an exact series solution of the corresponding Fokker-Planck equation in the steady-state. From this exact solution, we find that the radial probability distribution of the ABP has two “phases”: a “passive phase” described by a Boltzmann-like distribution with a peak at the trap center, and an “active phase” described by a non-Boltzmann distribution where the probability distribution is peaked away from the trap center. We find that increasing the trap stiffness induces a transition from the passive phase to the active phase as experimentally observed in Takatori et al. (2016). Our exact solution reveals a surprising prediction of a re-entrant transition to the passive phase upon stiffening the harmonic trap further. Note that these “transitions” are not thermodynamics phase-transitions, but are crossover behaviors in the probability distribution for the ABP. We estimate typical parameter values where such a re-entrant transition can be experimentally observed. Furthermore, the exact solution in our study unifies the various asymptotic limits of the steady-state distribution considered in earlier studies. Our numerical analysis is in good agreement with the analytical results.

The plan of the paper is as follows:- In Sec. (II) we define the precise model and then present a summary of our main analytical results for the steady state distribution, along with their numerical verification. The details of the analytical computations are presented in the appendices. In Sec. (III) we discuss some special limiting cases and finally we conclude with a discussion of our results in Sec. (IV).

II Trapped noninteracting ABPs: Steady state solution

The two-dimensional motion of an ABP in an isotropic harmonic potential ($k\rho^{2}/2$), with position coordinate ${\bm{\rho}}=(\rho\cos\varphi,\rho\sin\varphi)$ and internal angular degree of freedom $\theta$ at time $\tau$, are governed by the Langevin equations:

[TABLE]

where $\mu$ is the translational mobility. The ABP self-propels along the direction $\hat{\bf e}(\theta)=(\cos\theta,\sin\theta)$ with speed $u_{0}$. The Gaussian random variables ${\bm{\xi}}_{r}(\tau)$ and $\xi_{\theta}(\tau)$, with zero mean and unit variance, are uncorrelated in time, and represent translational and rotational noise terms. Defining dimensionless variables ${\bf r}={\bm{\rho}}\sqrt{{D_{\theta}}/{D_{t}}}$ and $t=\tau D_{\theta}$, the Langevin equations take the form

[TABLE]

where $\lambda={u_{0}}/{\sqrt{D_{\theta}D_{t}}}$ and $\beta={\mu k}/{D_{\theta}}$. Note that $\lambda$ is a Péclet number that indicates the relative importance of persistent motion compared to diffusion. Our aim is to calculate the steady state distribution of the trapped ABP. The Fokker-Planck equation for the probability distribution $\mathcal{P}(r,\varphi,\theta,t)$, corresponding to the above Langevin equation, can be obtained using standard methods Risken (1989). For the case of a free ABP in two dimension this was obtained in Sevilla and Sandoval (2015). In our case the corresponding Fokker-Planck equation is

[TABLE]

Note that the the radial symmetry of the problem implies that the steady state joint probability distribution $\mathcal{P}(r,\varphi,\theta)\equiv\mathcal{P}(r,\chi)$ depends only on the combination $\chi=\theta-\varphi$. In the steady state ($\partial_{t}\mathcal{P}=0$), the Fokker-Planck equation can be recast in the following form

[TABLE]

where $\mathcal{L}_{0}$ is an operator which is second-order in $r$ and $\chi$ and contains the terms arising from the confining potential. On the other hand $\mathcal{L}_{1}$ is a first-order operator that arises from the drift-terms in (3). Their explicit forms are given by

[TABLE]

We propose to solve Eq. (4) in the form of a power-series expansion in $\lambda$:

[TABLE]

Using this ansatz in (4) and equating powers of $\lambda$ on both sides, we get $\mathcal{L}_{0}\mathcal{P}^{(0)}=0$ and $\mathcal{L}_{0}\mathcal{P}^{(m)}=\mathcal{L}_{1}\mathcal{P}^{(m-1)}$ for $m\geq 1$. The zeroth-order solution is $\mathcal{P}^{(0)}(r,\chi)=e^{-\beta r^{2}/2}/{Z}$, where $Z=2\pi\int_{0}^{\infty}dr\,re^{-\beta r^{2}/2}$. Note that this is identical to the equilibrium distribution corresponding to the potential $V(r)=r^{2}/2$, if $D_{t}=\mu k_{B}T$. The solution at the next order satisfies $\mathcal{L}_{0}\mathcal{P}^{(1)}=\mathcal{L}_{1}\mathcal{P}^{(0)}$. Thus, if the eigensystem of $\mathcal{L}_{0}$ is explicitly known, then the source term $\mathcal{L}_{1}P^{(0)}$ (obtained from the previous order) can be expanded in the eigenbasis of $\mathcal{L}_{0}$. This allows one to solve for $\mathcal{P}^{(1)}$. Following this strategy at every order $m$, we find (see App. A) that the solutions $\mathcal{P}^{(m)}$ can be written as

[TABLE]

where the summation on $n$ and $l$ is constrained by $m=2n+|l|$ and

[TABLE]

is the right eigenfunction of $\mathcal{L}_{0}$ corresponding to the eigenvalue $\nu_{n,l}=-\beta(2n+|l|)-l^{2}$ with $n\geq 0$, $l$ being integers, $\Gamma$ is the gamma function, and $L_{n}^{|l|}(x)$ is the generalized Laguerre polynomial. The expansion coefficients $C_{n,l}^{(m)}$ satisfy the following recursion relations

[TABLE]

with $C^{(0)}_{0,0}=1/\sqrt{2\beta}$. Solving the above recursion relations, one can, in principle, evaluate $\mathcal{P}^{(m)}$ at any order $m$. For instance, the explicit form of $\mathcal{P}^{(m)}$ for the first few orders are,

[TABLE]

Expressions for few more orders are given in App. B. The analytical evaluation of $C^{(m)}_{n,l}$ becomes tedious with increasing order of $m$. However, they are easily computed numerically. Thus one can evaluate the joint distribution $\mathcal{P}(r,\chi)$ up to any power of $\lambda$. Integrating (7) over the angular variable $\chi$, only $l=0$ terms contribute and we find that the marginal radial distribution function is given by,

[TABLE]

where $L_{n}(z)$ is $n^{\rm th}$ Laguerre polynomial and the coefficients $C_{n,0}^{(2n)}$ can be computed numerically easily from the recursion relations in (9). From numerics, we found that the ratio of two consecutive $C^{(2n)}_{n,0}$ decreases as $\sim n^{-1}$ for large $n$, which ensures that the series in Eq. (11) is convergent.

We numerically solved Eqs. (1) using the Euler-Maruyama method Kloeden and Platen (1992) using a time-step $\Delta t=10^{-4}$ and averaged over $10^{6}$ realizations. Our analytical solution of the Fokker-Planck equation is in good agreement with direct numerical simulation results of Eqs. (1). This is shown in Fig. (1a) where we compare the analytical result for $P(r)$ in Eq. (11) (summed over a finite number of terms, $N$) with those from the simulations. The inset of Fig. (1a) shows the convergence of the series with increasing $M$. We find excellent agreement even for $M=5$.

The distribution $p(\chi)$, of the relative orientation $\chi$, can be obtained by integrating $r$ from the joint distribution $\mathcal{P}(r,\chi)$ (given in Eq. (7)). However, it is more interesting to look at the conditional distribution $p(\chi|r)=\mathcal{P}(r,\chi)/P(r)$ of the orientation $\chi$ at a given $r$. We observe from Fig. (1b) that $p(\chi|r)$ is peaked at $\chi=0$. With increasing $r$, the width of this peak decreases while its height increases. This indicates that the particles which successfully climb the potential must have their orientation $\hat{\bm{e}}(\theta)$ preferentially directed along the radial direction $\hat{\bm{r}}$.

In Fig. (1a), we observe that as $u_{0}$ is increased (at fixed $D_{\theta}$ and $D_{t}$), the position of the peak of $P(r)$ shifts away from the center ($r=0$) of the potential, while $P(r)|_{r=0}$ becomes a minima. We refer to this shift as the passive-to-active transition. The same transition can also be observed by changing $D_{\theta}$ at fixed $u_{0}$ and $D_{t}$. This transition is described by the curve given by

[TABLE]

In Fig. (2a), we plot this curve in the $\lambda-\beta$ plane. We see that at fixed $\beta$, increasing $\lambda$ induces a passive-to-active transition. However, for fixed $\lambda$ greater than a critical $\lambda^{*}$, on increasing $\beta$ from a small value, we first see a passive-to-active transition. Remarkably, upon further increasing $\beta$, we observe a re-entrant active-to-passive transition, as shown in Fig. (2b). The first transition was observed in a recent experiment Takatori et al. (2016), where $\beta$ was varied by changing the trap stiffness $k$. We here predict a re-entrant transition on increasing $k$ further. In Fig. (2c-e), we plot the steady-state distribution of particle positions $\mathbf{r}$ and angles $\theta$ at fixed $\lambda=4$ and varying $\beta$ to represent the three different regions of the phase-diagram in Fig. (2a). We observe that while the distribution is peaked at the center of the trap both in the passive phase (c) and the re-entrant phase (e), the distribution of $\theta$ in (e) is akin to that in the active phase (d). Note that the distribution is not peaked at the center of the trap in (d).

Below we provide a heuristic physical explanation of the two crossovers that we observe.

(i) Passive-to-active crossover — For a small stiffness $k$ we expect that the particle will effectively behave as a passive particle since it undergoes many orientation changes in the time scale ($\tau_{eq}=(\mu k)^{-1}$) required to reach the passive (Gaussian) steady state. The typical extent of this steady state distribution is $\ell_{th}\sim\sqrt{D_{\rm eff}/(\mu k)}$, where $D_{\rm eff}=D_{t}+u_{0}^{2}/(2D_{\theta})$. We note that the time taken by particles to reorder their orientation is given by $\tau_{rot}\sim D_{\theta}^{-1}$, in which time they travel a distance $\ell_{act}=u_{0}\tau_{rot}$. As we increase $k$, the length scale $\ell_{th}$ decreases. We then expect that the passive-like steady state will be stable provided $\ell_{act}\lesssim\ell_{th}$. Otherwise the particles tend to move out and the distribution starts having peaks away from the centre. This gives the passive-to-active crossover and the condition $\ell_{act}\gtrsim\ell_{th}$ for this gives

[TABLE]

Interestingly we note that even for large $\lambda$ one needs a finite value of $\beta\approx 1/2$ to see the passive-to-active crossover.

(ii) Active-to-passive crossover — this crossover is observed on further increasing $k$ (or equilvalentsly $\beta$). Note that this ‘reentrant passive phase’ is markedly different from the usual passive phase in that, the distribution of the orientation of the particle captures a notrivial signature of activity as explained in Figs. (2)(c)-(e). In this case we propose the following physical explanation. The activity pushes out the particle radially so that peaks appear at a distance given by $r_{peak}\sim u_{0}/\mu k$. If this distance is small compared to the thermal length scale $\ell^{\prime}_{th}\sim\sqrt{D_{t}/(\mu k)}$ (in this case the particles angular motion is restricted and does not contribute to $D_{\rm eff}$), then we will not see any off-centre peak, even though the particles do mainitain a strong radial orientation. This then leads to the re-entrant crossover occuring when $\ell^{\prime}_{th}\gtrsim r_{peak}$ which reduces to

[TABLE]

These two estimates, Eqs. (13,14), also indicate that there exists a critical $\lambda\approx 1.13$ below which one does not see any of the crossovers. These heuristic estimates for the crossovers are consistent with our numerical observations in Fig. (2)a, though a precise quantitative verification would require more work.

Moreover, in the limits $D_{\theta}\to 0$ and $D_{\theta}\to\infty$, we find closed form analytical expressions for this probability distribution. In particular, our approach offers a perturbative solution for the limit $D_{t}\to 0$, which is non-trivial to obtain for a problem without any translational diffusion ($D_{t}=0$) Basu et al. (2018). We discuss these limits in the next section.

III Special limits

We now discuss certain limiting cases of the radial probability distribution $\tilde{P}(\rho)=(D_{\theta}/D_{t})P(r)$. We first consider the case with zero translational noise, i.e., $D_{t}=0$. Figure (3a) shows a plot of $\tilde{P}(\rho)$ for $D_{t}=0$ obtained from our simulations. We observe that $\tilde{P}(\rho)$ has a finite support with a peak around $r=1$, in sharp contrast to the usual Boltzmann distribution. The singular nature of $\tilde{P}(\rho)$ for $D_{t}=0$ does not allow us to readily obtain an analytical expression in this case. However, we numerically evaluate our general series expression from Eq. (11) and compare $\tilde{P}(\rho)$ with numerical simulation results for small $D_{t}$. As shown in Fig. (3a), our approach thus allows a systematic way to attain the $D_{t}\to 0$ limit, although the computational cost of evaluating Eq. (11) escalates quickly for small values of $D_{t}$.

We next consider the limit $D_{\theta}\rightarrow\infty$, $u_{0}\rightarrow\infty$ such that $\beta\to 0$ with $u_{0}^{2}/D_{\theta}=\lambda^{2}D_{t}$ held constant. In this limit, the stochasticity in $u_{0}\hat{{\bf e}}(\theta)$ arising from the rotational diffusion in Eq. (1) approaches the limit of a zero-mean Gaussian white noise of strength $u_{0}^{2}/(2D_{\theta})$ Basu et al. (2018). This adds to the translational white noise ${\bm{\xi}}_{r}(t)$. As a result, we expect a Gaussian distribution of the form P(\rho)\simeq(\mu k)/D_{\rm eff}\exp\big{[}-\mu k\rho^{2}/(2D_{\rm eff})\big{]} with an effective diffusion constant $D_{\rm eff}=D_{t}+{u_{0}^{2}}/{(2D_{\theta})}$, as shown in Basu et al. (2018) for the case $D_{t}=0$. One can, however, find corrections to this result systematically at different orders in $1/D_{\theta}$. To accomplish this, we first solve the recursion relations in Eq. (9) perturbatively in powers of $\beta$ to get $C^{(2m)}_{m,0}\simeq(-1/2)^{m}\sqrt{{\beta}/{2\pi}}[1-\beta(7m^{2}+m)/8+O(\beta^{2})]$. Using this in Eq. (11) and performing some straightforward calculations we get $\tilde{P}(\rho)$ that takes the form,

[TABLE]

in the $D_{t}\rightarrow 0$ limit (calculation steps and the expression for finite $D_{t}$ are in App. C). Figure (3b) validates this asymptotic result with numerical simulations. Note that in the above expression one can take $D_{t}\to 0$ limit which would correspond to the case of an active Brownian particle without any translational noise.

Finally, we consider the opposite limit $D_{\theta}\rightarrow 0$, i.e., $\beta\to\infty$, in which the timescale for rotational diffusion of the orientation of $\theta$ diverges as $\sim D_{\theta}^{-1}$. The $\beta\to\infty$ limit can also be obtained by taking $k\to\infty$ limit. In this limit, the radial distribution, as shown in Pototsky and Stark (2012), is approximately

[TABLE]

where $Z$ is a normalization constant and $I_{0}(y)$ is the modified Bessel function of zeroth order. In order to arrive at this result from our general expression in Eq. (11), we look at the asymptotic behavour of the coefficients $C_{n,l}^{(m)}$ for large $\beta$. Note that for large $\beta$, Eq. (9) simplifies to the form $C^{(m)}_{n,l}=[C^{(m-1)}_{n,l-1}\sqrt{(n+|l|)}-C^{(m-1)}_{n-1,l+1}\sqrt{n}]/{[\sqrt{2\beta}(2n+|l|)]}$. It can be easily checked that this equation is satisfied by $C_{n,l}^{(2n+l)}={(-1)^{n}}\sqrt{\beta/(2\pi)}[{n!(n+l)!\beta^{2n+l}2^{2n+l}}]^{-1/2}$ for $n\geq 0$ and $\l\geq 0$. Using this result in Eq. (11) leads to Eq. (16) (see App. C for details).

IV Summary and conclusion

To conclude, we have obtained an exact series solution for the non-equilibrium steady state of an ABP, with translational diffusion, confined to a two dimensional harmonic trap. Our analytical results are in good agreement with explicit numerical simulations. We explain a recently observed passive-to-active transition of a harmonically confined ABP with increasing trap-stiffness. Furthermore, we predict a surprising active-to-passive reentrant transition upon a further increase in the confinement strength. The prediction of this reentrant transition is amenable to experimental validation. In fact, using the values from the experimental setup of Takatori et al. (2016), $D_{t}\approx 0.1\mu m^{2}/s$, $D_{\theta}^{-1}\approx 25s$, $u_{0}\approx 0.25\mu m/s$, i.e. Péclet number $~{}\lambda\approx 4$, we expect to observe the passive-to-active transition around $\beta_{1}\approx 0.35$, i.e., $\mu k_{1}\approx 0.01s^{-1}$ and the re-entrant active-to-passive transition around $\beta_{2}\approx 6.8$, i.e., $\mu k_{1}\approx 0.27s^{-1}$. These are plausible parameter values that could be accessed experimentally.

We thank S. C. Takatori for discussions on the experimental aspects. KM’s research is supported by Long Term Visiting Student’s Program-2018, ICTS-TIFR. A. Das acknowledges the CEFIPRA postdoctoral fellowship hosted at ICTS-TIFR. AK acknowledges support from DST grant under project No. ECR/2017/000634. KVK’s research is supported by the Department of Biotechnology, India, through a Ramalingaswami reentry fellowship and by the Max Planck Society and the Department of Science and Technology, India, through a Max Planck Partner Group at ICTS-TIFR. AK and AD acknowledge the support of the project 5604-2 of the Indo-French Centre for the Promotion of Advanced Research (IFCPAR). We acknowledge support of the Department of Atomic Energy, Government of India, under project no.12-R$\&$D-TFR-5.10-1100.

Appendix A Details of calculational of the steady state distribution

The dimensionless Langevin equation for the active Brownian particle in a two dimensional harmonic trap is given by,

[TABLE]

where the dimensionless parameters $\lambda$ and $\beta$ are defined as $\lambda=\frac{u_{0}}{\sqrt{D_{\theta}D_{t}}},\beta=\frac{\mu k}{D_{\theta}}.$ The corresponding Fokker-Planck equation is

[TABLE]

In the two dimensional polar coordinates the above equation takes the form

[TABLE]

where $\chi=\theta-\varphi$. As discussed in the main text, in the steady state given by ${\partial\mathcal{P}}/{\partial t}=0$, we expect the solution to have azimuthal symmetry, $\mathcal{P}(r,\varphi,\theta)=\mathcal{P}(r,\chi)$. Hence, in the steady state, after making the replacements $\frac{\partial}{\partial\theta}=\frac{\partial}{\partial\chi}$, $\frac{\partial}{\partial\varphi}=-\frac{\partial}{\partial\chi}$, we find that the equation satisfied by the steady state probability distribution $\mathcal{P}(r,\chi)$ can be written in the form (Eq. (4) of main text)

[TABLE]

where

[TABLE]

We try a series solution of the form

[TABLE]

Inserting this in Eq. (19) and equating powers of $\lambda$ on both sides, we get

[TABLE]

The first equation ($m=0$) above can be solved to give $\mathcal{P}^{(0)}(r,\chi)=e^{-\beta r^{2}/2}/{Z}$, where $Z=2\pi\int_{0}^{\infty}dr\,re^{-\beta r^{2}/2}$. As mentioned in the main text, if the eigenfunctions and eigenvalues of $\mathcal{L}_{0}$ are known explicitly and if they form a complete basis, then at any order $m$, the term on the right-hand-side of Eq. (23) calculated from the previous order ($m-1$) solution can be expanded in the eigenbasis of $\mathcal{L}_{0}$ and one obtains $\mathcal{P}^{(m)}(r,\chi)$. The procedure is detailed in the following.

Eigensystem of $\mathcal{L}_{0}$: Let $\phi(r,\chi)$ be the eigenvector of $\mathcal{L}_{0}$ corresponding to the eigenvalue $\nu$ satisfying the equation $\mathcal{L}_{0}\phi=\nu\phi$. With the transformation $\phi=e^{-\beta V(r)/2}\,\psi$, we find that $\psi$ satisfies the equation $H\psi=\nu\psi$ where $H=e^{-\beta V(r)/2}\mathcal{L}_{0}e^{\beta V(r)/2}$ is a Hermitian Schrödinger operator. Explicitly, the eigenvalue equation is

[TABLE]

The above equation without the last $\partial^{2}/\partial\chi^{2}$ term corresponds to the Schrödinger equation of an isotropic two-dimensional harmonic oscillator, for which the eigenfunctions ($\psi_{n,l}$) and eigenvalues ($\Lambda_{n,l}$) are known exactly Yáñez et al. (1994):

[TABLE]

with $\Lambda_{n,l}=-\beta(2n+|l|)~{}$ where $n\geq 0,~{}l$ are integers, $\Gamma$ is the gamma function, and $L_{n}^{|l|}(x)$ is the generalized Laguerre polynomial. Note that the presence of the last $\partial^{2}/\partial\chi^{2}$ term in Eq. (24) does not change the eigenfunctions while the eigenvalues get modified to

[TABLE]

The $\psi_{n,l}$ form a set of basis functions. Since the Fokker-Planck operator $\mathcal{L}_{0}$ is not a Hermitian operator, its right and left eigenfunctions are different. The right eigenfunctions are given by

[TABLE]

It can be shown that the left eigenfunctions are given by

[TABLE]

and these satisfy the orthogonality condition

[TABLE]

Next, we expand the solutions, at different order, in the $\mathcal{L}_{0}$-basis:

[TABLE]

Inserting this in Eq. (23), and using the orthogonality of the basis states, we find $C_{n,l}^{(m)}=({1}/{\nu_{n,l}})\int_{0}^{\infty}dr\,r\int_{0}^{2\pi}d\chi~{}\mathcal{L}_{1}\mathcal{P}^{(m-1)}\tilde{\phi}_{n,l}$. Using the specific form of $\mathcal{L}_{1}$ given in Eq. (21), we expand $\mathcal{L}_{1}\mathcal{P}^{(m-1)}\equiv f^{(m)}(r,\chi)$ in the same basis set and, on using the above mentioned orthonormality, we find that the coefficients satisfy the recursion relations ( Eqs. (9) in the main text)

[TABLE]

with $C^{(0)}_{0,0}=1/\sqrt{2\beta}$. Thus we have

[TABLE]

Inserting these expansions into Eq. (23) and using the orthonormality of $\psi_{n,l}$, we get

[TABLE]

Now, using the explicit forms of the eigenfunctions

[TABLE]

and defining N^{(m)}_{n,l}=C^{(m)}_{n,l}\sqrt{\frac{2n!\big{(}\frac{\beta}{2}\big{)}^{|l|+1}}{2\pi(n+|l|)!}}, we get

[TABLE]

and the source term

[TABLE]

where $x=\beta r^{2}/2$. Using the orthonormality of the Laguerre polynomials and the azimuthal eigenfunctions, we finally obtain the following recursion relation

[TABLE]

and from radial symmetry,

[TABLE]

The ‘boundary’ recursion equations are given by

[TABLE]

The recursion relations can be solved iteratively, starting from the “initial” condition

[TABLE]

and the boundary condition, $C^{(m)}_{-1,l}=0,~{}l\geq 0.$ The structure of the recursion relations is illustrated in Fig. (4). Here the red filled discs represent $C^{(m)}_{n,l}$ with the horizontal and the vertical axes representing $l$ and $n$ respectively, while the dashed lines correspond to constant $m=2n+l$. The coefficients $C_{n,l}^{(m)}$, at a given $m$, can be obtained from the coefficients $C_{n,l-1}^{(m-1)}$ and $C_{n-1,l+1}^{(m-1)}$ on the $(m-1)$ line, as indicated by the two arrows meeting at the point $(n,l)$. One can think of the index $m$ as a time index and $n,l$ as spatial indices. Thus, starting from the initial condition localized at $n=0,~{}l=0$, these recursions prescribe a time evolution in the $m$ direction. Following this procedure, we can in principle find all the coefficients for arbitrary $n,l$, and hence $\mathcal{P}^{(m)}$ at any order $m$.

So, starting with $C^{(0)}_{0,0}$, all other $C^{(m)}_{n,l}$ can be obtained using the above recursion relations, and hence we get an explicit series solution for $\mathcal{P}(r,\chi)$ given by Eqs. (22, 28). In the following we list first few coefficients,

[TABLE]

and so on. Integrating over $\chi$ we can obtain the radial distribution $P(r)$. It is easy to see that only the $l=0$ terms contribute to the radial distribution which gives,

[TABLE]

where $L_{n}$ is Laguerre polynomial of order $n$. Simplifying the above, we get,

[TABLE]

which is Eq. (11) of the main article. In the next section we demonstrate the results for first few $\mathcal{P}^{(m)}(r,\chi)$ calculated in this method.

Appendix B Explicit expressions for $\mathcal{P}^{(m)}(r,\chi)$ for $m=1,2,3,4$

Zeroth order

($\mathcal{P}^{(0)}$):

[TABLE]

This is the equilibrium distribution (for $\lambda=0$).

First order

($\mathcal{P}^{(1)}$):

[TABLE]

Second order

($\mathcal{P}^{(2)}$):

[TABLE]

Third order

($\mathcal{P}^{(3)}$):

[TABLE]

Fourth order

($\mathcal{P}^{(4)}$):

[TABLE]

Appendix C Special limiting cases

With the replacements, $\beta=\frac{\mu k}{D_{\theta}},~{}r=\rho\sqrt{\frac{D_{\theta}}{D_{t}}},~{}\lambda=\frac{u_{0}}{\sqrt{D_{\theta}D_{t}}},~{}t=\tau D_{\theta}$, the probability distribution $\tilde{P}(\rho)$ in the dimensionful variables will be,

[TABLE]

i.e., in equilibrium when activity is absent, $\tilde{P}^{(0)}(\rho)=\frac{\mu k}{D_{t}}e^{-\frac{\mu k}{2D_{t}}\rho^{2}}$, and in the non-equilibrium steady state in presence of nonzero activity,

[TABLE]

C.1 Limit $D_{\theta}\rightarrow\infty$, $u_{0}\rightarrow\infty$, i.e. $\beta\rightarrow 0$ and $\lambda=$ finite

In the limit $\beta\rightarrow 0$, the coefficients $C^{(m)}_{n,l}$ (see Fig. (4)) take the form :

[TABLE]

to leading order in $\beta$. It is evident from the table that the $l=0$ terms are proportional to $\sqrt{\beta}$, and the higher $l$ terms come with higher powers of $\beta$. Our aim is to evaluate $C^{(2m)}_{m,0}$, which, from Eq. (27) in the main text, depends on $C^{(2m-1)}_{m-1,l}$, that in turn depends on other $n~{}{\rm and}~{}l$ terms. However since the terms with $l$=2 and above appears with higher powers in $\beta$, to obtain the leading contribution in $\beta\rightarrow 0$ limit, we truncate the recursion equation at $l=1$,

[TABLE]

Combining these two, we obtain,

[TABLE]

Using the above in the series solution in dimensionful quantities, the steady state probability distribution takes the form,

[TABLE]

where we have replaced the sum using generating function of Laguerre polynomials, $\sum_{n=0}^{\infty}L_{n}(x)t^{n}=\frac{1}{1-t}\exp[-\frac{tx}{1-t}]$, with $t=-\frac{\lambda^{2}}{2}$ and $x=\frac{\mu k\rho^{2}}{2D_{t}}$. For $D_{t}\rightarrow 0$ , this reduces to the expression, $\tilde{P}(\rho)=\frac{2\mu kD_{\theta}}{u_{0}^{2}}\exp[-\frac{2\mu kD_{\theta}\rho^{2}}{u_{0}^{2}}]$ obtained in Basu et al. (2018).

Solution to $\mathcal{O}(\beta)$:

To find the leading correction to $\mathcal{O}(\beta)$, we keep one more order in the recursion for the $C^{(m)}_{n,l}$. The relevant recursion relations truncated at $l=2$ are,

[TABLE]

The above set of equations can be simplified as follows. Replacing $n$ by $n-1$ in Eq.(48), we get, $C^{(2n-1)}_{n-1,1}=\sqrt{\frac{\beta}{2}}~{}\frac{C^{(2n-2)}_{n-1,0}\sqrt{n}-C^{(2n-2)}_{n-2,2}\sqrt{n-1}}{1+\beta(2n-1)}$. $C^{(2n-2)}_{n-2,2}$ is obtained from Eq.(49) as, $C^{(2n-2)}_{n-2,2}=\frac{1}{2}\sqrt{\frac{\beta}{2}}~{}\frac{\sqrt{n}}{2+\beta(n-1)}C^{(2n-3)}_{n-2,1}$, which, using Eq.(47), gives, $C^{(2n-2)}_{n-2,2}=\frac{\beta}{2}\frac{\sqrt{n(n-1)}}{2+\beta(n-1)}C^{(2n-2)}_{n-1,0}$. Inserting the last expression in the equation for $C^{(2n-1)}_{n-1,1}$ and applying Eq.(47) once again, we finally obtain a simplified relation between the coefficients in the radial distribution function,

[TABLE]

Its solution correct to $\mathcal{O}(\beta)$ is,

[TABLE]

with $C^{(0)}_{0,0}=\sqrt{\frac{\beta}{2\pi}}$. Therefore, the radial distribution function takes the form,

[TABLE]

In the limit $D_{\theta}\rightarrow\infty$, i.e., $\beta\rightarrow 0$, the quantity $y^{2}=\beta r^{2}$ remains finite. So the correct radial distribution to $\mathcal{O}(\beta)$ in terms of the dimensionless quantities is,

[TABLE]

Using again the generating function of Laguerre polynomials, $G_{x}(t)=\sum_{n=0}^{\infty}t^{n}L_{n}(x)=\frac{1}{1-t}\exp(-\frac{tx}{1-t})$, we evaluate the terms in the sum above to get:

[TABLE]

and

[TABLE]

Inserting these in Eq.(LABEL:series), we get

[TABLE]

Reverting to dimensionful quantities,

[TABLE]

where $D_{\rm{eff}}=D_{t}+\frac{v^{2}}{2D_{\theta}}$. In the limit $D_{t}\rightarrow 0~{}(\mbox{i.e.}~{}\lambda\rightarrow\infty)$, the radial distribution function takes the form,

[TABLE]

which is Eq. (15) of the main text.

C.2 Limit $D_{\theta}\rightarrow 0~{}(\beta\rightarrow\infty)$

In this limit, the recursion relation takes the form

[TABLE]

solving which we get the following table for the coefficients

[TABLE]

This leads us to make an ansatz for the coefficients $C^{(m)}_{n,l}$

[TABLE]

which indeed satisfies the recursion relation (56). Therefore the coefficients corresponding to $l=0$ are given by

[TABLE]

Using Eq. (58) in the series solution of $\tilde{P}(\rho)$ we obtain

[TABLE]

Evaluating the double sum after the following rearrangement \sum_{m=0}^{\infty}\sum_{i=0}^{m}a_{m}b_{i}^{m}x^{i}=\sum_{m=0}^{\infty}\big{(}\sum_{i=m}^{\infty}a_{i}b_{m}^{i})x^{m} and using I_{\alpha}(x)=\sum_{m=0}^{\infty}\frac{1}{m!\,\Gamma(m+\alpha+1)}\big{(}\frac{x}{2}\big{)}^{2m+\alpha}, the probability distribution in the limit $D_{\theta}\rightarrow 0$ becomes,

[TABLE]

which is Eq. (16). This is in accordance with Eq. (8) of Pototsky and Stark (2012).