Extreme fluctuations of active Brownian motion

As a simple but fairly general description of active Brownian motion, we consider a system that is characterized by a continuous variable x for the displacement and a discrete set $\{i\}$ representing mesoscopic states, which could be inaccessible to direct observation. The switching between these states is modeled as a Markovian network with transition rates w_ij from state i to state j. Assuming the system to be homogeneous along the coordinate x, the rates w_ij do not depend on x. The system performs self-propelled motion along the coordinate x. As essential ingredient of the model, distinguishing active particles from merely externally pulled passive particles, the self-propulsion force f_i depends on the mesoscopic state i. At the present level of description, the microscopic mechanism that actually generates the propulsion force in not yet included in the model. As we will show below in section 4, the states $\{i\}$ and the forces f_i emerge via coarse-graining from a more detailed description. The switching between the states $\{i\}$ will typically be caused by thermal fluctuations with transition rates w_ij obeying detailed balance. In principle, however, it could involve an active component as well.

The Langevin equation for the displacement x(t) at time t is

$$\begin{eqnarray}&&\dot{x}(t)={\mu }_{i}{f}_{i}+{\zeta }_{i}(t).\end{eqnarray} \tag{ 1 }$$

For fixed i, this Langevin equation corresponds to the overdamped dynamics of a particle with mobility ${\mu }_{i}$ driven by a force f_i. In addition, the particle experiences thermal noise ${\zeta }_{i}(t)$ with the statistical averages $\langle {\zeta }_{i}(t)\rangle =0$ and $\langle {\zeta }_{i}(t){\zeta }_{i}({t}^{\prime })\rangle =2{D}_{i}\delta (t-{t}^{\prime })$. The diffusion coefficient D_i is related to the mobility via the Einstein relation ${D}_{i}={\mu }_{i}{k}_{{\rm{B}}}T$, with Boltzmann's constant k_B and the temperature T of the surrounding heat bath.

For a Janus particle, the variable i can label different geometrical alignments, while x denotes its position in a one-dimensional channel (or the projection of the position of a free Janus particle onto the x-axis). If the self-propulsion mechanism generates no torque on the particle, the internal transitions w_ij describe the rotational diffusion of the particle and satisfy detailed balance. For each state i a different component f_i of the propulsion force acts in the direction of the coordinate x. If the particle is not spherically symmetric, the mobility ${\mu }_{i}$ can depend on the state i. For biological microswimmers, the internal state can also encode different modes of motility, such as the 'run' and the 'tumble' mode in many types of bacteria.

The evolution of the joint probability distribution of x and i is governed by a combination of a Fokker-Planck equation for x and a master equation for i

$$\begin{eqnarray}&&{\partial }_{t}{p}_{i}(x,t)=(-{\mu }_{i}{f}_{i}{\partial }_{x}+{D}_{i}{\partial }_{x}^{2}){p}_{i}(x,t)+\displaystyle \sum _{j}{L}_{{ij}}{p}_{j}(x,t),\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}&&{L}_{{ij}}\equiv {w}_{{ji}}-{r}_{i}{\delta }_{{ij}}\end{eqnarray} \tag{ 3 }$$

and ${r}_{i}\equiv {\sum }_{{\ell }}{w}_{i{\ell }}$ is the exit rate from state i. With this joint probability distribution, the distributions for the internal variable i and the position variable x follow as

$$\begin{eqnarray}&&{p}_{i}(t)={\int }_{-\infty }^{\infty }{\rm{d}}x\;{p}_{i}(x,t)\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&p(x,t)=\displaystyle \sum _{i}{p}_{i}(x,t),\end{eqnarray} \tag{ 5 }$$

respectively. In the long time limit, the distribution of the internal variable attains a stationary distribution ${p}_{i}^{{\rm{s}}}$, which is independent of the initial distribution and satisfies

$$\begin{eqnarray}&&\displaystyle \sum _{j}{L}_{{ij}}{p}_{j}^{{\rm{s}}}=0.\end{eqnarray} \tag{ 6 }$$

The ensemble average of the velocity of the system in the long time limit can be calculated directly from the stationary distribution ${p}_{i}^{{\rm{s}}}$ as

$$\begin{eqnarray}&&v\equiv \underset{t\to \infty }{\mathrm{lim}}\langle x(t)/t\rangle =\underset{t\to \infty }{\mathrm{lim}}\langle \dot{x}(t)\rangle =\displaystyle \sum _{i}{p}_{i}^{{\rm{s}}}{\mu }_{i}{f}_{i}.\end{eqnarray} \tag{ 7 }$$

Along individual trajectories, the time averaged velocity $u\equiv x/t$ may differ from the ensemble average v. However, the probability of such fluctuations decreases exponentially with the length t of the trajectory,

$$\begin{eqnarray}&&p(x,t)\sim {{\rm{e}}}^{-{th}(x/t)},\end{eqnarray} \tag{ 8 }$$

characterized by of the so-called rate function or large deviation function [33]

$$\begin{eqnarray}&&h(u)\equiv -\underset{t\to \infty }{\mathrm{lim}}\displaystyle \frac{1}{t}\mathrm{ln}p(x={ut},t).\end{eqnarray} \tag{ 9 }$$

Expanding h(u) to second order around its minimum $h(v)=0$ yields a Gaussian approximation to the probability distribution, characterized by the effective diffusion coefficient

