Extreme fluctuations of active Brownian motion

In active Brownian motion, an internal propulsion mechanism interacts with translational and rotational thermal noise and other internal fluctuations to produce directed motion. We derive the distribution of its extreme fluctuations and identify its universal properties using large deviation theory. The limits of slow and fast internal dynamics give rise to a kink-like and parabolic behavior of the corresponding rate functions, respectively. For dipolar Janus particles in two and three dimensions interacting with a field, we predict a novel symmetry akin to, but different from, the one related to entropy production. Measurements of these extreme fluctuations could thus be used to infer properties of the underlying, often hidden, network of states.


Introduction
Stochastic thermal noise prevents predicting the dynamics of small systems exactly. In biological systems, one observes a persistent interplay between this noise and the goal of the system to keep vitally important processes running in a preferred direction. Such directed processes take place under non-equilibrium conditions, which have to be maintained by a permanent supply of energy. For instance, molecular motors burn chemical fuel molecules, typically ATP, to displace cargo in a preferred direction [1,2]. On larger length scales, many types of bacteria and some eucaryotic cells use flagella to generate propulsion forces that allow directed swimming [3,4]. Due to recent advances in nanotechnology, artificial systems with a variety of propulsion mechanisms have been built. For instance, Janus particles are colloidal particles with two differently coated halves, which react differently to external heating or to chemical non-equilibrium conditions and thereby generate self-phoretic motion [5][6][7][8]. The superposition of the directed motion in these systems with the ubiquitous Brownian motion is classified as active Brownian motion [9][10][11][12][13], for which phenomena like the interaction with confining boundaries [14][15][16][17] and external fields [18][19][20][21] and emergent hydrodynamic and collective properties like phase separation and swarming [22][23][24][25][26][27][28] have been intensely investigated recently.
An important characteristic of active Brownian motion is the probability distribution of fluctuations in the displacement of particles. As a consequence of the central limit theorem, typical fluctuations during long time intervals can be well described by a Gaussian distribution. In contrast, extreme fluctuations are usually non-Gaussian and depend crucially on the microscopic details of the underlying system. In this paper, we provide a full characterization of extreme fluctuations for active Brownian motion on large time scales. Earlier studies have analyzed the non-Gaussian behavior in terms of the low order cumulants of the distribution for the particular case of Janus particles [29][30][31][32].
As a mathematical tool, we make use of large deviation theory to calculate the so-called rate function, which captures the exponential contributions to a probability distribution [33]. In stochastic thermodynamics, this approach has led to the formulation of fluctuation theorems as a fundamental property of extreme fluctuations [34][35][36]. These fluctuation theorems, however, can only be applied when the system under consideration does not contain any slow hidden degrees of freedom [37]. If one is interested in the displacement of an active particle, different internal states affecting the propulsion force must often be considered as such a hidden degree of freedom.
So far, the only study of large deviations in the context of active Brownian motion was conducted experimentally for an asymmetric particle interacting with granular matter under nonequilibrium conditions, which revealed a fluctuation theorem for the integrated velocity of the particle [38,39] as an instance of the isometric fluctuation relation [40]. Recent experimental and numerical evidence suggests that this fluctuation relation holds for phoretic self-propelled particles as well [41]. Yet, in these studies the velocity was measured in a frame fixed to the particle. For active particles that undergo rotational diffusion, as for example Janus particles, this integrated velocity differs from the displacement perceived by an external observer. Our results are formulated in the lab frame and can therefore be applied even if the rotational diffusion of the particle cannot be observed or only with insufficient precision.
The paper is organized as follows: In section 2, we present a formalism using discrete internal states and biased diffusion for the propulsion. To illustrate this model, molecular motors are presented as a toy model before we analyze two analytically tractable limiting cases. In section 3, we discuss the extreme fluctuations of the displacement of a Janus particle with a dipole moment interacting with a homogeneous field, for which we derive a novel symmetry similar to a fluctuation theorem. In section 4, we generalize the formalism from section 2, allowing for a more detailed description of the propulsion mechanism. We conclude in section 5.