$$\begin{eqnarray}&&{D}_{\mathrm{eff}}\equiv \underset{t\to \infty }{\mathrm{lim}}\displaystyle \frac{1}{2t}(\langle {x}^{2}\rangle -{\langle x\rangle }^{2})\equiv \underset{t\to \infty }{\mathrm{lim}}{C}_{2}/(2t)=1{/h}^{\prime\prime }(v).\end{eqnarray} \tag{ 10 }$$

Higher derivatives of h(u) at u = v relate to the higher order cumulants C_n. The rate function captures the full range of non-Gaussian fluctuations in the long-time limit and is often easier to interpret than a large set of cumulants.

The standard approach to calculating the large deviation function [33, 35] starts with the transformation

$$\begin{eqnarray}&&{g}_{i}(\lambda ,t)\equiv {\int }_{-\infty }^{\infty }{\rm{d}}x\;{{\rm{e}}}^{\lambda x}\;{p}_{i}(x,t)\end{eqnarray} \tag{ 11 }$$

with a real variable λ. The Gaussian nature of the noise ensures that tails of ${p}_{i}(x,t)$ decay like a Gaussian, therefore the integral converges for all λ. The evolution equation (2) then transforms to

$$\begin{eqnarray}&&{\partial }_{t}{g}_{i}(\lambda ,t)=\displaystyle \sum _{j}{{ \mathcal L }}_{{ij}}(\lambda ){g}_{j}(\lambda ,t)\end{eqnarray} \tag{ 12 }$$

with the matrix

$$\begin{eqnarray}&&{{ \mathcal L }}_{{ij}}(\lambda )\equiv {L}_{{ij}}+({\mu }_{i}{f}_{i}\lambda +{D}_{i}{\lambda }^{2}){\delta }_{{ij}}.\end{eqnarray} \tag{ 13 }$$

Equation (12) can be solved using of an expansion in eigenvectors of the matrix ${{ \mathcal L }}_{{ij}}(\lambda )$. In the long-time limit, the cumulant generating function, defined as

$$\begin{eqnarray}&&\alpha (\lambda )\equiv \underset{t\to \infty }{\mathrm{lim}}\displaystyle \frac{1}{t}\mathrm{ln}{g}_{i}(\lambda ,t)=\underset{j}{\mathrm{max}}\;\mathrm{Re}[{\alpha }_{j}(\lambda )],\end{eqnarray} \tag{ 14 }$$

depends neither on i nor on the initial distribution and is given by the eigenvalue ${\alpha }_{j}(\lambda )$ with the largest real part. The Perron–Frobenius theorem ensures that this eigenvalue has no imaginary part. The cumulants C_n of the distribution of x in the long time limit follow as

$$\begin{eqnarray}&&\underset{t\to \infty }{\mathrm{lim}}{C}_{n}/t={{\partial }_{\lambda }^{n}\alpha (\lambda )| }_{\lambda =0}.\end{eqnarray} \tag{ 15 }$$

The cumulant generating function is convex, since it can be written as the Legendre transform

$$\begin{eqnarray}&&\alpha (\lambda )=\underset{u}{\mathrm{max}}\;[\lambda u-h(u)]\end{eqnarray} \tag{ 16 }$$

of the rate function h(u). If the latter is convex as well, $\alpha (\lambda )$ is differentiable for all λ and the Legendre transform can be inverted to

$$\begin{eqnarray}&&h(u)=\underset{\lambda }{\mathrm{max}}[\lambda u-\alpha (\lambda )].\end{eqnarray} \tag{ 17 }$$

Typical features of the cumulant generating function and the rate function can be illustrated in a two-state model as a simple example. The two states with drift velocities ${v}_{i}\equiv {\mu }_{i}{f}_{i}$ and diffusion constants D_i are connected via the transition rates $\eta {w}_{12}$ and $\eta {w}_{21}$. The dimensionless parameter η is introduced to study the effect of the time scale of the internal transitions. This system can be used, for example, to model a single molecular motor walking along a one-dimensional track, as shown in the top row of figure 1. For motors that can switch between forward and backward direction, like helicases [42] or complexes of kinesin and dynein [43], the parameters are chosen as ${v}_{1}\gt 0$, ${v}_{2}\lt 0$. and ${D}_{1}={D}_{2}$. For a unidirectional motor that can detach from and re-attach to the track [44, 45], ${v}_{1}\gt 0$, ${v}_{2}=0$ and ${D}_{2}\gt {D}_{1}$ would be an appropriate choice. The matrix (20) for these types of two state systems reads

$$\begin{eqnarray}{ \mathcal L }(\lambda ,\eta )=\left(\begin{array}{cc}-\eta {w}_{12}+{v}_{1}\lambda +{D}_{1}{\lambda }^{2} & \eta {w}_{21}\\ \eta {w}_{12} & -\eta {w}_{21}+{v}_{2}\lambda +{D}_{2}{\lambda }^{2}\end{array}\right).\end{eqnarray} \tag{ 18 }$$

Its maximal eigenvalue is, as in general for 2 × 2-matrices,

$$\begin{eqnarray}&&\alpha (\lambda ,\eta )=\mathrm{tr}{ \mathcal L }(\lambda ,\eta )/2+\sqrt{{(\mathrm{tr}{ \mathcal L }(\lambda ,\eta ))}^{2}/4-\mathrm{det}{ \mathcal L }(\lambda ,\eta )},\end{eqnarray} \tag{ 19 }$$

where $\mathrm{tr}{ \mathcal L }$ and $\mathrm{det}{ \mathcal L }$ denote the trace and the determinant of ${ \mathcal L }$, respectively. This cumulant generating function is plotted in figure 1 for the two simple motor models. As a universal characteristic, the cumulant generating function exhibits kink-like features for slow internal dynamics (small η) and approaches a smooth parabola for fast internal dynamics (large η). For the model with the dichotomous drift velocity, there is only one kink in the cumulant generating function around $\lambda =0$. The cumulant generating function for the model with stochastic detaching exhibits a second kink at $\lambda \simeq 0.5$. The corresponding rate functions h(u) in figure 1 have been calculated numerically via the Legendre transformation (17). Thereby, the kinks transform to intervals with constant (or zero) slope.

Zoom In Zoom Out Reset image size

The general properties of the cumulant generating function and the rate function can be understood by considering the limiting cases of slow and fast internal dynamics. In both cases, the timescales of the internal dynamics and the translational diffusion process get separated, allowing for an explicit calculation of the rate function. As in the two simple examples above, we rewrite the transition rates as ${w}_{{ij}}\mapsto \eta {w}_{{ij}}$ and ${r}_{i}\mapsto \eta {r}_{i}$, yielding the modified matrix (13)

$$\begin{eqnarray}&&{{ \mathcal L }}_{{ij}}(\lambda ,\eta )=\eta {L}_{{ij}}+({\mu }_{i}{f}_{i}\lambda +{D}_{i}{\lambda }^{2}){\delta }_{{ij}}\end{eqnarray} \tag{ 20 }$$

with L_ij from equation (3). The stationary distribution ${p}_{i}^{{\rm{s}}}$ is unaffected by this scaling, hence the velocity v in equation (7) is independent of η as well. In contrast, the effective diffusion coefficient (10) depends strongly on η. With an expansion around $\lambda =0$, it can be calculated exactly as

$$\begin{eqnarray}&&{D}_{\mathrm{eff}}(\eta )=\displaystyle \sum _{i}{p}_{i}^{{\rm{s}}}{D}_{i}+{\eta }^{-1}\displaystyle \sum _{i}{q}_{i}{v}_{i},\end{eqnarray} \tag{ 21 }$$

where the vector q_i (the first order correction to the eigenvector of ${ \mathcal L }(\lambda )$) is defined as the unique solution of the linear system of equations

$$\begin{eqnarray}&&\displaystyle \sum _{j}{L}_{{ij}}{q}_{j}=(v-{\mu }_{i}{f}_{i}){p}_{i}^{{\rm{s}}},\qquad \displaystyle \sum _{i}{q}_{i}=0.\end{eqnarray} \tag{ 22 }$$

The case of slow internal dynamics corresponds to $\eta \ll 1$. In the limit $\eta \to 0$, the eigenvalues of the matrix (20) are simply given by its diagonal elements ${\alpha }_{i}^{0}(\lambda )\equiv {\mu }_{i}{f}_{i}\lambda +{D}_{i}{\lambda }^{2}$. Thus, the limit of the cumulant generating function is given by

$$\begin{eqnarray}&&\underset{\eta \to 0}{\mathrm{lim}}\alpha (\lambda ,\eta )=\underset{i}{\mathrm{max}}\;{\alpha }_{i}^{0}(\lambda ).\end{eqnarray} \tag{ 23 }$$

This function can have kinks where different eigenvalues ${\alpha }_{i}^{0}(\lambda )$ intersect. Typically, if not all ${\mu }_{i}{f}_{i}$ are equal, there is at least one such kink at $\lambda =0$. In the limit $\eta \to 0$ the cumulant generating function is actually non-analytic at these kinks. In contrast, for all finite values $\eta \gt 0$, the Perron–Frobenius theorem ensures that all eigenvalues of ${{ \mathcal L }}_{{ij}}(\lambda ,\eta )$ are distinct, so that the intersections vanish and $\alpha (\lambda )$ is still analytic everywhere [46]. The increasing curvature of $\alpha (\lambda )$ with decreasing η is also reflected in the divergence of the effective diffusion coefficient ${D}_{\mathrm{eff}}(\eta )={\partial }_{\lambda }^{2}\alpha (\lambda ,\eta ){| }_{\lambda =0}$ in equation (21). For non-zero but still small values of η, corrections to the eigenvalue entering in the cumulant generating function (23) can be obtained via perturbation theory as

$$\begin{eqnarray}&&{\alpha }_{i}(\lambda ,\eta )={\alpha }_{i}^{0}(\lambda )-\eta {r}_{i}+{\eta }^{2}\displaystyle \sum _{i\ne j}\displaystyle \frac{{w}_{{ij}}{w}_{{ji}}}{{\alpha }_{i}^{0}(\lambda )-{\alpha }_{j}^{0}(\lambda )}+{ \mathcal O }({\eta }^{3}).\end{eqnarray} \tag{ 24 }$$

This approximation works well for most values of λ, but not in the region of the kinks, where the eigenvalues become degenerate.