Langevin formalism
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 Sec. 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 iṡ For fixed i, this Langevin equation corresponds to the overdamped dynamics of a particle with mobility µ i driven by a force f i . In addition, the particle experiences thermal noise ζ i (t) with the statistical averages ζ i (t) = 0 and ζ i (t)ζ i (t ) = 2D i δ(t − t ). The diffusion coefficient D i is related to the mobility via the Einstein relation D i = µ i k 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 µ 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, where and r i ≡ w i 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 and respectively. In the long time limit, the distribution of the internal variable attains a stationary distribution p s i , which is independent of the initial distribution and satisfies j L ij p s j = 0.
The ensemble average of the velocity of the system in the long time limit can be calculated directly from the stationary distribution p s i as v ≡ lim Along individual trajectories, the time averaged velocity u ≡ x/t may differ from the ensemble average v. However, the probability of such fluctuations decreases exponentially with the length t of the trajectory, characterized by of the so-called rate function or large deviation function [33] h(u) ≡ − lim 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 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 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 with the matrix Eq. (12) can be solved using of an expansion in eigenvectors of the matrix L ij (λ). In the long-time limit, the cumulant generating function, defined as depends neither on i nor on the initial distribution and is given by the eigenvalue α j (λ) 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 The cumulant generating function is convex, since it can be written as the Legendre transform of the rate function h(u). If the latter is convex as well, α(λ) is differentiable for all λ and the Legendre transform can be inverted to

Toy model: Two-state model for molecular motors
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 ≡ µ i f i and diffusion constants D i are connected via the transition rates ηw 12 and η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 Fig. 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 > 0, v 2 < 0. and D 1 = D 2 . For a unidirectional motor that can detach from and re-attach to the track [44,45], v 1 > 0, v 2 = 0 and D 2 > D 1 would be an appropriate choice. The matrix (20) for these types of two state systems reads Its maximal eigenvalue is, as in general for 2 × 2-matrices, where tr L and det L denote the trace and the determinant of L, respectively. This cumulant generating function is plotted in Fig. 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 λ = 0. The cumulant generating function for the model with stochastic detaching exhibits a second kink at λ 0.5. The corresponding rate functions h(u) in Fig. 1 have been calculated numerically via the Legendre transformation (17). Thereby, the kinks transform to intervals with constant (or zero) slope. . The left column refers to the motor model with dichotomous drift velocity and the right column to the motor model with stochastic detaching from the track. The functions are shown as thin curves with the internal (dimensionless) timescale η ranging from 0.01 to 10 in 10 steps with logarithmic spacing. The thick blue curves show the limit η → ∞ and the thick red curves the limit η → 0 in Eq. (23) with dashed continuations indicating the constituting parabolas α 0 i (λ) and h i (u), respectively. Parameters for the left column: w 12 = 1, w 21 = 2, v 1,2 = ±1, D 1 = 1 = D 2 . Right column:

Limiting cases: Slow and fast internal dynamics
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 → ηw ij and r i → ηr i , yielding the modified matrix (13) with L ij from Eq. (3). The stationary distribution p s i is unaffected by this scaling, hence the velocity v in Eq. (7) is independent of η as well. In contrast, the effective diffusion coefficient (10) depends strongly on η. With an expansion around λ = 0, it can be calculated exactly as where the vector q i (the first order correction to the eigenvector of L(λ)) is defined as the unique solution of the linear system of equations The case of slow internal dynamics corresponds to η 1. In the limit η → 0, the eigenvalues of the matrix (20) are simply given by its diagonal elements α 0 Thus, the limit of the cumulant generating function is given by This function can have kinks where different eigenvalues α 0 i (λ) intersect. Typically, if not all µ i f i are equal, there is at least one such kink at λ = 0. In the limit η → 0 the cumulant generating function is actually non-analytic at these kinks. In contrast, for all finite values η > 0, the Perron-Frobenius theorem ensures that all eigenvalues of L ij (λ, η) are distinct, so that the intersections vanish and α(λ) is still analytic everywhere [46]. The increasing curvature of α(λ) with decreasing η is also reflected in the divergence of the effective diffusion coefficient D eff (η) = ∂ 2 λ α(λ, η)| λ=0 in Eq. (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 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 Eq. (9) and the limit η → 0 are performed. If the limit η → 0 is taken first, the finite time probability distribution simply reads Provided that the initial distribution p i (x 0 , 0) is non-zero for all i, the corresponding rate function (9) is thenh 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 α(λ) for any finite η > 0. The limit η → 0 then leads to the convex envelope h(u) ofh(u) in Eq. (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 α(λ) at λ = λ k , the rate function h(u) depends linearly on u with the slope h (u) = λ k on the interval For the simple two-state models discussed above the limit (23) is shown in red in Fig. 1. In both models there is a kink in the cumulant generating function at λ = 0. In the model with the stochastic detaching, the parabolas for the η → 0 limit intersect a second time at λ = v 1 /(D 2 −D 1 ). In the Legendre transformed picture of the rate function, the kink at λ = 0 results in a flat plateau (slope 0) around u = 0. This plateau forms the convex envelope between the minima of the parabolas in Eq. (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 = 2v 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 λ → ±∞), 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 Figure 2. Model for a Janus particle. The internal state is characterized by the azimuthal angle θ between the propulsion force F 0 and the direction of the coordinate x. The propulsion force is aligned with the dipole moment, that interacts with the homogeneous field B in the x-direction. approach a parabolic shape shown in blue in Fig. 1. This behavior can be understood by considering the limit of fast internal dynamics in general. For η → ∞, the eigenvector associated with the largest eigenvalue of the matrix (20) approaches the stationary distribution p s i . Starting from this distribution, first order perturbation theory in 1/η yields the cumulant generating function and the rate function which corresponds effectively to the diffusion of a driven passive particle with drift velocity v and diffusion coefficientD ≡ p s i D i . This diffusion coefficient is equal to the limiting value D eff (η → ∞) in Eq. (21).

Model
Janus particles represent another important class of active particles that can be described within the formalism presented above. As an effective model, illustrated in Fig. 2, we assign to a particle a propulsion force with constant magnitude F 0 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 B T , is denoted as B in the following. We assume that the direction of the coordinate x is chosen such that B > 0.

Disk-like Janus particles
First we consider two dimensional, i.e. "disk-like", Janus particles, as in Ref. [29] (without field) and, more recently, Ref. [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 Eq. (1). The Langevin equation for the translation of the Janus particle readṡ with the mobility µ and correlations of the noise ζ(t) as in (1) with a uniform diffusion coefficient D. The rotational motion of the particle is described by the Langevin equatioṅ with the rotational mobility µ r ≡ D r /(k B T ). The rotational noise term obeys ζ r (t) = 0, ζ r (t)ζ r (t ) = 2D r δ(t − t ), and is uncorrelated with the translational noise ζ(t). The evolution of the joint probability p(x, θ, t) is governed by the Fokker-Plank equation which is a continuous version of Eq. (2) with the transition matrix L ij replaced by the operator for the rotational diffusion of the particle in the field. Analogously, the cumulant generating function α(λ) for the distribution of the displacement x is given by the maximal eigenvalue of the operator By rescaling to the dimensionless variable this maximal eigenvalue can be written as where a(B, z) is defined as the maximal value of a for which the differential equatioñ has a 2π-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] a(0, z) = 1 2 z 2 − 7 32 z 4 + 29 144 Due to the obvious symmetry all odd terms vanish in this expansion. From this expansion one can derive the low order cumulants of the distribution p(x, t) as in accordance with the results for C 2 and C 4 in Ref. [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.

A new symmetry
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 ψ(θ) ofL(B, z), the shifted function ψ(θ + π) is an eigenfunction of the adjoint operatorL † (B, −B − z), As a consequence, the eigenvalue a(B, z) satisfies the symmetry In terms of the cumulant generating function (36) the symmetry (42) can be expressed as with 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 with For the probability distribution p(x) of the displacement we obtain the fluctuation relation 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.

Numerical calculation of the rate function
For an efficient numerical calculation of the eigenvalue a(B, z) we expand the eigenfunction as ψ(θ) = c 0 + 2 ∞ k=1 c k cos(kθ), which leads to the representation of Eq. (37) as k L kk c k = ac k with the infinite matrix (up to the exceptionL 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 Fig. 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 Fig. 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 .   (52) and (61), respectively. The gray dotted line shows the obviously insufficient approximation a(B, z) = |z|.

Time-scale separation
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 r → 0 and D r → ∞, 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 Eq. (37). The relevant eigenvector of the term z cos θ for large positive z is simply δ(θ) with eigenvalue z. For large negative values of z the relevant eigenvector is δ(θ − π) with eigenvalue −z. Plugging this asymptotic behavior into Eq. (36) yields the cumulant generating function α(λ) = Dλ 2 + µF 0 |λ| (49) in the limit D r → 0, which is analogous to Eq. (23) in the discrete case. The resulting rate function is shown in red in Fig. 4. Although being asymptotically correct, in the sense that a(B, z)/|z| → 1 for large |z|, the approximation a(B, z) ∼ |z| is rather rough for realistic values of z. In particular, in Fig. 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 ψ(θ) 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 ψ(θ) for intermediate z somehow interpolates between the known eigenvector p s (θ) at z = 0 and δ(θ) for z → ∞, the ansatz Although this result is exact only for z = 0 and z → ∞, the numerical results in Fig. 3 show that there is a good agreement with the actual eigenvalue for all values of z > −B/2 for strong fields  (47) is marked by triangles. The limit Dr → 0 is shown as a thick red curve with dashed analytic continuations and the limit Dr → ∞ is shown as a thick blue curve.
where the eigenvectors are more concentrated around θ = 0. For z < −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 r → ∞. For an expansion in z, we use the stationary distribution p s (θ) = exp(B cos θ) 2πI 0 (B) (53) to obtain where the I n (B) denotes the modified Bessel functions of the first kind. Plugging this approximation into Eq. (36) yields the parabolic approximation of the cumulant generating function analogous to Eq. (28) and the corresponding rate function (29) with the velocity v = µF 0 I 1 (B)/I 0 (B) and diffusion coefficientD = D. In Fig. 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 α(λ) = vλ + Dλ 2 . As visible in the plots of the function a(B, z) in Fig. 3, the location of the kink coincides with the center of the symmetry (43) at λ = λ 0 . As explained for the discrete model in Sec. 2.3, this kink leads in the Legendre transformed rate function to an interval with constant slope λ 0 . Based on Eq. (27) and the limits −1 < ∂ z a(B, z) < 1, the range where this linear slope occurs can be limited to u 0 − µF 0 < u < u 0 + µ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.

Spherical Janus particles
For three dimensional, spherical Janus particles we write the orientation in spherical coordinates (θ, φ). A distribution over these coordinates is to be understood as a probability per solid angle, hence integrals over the azimuthal angels θ require sin θ as Jacobi determinant. For the joint probability distribution p(x, θ, φ, t), the rotational operator in the Fokker-Planck equation (32) gets modified to Eigenfunctions of the accordingly modified operator L(λ) in Eq. (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 ψ(θ) and thus drop the part of L r (B) acting on φ. The steps leading to the operatorL(B, z, θ) in Eq. (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 relatioñ As a consequence, the eigenvalue a(B, z) satisfies the symmetry and the property (43) of the cumulant generating function still holds, but with a different center of symmetry λ 0 = −D r B/(µF 0 ). A numerical calculation of the function a(B, z) is possible via an expansion of the eigenfunctions in Legendre polynomials P (cos θ). The matrix form ofL(B, z, θ), analogous to (48), then reads with , ≥ 0. The dominant eigenvalue can be determined numerically in a truncated approximation of this matrix. The result of the numerics shown in Fig. 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.
Accordingly, the average velocity of the Janus particle is given by v = µF 0 (coth B − 1/B). The effective field approximation analogous to Eq. (50) yields for z > −B which can be continued to z < −B via the symmetry (57).

General model for the active mechanism
In Sec. 2.1 we have used diffusion with drift as the arguably simplest possible dynamics for fixed internal state i. Molecular motors that proceed in a stepwise fashion do usually not exert a constant force as we have assumed in the simple model of Sec. 2.
2. An effective description in terms of the average velocity v i = µ i f i and effective diffusion coefficient D i covers moderate fluctuations quite well. For extreme fluctuations, however, the approximation of the rate function for fixed i as a parabola does not capture the previously reported rich structure of the rate function for molecular motors [51][52][53] modeled via ratchet models [1] or discrete stepping processes [54]. In a more detailed and general description, we assign to each of the internal mesostates (labeled by the Latin letters i, j, . . .) further microstates (labeled by Greek letters γ, δ, . . .), as shown in enzyme. For microswimmers, the mesostate i labels the swimming direction, while microstates γ label an arbitrarily fine description of the hydrodynamic flow field around the particle and possibly the geometric configuration of flagella. In a less detailed, effective description of active Brownian motion via a velocity dependent force [9,55,56], the microstates γ correspond to a finely discretized set of velocity states. We assume Markovian transition rates k i,δγ from microstate δ to microstate γ within the mesostate i. Associated with this transition is the displacement d i,δγ = −d i,γδ of the whole system along the coordinate x. The transitions between different mesostates are modeled via the Markovian rates ηw iγ,jδ from microstate γ of mesostate i to microstate δ of mesostate j. For simplicity, we assume these transitions to come without displacement. The dynamics of the system is governed by the master equation and The stationary distribution p s (i, γ) follows upon setting the left hand side of Eq. (62) to zero and can be used to calculate the average velocity The cumulant generating function α(λ) for the displacement is given by the maximal eigenvalue of L iγ,jδ (λ).
In the limit of slow transitions between mesostates, i.e. for η → 0, the cumulant generating function is given by where α i (λ) is the maximal eigenvalue of L i,γδ (λ). Kinks in the cumulant generating function α(λ) occur when the two maximal eigenvalues α i (λ) intersect. In the Legendre transformed rate function h(u), these kinks result in ranges with linear slope, forming the convex envelope of the These features are analogous to what we have described for the simple model with diffusion and drift in Sec. 2. Merely the parabolas in Eqs. (23) and (26) are replaced by the less trivial functions α i (λ) and h i (u), respectively, corresponding to the dynamics for fixed i. The simple model is recovered if the parabolic approximation of h i (u) around u = v extends to a wide range of the parameter u, which is the case in the linear response regime where the driving chemical and/or mechanical forces are weak and for small step-sizes d i,γδ .
For small non-zero values of η, the kinks in the cumulant generating function and the linear slope in the rate function smear out. For larger η, the rates k i,γδ and ηw iγ,jδ become comparable. Then, the distinction between meso-and microstates becomes pointless and one has to diagonalize the full matrix (63). In the limit η → ∞, one can simplify the network by coarse graining microstates that are connected via the transitions ηw iγ,jδ . Among these, a local stationary distribution is reached quickly, which can be used to calculate the effective transition rates in the coarse grained network based on the rates k i,γδ .
This generalized formalism can be illustrated by reconsidering the motor model with dichotomous drift velocity in Sec. 2.2. To account for the discrete stepping of the motor with step length d (e.g., d = 8 nm for kinesin [57]), we here model the motor as a random walker. In the simplest description, each of the two mesostates i ∈ {1, 2} contains only a single transition with forward rate k + i and backward rate k − i . These rates usually follow through coarse graining from a more detailed motor model [58,59]. Switching between the two mesostates occurs at rates w 12 and w 21 . The functions relevant for the η → 0 limit then read [35] and with , and a i ≡ sinh(A i /2). For η 1, we can make use of the analytical solution for the cumulant generating function as the maximal eigenvalue (19) of the matrix Finally, in the limit η → ∞, the cumulant generating function becomes with coarse-grained rates and the rate function analogous to Eq. (68). The cumulant generating function and rate function for selected values of the transition rates are shown in Fig. 5.

Summary and perspectives
In summary, we have shown that the rate function characterizing extreme fluctuations of the displacement in active Brownian motion displays a rich structure, which can be understood in a simple model describing the propulsion mechanism as a biased diffusion process. While the rate function for passive Brownian motion is globally parabolic, corresponding to the Gaussian shape of the distribution, we find for active Brownian motion an interplay between parabolic and linear sections. This structure of the rate function can be attributed to the dependence of the propulsion force on the mesoscopic state of the system. For fluctuations contributing to the parabolic sections, the internal mechanism typically rests in an individual mesostate, so that the propulsion force is constant and fluctuations are dominated by translational diffusion. In contrast, the fluctuations in the linear sections are dominated by the noise in the internal process, which results in a fluctuating propulsion force. If the internal process is slow compared to the time scale of the translational diffusion, the borders between the parabolic and linear sections of the rate function become sharp.
With increasing speed of the internal process the shape of the rate function gets more smeared out, until it approaches the globally parabolic rate function of a particle that diffuses with constant drift and an effective diffusion coefficient. For experimentally measured distributions, the identification of the different sections in the rate function could help to elucidate aspects of the hidden internal states of active systems. For example, fitting of the parabolic sections would yield information about the propulsion forces f i and diffusion coefficients D i of the corresponding dominant internal states.
For a propulsion mechanism that is more intricate than biased diffusion, like for stepping motors, we have shown that the general features of the extreme fluctuations remain the same except that the parabolic sections get replaced by the more complex rate functions corresponding to the different internal mesostates.
Although being characteristic for active Brownian motion, the occurrence of linear sections in the rate function is not restricted to this type of Brownian motion. For example, a passive but asymmetric particle pulled with a constant force has different mobilities in different geometrical alignments and hence also different drift velocities. This fact results in a similar rate function as we have described for active particles. Similarly, passive particles in a Poiseuille flow have different drift velocities depending on the radial position.
For the example of Janus particles with a dipole moment interacting with a homogeneous field, we have found a new fluctuation relation-like symmetry for the rate function. Unlike fluctuation relations related to entropy production, the center of this symmetry is not at velocity zero but at a negative, field-dependent velocity u 0 . In experiments, one could make use of this relation to extract the dipole moment of a Janus particle from the slope 2λ 0 . While this information could in principle be extracted from a measurement of the average velocity v as well, using the fluctuation relation would have the advantage of being robust to a spurious drift from possible external forces. These would only affect the center of the symmetry u 0 but not the slope. This distinction is relevant in experiments assessing gravitaxis, where gravity and buoyancy jointly cause both a torque and a translational force.
The fluctuation relation for the Janus particles is complementary to the fluctuation relation for entropy production also in the respect that it comes with a linear section of the rate function around the center of symmetry. In contrast, for the entropy production there is a kink at this center of symmetry [60,61].
Based on the insight from the discrete model, we expect that the occurrence of linear sections in the rate function for the displacement is a universal feature of Janus particles and independent of the geometry and the type of interaction with external fields. However, the new fluctuation relation we have proven for the dipolar interaction would not extend to Janus particles with broken spatial symmetry, as a numerical check for Janus particles in an asymmetric periodic potential reveals. In future work, it will be a challenging task to figure out universal conditions under which currents in arbitrary systems exhibit fluctuation relations with a non-zero center of symmetry.