For the calculation of the rate function it plays a crucial role in which order the long-time limit in equation (9) and the limit $\eta \to 0$ are performed. If the limit $\eta \to 0$ is taken first, the finite time probability distribution simply reads

$$\begin{eqnarray}&&{p}_{i}(x,t)=\displaystyle \sum _{i}{\int }_{-\infty }^{\infty }{\rm{d}}{x}_{0}\;{p}_{i}({x}_{0},0){(4\pi {D}_{i}t)}^{-1/2}\mathrm{exp}\left[\displaystyle \frac{{-(x-{x}_{0}-{\mu }_{i}{f}_{i}t)}^{2}}{4{D}_{i}t}\right].\end{eqnarray} \tag{ 25 }$$

Provided that the initial distribution ${p}_{i}({x}_{0},0)$ is non-zero for all i, the corresponding rate function (9) is then

$$\begin{eqnarray}&&\tilde{h}(u)=\underset{i}{\mathrm{min}}\;{h}_{i}(u)\equiv \;\underset{i}{\mathrm{min}}\;\displaystyle \frac{{(u-{\mu }_{i}{f}_{i})}^{2}}{4{D}_{i}},\end{eqnarray} \tag{ 26 }$$

which is typically a non-convex function. In contrast, if the long-time limit is taken first, the rate function is uniquely determined by the Legendre transform (17) of the analytic function $\alpha (\lambda )$ for any finite $\eta \gt 0$. The limit $\eta \to 0$ then leads to the convex envelope h(u) of $\tilde{h}(u)$ in equation (26) as rate function. This order of the limits is the one relevant for approximating the rate function at small, non-zero values of η. For every kink in the cumulant generating function $\alpha (\lambda )$ at $\lambda ={\lambda }_{{\rm{k}}}$, the rate function h(u) depends linearly on u with the slope ${h}^{\prime }(u)={\lambda }_{{\rm{k}}}$ on the interval

$$\begin{eqnarray}&&\alpha ({\lambda }_{{\rm{k}}}^{-})\lt u\lt \alpha ({\lambda }_{{\rm{k}}}^{+}).\end{eqnarray} \tag{ 27 }$$

For the simple two-state models discussed above the limit (23) is shown in red in figure 1. In both models there is a kink in the cumulant generating function at $\lambda =0$. In the model with the stochastic detaching, the parabolas for the $\eta \to 0$ limit intersect a second time at $\lambda ={v}_{1}/({D}_{2}-{D}_{1})$. In the Legendre transformed picture of the rate function, the kink at $\lambda =0$ results in a flat plateau (slope 0) around u = 0. This plateau forms the convex envelope between the minima of the parabolas in equation (26) at $u={v}_{1}$ and $u={v}_{2}$, respectively. For the model with stochastic detaching, there is a further interval with linear slope ${v}_{1}/({D}_{2}-{D}_{1})$ between $u={v}_{1}({D}_{2}+{D}_{1})/({D}_{2}-{D}_{1})$ and $u=2{v}_{1}{D}_{2}/({D}_{2}-{D}_{1})$. Within these intervals, fluctuations in the displacement of the system are dominated by the fluctuations of the internal variable. Beyond these ranges (especially for $\lambda \to \pm \infty$), fluctuations are dominated by translational noise and occur almost exclusively if the internal variable happens to be constant most of the time.

With increasing internal speed η of the internal dynamics, the kinks in the cumulant generating function and the linear regions in the rate function are softened until both functions approach a parabolic shape shown in blue in figure 1. This behavior can be understood by considering the limit of fast internal dynamics in general. For $\eta \to \infty$, the eigenvector associated with the largest eigenvalue of the matrix (20) approaches the stationary distribution ${p}_{i}^{{\rm{s}}}$. Starting from this distribution, first order perturbation theory in $1/\eta$ yields the cumulant generating function

$$\begin{eqnarray}&&\alpha (\lambda )=v\lambda +\bar{D}{\lambda }^{2}+{ \mathcal O }(1/\eta )\end{eqnarray} \tag{ 28 }$$

and the rate function

$$\begin{eqnarray}&&h(u)=\displaystyle \frac{{(u-v)}^{2}}{4\bar{D}}+{ \mathcal O }(1/\eta ),\end{eqnarray} \tag{ 29 }$$

which corresponds effectively to the diffusion of a driven passive particle with drift velocity v and diffusion coefficient $\bar{D}\equiv {p}_{i}^{{\rm{s}}}{D}_{i}$. This diffusion coefficient is equal to the limiting value ${D}_{\mathrm{eff}}(\eta \to \infty )$ in equation (21).

Janus particles represent another important class of active particles that can be described within the formalism presented above. As an effective model, illustrated in figure 2, we assign to a particle a propulsion force with constant magnitude F₀ acting in the direction of the instantaneous orientation of the particle. Since the particle is subject to rotational diffusion, the component of the propulsion force in the direction of an individual coordinate x varies with time. We assume that both the rotational and the translational diffusion coefficient are independent of the orientation, which should be a good approximation for spherically symmetric particles. In the absence of external potentials, the symmetry of the configuration leads to a vanishing average velocity in the x-direction. Since the propulsion mechanism drives the system out of equilibrium, breaking this symmetry leads to a non-vanishing velocity. For example, such a symmetry breaking could be achieved by an asymmetrically patterned substrate, modeled as an asymmetric periodic potential [14, 47] or through interaction with an external field, for example the gravitational field in gravitaxis [18, 19] or a magnetic field in magnetotaxis [21]. A symmetry breaking that can easily be implemented in the present formalism is a dipolar interaction of the Janus particle with a homogeneous external field. We identify the coordinate x with the direction of the field and assume that the dipole moment of the Janus particle is aligned with the orientation of the propulsion force. The strength of the interaction, i.e., the dipole moment multiplied by the field strength in units of ${k}_{{\rm{B}}}T$, is denoted as B in the following. We assume that the direction of the coordinate x is chosen such that $B\gt 0$.

Zoom In Zoom Out Reset image size

First we consider two-dimensional, i.e. 'disk-like', Janus particles, as in [29] (without field) and, more recently, [48]. In this case the angle θ between the propulsion force and the direction of the coordinate x fully characterizes the internal state of the particle, replacing the discrete variable i in equation (1). The Langevin equation for the translation of the Janus particle reads

$$\begin{eqnarray}&&\dot{x}(t)=\mu {F}_{0}\mathrm{cos}\theta +\zeta (t)\end{eqnarray} \tag{ 30 }$$

with the mobility μ and correlations of the noise $\zeta (t)$ as in (1) with a uniform diffusion coefficient D. The rotational motion of the particle is described by the Langevin equation

$$\begin{eqnarray}&&\dot{\theta }(t)=-{\mu }_{{\rm{r}}}\;{k}_{{\rm{B}}}{TB}\mathrm{sin}\theta +{\zeta }_{{\rm{r}}}(t)\end{eqnarray} \tag{ 31 }$$

with the rotational mobility ${\mu }_{{\rm{r}}}\equiv {D}_{{\rm{r}}}/({k}_{{\rm{B}}}T)$. The rotational noise term obeys $\langle {\zeta }_{{\rm{r}}}(t)\rangle =0$, $\langle {\zeta }_{{\rm{r}}}(t){\zeta }_{{\rm{r}}}({t}^{\prime })\rangle =2{D}_{{\rm{r}}}\delta (t-{t}^{\prime })$, and is uncorrelated with the translational noise $\zeta (t)$. The evolution of the joint probability $p(x,\theta ,t)$ is governed by the Fokker–Plank equation

$$\begin{eqnarray}&&{\partial }_{t}p(x,\theta ,t)=[{D}_{{\rm{r}}}{L}_{{\rm{r}}}(B)-\mu {F}_{0}\mathrm{cos}\theta \;{\partial }_{x}+D{\partial }_{x}^{2}]p(x,\theta ,t),\end{eqnarray} \tag{ 32 }$$

which is a continuous version of equation (2) with the transition matrix L_ij replaced by the operator

$$\begin{eqnarray}&&{L}_{{\rm{r}}}(B)\equiv {\partial }_{\theta }^{2}+B{\partial }_{\theta }\mathrm{sin}\theta \end{eqnarray} \tag{ 33 }$$

for the rotational diffusion of the particle in the field. Analogously, the cumulant generating function $\alpha (\lambda )$ for the distribution of the displacement x is given by the maximal eigenvalue of the operator

$$\begin{eqnarray}&&{ \mathcal L }(\lambda )={D}_{{\rm{r}}}{L}_{{\rm{r}}}(B)+\mu {F}_{0}\mathrm{cos}\theta \;\lambda +D{\lambda }^{2}.\end{eqnarray} \tag{ 34 }$$

By rescaling to the dimensionless variable

$$\begin{eqnarray}&&z\equiv \lambda \mu {F}_{0}/{D}_{{\rm{r}}},\end{eqnarray} \tag{ 35 }$$

this maximal eigenvalue can be written as

$$\begin{eqnarray}&&\alpha (\lambda )=D{\lambda }^{2}+{D}_{{\rm{r}}}\;a(B,\lambda \mu {F}_{0}/{D}_{{\rm{r}}}),\end{eqnarray} \tag{ 36 }$$

where $a(B,z)$ is defined as the maximal value of a for which the differential equation

$$\begin{eqnarray}&&\tilde{{ \mathcal L }}(B,z,\theta )\psi (\theta )\equiv [{L}_{{\rm{r}}}(B)+z\mathrm{cos}\theta ]\psi (\theta )=a\psi (\theta )\end{eqnarray} \tag{ 37 }$$

has a $2\pi$-periodic solution.

Without external field, i.e. for B = 0, this equation is identical to the Mathieu equation and $a(0,z)$ can be expanded as [49]

$$\begin{eqnarray}&&a(0,z)=\displaystyle \frac{1}{2}{z}^{2}-\displaystyle \frac{7}{32}{z}^{4}+\displaystyle \frac{29}{144}{z}^{6}+{ \mathcal O }({z}^{8}).\end{eqnarray} \tag{ 38 }$$

Due to the obvious symmetry

$$\begin{eqnarray}&&a(0,z)=a(0,-z),\end{eqnarray} \tag{ 39 }$$

all odd terms vanish in this expansion. From this expansion one can derive the low order cumulants of the distribution $p(x,t)$ as

$$\begin{eqnarray}&&{C}_{2}=[D+{D}_{{\rm{r}}}{(\mu {F}_{0}/{D}_{{\rm{r}}})}^{2}]\;t,{C}_{4}=-\displaystyle \frac{21}{4}\;{D}_{{\rm{r}}}{(\mu {F}_{0}/{D}_{{\rm{r}}})}^{4}\;t,{C}_{6}=145\;{D}_{{\rm{r}}}{(\mu {F}_{0}/{D}_{{\rm{r}}})}^{6}\;t,\end{eqnarray} \tag{ 40 }$$

in accordance with the results for C₂ and C₄ in [29]. It should be noted, however, that this Taylor expansion has a finite radius of convergence beyond which one has to rely on a numerical evaluation of the eigenvalue.

Remarkably, the trivial symmetry (39) persists even in the more interesting case of a non-zero interaction parameter B, though with a shifted center of symmetry. For every eigenfunction $\psi (\theta )$ of $\tilde{{ \mathcal L }}(B,z)$, the shifted function $\psi (\theta +\pi )$ is an eigenfunction of the adjoint operator ${\tilde{{ \mathcal L }}}^{\dagger }(B,-B-z)$, since

$$\begin{eqnarray}{\tilde{{ \mathcal L }}}^{\dagger }(B,-B-z,\theta )\psi (\theta +\pi ) & = & [{\partial }_{\theta }^{2}-B{\partial }_{\theta }\mathrm{sin}\theta -z\mathrm{cos}\theta ]\psi (\theta +\pi )\\ & = & \tilde{{ \mathcal L }}(B,z,\theta +\pi )\psi (\theta +\pi ).\end{eqnarray} \tag{ 41 }$$

As a consequence, the eigenvalue $a(B,z)$ satisfies the symmetry

$$\begin{eqnarray}&&a(B,z)=a(B,-B-z).\end{eqnarray} \tag{ 42 }$$

In terms of the cumulant generating function (36) the symmetry (42) can be expressed as

$$\begin{eqnarray}&&\alpha ({\lambda }_{0}+\lambda )=\alpha ({\lambda }_{0}-\lambda )+4D{\lambda }_{0}\lambda \end{eqnarray} \tag{ 43 }$$

with

$$\begin{eqnarray}&&{\lambda }_{0}\equiv -{D}_{{\rm{r}}}B/(2\mu {F}_{0}).\end{eqnarray} \tag{ 44 }$$

This new symmetry is different from the well-known Gallavotti–Cohen symmetry [35, 50]. The latter holds if the observable of interest is proportional to the entropy production in the medium along individual trajectories. However, this is not the case for the displacement of Janus particles. For example, such a particle could move back and forth with parallel propulsion force for each direction, thus producing entropy without net displacement. After the Legendre transformation (17), the symmetry (43) can be expressed in terms of the rate function as

$$\begin{eqnarray}&&h({u}_{0}+u)=2{\lambda }_{0}u+h({u}_{0}-u)\end{eqnarray} \tag{ 45 }$$

with

$$\begin{eqnarray}&&{u}_{0}\equiv 2D{\lambda }_{0}=-{k}_{{\rm{B}}}{{TD}}_{{\rm{r}}}B/{F}_{0}.\end{eqnarray} \tag{ 46 }$$

For the probability distribution p(x) of the displacement we obtain the fluctuation relation

$$\begin{eqnarray}&&\underset{t\to \infty }{\mathrm{lim}}\mathrm{ln}\displaystyle \frac{p({u}_{0}t+{\rm{\Delta }}x)}{p({u}_{0}t-{\rm{\Delta }}x)}=-2{\lambda }_{0}{\rm{\Delta }}x.\end{eqnarray} \tag{ 47 }$$

In contrast to the well-established fluctuation theorems that are related to entropy production, this fluctuation relation refers to a center of symmetry that is not located at x = 0. Instead, the center of symmetry at $x={u}_{0}t$ moves in backward direction at constant speed.

For an efficient numerical calculation of the eigenvalue $a(B,z)$ we expand the eigenfunction as $\psi (\theta )={c}_{0}+2{\sum }_{k=1}^{\infty }{c}_{k}\mathrm{cos}(k\theta )$, which leads to the representation of equation (37) as ${\sum }_{{k}^{\prime }}{\tilde{{ \mathcal L }}}_{{{kk}}^{\prime }}{c}_{{k}^{\prime }}={{ac}}_{k}$ with the infinite matrix

$$\begin{eqnarray}&&{\tilde{{ \mathcal L }}}_{{{kk}}^{\prime }}=-{k}^{2}{\delta }_{k,{k}^{\prime }}+(z+{Bk}){\delta }_{k-1,{k}^{\prime }}/2+(z-{Bk}){\delta }_{k+1,{k}^{\prime }}/2\end{eqnarray} \tag{ 48 }$$

(up to the exception ${\tilde{{ \mathcal L }}}_{01}=z$). Truncating this matrix after about 10 to 20 rows and columns already gives very good approximations of $a(B,z)$. The result of this numerical calculation is shown in figure 3 for selected fields B. These plots display the symmetry (42). At the center of this symmetry, the function $a(B,z)$ exhibits a kink-like feature, which becomes more pronounced with increasing B. Further away from this kink, the function $a(B,z)$ has a rather small curvature. In figure 4, the numerical result for $a(B=5,z)$ has been used to generate a family of rate functions for varying rotational diffusion coefficients D_r.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The two limiting cases with time scale separation can be analyzed for Janus particles in a similar fashion as for the discrete systems in section 2.3. The operator (34) is a continuous version of the matrix (20), with the rotational diffusion coefficient D_r playing the role of the parameter η. For the extreme cases ${D}_{{\rm{r}}}\to 0$ and ${D}_{{\rm{r}}}\to \infty$, we can again make use of the separation of timescales to calculate the cumulant generating function analytically.

The limit of slow internal dynamics corresponds to small D_r and thus large values of $z$ in equation (37). The relevant eigenvector of the term $z\mathrm{cos}\theta$ for large positive $z$ is simply $\delta (\theta )$ with eigenvalue $z$. For large negative values of $z$ the relevant eigenvector is $\delta (\theta -\pi )$ with eigenvalue $-z$. Plugging this asymptotic behavior into equation (36) yields the cumulant generating function

$$\begin{eqnarray}&&\alpha (\lambda )=D{\lambda }^{2}+\mu {F}_{0}| \lambda | \end{eqnarray} \tag{ 49 }$$

in the limit ${D}_{{\rm{r}}}\to 0$, which is analogous to equation (23) in the discrete case. The resulting rate function is shown in red in figure 4.

Although being asymptotically correct, in the sense that $a(B,z)/| z| \to 1$ for large $| z|$, the approximation $a(B,z)\sim | z|$ is rather rough for realistic values of $z$. In particular, in figure 3, this approximation captures neither the shift of the kink nor the slope in the wings of the function $a(B,z)$. The problem is that for these values of $z$ the eigenfunction $\psi (\theta )$ is not yet delta-like but still has a finite width. Thus, for a better approximation of the eigenvalue for intermediate values of $z,$ we choose a representation of the delta-function with finite width as eigenvector. Since the eigenvector $\psi (\theta )$ for intermediate $z$ somehow interpolates between the known eigenvector ${p}^{{\rm{s}}}(\theta )$ at $z=0$ and $\delta (\theta )$ for $z\to \infty$, the ansatz

$$\begin{eqnarray}&&\tilde{\psi }(\theta )=\mathrm{exp}[{B}_{\mathrm{eff}}(B,z)\mathrm{cos}\theta ]\end{eqnarray} \tag{ 50 }$$

with an effective field ${B}_{\mathrm{eff}}(B,z)$ that increases with $z$ should be suitable. This parameter can be fixed along with the approximate eigenvalue $\tilde{a}(B,z)$ by requiring $[\tilde{{ \mathcal L }}(B,z,\theta )\tilde{\psi }(\theta )]/\tilde{\psi }(\theta )=\tilde{a}(B,z)+{ \mathcal O }({\theta }^{4})$ for the expansion around the maximum of the eigenvector at $\theta =0$. The resulting effective field is

$$\begin{eqnarray}&&{B}_{\mathrm{eff}}(B,z)=[(B-1/2)+\sqrt{{(B+1/2)}^{2}+2z}]/2\end{eqnarray} \tag{ 51 }$$

and the approximate eigenvalue becomes

$$\begin{eqnarray}&&\tilde{a}(B,z)=z+B-{B}_{\mathrm{eff}}(B,z).\end{eqnarray} \tag{ 52 }$$

Although this result is exact only for $z=0$ and $z\to \infty$, the numerical results in figure 3 show that there is a good agreement with the actual eigenvalue for all values of $z\gt -B/2$ for strong fields where the eigenvectors are more concentrated around $\theta =0$. For $z\lt -B/2$, we can make use of the symmetry (42) to obtain a global approximation for $a(B,z)$.

The limit of fast internal dynamics is obtained for ${D}_{{\rm{r}}}\to \infty$. For an expansion in z, we use the stationary distribution

$$\begin{eqnarray}&&{p}^{{\rm{s}}}(\theta )=\displaystyle \frac{\mathrm{exp}(B\mathrm{cos}\theta )}{2\pi {I}_{0}(B)}\end{eqnarray} \tag{ 53 }$$

to obtain

$$\begin{eqnarray}&&a(B,z)={\int }_{0}^{2\pi }{\rm{d}}\theta \;{p}^{{\rm{s}}}(\theta )\mathrm{cos}\theta \;z+{ \mathcal O }({z}^{2})=\displaystyle \frac{{I}_{1}(B)}{{I}_{0}(B)}z+{ \mathcal O }({z}^{2}),\end{eqnarray} \tag{ 54 }$$

where the I_n(B) denotes the modified Bessel functions of the first kind. Plugging this approximation into equation (36) yields the parabolic approximation of the cumulant generating function analogous to equation (28) and the corresponding rate function (29) with the velocity $v=\mu {F}_{0}{I}_{1}(B)/{I}_{0}(B)$ and diffusion coefficient $\bar{D}=D$. In figure 4, this limiting case is shown as a blue parabola.

Similar to the case of discrete internal degrees of freedom, the cumulant generating function (36) of a Janus particle exhibits a kink-like feature. This kink gets smeared out with increasing speed of the internal dynamics, i.e., with increasing D_r, until the function approaches the parabola $\alpha (\lambda )=v\lambda +D{\lambda }^{2}$. As visible in the plots of the function $a(B,z)$ in figure 3, the location of the kink coincides with the center of the symmetry (43) at $\lambda ={\lambda }_{0}$. As explained for the discrete model in section 2.3, this kink leads in the Legendre transformed rate function to an interval with constant slope ${\lambda }_{0}$. Based on equation (27) and the limits $-1\lt {\partial }_{z}a(B,z)\lt 1$, the range where this linear slope occurs can be limited to ${u}_{0}-\mu {F}_{0}\lt u\lt {u}_{0}+\mu {F}_{0}$. In this range, the probability distribution of the displacement does not decay like a Gaussian, as it would be common for diffusion processes, but rather exponentially. Loosely speaking, this slower decay in the probability is due to the fact that extreme fluctuations in this range can be assisted by internal fluctuations of the rotational degree of freedom.

For three-dimensional, spherical Janus particles we write the orientation in spherical coordinates $(\theta ,\phi )$. A distribution over these coordinates is to be understood as a probability per solid angle, hence integrals over the azimuthal angels θ require $\mathrm{sin}\theta$ as Jacobi determinant. For the joint probability distribution $p(x,\theta ,\phi ,t)$, the rotational operator in the Fokker–Planck equation (32) gets modified to

$$\begin{eqnarray}&&{L}_{{\rm{r}}}(B)=\displaystyle \frac{1}{\mathrm{sin}\theta }\displaystyle \frac{\partial }{\partial \theta }\left(\mathrm{sin}\theta \displaystyle \frac{\partial }{\partial \theta }\right)+\displaystyle \frac{1}{{\mathrm{sin}}^{2}\theta }\displaystyle \frac{{\partial }^{2}}{\partial {\phi }^{2}}+\displaystyle \frac{B}{\mathrm{sin}\theta }\displaystyle \frac{\partial }{\partial \theta }{\mathrm{sin}}^{2}\theta ,\end{eqnarray} \tag{ 55 }$$

Eigenfunctions of the accordingly modified operator ${ \mathcal L }(\lambda )$ in equation (34) can be written as products of a θ-dependent and a ϕ-dependent part. However, since any dependence on ϕ decreases the corresponding eigenvalue, we can focus on eigenfunctions of the form $\psi (\theta )$ and thus drop the part of ${L}_{{\rm{r}}}(B)$ acting on ϕ. The steps leading to the operator $\tilde{{ \mathcal L }}(B,z,\theta )$ in equation (37) then remain unaltered in three-dimensions. The adjoint operator, defined with respect to the natural scalar product on the unit sphere, can be shown to satisfy the relation

$$\begin{eqnarray}&&{\tilde{{ \mathcal L }}}^{\dagger }(B,-2B-z,\theta )=\tilde{{ \mathcal L }}(B,z,\theta +\pi ).\end{eqnarray} \tag{ 56 }$$

As a consequence, the eigenvalue $a(B,z)$ satisfies the symmetry

$$\begin{eqnarray}&&a(B,z)=a(B,-2B-z)\end{eqnarray} \tag{ 57 }$$

and the property (43) of the cumulant generating function still holds, but with a different center of symmetry ${\lambda }_{0}=-{D}_{{\rm{r}}}B/(\mu {F}_{0})$.

A numerical calculation of the function $a(B,z)$ is possible via an expansion of the eigenfunctions in Legendre polynomials ${P}_{{\ell }}(\mathrm{cos}\theta )$. The matrix form of $\tilde{{ \mathcal L }}(B,z,\theta )$, analogous to (48), then reads

$$\begin{eqnarray}&&{\tilde{{ \mathcal L }}}_{{\ell }{{\ell }}^{\prime }}=-{\ell }({\ell }+1){\delta }_{{\ell },{{\ell }}^{\prime }}+\displaystyle \frac{{\ell }}{2{\ell }-1}[z+B({\ell }+1)]{\delta }_{{\ell }-1,{{\ell }}^{\prime }}+\displaystyle \frac{{\ell }+1}{2{\ell }+3}[z-B{\ell }]{\delta }_{{\ell }+1,{{\ell }}^{\prime }}\end{eqnarray} \tag{ 58 }$$

with ${\ell },{{\ell }}^{\prime }\geqslant 0$. The dominant eigenvalue can be determined numerically in a truncated approximation of this matrix. The result of the numerics shown in figure 3 reveals that the kink in the cumulant generating function at the center of symmetry is somewhat less pronounced than for the two-dimensional model.

For small values of $z,$ the function $a(B,z)$ can be expanded similarly to equation (54) as

$$\begin{eqnarray}&&a(B,z)=(\mathrm{coth}B-1/B)z+{ \mathcal O }({z}^{2}).\end{eqnarray} \tag{ 59 }$$

Accordingly, the average velocity of the Janus particle is given by $v=\mu {F}_{0}(\mathrm{coth}B-1/B)$. The effective field approximation analogous to equation (50) yields for $z\gt -B$

$$\begin{eqnarray}&&{B}_{\mathrm{eff}}=[B-1+\sqrt{{(B+1)}^{2}+2z}]/2\end{eqnarray} \tag{ 60 }$$

and

$$\begin{eqnarray}&&\tilde{a}(B,z)=B+z+1-\sqrt{{(B+1)}^{2}+2z},\end{eqnarray} \tag{ 61 }$$

which can be continued to $z\lt -B$ via the symmetry (57).